Operators

FORGE treats an MRI scan as a linear system. The image x produces the measured k-space data y through an encoding operator A:

\[\mathbf{y} = A\mathbf{x} + \text{noise}\]

A holds everything about the acquisition: the trajectory, the sample times, the coil sensitivities, the field map, and the shot phases. Reconstruction means finding the x that best explains y under A (see The solver and regularization). The solver never needs A as a matrix, which would be far too large. It only needs to apply A and its adjoint A^H to vectors, and every operator class in FORGE provides exactly those two things.

The operator convention

Expression	Meaning
A * x	Forward transform, image to k-space
A / y	Adjoint transform, k-space to image

/ here is the adjoint, not a division or an inverse. Each operator has overloads for arma::Col<std::complex<T>> and for forgeCol<forgeComplex<T>>. The solver uses the forgeCol versions where they exist, which keeps the data on the GPU between iterations.

Operators wrap operators

Each part of the physics is a separate class, and they nest. A class that adds coils, field correction or shot phases takes the next operator in as a template parameter and calls it:

Class	Adds
Gnufft	Non-uniform FFT by Kaiser-Bessel gridding (see The NUFFT)
Gdft	Exact non-uniform DFT, with field correction built in
GdftR2	Gdft plus field-map gradients within each voxel
Gfft	Cartesian FFT
TimeSegmentation<T, Tobj>	Field correction around any of the above (see Field correction)
SENSE<T, Tobj>	Coil sensitivities: one application of Tobj per coil
pcSENSE	Coil sensitivities and a phase map per shot, with a Gdft per shot
pcSenseTimeSeg	The same, with a time-segmented Gnufft per shot

So the command-line programs build these stacks:

Program	Operator
forgeSense -F NUFFT	SENSE<float, TimeSegmentation<float, Gnufft<float>>>
forgeSense -F DFT	SENSE<float, Gdft<float>>
forgeSense -F DFTGrads	SENSE<float, GdftR2<float>>
forgePcSense	pcSENSE<float>
forgePcSenseTimeSeg	pcSenseTimeSeg<float>
oscillateRecon	LRobj around one pcSenseTimeSeg<float> per repetition (see Low-rank reconstruction)

Because each layer only calls the one inside it, a new trajectory type or a new correction can be added as one class without touching the others.

Shot phase correction

In multi-shot imaging, motion during diffusion or MRE encoding gives each shot its own phase across the image. SENSE assumes all shots see the same object, so these phases become ghosts. The pcSENSE operators multiply the image by each shot's phase before encoding it:

\[\mathbf{y}_{c,s} = G_s\left(e^{\,i P_s} \circ S_c \circ \mathbf{x}\right)\]

for coil c and shot s, where S_c is the coil sensitivity, P_s the shot's phase map, and G_s the shot's Fourier operator. The phase maps come from navigators or a low-resolution reconstruction of each shot, computed before FORGE runs.

Precision and backends

The operators are templates on the floating-point type, with float and double instantiations. The command-line programs use float. GPU code exists only for float; the double versions always run on the CPU. GPU backends describes which parts of each operator run on the GPU.