1.

Introduction

Interest has been growing rapidly in various imaging modalities of optical tomography (OT) since computed tomography (CT) apparatus was first introduced in the 1970s. OT systems using nonradioactive sources can be built at relatively low cost. Diffuse optical tomography (DOT) providing functional information related to tissues has drawn great attention for the last two decades.^{1, 2, 3, 4, 5, 6} Developed techniques include mainly estimating in the near-infrared [(NIR), $650–950\mathrm{nm}$ ] the distribution of optical properties and their changes within a tissue volume, and then relating to spatial variations of physiological parameters, such as hemoglobin concentration and oxygen saturation. A wide range of applications has been explored, such as brain, breast, forearm, and hand joint, etc.^{7, 8, 9, 10, 11, 12, 13, 14} However, the NIR DOT imaging techniques suffer from low spatial resolution owing to the diffusive nature of scattered light. Additionally, optical property image reconstruction through inverse problems is computed off line through an iterative process that is usually computationally expensive.

To date, various reconstruction algorithms^{15, 16, 17, 18, 19} have been developed and evaluated by using considerable experimental data acquired from different hardware technology, including continuous wave (cw), frequency domain, and time domain. Although a large number of algorithms have been developed, a keen need to explore and improve optical imaging performance still remains. Thus, further studies are being conducted by employing a priori information from other modalities or implementing an inverse solution regularized with spatial or spectral constraints.

With the use of a priori information, a reconstructed optical property image can be improved. Some studies^{20, 21, 22} used gradient-based iterative image reconstruction schemes consisting of the minimization of an appropriately defined objective function separated into both least squares of errors and additional penalty terms containing a priori information. This gradient-based iterative image reconstruction method uses the gradient of the objective function in a line-minimization scheme to provide subsequent guesses of the spatial distribution of the optical properties for the forward model, and the reconstruction of these properties is completed once a minimum of this objective function is found. With the regions of interest (ROI) a priori known, a model-based iterative reconstruction procedure²³ was used to reconstruct the optical properties of the region. Also, we find that the reconstruction of a priori specified ROI converged faster than without prior information. Furthermore, one potential method²⁴ improved resolution and contrast consisting in the application of a spatially variant regularization parameter during the NIR DOT image reconstruction. Other approaches proposed a regularization parameter $\left(\lambda \right)$ selection method, such as using the generalized cross-validation or the quasi-optimality criterion,²⁵ the $f$ slope plotted in the curve of the solution norm versus $\ln (1∕\lambda )$ ,²⁶ and fixed noise figures defined as the ratio between the signal-to-noise ratios (SNRs) in measurements, and the SNRs in the image or the blur radius defined as a measure of the resolution.²⁷ The value-preserved images with adopting high-pass filtering were performed by the derivation of the Poisson maximum a posteriori superresolution algorithm;²⁸ that previous work enhancing image resolution was evaluated with quantitative measures for all tested cases.

Besides incorporating a priori information into the inverse solution, a common way to implement regularization in optical tomography is using Tikhonov regularization. Then, the inverse problem is formulated as a minimization of the mismatch of the measurements against the model and the regularization function. Usually, differential operators, such as the Laplacian- or Gaussian-type are used as regularization functions, leading to smooth solutions.^{29, 30, 31} Assuming the breast is segmented into distinct glandular and adipose regions, the overall chromophore concentrations and the inhomogeneities were estimated using the Gauss-Newton optimization approach and a linear perturbation algorithm, respectively.³²

Recently, there have been significant breakthroughs in the development of algorithms that incorporate multispectral constraints into the image-reconstruction process.^{33, 34, 35, 36} These studies have demonstrated the potential benefits of this spectrally constrained inverse resolution method. The advantages of using direct spectral reconstruction include distinct estimations of oxygen saturation, water fraction, and scattering parameters. Spectral priors generated superior quantification of all estimated NIR optical parameters compared to only using spatial priors.³⁷ A general framework for incorporating a single or multiple priors in diffuse OT was studied as well.³⁸ This was performed by both utilizing spatial and spectral priors in the context of imaging breast cancer.

