API documentation

NUFFT class

class pynufft.nufft.__init__.NUFFT(device_indx=None, legacy=None)

NUFFT class

A super class of cpu and gpu NUFFT functions.

Note: NUFFT does not inherit NUFFT_cpu (deprecated) and NUFFT_hsa (deprecated).

__init__(device_indx=None, legacy=None)

Constructor.

Parameters:

None (Python NoneType)

Returns:

NUFFT: the pynufft.NUFFT instance

Return type:

NUFFT: the pynufft.NUFFT class

Example:

>>> from pynufft import NUFFT
>>> NufftObj = NUFFT()

or

>>> from pynufft import NUFFT, helper
>>> device = helper.device_list()[0]
>>> NufftObj = NUFFT(device) # for first acceleration device in the system
__weakref__

list of weak references to the object (if defined)

_adjoint_cpu(y)

Adjoint NUFFT on CPU

Parameters:

y – The input numpy array, with the size of (M,)

Type:

numpy array with the dtype of numpy.complex64

Returns:

x: The output numpy array, with the size of Nd or Nd

Return type:

numpy array with the dtype of numpy.complex64

_adjoint_device(gy)

Adjoint NUFFT on the heterogeneous device

Parameters:

gy – The input gpu array, with size=(M,)

Type:

reikna gpu array with dtype =numpy.complex64

Returns:

gx: The output gpu array, with size=Nd

Return type:

reikna gpu array with dtype =numpy.complex64

_adjoint_legacy(gy)

Adjoint NUFFT on the heterogeneous device

Parameters:

gy – The input gpu array, with size=(M,)

Type:

reikna gpu array with dtype =numpy.complex64

Returns:

gx: The output gpu array, with size=Nd

Return type:

reikna gpu array with dtype =numpy.complex64

_forward_cpu(x)

Forward NUFFT on CPU

Parameters:

x – The input numpy array, with the size of Nd

Type:

numpy array with the dtype of numpy.complex64

Returns:

y: The output numpy array, with the size of (M,)

Return type:

numpy array with the dtype of numpy.complex64

_forward_device(gx)

Forward NUFFT on the heterogeneous device

Parameters:

gx (reikna gpu array with dtype = numpy.complex64) – The input gpu array, with size = Nd

Returns:

gy: The output gpu array, with size = (M,)

Return type:

reikna gpu array with dtype = numpy.complex64

_forward_legacy(gx)

Forward NUFFT on the heterogeneous device

Parameters:

gx (reikna gpu array with dtype = numpy.complex64) – The input gpu array, with size = Nd

Returns:

gy: The output gpu array, with size = (M,)

Return type:

reikna gpu array with dtype = numpy.complex64

_init__cpu()

Constructor.

Parameters:

None (Python NoneType)

Returns:

NUFFT: the pynufft_hsa.NUFFT instance

Return type:

NUFFT: the pynufft_hsa.NUFFT class

Example:

>>> from pynufft import NUFFT
>>> NufftObj = NUFFT()
_init__device(device_indx=None)

Constructor.

Parameters:
  • API (string) – The API for the heterogeneous system. API=’cuda’ or API=’ocl’

  • platform_number (integer) – The number of the platform found by the API.

  • device_number (integer) – The number of the device found on the platform.

  • verbosity (integer) – Defines the verbosity level, default value is 0

Returns:

0

Return type:

int, float

Example:

>>> from pynufft import NUFFT
>>> NufftObj = NUFFT_device(API='cuda', platform_number=0,
                                 device_number=0, verbosity=0)
_k2xx_cpu(k)
Private: the inverse FFT and image cropping (which is the reverse of

_xx2k() method)

_k2xx_device(k)

Private: the inverse FFT and image cropping (which is the reverse of _xx2k() method)

_k2xx_one2one_cpu(k)
Private: the inverse FFT and image cropping

(which is the reverse of _xx2k() method)

_k2y2k_cpu(k)
Private: the integrated interpolation-gridding by the Sparse

Matrix-Vector Multiplication

_k2y_cpu(k)

Private: interpolation by the Sparse Matrix-Vector Multiplication

_k2y_device(k)

Private: interpolation by the Sparse Matrix-Vector Multiplication

_k2y_legacy(k)

Private: interpolation by the Sparse Matrix-Vector Multiplication

_offload_device()

self.offload():

Off-load NUFFT to the opencl or cuda device(s)

Parameters:
  • API (string) – define the device type, which can be ‘cuda’ or ‘ocl’

  • platform_number (int) – define which platform to be used. The default platform_number = 0.

  • device_number (int) – define which device to be used. The default device_number = 0.

