A bilinear inverse problem with forward operator inaccuracy applied to neonatal atlas-based diffuse optical tomography

1. Problem Statement

The paper addresses a critical challenge in Diffuse Optical Tomography (DOT) for neonates: the inaccuracy of the forward operator caused by unknown or approximated subject anatomy.

• Context: DOT is a functional neuroimaging technique used to monitor hemodynamic changes (HbO2, HbR, HbT) in the brain. Accurate reconstruction requires a precise model of the head's anatomy (tissue segmentation, geometry) to compute the sensitivity profile (Jacobian matrix).
• The Issue: In neonatal imaging, acquiring individual MRI scans is often impractical. Instead, atlas-based models (average or population-level anatomies) are used. However, linear registration of an atlas to a specific neonate's head surface introduces inaccuracies in the internal anatomy (e.g., sulci location, tissue thickness), leading to errors in the Jacobian matrix.
• Mathematical Formulation: Standard approaches substitute a "mean" forward operator ($A_0$) into the linear inverse problem $Ax=b$. This ignores the statistical variation of the operator. The authors propose treating the forward operator $A$ as an unknown variable dependent on expansion coefficients, transforming the problem from a linear inverse problem into a bilinear inverse problem:
  $$A \approx A_0 + \sum_{i=1}^p y_i V_i$$
  where $V_i$ are basis matrices (derived from Principal Component Analysis) and $y_i$ are unknown coefficients. The goal is to simultaneously solve for the target image $x$ (absorption changes) and the operator coefficients $y$.

2. Methodology

2.1. Bayesian Framework and Bilinear Formulation

The authors formulate the problem within a Bayesian inversion framework.

• Model: The measurement equation is modeled as $b = A_0 x + A \cdot (y, x) + E$, where $E$ is noise.
• Priors: Mutually independent Gaussian priors are assumed for the noise ($E$), the operator coefficients ($y$), and the target ($x$).
• Objective: The goal is to find the Maximum A Posteriori (MAP) estimate or the Conditional Mean (CM) estimate of the joint posterior distribution $\pi(x, y | b)$.
• Challenge: The resulting objective function (negative log-posterior) is non-convex and possesses multiple local minima, making global optimization difficult.

2.2. Algorithms for Solution

Three algorithms are proposed to approximate the solution:

1. Block Coordinate Descent (BCD): Alternates between minimizing the objective with respect to $x$ (fixing $y$) and $y$ (fixing $x$). Both sub-problems are linear least squares solvable via the Woodbury matrix identity.
2. Gauss–Newton (G-N) Method: Linearizes the bilinear term around the current iterate and solves the resulting linear least squares problem. This method utilizes the low-rank structure of the Jacobian to remain computationally efficient.
3. Gibbs Sampling: A Markov Chain Monte Carlo (MCMC) method that samples from the joint posterior distribution to estimate the Conditional Mean (CM). This explores the full posterior landscape rather than seeking a single local minimum.

2.3. Application to Neonatal DOT

• Data Source: 215 segmented neonatal head models from the UCL Human Connectome Project (gestational ages 29–44 weeks).
• Jacobian Atlas: Monte Carlo simulations (using MCX software) generated 10 billion photon packets per model to create a database of Jacobians.
• PCA: Principal Component Analysis was performed row-wise on the Jacobian matrices (sensitivity profiles of individual source-detector pairs) to capture anatomical variations. The first 10 principal components were used to model the operator uncertainty.
• Registration: A linear registration (rotation, translation, scaling) aligned the atlas models to a target neonate's surface landmarks.
• Simulation: Synthetic hemodynamic perturbations (changes in absorption coefficients) were introduced in brain and scalp regions to generate "ground truth" data for testing.

3. Key Contributions

1. Bilinear Tensor Formulation: The paper introduces a novel approach to handle forward operator inaccuracy by embedding the operator uncertainty directly into the inverse problem as a bilinear tensor equation, rather than treating it as a fixed error or using simple regularization.
2. PCA-Based Operator Modeling: Instead of relying on a single "best-fit" atlas, the method utilizes a statistical basis (PCA) of Jacobian variations to allow the reconstruction algorithm to "tune" the forward operator during the solution process.
3. Algorithmic Comparison: The study provides a comprehensive comparison of optimization (BCD, G-N) and sampling (Gibbs) methods for this specific class of bilinear problems in medical imaging.
4. Real-World Validation: The methods are tested on a realistic, high-dimensional neonatal DOT setup, demonstrating improvements over standard mean-Jacobian approaches.

4. Results

The algorithms were tested on simulated data with various perturbation patterns and neonatal ages.

• Performance Metrics: The study evaluated Contrast-to-Noise Ratio (CNR) and Root-Mean-Squared Error (RMSE).
• Comparison with Mean Jacobian:
  • Reconstructions using the standard mean Jacobian suffered from poor spatial localization, severe background oscillations, and low CNR.
  • All three proposed algorithms (G-N, BCD, Gibbs) significantly improved CNR and reduced RMSE compared to the mean Jacobian.
• Algorithm Specifics:
  • Gauss–Newton (G-N): Converged quickly (seconds on a laptop) and provided substantial improvements in localization and contrast. It effectively reduced artifacts caused by anatomical mismodeling, though it did not perfectly recover the true operator coefficients ($y$).
  • Block Coordinate Descent (BCD): Produced similar qualitative results to G-N but required significantly more iterations (slower convergence).
  • Gibbs Sampler: Produced the best visual and numerical results (highest CNR, lowest error), closely matching the "true Jacobian" reconstruction. However, it was computationally expensive (orders of magnitude slower than G-N) and required a large number of samples.
• Local Minima: The study confirmed the non-convex nature of the problem. While the algorithms did not always recover the exact true operator coefficients (due to local minima), the resulting image reconstructions were superior to those obtained with a fixed mean operator.

5. Significance and Conclusion

• Clinical Impact: This work offers a practical solution for neonatal DOT where individual anatomical models are unavailable. By statistically modeling the uncertainty of atlas-based head models, the method enhances the reliability of functional brain imaging in premature and full-term infants without requiring invasive or time-consuming MRI scans for every subject.
• Methodological Insight: The paper demonstrates that even if the exact forward operator cannot be recovered (due to the ill-posedness of the bilinear problem), optimizing the operator coefficients simultaneously with the image reconstruction yields significantly better diagnostic images.
• Future Directions: The authors suggest that while Gibbs sampling provides the best results, the computational cost limits its real-time use. Future work may focus on hybrid approaches (e.g., Expectation-Maximization) or utilizing the marginalized posterior to improve the efficiency of the G-N method.

In summary, the paper successfully bridges the gap between theoretical bilinear inverse problems and practical medical imaging, showing that accounting for forward operator inaccuracy via a PCA-based bilinear framework significantly improves the quality of neonatal brain imaging.