As mentioned before, the reconstruction of optical-property images with an iteration procedure is usually computed off line and computationally expensive. The above studies, however, focused mainly on improving the spatial and spatial-frequency resolutions of reconstructed optical-property images with the introduction of the regularization term or varied priors in the procedure. If real-time resolution is required, then dedicated reconstruction hardwares or specialized computers are mandatory. Moreover, fast reconstruction algorithms should also be considered to reduce the computation load. It is worth emphasizing that our method can reduce computation time with the regularization term, which is designed on the viewpoint of the update characteristics in the iteration procedure but not utilizing any spatial/spectral a priori knowledge or constraints; some results can be found in the previous work.³⁹ In this paper, we show how to speed up the computation to find an inverse solution by using regularization with an iteration domain technique; additionally, the convergence rate of proposed soft priors superior to that of the primary method is illustrated.

2.

Theory

Image-reconstruction tasks contain forward modeling and inverse problem resolution. The forward computation consists of obtaining the intensity out of a subject under investigation for a given source, and the initial guess (or iterated result) on scattering and absorption coefficients. The inverse computation is to compute the scattering and absorption coefficients for a known light source and measured intensities in an iterative manner. The modeling described here is using light propagation through diffusive media via the diffusion approximation to the Boltzmann equation. The image reconstruction is sought from boundary value measurements based on the diffusion equation and the finite element method.

Because we utilize cw light illumination or dc data, the physical process of NIR light illuminating through a highly scattering medium can be approximated by the steady state diffusion equation,

With an approximation to applying the Newton-Raphson method and ignoring higher order terms, we obtain

It is known that to solve Eq. 4 is an ill-posed problem. Tikhonov regularization is a method stabilizing the inverse problem through incorporating a priori assumptions to constrain the desired solution. It is able to convert an ill-posed problem into a well-posed one and, further, to improve an ill-conditioned problem. The regularization term (penalty term) introduced in the process regularizes the problem and makes the update stable. It also strengthens the robustness of algorithm to noisy data with the adequate design of the regularization term. Generally, Tikhonov regularization is to optimize this ill-conditioned problem as

Following the general derivation above, further discussion is divided into Sections 2.1, 2.2. Most usual concerns for the constraints on the spatial domain are described in Section 2.1; subsequently, the novel viewpoint concerning the constraint on the iteration domain is delivered in Section 2.2.

2.2.

Constraints on the Iteration Domain

In NIR DOT, it is also crucial to accelerate the computation. But, to date, speeding up the computation in the iteration domain has not yet been explored. Here, we consider this issue through the use of a Lorentzian distributed function taking a natural logarithm computation as a constraint, i.e.,

Eq. 13

$$\Psi \left(\Delta \chi \right)=\sum _{p=1}^{K+L}\ln \frac{\gamma ∕\pi }{{\left(\Delta \chi \right)}_p^2+{\gamma }^2},$$

where $p$ is the calculated nodes in the subject under investigation and $\gamma$ is a user-defined positive parameter. As can be seen, $\Psi \left(\Delta \chi \right)⩽{\sum }_p\ln (1∕\pi \gamma )$ , $\forall \Delta \chi$ , meets the requirement of Eq. 5. Performing the differentiation indicated in Eq. 9, we can obtain the solution in an iterative formality

Eq. 14

$${\left(\Delta \chi \right)}_n={(J^{\mathrm{T}}J+\frac{\lambda I}{{\left(\Delta \chi \right)}_{n-1}^2+{\gamma }_n^{2}I})}^{-1}{J}^{\mathrm{T}}\Delta \Phi ,$$

and some parameters are chosen as

$$\lambda =0.75\max \left\{J^{T}J\right\},{\left(\Delta \chi \right)}_0^2=I,{\gamma }_n=\lambda \mathit{\text{or}}2.5e^{-2n},$$