Returns:

self: instance

_offload_legacy()

self.offload():

Off-load NUFFT to the opencl or cuda device(s)

Parameters:
  • API (string) – define the device type, which can be ‘cuda’ or ‘ocl’

  • platform_number (int) – define which platform to be used. The default platform_number = 0.

  • device_number (int) – define which device to be used. The default device_number = 0.

Returns:

self: instance

_plan_cpu(om, Nd, Kd, Jd, ft_axes=None)

Plan the NUFFT object with the geometry provided.

Parameters:
  • om (numpy.float array, matrix size = M * ndims) – The M off-grid locates in the frequency domain, which is normalized between [-pi, pi]

  • Nd (tuple, ndims integer elements.) –

    The matrix size of the equispaced image. Example: Nd=(256,256) for a 2D image;

    Nd = (128,128,128) for a 3D image

  • Kd (tuple, ndims integer elements.) –

    The matrix size of the oversampled frequency grid. Example: Kd=(512,512) for 2D image;

    Kd = (256,256,256) for a 3D image

  • Jd (tuple, ndims integer elements.) –

    The interpolator size. Example: Jd=(6,6) for 2D image;

    Jd = (6,6,6) for a 3D image

  • ft_axes (None, or tuple with optional integer elements.) – (Optional) The axes for Fourier transform. The default is all axes if ‘None’ is given.

  • batch – (Optional) Batch mode. If the batch is provided, the last appended axis is the number of identical NUFFT to be transformed. The default is ‘None’.

Returns:

0

Return type:

int, float

Variables:
  • Nd – initial value: Nd

  • Kd – initial value: Kd

  • Jd – initial value: Jd

  • ft_axes – initial value: None

Example:

>>> from pynufft import NUFFT
>>> NufftObj = NUFFT()
>>> NufftObj.plan(om, Nd, Kd, Jd)

or

>>> NufftObj.plan(om, Nd, Kd, Jd, ft_axes)
_plan_device(om, Nd, Kd, Jd, ft_axes=None, radix=None)

Design the multi-coil or single-coil memory reduced interpolator.

Parameters:
  • om (numpy.float array, matrix size = (M, ndims)) – The M off-grid locations in the frequency domain. Normalized between [-pi, pi]

  • Nd (tuple, ndims integer elements.) –

    The matrix size of equispaced image. Example: Nd=(256, 256) for a 2D image;

    Nd = (128, 128, 128) for a 3D image

  • Kd (tuple, ndims integer elements.) – The matrix size of the oversampled frequency grid. Example: Kd=(512,512) for 2D image; Kd = (256,256,256) for a 3D image

  • Jd (tuple, ndims integer elements.) – The interpolator size. Example: Jd=(6,6) for 2D image; Jd = (6,6,6) for a 3D image

  • ft_axes (tuple, selected axes to be transformed.) – The dimensions to be transformed by FFT. Example: ft_axes = (0, 1) for 2D, ft_axes = (0, 1, 2) for 3D; ft_axes = None for all dimensions.

  • radix – expert mode. If provided, the shape is Nd. The last axis is the number of parallel coils.

Returns:

0

Return type:

int, float

Example:

>>> import pynufft
>>> device=pynufft.helper.device_list()[0]
>>> NufftObj = pynufft.NUFFT(device)
>>> NufftObj.plan(om, Nd, Kd, Jd)
_plan_legacy(om, Nd, Kd, Jd, ft_axes=None)

Design the min-max interpolator.

Parameters:
  • om (numpy.float array, matrix size = M * ndims) – The M off-grid locations in the frequency domain. Normalized between [-pi, pi]

  • Nd (tuple, ndims integer elements.) – The matrix size of equispaced image. Example: Nd=(256,256) for a 2D image; Nd = (128,128,128) for a 3D image

  • Kd (tuple, ndims integer elements.) – The matrix size of the oversampled frequency grid. Example: Kd=(512,512) for 2D image; Kd = (256,256,256) for a 3D image

  • Jd (tuple, ndims integer elements.) – The interpolator size. Example: Jd=(6,6) for 2D image; Jd = (6,6,6) for a 3D image

Returns:

0

Return type:

int, float

Example:

>>> import pynufft
>>> NufftObj = pynufft.NUFFT_cpu()
>>> NufftObj.plan(om, Nd, Kd, Jd) 
_precompute_sp_cpu()
Private: Precompute adjoint (gridding) and Toepitz interpolation

matrix.

Parameters:

None (Python Nonetype)

Returns:

self: instance

