Model implementation

This page collects implementation notes that complement the auto-generated API documentation (see nifreeze.model package). It documents design decisions and empirical characterizations that are specific to NiFreeze’s code — not the underlying domain theory, which lives in the project’s grounding knowledge base.

Generalized q-Sampling Imaging (GQI)

GeneralizedQSamplingModel is a vendored, lightly adapted copy of the GQI reconstruction of Yeh et al., 2010 as implemented in DIPY. The forward kernel (gqi_kernel(), method="standard") implements the paper’s spin-distribution-function reconstruction (Eq. 6/9) verbatim — it matches DIPY’s standard kernel and a direct transcription of Eq. 9 to machine precision.

The model and its fit descend from DIPY’s ODF-model family (OdfModel / OdfFit), and odf() provides the forward signal → SDF transform expected of that family. Unlike DIPY’s ODF models (and NiFreeze’s DTI/DKI), GQI deliberately does not use the @multi_voxel_fit decorator: its fit is a linear operator that merely stores the data, and prediction is a single vectorized matmul over all voxels, so per-voxel looping would only discard that vectorization.

GQI as a signal predictor (a NiFreeze extension)

Yeh (2010) defines only the forward map signal → SDF. NiFreeze additionally needs the inverse to predict the diffusion signal given a b-vector (and -value). predict() composes the forward GQI kernel with a Tikhonov-regularized reconstruction kernel (prediction_kernel()):

\[\mathbf{K}^{+} = (\mathbf{K}\mathbf{K}^{\mathsf T} + \lambda_0 \mathbf{I})^{-1} \mathbf{K}, \qquad \lambda_0 = 10^{-6}.\]

This is the regularized least-squares signal whose forward GQI transform reproduces the fitted SDF. It is a NiFreeze modeling choice and is not part of Yeh (2010).

Reconstruction fidelity and the intercept behaviour

Applied to the signal, the round-trip operator is \(\mathbf{P} = \mathbf{A}(\mathbf{A}+\lambda_0\mathbf{I})^{-1}\) with \(\mathbf{A} = \mathbf{K}\mathbf{K}^{\mathsf T}\). Empirically, on diffusion-weighted (b > 0) volumes:

  • \(\mathbf{P}\) preserves the mean/DC exactly (\(\mathbf{P}\cdot\mathbf 1 = \mathbf 1\)) and is scale-homogeneous (a per-voxel S0 cancels). So there is no missing additive baseline/S0 term, and normalizing by S0 changes nothing.

  • What it does lose is angular signal energy: \(\mathbf{P}\) is a contraction that projects out the part of the signal lying in the null space of the GQI basis. The per-voxel regression slope of predicted vs. observed is the fraction of angular variance the basis can represent (≈ 0.76 on single-shell HARDI, ≈ 0.99 on rich grid sampling).

  • The apparent non-zero intercept is not an independent term because the mean is preserved: intercept = mean · (1 − slope) exactly. There is a single effect (angular shrinkage), not two.

Consequently the missing degree of freedom is a representation gap, worse for q-space-poor acquisitions (single-shell) than for grid/multi-shell data, not an S0 intercept. Two practical corollaries:

  • b=0 volumes are excluded from GQI fitting/prediction in NiFreeze. Feeding raw b=0 signal into the fit without an explicit constant term degrades the diffusion-weighted reconstruction (the large b=0 amplitude is not representable by the diffusion-weighted sinc basis); excluding b=0 avoids this cleanly.

  • method="gqi2" is a weaker signal predictor than the default "standard": its kernel is more oscillatory and ill-conditioned, so its round-trip does not preserve signal scale and its leave-one-volume-out correlation falls below that of "standard". This is why NiFreeze defaults to "standard" even though DIPY’s GQI defaults to "gqi2".

For motion estimation this amplitude shrinkage is benign (registration keys on relative contrast, which the high correlation preserves); it would matter for a downstream quantitative use of the predicted signal.

Sphere density and the default recursion level

The SDF is sampled on an icosahedral sphere whose subdivision recursion_level (a parameter of GeneralizedQSamplingModel) sets its vertex count. The default, DEFAULT_SPHERE_RECURSION_LEVEL = 5 (1026 vertices), was chosen from the following experiment.

Settings. method="standard", sampling_length \(\sigma = 1.2\), INVERSE_LAMBDA \(\lambda_0 = 10^{-6}\). Metric: per-voxel round-trip slope (fit the model on the diffusion-weighted set, reconstruct that same set, then take the least-squares slope of predicted vs. observed, averaged over voxels) and the mean per-voxel correlation. No direction is held out: this is a self-consistency measure that isolates the reconstruction, not a leave-one-out prediction.

Data. (1) Grid / multi-shell: DIPY’s dsi_voxels() phantom, b=0 excluded → 101 diffusion-weighted directions, 600 voxels. (2) Single-shell: the Stanford HARDI dataset (DIPY get_fnames("stanford_hardi")), central slab [20:60, 40:80, 28:45], brain = b=0 reference above its 60th percentile → 150 diffusion-weighted directions, 10 880 voxels. b=0 excluded in both.

Table 1 Round-trip fidelity vs. sphere subdivision (standard method)

recursion_level

vertices

grid slope

grid corr

single-shell slope

single-shell corr

2

18

0.759

0.839

0.339

0.377

3

66

0.829

0.952

0.519

0.700

4

258

0.931

0.974

0.750

0.869

5 (default)

1026

0.965

0.980

0.765

0.879

6

4098

0.992

0.983

0.788

0.898

Reading. Both curves rise steeply up to a knee at recursion_level 4–5; the operator’s effective rank saturates there (beyond it, a denser sphere only refines the quadrature). Past the knee the two data types diverge: the grid acquisition keeps improving toward a near-perfect round-trip (slope → 0.99) — it is sphere-limited — whereas the single-shell acquisition plateaus near 0.78 — it is q-space-limited, and no sphere density recovers the angular signal a single shell cannot encode.

Decision. The default recursion_level = 5 sits just past the knee for both regimes, so it is adequate for the single-shell data NiFreeze typically sees while leaving headroom for grid/multi-shell acquisitions, where a higher value (e.g. 6) keeps paying off. Denser spheres cost compute (O(vertices) per kernel build) for no single-shell benefit, so they are opt-in rather than the default.

Reproduction.

import numpy as np
from dipy.core.gradients import gradient_table
from dipy.core.subdivide_octahedron import create_unit_sphere
from dipy.data import dsi_voxels
from nifreeze.model.gqi import gqi_kernel

L, lam = 1.2, 1e-6
data, gtab = dsi_voxels()
dw = np.where(~gtab.b0s_mask)[0]
g = gradient_table(bvals=gtab.bvals[dw], bvecs=gtab.bvecs[dw])
d2d = data.reshape(-1, data.shape[-1])[:, dw]
for r in (2, 3, 4, 5, 6):
    K = gqi_kernel(g, L, create_unit_sphere(recursion_level=r), method="standard")
    A = K @ K.T
    P = np.linalg.inv(A + lam * np.eye(A.shape[0])) @ A
    p = np.maximum(d2d @ P.T, 0)
    dc, pc = d2d - d2d.mean(1, keepdims=True), p - p.mean(1, keepdims=True)
    slope = ((dc * pc).sum(1) / ((dc**2).sum(1) + 1e-12)).mean()
    print(r, len(create_unit_sphere(recursion_level=r).vertices), round(float(slope), 3))