Operators & Solvers¶

Operator Convention¶

All encoding operators in forge follow a unified interface:

Expression	Meaning
`G * x`	Forward transform: image → k-space
`G / y`	Adjoint transform: k-space → image

Both operators accept and return arma::Col<std::complex<T1>> vectors at the Armadillo interface boundary. Internally, hot paths use forgeCol/forgeMat for GPU dispatch (see forge Types).

Encoding Operators¶

Class	Description
`Gdft<T1>`	Direct non-uniform DFT with optional off-resonance correction
`GdftR2<T1>`	Gdft extended with R2* field map gradient support
`Gnufft<T1>`	Fast NUFFT via Kaiser-Bessel gridding (CPU + GPU)
`Gfft<T1>`	Uniform Cartesian FFT (wraps FFTW / cuFFT)
`SENSE<T1, Tobj>`	Multi-coil sensitivity encoding wrapping any encoding operator
`pcSENSE<T1>`	Phase-corrected SENSE with per-shot phase maps
`pcSenseTimeSeg<T1>`	pcSENSE combined with time-segmented field correction
`TimeSegmentation<T1, Tobj>`	Off-resonance correction via time-segmented interpolation

Gnufft Pipeline¶

The NUFFT forward and adjoint transforms are multi-step pipelines:

graph LR
    subgraph "Forward (image → k-space)"
        A[Image] --> B[Deapodize]
        B --> C[Zero-pad]
        C --> D[Oversampled FFT]
        D --> E[KB Interpolation]
        E --> F[k-space samples]
    end

graph RL
    subgraph "Adjoint (k-space → image)"
        G[k-space samples] --> H[KB Gridding]
        H --> I[Inverse FFT]
        I --> J[Crop]
        J --> K[Deapodize]
        K --> L[Image]
    end

Density Compensation¶

Non-uniform k-space sampling requires density compensation weights (DCF) to correct for variable sampling density before adjoint gridding. forge uses the Pipe iterative algorithm:

\[W_{n+1} = \frac{W_n}{G^H G \, W_n}\]

where \(G^H G\) is the forward-then-adjoint gridding operation. This converges in ~15-20 iterations. Implemented in calc3DDensityCompensation().

Regularization¶

All regularization operators inherit from Robject<T1>:

Class	Description
`Robject<T1>`	Abstract base class for penalty functions
`QuadPenalty<T1>`	Quadratic (Tikhonov) penalty: \(R(\mathbf{x}) = \frac{\beta}{2}\\|C\mathbf{x}\\|^2\)
`TVPenalty<T1>`	Smooth TV (Charbonnier) penalty: \(\psi(d) = \delta^2(\sqrt{1+(d/\delta)^2}-1)\)

Solvers¶

Function	Description
`solve_pwls_pcg`	Preconditioned weighted least-squares via conjugate gradient (Polak-Ribiere-Polyak)
`solve_grad_desc`	Gradient descent solver
`reconSolve`	High-level reconstruction entry point

The primary solver minimizes:

\[\hat{\mathbf{x}} = \arg\min_{\mathbf{x}} \|\mathbf{W}^{1/2}(\mathbf{y} - A\mathbf{x})\|^2 + R(\mathbf{x})\]

Quick Example¶

#include "forge/Operators/Gdft.h"
#include "forge/Penalties/QuadPenalty.h"
#include "forge/Solvers/solve_pwls_pcg.hpp"
using namespace arma;

// Build a Cartesian Gdft operator (8x8 image, 64 k-space samples)
uword Nx = 8, Ny = 8, n1 = Nx*Ny, n2 = Nx*Ny;
Col<float> kx(n2), ky(n2), kz(n2, fill::zeros);
Col<float> ix(n1), iy(n1), iz(n1, fill::zeros);
// ... fill kx/ky with Cartesian grid, ix/iy with pixel positions ...
Col<float> FM(n1, fill::zeros), t(n2, fill::zeros);
Gdft<float> G(n1, n2, kx, ky, kz, ix, iy, iz, FM, t);

// Forward transform
Col<cx_float> image(n1, fill::ones);
Col<cx_float> kspace = G * image;

// Adjoint transform
Col<cx_float> back = G / kspace;

// Iterative reconstruction via PCG
Col<float> W(n2, fill::ones);
QuadPenalty<float> R(Nx, Ny, 1, 0.01f);
Col<cx_float> x0(n1, fill::zeros);
Col<cx_float> xhat = solve_pwls_pcg<float>(x0, G, W, kspace, R, 50);