_selfadjoint_cpu(x)

selfadjoint NUFFT on CPU

Parameters:

x – The input numpy array, with size=Nd

Type:

numpy array with dtype =numpy.complex64

Returns:

x: The output numpy array, with size=Nd

Return type:

numpy array with dtype =numpy.complex64

_selfadjoint_device(gx)

selfadjoint NUFFT on the heterogeneous device

Parameters:

gx (reikna gpu array with dtype =numpy.complex64) – The input gpu array, with size=Nd

Returns:

gx: The output gpu array, with size=Nd

Return type:

reikna gpu array with dtype =numpy.complex64

_selfadjoint_legacy(gx)

selfadjoint NUFFT on the heterogeneous device

Parameters:

gx (reikna gpu array with dtype =numpy.complex64) – The input gpu array, with size=Nd

Returns:

gx: The output gpu array, with size=Nd

Return type:

reikna gpu array with dtype =numpy.complex64

_solve_cpu(y, solver=None, *args, **kwargs)

Solve NUFFT. :param y: data, numpy.complex64. The shape = (M,) or (M, batch) :param solver: ‘cg’, ‘L1TVOLS’, ‘lsmr’, ‘lsqr’, ‘dc’, ‘bicg’,

‘bicgstab’, ‘cg’, ‘gmres’,’lgmres’

Parameters:

maxiter (int) – the number of iterations

Returns:

numpy array with size. The shape = Nd (‘L1TVOLS’) or Nd (‘lsmr’, ‘lsqr’, ‘dc’,’bicg’,’bicgstab’,’cg’, ‘gmres’,’lgmres’)

_solve_device(gy, solver=None, *args, **kwargs)

The solver of NUFFT

Parameters:
  • gy (reikna array, dtype = numpy.complex64) – data, reikna array, (M,) size

  • solver (string) – could be ‘cg’, ‘L1TVOLS’, ‘L1TVLAD’

  • maxiter (int) – the number of iterations

Returns:

reikna array with size Nd

_solve_legacy(gy, solver=None, *args, **kwargs)

The solver of NUFFT

Parameters:
  • gy (reikna array, dtype = numpy.complex64) – data, reikna array, (M,) size

  • solver (string) – could be ‘cg’, ‘L1TVOLS’, ‘L1TVLAD’

  • maxiter (int) – the number of iterations

Returns:

reikna array with size Nd

_vec2k_cpu(k_vec)

Sorting the vector to k-spectrum Kd array

_vec2y_cpu(k_vec)

gridding:

_x2xx_cpu(x)

Private: Scaling on CPU Inplace multiplication of self.x_Nd by the scaling factor self.sn.

_xx2k_cpu(xx)

Private: oversampled FFT on CPU

Firstly, zeroing the self.k_Kd array Second, copy self.x_Nd array to self.k_Kd array by cSelect Third, inplace FFT

_xx2k_device(xx)

Private: oversampled FFT on the heterogeneous device

First, zeroing the self.k_Kd array Second, copy self.x_Nd array to self.k_Kd array by cSelect Third, inplace FFT

_xx2k_one2one_cpu(xx)

Private: oversampled FFT on CPU

First, zeroing the self.k_Kd array Second, copy self.x_Nd array to self.k_Kd array by cSelect Third, inplace FFT

_xx2x_cpu(xx)

Private: rescaling, which is identical to the _x2xx() method

_y2k_cpu(y)

Private: gridding by the Sparse Matrix-Vector Multiplication

_y2k_device(y)

Private: gridding by the Sparse Matrix-Vector Multiplication Atomic_twosum together provide better accuracy than generic atomic_add. See: ocl_add and cuda_add code-strings in atomic_add(), inside the re_subroutine.py.

_y2k_legacy(y)

Private: gridding by the Sparse Matrix-Vector Multiplication

_y2vec_cpu(y)

regridding non-uniform data (unsorted vector)

adjoint(*args, **kwargs)

Adjoint NUFFT (host code)

Parameters:

y – The input numpy array, with the size of (M,)

Type:

numpy array with the dtype of numpy.complex64

Returns:

x: The output numpy array, with the size of Nd or Nd

Return type:

numpy array with the dtype of numpy.complex64

forward(*args, **kwargs)

Forward NUFFT (host code)

Parameters:

x – The input numpy array, with the size of Nd

Type:

numpy array with the dtype of numpy.complex64

Returns:

y: The output numpy array, with the size of (M,)

Return type:

numpy array with the dtype of numpy.complex64

plan(*args, **kwargs)

Plan the NUFFT object with the geometry provided.

