MRICufiNUFFT#

class mrinufft.operators.MRICufiNUFFT(samples, shape, density=False, n_coils=1, n_batchs=1, smaps=None, smaps_cached=False, verbose=False, squeeze_dims=False, n_trans=1, async_transfer=False, **kwargs)[source]#

Bases: _GramOpGpuMixin, FourierOperatorBase, _ToggleGradPlanMixin

MRI Transform operator, build around cufinufft.

This operator adds density estimation and compensation (preconditioning) and multicoil support.

Parameters:

samples (np.ndarray or GPUArray.) – The samples location of shape Nsamples x N_dimensions.
shape (tuple) – Shape of the image space.
n_coils (int) – Number of coils.
n_batchs (int) – Size of the batch dimension.
density (bool or array) –
Density compensation support.
- If array, use this for density compensation
- If True, the density compensation will be automatically estimated, using the fixed point method.
- If False, density compensation will not be used.
smaps (np.ndarray or GPUArray , optional) –
- If None: no Smaps wil be used.
- If np.ndarray: Smaps will be copied on the device, according to smaps_cached.
- If GPUArray, the smaps are already cached.
smaps_cached (bool, default False) –
- If False the smaps are copied on device and free at each iterations.
- If True, the smaps are copied on device and stay on it.
squeeze_dims (bool, default False) – If True, will try to remove the singleton dimension for batch and coils.
n_trans (int, default 1) – Number of transform to perform in parallel by cufinufft.
async_transfer (bool, default False) – If True, pipeline host<->device transfers with compute (double buffering) in the host-input code paths, using dedicated non-blocking CUDA streams. This overlaps H2D/D2H copies with cufinufft compute across batches, at the cost of extra pinned host buffers and device buffers (roughly 2x the per-batch memory).
kwargs – Extra kwargs for the raw cufinufft operator

Notes

Cufinufft is able to run multiple transform in parallel, this is controlled by the n_trans parameter. The data provided should be of shape, (n_batch, n_coils, img_shape) for op (type2) and (n_batch, n_coils, n_samples) for adjoint (type1). and in contiguous memory order.

For now only single precision (float32 and complex64) is supported

Examples using `mrinufft.operators.MRICufiNUFFT`#

Density Compensation Routines

Least Squares Image Reconstruction

Model-based iterative reconstruction

Sensitivity maps estimation

backend: ClassVar[str] = 'cufinufft'[source]#

available: ClassVar[bool] = True[source]#

property shape: tuple[int, ...][source]#: Shape of the image space of the operator.

property n_batchs[source]#: Number of coils for the operator.

property n_coils: int[source]#: Number of coils for the operator.

autograd_available = True[source]#

property dtype[source]#: Return floating precision of the operator.

property smaps[source]#: Sensitivity maps of the operator.

update_samples(new_samples, *, unsafe=False)[source]#

Update the samples of the NUFFT operator.

Parameters:

new_samples (np.ndarray or GPUArray) – The new samples location of shape Nsamples x N_dimensions.
unsafe (bool, default False) – If True, the original array is used directly without any checks. This should be used with caution as it might lead to unexpected behavior.

Notes

If unsafe is True, the new_samples should be of shape (Nsamples, N_dimensions), F-ordered (column-major) and in the range [-pi, pi]. If not, this will lead to unexpected behavior. You have been warned.

If unsafe is False, this is automatically handled.

property density: ndarray[tuple[Any, ...], dtype[floating]] | None[source]#: Density compensation of the operator.

op(data, ksp_d=None)[source]#

Non Cartesian MRI forward operator.

Parameters:

data (np.ndarray or GPUArray)
space. (The uniform (2D or 3D) data in image)

Return type:

Results array on the same device as data.

Notes

this performs for every coil ell: ..math:: mathcal{F}mathcal{S}_ell x

Note

This function uses numpy for all CPU arrays, and cupy for all on-gpu array. It will convert all its array argument to the respective array library. The outputs will be converted back to the original array module and device.

_get_async_streams()[source]#

Return (and lazily create) the dedicated H2D/D2H transfer streams.

