diff --git a/examples/GPU_Examples/NonCartesian_gpuNUFFT_DensityCompensation.py b/examples/GPU_Examples/NonCartesian_gpuNUFFT_DensityCompensation.py
deleted file mode 100644
index 2f219c70..00000000
--- a/examples/GPU_Examples/NonCartesian_gpuNUFFT_DensityCompensation.py
+++ /dev/null
@@ -1,78 +0,0 @@
-"""
-Neuroimaging non-cartesian reconstruction
-=========================================
-
-Author: Chaithya G R
-
-In this tutorial we will reconstruct an MR Image directly with density
-compensation
-
-Import neuroimaging data
-------------------------
-
-We use the toy datasets available in pysap, more specifically a 2D brain slice
-and the radial acquisition scheme (non-cartesian).
-"""
-
-# Package import
-from mri.operators import NonCartesianFFT, WaveletUD2
-from mri.operators.utils import convert_locations_to_mask, \
-    gridded_inverse_fourier_transform_nd
-from mri.operators.fourier.utils import estimate_density_compensation
-from mri.reconstructors import SingleChannelReconstructor
-import pysap
-from pysap.data import get_sample_data
-
-# Third party import
-from modopt.math.metrics import ssim
-from modopt.opt.linear import Identity
-from modopt.opt.proximity import SparseThreshold
-import numpy as np
-
-# Loading input data
-image = get_sample_data('2d-mri')
-
-# Obtain MRI non-cartesian mask and estimate the density compensation
-radial_mask = get_sample_data("mri-radial-samples")
-kspace_loc = radial_mask.data
-density_comp = estimate_density_compensation(kspace_loc, image.shape)
-
-#############################################################################
-# Generate the kspace
-# -------------------
-#
-# From the 2D brain slice and the acquisition mask, we retrospectively
-# undersample the k-space using a radial acquisition mask
-# We then reconstruct using adjoint with and without density compensation
-
-# Get the locations of the kspace samples and the associated observations
-fourier_op = NonCartesianFFT(
-    samples=kspace_loc,
-    shape=image.shape,
-    implementation='gpuNUFFT',
-)
-fourier_op_density_comp = NonCartesianFFT(
-    samples=kspace_loc,
-    shape=image.shape,
-    implementation='gpuNUFFT',
-    density_comp=density_comp
-)
-# Get the kspace data retrospectively. Note that this can be done with
-# `fourier_op_density_comp` as the forward operator is the same
-kspace_obs = fourier_op.op(image.data)
-
-# Simple adjoint
-image_rec0 = pysap.Image(data=np.abs(fourier_op.adj_op(kspace_obs)))
-# image_rec0.show()
-base_ssim = ssim(image_rec0, image)
-print('The SSIM from Adjoint is : ' + str(base_ssim))
-
-# Density Compensation adjoint:
-# This preconditions k-space giving a result closer to inverse
-image_rec1 = pysap.Image(data=np.abs(
-    fourier_op_density_comp.adj_op(kspace_obs))
-)
-# image_rec1.show()
-new_ssim = ssim(image_rec1, image)
-print('The SSIM from Density '
-      'compensated Adjoint is : ' + str(new_ssim))
diff --git a/mri/operators/fourier/non_cartesian.py b/mri/operators/fourier/non_cartesian.py
index c40079c3..63f988ba 100644
--- a/mri/operators/fourier/non_cartesian.py
+++ b/mri/operators/fourier/non_cartesian.py
@@ -56,6 +56,7 @@ def __init__(self, samples, shape, implementation='finufft', n_coils=1,
         self.impl = get_operator(self.implementation)(
             self.samples,
             self.shape,
+            density=self.density_comp,
             **self.kwargs
         )
         
diff --git a/mri/operators/fourier/utils/processing.py b/mri/operators/fourier/utils/processing.py
index 397e862b..a2f40a8a 100644
--- a/mri/operators/fourier/utils/processing.py
+++ b/mri/operators/fourier/utils/processing.py
@@ -285,4 +285,4 @@ def estimate_density_compensation(kspace_loc, volume_shape, implementation='pipe
         volume_shape,
         **kwargs
     )
-    return density_comp
+    return np.squeeze(density_comp)