Sparsity promoting reconstructions via hierarchical prior models in diffuse optical tomography

Anssi Manninen; Meghdoot Mozumder; Tanja Tarvainen; Andreas Hauptmann

doi:10.3934/ipi.2023025

Article Contents

2024, Volume 18, Issue 1: 113-137. Doi: 10.3934/ipi.2023025

This issue Previous Article A uniqueness theorem for inverse problems in quasilinear anisotropic media II Next Article Linear inverse problems with nonnegativity constraints: Singularity of optimisers

Sparsity promoting reconstructions via hierarchical prior models in diffuse optical tomography

1.
Research Unit of Mathematical Sciences, University of Oulu, Oulu 90014, Finland
2.
Department of Technical Physics, University of Eastern Finland, Kuopio 70211, Finland
3.
Department of Computer Science, University College London, London WC1E 6BT, United Kingdom

^*Corresponding author: Anssi Manninen
^*Corresponding author: Anssi Manninen

Received: July 2022

Revised: April 2023

Early access: May 10, 2023

Published online: May 10, 2023

Abstract / Introduction Full Text(HTML) Figure(7) / Table(3) Related Papers Cited by

Abstract

Diffuse optical tomography (DOT) is a severely ill-posed nonlinear inverse problem that seeks to estimate optical parameters from boundary measurements. In the Bayesian framework, the ill-posedness is diminished by incorporating a priori information of the optical parameters via the prior distribution. In case the target is sparse or sharp-edged, the common choice as the prior model are non-differentiable total variation and $ \ell^1 $ priors. Alternatively, one can hierarchically extend the variances of a Gaussian prior to obtain differentiable sparsity promoting priors. By doing this, the variances are treated as unknowns allowing the estimation to locate the discontinuities.

In this work, we formulate hierarchical prior models for the nonlinear DOT inverse problem using exponential, standard gamma and inverse-gamma hyperpriors. Depending on the hyperprior and the hyperparameters, the hierarchical models promote different levels of sparsity and smoothness. To compute the MAP estimates, the previously proposed alternating algorithm is adapted to work with the nonlinear model. We then propose an approach based on the cumulative distribution function of the hyperpriors to select the hyperparameters. We evaluate the performance of the hyperpriors with numerical simulations and show that the hierarchical models can improve the localization, contrast and edge sharpness of the reconstructions.

Keywords:

Mathematics Subject Classification: Primary: 65N21, 94A08; Secondary: 62F15.

Citation:

Full Text(HTML)

Figure 1. Computed MAP estimates. The two top-most rows show the true values and the estimates when the uncorrelated Gaussian prior (5) with fixed variances of $ 0.25^2 $ and $ 0.0025^2 $ for reduced scattering and absorption were used. Other rows show the estimates when variances were a) inverse-gamma, b) standard gamma, and c) exponentially distributed. The columns show the estimates with small (left), intermediate (mid) and large (right) hyperparameter $ \vartheta $ or $ \gamma $ values as listed in Table 1. The colorbars exclude some of the highest and smallest values.

Download: Full-size image PowerPoint slide

Figure 2. Reconstructions of smaller inclusions using uncorrelated Gaussian prior (5) with standard gamma hyperpriors. The leftmost column shows the used targets, where the radiuses of the inclusions were 2.5 mm, 1.1 mm (lower) and 2.2 mm, 0.8mm (upper). The remaining columns from left to right show the reconstructions with fixed variance values and variances with low, intermediate and high hyperparameter values. The fixed values and the hyperparameter values are the same as used in Figure 1 row b).

Download: Full-size image PowerPoint slide

Figure 3. Computed MAP estimates. The two top-most rows show the true values and the estimates when the difference prior (5) with fixed variances of $ 0.1^2 $ and $ 0.001^2 $ for reduced scattering and absorption were used. Other rows show the estimates, when variances were a) inverse-gamma, b) standard gamma and c) exponentially distributed. The columns show the estimates with small (left), intermediate (mid) and large (right) hyperparameter $ \vartheta $ or $ \gamma $ values as listed in Table 1. The colorbars exclude some of the highest and smallest values.

