Diffuse Optical Tomography

Diffuse Optical Tomography

Diffuse Optical Tomography (DOT) in pf2 reconstructs an image-space map of hemoglobin change on the cortical surface by inverting an atlas forward sensitivity model, rather than interpolating channel values for on-screen display. This page is for researchers who want to move beyond channel-space topography toward a physically grounded cortical image, and covers the forward model, image reconstruction, rendering, and the fidelity diagnostics that keep expectations honest for sparse montages.

DOT vs. channel-space topography

pf2.probe.plot.topo and pf2.probe.project.biomarker interpolate channel values across the scalp/cortex with a lateral (Gaussian) kernel — fast and useful for a first look, but not physically grounded in how light actually samples tissue. DOT instead:

1. Builds a forward sensitivity matrix describing how light travels between each source-detector pair and reaches each point on the cortex (the "banana"-shaped photon-measurement-density, or PMDF).
2. Inverts that model to estimate a vertex-space hemoglobin map from the channel measurements — an actual image reconstruction, not a interpolation.

Resolution is coarse (~2-3 cm) for a sparse single-distance montage — most standard fNIRS probes. The infrastructure is valuable regardless (honest projections, coverage masking, reconstruction diagnostics), and resolution sharpens with high-density / multi-distance / short-separation montages (see Fidelity and high-density montages below).

proc = processFNIRS2(pf2.import.sampleData.fNIR2000());   % needs 3D geometry

All DOT functions require the recording's device to carry 3D optode coordinates (dev.hasMNI()); layout-only 2D-schematic devices cannot be reconstructed.

Forward model

pf2.probe.forward.sensitivity builds the channel-by-vertex PMDF sensitivity matrix A(c, v) on the bundled MNI cortical surface — the operator every other DOT function is built on:

[A, mesh] = pf2.probe.forward.sensitivity(proc);
% A [nChannels x nVertices] sparse
[A, mesh, info] = pf2.probe.forward.sensitivity(proc, 'Wavelength', 850);

Key name-value options: 'Wavelength' ('mean' default, a numeric value, or 'all' for a per-wavelength cell array); 'ScalpOffset' (mm, default 12); 'MaxDistance' (mm, default 50); 'Prune' (sparsity threshold); 'mua', 'musp' optical-property overrides; 'HeadModel' ('homogeneous' default or 'layered'). Results are memoised per (device, options), so repeated calls in a batch loop are cheap.

pf2.probe.forward.coverage collapses the forward model to a single per-vertex sensitivity field — where the montage can actually "see":

cov = pf2.probe.forward.coverage(proc);
supported = cov > 0.1;

This coverage field is what gates masking in reconstruction and rendering — no signal is ever invented outside the montage's optical support.

Honest channel projection: the PMDF "banana"

For a quick, physically grounded display of raw channel values (not a full inverse), pf2.probe.project.pmdf backprojects through the real sensitivity footprint instead of a Gaussian lateral kernel:

meanHbO = mean(proc.HbO(100:200, :), 1, 'omitnan');
pf2.probe.project.pmdf(meanHbO, proc, ...
    'initCamPosition', 'front-left', 'savePath', 'banana.png');

This is a display smoothing through the true forward operator, not a deconvolution — for that, reconstruct (next section).

Image reconstruction

pf2.probe.dot.reconstruct inverts the forward model with a depth-weighted, channel-whitened Tikhonov minimum-norm estimator, masked to the montage's coverage support:

recon = pf2.probe.dot.reconstruct(proc);                      % whole-recording mean
recon = pf2.probe.dot.reconstruct(proc, 'Time', [10 30]);     % window mean
recon = pf2.probe.dot.reconstruct(proc, 'Time', 30);           % nearest sample
recon = pf2.probe.dot.reconstruct(proc, 'Biomarkers', {'HbO'}, 'AllTimes', true);  % [T x nV]
recon = pf2.probe.dot.reconstruct(proc, 'Lambda', 1e9, 'DepthWeight', false);       % fixed reg, no depth prior

Key name-value options:

Parameter	Purpose
'Biomarkers'	Fields to reconstruct (default {'HbO','HbR'})
'Time'	[] reconstructs the recording time-mean; a scalar t the nearest sample; [t0 t1] the window mean. Ignored when 'AllTimes' is true
'AllTimes'	Reconstruct every sample → [T x nV] per biomarker (memory-heavy; warns above ~5e7 elements)
'Mask' / 'MaskThreshold'	NaN out vertices outside the coverage support (default true, cutoff 0.05)
'Prior'	'minnorm' (default, depth-weighted min-norm), 'laplacian' (graph-Laplacian smoothness), or 'parcel' (penalizes within-Brodmann-area variance)
'DepthWeight', 'Whiten', 'Lambda', 'LambdaFraction', 'RegMethod'	Forwarded to pf2_base.dot.reconstructImage; 'Lambda' fixes regularization, otherwise auto-chosen by GCV ('RegMethod','gcv', default) or L-curve
'ScalpRegression', 'ScalpMethod', 'ShortSepThreshold'	Remove superficial signal via short-separation channel regression before inversion (see below)
'HighRes', 'ScalpOffset', 'MaxDistance', 'Prune'	Forwarded to pf2.probe.forward.sensitivity

The returned recon struct:

Field	Contents
.<Biomarker>	[1 x nV] (time-mean/point) or [T x nV] (AllTimes)
.vertices, .faces, .brodmann	Cortical mesh geometry
.coverage	[1 x nV] normalized sensitivity field
.mask	[1 x nV] logical support used
.time	Scalar / [1x2] / [T x 1] time reconstructed
.units	Biomarker units, suffixed (relative)
.lambda	Regularization used, per biomarker
.meta	Forward-model .info plus inverse metadata