Compute keeps running on the default stream (where cufinufft’s Plans already run, since a Plan’s CUDA stream cannot be changed after construction). These streams are created non-blocking so they do not implicitly synchronize with the default stream, allowing real overlap of transfers with compute.

_host_register(arr, anchor=None)[source]#

Page-lock ARR’s memory in place (no copy), for async H2D/D2H.

cudaHostRegister/cudaHostUnregister are expensive (roughly size-proportional: ~230ms combined for a 4GiB array), so the registration is cached by memory address and reused across calls as long as the underlying buffer is still alive, instead of paying that cost on every op/adj_op call. It is released automatically via a weakref.finalize callback once anchor (the array whose lifetime owns that memory – arr itself if not given explicitly) is garbage collected.

Returns True on success. On failure (e.g. the memory cannot be pinned), returns False so the caller can fall back to the synchronous path rather than silently doing unsafe transfers.

_register_contiguous(arr: ndarray[tuple[Any, ...], dtype[_ScalarT]], shape: tuple[int, ...] | None = None) → ndarray[tuple[Any, ...], dtype[_ScalarT]] | None[source]#

Page-lock a contiguous view of arr, reshaped if shape given.

Fuses making arr contiguous with _host_register: registration is keyed off the contiguous buffer’s own address, anchored on whichever array owns that memory – arr itself if it was already contiguous (no copy made), or the freshly made contiguous copy otherwise. Returns None (caller falls back to the synchronous path) if the memory cannot be page-locked.

Return type:: ndarray[tuple[Any, …], dtype[_ScalarT]] | None

_op_sense_device(data, ksp_d=None)[source]#

_op_sense_host(data, ksp=None)[source]#

_op_sense_host_async(data, ksp)[source]#

Pipelined (double-buffered) forward SENSE op for host data.

Registers data/ksp in place (no staging copy) and issues H2D/compute/D2H in an overlapped double-buffered pipeline. Returns None (caller falls back to the synchronous loop) if the arrays cannot be page-locked.

_op_calibless_device(data, ksp_d=None)[source]#

_op_calibless_host(data, ksp=None)[source]#

_op_calibless_host_async(data, ksp)[source]#

Pipelined (double-buffered) forward calibrationless op for host data.

_op(image_d, coeffs_d)[source]#: Low level operator implementation.

adj_op(coeffs, img_d=None)[source]#

Non Cartesian MRI adjoint operator.

Parameters:: coeffs (np.array or GPUArray)
Return type:: Array in the same memory space of coeffs. (ie on cpu or gpu Memory).

Note

_adj_op_sense_device(coeffs, img_d=None)[source]#: Perform sense reconstruction when data is on device.

_adj_op_sense_host(coeffs, img_d=None)[source]#

Perform sense reconstruction when data is on host.

On device the following array are involved: - coil_img(S, T, 1, X,Y,Z) - ksp_batch(B, 1, X,Y,Z) - smaps_batched(S, T, X,Y,Z)

_adj_op_sense_host_async(coeffs, img_d)[source]#

Pipelined (double-buffered) adjoint SENSE op for host data.

Only the input (ksp + smaps) side is pipelined: the output accumulates in-place into the shared img_d across every coil, so a single readback after the loop is used instead of a per-batch D2H copy. Registers coeffs in place (no staging copy). Returns None (caller falls back to the synchronous loop) if it cannot be page-locked.

_adj_op_calibless_device(coeffs, img_d=None)[source]#

_adj_op_calibless_host(coeffs, img_batched=None)[source]#

_adj_op_calibless_host_async(coeffs, img)[source]#

Pipelined (double-buffered) adjoint calibrationless op for host data.

Registers coeffs/img in place (no staging copy) and issues H2D/compute/D2H in an overlapped double-buffered pipeline. Returns None (caller falls back to the synchronous loop) if the arrays cannot be page-locked.

_adj_op(coeffs_d, image_d)[source]#: Low level adjoint operator implementation.

data_consistency(image_data, obs_data)[source]#

Compute the data consistency estimation directly on gpu.

This mixes the op and adj_op method to perform F_adj(F(x-y)) on a per coil basis. By doing the computation coil wise, it uses less memory than the naive call to adj_op(op(x)-y)

