diff --git a/dev/abstract/index.html b/dev/abstract/index.html index 9df5653..9e4fa5b 100644 --- a/dev/abstract/index.html +++ b/dev/abstract/index.html @@ -9,4 +9,4 @@ ... end

We next outline all of the aforementioned functions and describe their behavior:

    size_in(p)

Size of the input array for an NFFT operation. The returned tuple has D entries. Note that this will be the output array for an adjoint NFFT.

    size_out(p)

Size of the output array for an NFFT operation. The returned tuple has R entries. Note that this will be the input array for an adjoint NFFT.

    mul!(fHat, p, f) -> fHat

Inplace NFFT transforming the D dimensional array f to the R dimensional array fHat. The transformation is applied along D-R+1 dimensions specified in the plan p. Both f and fHat must be complex arrays of element type Complex{T}.

    mul!(f, p::Adjoint{Complex{T},<:AbstractNFFTPlan{T}}, fHat) -> f

Inplace adjoint NFFT transforming the R dimensional array fHat to the D dimensional array f. The transformation is applied along D-R+1 dimensions specified in the plan p. Both f and fHat must be complex arrays of element type Complex{T}.

    nodes!(p, k)

Exchange the nodes k in the plan p and return the plan. The implementation of this function is optional.

Plan Interface

The constructor for a plan also has a defined interface. It should be implemented in this way:

function MyNFFTPlan(k::Matrix{T}, N::NTuple{D,Int}; kwargs...) where {T,D}
   ...
-end

All parameters are put into keyword arguments that have to match as well. We describe the keyword arguments in more detail in the overview page. Using the same plan interface allows to load several NFFT libraries simultaneously and exchange the constructor dynamically by storing the constructor in a function object. This is how the unit tests of NFFT.jl run.

Additionally, to the type-specific constructor one can provide the factory

plan_nfft(Q::Type, k::Matrix{T}, N::NTuple{D,Int}; kargs...) where {D}

where Q is the Array type, e.g. Array. The reason to require the array type is, that this allows for GPU implementations, which would use for instance CuArray here.

The package AbstractNFFTs provides a convenient constructor

plan_nfft(k::Matrix{T}, N::NTuple{D,Int}; kargs...) where {D}

defaulting to the Array type.

Note

Different packages implementing plan_nfft will conflict if the same Q is implemented. In case of NFFT.jl and CuNFFT.jl there is no conflict since the array type is different.

Derived Interface

Based on the core low-level interface that an AbstractNFFTPlan needs to provide, the package AbstractNFFT.jl also provides high-level functions like *, nfft, and nfft_adjoint, which internally use the low-level interface. Thus, the implementation of high-level function is shared among all AbstractNFFT.jl implementations.

+end

All parameters are put into keyword arguments that have to match as well. We describe the keyword arguments in more detail in the overview page. Using the same plan interface allows to load several NFFT libraries simultaneously and exchange the constructor dynamically by storing the constructor in a function object. This is how the unit tests of NFFT.jl run.

Additionally, to the type-specific constructor one can provide the factory

plan_nfft(Q::Type, k::Matrix{T}, N::NTuple{D,Int}; kargs...) where {D}

where Q is the Array type, e.g. Array. The reason to require the array type is, that this allows for GPU implementations, which would use for instance CuArray here.

The package AbstractNFFTs provides a convenient constructor

plan_nfft(k::Matrix{T}, N::NTuple{D,Int}; kargs...) where {D}

defaulting to the Array type.

Note

Different packages implementing plan_nfft will conflict if the same Q is implemented. In case of NFFT.jl and CuNFFT.jl there is no conflict since the array type is different.

Derived Interface

Based on the core low-level interface that an AbstractNFFTPlan needs to provide, the package AbstractNFFT.jl also provides high-level functions like *, nfft, and nfft_adjoint, which internally use the low-level interface. Thus, the implementation of high-level function is shared among all AbstractNFFT.jl implementations.

diff --git a/dev/api/index.html b/dev/api/index.html index c1d4b94..b1c2606 100644 --- a/dev/api/index.html +++ b/dev/api/index.html @@ -1,5 +1,5 @@ -API · NFFT
Edit on GitHub

API

In the following, you find the documentation of some, and hopefully soon all, exported functions of the NFFT.jl package family:

AbstractNFFTs.jl

NFFT.jl

NFFTTools.jl

AbstractNFFTs.AbstractComplexFTPlan — Type

AbstractComplexFTPlan{T,D,R}

Abstract type for either an NFFT plan or an NNFFT plan.

  • T is the element type (Float32/Float64)
  • D is the number of dimensions of the input array.
  • R is the number of dimensions of the output array.
