TimeSegmentation¶

template <typename T1, typename Tobj> class TimeSegmentation

Off-resonance correction via time-segmented interpolation. Off-resonance correction via time-segmented interpolation.

Wraps any encoding operator Tobj (typically Gnufft) with field-map-based off-resonance correction, replacing the exact but \(O(N^2)\) field-corrected DFT (Gdft) with \(L\) fast NUFFT evaluations.

The problem: The field-corrected MRI signal model is:

\[d_j = \sum_k x_k \, e^{-i\omega_k t_j} \, e^{-i 2\pi \mathbf{k}_j \cdot \mathbf{r}_k}\]

The spatially-varying phase \(e^{-i\omega_k t_j}\) couples every k-space sample \(j\) to every image pixel \(k\) through the field map, preventing the use of an FFT. Computing this directly (as Gdft does) costs \(O(N_d \cdot N_i)\).

The approximation: Time segmentation factors the off-resonance phase into a temporal interpolation and a spatially-varying modulation:

\[e^{-i\omega(\mathbf{r}) t} \approx \sum_{l=1}^{L} b_l(t) \, e^{-i\omega(\mathbf{r}) \tau_l}\]

where \(\tau_l\) are \(L\) segment center times and \(b_l(t)\) are temporal interpolation coefficients (stored in the matrix AA).

Substituting into the signal model:

\[d_j \approx \sum_{l=1}^{L} b_l(t_j) \underbrace{\sum_k \left[ x_k \, e^{-i\omega_k \tau_l} \right] e^{-i 2\pi \mathbf{k}_j \cdot \mathbf{r}_k}}_{\text{NUFFT of phase-modulated image}}\]

The inner sum is a standard NUFFT of the image modulated by a time-independent phase \(e^{-i\omega_k \tau_l}\). This replaces one \(O(N^2)\) DFT with \(L\) fast \(O(N \log N)\) NUFFTs plus element-wise phase modulations and temporal interpolation.

Forward (image → k-space with field correction):

graph LR
    X["$$\mathbf{x}$$"] --> PHI1["$$e^{-i\omega\tau_1} \circ \mathbf{x}$$"]
    X --> PHI2["$$e^{-i\omega\tau_2} \circ \mathbf{x}$$"]
    X --> PHIL["$$e^{-i\omega\tau_L} \circ \mathbf{x}$$"]
    PHI1 --> G1["NUFFT"]
    PHI2 --> G2["NUFFT"]
    PHIL --> GL["NUFFT"]
    G1 --> B1["$$b_1(t) \circ$$"]
    G2 --> B2["$$b_2(t) \circ$$"]
    GL --> BL["$$b_L(t) \circ$$"]
    B1 --> SUM["$$\sum_l$$"]
    B2 --> SUM
    BL --> SUM
    SUM --> D["$$\mathbf{d}$$"]

Adjoint (k-space → image with field correction):

graph LR
    D["$$\mathbf{d}$$"] --> B1["$$b_1^*(t) \circ$$"]
    D --> B2["$$b_2^*(t) \circ$$"]
    D --> BL["$$b_L^*(t) \circ$$"]
    B1 --> GH1["NUFFT$$^H$$"]
    B2 --> GH2["NUFFT$$^H$$"]
    BL --> GHL["NUFFT$$^H$$"]
    GH1 --> PHI1["$$e^{+i\omega\tau_1} \circ$$"]
    GH2 --> PHI2["$$e^{+i\omega\tau_2} \circ$$"]
    GHL --> PHIL["$$e^{+i\omega\tau_L} \circ$$"]
    PHI1 --> SUM["$$\sum_l$$"]
    PHI2 --> SUM
    PHIL --> SUM
    SUM --> X["$$\hat{\mathbf{x}}$$"]

Interpolation types: - Hanning window (type 1): fast to precompute, moderate accuracy - Min-max / Fessler (type 2): more accurate but slower precomputation - Approximate min-max (type 3): balanced accuracy and speed

\(L\) controls the accuracy-vs-speed tradeoff. More segments = better approximation of the off-resonance phase but more NUFFT evaluations per iteration. Rules of thumb for choosing \(L\): - Hanning interpolation: ~1 segment per 1 ms of readout length - Min-max interpolation: ~1 segment per 2--3 ms of readout length (more accurate interpolation allows fewer segments)

For example, a 10 ms spiral readout would need ~10 segments with Hanning or ~4--5 segments with min-max.

Gnufft<float> G(n1, 2.0f, Nx, Ny, Nz, kx, ky, kz, ix, iy, iz);
TimeSegmentation<float, Gnufft<float>> Gts(G, FM, t, L, interptype, gridos);
SENSE<float, TimeSegmentation<float, Gnufft<float>>> S(Gts, SMap, nc);
auto xhat = solve_pwls_pcg<float>(x0, S, W, yi, R, niter);

See : Sutton et al., "Fast, Iterative Image Reconstruction for MRI in the Presence of Field Inhomogeneities," IEEE TMI, 2003. https://doi.org/10.1109/TMI.2002.808360

See : Man et al., "Multifrequency Interpolation for Fast Off-Resonance Correction," IEEE TMI, 1997. https://doi.org/10.1109/42.611354