where the subscript $n$ is the $n$ ’th iteration, “max” means the maximum value, and the superscript T denotes a transposition operation. One way to improve the convergence rate is using ${\gamma }_n=\lambda$ as the type 1 soft prior and using ${\gamma }_n=2.5e^{-2n}$ , an exponentially decreasing form, as the type 2 soft prior, where Type 1 is a parameter related to the system function (Jacobian matrix) and Type 2 is a user-defined parameter. Both values of ${\gamma }_n$ have been respectively employed to seek an inverse solution for comparison. It is noted that the minus sign in Eq. 6, the objective function, corresponds to the regularization term proposed here, because the term is constrained to an energy bound. For further inspection in Eqs. 13, 14, as known, ${\mu }_{\mathrm{a}}$ and $D$ are generally searched in a range of $[{10}^{-3}{\mathrm{mm}}^{-1}–{10}^{-1}\mathrm{mm}]$ and thus, $\Delta \chi$ is much smaller than a unit. It can be proven that even the use of the natural logarithm in the constraint $\Psi \left(\Delta \chi \right)$ still makes it a positive and finite value. The other reason to use “ln” is because the regularization term in Eq. 14 still remains in a form of the Lorentzian distributed function derived from the constraint associated with the Lorentzian distributed function in Eq. 13.

The Lorentzian distributed function, as depicted in Fig. 1 , is employed here owing to its following two characteristics:

• 1. Lorentzian distributed function has a sharp peak with a long tail, describing the histogram distribution of $\Delta \chi$ , many of $\Delta \chi$ $(\sim 0)$ at its peak and a small rest of $\Delta \chi$ distributing along its long tail.

• 2. Its histogram distribution can be further tuned with the parameter $\left(\gamma \right)$ as iteration increasing. Related to the consideration in convergence, the updated quantity, $\Delta \chi$ , decreases, ranging from the peak to the tail, as the iteration increases; whereas it has a smooth distribution in the beginning stage of iteration.

Additionally, as the shape of the histogram would be affected, it is smooth with a big value of $\gamma$ and sharp with a small value of $\gamma$ . Thus, Lorentzian distributed function can characterize the nature of $\Delta \chi$ in the iterative process as the distribution from a smooth to a sharp distribution to be used as a constraint for the purpose of speeding up computation.

Fig. 1

Charts of the Lorentzian distributed functions $[(\gamma ∕\pi )∕{\left(\Delta \chi \right)}^2+{\gamma }^2]$ at various $\gamma$ . As can be seen, it has a smooth distribution for a big $\gamma$ and a sharp distribution as $\gamma$ is small.

3.

Results and Discussion

To validate the efficiency and effectiveness of the proposed method, numerical simulation was performed under the condition of measured intensity with or without noise. A wide range of phantoms was tested to investigate the convergence rate. The test phantoms include a circular background ( $D=100\mathrm{mm}$ in diam) associated with varied numbers, sizes, and locations of inclusions, where the homogeneous background ( ${\mu }_{\mathrm{a}}=0.025{\mathrm{m}}^{-1}$ , ${\mu }_{\mathrm{s}}^{\prime }=2{\mathrm{mm}}^{-1}$ ) as a normal tissue was inserted with inclusions as “tumors” ( ${\mu }_{\mathrm{a}}=0.1{\mathrm{mm}}^{-1}$ , ${\mu }_{\mathrm{s}}^{\prime }=8{\mathrm{mm}}^{-1}$ ). Of the various test examples, two are chosen to be as representative as possible, where the phantoms containing two small and three big inclusions are demonstrated in Fig. 2 . These two examples have totally unlike features, resulting in different response to noise. The homogeneous background was adopted as an initial guess during the computation. The following results using the proposed and primary [obtained through Eq. 11] methods are shown by using 6 and 30 iterations, respectively, to justify the scheme on rapid convergence. It should be noted that no significant differences using more than 6 iterations of computation through the soft priors, and more than 30 iterations through the primary. To compare the results in details, as shown in Fig. 2 we use 1-D circular transections passing through the centers of inclusions from $0\text{to}360\mathrm{deg}$ .