Download: Full-size image PowerPoint slide

Figure 4. Computed MAP estimates of alternative target type with only positive inclusions. The left-most column show the target. The second column shows the estimates when the difference prior (6) with fixed variances of $ 0.1^2 $ and $ 0.001^2 $ for reduced scattering and absorption was used. Other columns show the estimates, when variances were inverse-gamma, standard gamma and exponentially distributed. The hyperparameters $ \vartheta $ were set as the intermediate values, also used for recontructions in Figure 3.

Download: Full-size image PowerPoint slide

Figure 5. Computed inner (Gauss-Newton) iterations cumulatively as a function of the IAS iterations. The computed iterations were performed during the computation of the reconstruction with the uncorrelated Gaussian model as shown in Figure 1.

Download: Full-size image PowerPoint slide

Figure 6. Convergence of the MAP estimates of the hierarchical models. The vertical axis shows the error $ ||x^i-x_{\rm MAP}||_2 $ of the $ i $th iteration and the horizontal axis the error of the $ (i+1) $th iteration. The MAP estimates $ x_{\rm MAP} $ were computed by accurately solving the corresponding optimization problems. The bottom row shows the convergence of the difference prior model and the top row convergence of the uncorrelated Gaussian prior. The left, middle and right columns show the convergence for standard gamma, exponential and inverse-gamma hyperpriors, respectively. The red lines show the convergence of the used nonlinear IAS algorithm and the dashed blue line shows a linear convergence with convergence constant of $ \mu = 0.6 $.

Download: Full-size image PowerPoint slide

Figure 7. The absorption and reduced scattering estimates after three IAS iterations, when uncorrelated Gaussian prior was used with the intermediate hyperparameter values as in Figure 1.

Download: Full-size image PowerPoint slide

Table 1. Used hyperparameter ($ \vartheta $ and $ \gamma $) values for the hierarchical uncorrelated Gaussian prior (unc.) and structural difference prior (dif.). The hyperparameter $ \vartheta $ values were computed from the CDF (30). The hyperparameter values for reduced scattering (scat.) were computed with $ M\in\{0.3, \, 1, \, 10\} $ (unc.) and $ M\in\{1, \, 5, \, 10\} $ (dif.) for standard gamma and $ M\in\{0.3, \, 1, \, 5\} $ (unc.) and $ M\in\{0.25, \, 1, \, 4\} $ (dif.) for inverse-gamma. For absorption (abs.), the assumed $ M $ was 0.01 times the corresponding reduced scattering values.

		Low		Intermediate		High
		Unc.	Dif.	Unc.	Dif.	Unc.	Dif.
Exponential ($ \gamma) $	Scat.	$ 10^{-{10}} $	$ 9\cdot10^{-4} $	2.5$ \cdot 10^{-3} $	3.6$ \cdot 10^{-3} $	0.25	0.14
Exponential ($ \gamma) $	Abs.	$ 10^{-{14}} $	$ 9\cdot10^{-8} $	2.5$ \cdot 10^{-7} $	3.6$ \cdot 10^{-7} $	2.5$ \cdot10^{-3} $	1.4$ \cdot10^{-3} $
Standard gamma ($ \vartheta $)	Scat.	5.8$ \cdot 10^{-3} $	6.4$ \cdot 10^{-2} $	6.4$ \cdot 10^{-2} $	1.6	6.4	6.4
Standard gamma ($ \vartheta $)	Abs.	5.8$ \cdot 10^{-7} $	6.4$ \cdot 10^{-6} $	6.4$ \cdot 10^{-6} $	1.6 $ \cdot 10^{-4} $	6.4 $ \cdot 10^{-4} $	6.4 $ \cdot 10^{-4} $
Inverse-gamma ($ \vartheta $)	Scat.	$ 4\cdot10^{-3} $	5.5$ \cdot10^{-3} $	1.1$ \cdot 10^{-2} $	8.8$ \cdot 10^{-2} $	0.4	1.4
Inverse-gamma ($ \vartheta $)	Abs.	$ 4\cdot10^{-7} $	5.5$ \cdot10^{-7} $	1.1$ \cdot 10^{-6} $	8.8$ \cdot 10^{-6} $	4 $ \cdot 10^{-3} $	1.4 $ \cdot 10^{-4} $

