From d6acae5330c684b66a99d0fbf13e44b0b66046cf Mon Sep 17 00:00:00 2001 From: Chaithya G R Date: Thu, 9 Jul 2020 14:21:24 +0200 Subject: [PATCH] Add support for density compensation (#109) * Add support for density compensation * Fix Documentation and PEP8 * Update to utils.py and README * Update tests for sensitivity extraction * remove warning for pyNUFFT and make it an error * Update based on comments from Zac * Revert a single change to prevent unwanted changes * Final changes * Make den_c None in case not used * Fix Doc Co-authored-by: chaithyagr --- README.rst | 2 +- ...nCartesian_gpuNUFFT_DensityCompensation.py | 78 +++++++++++++++++ ...sian_gpuNUFFT_SENSE_DensityCompensation.py | 86 +++++++++++++++++++ mri/operators/fourier/non_cartesian.py | 57 ++++++++---- mri/operators/fourier/utils.py | 33 +++++++ .../utils/extract_sensitivity_maps.py | 55 +++++++++--- mri/test_local/test_fourier_adjoint_gpu.py | 15 ++-- mri/tests/test_fourier_adjoint.py | 15 ++-- mri/tests/test_sensitivity_extraction.py | 41 ++++++--- 9 files changed, 332 insertions(+), 50 deletions(-) create mode 100644 examples/GPU_Examples/NonCartesian_gpuNUFFT_DensityCompensation.py create mode 100644 examples/GPU_Examples/NonCartesian_gpuNUFFT_SENSE_DensityCompensation.py diff --git a/README.rst b/README.rst index ebabc211..2d5f0922 100644 --- a/README.rst +++ b/README.rst @@ -38,7 +38,7 @@ Special Installations For faster NUFFT operation, pysap-mri uses gpuNUFFT, to run the NUFFT on GPU. To install gpuNUFFT, please use: -``pip install git+https://github.com/andyschwarzl/gpuNUFFT`` +``pip install gpuNUFFT`` We are still in the process of merging this and getting a pip release for gpuNUFFT. Note, that you can still use CPU NFFT without installing the package. diff --git a/examples/GPU_Examples/NonCartesian_gpuNUFFT_DensityCompensation.py b/examples/GPU_Examples/NonCartesian_gpuNUFFT_DensityCompensation.py new file mode 100644 index 00000000..2f219c70 --- /dev/null +++ b/examples/GPU_Examples/NonCartesian_gpuNUFFT_DensityCompensation.py @@ -0,0 +1,78 @@ +""" +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/examples/GPU_Examples/NonCartesian_gpuNUFFT_SENSE_DensityCompensation.py b/examples/GPU_Examples/NonCartesian_gpuNUFFT_SENSE_DensityCompensation.py new file mode 100644 index 00000000..702272bc --- /dev/null +++ b/examples/GPU_Examples/NonCartesian_gpuNUFFT_SENSE_DensityCompensation.py @@ -0,0 +1,86 @@ +""" +Neuroimaging non-cartesian reconstruction +========================================= + +Author: Chaithya G R + +In this tutorial we will reconstruct an MR Image directly with density +compensation and SENSE from gpuNUFFT + +Import neuroimaging data +------------------------ + +We use the toy datasets available in pysap, more specifically a 3D orange data +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 +from mri.reconstructors.utils.extract_sensitivity_maps import get_Smaps +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('3d-pmri') +cartesian = np.linalg.norm(image, axis=0) + +# Obtain MRI non-cartesian mask and estimate the density compensation +radial_mask = get_sample_data("mri-radial-3d-samples") +kspace_loc = radial_mask.data +density_comp = estimate_density_compensation(kspace_loc, cartesian.shape) + +############################################################################# +# Generate the kspace +# ------------------- +# +# From the 3D orange slice and 3D radial acquisition mask, we retrospectively +# undersample the k-space +# 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=cartesian.shape, + n_coils=image.shape[0], + implementation='gpuNUFFT', +) +kspace_obs = fourier_op.op(image.data) +Smaps, SOS = get_Smaps( + k_space=kspace_obs, + img_shape=fourier_op.shape, + samples=kspace_loc, + thresh=(0.05, 0.05, 0.05), # The cutoff threshold in each kspace + # direction between 0 and kspace_max (0.5) + min_samples=kspace_loc.min(axis=0), + max_samples=kspace_loc.max(axis=0), + density_comp=density_comp, + mode='NFFT', +) +fourier_op_sense_dc = NonCartesianFFT( + samples=kspace_loc, + shape=cartesian.shape, + implementation='gpuNUFFT', + n_coils=image.shape[0], + density_comp=density_comp, + smaps=Smaps, +) + +# Density Compensation SENSE adjoint: +# This preconditions k-space giving a result closer to inverse +image_rec1 = pysap.Image(data=np.abs( + fourier_op_sense_dc.adj_op(kspace_obs)) +) +# image_rec1.show() +base_ssim = ssim(image_rec1, cartesian) +print('The SSIM for simple Density ' + 'compensated SENSE Adjoint is : ' + str(base_ssim)) diff --git a/mri/operators/fourier/non_cartesian.py b/mri/operators/fourier/non_cartesian.py index cd1441b8..d4e04627 100644 --- a/mri/operators/fourier/non_cartesian.py +++ b/mri/operators/fourier/non_cartesian.py @@ -27,12 +27,13 @@ warnings.warn("pynfft python package has not been found. If needed use " "the master release.") pass +pynufft_available = False try: from pynufft import NUFFT_hsa, NUFFT_cpu except Exception: - warnings.warn("pynufft python package has not been found. If needed use " - "the master release. Till then you cannot use NUFFT on GPU") pass +else: + pynufft_available = True gpunufft_available = False try: @@ -236,6 +237,10 @@ def __init__(self, samples, shape, platform='cuda', Kd=None, Jd=None, """ if (n_coils < 1) or (type(n_coils) is not int): raise ValueError('The number of coils should be an integer >= 1') + if not pynufft_available: + raise ValueError('PyNUFFT Package is not installed, please ' + 'consider using `gpuNUFFT` or install the ' + 'PyNUFFT package') self.nb_coils = n_coils self.shape = shape self.platform = platform @@ -464,7 +469,7 @@ def __init__(self, samples, shape, n_coils=1, density_comp=None, balance_workload ) - def op(self, image): + def op(self, image, interpolate_data=False): """ This method calculates the masked non-cartesian Fourier transform of a 2D / 3D image. @@ -472,6 +477,9 @@ def op(self, image): ---------- image: np.ndarray input array with the same shape as shape. + interpolate_data: bool, default False + if set to True, the image is just apodized and interpolated to + kspace locations. This is used for density estimation. Returns ------- @@ -484,16 +492,19 @@ def op(self, image): if self.n_coils > 1 and not self.uses_sense: coeff = self.operator.op(np.asarray( [np.reshape(image_ch.T, image_ch.size) for image_ch in image] - ).T) + ).T, interpolate_data) else: - coeff = self.operator.op(np.reshape(image.T, image.size)) + coeff = self.operator.op( + np.reshape(image.T, image.size), + interpolate_data + ) # Data is always returned as num_channels X coeff_array, # so for single channel, we extract single array if not self.uses_sense: coeff = coeff[0] return coeff - def adj_op(self, coeff): + def adj_op(self, coeff, grid_data=False): """ This method calculates adjoint of non-uniform Fourier transform of a 1-D coefficients array. @@ -501,14 +512,16 @@ def adj_op(self, coeff): ---------- coeff: np.ndarray masked non-uniform Fourier transform 1D data. - + grid_data: bool, default False + if True, the kspace data is gridded and returned, + this is used for density compensation Returns ------- np.ndarray adjoint operator of Non Uniform Fourier transform of the input coefficients. """ - image = self.operator.adj_op(coeff) + image = self.operator.adj_op(coeff, grid_data) if self.n_coils > 1 and not self.uses_sense: image = np.asarray( [image_ch.T for image_ch in image] @@ -523,7 +536,7 @@ def adj_op(self, coeff): class NonCartesianFFT(OperatorBase): """This class wraps around different implementation algorithms for NFFT""" def __init__(self, samples, shape, implementation='cpu', n_coils=1, - **kwargs): + density_comp=None, **kwargs): """ Initialize the class. Parameters @@ -548,9 +561,11 @@ def __init__(self, samples, shape, implementation='cpu', n_coils=1, self.samples = samples self.n_coils = n_coils if implementation == 'cpu': + self.density_comp = density_comp self.implementation = NFFT(samples=samples, shape=shape, n_coils=self.n_coils) elif implementation == 'cuda' or implementation == 'opencl': + self.density_comp = density_comp self.implementation = NUFFT(samples=samples, shape=shape, platform=implementation, n_coils=self.n_coils) @@ -559,14 +574,19 @@ def __init__(self, samples, shape, implementation='cpu', n_coils=1, raise ValueError('gpuNUFFT library is not installed, ' 'please refer to README' 'or use cpu for implementation') - self.implementation = gpuNUFFT(samples=samples, shape=shape, - n_coils=self.n_coils, **kwargs) + self.implementation = gpuNUFFT( + samples=samples, + shape=shape, + n_coils=self.n_coils, + density_comp=density_comp, + **kwargs + ) else: raise ValueError('Bad implementation ' + implementation + ' chosen. Please choose between "cpu" | "cuda" |' '"opencl" | "gpuNUFFT"') - def op(self, data): + def op(self, data, *args): """ This method calculates the masked non-cartesian Fourier transform of an image. @@ -579,9 +599,9 @@ def op(self, data): ------- masked Fourier transform of the input image. """ - return self.implementation.op(data) + return self.implementation.op(data, *args) - def adj_op(self, coeffs): + def adj_op(self, coeffs, *args): """ This method calculates inverse masked non-uniform Fourier transform of a 1-D coefficients array. @@ -594,7 +614,14 @@ def adj_op(self, coeffs): ------- inverse discrete Fourier transform of the input coefficients. """ - return self.implementation.adj_op(coeffs) + if not isinstance(self.implementation, gpuNUFFT) and \ + self.density_comp is not None: + return self.implementation.adj_op( + coeffs * self.density_comp, + *args + ) + else: + return self.implementation.adj_op(coeffs, *args) class Stacked3DNFFT(OperatorBase): diff --git a/mri/operators/fourier/utils.py b/mri/operators/fourier/utils.py index 1541a373..d234575d 100644 --- a/mri/operators/fourier/utils.py +++ b/mri/operators/fourier/utils.py @@ -308,3 +308,36 @@ def check_if_fourier_op_uses_sense(fourier_op): return fourier_op.implementation.uses_sense else: return False + + +def estimate_density_compensation(kspace_loc, volume_shape, num_iterations=10): + """ Utils function to obtain the density compensator for a + given set of kspace locations. + + Parameters + ---------- + kspace_loc: np.ndarray + the kspace locations + volume_shape: np.ndarray + the volume shape + num_iterations: int default 10 + the number of iterations for density estimation + """ + from .non_cartesian import NonCartesianFFT + from .non_cartesian import gpunufft_available + if gpunufft_available is False: + raise ValueError("gpuNUFFT is not available, cannot " + "estimate the density compensation") + grid_op = NonCartesianFFT( + samples=kspace_loc, + shape=volume_shape, + implementation='gpuNUFFT', + osf=1, + ) + density_comp = np.ones(kspace_loc.shape[0]) + for _ in range(num_iterations): + density_comp = ( + density_comp / + np.abs(grid_op.op(grid_op.adj_op(density_comp, True), True)) + ) + return density_comp diff --git a/mri/reconstructors/utils/extract_sensitivity_maps.py b/mri/reconstructors/utils/extract_sensitivity_maps.py index 53bb1677..98be1525 100644 --- a/mri/reconstructors/utils/extract_sensitivity_maps.py +++ b/mri/reconstructors/utils/extract_sensitivity_maps.py @@ -15,6 +15,7 @@ # Package import from mri.operators import NonCartesianFFT +from mri.operators.fourier.non_cartesian import gpunufft_available from mri.operators.utils import get_stacks_fourier, \ gridded_inverse_fourier_transform_nd, \ gridded_inverse_fourier_transform_stack, convert_locations_to_mask @@ -27,7 +28,7 @@ def extract_k_space_center_and_locations(data_values, samples_locations, thr=None, img_shape=None, - is_fft=False): + is_fft=False, density_comp=None): """ This class extract the k space center for a given threshold and extracts the corresponding sampling locations @@ -45,9 +46,15 @@ def extract_k_space_center_and_locations(data_values, samples_locations, is_fft: bool default False Checks if the incoming data is from FFT, in which case, masking can be done more directly + density_comp: np.ndarray default None + The density compensation for kspace data in case it exists and we + use density compensated adjoint for Smap estimation Returns ------- - The extracted center of the k-space + The extracted center of the k-space, i.e. both the kspace locations and + kspace values. If the density compensators are passed, the corresponding + compensators for the center of k-space data will also be returned. The + return stypes for density compensation and kspace data is same as input """ if thr is None: if img_shape is None: @@ -76,12 +83,17 @@ def extract_k_space_center_and_locations(data_values, samples_locations, index = np.extract(condition, index) center_locations = samples_locations[index, :] data_thresholded = data_ordered[:, index] - return data_thresholded, center_locations + if density_comp is not None: + density_comp = density_comp[index] + return data_thresholded, center_locations, density_comp + else: + return data_thresholded, center_locations def get_Smaps(k_space, img_shape, samples, thresh, min_samples, max_samples, mode='Gridding', - method='linear', n_cpu=1): + method='linear', density_comp=None, n_cpu=1, + fourier_op_kwargs=None): """ This method estimate the sensitivity maps information from parallel mri acquisition and for variable density sampling scheme where teh k-space @@ -111,8 +123,15 @@ def get_Smaps(k_space, img_shape, samples, thresh, been sampled on the grid method: string 'linear' | 'cubic' | 'nearest', default='linear' For gridding mode, it defines the way interpolation must be done - n_cpu: intm default=1 + density_comp: np.ndarray default None + The density compensation for kspace data in case it exists and we + use density compensated adjoint for Smap estimation + n_cpu: int default=1 Number of parallel jobs in case of parallel MRI + fourier_op_kwargs: dict, default None + The keyword arguments given to fourier_op initialization if + mode == 'NFFT'. If None, we choose implementation of fourier op to + 'gpuNUFFT' if library is installed. Returns ------- @@ -127,13 +146,19 @@ def get_Smaps(k_space, img_shape, samples, thresh, or len(thresh) != len(img_shape): raise NameError('The img_shape, max_samples, ' 'min_samples and thresh must be of same length') - k_space, samples = \ + k_space, samples, *density_comp = \ extract_k_space_center_and_locations( data_values=k_space, samples_locations=samples, thr=thresh, img_shape=img_shape, - is_fft=mode == 'FFT') + is_fft=mode == 'FFT', + density_comp=density_comp, + ) + if density_comp: + density_comp = density_comp[0] + else: + density_comp = None if samples is None: mode = 'FFT' L, M = k_space.shape @@ -148,9 +173,19 @@ def get_Smaps(k_space, img_shape, samples, thresh, for l in range(Smaps_shape[2]): Smaps[l] = pfft.ifftshift(pfft.ifft2(pfft.fftshift(k_space[l]))) elif mode == 'NFFT': - fourier_op = NonCartesianFFT(samples=samples, shape=img_shape, - implementation='cpu') - Smaps = np.asarray([fourier_op.adj_op(k_space[l]) for l in range(L)]) + if fourier_op_kwargs is None: + if gpunufft_available: + fourier_op_kwargs = {'implementation': 'gpuNUFFT'} + else: + fourier_op_kwargs = {} + fourier_op = NonCartesianFFT( + samples=samples, + shape=img_shape, + density_comp=density_comp, + n_coils=L, + **fourier_op_kwargs, + ) + Smaps = fourier_op.adj_op(np.ascontiguousarray(k_space)) elif mode == 'Stack': grid_space = [np.linspace(min_samples[i], max_samples[i], diff --git a/mri/test_local/test_fourier_adjoint_gpu.py b/mri/test_local/test_fourier_adjoint_gpu.py index e392c825..5d79bc73 100644 --- a/mri/test_local/test_fourier_adjoint_gpu.py +++ b/mri/test_local/test_fourier_adjoint_gpu.py @@ -62,14 +62,17 @@ def test_NUFFT_2D(self): for platform in self.platforms: _mask = np.random.randint(2, size=(self.N, self.N)) _samples = convert_mask_to_locations(_mask) - fourier_op_adj = NonCartesianFFT(samples=_samples, - shape=(self.N, self.N), - implementation=platform, - n_coils=num_channels) + fourier_op_adj = NonCartesianFFT( + samples=_samples, + shape=(self.N, self.N), + implementation=platform, + n_coils=num_channels, + density_comp=np.ones((num_channels, _samples.shape[0])) + ) Img = (np.random.randn(num_channels, self.N, self.N) + 1j * np.random.randn(num_channels, self.N, self.N)) - f = (np.random.randn(num_channels, _samples.shape[0], 1) + - 1j * np.random.randn(num_channels, _samples.shape[0], 1)) + f = (np.random.randn(num_channels, _samples.shape[0]) + + 1j * np.random.randn(num_channels, _samples.shape[0])) f_p = fourier_op_adj.op(Img) I_p = fourier_op_adj.adj_op(f) x_d = np.vdot(Img, I_p) diff --git a/mri/tests/test_fourier_adjoint.py b/mri/tests/test_fourier_adjoint.py index 253e9985..e3334a4d 100644 --- a/mri/tests/test_fourier_adjoint.py +++ b/mri/tests/test_fourier_adjoint.py @@ -113,7 +113,9 @@ def test_FFT(self): print(" FFT adjoint test passes") def test_NFFT_2D(self): - """Test the adjoint operator for the 2D non-Cartesian Fourier transform + """Test the adjoint operator for the 2D non-Cartesian Fourier + transform, with density compensator set to 1, to vet the code + path, the test is unchanged otherwise """ for num_channels in self.num_channels: print("Testing with num_channels=" + str(num_channels)) @@ -121,10 +123,13 @@ def test_NFFT_2D(self): _mask = np.random.randint(2, size=(self.N, self.N)) _samples = convert_mask_to_locations(_mask) print("Process NFFT in 2D test '{0}'...", i) - fourier_op = NonCartesianFFT(samples=_samples, - shape=(self.N, self.N), - n_coils=num_channels, - implementation='cpu') + fourier_op = NonCartesianFFT( + samples=_samples, + shape=(self.N, self.N), + n_coils=num_channels, + implementation='cpu', + density_comp=np.ones((num_channels, _samples.shape[0])) + ) Img = np.random.randn(num_channels, self.N, self.N) + \ 1j * np.random.randn(num_channels, self.N, self.N) f = np.random.randn(num_channels, _samples.shape[0]) + \ diff --git a/mri/tests/test_sensitivity_extraction.py b/mri/tests/test_sensitivity_extraction.py index 9f87d921..38aaf4dc 100644 --- a/mri/tests/test_sensitivity_extraction.py +++ b/mri/tests/test_sensitivity_extraction.py @@ -63,9 +63,10 @@ def test_extract_k_space_center_3D(self): data_thresholded) def test_extract_k_space_center_2D(self): - """ Ensure that the extracted k-space center is right""" - _mask = np.ones((self.N, self.N)) - _samples = convert_mask_to_locations(_mask) + """ Ensure that the extracted k-space center is right, + send density compensation also and vet the code path""" + mask = np.ones((self.N, self.N)) + samples = convert_mask_to_locations(mask) Img = (np.random.randn(self.num_channel, self.N, self.N) + 1j * np.random.randn(self.num_channel, self.N, self.N)) Nby2_percent = self.N * self.percent / 2 @@ -73,13 +74,15 @@ def test_extract_k_space_center_2D(self): high = int(self.N / 2 + Nby2_percent + 1) center_Img = Img[:, low:high, low:high] thresh = self.percent * 0.5 - data_thresholded, samples_thresholded = \ + data_thresholded, samples_thresholded, dc = \ extract_k_space_center_and_locations( data_values=np.reshape(Img, (self.num_channel, self.N * self.N)), - samples_locations=_samples, + samples_locations=samples, thr=(thresh, thresh), - img_shape=(self.N, self.N)) + img_shape=(self.N, self.N), + density_comp=np.ones(samples.shape[0]) + ) np.testing.assert_allclose( center_Img.reshape(data_thresholded.shape), data_thresholded) @@ -114,9 +117,9 @@ def test_sensitivity_extraction_2D(self): """ This test ensures that the output of the non cartesian kspace extraction is same a that of mimicked cartesian extraction in 2D """ - _mask = np.ones((self.N, self.N)) - _samples = convert_mask_to_locations(_mask) - fourier_op = NonCartesianFFT(samples=_samples, shape=(self.N, self.N)) + mask = np.ones((self.N, self.N)) + samples = convert_mask_to_locations(mask) + fourier_op = NonCartesianFFT(samples=samples, shape=(self.N, self.N)) Img = (np.random.randn(self.num_channel, self.N, self.N) + 1j * np.random.randn(self.num_channel, self.N, self.N)) F_img = np.asarray([fourier_op.op(Img[i]) @@ -124,20 +127,32 @@ def test_sensitivity_extraction_2D(self): Smaps_gridding, SOS_Smaps = get_Smaps( k_space=F_img, img_shape=(self.N, self.N), - samples=_samples, + samples=samples, thresh=(0.4, 0.4), mode='gridding', min_samples=(-0.5, -0.5), max_samples=(0.5, 0.5), n_cpu=1) + Smaps_NFFT_dc, SOS_Smaps_dc = get_Smaps( + k_space=F_img, + img_shape=(self.N, self.N), + thresh=(0.4, 0.4), + samples=samples, + min_samples=(-0.5, -0.5), + max_samples=(0.5, 0.5), + mode='NFFT', + density_comp=np.ones(samples.shape[0]) + ) Smaps_NFFT, SOS_Smaps = get_Smaps( k_space=F_img, img_shape=(self.N, self.N), thresh=(0.4, 0.4), - samples=_samples, + samples=samples, min_samples=(-0.5, -0.5), max_samples=(0.5, 0.5), - mode='NFFT') + mode='NFFT', + ) + np.testing.assert_allclose(Smaps_gridding, Smaps_NFFT_dc) np.testing.assert_allclose(Smaps_gridding, Smaps_NFFT) # Test that we raise assert for bad mode np.testing.assert_raises( @@ -146,7 +161,7 @@ def test_sensitivity_extraction_2D(self): k_space=F_img, img_shape=(self.N, self.N), thresh=(0.4, 0.4), - samples=_samples, + samples=samples, min_samples=(-0.5, -0.5), max_samples=(0.5, 0.5), mode='test'