AbstractNFFTs.AbstractFTPlan — Type

AbstractFTPlan{T,D,R}

Abstract type for any NFFT-like plan (NFFT, NFFT, NFCT, NFST).

  • T is the element type (Float32/Float64)
  • D is the number of dimensions of the input array.
  • R is the number of dimensions of the output array.
AbstractNFFTs.AbstractNFCTPlan — Type

AbstractNFCTPlan{T,D,R}

Abstract type for an NFCT plan.

  • T is the element type (Float32/Float64)
  • D is the number of dimensions of the input array.
  • R is the number of dimensions of the output array.
AbstractNFFTs.AbstractNFFTPlan — Type

AbstractNFFTPlan{T,D,R}

Abstract type for an NFFT plan.

  • T is the element type (Float32/Float64)
  • D is the number of dimensions of the input array.
  • R is the number of dimensions of the output array.
AbstractNFFTs.AbstractNFSTPlan — Type

AbstractNFSTPlan{T,D,R}

Abstract type for an NFST plan.

  • T is the element type (Float32/Float64)
  • D is the number of dimensions of the input array.
  • R is the number of dimensions of the output array.
AbstractNFFTs.AbstractNNFFTPlan — Type

AbstractNNFFTPlan{T,D,R}

Abstract type for an NNFFT plan.

  • T is the element type (Float32/Float64)
  • D is the number of dimensions of the input array.
  • R is the number of dimensions of the output array.
AbstractNFFTs.AbstractRealFTPlan — Type

AbstractRealFTPlan{T,D,R}

Abstract type for either an NFCT plan or an NFST plan.

  • T is the element type (Float32/Float64)
  • D is the number of dimensions of the input array.
  • R is the number of dimensions of the output array.
AbstractNFFTs.accuracyParams — Method

accuracyParams(; [m, σ, reltol]) -> m, σ, reltol

Calculate accuracy parameters m, σ, reltol based on either

  • reltol

or

  • m, σ

TODO: Right now, the oversampling parameter is not taken into account, i.e. σ=2.0 is assumed

AbstractNFFTs.nfct — Method

nfft(k, f, rest...; kwargs...)

calculates the nfft of the array f for the nodes contained in the matrix k The output is a vector of length M=size(nodes,2)

AbstractNFFTs.nfct_transpose — Method

nfft_adjoint(k, N, fHat, rest...; kwargs...)

calculates the adjoint nfft of the vector fHat for the nodes contained in the matrix k. The output is an array of size N

AbstractNFFTs.nfft — Method

nfft(k, f, rest...; kwargs...)

calculates the nfft of the array f for the nodes contained in the matrix k The output is a vector of length M=size(nodes,2)

AbstractNFFTs.nfft_adjoint — Method

nfft_adjoint(k, N, fHat, rest...; kwargs...)

calculates the adjoint nfft of the vector fHat for the nodes contained in the matrix k. The output is an array of size N

AbstractNFFTs.nfst — Method

nfft(k, f, rest...; kwargs...)

calculates the nfft of the array f for the nodes contained in the matrix k The output is a vector of length M=size(nodes,2)

AbstractNFFTs.nfst_transpose — Method

nfft_adjoint(k, N, fHat, rest...; kwargs...)

calculates the adjoint nfft of the vector fHat for the nodes contained in the matrix k. The output is an array of size N

AbstractNFFTs.nodes! — Method
nodes!(p, k) -> p

Change nodes k in the plan p operation and return the plan.

AbstractNFFTs.size_in — Method
size_in(p)

Size of the input array for the plan p (NFFT/NFCT/NFST/NNFFT). The returned tuple has R entries. Note that this will be the output size for the transposed / adjoint operator.

AbstractNFFTs.size_out — Method
size_out(p)

Size of the output array for the plan p (NFFT/NFCT/NFST/NNFFT). The returned tuple has R entries. Note that this will be the input size for the transposed / adjoint operator.

Base.:* — Method
    *(p, f) -> fHat

For a non-directional D dimensional plan p this calculates the NFFT/NNFFT of a D dimensional array f of size N. fHat is a vector of length M. (M and N are defined in the plan creation)

For a directional D dimensional plan p both f and fHat are D dimensional arrays, and the dimension specified in the plan creation is affected.

Base.:* — Method
    *(p, f) -> fHat

For a non-directional D dimensional plan p this calculates the NFCT/NFST of a D dimensional array f of size N. fHat is a vector of length M. (M and N are defined in the plan creation)