Fig. 2

Schematic diagram of test phantoms: (a) Phantom with two small inclusions ( $15\mathrm{mm}$ diam) centered at $(35,135\mathrm{deg})$ and $(35,180\mathrm{deg})$ and (b) phantom with three big inclusions ( $30\mathrm{mm}$ diam) centered at $(25,90\mathrm{deg})$ , $(25,180\mathrm{deg})$ and $(25,270\mathrm{deg})$ , where the phantom ( $100\mathrm{mm}$ diam) is centered at (0,0). The dashed-line circles indicate 1-D circular transections for comparison among varied methods.

To justify the proposed method robust to noisy data, the measured intensity signals through the forward process were added with Gaussian white noise, which is defined by a zero mean and a variance $\left({\sigma }^2\right)$ . Here, the measured (or detected) SNR is specified by

Figure 3 demonstrates the example of image reconstruction performed on a phantom with two inclusions, where the reconstructed absorption images were obtained respectively from the primary inverse solution [Eq. 11], as shown in Figs. 3, 3, 3 and the constrained inverse solutions, Eq. 14, with types 1 and 2 soft priors, as shown in Figs. 3, 3, 3, 3, 3, 3, respectively. It is noted that there are a few differences among these reconstructed images through visual perception, where more fluctuations exist in the reconstructed absorption images in Figs. 3, 3, 3.

Fig. 3

Reconstructed absorption images for the phantom with two small inclusions, using (a–c) the primary inverse solution, (d–f) the constrained inverse solution (type 1 soft prior), and (g–i) the constrained inverse solution (type 2 soft prior). Further, left column is noise free, and middle and right columns are with Gaussian white noise (30- and $25\text{-}\mathrm{dB}$ SNRs, respectively).

As shown in Fig. 4 , the 1-D circular transection profile of each optical-property image is plotted for the inspection in detail. Obviously, the reconstructed 1-D circular transection profiles [Figs. 4, 4, 4, 4, 4, 4], using the constrained inverse solutions have better qualities showing two reconstructed peaks with higher contrast than those [Fig. 4, 4, 4], using the primary inverse solution. As to further inspect the influence of noise on the reconstructed, the cases with Gaussian white noise (30- and $25\text{-}\mathrm{dB}$ SNRs), Figs. 4 and 4 shows distorted reconstruction values using the primary method, whereas Figs. 4, 4, 4, 4 still show good reconstruction values using the constrained inverse solutions.

Fig. 4

1-D circular transection profiles corresponding to Fig. 3, where dotted lines are the original and solid lines are the reconstructed.

Concerning on the phantom with three big inclusions, as shown in Fig. 5 there are few differences among these reconstructed absorption images through visual perception, whether they were obtained either from the primary inverse solution or from the constrained inverse solutions. But, numerical figures via 1-D circular transection profiles help highlight their differences. Likewise, Fig. 6 depicts the 1-D circular transection profiles corresponding to Fig. 5. On the basis of these charts, it is apparently noted that the profiles using the constrained inverse solution [Figs. 6, 6, 6, 6, 6, 6] exhibit better reconstruction (three reconstructed peaks with higher contrast) than those [Figs. 6, 6, 6] using the primary inverse solution. It is worth noting that we can see the primary inverse solution remain a good reconstructed profile, and not be distorted in the case of three big inclusions, with even noise added. It can be explained that the same level of noise in boundary values has different effects on the phantoms with respectively two small or three big inclusions. Noise in the measured intensity affects more seriously on the reconstruction task for a phantom with small inclusions than that for a phantom with big inclusions.

Fig. 5

Same as in Fig. 3 except for the phantom with three big inclusions.