Parameters:
  • om (numpy.float array, matrix size = M * ndims) – The M off-grid locates in the frequency domain, which is normalized between [-pi, pi]

  • Nd (tuple, ndims integer elements.) –

    The matrix size of the equispaced image. Example: Nd=(256,256) for a 2D image;

    Nd = (128,128,128) for a 3D image

  • Kd (tuple, ndims integer elements.) –

    The matrix size of the oversampled frequency grid. Example: Kd=(512,512) for 2D image;

    Kd = (256,256,256) for a 3D image

  • Jd (tuple, ndims integer elements.) –

    The interpolator size. Example: Jd=(6,6) for 2D image;

    Jd = (6,6,6) for a 3D image

  • ft_axes (None, or tuple with optional integer elements.) – (Optional) The axes for Fourier transform. The default is all axes if ‘None’ is given.

Returns:

0

Return type:

int, float

Variables:
  • Nd – initial value: Nd

  • Kd – initial value: Kd

  • Jd – initial value: Jd

  • ft_axes – initial value: None

Example:

>>> from pynufft import NUFFT
>>> NufftObj = NUFFT()
>>> NufftObj.plan(om, Nd, Kd, Jd)

or

>>> NufftObj.plan(om, Nd, Kd, Jd, ft_axes)
selfadjoint(*args, **kwargs)

selfadjoint NUFFT (host code)

Parameters:

x – The input numpy array, with size=Nd

Type:

numpy array with dtype =numpy.complex64

Returns:

x: The output numpy array, with size=Nd

Return type:

numpy array with dtype =numpy.complex64

solve(*args, **kwargs)

Solve NUFFT (host code) :param y: data, numpy.complex64. The shape = (M,) :param solver: ‘cg’, ‘L1TVOLS’, ‘lsmr’, ‘lsqr’, ‘dc’, ‘bicg’,

‘bicgstab’, ‘cg’, ‘gmres’,’lgmres’

Parameters:

maxiter (int) – the number of iterations

Returns:

numpy array with size Nd.

class pynufft.nufft.__init__.NUFFT_cupy

A tentative torch interface

Jd

initial value: ()

__init__()
__weakref__

list of weak references to the object (if defined)

debug

initial value: 0

plan(om, Nd, Kd, Jd)

Plan the NUFFT object with the geometry provided.

Parameters:
  • om (numpy.float array, matrix size = M * ndims) – The M off-grid locates in the frequency domain, which is normalized between [-pi, pi]

  • Nd (tuple, ndims integer elements.) –

    The matrix size of the equispaced image. Example: Nd=(256,256) for a 2D image;

    Nd = (128,128,128) for a 3D image

  • Kd (tuple, ndims integer elements.) –

    The matrix size of the oversampled frequency grid. Example: Kd=(512,512) for 2D image;

    Kd = (256,256,256) for a 3D image

  • Jd (tuple, ndims integer elements.) –

    The interpolator size. Example: Jd=(6,6) for 2D image;

    Jd = (6,6,6) for a 3D image

  • ft_axes (None, or tuple with optional integer elements.) – (Optional) The axes for Fourier transform. The default is all axes if ‘None’ is given.

  • batch – (Optional) Batch mode. If the batch is provided, the last appended axis is the number of identical NUFFT to be transformed. The default is ‘None’.

Returns:

0

Return type:

int, float

Variables:
  • Nd – initial value: Nd

  • Kd – initial value: Kd

  • Jd – initial value: Jd

  • ft_axes – initial value: None

Example:

>>> from pynufft import NUFFT
>>> NufftObj = NUFFT()
>>> NufftObj.plan(om, Nd, Kd, Jd)

or

>>> NufftObj.plan(om, Nd, Kd, Jd, ft_axes)
class pynufft.nufft.__init__.NUFFT_tf_eager

A tentative torch interface

Jd

initial value: ()

__init__()
__weakref__

list of weak references to the object (if defined)

debug

initial value: 0

plan(om, Nd, Kd, Jd)

Plan the NUFFT object with the geometry provided.

Parameters:
  • om (numpy.float array, matrix size = M * ndims) – The M off-grid locates in the frequency domain, which is normalized between [-pi, pi]

  • Nd (tuple, ndims integer elements.) –

    The matrix size of the equispaced image. Example: Nd=(256,256) for a 2D image;

    Nd = (128,128,128) for a 3D image

  • Kd (tuple, ndims integer elements.) –

    The matrix size of the oversampled frequency grid. Example: Kd=(512,512) for 2D image;

    Kd = (256,256,256) for a 3D image

  • Jd (tuple, ndims integer elements.) –

    The interpolator size. Example: Jd=(6,6) for 2D image;

    Jd = (6,6,6) for a 3D image

  • ft_axes (None, or tuple with optional integer elements.) – (Optional) The axes for Fourier transform. The default is all axes if ‘None’ is given.

  • batch – (Optional) Batch mode. If the batch is provided, the last appended axis is the number of identical NUFFT to be transformed. The default is ‘None’.