Parameters:

image (array) – Image on which the gradient operation will be evaluated. N_coil x Image shape is not using sense.
obs_data (array) – Observed data.

_dc_sense_host(image_data, obs_data)[source]#: Gradient computation when all data is on host.

_dc_sense_device(image_data, obs_data)[source]#: Gradient computation when all data is on device.

_dc_calibless_host(image_data, obs_data)[source]#: Calibrationless Gradient computation when all data is on host.

_dc_calibless_device(image_data, obs_data)[source]#: Calibrationless Gradient computation when all data is on device.

property eps[source]#: Return the underlying precision parameter.

property bsize_ksp[source]#: Size in Bytes of the compute batch of samples.

property bsize_img[source]#: Size in Bytes of the compute batch of images.

property img_size[source]#: Image size in bytes.

property ksp_size[source]#: k-space size in bytes.

property norm_factor[source]#: Norm factor of the operator.

_abc_impl = <_abc._abc_data object>[source]#

_accumulate_coil_combine(img_d: ndarray[tuple[Any, ...], dtype[_ScalarT]], i: int, data_batched: ndarray[tuple[Any, ...], dtype[_ScalarT]], smaps_batched: ndarray[tuple[Any, ...], dtype[_ScalarT]])[source]#

img_d[b] += sum_t data_batched[t] * conj(smaps_batched[t]).

A single custom kernel: conj, multiply and the reduction over the T coils all happen in one pass over memory, with no intermediate arrays and no extra kernel launches.

Since n_trans is required to divide n_coils, a chunk never straddles a batch boundary: the batch index for every one of the T coils in chunk i is the same single value (i*T)//n_coils, computed on the host as a plain Python int (no device array/upload needed at all).

_coil_slice(i)[source]#

Return the coil slice for chunk i of the (B*C)//T loop.

n_trans (T) is required to divide n_coils (C), so a chunk’s T coils are always a contiguous, non-wrapping range within a single batch: a plain slice into self.smaps is a zero-copy view (unlike fancy/advanced indexing with an int array, which always allocates and copies).

_gram_op_calibless(data)[source]#: Compute the Gram operator without sensitivity maps.

_gram_op_calibless_device(data, img_d)[source]#

_gram_op_calibless_host(data, img_d)[source]#

_gram_op_raw_device(in_d, out_d, padded_array=None)[source]#: Apply the Toeplitz Gram operator on device to a (batched) image.

_gram_op_sense(data)[source]#: Compute the Gram operator with sensitivity maps.

_gram_op_sense_device(data, img_d)[source]#

_gram_op_sense_host(data, img_d)[source]#

_gram_op_toeplitz_raw(data) → ndarray[tuple[Any, ...], dtype[_ScalarT]][source]#

Return type:: ndarray[tuple[Any, …], dtype[_ScalarT]]

_make_plan_grad(**kwargs)[source]#

_safe_squeeze(arr)[source]#: Squeeze the first two dimensions of shape of the operator.

check_shape(*, image=None, ksp=None)[source]#

Validate the shapes of the image or k-space data against operator shapes.

Parameters:

image (NDArray, optional) – If passed, the shape of image data will be checked.
ksp (NDArray or object, optional) – If passed, the shape of the k-space data will be checked.

Raises:

ValueError – If the shape of the provided image does not match the expected operator shape, or if the number of k-space samples does not match the expected number of samples.

compute_density(method: Callable[[...], ndarray[tuple[Any, ...], dtype[_ScalarT]]] | bool | None | str | dict[str, Any] = None)[source]#

Compute the density compensation weights and set it.

Parameters:

method (str or callable or array or dict or bool) –

The method to use to compute the density compensation.

If a string, the method should be registered in the density registry.
If a callable, it should take the samples and the shape as input.
If a dict, it should have a key ‘name’, to determine which method to use. other items will be used as kwargs.
If an array, it should be of shape (Nsamples,) and will be used as is.
If True, the method pipe is chosen as default estimation method.

Notes

The “pipe” method is only available for the following backends: tensorflow, finufft, cufinufft, gpunufft, torchkbnufft-cpu and torchkbnufft-gpu.