For a directional D dimensional plan p both f and fHat are D dimensional arrays, and the dimension specified in the plan creation is affected.

LinearAlgebra.mul! — Method
mul!(fHat, p, f) -> fHat

Inplace NFFT/NFCT/NFST/NNFFT transforming the D dimensional array f to the R dimensional array fHat. The transformation is applied along D-R+1 dimensions specified in the plan p.

LinearAlgebra.mul! — Method
mul!(f, p, fHat) -> f

Inplace adjoint NFFT/NNFFT transforming the R dimensional array fHat to the D dimensional array f. The transformation is applied along D-R+1 dimensions specified in the plan p.

LinearAlgebra.mul! — Method
mul!(f, p, fHat) -> f

Inplace transposed NFCT/NFST transforming the R dimensional array fHat to the D dimensional array f. The transformation is applied along D-R+1 dimensions specified in the plan p.

AbstractNFFTs.plan_nfft — Method
    plan_nfft(k::Matrix{T}, N::NTuple{D,Int}, rest...;  kargs...)

compute a plan for the NFFT of a size-N array at the nodes contained in k.

source
NFFT.precomputeLinInterp — Method

Precompute the look up table for the window function φ.

Remarks:

  • Only the positive half is computed
  • The window is computed for the interval [0, (m)/Ñ]. The reason for the +2 is that we do evaluate the window function outside its interval, since x does not necessary match the sampling points
  • The window has K+1 entries and during the index calculation we multiply with the factor K/(m).
  • It is very important that K/(m) is an integer since our index calculation exploits this fact. We therefore always use Int(K/(m))instead of K÷(m) since this gives an error while the later variant would silently error.
source
NFFTTools.calculateToeplitzKernel! — Method
calculateToeplitzKernel!(f::Array{Complex{T}}, p::AbstractNFFTPlan, tr::Matrix{T}, fftplan)

Calculate the kernel for an implementation of the Gram matrix that utilizes its Toeplitz structure and writes it in-place in f, which has to be twice the size of the desired image matrix, as the Toeplitz trick requires an oversampling factor of 2 (cf. Wajer, F. T. A. W., and K. P. Pruessmann. "Major speedup of reconstruction for sensitivity encoding with arbitrary trajectories." Proc. Intl. Soc. Mag. Res. Med. 2001.). The type of the kernel f has to be Complex{T}, i.e. the complex of the k-space trajectory's type; for speed and memory efficiecy, call this function with Float32.(tr), and set the type of f accordingly.

Required Arguments

  • f::Array{T}: Array in which the kernel will be written.
  • p::AbstractNFFTPlan: NFFTPlan with the same dimensions as tr, which will be overwritten in place.
  • tr::Matrix{T}: non-Cartesian k-space trajectory in units revolutions/voxel, i.e. tr[i] ∈ [-0.5, 0.5] ∀ i. The matrix has the size 2 x Nsamples for 2D imaging with a trajectory length Nsamples, and 3 x Nsamples for 3D imaging.
  • fftplan: plan for the final FFT of the kernel from image to k-space. Therefore, it has to have twice the size of the original image. Calculate, e.g., with fftplan = plan_fft(zeros(Complex{T}, 2 .* shape); flags=FFTW.MEASURE), where shape is the size of the reconstructed image.

Examples

julia> using FFTW
+API · NFFT
Edit on GitHub
API
In the following, you find the documentation of some, and hopefully soon all, exported functions of the NFFT.jl package family:
AbstractNFFTs.jl
NFFT.jl
NFFTTools.jl
AbstractNFFTs.AbstractComplexFTPlan — Type
AbstractComplexFTPlan{T,D,R}
Abstract type for either an NFFT plan or an NNFFT plan.
  • T is the element type (Float32/Float64)
  • D is the number of dimensions of the input array.
  • R is the number of dimensions of the output array.
AbstractNFFTs.AbstractFTPlan — Type
AbstractFTPlan{T,D,R}
Abstract type for any NFFT-like plan (NFFT, NFFT, NFCT, NFST).
  • T is the element type (Float32/Float64)
  • D is the number of dimensions of the input array.
  • R is the number of dimensions of the output array.
AbstractNFFTs.AbstractNFCTPlan — Type
AbstractNFCTPlan{T,D,R}
Abstract type for an NFCT plan.
  • T is the element type (Float32/Float64)
  • D is the number of dimensions of the input array.
  • R is the number of dimensions of the output array.
AbstractNFFTs.AbstractNFFTPlan — Type
AbstractNFFTPlan{T,D,R}
Abstract type for an NFFT plan.
  • T is the element type (Float32/Float64)
  • D is the number of dimensions of the input array.
  • R is the number of dimensions of the output array.