| Show Table

DownLoad: CSV

Table 2. Relative errors (%) of the MAP estimates when the uncorrelated Gaussian prior (5) was used with low, intermediate and high hyperparameter $ \vartheta $ (or $ \gamma $) values. The left-most column indicates the type of the hyperprior used for the unknown variances. The non-hierarchical model with the fixed variances is denoted as "no hyperprior". For reduced scattering and absorption, the smallest relative error is bold.

	Low		Intermediate		High
	RE ($ \mu_{\rm{a}} $)	RE ($ \mu_{\rm{s}}^\prime $)	RE ($ \mu_{\rm{a}} $)	RE ($ \mu_{\rm{s}}^\prime $)	RE ($ \mu_{\rm{a}} $)	RE ($ \mu_{\rm{s}}^\prime $)
No hyperprior	7.35	6.80	7.35	6.80	7.35	6.80
Exponential	7.87	16.46	6.55	5.42	6.81	7.01
Standard gamma	6.63	6.51	6.29	5.81	5.57	5.28
Inverse-gamma	7.94	7.69	7.73	7.34	7.06	6.74

| Show Table

DownLoad: CSV

Table 3. Relative errors (%) of the MAP estimates when the difference prior (6) was used with low, intermediate and high hyperparameter $ \vartheta $ (or $ \gamma $) values. The left-most column indicates the type of the hyperprior used for the unknown variances. The non-hierarchical model with the fixed variances is denoted as "no hyperprior". For reduced scattering and absorption, the smallest relative error is bold.

	Low		Intermediate		High
	RE ($ \mu_{\rm{a}} $)	RE ($ \mu_{\rm{s}}^\prime $)	RE ($ \mu_{\rm{a}} $)	RE ($ \mu_{\rm{s}}^\prime $)	RE ($ \mu_{\rm{a}} $)	RE ($ \mu_{\rm{s}}^\prime $)
No hyperprior	9.78	7.17	9.78	7.17	9.78	7.17
Exponential	11.14	9.00	9.60	6.89	9.44	6.67
Standard gamma	10.55	6.54	8.74	4.28	7.80	4.40
Inverse-gamma	11.08	6.14	10.26	7.67	9.12	6.34

| Show Table

DownLoad: CSV

Related Papers

Cited by

References