See : Gdft for brute-force field correction (exact but slow)

T1 : Floating-point precision type (float or double).

Tobj : Underlying encoding operator type (e.g., Gnufft<T1>).

Variables¶

Name	Description
n1	Number of k-space samples (output length of forward transform).
n2	Number of image pixels (input length of forward transform).
L	Number of time segments L.
type	Interpolation type: 1 = Hanning, 2 = exact min-max, 3 = approx min-max.
Nshots	Number of shots; used to reduce complexity of interpolator calculation.
tau	Time-segment length tau (seconds).
T_min	Minimum time in the time vector (e.g., TE for spiral-out), in seconds.
obj	Pointer to the underlying single-coil encoding operator.
fieldMap	Off-resonance field map (rad/s), length n2.
timeVec	Per-sample readout time vector (s), length n1.
AA	Interpolation coefficient matrix, size L x n1.
conjAA	Cached conj(AA) for adjoint operator — avoids recomputing per call.
Wo	Phase modulation matrices Wo (n2 x L) for forward transform.
WoH	Conjugate phase modulation matrices WoH (n2 x L) for adjoint transform.

Operators¶

Name	Description
operator*	Forward transform: image -> k-space with off-resonance correction. `d` : Input image vector of length `n2.` Return : Output k-space vector of length `n1.`
operator/	Adjoint transform: k-space -> image with off-resonance correction. `d` : Input k-space vector of length `n1.` Return : Output image vector of length `n2.`

Functions¶

Name	Description
TimeSegmentation	Default constructor.
TimeSegmentation	Construct a TimeSegmentation operator.
calcFFT1D	Compute 1-D FFT along the time dimension for interpolator construction.

Variable Details¶

AA ¶

Mat<CxT1> AA

Interpolation coefficient matrix, size L x n1.

L ¶

int L

Number of time segments L.

Nshots ¶

uword Nshots

Number of shots; used to reduce complexity of interpolator calculation.

T_min ¶

T1 T_min

Minimum time in the time vector (e.g., TE for spiral-out), in seconds.

Wo ¶

Mat<CxT1> Wo

Phase modulation matrices Wo (n2 x L) for forward transform.

WoH ¶

Mat<CxT1> WoH

Conjugate phase modulation matrices WoH (n2 x L) for adjoint transform.

conjAA ¶

Mat<CxT1> conjAA

Cached conj(AA) for adjoint operator — avoids recomputing per call.

fieldMap ¶

Col<T1> fieldMap

Off-resonance field map (rad/s), length n2.

n1 ¶

uword n1

Number of k-space samples (output length of forward transform).

n2 ¶

uword n2

Number of image pixels (input length of forward transform).

obj ¶

Tobj* obj

Pointer to the underlying single-coil encoding operator.

tau ¶

T1 tau

Time-segment length tau (seconds).

timeVec ¶

Col<T1> timeVec

Per-sample readout time vector (s), length n1.

type ¶

uword type

Interpolation type: 1 = Hanning, 2 = exact min-max, 3 = approx min-max.

Operator Details¶

operator*¶

Col<CxT1> operator*(const Col<CxT1>& d) const

Forward transform: image -> k-space with off-resonance correction.

d : Input image vector of length n2.

Return : Output k-space vector of length n1.

operator/¶

Col<CxT1> operator/(const Col<CxT1>& d) const

Adjoint transform: k-space -> image with off-resonance correction.

d : Input k-space vector of length n1.

Return : Output image vector of length n2.

Function Details¶

TimeSegmentation ¶

TimeSegmentation()

Default constructor.

TimeSegmentation( Tobj& G, Col<T1> map_in, Col<T1> timeVec_in, uword a, uword b, uword c, uword interptype = 1, uword shots = 1)

Construct a TimeSegmentation operator.

G : Reference to the underlying encoding operator.

map_in : Off-resonance field map in rad/s, length n2.

timeVec_in : Per-sample readout times in seconds, length n1.

a : Number of k-space samples (n1).

b : Number of image pixels (n2).

c : Number of time segments (L).

interptype : Interpolation method: 1=Hanning, 2=exact min-max, 3=approx min-max (default 1).

shots : Number of shots for segmentation (default 1).

calcFFT1D ¶

Col<CxT1> calcFFT1D(const Col<CxT1>& d, uword KK) const

Compute 1-D FFT along the time dimension for interpolator construction.

d : Input complex vector.

KK : FFT length.

Return : FFT output vector of length KK.

TimeSegmentation¶

Variables¶

Operators¶

Functions¶

Variable Details¶

AA¶

L¶

Nshots¶

T_min¶

Wo¶

WoH¶

conjAA¶

fieldMap¶

n1¶

n2¶

obj¶

tau¶

timeVec¶

type¶

Operator Details¶

operator*¶

operator/¶

Function Details¶

TimeSegmentation¶

calcFFT1D¶

AA ¶

L ¶

Nshots ¶

T_min ¶

Wo ¶

WoH ¶

conjAA ¶

fieldMap ¶

n1 ¶

n2 ¶

obj ¶

tau ¶

timeVec ¶

type ¶

TimeSegmentation ¶

calcFFT1D ¶