AbstractNFFTs.AbstractNFSTPlan — Type
AbstractNFSTPlan{T,D,R}
Abstract type for an NFST plan.
  • T is the element type (Float32/Float64)
  • D is the number of dimensions of the input array.
  • R is the number of dimensions of the output array.
AbstractNFFTs.AbstractNNFFTPlan — Type
AbstractNNFFTPlan{T,D,R}
Abstract type for an NNFFT plan.
  • T is the element type (Float32/Float64)
  • D is the number of dimensions of the input array.
  • R is the number of dimensions of the output array.
AbstractNFFTs.AbstractRealFTPlan — Type
AbstractRealFTPlan{T,D,R}
Abstract type for either an NFCT plan or an NFST plan.
  • T is the element type (Float32/Float64)
  • D is the number of dimensions of the input array.
  • R is the number of dimensions of the output array.
AbstractNFFTs.accuracyParams — Method
accuracyParams(; [m, σ, reltol]) -> m, σ, reltol
Calculate accuracy parameters m, σ, reltol based on either
  • reltol
or
  • m, σ
TODO: Right now, the oversampling parameter is not taken into account, i.e. σ=2.0 is assumed
AbstractNFFTs.nfct — Method
nfft(k, f, rest...; kwargs...)
calculates the nfft of the array f for the nodes contained in the matrix k The output is a vector of length M=size(nodes,2)
AbstractNFFTs.nfct_transpose — Method
nfft_adjoint(k, N, fHat, rest...; kwargs...)
calculates the adjoint nfft of the vector fHat for the nodes contained in the matrix k. The output is an array of size N
AbstractNFFTs.nfft — Method
nfft(k, f, rest...; kwargs...)
calculates the nfft of the array f for the nodes contained in the matrix k The output is a vector of length M=size(nodes,2)
AbstractNFFTs.nfft_adjoint — Method
nfft_adjoint(k, N, fHat, rest...; kwargs...)
calculates the adjoint nfft of the vector fHat for the nodes contained in the matrix k. The output is an array of size N
AbstractNFFTs.nfst — Method
nfft(k, f, rest...; kwargs...)
calculates the nfft of the array f for the nodes contained in the matrix k The output is a vector of length M=size(nodes,2)
AbstractNFFTs.nfst_transpose — Method
nfft_adjoint(k, N, fHat, rest...; kwargs...)
calculates the adjoint nfft of the vector fHat for the nodes contained in the matrix k. The output is an array of size N
AbstractNFFTs.nodes! — Method
nodes!(p, k) -> p
Change nodes k in the plan p operation and return the plan.
AbstractNFFTs.size_in — Method
size_in(p)
Size of the input array for the plan p (NFFT/NFCT/NFST/NNFFT).  The returned tuple has R entries.  Note that this will be the output size for the transposed / adjoint operator.
AbstractNFFTs.size_out — Method
size_out(p)
Size of the output array for the plan p (NFFT/NFCT/NFST/NNFFT).  The returned tuple has R entries.  Note that this will be the input size for the transposed / adjoint operator.
Base.:* — Method
    *(p, f) -> fHat
For a non-directional D dimensional plan p this calculates the NFFT/NNFFT of a D dimensional array f of size N. fHat is a vector of length M. (M and N are defined in the plan creation)
For a directional D dimensional plan p both f and fHat are D dimensional arrays, and the dimension specified in the plan creation is affected.
Base.:* — Method
    *(p, f) -> fHat
For a non-directional D dimensional plan p this calculates the NFCT/NFST of a D dimensional array f of size N. fHat is a vector of length M. (M and N are defined in the plan creation)
For a directional D dimensional plan p both f and fHat are D dimensional arrays, and the dimension specified in the plan creation is affected.
LinearAlgebra.mul! — Method
mul!(fHat, p, f) -> fHat
Inplace NFFT/NFCT/NFST/NNFFT transforming the D dimensional array f to the R dimensional array fHat. The transformation is applied along D-R+1 dimensions specified in the plan p.
LinearAlgebra.mul! — Method
mul!(f, p, fHat) -> f
Inplace adjoint NFFT/NNFFT transforming the R dimensional array fHat to the D dimensional array f. The transformation is applied along D-R+1 dimensions specified in the plan p.
LinearAlgebra.mul! — Method
mul!(f, p, fHat) -> f
Inplace transposed NFCT/NFST transforming the R dimensional array fHat to the D dimensional array f. The transformation is applied along D-R+1 dimensions specified in the plan p.
AbstractNFFTs.plan_nfft — Method
    plan_nfft(k::Matrix{T}, N::NTuple{D,Int}, rest...;  kargs...)