Units caveat: the reconstruction is concentration-domain through a geometric (unnormalized) sensitivity operator, so its magnitude is relative, not calibrated µM — recon.units is suffixed (relative) to make this explicit. A spectral ΔOD→chromophore reconstruction (absolute units) is a planned extension.

Bad channels (fchMask == 0) are automatically excluded from the inverse with a warning, so a masked-out channel cannot bias the image or the auto-lambda selection.

Rendering the reconstruction

pf2.probe.project.tomography paints a recon struct onto the cortical surface, masked to the sensitivity field so unsupported cortex stays gray:

pf2.probe.project.tomography(recon, 'Biomarker', 'HbO', ...
    'initCamPosition', 'front-left', 'savePath', 'dot_HbO.png');
pf2.probe.project.tomography(recon, 'Biomarker', 'HbR', 'savePath', 'dot_HbR.png');

'Range' sets colorbar limits (default symmetric auto), 'CoverageAlpha' fades opacity by sensitivity (default true, honest-confidence display), and 'TimeIndex' selects a frame from an 'AllTimes' reconstruction (default: the peak-absolute-value frame). Other options ('savePath', 'Style', 'initCamPosition') forward to interpolateValues3D — as with all 3D renders, use 'savePath' for headless output rather than figure('Visible','off') + saveas.

For a time-resolved view, pf2.probe.plot.tomographyMovie reconstructs a sequence of frames sharing one regularization and one global color scale:

pf2.probe.plot.tomographyMovie(proc, 'Biomarker', 'HbO', ...
    'TimeRange', [0 30], 'NFrames', 40, 'FPS', 15, 'savePath', 'dot_HbO_movie.mp4');

'Prior', 'ScalpRegression', 'DepthWeight', 'Lambda', etc. forward to pf2.probe.dot.reconstruct; 'Style', 'initCamPosition', 'cmap' forward to pf2.probe.project.tomography.

Fidelity and high-density montages

pf2.probe.dot.montageInfo characterizes a montage's DOT regime before you invest in reconstruction — sparse single-distance montages get coarse, topography-like images; high-density, overlapping multi-distance montages support genuine depth-resolved reconstruction:

info = pf2.probe.dot.montageInfo(proc, 'Print', true);
if info.isHighDensity
    recon = pf2.probe.dot.reconstruct(proc);
end

info reports .nChannels/.nSources/.nDetectors, .separations and .sepClasses (distinct SD-distance groups), .hasShortSep, .meanOverlap (sensitivity redundancy per covered vertex), .isHighDensity (heuristic: multi-distance AND meanOverlap > 2), and a plain-English .recommendation.

pf2.probe.dot.resolution quantifies the honest point-spread of the reconstruction operator — how well a focal cortical change can actually be localized — by probing seed vertices across the coverage support:

res = pf2.probe.dot.resolution(proc, 'NSeeds', 60);
fprintf('median localization %.1f mm, spread %.1f mm, FWHM %.1f mm\n', ...
    res.summary.medianLocalization, res.summary.medianSpread, res.summary.medianFWHM);

Pass 'Lambda' (e.g. a reconstruction's recon.lambda.HbO) to make the reported resolution match the operator a specific reconstruction actually used; otherwise it reflects 'LambdaFraction' (default 0.05), not a data-driven lambda. Reference: Arridge (1999), Optical tomography in medical imaging, Inverse Problems, 15(2), R41-R93. DOI: 10.1088/0266-5611/15/2/022.

For higher-fidelity reconstructions on high-density/multi-distance montages, combine a layered head model, short-separation scalp regression, and a spatial-smoothness prior:

reconHD = pf2.probe.dot.reconstruct(proc, 'Time', [5 20], ...
    'HeadModel', 'layered', 'ScalpRegression', true, 'Prior', 'laplacian');
pf2.probe.project.tomography(reconHD, 'Biomarker', 'HbO', ...
    'savePath', 'dot_HbO_layered_smooth.png');

'HeadModel','layered' uses a scalp/skull/CSF/gray-matter ray-segmented model instead of the default homogeneous semi-infinite medium. 'ScalpRegression' removes superficial signal via the montage's short-separation channels before inversion (no-op with a warning if the montage has none); 'ScalpMethod' selects 'nearest', 'all', or 'pca'. 'Prior','laplacian' favors spatially smooth images; 'Prior','parcel' penalizes within-Brodmann-area variance for a piecewise-per-region result.

Advanced internals

The functions above are built on lower-level primitives in pf2_base.+dot, useful when extending the forward/inverse model directly:

Function	Role
channelGeometry	Per-channel source/detector positions, midpoints, SD distances
corticalMesh	Loads the bundled MNI cortical surface (high-res/low-res)
greensFunction / greensFunctionLayered	Homogeneous / layered semi-infinite Green's function for light propagation
layeredHeadModel	Scalp/skull/CSF/gray-matter layer definitions for 'HeadModel','layered'
opticalProperties	Wavelength-dependent mua/musp defaults
sensitivityMatrix	Assembles the PMDF sensitivity matrix from geometry + Green's functions
meshLaplacian	Graph-Laplacian operator over the cortical mesh (for 'Prior','laplacian')
reconstructImage	The regularized (Tikhonov) inverse solver itself
regParam	Regularization parameter selection (GCV/L-curve)
surfaceNormals	Per-vertex/per-optode surface normals for 'NormalMode','surface'

See examples/scripts/example_dot_reconstruction.m for a complete run through the forward model, coverage, PMDF projection, reconstruction, rendering, and the fidelity/HD-DOT extras above.

See also

• Visualization
• Processing Pipeline
• Export and Interoperability
• API Reference