Returns:

0

Return type:

int, float

Variables:
  • Nd – initial value: Nd

  • Kd – initial value: Kd

  • Jd – initial value: Jd

  • ft_axes – initial value: None

Example:

>>> from pynufft import NUFFT
>>> NufftObj = NUFFT()
>>> NufftObj.plan(om, Nd, Kd, Jd)

or

>>> NufftObj.plan(om, Nd, Kd, Jd, ft_axes)
class pynufft.nufft.__init__.NUFFT_torch

A tentative torch interface

Jd

initial value: ()

__init__()
__weakref__

list of weak references to the object (if defined)

debug

initial value: 0

plan(om, Nd, Kd, Jd)

Plan the NUFFT object with the geometry provided.

Parameters:
  • om (numpy.float array, matrix size = M * ndims) – The M off-grid locates in the frequency domain, which is normalized between [-pi, pi]

  • Nd (tuple, ndims integer elements.) –

    The matrix size of the equispaced image. Example: Nd=(256,256) for a 2D image;

    Nd = (128,128,128) for a 3D image

  • Kd (tuple, ndims integer elements.) –

    The matrix size of the oversampled frequency grid. Example: Kd=(512,512) for 2D image;

    Kd = (256,256,256) for a 3D image

  • Jd (tuple, ndims integer elements.) –

    The interpolator size. Example: Jd=(6,6) for 2D image;

    Jd = (6,6,6) for a 3D image

  • ft_axes (None, or tuple with optional integer elements.) – (Optional) The axes for Fourier transform. The default is all axes if ‘None’ is given.

  • batch – (Optional) Batch mode. If the batch is provided, the last appended axis is the number of identical NUFFT to be transformed. The default is ‘None’.

Returns:

0

Return type:

int, float

Variables:
  • Nd – initial value: Nd

  • Kd – initial value: Kd

  • Jd – initial value: Jd

  • ft_axes – initial value: None

Example:

>>> from pynufft import NUFFT
>>> NufftObj = NUFFT()
>>> NufftObj.plan(om, Nd, Kd, Jd)

or

>>> NufftObj.plan(om, Nd, Kd, Jd, ft_axes)
sparse_mx_to_torch_sparse_tensor(sparse_mx)

Convert a scipy sparse matrix to a torch sparse tensor.

https://github.com/DSE-MSU/DeepRobust

pynufft.nufft.__init__.push_cuda_context(hsa_method)

Decorator: Push cude context to the top of the stack for current use Add @push_cuda_context before the methods of NUFFT_device()

CPU solvers

pynufft.linalg.solve_cpu.L1TVOLS(nufft, y, maxiter, rho)

L1-total variation regularized ordinary least square

pynufft.linalg.solve_cpu._create_kspace_sampling_density(nufft)

Compute k-space sampling density

pynufft.linalg.solve_cpu._pipe_density(nufft, maxiter)

Create the density function by iterative solution Generate pHp matrix

pynufft.linalg.solve_cpu.cDiff(x, d_indx)

Compute image gradient, which needs the results of indxmap_diff(Nd) :param x: The image array :param d_indx: The index of the shifted image :type x: numpy.float array, matrix size = Nd :type d_indx: int32 :returns: x_diff: Image gradient determined by d_indx :rtype: x_diff: numpy.complex64

pynufft.linalg.solve_cpu.solve(nufft, y, solver=None, *args, **kwargs)

Solve NUFFT. The current version supports solvers = ‘cg’ or ‘L1TVOLS’ or ‘L1TVLAD’.

Parameters:
  • nufft – NUFFT_cpu object

  • y – (M,) array, non-uniform data

Returns:

x: image

HSA solvers

pynufft.linalg.solve_hsa.L1TVLAD(nufft, gy, maxiter, rho)

(testing) L1-total variation regularized least absolute deviation

pynufft.linalg.solve_hsa.L1TVOLS(nufft, gy, maxiter, rho)

L1-total variation regularized ordinary least square

pynufft.linalg.solve_hsa._create_kspace_sampling_density(nufft)

(stable) Compute k-space sampling density from the nufft object