compute a plan for the NFFT of a size-N array at the nodes contained in k.
source
NFFT.precomputeLinInterp — Method
Precompute the look up table for the window function φ.
Remarks:
  • Only the positive half is computed
  • The window is computed for the interval [0, (m)/Ñ]. The reason for the +2 is that we do evaluate the window function outside its interval, since x does not necessary match the sampling points
  • The window has K+1 entries and during the index calculation we multiply with the factor K/(m).
  • It is very important that K/(m) is an integer since our index calculation exploits this fact. We therefore always use Int(K/(m))instead of K÷(m) since this gives an error while the later variant would silently error.
source
NFFTTools.calculateToeplitzKernel! — Method
calculateToeplitzKernel!(f::Array{Complex{T}}, p::AbstractNFFTPlan, tr::Matrix{T}, fftplan)
Calculate the kernel for an implementation of the Gram matrix that utilizes its Toeplitz structure and writes it in-place in f, which has to be twice the size of the desired image matrix, as the Toeplitz trick requires an oversampling factor of 2 (cf. Wajer, F. T. A. W., and K. P. Pruessmann. "Major speedup of reconstruction for sensitivity encoding with arbitrary trajectories." Proc. Intl. Soc. Mag. Res. Med. 2001.). The type of the kernel f has to be Complex{T}, i.e. the complex of the k-space trajectory's type; for speed and memory efficiecy, call this function with Float32.(tr), and set the type of f accordingly.
Required Arguments
  • f::Array{T}: Array in which the kernel will be written.
  • p::AbstractNFFTPlan: NFFTPlan with the same dimensions as tr, which will be overwritten in place.
  • tr::Matrix{T}: non-Cartesian k-space trajectory in units revolutions/voxel, i.e. tr[i] ∈ [-0.5, 0.5] ∀ i. The matrix has the size 2 x Nsamples for 2D imaging with a trajectory length Nsamples, and 3 x Nsamples for 3D imaging.
  • fftplan: plan for the final FFT of the kernel from image to k-space. Therefore, it has to have twice the size of the original image. Calculate, e.g., with fftplan = plan_fft(zeros(Complex{T}, 2 .* shape); flags=FFTW.MEASURE), where shape is the size of the reconstructed image.
Examples
julia> using FFTW
 
 julia> Nx = 32;
 
@@ -108,4 +108,4 @@
   1685.44-1243.29im   1911.16-1157.72im     -991.639-42.8214im
   1054.11-1282.41im   66.9358-588.991im  …  -952.238+1026.35im
  -417.276-273.367im  -946.698+1971.77im     -890.339-882.05im
-

+
diff --git a/dev/background/index.html b/dev/background/index.html index 499966b..46b9046 100644 --- a/dev/background/index.html +++ b/dev/background/index.html @@ -3,4 +3,4 @@ \hat{\bm{f}} &:= \left( \hat{f}_j \right)_{j=1}^{J} \in \mathbb{C}^J \\ \bm{f} &:= \left( f_{\bm{n}} \right)_{\bm{n} \in I_\mathbf{N}} \in \mathbb{C}^{|\bm{N}|}\\ \bm{A} &:= \left( \mathrm{e}^{-2 \pi \mathrm{i} \, \bm{n} \cdot \bm{k}_j} \right)_{j=1,\dots,J; \bm{n} \in I_{\mathbf{N}}} \in \mathbb{C}^{J \times |\bm{N}|} -\end{aligned}\]

and where $|\bm{N}| := \text{prod}(\bm{N})$. The adjoint can then be written as

\[ \bm{y} = \bm{A}^{\mathsf{H}} \hat{\bm{f}}\]

where $\bm{y} \in \mathbb{C}^{|\bm{N}|}$.

NFFT

The NFFT is an approximative algorithm that realizes the NDFT in just ${\mathcal O}(|\bm{N}| \log |\bm{N}| + J)$ steps. This is at the same level as the ordinary FFT with the exception of the additional linear term $J$, which is unavoidable since all nodes need to be touched as least once.

The NFFT has two important parameters that influence its accuracy:

From the latter we can derive $\tilde{\bm{N}} = \sigma \bm{N} \in (2\mathbb{N})^D$. As the definition indicates, the oversampling factor $\sigma$ is usually adjusted such that $\tilde{\bm{N}}$ consists of even integers.

The NFFT now approximates $\bm{A}$ by the product of three matrices

