mrpro.operators.FastFourierOp

class mrpro.operators.FastFourierOp(dim: Sequence[int] = (-3, -2, -1), recon_matrix: SpatialDimension[int] | Sequence[int] | None = None, encoding_matrix: SpatialDimension[int] | Sequence[int] | None = None)[source]

Bases: LinearOperator

Fast Fourier operator class.

Applies a Fast Fourier Transformation along selected dimensions with cropping/zero-padding along these selected dimensions

The transformation is done with ‘ortho’ normalization, i.e. the normalization constant is split between forward and adjoint [FFT].

Remark regarding the fftshift/ifftshift: fftshift shifts the zero-frequency point to the center of the data, ifftshift undoes this operation. The input to both forward and ajoint assumes that the zero-frequency is in the center of the data. Torch.fft.fftn and torch.fft.ifftn expect the zero-frequency to be the first entry in the tensor. Therefore for forward and ajoint first ifftshift needs to be applied, then fftn or ifftn and then ifftshift.

References

__init__(dim: Sequence[int] = (-3, -2, -1), recon_matrix: SpatialDimension[int] | Sequence[int] | None = None, encoding_matrix: SpatialDimension[int] | Sequence[int] | None = None) → None[source]

Initialize a Fast Fourier Operator.

If both recon_matrix and encoding_matrix are set, the operator will perform padding/cropping before and after the transforms to match the shape in image space (recon_matrix) and k-shape (encoding_matrix). If both are set to None, no padding or cropping will be performed. If these are SpatialDimension, the transform dimensions must be within the last three dimensions, typically corresponding to the (k2,k1,k0) and (z,y,x) axes of KData and IData, respectively.

Parameters:

  • dim – dim along which FFT and IFFT are applied, by default last three dimensions (-3, -2, -1), as these correspond to k2, k1, and k0 of KData.

  • encoding_matrix – shape of encoded k-data along the axes in dim. Must be set if recon_matrix is set. If encoding_matrix and recon_matrix are None, no padding or cropping will be performed. If all values in dim are -3, -2 or -1, this can also be a SpatialDimension describing the k-space shape in all 3 dimensions (k2, k1, k0), but only values in the dimensions in dim will be used. Otherwise, it should be a Sequence of the same length as dim.

  • recon_matrix – shape of reconstructed image data. Must be set if encoding_matrix is set. If encoding_matrix and recon_matrix are None, no padding or cropping will be performed. If all values in dim are -3, -2 or -1, this can also be a SpatialDimension describing the image-space shape in all 3 dimensions (z,y,x), but only values in the dimensions in dim will be used. Otherwise, it should be a Sequence of the same length as dim.

adjoint(y: Tensor) → tuple[Tensor][source]

IFFT from k-space to image space.

Parameters:

y – k-space data on Cartesian grid

Return type:

IFFT of y

forward(x: Tensor) → tuple[Tensor][source]

FFT from image space to k-space.

Parameters:

x – image data on Cartesian grid

Return type:

FFT of x

operator_norm(initial_value: Tensor, dim: Sequence[int] | None, max_iterations: int = 20, relative_tolerance: float = 0.0001, absolute_tolerance: float = 1e-05, callback: Callable[[Tensor], None] | None = None) → Tensor

Power iteration for computing the operator norm of the linear operator.

Parameters:

  • initial_value – initial value to start the iteration; if the initial value contains a zero-vector for one of the considered problems, the function throws an value error.

  • dim – the dimensions of the tensors on which the operator operates. For example, for a matrix-vector multiplication example, a batched matrix tensor with shape (4,30,80,160), input tensors of shape (4,30,160) to be multiplied, and dim = None, it is understood that the matrix representation of the operator corresponds to a block diagonal operator (with 4*30 matrices) and thus the algorithm returns a tensor of shape (1,1,1) containing one single value. In contrast, if for example, dim=(-1,), the algorithm computes a batched operator norm and returns a tensor of shape (4,30,1) corresponding to the operator norms of the respective matrices in the diagonal of the block-diagonal operator (if considered in matrix representation). In any case, the output of the algorithm has the same number of dimensions as the elements of the domain of the considered operator (whose dimensionality is implicitly defined by choosing dim), such that the pointwise multiplication of the operator norm and elements of the domain (to be for example used in a Landweber iteration) is well-defined.

  • max_iterations – maximum number of iterations

  • relative_tolerance – absolute tolerance for the change of the operator-norm at each iteration; if set to zero, the maximal number of iterations is the only stopping criterion used to stop the power iteration

  • absolute_tolerance – absolute tolerance for the change of the operator-norm at each iteration; if set to zero, the maximal number of iterations is the only stopping criterion used to stop the power iteration

  • callback – user-provided function to be called at each iteration

Return type:

an estimaton of the operator norm

property H: LinearOperator

Adjoint operator.

property gram: LinearOperator

Gram operator.

For a LinearOperator \(A\), the self-adjoint Gram operator is defined as \(A^H A\).

Note: This is a default implementation that can be overwritten by subclasses for more efficient implementations.