pynufft.linalg.solve_hsa._pipe_density(nufft, maxiter)

Private: create the density function in the data space by a iterative solution Pipe et al. 1999

pynufft.linalg.solve_hsa.cDiff(x, d_indx)

(stable) Compute image gradient Work with indxmap_diff(Nd). …

pynufft.linalg.solve_hsa.solve(nufft, gy, solver=None, maxiter=30, *args, **kwargs)

The solve function of NUFFT_hsa. The current version supports solvers = ‘cg’ or ‘L1TVOLS’.

Parameters:
  • nufft – NUFFT_hsa object

  • y (numpy.complex64 reikna array) – (M,) array, non-uniform data. If batch is provided, ‘cg’ and ‘L1TVOLS’ returns different image shape.

Returns:

x: Nd image. L1TVOLS always returns Nd. ‘cg’ returns Nd.

Return type:

x: reikna array, complex64.

NUDFT class

class pynufft.linalg.nudft_cpu.NUDFT

Class NUDFT

The non-uniform DFT operator

__init__()

Constructor.

Parameters:

None (Python NoneType)

Returns:

NUFFT: the pynufft.NUDFT instance

Return type:

NUFFT: the pynufft.NUFFT class

Example:

>>> from pynufft import NUDFT
>>> NufftObj = NUDFT()
__weakref__

list of weak references to the object (if defined)

debug

initial value: 0

Helper functions

class pynufft.src._helper.helper.ELL(elldata, ellcol)

ELL is slow on a single core CPU

__init__(elldata, ellcol)
__weakref__

list of weak references to the object (if defined)

pynufft.src._helper.helper.OMEGA_k(J, K, omd, Kd, dimid, dd, ft_flag)

Compute the index of k-space k_indx

pynufft.src._helper.helper.QR_process(om, N, J, K, sn)

1D QR method for generating min-max interpolator

Parameters:
  • om (numpy.float32) – non-Cartesian coordinate. shape = (M, dims)

  • N (int) – length

  • J (int) – size of interpolator

  • K (int) – size of oversampled grid

  • sn (numpy.float32 shape = (N,)) – scaling factor as a length-N vector

class pynufft.src._helper.helper.Tensor_sn(snd, radix)

Not implemented:

__init__(snd, radix)
__weakref__

list of weak references to the object (if defined)

pynufft.src._helper.helper.block_outer_prod(x1, x2)

Multiply x1 (J1 x M) and x2 (J2xM) and extend the dimensions to 3D (J1xJ2xM)

pynufft.src._helper.helper.block_outer_sum(x1, x2)

Update the new index after adding a new axis

pynufft.src._helper.helper.block_outer_sum0(x1, x2)

Multiply x1 (J1 x M) and x2 (J2xM) and extend the dimensions to 3D (J1xJ2xM)

pynufft.src._helper.helper.cat_snd(snd)
Parameters:

snd (tuple) – tuple of input 1D vectors

Returns:

tensor_sn: vector of concatenated scaling factor, shape = (numpy.sum(Nd), )

Return type:

tensor_sn: numpy.float32

pynufft.src._helper.helper.create_laplacian_kernel(nufft)

Create the multi-dimensional laplacian kernel in k-space

Parameters:

nufft – the NUFFT object

Returns:

uker: the multi-dimensional laplacian kernel in k-space (no fft shift used)

Return type:

numpy ndarray

pynufft.src._helper.helper.create_partialELL(ud, kd, Jd, M)
Parameters:
  • ud (tuple of numpy.complex64 arrays) – tuple of all 1D interpolators

  • kd (tuple of numpy.int32 arrays) – tuple of all 1D indices

  • Jd (tuple of int32) – tuple of interpolation sizes

  • M (int) – number of samples

Returns:

partialELL:

Return type:

partialELL: pELL instance

pynufft.src._helper.helper.crop_slice_ind(Nd)

(Deprecated in v.0.3.4) Return the “slice” of Nd size to index multi-dimensional array. “Slice” functions as the index of the array. This function is superseded by preindex_copy(), which avoid run-time indexing.

pynufft.src._helper.helper.device_list()

device_list() returns available devices for acceleration as a tuple. If no device is available, it returns an empty tuple.

pynufft.src._helper.helper.diagnose(verbosity=0)

Diagnosis function Find available devices when NUFFT.offload() fails.

pynufft.src._helper.helper.get_sn(J, K, N)

Compute the 1D scaling factor for the given J, K, N

Parameters:
  • J (int) – size of interpolator

  • K (int) – size of oversampled grid

  • N (int) – length