\[\bm{A} \approx \bm{B} \bm{F} \bm{D}\]

where

The NFFT is based on the convolution theorem. It applies a convolution in the non-equidistant domain, which is evaluated at equidistant sampling nodes. This convolution is then corrected in the equidistant domain by division with the inverse Fourier transform $\hat{\varphi}$.

The adjoint NFFT matrix approximates $\bm{A}^{\mathsf{H}}$ by

\[\bm{A}^{\mathsf{H}} \approx \bm{D}^{\mathsf{H}} \bm{F}^{\mathsf{H}} \bm{B}^{\mathsf{H}} \]

Implementation-wise, the matrix-vector notation illustrates that the NFFT consists of three independent steps that are performed successively.

Since in practice the multiplication with $\bm{B}$ is also the most expensive step, an NFFT library needs to pay special attention to optimizing it appropriately.

Directional NFFT

In many cases one not just needs to apply a single NFFT but needs to apply many on different data. This leads us to the directional NFFT. The directional NFFT is defined as

\[ f_{\bm{l},j,\bm{r}} := \sum_{ \bm{n} \in I_{\bm{N}_\text{sub}}} \hat{f}_{\bm{l},\bm{n},\bm{r}} \, \mathrm{e}^{-2\pi\mathrm{i}\,\bm{n}\cdot\bm{k}_j}\]

where now $(\bm{l}, \bm{k}, \bm{r}) \in I_\mathbf{N}$ and $I_{\bm{N}_\text{sub}}$ is a subset of $I_{\bm{N}}$. The transform thus maps a $D$-dimensional tensor $\hat{f}_{\bm{l},\bm{n},\bm{r}}$ to an $R$-dimensional tensor $f_{\bm{l},j,\bm{r}}$. $\bm{N}_\text{sub}$ is thus a vector of length $D-R+1$. The indices $\bm{l}$ and $\bm{r}$ can also have length zero. Thus, for $R=1$, the conventional NFFT arises as a special case of the directional NFFT.

Note

The directional NFFT can also be considered to be a slicing of a tensor with subsequent application of a regular NFFT. But the aforementioned formulation can be used to implement a much more efficient algorithm than can be achieved with slicing.

+\end{aligned}\]

and where $|\bm{N}| := \text{prod}(\bm{N})$. The adjoint can then be written as

\[ \bm{y} = \bm{A}^{\mathsf{H}} \hat{\bm{f}}\]

where $\bm{y} \in \mathbb{C}^{|\bm{N}|}$.

NFFT

The NFFT is an approximative algorithm that realizes the NDFT in just ${\mathcal O}(|\bm{N}| \log |\bm{N}| + J)$ steps. This is at the same level as the ordinary FFT with the exception of the additional linear term $J$, which is unavoidable since all nodes need to be touched as least once.

The NFFT has two important parameters that influence its accuracy:

From the latter we can derive $\tilde{\bm{N}} = \sigma \bm{N} \in (2\mathbb{N})^D$. As the definition indicates, the oversampling factor $\sigma$ is usually adjusted such that $\tilde{\bm{N}}$ consists of even integers.

The NFFT now approximates $\bm{A}$ by the product of three matrices

\[\bm{A} \approx \bm{B} \bm{F} \bm{D}\]

where

The NFFT is based on the convolution theorem. It applies a convolution in the non-equidistant domain, which is evaluated at equidistant sampling nodes. This convolution is then corrected in the equidistant domain by division with the inverse Fourier transform $\hat{\varphi}$.

The adjoint NFFT matrix approximates $\bm{A}^{\mathsf{H}}$ by

\[\bm{A}^{\mathsf{H}} \approx \bm{D}^{\mathsf{H}} \bm{F}^{\mathsf{H}} \bm{B}^{\mathsf{H}} \]

Implementation-wise, the matrix-vector notation illustrates that the NFFT consists of three independent steps that are performed successively.

Since in practice the multiplication with $\bm{B}$ is also the most expensive step, an NFFT library needs to pay special attention to optimizing it appropriately.

Directional NFFT

In many cases one not just needs to apply a single NFFT but needs to apply many on different data. This leads us to the directional NFFT. The directional NFFT is defined as

\[ f_{\bm{l},j,\bm{r}} := \sum_{ \bm{n} \in I_{\bm{N}_\text{sub}}} \hat{f}_{\bm{l},\bm{n},\bm{r}} \, \mathrm{e}^{-2\pi\mathrm{i}\,\bm{n}\cdot\bm{k}_j}\]