Fig. 6

Same as in Fig. 4 except for it being obtained from Fig. 5.

As what has been described above, ten cases at least were actually tested in our simulation. All the reconstructed images using the proposed schemes have a good convergence when maintaining a good quality as using the primary method. With the above discussion, it is concluded that the inverse solution regularized with the Lorentzian distributed function can be more robust to the noise in measured boundary intensities.

As described previously, the reconstructed images using the constrained inverse solution have superior qualities to those using the primary inverse solution. Furthermore, the proposed method can reconstruct optical-property images better than the primary method in the cases of noisy data, as shown in Figs. 4 and 6 (1D circular transection profiles). This section is to investigate the convergence rate as well as the performance measure associated with the above evaluated results. A computation scheme having a rapid convergence rate means that it can accelerate an implementation and reduce evaluation time. This feature is rather important for an NIR DOT imaging system to seek an inverse solution during the image reconstruction.

For the purpose of objectively evaluating the performance of image reconstruction on the convergence rate, a relative mean-square error (MSE) on the calculated intensity is adopted for the performance measure, given by the expression as

Figures 7 and 8 show the convergence rates for both examples discussed above, phantoms with two small inclusions and three big inclusions, where the cases with noise free, 30-, and $25\text{-}\mathrm{dB}$ SNRs are depicted. In each figure, the dotted line characterizes the result obtained from the primary inverse solution and the solid lines show the constrained inverse solutions with types 1 and 2 soft priors. To zoom in the MSE versus iterations, only MSEs ranging between 0 and 0.01% are charted. As will be seen, the convergence rate reaches a considerably small value of $<{10}^{-3}\mathrm{\%}$ $\left({10}^{-5}\right)$ .

Fig. 7

Convergence rates of the phantom with two small inclusions at various conditions: (a) noise free, Gaussian white noise with (b) $30\text{-}\mathrm{dB}$ SNR, and (c) $25\text{-}\mathrm{dB}$ SNR.

Fig. 8

Same as in Fig. 7 except for the phantom with three big inclusions.

For a phantom with two small inclusions, Fig. 7 explicitly demonstrates rapid convergence rates as using the constrained inverse solution to the cases of noise and noise free, whereas a slow convergence rate is obtained by using the primary inverse solution only for the noise-free case. In this example, the constrained inverse solutions with types 1 and 2 soft priors obtain similar convergence rates. For the noisy-data case, however, when using the primary inverse solution the convergence performance is less competitive, especially for a low value of SNR. This is predictable but unwanted as it may lead to a false result during minimizing the errors of calculated data against measurements.

For the phantom with three big inclusions, Fig. 8 illustrates rapid convergence rates for the constrained inverse solution to all the cases of noise and noise free, whereas a slow convergence rate is obtained in the primary inverse solution. In this example, the constrained inverse solutions with soft priors reach similar convergences. Additionally, the primary inverse solution can also achieve a slow convergence for the noisy-data case.

As shown in Figs. 7 and 8, it is found that only five iterations are required for a good convergence as using the constrained inverse solution, whereas the primary inverse solution takes 25 iterations or so to converge. Comparing Fig. 7 with Fig. 8, it is noted that noisy data affect the optical-property reconstruction of a phantom with small inclusions more severely than that of a phantom with big ones. For further observation in each panel within Figs. 7 and 8, constrained inverse solutions proposed in the study have similarly rapid convergence rates, implying that the choice of $\gamma$ is flexible. This shows the robustness of the proposed scheme in noisy-data cases and various inclusion conditions.

To further study the performance of prior constraints, two more noisy-data cases with 20- and $15\text{-}\mathrm{dB}$ SNRs are conducted to examine the robustness of the scheme. In Eq. 14, let

Fig. 9

Reconstruction data through various priors with intensity signals corrupted by Gaussian white noise $(\mathrm{SNR}=20\mathrm{dB})$ . Left column: constrained inverse solution with soft prior 1; middle column: constrained inverse solution with soft prior 2; right column: constrained inverse solution with hard prior.