Returns:

sn: scaling factor as a length-N vector

Rtype sn:

numpy.float32 shape = (N,)

pynufft.src._helper.helper.indxmap_diff(Nd)

Preindixing for rapid image gradient. Diff(x) = x.flat[d_indx[0]] - x.flat Diff_t(x) = x.flat[dt_indx[0]] - x.flat

Parameters:

Nd (tuple with integers) – the dimension of the image

Returns:

d_indx: image gradient

Returns:

dt_indx: the transpose of the image gradient

Return type:

d_indx: lists with numpy ndarray

Return type:

dt_indx: lists with numpy ndarray

pynufft.src._helper.helper.kaiser_bessel_ft(u, J, alpha, kb_m, d)

Interpolation weight for given J/alpha/kb-m

pynufft.src._helper.helper.kronecker_scale(snd)

Compute the Kronecker product of the scaling factor.

Parameters:
  • snd (tuple of 1D numpy.array) – Tuple of 1D scaling factors

  • dd – Number of dimensions

Returns:

sn: The multi-dimensional Kronecker of the scaling factors

Return type:

Nd array

pynufft.src._helper.helper.nufft_T(N, J, K, alpha, beta)

Equation (29) and (26) in Fessler and Sutton 2003. Create the overlapping matrix CSSC (diagonal dominant matrix) of J points, then find the pseudo-inverse of CSSC

pynufft.src._helper.helper.nufft_alpha_kb_fit(N, J, K)

Find parameters alpha and beta for scaling factor st[‘sn’] The alpha is hardwired as [1,0,0…] when J = 1 (uniform scaling factor)

Parameters:
  • N (int) – size of image

  • J (int) – size of interpolator

  • K (int) – size of oversampled k-space

Returns:

alphas:

Returns:

beta:

Return type:

alphas: list of float

Return type:

beta:

pynufft.src._helper.helper.nufft_offset(om, J, K)

For every om point (outside regular grids), find the nearest central grid (from Kd dimension)

pynufft.src._helper.helper.nufft_r(om, N, J, K, alpha, beta)

Equation (30) of Fessler & Sutton’s paper

pynufft.src._helper.helper.nufft_scale1(N, K, alpha, beta, Nmid)

Calculate image space scaling factor

pynufft.src._helper.helper.outer_sum(xx, yy)

Superseded by numpy.add.outer() function

class pynufft.src._helper.helper.pELL(M, Jd, curr_sumJd, meshindex, kindx, udata)

class pELL: partial ELL format

__init__(M, Jd, curr_sumJd, meshindex, kindx, udata)

Constructor

Parameters:
  • M (int) – Number of samples

  • Jd (tuple of int) – Interpolator size

  • curr_sumJd (tuple of int) – Summation of Jd[0:d-1], for fast shift computing

  • meshindex (numpy.uint32, shape = (numpy.prod(Jd), dd)) – The tensor indices to all interpolation points

  • kindx (numpy.uint32, shape = (M, numpy.sum(Jd))) – Premixed k-indices to be combined

  • udata (numpy.complex64, shape = (M, numpy.sum(Jd))) – Premixed interpolation data values

Returns:

pELL: partial ELLpack class with the given values

Return type:

pELL: partial ELLpack class

__weakref__

list of weak references to the object (if defined)

pynufft.src._helper.helper.plan(om, Nd, Kd, Jd, ft_axes=None, format='CSR', radix=None)

Plan for the NUFFT object.

Parameters:
  • om (numpy.float) – Coordinate

  • Nd (tuple of int) – Image shape

  • Kd (tuple of int) – Oversampled grid shape

  • Jd (tuple of int) – Interpolator size

  • ft_axes (tuple of int) – Axes where FFT takes place

  • format (string, 'CSR' or 'pELL') – Output format of the interpolator. ‘CSR’: the precomputed Compressed Sparse Row (CSR) matrix. ‘pELL’: partial ELLPACK which precomputes the concatenated 1D interpolators.

Return st:

dictionary for NUFFT

pynufft.src._helper.helper.preindex_copy(Nd, Kd)

Building the array index for copying two arrays of sizes Nd and Kd. Only the front parts of the input/output arrays are copied. The oversize parts of the input array are truncated (if Nd > Kd), and the smaller size are zero-padded (if Nd < Kd)

Parameters:
  • Nd (tuple with integer elements) – tuple, the dimensions of array1

  • Kd (tuple with integer elements) – tuple, the dimensions of array2

Returns:

inlist: the index of the input array

Returns:

outlist: the index of the output array

Returns:

nelem: the length of the inlist and outlist (equal length)

Return type:

inlist: list with integer elements

Return type:

outlist: list with integer elements

Return type:

nelem: int

pynufft.src._helper.helper.rdx_kron(ud, kd, Jd, radix=None)

Radix-n Kronecker product of multi-dimensional array

Parameters:
  • ud (tuple of (M, Jd[d]) numpy.complex64 arrays) – 1D interpolators

  • kd (tuple of (M, Jd[d]) numpy.uint arrays) – 1D indices to interpolators

  • Jd (tuple of int) – 1D interpolator sizes

  • radix (int) – radix of Kronecker product

  • kk (tuple of (M, Jd[d]) numpy.uint arrays) – 1D indices to interpolators

  • JJ (tuple of int) – 1D interpolator sizes

Returns:

uu: 1D interpolators

pynufft.src._helper.helper.strides_divide_itemsize(Nd)

strides_divide_itemsize function computes the step_size (strides/itemsize) along different axes, and its inverse as float32. For fast GPU computing, preindexing allows for fast Hadamard product and copy. However preindexing costs some memory. strides_divide_itemsize aims to replace preindexing by run-time calculation of the index, given the invNd_elements.

Parameters:

Nd (tuple of int) – Input shape

Returns:

Nd_elements: strides/itemsize of the Nd.

Returns:

invNd_elements: (float32)(1/Nd_elements). Division on GPU is slow but multiply is fast. Thus we can precompute the inverse and then multiply the inverse on GPU.

Return type:

Nd_elements: tuple of int

Return type:

invNd_elements: tuple of float32

See also

pynufft.NUFFT_hsa

Metaprogramming subroutines (using reikna, pyopencl, pycuda)

pynufft.src.re_subroutine.atomic_add(API)

Return the atomic_add for the given API. Overcome the missing atomic_add_float for OpenCL-1.2. Note: will be checked if OpenCL 2.0 provided by all GPU vendors.

pynufft.src.re_subroutine.cAddScalar()

Return the kernel source for cAddScalar.

pynufft.src.re_subroutine.cAddVec()

Return the kernel source for cAddVec.

pynufft.src.re_subroutine.cAnisoShrink()

Return the kernel source of cAnisoShrink

pynufft.src.re_subroutine.cCopy()

Return the kernel source for cCopy

pynufft.src.re_subroutine.cDiff()

Return the kernel source of cDiff.

pynufft.src.re_subroutine.cHadamard()

Return the Hadamard operations related kernel sources.

pynufft.src.re_subroutine.cHypot()

Return the kernel code for hypot, which computes the sqrt(x*x + y*y) without intermediate overflow.

pynufft.src.re_subroutine.cMultiplyConjVec()

Return the kernel source of cMultiplyConjVec.

pynufft.src.re_subroutine.cMultiplyConjVecInplace()

Return the kernel source of cMultiplyConjVecInplace

pynufft.src.re_subroutine.cMultiplyRealInplace()

Return the kernel source of cMultiplyRealInplace.

pynufft.src.re_subroutine.cMultiplyScalar()

Return the kernel source for cMultiplyScalar.

pynufft.src.re_subroutine.cMultiplyVec()

Return the kernel source of cMultiplyVec

pynufft.src.re_subroutine.cMultiplyVecInplace()

Return the kernel source of cMultiplyVecInplace.

pynufft.src.re_subroutine.cRealShrink()

Return the kernel source of xAnisoShrink

pynufft.src.re_subroutine.cSelect()

Return the kernel source of cSelect.

pynufft.src.re_subroutine.cSpmv()

Return the kernel sources for cSpmv related operations, providing cCSR_spmv_vector and cpELL_spmv_mCoil.

pynufft.src.re_subroutine.cSpmvh()

Return the cSpmvh related kernel source. Only pELL_spmvh_mCoil is provided for Spmvh. NUFFT_hsa_legacy reuses the cCSR_spmv() function, which doubles the storage.

pynufft.src.re_subroutine.cSqrt()

Return the kernel source of cSqrt.

pynufft.src.re_subroutine.cTensorCopy()

Return the kernel source for cTensorCopy.

pynufft.src.re_subroutine.cTensorMultiply()

Return the kernel source for cTensorMultiply

pynufft.src.re_subroutine.create_kernel_sets(API)

Create the kernel from the kernel sets. Note that in some tests (Benoit’s and my tests) CUDA shows some degraded accuracy. This loss of accuracy was due to undefined shared memory behavior, which I don’t fully understand. This has been fixed in 2019.2.0 as the operations are moved to global memory.