where now $(\bm{l}, \bm{k}, \bm{r}) \in I_\mathbf{N}$ and $I_{\bm{N}_\text{sub}}$ is a subset of $I_{\bm{N}}$. The transform thus maps a $D$-dimensional tensor $\hat{f}_{\bm{l},\bm{n},\bm{r}}$ to an $R$-dimensional tensor $f_{\bm{l},j,\bm{r}}$. $\bm{N}_\text{sub}$ is thus a vector of length $D-R+1$. The indices $\bm{l}$ and $\bm{r}$ can also have length zero. Thus, for $R=1$, the conventional NFFT arises as a special case of the directional NFFT.

Note

The directional NFFT can also be considered to be a slicing of a tensor with subsequent application of a regular NFFT. But the aforementioned formulation can be used to implement a much more efficient algorithm than can be achieved with slicing.

diff --git a/dev/index.html b/dev/index.html index 2ec1418..87e0a98 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,4 +1,4 @@ Home · NFFT
Edit on GitHub

NFFT.jl

Julia package for the Non-equidistant Fast Fourier Transform

Introduction

This package provides a Julia implementation of the Non-equidistant Fast Fourier Transform (NFFT). For a detailed introduction into the NFFT and its application please have a look at the software paper on the NFFT.jl. Further resources are nfft.org and finufft.readthedocs.io. You

The NFFT is a fast implementation of the Non-equidistant Discrete Fourier Transform (NDFT) that is basically a Discrete Fourier Transform (DFT) with non-equidistant sampling nodes in either Fourier or time/space domain. In contrast to the Fast Fourier Transform (FFT), the NFFT is an approximative algorithm whereas the accuracy can be controlled by two parameters: the window width parameter m and the oversampling factor σ.

The NFFT.jl project serves two different purposes:

  • Provide a package AbstractNFFTs.jl that allows to use any NFFT Julia package such as NFFT3.jl or FINUFFT.jl using the same interface
  • Provide a high-performance, multi-threaded reference implementation in pure Julia. This is realized in the Julia package NFFT.jl.

The term NFFT.jl thus may either mean the entire Github project consisting of several packages or the concrete reference implementation.

Installation

Start julia and open the package mode by entering ]. Then enter

add NFFT

This will install the packages NFFT.jl and all its dependencies. Most importantly it installs the abstract interface package AbstractNFFTs.jl, which NFFT.jl implements.

Additional NFFT related tools can be obtained by adding the package NFFTTools.jl. If you need support for CUDA you also need to install the package CuNFFT.jl.

In case you want to use an alternative NFFT implementation such as NFFT3.jl or FINUFFT.jl we provide wrapper types allowing to use them as AbstractNFFTs implementations. They can be used like this:

julia> using AbstractNFFTs
 julia> include(joinpath(dirname(pathof(AbstractNFFTs)), "..", "..", "Wrappers", "FINUFFT.jl"))
-julia> include(joinpath(dirname(pathof(AbstractNFFTs)), "..", "..", "Wrappers", "NFFT3.jl"))

This requires that you first add the package you want to use.

Guide

  • The documentation starts with the Mathematical Background that properly defines the NDFT, the NFFT and its directional variants. You might want to skip this part if you are familiar with the notation and concepts of the NFFT.
  • Then, an Overview about the usage of the NFFT functions is given in a tutorial style manner.
  • Then, an overview about Accuracy and Performance is given.
  • The section about Tools introduced some high-level functions that build upon the NFFT. For instance NFFT inversion is discussed in that section.
  • In the section about the Abstract Interface for NFFTs we outline how the package is divided into an interface package and implementation packages. This part is useful if you plan to use different NFFT implementations, e.g. one for the CPU and another for the GPU and would like to switch.
  • Finally, the documentation contains an API index.

License / Terms of Usage

The source code of this project is licensed under the MIT license. This implies that you are free to use, share, and adapt it. However, please give appropriate credit by citing the project. You can do so by citing the publication

T. Knopp, M. Boberg and M. Grosser, NFFT.jl: Generic and Fast Julia Implementation of the Nonequidistant Fast Fourier Transform, 2022 arXiv:2208.00049

A BibTeX file NFFT.bib can be found in the root folder of the Github repository.

Contact

If you have problems using the software, find bugs or have ideas for improvements please use the issue tracker. For general questions please use the discussions section on Github.

Contributors

A complete list of contributors can be found on the Github page.

+julia> include(joinpath(dirname(pathof(AbstractNFFTs)), "..", "..", "Wrappers", "NFFT3.jl"))

This requires that you first add the package you want to use.

Guide

License / Terms of Usage