Fig. 10

Same as in Fig. 9 except for $15\text{-}\mathrm{dB}$ SNR.

Figure 9 illustrates the comparisons between constrained solutions using soft priors (types 1 and 2) and a hard prior, where the left, middle, and right columns are the constrained inverse solutions with soft prior 1, soft prior 2, and hard prior, respectively. Figures 9, 9, 9, 9, 9, 9 show the 2-D reconstructions of phantoms with two and three inclusions, where slight discrepancy can be observed. Figures 9, 9, 9, 9, 9, 9 depict their corresponding 1-D circular transection profiles to reveal noticeable differences. Basically, there is a better separation resolution but a lower intensity owing to a highly suppressed signal by a hard prior rather than a soft prior. Additionally, Figs. 9, 9, 9 exhibit good convergences obtained by using both soft and hard priors.

Figure 10 evidently supports our thought that a hard prior to obtain a constrained inverse solution can be immune to more noisy $\left(15\mathrm{dB}\right)$ data than a soft prior is. In Fig. 10, first two columns obtained through soft-prior reconstructions, show a few distortions and fluctuations although a convergence remains reached; whereas the right column obtained by using a hard prior is more robust to noise, but traded off with contrast.

As has been shown, a hard prior can be immune to noise better than soft priors. But, according to our experience it is rather case sensitive to choose an adequate hard prior. To choose an appropriate hard prior is worthy of further study.

In addition to the investigation of convergence rates, global reconstruction qualities by the MSEs versus different noise levels are demonstrated in Fig. 11 . Figures 11 and 11 show the evaluation corresponding to Figs. 3 and 5, respectively, where the image reconstruction using priors with smaller MSEs means better results than that without a prior; whereas Figs. 11 and 11 show the evaluation corresponding to Figs. 9 and 10, where the reconstruction using a hard prior even has the smallest MSEs at different noise levels.

Fig. 11

Evaluation of the global reconstruction quality by the MSEs versus different noise levels: (a) Reconstruction for two inclusions through various priors corresponding to Fig. 3, (b) reconstruction for three inclusions through various priors corresponding to Fig. 5, (c) reconstruction for two inclusions through various priors corresponding to Figs. 9 and 10, and (d) reconstruction for three inclusions through various priors corresponding to Figs. 9 and 10.

4.

Conclusions

In this study, we have developed and realized the schemes for expediting NIR DOT image reconstruction through the inverse solution regularized with the constraint of a Lorentzian distributed function. Substantial improvements in reconstruction have been achieved without incurring additional hardware cost. It is valuable for the evaluation of NIR DOT being used as a medical imaging modality since promising results have revealed this probability.

With the introduction of constraints having a form of the Lorentzian distributed function, rapid convergence can be achieved owing to the fact that decreasing $\Delta \chi$ results in the increase of ${\lambda }^{\prime }$ as the iteration process proceeds and vice versa. It behaves like a criterion in the sense of a rapid convergence that the optimal iteration number is founded when seeking an inverse solution regularized with the Lorentzian distributed function. In the study, five or six iterations are recommended for the use of finding a constrained inverse solution [Eq. 14]. Besides, both the soft priors and a hard prior are introduced. The constrained inverse solution through soft priors can approach an expected value in a self-convergent manner based on expected decreasing updates $\left(\Delta \chi \right)$ . The solution through a hard prior, especially feasible for the case of severe noise, reaches the value in a compulsorily convergent manner.

We believe that the proposed method can be generalized as constraints in the regularization algorithm to lead to rapid convergence. There are, however, tremendous probabilities for further development in NIR DOT imaging systems, many of which have become current research topics such as the improvement in the iteration domain, the spatial, or spatial frequency domain with a refined mesh, structural, or spectral priors. Our further studies will explore these by means of either hardware design or algorithm exploitation.