[1]	C. Aguerrebere, A. Almansa, J. Delon, Y. Gousseau and P. Musé, A Bayesian hyperprior approach for joint image denoising and interpolation, with an application to hdr imaging, IEEE Transactions on Computational Imaging, 3 (2017), 633-646. doi: 10.1109/TCI.2017.2704439.
[2]	S. R. Arridge, Optical tomography in medical imaging, Inverse Problems, 15 (1999), R41-R93. doi: 10.1088/0266-5611/15/2/022.
[3]	S. R. Arridge and J. C. Schotland., Optical tomography: Forward and inverse problems, Inverse Problems, 25 (2009), 123010, 59 pp. doi: 10.1088/0266-5611/25/12/123010.
[4]	S. R. Arridge, M. Schweiger, M. Hiraoka and D. Delpy, Performance of an iterative reconstruction algorithm for near-infrared absorption and scatter imaging, Photon Migration and Imaging in Random Media and Tissues, International Society for Optics and Photonics, 1888 (1993), 360-371. doi: 10.1117/12.154654.
[5]	J. M. Bardsley, A. Seppanen, A. Solonen, H. Haario and J. P. Kaipio, Randomize-then-optimize for sampling and uncertainty quantification in electrical impedance tomography, SIAM/ASA Journal on Uncertainty Quantification, 3 (2015), 1136-1158. doi: 10.1137/140978272.
[6]	D. Calvetti, H. Hakula, S. Pursiainen and E. Somersalo, Conditionally Gaussian hypermodels for cerebral source localization, SIAM Journal on Imaging Sciences, 2 (2009), 879-909. doi: 10.1137/080723995.
[7]	D. Calvetti, A. Pascarella, F. Pitolli, E. Somersalo and B. Vantaggi, A hierarchical Krylov-Bayes iterative inverse solver for MEG with physiological preconditioning, Inverse Problems, 31 (2015), 125005, 23 pp. doi: 10.1088/0266-5611/31/12/125005.
[8]	D. Calvetti, M. Pragliola, E. Somersalo and A. Strang, Sparse reconstructions from few noisy data: Analysis of hierarchical {B}ayesian models with generalized gamma hyperpriors, Inverse Problems, 36 (2020), 025010, 29 pp. doi: 10.1088/1361-6420/ab4d92.
[9]	D. Calvetti, F. Sgallari and E. Somersalo, Image inpainting with structural bootstrap priors, Image and Vision Computing, 24 (2006), 782-793. doi: 10.1016/j.imavis.2006.01.015.
[10]	D. Calvetti and E. Somersalo, Local regularization and Bayesian hypermodels, Proc. SPIE 5910, Advanced Signal Processing Algorithms, Architectures, and Implementations XV, 5910 (2005).
[11]	D. Calvetti and E. Somersalo., A Gaussian hypermodel to recover blocky objects, Inverse Problems, 23 (2007), 733-754. doi: 10.1088/0266-5611/23/2/016.
[12]	D. Calvetti and E. Somersalo, Hypermodels in the Bayesian imaging framework, Inverse Problems, 24 (2008), 034013, 20 pp. doi: 10.1088/0266-5611/24/3/034013.
[13]	D. Calvetti and E. Somersalo, Inverse problems: From regularization to Bayesian inference, WIREs Computatinal Statistics, 10 (2018), e1427, 19 pp. doi: 10.1002/wics.1427.
[14]	D. Calvetti, E. Somersalo and A. Strang, Hierachical Bayesian models and sparsity: $\ell^2$-magic, Inverse Problems, 35 (2019), 035003, 26 pp. doi: 10.1088/1361-6420/aaf5ab.
[15]	J. Chen, Optical tomography in small animals with time-resolved Monte Carlo methods, Dissertation Abstracts International, 74 (2012).
[16]	S. G. Diamond, T. J. Huppert, V. Kolehmainen, M. A. Franceschini, J. P. Kaipio, S. R. Arridge and D. A. Boas, Dynamic physiological modeling for functional diffuse optical tomography, NeuroImage, 30 (2006), 88-101. doi: 10.1016/j.neuroimage.2005.09.016.
[17]	T. Durduran, R. Choe, W. B. Baker and A. G. Yodh, Diffuse optics for tissue monitoring and tomography, Reports on Progress in Physics, 73 (2015), 076701. doi: 10.1088/0034-4885/73/7/076701.
[18]	A. Gibson, J. Hebden and S. R. Arridge, Recent advances in diffuse optical imaging, Phys. Med. Biol., 50 (2005), R1. doi: 10.1088/0031-9155/50/4/R01.
[19]	D. Grosenick, H. Rinneberg, R Cubeddu and P. Taroni, Review of optical breast imaging and spectroscopy, Journal of Biomedical Optics, 21 (2016), 091311. doi: 10.1117/1.JBO.21.9.091311.
[20]	M. Guven, B. Yazici, X. Intes and B. Chance, Diffuse optical tomography with a priori anatomical information, Physics in Medicine & Biology, 50 (2005), 2837. doi: 10.1117/12.479799.
[21]	M. Guven, B. Yazici, X. Intes and B. Chance, Hierarchical Bayesian algorithm for diffuse optical tomography, 34th Applied Imagery and Pattern Recognition Workshop (AIPR'05), (2005), 6 pp. doi: 10.1109/AIPR.2005.30.
[22]	J. C. Hebden, A. Gibson, R. Md. Yusof, N. Everdell, E. M. C. Hillman, D. T. Delpy, S. R. Arridge, T. Austin, J. H. Meek and J. S. Wyatt, Three-dimensional optical tomography of the premature infant brain, Physics in Medicine & Biology, 47 (2002), 4155. doi: 10.1088/0031-9155/47/23/303.
[23]	P. Hiltunen, S. Prince and S. R. Arridge, A combined reconstruction-classification method for diffuse optical tomography, Physics in Medicine & Biology, 54 (2009), 6457. doi: 10.1088/0031-9155/54/21/002.
[24]	Y. Hoshi and Y. Yamada, Overview of diffuse optical tomography and its clinical applications, Journal of Biomedical Optics, 21 (2016), 091312. doi: 10.1117/1.JBO.21.9.091312.
[25]	D. Isaacson, Distinguishability of conductivities by electric current computed tomography, IEEE Transactions on Medical Imaging, 5 (1986), 91-95. doi: 10.1109/TMI.1986.4307752.
[26]	A. Ishimaru, Wave Propagation and Scattering in Random Media, IEEE/OUP Series on Electromagnetic Wave Theory, An IEEE/OUP Classic Reissue, IEEE Press, New York, 1997.
[27]	J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems, Applied Mathematical Sciences, 160. Springer-Verlag, New York, 2005.
[28]	F. Lin, T. Witzel, S. P. Ahlfors, S. M. Stufflebeam, J. W. Belliveau and M. S. H$\ddot{\text{a}}$m$\ddot{\text{a}}$l$\ddot{\text{a}}$inen, Assessing and improving the spatial accuracy in MEG source localization by depth-weighted minimum-norm estimates, NeuroImage, 31 (2006), 160-171. doi: 10.1016/j.neuroimage.2005.11.054.
[29]	S. Liu, J. Jia, D. Zhang, Y. Yimi and Y. Yang, Image reconstruction in electrical impedance tomography based on structure-aware sparse Bayesian learning, IEEE Transactions on Medical Imaging, 37 (2018), 2090-2102. doi: 10.1109/TMI.2018.2816739.
[30]	J. Mattout, C. Phillips, W. Penny, M. Rugg and K. Friston, MEG source localization under multiple constraints: An extended Bayesian framework, NeuroImage, 30 (2006), 753-767. doi: 10.1007/978-4-431-30962-8_33.
[31]	A. Miyamoto, K. Watanabe, K. Ikeda and M. Sato, Phase diagrams of a variational Bayesian approach with ARD prior in NIRS-DOT, The 2011 International Joint Conference on Neural Networks, (2011), 1230-1236. doi: 10.1109/IJCNN.2011.6033364.
[32]	M. Mozumder, A. Hauptmann, I. Nissilä, S. R. Arridge and T. Tarvainen, A model-based iterative learning approach for diffuse optical tomography, IEEE Transactions on Medical Imaging, 41 (2021), 1289-1299. doi: 10.1109/TMI.2021.3136461.
[33]	M. Mozumder, T. Tarvainen, J.P. Kaipio, S. R. Arridge and V. Kolehmainen, Compensation of modeling errors due to unknown domain boundary in diffuse optical tomography, JOSA A, 31 (2014), 1847-1855. doi: 10.1364/JOSAA.31.001847.
[34]	A. Nummenmaa, T. Auranen, M. Hämäläinen, I. Jääskeläinen, J. Lampinen, M. Sams and A. Vehtari, Hierarchical Bayesian estimates of distributed MEG sources: Theoretical aspects and comparison of variational and MCMC methods, NeuroImage, 35 (2007), 669-685. doi: 10.1016/j.neuroimage.2006.05.001.
[35]	K. Paulsen and H. Jiang, Enhanced frequency-domain optical image reconstruction in tissues through total-variation minimization, Applied Optics, 35 (1996), 3447-3458. doi: 10.1364/AO.35.003447.
[36]	P. Perona and J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Transactions on Pattern Analysis and Machine Intelligence, 12 (1990), 629-639. doi: 10.1109/34.56205.
[37]	L. Roininen, M. Girolami, S. Lasanen and M. Markkanen, Hyperpriors for matérn fields with applications in Bayesian inversion, Inverse Probl. Imaging, 13 (2019), 1-29, arXiv: 1612.02989. doi: 10.3934/ipi.2019001.
[38]	M. Schweiger and S. R. Arridge, The Toast++ software suite for forward and inverse modeling in optical tomography, Journal of Biomedical Optics, 19 (2014), 040801. doi: 10.1117/1.JBO.19.4.040801.
[39]	M. Schweiger, S. R. Arridge and I. Nissilä, Gauss-Newton method for image reconstruction in diffuse optical tomography, Physics in Medicine & Biology, 50 (2005), 2365. doi: 10.1088/0031-9155/50/10/013.
[40]	C. B. Shaw and P. K. Yalavarthy, Performance evaluation of typical approximation algorithms for nonconvex $\ell_p$-minimization in diffuse optical tomography, Journal of the Optical Society of America A, 31 (2014), 852-862. doi: 10.1364/JOSAA.31.000852.
[41]	T. Shimokawa, T. Kosaka, O. Yamashita, N. Hiroe, T. Amita, Y. Inoue and M. Sato, Hierarchical Bayesian estimation improves depth accuracy and spatial resolution of diffuse optical tomography, Optics Express, 20 (2012), 20427-20446. doi: 10.1364/OE.20.020427.
[42]	T. Shimokawa, T. Kosaka, O. Yamashita, N. Hiroe, T. Amita, Y. Inoue and M. Sato, Extended hierarchical Bayesian diffuse optical tomography for removing scalp artifact, Biomedical Optics Express, 4 (2013), 2411-2432. doi: 10.1364/BOE.4.002411.
[43]	D. Vidaurre, C. Bielza and P. Larra naga, A survey of L1 regression, International Statistical Review, 81 (2013), 361-387. doi: 10.1111/insr.12023.
[44]	M. Wheelock, J. Culver and A. Eggebrecht, High-density diffuse optical tomography for imaging human brain function, Review of Scientific Instruments, 90 (2019), 051101. doi: 10.1063/1.5086809.
[45]	O. Yamashita, T. Shimokawa, R. Aisu, T. Amita, Y. Inoue and M. Sato, Multi-subject and multi-task experimental validation of the hierarchical Bayesian diffuse optical tomography algorithm, NeuroImage, 135 (2016), 287-299. doi: 10.1016/j.neuroimage.2016.04.068.
[46]	J. Yoo, S. Sabir, D. Heo, K. H. Kim, A. Wahab, Y. Choi, S. Lee, E. Y. Chae, H. H. Kim, Y. M. Bae, Y.-W. Choi, S. Cho and J. C. Ye, Deep learning diffuse optical tomography, IEEE Transactions on Medical Imaging, 39 (2019), 877-887. doi: 10.1109/TMI.2019.2936522.
[47]	G. Zhang, X. Cao, B. Zhang, F. Liu, J. Luo and J. Bai, MAP estimation with structural priors for fluorescence molecular tomography, Physics in Medicine & Biology, 58 (2012), 351. doi: 10.1088/0031-9155/58/2/351.