The source code of this project is licensed under the MIT license. This implies that you are free to use, share, and adapt it. However, please give appropriate credit by citing the project. You can do so by citing the publication

T. Knopp, M. Boberg and M. Grosser, NFFT.jl: Generic and Fast Julia Implementation of the Nonequidistant Fast Fourier Transform, 2022 arXiv:2208.00049

A BibTeX file NFFT.bib can be found in the root folder of the Github repository.

Contact

If you have problems using the software, find bugs or have ideas for improvements please use the issue tracker. For general questions please use the discussions section on Github.

Contributors

A complete list of contributors can be found on the Github page.

diff --git a/dev/overview/index.html b/dev/overview/index.html index 5b63a68..c8312d8 100644 --- a/dev/overview/index.html +++ b/dev/overview/index.html @@ -21,4 +21,4 @@ p1 = plan_nfft(y, N, dims=1) f = randn(ComplexF64, N) fHat = p1 * f

Here size(f) = (16,20) and size(fHat) = (11,20) since we compute an NFFT along the first dimension. To compute the NFFT along the second dimension

p2 = plan_nfft(y, N, dims=2)
-fHat = p2 * f

Now size(fHat) = (16,11).

+fHat = p2 * f

Now size(fHat) = (16,11).

diff --git a/dev/performance/index.html b/dev/performance/index.html index cbdad60..ba4e76e 100644 --- a/dev/performance/index.html +++ b/dev/performance/index.html @@ -1,2 +1,2 @@ -Performance · NFFT
Edit on GitHub

Accuracy and Performance

On this page, the accuracy and the performance of NFFT.jl are investigated. For comparison we use the C library NFFT3 and the C++ library FINUFFT. The parameters for the benchmark are

  • $\bm{N}_\text{1D}=(512*512,), \bm{N}_\text{2D}=(512,512), \bm{N}_\text{3D}=(64,64,64)$
  • $J= \vert \bm{N} \vert$
  • $m=3, \dots, 8$
  • $\sigma = 2$
  • POLYNOMIAL and TENSOR precomputation
  • 1 thread
  • random nodes

All benchmarks are performed with @belapsed from BenchmarkTools.jl which takes the minimum of several runs (120 s upper benchmark time). The benchmark is run on a computer with 2 AMD EPYC 7702 CPUs running at 2.0 GHz (256 cores in total) and a main memory of 1024 GB. The benchmark suite is described here.

The results for $D=1,\dots,3$ are shown in the following graphic illustrating the accuracy (x-axis) versus the performance (y-axis) for various $m$.

Performance vs Accurracy 1D

The results show that NFFT.jl is one of the fastest NFFT libraries. One can chose between shorter precomputation time using POLYNOMIAL precomputation or faster transforms using TENSOR precomputation.

+Performance · NFFT
Edit on GitHub

Accuracy and Performance

On this page, the accuracy and the performance of NFFT.jl are investigated. For comparison we use the C library NFFT3 and the C++ library FINUFFT. The parameters for the benchmark are

  • $\bm{N}_\text{1D}=(512*512,), \bm{N}_\text{2D}=(512,512), \bm{N}_\text{3D}=(64,64,64)$
  • $J= \vert \bm{N} \vert$
  • $m=3, \dots, 8$
  • $\sigma = 2$
  • POLYNOMIAL and TENSOR precomputation
  • 1 thread
  • random nodes

All benchmarks are performed with @belapsed from BenchmarkTools.jl which takes the minimum of several runs (120 s upper benchmark time). The benchmark is run on a computer with 2 AMD EPYC 7702 CPUs running at 2.0 GHz (256 cores in total) and a main memory of 1024 GB. The benchmark suite is described here.

The results for $D=1,\dots,3$ are shown in the following graphic illustrating the accuracy (x-axis) versus the performance (y-axis) for various $m$.

Performance vs Accurracy 1D

The results show that NFFT.jl is one of the fastest NFFT libraries. One can chose between shorter precomputation time using POLYNOMIAL precomputation or faster transforms using TENSOR precomputation.

diff --git a/dev/search/index.html b/dev/search/index.html index bb5ff50..b6328ef 100644 --- a/dev/search/index.html +++ b/dev/search/index.html @@ -1,2 +1,2 @@ -Search · NFFT

Loading search...

    +Search · NFFT

    Loading search...

      diff --git a/dev/tools/index.html b/dev/tools/index.html index 92fd515..653f11d 100644 --- a/dev/tools/index.html +++ b/dev/tools/index.html @@ -13,4 +13,4 @@ y = randn(ComplexF32, Ñ) convolveToeplitzKernel!(y, λ) # multiply with Gram matrix - +