An L_{p} (0 ≤ p ≤ 1)norm regularized image reconstruction scheme for breast DOT with nonnegativeconstraint
 Bingyuan Wang^{1},
 Wenbo Wan^{1},
 Yihan Wang^{1},
 Wenjuan Ma^{2},
 Limin Zhang^{1},
 Jiao Li^{1},
 Zhongxing Zhou^{1},
 Huijuan Zhao^{1, 3}Email author and
 Feng Gao^{1, 3}Email author
DOI: 10.1186/s129380170318y
© The Author(s) 2017
Received: 20 September 2016
Accepted: 30 January 2017
Published: 3 March 2017
Abstract
Background
In diffuse optical tomography (DOT), the image reconstruction is often an illposed inverse problem, which is even more severe for breast DOT since there are considerably increasing unknowns to reconstruct with regard to the achievable number of measurements. One common way to address this illposedness is to introduce various regularization methods. There has been extensive research regarding constructing and optimizing objective functions. However, although these algorithms dramatically improved reconstruction images, few of them have designed an essentially differentiable objective function whose full gradient is easy to obtain to accelerate the optimization process.
Methods
This paper introduces a new kind of nonnegative prior information, designing differentiable objective functions for cases of L_{1}norm, L_{p} (0 < p < 1)norm and L_{0}norm. Incorporating this nonnegative prior information, it is easy to obtain the gradient of these differentiable objective functions, which is useful to guide the optimization process.
Results
Performance analyses are conducted using both numerical and phantom experiments. In terms of spatial resolution, quantitativeness, gray resolution and execution time, the proposed methods perform better than the conventional regularization methods without this nonnegative prior information.
Conclusions
The proposed methods improves the reconstruction images using the introduced nonnegative prior information. Furthermore, the nonnegative constraint facilitates the gradient computation, accelerating the minimization of the objective functions.
Keywords
Diffuse optical tomography Inverse problem Sparsity regularization NonnegativeBackground
Diffuse optical tomography (DOT) has attracted widespread interest in recent years due to its noninvasive [1, 2] and sensitive [3] properties, which offers huge clinical potential. This biomedical imaging modality presents lower radiation risk than those methods using Xray [3] and has extensive applications including optical mammography and functional brain imaging [1, 4]. DOT reconstructs the spatial distribution of the optical properties such as absorption coefficients, which provide useful functional information about blood volume and oxygenation. However, the reconstructed images of DOT often suffer from the low spatial resolution and are significantly affected by measurement noise because the inverse problem of DOT is highly illposed [1–5]. The inverse problem of DOT is both undetermined because the available measurements are far fewer than the unknown variables to be reconstructed and illconditioned because of the dominance of scattering during the propagation of light in human tissues.
To alleviate this illposedness, researchers introduced various useful prior information to regularize the DOT inverse problem. The first introduced prior information is smoothness constrain incorporated in the wellestablished Tikhonov regularization method, which to some extent improved the reconstruction image. The involved optimization problem can be effectively solved because its objective function is differentiable. However, the reconstructed image of this method is often oversmoothed, due to the incorporated smoothness constraint, discouraging sharp edges in the reconstructed images. The second important prior information, used to alleviate the illposedness, is that the real solution is sparse. Its corresponding regularization term is based on the L_{p}norm (0 ≤ p ≤ 1) instead of L_{2}norm. These sparsity regularization methods can facilitate the recovery of the sharp edges and are robust with noise [4, 5]. This is based on the clinical fact that the breast tumor usually accounts for a small part of the overall breast, and the remaining part is healthy [2, 4]. Thus, the changes of the optical coefficients, which can be caused by the regional blood blow changes or the early stage of breast cancers are expected to be localized. That is to say, the original reconstruction problem itself is sparse. Up to now, researchers have devised a variety of sparsity regularization methods using different L_{p} (0 ≤ p ≤ 1)norms [6, 7]. For example, in case of L_{1}norm, Amir Beck et al. proposed a practical method named FISTA [8] using the gradient information; Figueiredo et al. proposed a method called GPSR [9] to split the L_{1}norm into two parts, from which the gradient is easy to obtain. To use L_{p} (0 < p < 1)norm regularization, several approaches [4, 10, 11] were proposed to obtain the optimal solution with the help of majorizationminimization framework. When p = 0, a wellestablished smooth L_{0}norm based regularization method [4, 12] was designed, which had proven to be effective. The total variation methods of special sparsity regularization [13] also improved the reconstruction image by using the L_{1}norm of the difference of neighboring pixels rather than that of the vector in the reconstructed image as the constraint term.
However, image reconstruction using these algorithms is a rather lengthy process. Their objective functions are essentially nondifferentiable, unlike those of the Tikhonov method, making it impossible to use their gradient information directly to minimize them. More specifically, using these algorithms requires more computational work and time seeing as the full gradient information cannot be easily obtained to guide the optimization process.
In this paper, a new kind of nonnegative prior information is introduced for the first time, designing regularization methods with differentiable objective functions for L_{1}norm, L_{p} (0 < p < 1)norm and L_{0}norm. These regularization methods not only make full use of the gradient information of the objective functions, like the Tikhonov method, but also retain the advantages of the sparsity regularizations in improving image quality. To investigate the performances of these proposed methods, the methods with the nonnegative constraint and those without the nonnegative constraint are compared, using numerical simulation and phantom experiments.
Methods
The nonnegative prior information
To optimize Eq. (2) more efficiently and more precisely, a new kind of nonnegative prior information is incorporated to constrain the solution. Specially speaking, from the viewpoint of physiology, the abnormal absorption coefficients caused by the breast tumor with angiogenesis must have bigger absorption coefficients than those of the normal region because the tumor has more hemoglobin [5]. Thus the changes of the absorption coefficients compared with the normal region are nonnegative; that is to say, \(\delta {\varvec{\upmu}}_{a} \ge 0\) can be used as an inequality constraint. This is the nonnegative prior information, which will play important role in designing differentiable objective functions for sparsity regularization methods and reconstructing images of higher quality.
L_{1}norm regularized reconstruction scheme with the nonnegative constraint
The detailed gradient projection method for L_{1}norm regularized reconstruction scheme with nonnegative constraint is shown below, according to the optimization theory. Hereinafter, this method will be referred to as the NL_{1} (NonNegative L_{1}) method for short.
L_{p}norm regularized reconstruction scheme with the nonnegative constraint
In this scenario, the gradient of the objective function is \({\mathbf{g}} = {\mathbf{J}}^{T} *({\mathbf{J}}\delta {\varvec{\upmu}}_{a}  \delta {\varvec{\Gamma}}){ + }\lambda \sum\nolimits_{i = 1}^{{N_{n} }} {\frac{1}{{\delta {\varvec{\upmu}}_{ai}^{1  p} + C}}}\), where C is a small positive constant that provides stability by ensuring that the zero components of δ μ _{a} do not prohibit a nonzero estimate in the next step. Like the NL_{1} method algorithm, the NL_{p} (NonNegative L_{1}) algorithm is designed by replacing the gradient term with \({\mathbf{g}} = {\mathbf{J}}^{T} *({\mathbf{J}}\delta {\varvec{\upmu}}_{a}  \delta {\varvec{\Gamma}}){ + }\lambda \sum\nolimits_{i = 1}^{{N_{n} }} {\frac{1}{{\delta {\varvec{\upmu}}_{ai}^{1  p} + C}}}\). The remaining part of this algorithm is the same as that of the NL_{1} algorithm. In this paper, the values of p and C are fixed at 1/2 and 1e6, respectively.
L_{0}norm regularized reconstruction scheme with the nonnegative constraint
The whole algorithm for the NL_{0} (NonNegative L_{0}) is presented as follows using the pseudo code:
Numerical and phantom experiments
To evaluate the performances of the proposed regularization methods, several numerical and phantom experiments are conducted. Their reconstructed results are evaluated and then compared with those of the conventional regularization methods based on the same norm.
Conditions
The optical properties of healthy and diseased breast tissues have been extensively investigated using both in vitro and in vivo methods [14]. Most investigations have measured the absorption and (reduced) scattering coefficients of healthy breast tissues in the range of 0.002–0.008 mm^{−1} and 0.6–1.3 mm^{−1}, respectively, in the nearinfrared wavelength range. The in vivo measurements reported by Tromberg et al. indicate that some tumors exhibit 1.25 to threefold higher absorption than normal breast tissue [15]. Overall, most evidence suggests that the ratio between the absorption for healthy and diseased tissues are of the order of a factor of 2 [14]. In this paper, to simulate the clinical cases, the absorption and (reduced) scattering coefficients of the background medium are set to 0.004 and 1 mm^{−1}, respectively. For the target regions, the absorption coefficients vary with the purposes of the experiments while the scattering coefficients are the same as that of the background: for the spatial resolution, the absorption coefficient of the targets are set as 0.008 mm^{−1}; for the quantitativeness, both the two targets have the same absorption coefficients of 0.005, 0.006 mm, and 0.008 mm^{−1}, mimicking the tumors with contrasts of 1.25, 1.5 and 2.0, respectively. For the grayscale resolution, the two targets have different absorption coefficients, paired as (0.0072, 0.0088 mm^{−1}), (0.0064, 0.0096 mm^{−1}) and (0.0056, 0.0104 mm^{−1}), corresponding to a grayscale difference of 10, 20 and 30%. It is worth noting that the same experimental setup, including the geometry and optical coefficients, are also used in the phantom experiments, facilitating the mutual collation between the numerical and phantom experiments.
Thirtytwo source positions (red regions) evenly located one mean transport length inside the boundary are illuminated by the continuouswave source one by one. For each illumination, except the 15 detectors nearest to the source position, the remaining 17 detectors evenly located on the surface of the medium are used to measure the intensities of the reemitted light, leading to a total of 32 × 17 measurements.
When optimizing the regularization parameter λ, TE is calculated as a function of the regularization parameter, λ. The λ producing the smallest TE is selected because it provides the best balance across the above 3 metrics.
Generally speaking, better image reconstruction can be identified by: smaller RMSE, TE, and SR values; a larger CNR value; an AR value closer to 1; and AC and GR values closer to the true values.
Numerical experiments
To more closely simulate the real case, we assume that the main noises in the measurements are Gaussian noise [4–6]. Gaussian of different levels are added to the pure simulated measurements. The signaltonoise ratio (SNR) of the detecting location with the weakest light is set to SNR _{min} = 20 and 30 dB, respectively. The best signaltonoise ratio for a number of photons, N, reaching the detector in a given time interval is \(SNR = \sqrt N\) [18]. So, the SNRs for the other detecting locations are set according to \(SNR = SNR_{\rm{min} } \sqrt {I/I_{\rm{min} } }\), with I and I _{min} being the intensities of the current detecting locations and the one with the weakest light, respectively. Different meshes are adopted in the forward and inverse problems to avoid tricky problems. Specifically, a fine mesh with 12,290 nodes and 24,098 triangle elements and a coarse mesh with 4526 nodes and 8762 triangle elements are used for the forward and inverse calculations, respectively. To comprehensively analyze and compare the methods, we repeat each experiment 10 times and plot the average value and 95% confidence interval of the different metrics.
Spatial resolution
Furthermore, to more clearly and directly demonstrate the results, we plot the reconstructed images for the cases of CCS = 22 mm in Fig. 3b. These are the most difficult cases to separate the two nearest targets and the best cases to highlight the differences of spatial resolution abilities between different methods. The first and second columns present the images reconstructed using the nonnegative methods and the conventional methods without the nonnegative constraint, respectively. The third column contains the profiles of the reconstructed absorption coefficients along the Xaxis. The first, second, and third rows contain the results corresponding to the L_{1}, L_{1/2} and L_{0}norm, respectively. To facilitate the comparison between the methods with and without the nonnegative prior information, all the images are shown using the same colormap. The profiles for the same normbased regularization methods are plotted in the same subfigures. The left and right parts contain the metrics for SNR _{min} = 20 dB and SNR _{min} = 30 dB, respectively. All the following figures presenting the images are likewise presented. From Fig. 3b, it is clear that the targets reconstructed using the proposed methods have smaller sizes and are better separated than those of the conventional methods based on the same norm. At the same time, the results of the nonnegative regularization methods have less artifacts than those of the conventional methods. This validates the bigger CNR of the nonnegative regularization methods compared to those of the conventional methods. This is because the proposed methods better suppress the influences of the noise and of the difference in the meshes used in the forward and inverse problems. Figure 3a, b jointly demonstrate that the nonnegative regularization methods have better spatial resolution while retaining better RMSE and CNR.
Quantitativeness
Furthermore, to more clearly and directly demonstrate the results, we plot the reconstructed images for the cases of AC = 1.25 in Fig. 4b. These are the most difficult cases to distinguish the targets from the background, due to the low contrast, and the best cases to highlight the difference of quantitativeness between different methods. From Fig. 4b, we observe that the images of the nonnegative regularization methods have less artifacts and are better distinguished from the background than those of the conventional methods based on the same norm. Figure 4a, b jointly demonstrate that the nonnegative regularization methods have better quantitativeness abilities while retaining better RMSE and CNR.
Grayscale resolution
Furthermore, to more clearly and directly demonstrate the results, we plot the reconstructed images for the cases with a GR being 10% in Fig. 5b. These are the most difficult cases to reconstruct the small grayscale difference, and the best cases to highlight the difference of grayscale resolution between different methods. It is easy to find that targets reconstructed using the nonnegative methods have less artifacts and closer grayscale difference to the real values. Figure 5a, b jointly demonstrate that the nonnegative regularization methods have better grayscale resolution ability while retaining better RMSE and CNR.
Execution time
Phantom experiment
Along the boundary of cylindrical phantom, 32 optodes with a coaxial structure of source and detector are placed equidistantly from one another (as shown in Fig. 7b). A steady state laser operating at a wavelength of 675 nm stimulates the phantom boundary through a 32:1 mechanical optical switch. Meanwhile, the reemitted light on the phantom boundary is collected and sent to 4 photomultiplier tubes (PMT) through four 8:1 mechanical optical switches. Then the single electron responses from the PMTs are detected and accumulated by the field programmable gate array (FPGA). The total number accumulated in 250 ms is used as the intensity of the light on the detecting locations on the phantom boundary. It is worth noting that the light of the locations near the sources is so strong that it possibly induces the stack effect, making the measurement imprecise. Therefore, only the remaining opposite 17 measurements farthest away from the sources are used. A total of 32 × 17 = 544 measurements, equal to that in the numerical experiments, are then used to reconstruct the spatial distribution of absorption coefficients.
All phantom experiments are conducted according to the following process. First, the holes embedded in the phantom are filled with the mixed solution having absorption and (reduced) scattering coefficients of 0.004 and 1 mm^{−1}, identical to those of the healthy breast tissue. One of the optical fibers works as a light source and the remaining opposite 17 optical fibers are used to measure the light intensity reemitted from the surface of the phantom. Then, the remaining optical fibers work as the light source one by one, and the corresponding light intensities are measured. The mixed solution in the target holes is then replaced with the mixed solution with larger absorption coefficients, mimicking the breast tumors of different absorption contrasts. Repeating the measurement process, another 544 measurements are obtained and then are utilized in the reconstruction process.
Spatial resolution analysis
Furthermore, to more clearly and directly demonstrate the results, we plot the reconstructed images for the cases of CCS = 22 mm in Fig. 8b. These are the most difficult cases to separate the two nearest targets and the best cases to highlight the differences of spatial resolution abilities between different methods. Figure 8b tells us that the images of the nonnegative regularization methods have less artifacts and the two targets are better separated than those of the conventional methods. Figure 8a, b jointly demonstrate that the nonnegative regularization methods have better spatial resolution ability while retaining better RMSE and CNR.
Quantitativeness
Furthermore, to more clearly and directly demonstrate the results, we plot the reconstructed images for the cases of AC = 1.25 in Fig. 9b. These are the most difficult cases to distinguish the targets from the background, due to the low contrast, and the best cases to highlight the difference of quantitativeness between different methods. From Fig. 9b, we observe that the images of the nonnegative regularization methods have less artifacts and are better distinguished from the background than those of the conventional methods. Figure 9a, b jointly demonstrate that the nonnegative regularization methods have better quantitativeness ability while retaining better RMSE and CNR.
Grayscale resolution
Furthermore, to more clearly and directly demonstrate the results, we present the reconstructed images for the cases of a grayscale difference of 10% in Fig. 10b. These are the most difficult cases to reconstruct the small grayscale difference, and the best cases to highlight the difference of grayscale resolution between different methods. Figure 10b shows that the targets reconstructed using the nonnegative methods have less artifacts and closer grayscale difference to the real values than those of the conventional methods. Figure 10a, b jointly demonstrate that the nonnegative regularization methods have better grayscale resolution abilities while retaining better RMSE and CNR.
Execution time
The same stopping criteria and computer as those of the numerical simulation experiments are also adopted. Execution time of all the regularization methods for the phantom experiments analyzing spatial resolution for CCS = 22 mm are recorded and presented in Fig. 6c. It is obvious that the nonnegative regularization methods require less time to obtain the optical solution than those of the conventional regularization methods based on the same norm.
Discussion and conclusion
Though the investigations described in this paper are implemented on the continuous wave breast diffuse optical tomography, they can easily be extended to diffuse fluorescence tomography, because they share the same sparse and nonnegative prior information.
In conclusion, we propose several new regularized methods based on L_{1}norm, L_{p}norm (0 < p < 1) and approximate L_{0}norm with nonnegative prior information. Both numerical and phantom experiments demonstrate that, excepting AR, all metrics are improved by these proposed methods. Furthermore, the proposed methods require less time than those of the conventional methods based on the same norm.
Abbreviations
 DOT:

diffuse optical tomography
 CCS:

centertocenter separation
 RMSE:

root mean square error
 AR:

area ratio
 CNR:

contrasttonoise ratio
 TE:

total error
 SR:

spatial resolution
 AC:

absorption contrast
 GR:

grayscale resolution
 SNR _{min} :

the minimum signal noise ratio
 PMT:

photomultiplier tubes
 FPGA:

field programmable gate array
 NonNegL_{1} :

L_{1}norm regularized reconstruction scheme with nonnegative constraint
Declarations
Authors’ contributions
BW carried out the design of the algorithms and the experiments. WW and YW participated in the implementation of the experiments. WM, LZ, JL and ZZ participated in providing guidance for designing the algorithm. HZ and FG carried out the design of comparison of the different regularization methods. All authors read and approved the final manuscript.
Acknowledgements
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Data availability statement
Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.
Funding
The authors acknowledge the funding supports from the National Natural Science Foundation of China (81371602, 61475115, 61475116, 81401453, 61575140, 81571723, 81671728), and Tianjin Municipal Government of China (14JCQNJC14400, 15JCZDJC31800, 15JCQNJC14500, 16JCZDJC31200).
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Cao N, Nehorai A, Jacob M. Image reconstruction for diffuse optical tomography using sparsity regularization and expectationmaximization algorithm. Opt Express. 2007;15:13695–708.View ArticleGoogle Scholar
 Chen C, Tian F, Liu H, Huang J. Diffuse optical tomography enhanced by clustered sparsity for functional brain imaging. IEEE Trans Med Imaging. 2014;33:2323–31.View ArticleGoogle Scholar
 Lee O, Kim JM, Bresler Y, Ye JC. Compressive diffuse optical tomography: noniterative exact reconstruction using joint sparsity. IEEE Trans Med Imaging. 2011;30:1129–42.View ArticleGoogle Scholar
 Prakash J, Shaw CB, Manjappa R, Kanhirodan R, Yalavarthy PK. Sparse recovery methods hold promise for diffuse optical tomographic image reconstruction. IEEE J Sel Top Quantum Electron. 2014;20:6800609.View ArticleGoogle Scholar
 Okawa S, Yoko H, Yamada Y. Improvement of image quality of timedomain diffuse optical tomography with lp sparsity regularization. Biomed Opt Express. 2011;2:3334–48.View ArticleGoogle Scholar
 Han D, Tian J, Zhu S, Feng J, Qin C, Zhang B, Yang X. A fast reconstruction algorithm for fluorescence molecular tomography with sparsity regularization. Opt Express. 2010;18:8630–46.View ArticleGoogle Scholar
 Candès EJ, Wakin MB. An Introduction to Compressive Sampling. IEEE Signal Process Mag. 2008;2:21–30.View ArticleGoogle Scholar
 Beck A, Teboulle M. A fast iterative shrinkagethresholding algorithm for linear inverse problems. SIAM J Imaging Sci. 2009;2:183–202.MathSciNetView ArticleMATHGoogle Scholar
 Figueiredo MAT, Nowak RD, Wright SJ. Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems. IEEE J Sel Top Signal Process. 2007;1:586–97.View ArticleGoogle Scholar
 Guo H, Yu J, He X, Hou Y, Dong F, Zhang S. Improved sparse reconstruction for fluorescence molecular tomography with L_{1/2} regularization. Biomed Opt Express. 2015;6:1648–64.View ArticleGoogle Scholar
 Zhao L, Yang H, Cong W, Wang G, Intes X. L_{p} regularization for early gate fluorescence molecular tomography. Opt Lett. 2014;39:4156–9.View ArticleGoogle Scholar
 Mohimani H, Zadeh MB, Jutten C. A fast approach for overcomplete sparse decomposition based on smoothed ℓ_{0} norm. IEEE Trans Signal Process. 2009;57:289–301.MathSciNetView ArticleGoogle Scholar
 Li Y, Santosa F. A computational algorithm for minimizing total variation in image restoration. IEEE Trans Image Process. 1996;5:987–95.View ArticleGoogle Scholar
 Hebden JC, Veenstra H, Dehghani H, Hillman EMC, Schweiger M, Arridge SR, et al. Threedimensional timeresolved optical tomography of a conical breast phantom. Appl Opt. 2001;2001(40):3278–87.View ArticleGoogle Scholar
 Tromberg BJ, Shah N, Lanning R, et al. Noninvasive in vivo characterization of breast tumors using photon migration spectroscopy. Neoplasia. 2000;2:26–40.View ArticleGoogle Scholar
 Zhu D, Li C. Nonconvex regularizations in fluorescence molecular tomography for sparsity enhancement. Phys Med Biol. 2014;59:2901–12.View ArticleGoogle Scholar
 Gao F, Zhao H, Zhang L, Tanikawa Y, Marjono A, Yamada Y. A selfnormalized, full time resolved method for fluorescence diffuse optical tomography. Opt Express. 2008;16:13104–21.View ArticleGoogle Scholar
 Vojnovic B. Advanced timecorrelated single photon counting techniques. Berlin: Springer; 2006.Google Scholar
 Tehrani JN, Mcewan A, Jin C, et al. L_{1} regularization method in electrical impedance tomography by using the L_{1}curve (Pareto frontier curve). Appl Math Model. 2011;36:1095–105.MathSciNetView ArticleMATHGoogle Scholar
 Chen W, Wang X, Wang B, Wang Y, Zhang Y, Zhao H, Gao F. Lockinphotoncountingbased highlysensitive and largedynamic imaging system for continuouswave diffuse optical tomography. Biomed Opt Express. 2016;7:499–511.View ArticleGoogle Scholar
 Qin D, Ma Z, Gao F, Zhao H. Determination of optical properties in turbid medium based on time resolved determination. Proc SPIE. 2007;6534:65340T.View ArticleGoogle Scholar