compute_smaps(method: ndarray[tuple[Any, ...], dtype[_ScalarT]] | Callable[[...], ndarray[tuple[Any, ...], dtype[_ScalarT]]] | str | dict[str, Any] | None = None)[source]#

Compute the sensitivity maps and set it.

Parameters:: method (callable or dict or array) – The method to use to compute the sensitivity maps. If an array, it should be of shape (NCoils,XYZ) and will be used as is. If a dict, it should have a key ‘name’, to determine which method to use. other items will be used as kwargs. If a callable, it should take the samples and the shape as input. Note that this callable function should also hold the k-space data (use funtools.partial)

compute_toeplitz_kernel() → ndarray[tuple[Any, ...], dtype[_ScalarT]][source]#

Compute the Toeplitz kernel and set it.

Return type:: ndarray[tuple[Any, …], dtype[_ScalarT]]

property cpx_dtype[source]#: Return complex floating precision of the operator.

grad_traj_plan()[source]#: Context manager to enable gradient computation with respect to trajectory.

gram_op(data, img_d=None, toeplitz=True)[source]#

Compute the Gram operator of the NUFFT.

Parameters:

data (array) – Input data array.
img_d (array, optional) – Preallocated output array.
toeplitz (bool, default True) – If True, use the Toeplitz method to compute the Gram operator. If False, use the direct method.

Returns:

Array with the Gram operator applied.

Return type:

NDArray

Note

property img_full_shape: tuple[int, ...][source]#: Full image shape with batch and coil dimensions.

interfaces: dict[str, tuple[bool, type[FourierOperatorBase]]] = {'bart': (False, <class 'mrinufft.operators.interfaces.bart.MRIBartNUFFT'>), 'cartesian': (True, <class 'mrinufft.operators.cartesian.MRICartesianOperator'>), 'cufinufft': (True, <class 'mrinufft.operators.interfaces.cufinufft.MRICufiNUFFT'>), 'ducc0': (False, <class 'mrinufft.operators.interfaces.ducc0.MRIDUCC0'>), 'finufft': (True, <class 'mrinufft.operators.interfaces.finufft.MRIfinufft'>), 'gpunufft': (True, <class 'mrinufft.operators.interfaces.gpunufft.MRIGpuNUFFT'>), 'numpy': (True, <class 'mrinufft.operators.interfaces.nudft_numpy.MRInumpy'>), 'pynfft': (False, <class 'mrinufft.operators.interfaces.nfft.MRInfft'>), 'pynufft-cpu': (False, <class 'mrinufft.operators.interfaces.pynufft_cpu.MRIPynufft'>), 'sigpy': (True, <class 'mrinufft.operators.interfaces.sigpy.MRISigpyNUFFT'>), 'stacked': (True, <class 'mrinufft.operators.stacked.MRIStackedNUFFT'>), 'stacked-cufinufft': (True, <class 'mrinufft.operators.stacked.MRIStackedNUFFTGPU'>), 'tensorflow': (False, <class 'mrinufft.operators.interfaces.tfnufft.MRITensorflowNUFFT'>), 'torchkbnufft-cpu': (False, <class 'mrinufft.operators.interfaces.torchkbnufft.TorchKbNUFFTcpu'>), 'torchkbnufft-gpu': (False, <class 'mrinufft.operators.interfaces.torchkbnufft.TorchKbNUFFTgpu'>)}[source]#

property inv_norm_factor: floating[source]#

Reciprocal of norm_factor, cached for the operator lifetime.

shape is fixed at construction, so the normalization is constant; caching turns the per-call 1 / norm_factor (a property recompute plus a division) into a single stored reciprocal-multiply.

property ksp_full_shape: tuple[int, int, int][source]#: Full kspace shape with batch and coil dimensions.

property log: Logger[source]#: Logger for this operator instance, named after its defining module.

make_autograd(*, wrt_data: bool = True, wrt_traj: bool = False, paired_batch: bool = False) → MRINufftAutoGrad[source]#

Make a new Operator with autodiff support.

Parameters:

wrt_data (bool, optional) – If the gradient with respect to the data is computed, default is true
wrt_traj (bool, optional) – If the gradient with respect to the trajectory is computed, default is false
paired_batch (int, optional) – If provided, specifies batch size for varying data/smaps pairs. Default is None, which means no batching

Returns:

A NUFFT operator with autodiff capabilities.

Return type:

torch.nn.module

Raises:

ValueError – If autograd is not available.

make_deepinv_phy(*args, **kwargs) → DeepInvPhyNufft[source]#

Make a new DeepInv Physics with NUFFT operator.

Parameters:

wrt_data (bool, optional) – If the gradient with respect to the data is computed, default is true
wrt_traj (bool, optional) – If the gradient with respect to the trajectory is computed, default is false
paired_batch (int, optional) – If provided, specifies batch size for varying data/smaps pairs. Default is None, which means no batching
viewed_as_real (bool, optional) – If True, the DeepInverse physics wrapper accepts and returns real-view tensors with a final dimension of size 2 representing the real and imaginary parts. Default is False.

Returns:

A NUFFT operator with autodiff capabilities.

Return type:

torch.nn.module

Raises:

ValueError – If autograd is not available.

make_linops(*, cupy: bool = False)[source]#

Create a Scipy Linear Operator from the NUFFT operator.

We add a _nufft private attribute with the current operator.

Parameters:: cupy (bool, default False) – If True, create a Cupy Linear Operator

`__init__`
`adj_op`	Non Cartesian MRI adjoint operator.
`check_shape`	Validate the shapes of the image or k-space data against operator shapes.
`compute_density`	Compute the density compensation weights and set it.
`compute_smaps`	Compute the sensitivity maps and set it.
`compute_toeplitz_kernel`	Compute the Toeplitz kernel and set it.
`data_consistency`	Compute the data consistency estimation directly on gpu.
`get_lipschitz_cst`	Return the Lipschitz constant of the operator.
`grad_traj_plan`	Context manager to enable gradient computation with respect to trajectory.
`gram_op`	Compute the Gram operator of the NUFFT.
`make_autograd`	Make a new Operator with autodiff support.
`make_deepinv_phy`	Make a new DeepInv Physics with NUFFT operator.
`make_linops`	Create a Scipy Linear Operator from the NUFFT operator.
`op`	Non Cartesian MRI forward operator.
`pinv_solver`	Solves the linear system Ax = y.
`pipe`	Compute the density compensation weights for a given set of kspace locations.
`toggle_grad_traj`	Toggle gradient computation with respect to trajectory.
`update_samples`	Update the samples of the NUFFT operator.
`with_autograd`	Return a Fourier operator with autograd capabilities.
`with_off_resonance_correction`	Return a new operator with Off Resonnance Correction.

`autograd_available`
`available`
`backend`
`bsize_img`	Size in Bytes of the compute batch of images.
`bsize_ksp`	Size in Bytes of the compute batch of samples.
`cpx_dtype`	Return complex floating precision of the operator.
`density`	Density compensation of the operator.
`dtype`	Return floating precision of the operator.
`eps`	Return the underlying precision parameter.
`img_full_shape`	Full image shape with batch and coil dimensions.
`img_size`	Image size in bytes.
`interfaces`
`inv_norm_factor`	Reciprocal of `norm_factor`, cached for the operator lifetime.
`ksp_full_shape`	Full kspace shape with batch and coil dimensions.
`ksp_size`	k-space size in bytes.
`log`	Logger for this operator instance, named after its defining module.
`n_batchs`	Number of coils for the operator.
`n_coils`	Number of coils for the operator.
`n_samples`	Return the number of samples used by the operator.
`ndim`	Number of dimensions in image space of the operator.
`norm_factor`	Norm factor of the operator.
`samples`	Return the samples used by the operator.
`shape`	Shape of the image space of the operator.
`smaps`	Sensitivity maps of the operator.
`uses_density`	Return True if the operator uses density compensation.
`uses_sense`	Return True if the operator uses sensitivity maps.

MRICufiNUFFT#

Examples using mrinufft.operators.MRICufiNUFFT#

Examples using `mrinufft.operators.MRICufiNUFFT`#