Improvement of image quality of time-domain diffuse optical tomography with lp sparsity regularization

doi:10.1364/BOE.2.003334

Biomed Opt Express. 2011 December 1; 2(12): 3334–3348.

Published online 2011 November 21. doi: 10.1364/BOE.2.003334

PMCID: PMC3233252

Improvement of image quality of time-domain diffuse optical tomography with l_p sparsity regularization

Shinpei Okawa,¹^* Yoko Hoshi,² and Yukio Yamada¹

¹Department of Mechanical Engineering and Intelligent Systems, Graduate School of Informatics and Engineering, The University of Electro-Communications, 1-5-1 Chofugaoka, Chofu, Tokyo 182-8585 Japan

²Integrated Neuroscience Research Project, Tokyo Metropolitan Institute of Medical Science, 2-1-6 Kami-kitazawa, Setagaya, Tokyo 156-8506, Japan

^*Email: okawa/at/mce.uec.ac.jp

Received September 23, 2011; Revised November 14, 2011; Accepted November 14, 2011.

This is an open-access article distributed under the terms of the Creative Commons Attribution-Noncommercial-No Derivative Works 3.0 Unported License, which permits download and redistribution, provided that the original work is properly cited. This license restricts the article from being modified or used commercially.

This article has been cited by other articles in PMC.

Abstract

An l_p (0 < p ≤ 1) sparsity regularization is applied to time-domain diffuse optical tomography with a gradient-based nonlinear optimization scheme to improve the spatial resolution and the robustness to noise. The expression of the l_p sparsity regularization is reformulated as a differentiable function of a parameter to avoid the difficulty in calculating its gradient in the optimization process. The regularization parameter is selected by the L-curve method. Numerical experiments show that the l_p sparsity regularization improves the spatial resolution and recovers the difference in the absorption coefficients between two targets, although a target with a small absorption coefficient may disappear due to the strong effect of the l_p sparsity regularization when the value of p is too small. The l_p sparsity regularization with small p values strongly localizes the target, and the reconstructed region of the target becomes smaller as the value of p decreases. A phantom experiment validates the numerical simulations.

OCIS codes: (170.3880) Medical and biological imaging, (100.3190) Inverse problems

1. Introduction

Diffuse optical tomography (DOT) reconstructs the distributions of the optical properties, such as the scattering and absorption coefficients in the biological media. Near infrared (NIR) light incident on the surface of the biological medium propagates diffusively inside the biological medium, and the fractions of the NIR light are reemitted from and detected at the surface of the medium. By solving the inverse problem based on the photon diffusion equation with the measurement data of the detected light, DOT images are obtained noninvasively [1, 2]

The reconstructed optical properties contain the information about not only the structure but also the constituents in the measured object. By using the spectroscopic technique, the concentrations of hemoglobin, lipid etc. can be calculated from the optical properties. The abnormal absorption coefficients which can be caused by the tumor with angiogenesis are detected by DOT [3–7]. Therefore DOT attracts attentions as a new imaging modality for breast cancer diagnoses. DOT also can monitor the brain activities. The large change in the regional blood flow due to the higher brain activities is detectable by DOT [8, 9].

The advantage of the near infrared optical imaging modalities is that the measuring instruments are smaller in size and easier to use than the other modalities. Therefore, it is desirable that the number of the detectors is small. However, when the number of the detectors is limited, the inverse problem of DOT is highly ill-posed. The reconstructed optical properties usually have a low spatial resolution, and measurement noises significantly affect the reconstructed images. One of the solutions of these problems is to use a regularization method in the image reconstruction. So, recently the regularization techniques are actively studied for improvement of the DOT image quality.

The important approach is to incorporate the structural prior information. Pogue et al [10] uses a spatially variant Tikhonov regularization to reduce high frequency noises in the reconstructed images. Boverman et al employed prior segmentation of breast into glandular and adipose tissues [11]. Yalavarthy et al uses MRI-derived breast geometry to regularize the inverse solution [12]. Douiri et al applied anisotropic diffusion regularization with a priori edge information of the object to preserve the edge of the inner structure [13]. An algorithm assuming that the targets have blocky structures for regularization is introduced by Hiltunen et al [14]. Mutual information and joint entropy are used by Panagiotou et al to reflect the structure obtained from an alternative high resolution modality in the DOT images [15].

Another possible approach to improve the spatial resolution of the DOT image is to use a sparsity constraint. Cao et al implements L₁ norm minimization by use of an expectation maximization algorithm for a linearized DOT inverse problem, and show that the reconstructed region with abnormal optical properties are localized more than the other methods they use [16].

Sparsity regularization with L₁ norm and other methods to obtain sparse solution are used to several applications of inverse problems. Image restoration with sparsity constrained regularization is proposed by Shankar et al [17]. The L₁ sparsity constraint is applied to fluorescence/bioluminescence diffuse optical tomography (F/BDOT) [18, 19]. Okawa et al used recursive spatial filtering for F/BDOT [20]. And other applications are found in the functional brain imaging such as electroencephalography [21, 22].

When the changes in the optical properties, which can be caused by the regional blood flow changes or the breast cancer in early stage, are expected to be localized, the sparsity regularization will provide a good reconstruction in DOT inverse problem.

In this paper, we apply the l_p (0 < p ≤ 1) sparsity regularization by partial use of the focal underdetermined system solver (FOCUSS) algorithm, which was introduced by He et al. [23] to solve linear underdetermined inverse problems, and investigate the effects of the regularization on the non-linear time-domain DOT image reconstruction.

The changes in the absorption coefficients are described by a parameter to solve the difficulty in l_p minimization. And the DOT images are reconstructed by minimizing the residual error between the measurement and predicted data sets and the l_p norm of the changes in the absorption coefficients simultaneously. Numerical simulations and a phantom experiment demonstrate that the regularization improves the localization of the changes in the distribution of the absorption coefficient.

2. Methods

2.1. Forward problem in DOT

Light propagation in biological media is a phenomenon of the transport of radiative energy. Therefore, the radiative transport equation (RTE) rigorously describes the phenomenon. In DOT, NIR light illuminates the surface of the medium. The biological tissues strongly scatter and weakly absorb NIR light while propagating. The fraction of the light is reemitted from and detected at the surface of the medium. Because RTE is an integro-differential equation, solving and using RTE directly for the inverse problem is not an easy task. So the photon diffusion equation (PDE) is usually used [1]. PDE is a partial differential equation approximating RTE, and this approximation holds when the media are thicker than several centimeters and time after the pulse light incidence is longer than several hundred picoseconds.

Under the diffusion approximation, and by use of the spherical harmonic expansion of the quantities in RTE, following PDE for the fluence rate of light at the position r and time t, Φ(r,t), is obtained,

(1)

where D = 1/(3μ′_s) represents the diffusion coefficient with the reduced scattering coefficient μ′_s, μ_a the absorption coefficient, c the speed of light, q₀ the light source. The boundary condition is given as −n · D [nabla]

Φ = 1/(2A)Φ where n is the vector normal to the surface of the medium, and A is the parameter depending on the internal reflection ratio. The forward solution is obtained by solving Eq.(1) by use of the finite element method (FEM).

The fluxes of the fluence rate, Γ(r,t) = −n · D [nabla]

Φ, are the quantities measured by the detectors located at various positions. The mean time of flight (MTF) is calculated from Γ as equation m2

. In this study, a set of the measured MTFs, M, is used as the input data for reconstruction of the distribution of μ_a. An efficient method of solving the forward problem introduced by Schweiger et al [24] is also employed in this study to calculate MTFs using the moments of the time-resolved profiles.

2.2. Inverse problem with l_p sparsity regularization

Reconstruction of the distribution of μ_a is carried out by minimizing the residual error between the measured MTF data and the MTF data calculated by solving the forward problem. Reconstruction to estimate the μ_a distribution is conventionally achieved by solving the following optimization,

(2)

where M and are the sets of the measured and the calculated MTFs, respectively. λ is a regularization parameter, and f is a regularization function depending on the regularization method such as Tikhonov regularization method, the total variation method, etc [25]. The nonlinear optimization method such as Gauss-Newton method etc. is used.

To calculate the gradient of the residual error in the first term in Eq. (2), Jacobian matrix which represents the sensitivity of the optical properties to the measured data, _j/μ_{a_i}, is calculated by using FEM, where _j is the j-th calculated datum and μ_{a_i} is the absorption coefficient at the i-th FEM node.

To reconstruct the changes in μ_a localized in the small regions, we apply an l_p (0 < p ≤ 1) sparsity regularization which minimizes equation m4

where Δμ_{a_i} is a change in μ_a from the baseline at the i-th FEM node, and I is the number of the FEM nodes. However, the expression of the l_p (0 < p < 1) sparsity regularization has a difficulty in calculating the gradient in the optimization by use of the gradient based method, because |Δμ_{a_i}|^p⁻¹ tends to infinity when Δμ_{a_i} is minimized. To avoid this difficulty, Δμ_{a_i} is reformulated with a parameter z_i as follows [23],

(3)

Thus μ_{a_i} is described as

(4)

where _a is the constant baseline. Then the DOT reconstruction with the l_p sparsity regularization is represented as follows,

(5)

where z represents a vector consisting of z_i. By solving this optimization problem, a solution which selects the changes in μ_a in small localized region is obtained. As p becomes smaller, the localization is expected to be improved.

We use the nonlinear conjugate gradient method [26] for the optimization which requires the gradient of the cost function. Therefore, Jacobian matrix, J_ji(z_i) = _j/z_i, is calculated with the perturbation of z_i which is equivalent to 0.01 percent of the baseline _a in this study. The regularization parameter λ is selected with the L-curve method [27].

3. Numerical experiments

3.1. Conditions

Numerical experiments are conducted to investigate the effect of the l_p sparsity regularization with p = 1, 1/2, 1/4 on the image reconstruction of time-domain diffuse optical tomography and the results are compared with Tikhonov regularization which is identical to the l_p sparsity regularization with p = 2.

A μ_a distribution in a 2D circular medium with a radius of 40 mm is reconstructed with the regularizations. The medium has strongly absorbing targets. The sizes and the number of the targets are varied in the simulation. μ_a of the targets are changed, depending on the purpose of the experiments. Except the targets, the medium has μ_a = 0.007 mm⁻¹ and μ′_s = 0.8 mm⁻¹ homogeneously. These optical properties are given as those of breast in the NIR wavelength range by referring a literature [6]. The parameter A in the boundary condition for PDE is equal to unity because no reflection is assumed at the boundary.

16 positions working for both light sources and detectors are located at the boundary of the medium with an equal spacing. An ultra-short pulse light illuminates one of the 16 positions and the reemitted light after propagating inside the medium is detected at the rest of the 16 positions. Each position works as the light source one by one. Thus we obtain 16×15 MTF data.

By use of the same forward solver for generating simulated data and for solving inverse problem, the reconstruction can be fairly successful [28]. To avoid this ‘inverse crime’ problem, the measured MTF data sets are generated by solving PDE using FEM with 12800 elements and 6561 nodes, and the inverse problem is solved using FEM with 3200 elements and 1681 nodes. Gaussian noises with the standard deviation of one percent of MTF are added to the measured data, and the noisy data are used for reconstruction.

3.2. Results and discussions

3.2.1. Effect on the spatial resolution and localization of the reconstructed image

The spatial resolution of the reconstructed images is investigated. Two circular targets are placed in the medium. The centers of the targets are at (x,y) = (20 mm,10 mm) and (20 mm,−10 mm) with the origin of the coordinate located at the center of the medium. The radius and μ_a of both of the targets are 5 mm and 0.014 mm⁻¹ (Δμ_a = 0.007 mm⁻¹), respectively. μ′_s of the targets is identical to that of the background.

Figure 1 shows the L-curves with p = 2, 1, 1/2 and 1/4. Reconstructions are carried out with λ = 10⁻¹⁰, 10⁻⁹, 10⁻⁸,···,10⁻² for each p, and the L-curves are plotted with R = log₁₀ ||M – ||² for the abscissa and F = log₁₀f(z) for the ordinate. The corner of the L-curve, in which both terms in Eq. (5) are minimized with a good balance, is visually judged from the plots. The optimum λ determined at the corner decreases with the decrease in p.

Fig. 1

L-curves for the reconstructions with (a) p = 2, (b) p = 1, (c) p = 1/2, (d) p = 1/4 in the simulations in section 3.2.1. R = log₁₀ ||M–||² for the abscissa and F = log₁₀f(z) for the ordinate.

The reconstructed μ_a distributions are shown in Fig. 2 with the same color scale for all images using (a) Tikhonov regularization with p = 2 and (b) to (d) the l_p sparsity regularization with p = 1, 1/2 and 1/4, respectively. The maximum values of the reconstructed Δμ_a are 0.0015 mm⁻¹, 0.0020 mm⁻¹, 0.0036 mm⁻¹ and 0.0050 mm⁻¹ for p = 2,1,1/2 and 1/4, respectively. Figure 3(a) shows the profiles of Δμ_a at the FEM nodes along the lines in the y-direction passing through the two peaks of reconstructed μ_a. The two targets are observable and almost equally reconstructed with their centers at the correct positions. Especially with p = 1/2 and 1/4, the two targets are clearly separated.

Fig. 2

The μ_a distributions reconstructed using (a) Tikhonov regularization with p = 2 and using the l_p sparsity regularization with (b) p = 1, (c) p = 1/2 and (d) p = 1/4 in the simulations in section 3.2.1. The black open circles indicate the true (more ...)

Fig. 3

The quantitative evaluations of the reconstructed images obtained by the simulations in section 3.2.1: (a) The profiles of Δμ_a along the lines in the y-direction passing through the two peaks of μ_a, reconstructed with p = 2 (green (more ...)

Figure 3(b) plots ζ = Δμ_{a_min}/Δμ_{a_max} as a function of p, where Δμ_{a_min} is the minimum of Δμ_a reconstructed between the peaks, and Δμ_{a_max} is the average of the two peaks of Δμ_a. ζ is an index for evaluation of spatial resolution. ζ = 0 indicates that the two reconstructed targets are perfectly separated while ζ = 1 indicates that the two targets are not separated but combined into one large target. Figure 3 (b) shows that the degree of the separation is improved as p decreases.

The average of the peak values of Δμ_a, Δμ_{a_max}, is plotted as a function of p in Fig. 3(c). The absolute value of Δμ_a of the reconstructed target with small p becomes closer to the true Δμ_a of 0.007 mm⁻¹ in the targets than that with large p.

The area of the peak, S, is defined as the sum of the area of the FEM elements having Δμ_a ≥ Δμ_{a_max}/2, where Δμ_{a_max} is the maximum of Δμ_a, and is plotted as a function of p in Fig. 3 (d). S is an index for evaluation of localization by comparing with the true value of S. S decreases with the decrease in p, and localization of the reconstructed targets is improved.

The μ_a distribution reconstructed using Tikhonov regularization (Fig. 2(a)) has small undulations around the targets, which are caused by the random noise, and the difference in the FEM meshing between that for generating the measurement data and that for solving the inverse problem.

On the other hand, undulations are hardly observed in the μ_a distributions reconstructed with the l_p sparsity regularization. Therefore, the l_p sparsity regularization reduces the influences of noise and of the difference in the FEM meshing. From these results, it can be said that the l_p sparsity regularization achieves a high spatial resolution. When the true changes in μ_a are localized in small regions, the l_p sparsity regularization provides preferable reconstruction with the robustness to noise.

3.2.2. Effect on the sensitivity to small changes in the absorption coefficient

In this subsection we investigate the sensitivity of reconstruction to small changes in the absorption coefficient by reconstructing μ_a of a medium having two targets with different μ_a. The positions and the sizes of the targets are the same as those in the previous subsection. One of the targets with its center at (x,y) = (20 mm, 10 mm) has μ_a = 0.0105 mm⁻¹, and the other target with its center at (x,y) = (20mm,−10mm) has μ_a = 0.0140 mm⁻¹. The difference in μ_a from the background μ_a = 0.0070 mm⁻¹ of the former target (Δμ_a₁ = 0.0035 mm⁻¹) is half of that of the latter one (Δμ_a₂ = 0.0070 mm⁻¹). Therefore, the ratio of the smaller changes in μ_a from the background to larger one, defined as γ = Δμ_a₁/Δμ_a₂, is 0.5. γ is an index for evaluating sensitivity to small changes in μ_a by comparing with the true value of γ.

Reconstructions are carried out with the same manner as in the previous section. λ is selected by the L-curve method.

Figure 4 shows the reconstructed images. The μ_a image reconstructed using Tikhonov regularization shown in Fig. 4(a) reveals weak two peaks with Δμ_{a_max} = 0.0009 mm⁻¹ and 0.0015 mm⁻¹ at (x,y) = (23.0 mm, +12.1 mm) and (25.2 mm, −12.1 mm), respectively. The positions of the targets are well reconstructed. However, γ calculated from the reconstructed peaks of Δμ_a is 0.60, and the reconstructed μ_a are underestimated in this case. As well as in the case of the previous subsection, the spatial resolutions of the reconstructed images are low, and the small undulations due to noise are seen in Fig. 4(a).

Fig. 4

Reconstructed μ_a distributions using (a) Tikhonov regularization with p = 2 and the l_p sparsity regularization with (b) p = 1, (c) p = 1/2 and (d) p = 1/4 in the simulations in section 3.2.2. The black open circles indicate the true positions (more ...)

When the l_p sparsity regularization with p = 1 is used, undulation due to noise is suppressed, and localization of the reconstructed targets is improved. Two peaks of the reconstructed Δμ_a are 0.0013 mm⁻¹ and 0.0027 mm⁻¹, respectively. γ is 0.48 which is close to the true value of 0.50. This means that the l₁ sparsity regularization recovers the difference in μ_a between the two targets better than Tikhonov regularization.

Nevertheless the l_p sparsity regularization is not always successful. When p = 1/2, one of the target is lost in the reconstructed image in Fig. 4(c). The changes in μ_a are so excessively localized that the weakly absorbing target disappears. Two peaks are found in the μ_a distribution reconstructed with p = 1/4 in Fig. 4(d), but γ is 0.23 which is much smaller than the true value of 0.50. Figures 5(a) and (b) show the reconstructed peaks of Δμ_a of the two targets, Δμ_{a1_max} and Δμ_{a2_max}, and γ = Δμ_{a1_max}/Δμ_{a2_max} calculated from the reconstructed peaks of Δμ_a as a function of p, respectively. The l_p sparsity regularization is effective to enhance the localization, to reduce the artifacts and to recover the difference in μ_a between the large and small changes, although small changes in μ_a may be eliminated by the l_p sparsity regularization when p is too small because it strongly localizes the changes in μ_a.

Fig. 5

The quantitative evaluations of the reconstructed images obtained by the simulations in section 3.2.2: (a) The reconstructed peaks of Δμ_a of the two targets at (x,y) = (20 mm, 10 mm) (times symbols, true Δμ_{a1_max} = 0.0035 (more ...)

3.2.3. Effect on the reconstruction of broad target

The l_p sparsity regularization is found to be effective for reconstructing localized targets as shown in the previous subsections. However, the targets are not always localized well in practical applications. We demonstrate reconstructions of broad targets in this subsection.

The circular target has a radius of 10 mm, a center at (x,y) = (20 mm, 0 mm) and μ_a of 0.014 mm⁻¹ in the background medium having μ_a = 0.0070 mm⁻¹. Reconstructions are carried out with the manner mentioned above.

Figure 6 shows the reconstructed μ_a distributions using Tikhonov and the l_p sparsity regularizations with p = 1, 1/2, and 1/4. The maximum μ_a reconstructed using Tikhonov regularization shown in Fig. 6(a) is 0.0122 mm⁻¹. Undulations in the μ_a distribution appear due to the noise in the input data.

Fig. 6

Reconstructed μ_a distributions using (a) Tikhonov regularization with p = 2 and the l_p sparsity regularization with (b) p = 1, (c) p = 1/2 and (d) p = 1/4 in the simulations in section 3.2.3. The black open circle indicates the true position of (more ...)

On the other hand, the l_p sparsity regularizations remove the undulations. The shape and size of the target are well reconstructed when p = 1 as shown in Fig. 6(b). In Fig. 6(c), the target reconstructed with p = 1/2 has a smaller area than that with p = 1, and the maximum μ_a value of 0.0130 mm⁻¹ with p = 1/2 is larger than that of 0.0117 mm⁻¹ with p = 1. The maximum μ_a value further increases to 0.0184 mm⁻¹ in Fig. 6(d) as p decreases to 1/4 although the reconstructed target is localized excessively and smaller in size than the true target. Figures 7(a) and (b) plot the reconstructed peak value of Δμ_a, Δμ_{a_max}, and the area of the peak of Δμ_a, S, as a function of p, respectively. S is defined as the sum of the area of the FEM elements having Δμ_a ≥ Δμ_{a_max}/2, and S is the index for evaluating localization of the reconstructed target by comparing with the true value of S as mentioned before.

Fig. 7

The quantitative evaluations of the reconstructed images obtained by the simulations in section 3.2.3: (a) Reconstructed peaks of Δμ_a with the true Δμ_a = 0.0070 mm⁻¹, and (b) the area at half maximum of Δ (more ...)

From the results of Figs. 6 and and7,7, it can be said that the l_p sparsity regularization is effective to improve the quality of the reconstructed images of broad targets by reducing the influence of noise. However, the quality of the reconstructed images highly depends on the parameter p, and the size of the reconstructed target becomes smaller than that of the true one as p decreases too much.

3.2.4. Determination of optimum p value

In the simulations above, we investigated the effects of the p value on the quality of the reconstructed images in the cases of various sizes and μ_a of the targets. According to the simulation results, the optimum p value depends on the cases. For the case in the section 3.2.1, the size and μ_a of the small targets with a radius of 5 mm is reconstructed clearly when p is small, and p = 1/4 is a good choice in this case as shown in Figs. 3(c) and (d).

However, other simulation results indicate that smaller p does not always lead to better reconstruction. For the case of multiple targets with different μ_a, too small p values underestimate μ_a of the target with smaller μ_a. From the results in the subsection 3.2.2, p = 1/4 may be a good choice from the Δμ_a point of view (Fig. 5(a)), but p = 1 may be a good choice from the γ point of view (Fig. 5(b)). When the target has a radius of 10 mm, the size of the target is well reconstructed with p = 1, although p = 1/2 reconstructs the peak of Δμ_a of the target better than p = 1 as plotted in Fig. 7. As a whole it can be said that p ≤ 1 is preferable to reduce the measurement noise and the error between the true and reconstructed μ_a values.

A criterion will be required for determination of the optimum p value. Prior information provided by other imaging modalities may be useful for that purpose. Prior information of the target size may help determine the optimum p value, for example.

4. Phantom experiment

4.1. Conditions

The effect of the l_p sparsity regularization is validated by a phantom experiment in this section. The time-resolved measurement system consisted of an ultra-short pulse laser operating at the wavelength of 759 nm, time-correlated single-photon counting units and 16 source/detector optical fiber bundles. The ultra-short pulse light with a duration of about 100 ps, a mean power of 0.25 mW and a repetition rate of 5 MHz illuminated the measured object. The details of the time-resolved measurement system is found in the literature [29].

One of the optical fiber bundles worked as a light source and the others detected light reemitted from the surface of the measured object. Therefore, we acquired 16×15 time-resolved data. The MTF dataset obtained from the measured time-resolved data at the wavelength of 759 nm is used as the input for reconstruction of the μ_a distribution.

In the phantom experiment, we reconstructed μ_a in a 2D-like tissue simulating phantom which was a cylinder made of polyacetal resin with a height of 240 mm and a radius of 40mm. The phantom had the background optical properties of μ_a = 0.0006 mm⁻¹ and μ′_s = 0.863 mm⁻¹. In this phantom, there existed a target of a cylindrical hole with a radius of 10 mm and the center at (x,y) = (20 mm, 0 mm). The cylindrical hole was filled with 1.0 percent Intralipid solution having μ_a = 0.0026 mm⁻¹ and μ′_s = 1.054 mm⁻¹.

The optical fiber bundles were circularly attached on the surface of the phantom with an equal spacing in a 2D plane perpendicular to the cylinder axis. The μ_a distribution in the 2D plane was reconstructed with Tikhonov or the l_p sparsity regularizations. We used the FEM meshing identical to that used for reconstructions in the previous simulations. We set _a = 0.0006 mm⁻¹ and selected λ based on the L-curve method. The refractive index of the polyacetal resin was given as 1.54 leading to the value of A = 4.26 in the boundary condition of PDE.

4.2. Results and discussions

Figure 8 shows the L-curves plotted with the evaluating points with λ = 10⁻¹⁰,10⁻⁹,10⁻⁸, (...)

,10⁻². λ values at the corners are determined as 10⁻² and 10⁻⁵ for p =2 and 1, respectively. However, it was difficult to obtain typical L-curves for p = 1/2 and 1/4. The evaluating points for p = 1/2 and for 1/4 in Fig. 8(c) and (d) subtly form the L-curves, and λ = 10⁻⁷ is at the corners in both cases. More points on the L-curve for p = 1/4 were evaluated to confirm the corner as shown in Fig. 8(d). In the simulation sections and this phantom experiment, the regularized solution with small p tends to change drastically within a narrow range of λ. So it is better to make the L-curve with fine evaluating points of λ to find an appropriate corner of the L-curve.

Fig. 8

L-curves for the reconstructions with (a) p = 2, (b) p = 1, (c) p = 1/2, (d) p = 1/4 in the phantom experiment. R = log₁₀ ||M – ||² for the abscissa and F = log₁₀f(z) for the ordinate

The reconstructed μ_a distributions are shown in Fig. 9. As expected from the results in the previous section, the smaller the value of p is, the more localized the target is. Tikhonov regularization reconstructed the target clearly as shown in Fig. 9(a), but there are undesirable changes in μ_a in the wide area out of the true target position which were due to unknown noise factors and took the values from 38 percent to 63 percent of the reconstructed μ_a of the target. The area with the undesirable changes in μ_a became smaller by the l_p sparsity regularization, although some of the undesirable large changes still remained near the position of (x,y) = (−20 mm, −10 mm) and (−10 mm, 25 mm) in Figs. 9 (b), (c) and (d). The l_p sparsity regularizations with p = 1/2 and 1/4 provide highly localized targets. The maximum value of μ_a became larger as p decreased.

Fig. 9

The μ_a distributions reconstructed with (a) Tikhonov regularization with p = 2 and the l_p sparsity regularization with (b) p = 1, (c) p = 1/2 and (d) p = 1/4 in the phantom experiment. The black open circle indicates the true position of the target. (more ...)

The results of the phantom experiment indicates that the l_p sparsity regularization localizes the change in the μ_a distribution and that the regularization is effective to reduce the influence of noise. Figure 10 shows Δμ_{a_max} and S in this phantom experiment. Small p enhances Δμ_{a_max}, and the size of the reconstructed target becomes smaller than that of the true one as p decreases. These results validate the simulations in the previous section.

Fig. 10

The quantitative evaluations of the reconstructed images obtained by the phantom experiment: (a) Reconstructed peaks of Δμ_a with the true Δμ_a = 0.0020 mm⁻¹, and (b) the area at half maximum of Δμ (more ...)

The maximum values of the reconstructed Δμ_a are 0.0047 mm⁻¹, 0.0068 mm⁻¹ and 0.0132 mm⁻¹ for p = 1, 1/2 and 1/4, respectively, and all of these values are larger than the true value of 0.002 mm⁻¹. These differences may be caused by the small changes in μ′_s of the target from that of the background because only μ_a is reconstructed under the assumption of homogeneous μ′_s. Simultaneous reconstruction of μ_a and μ′_s is an interesting topic for future study. Unknown noise factor must prevent a good reconstruction also for the cases of multiple targets. Although the phantom experiment was conducted only for single target, the effects of the proposed method, i.e. improvement of localization and alleviation of influence of noise, are validated and similar results will be obtained for the cases of multiple targets.

5. Conclusion

To improve the spatial resolution and the robustness to noise, the l_p (0 < p = 1) sparsity regularization is applied to the time-domain diffuse optical tomography with the gradient-based nonlinear optimization scheme. The expression of the l_p sparsity regularization is reformulated as a differentiable function of a parameter to avoid the difficulty in calculating its gradient in the optimization process. The regularization parameter is selected by the L-curve method.

Numerical experiments show that the l_p sparsity regularization improves the spatial resolution. Two localized targets with the absorption coefficients higher than that of the background are clearly separated in the images reconstructed with the l_p sparsity regularization, while Tikhonov regularization reconstructs the two targets as a broad single target. The reconstructed values of the absorption coefficients approaches the correct value by use of the l_p sparsity regularization.

The difference in the absorption coefficient between targets at different positions is also recovered well by the l_p sparsity regularization. However, it is demonstrated that a target with a small differences in the absorption coefficient from the background may disappear due to the excessive effect of the l_p sparsity regularization.

The l_p sparsity regularization is also effective for reconstructing broad targets. The influence of noise such as non-targeted small undulations in the reconstructed image is reduced. The size of the reconstructed target is reconstructed well when p = 1 and 1/2 in the simulation. The l_p sparsity regularization with a small p strongly localizes the target, and the size of the reconstructed target decreases with the decrease in p. A phantom experiment validates the numerical simulations.

The l_p sparsity regularization can be useful to reconstruct the localized changes in the absorption coefficient. The criterion to determine the optimum p value can be discussed in a future work.

References and links

1. Arridge S. R., “Optical tomography in medical imaging,” Inverse Prob. 15, R41–R93 (1999). doi: 10.1088/0266-5611/15/2/022. [Cross Ref]

2. Gibson A. P., Hebden J. C., Arridge S. R., “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005). doi: 10.1088/0031-9155/50/4/R01. [PubMed] [Cross Ref]

3. Grosenick D., Wabnitz H., Rinneberg H. H., Moesta T., Schlag P. M., “Development of a time-domain optical mammography and first in vivo applications,” Appl. Opt. 38(13), 2927–2943 (1999). doi: 10.1364/AO.38.002927. [PubMed] [Cross Ref]

4. Grosenick D., Moesta K. T., Möller M., Mucke J., Wabnitz H., Gebauer B., Stroszczynski C., Wassermann B., Schlag P. M., Rinneberg H., “Time-domain scanning optical mammography: I. Recording and assessment of mammograms of 154 patients,” Phys. Med. Biol. 50, 2429–2449 (2005). doi: 10.1088/0031-9155/50/11/001. [PubMed] [Cross Ref]

5. Grosenick D., Wabnitz H., Moesta K. T., Mucke J., Schlag P. M., Rinneberg H., “Time-domain scanning optical mammography: II. Optical properties and tissue parameters of 87 carcinomas,” Phys. Med. Biol. 50, 2451–2468 (2005). doi: 10.1088/0031-9155/50/11/002. [PubMed] [Cross Ref]

6. Hebden J. C., Veenstra H., Dehghani H., Hillman E. M. C., Schweiger M., Arridge S. R., Delpy D. T., “Three-dimensional time-resolved optical tomography of a conical breast phantom,” Appl. Opt. 40 (19), 3278–32887 (2001). doi: 10.1364/AO.40.003278. [PubMed] [Cross Ref]

7. Yates T., Hebdan C., Gibson A., Everdell N., Arridge S. R., Douek M., “Optical tomography of the breast using a multi-channel time-resolved imager,” Phys. Med. Biol. 50, 2503–2517 (2005). doi: 10.1088/0031-9155/50/11/005. [PubMed] [Cross Ref]

8. Gibson A. P., Austin T., Everdell N. L., Schweiger M., Arridge S.R., Meek J. H., Wyatt J. S., Delpy D. T., Hebden J. C., “Three-dimensional whole-head optical tomography for passive motor evoked responses in the neonate,” NueroImage 30, 521–528 (2006). doi: 10.1016/j.neuroimage.2005.08.059. [PubMed] [Cross Ref]

9. Hebden J. C., Gibson A., Yusof R. M., Everdell N., Hillman E. M. C., Delpy D. T., Arridge S. R., Austin T., Meek J. H., Wyatt J. S., “Three-dimensional optical tomography of the premature infant brain,” Phys. Med. Biol. 47, 4155–4166 (2002) doi: 10.1088/0031-9155/47/23/303. [PubMed] [Cross Ref]

10. Pogue B. W., McBride T. O., Prewitt J., Lösterberg U., Paulsen K. D., “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt. 38(13), 2950–2961, (1999). doi: 10.1364/AO.38.002950. [PubMed] [Cross Ref]

11. Boverman G., Miller E. L., Li A., Zhang Q., Chaves T., Brooks D. H., Boas D. A., “Quantitative spectroscopic optical tomography of the breast guided by imperfect a priori structural information,” Phys. Med. Biol. 50, 3941–3956 (2005). doi: 10.1088/0031-9155/50/17/002. [PubMed] [Cross Ref]

12. Yalavarthy P. K., Pogue B. W., Dehghani H., Carpenter C. M., Jiang S., Paulsen K. D., “Structural information within regularization matrices improves near infrared diffuse optical tomography,” Opt. Express 15(13), 8043–8058, (2007). doi: 10.1364/OE.15.008043. [PubMed] [Cross Ref]

13. Douiri A., Schweiger M., Riley J., Arridge S. R., “Anisotropic diffusion regularization methods for diffuse optical tomography using edge prior information,” Meas. Sci. Technol. 18, 87–95 (2007). doi: 10.1088/0957-0233/18/1/011. [Cross Ref]

14. Hiltunen P., Calvetti D., Somersalo E., “An adaptive smoothness regularization algorithm for optical tomography,” Opt, Express 16(24), 19957–19977, (2008). doi: 10.1364/OE.16.019957. [PubMed] [Cross Ref]

15. Panagiotou C., Somayajula S., Gibson A. P., Schweiger M., Leahy R. M., Arridge S. R., “Information theoretic regularization in diffuse optical tomography,” J. Opt. Soc. Am. A 26(5), 1277–1290 (2009). doi: 10.1364/JOSAA.26.001277. [PubMed] [Cross Ref]

16. Cao N., Nehorai A., Jacob M., “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Opt. Express, 15(21), 13695–13708 (2007). doi: 10.1364/OE.15.013695. [PubMed] [Cross Ref]

17. Shankar P. M., Neifeld M. A., “Sparsity constrained regularization for multiframe image restoration,” J. Opt. Soc. Am. A 25(5), 1199–1214 (2008). doi: 10.1364/JOSAA.25.001199. [PubMed] [Cross Ref]

18. Mohajerani P., Eftekhar A. A., Huang J., Adibi A., “Optimal sparse solution for fluorescent diffuse optical tomography: theory and phantom experimental results,” Appl. Opt. 46(10), 1679–1685 (2007). doi: 10.1364/AO.46.001679. [PubMed] [Cross Ref]

19. Lu Y., Zhang X., Douraghy A., Stout D., Tian J., Chan T. F., Chatziioannou A. F., “Source reconstruction for spectrally-resolved bioluminescence tomography with sparse a priori information,” Opt. Express 17(10), 8062–8088 (2009). doi: 10.1364/OE.17.008062. [PubMed] [Cross Ref]

20. Okawa S., Yamada Y., “Reconstruction of fluorescence/bioluminescence sources in biological medium with spatial filter,” Opt. Express 18(12), 13151–13172 (2010). doi: 10.1364/OE.18.013151. [PubMed] [Cross Ref]

21. Baillet S., Mosher J. C., Leahy R. M., “Electromagnetic brain mapping,” IEEE Signal Process. Mag. 18, 14–30 (2001). doi: 10.1109/79.962275. [Cross Ref]

22. Xu P., Tian Y., Chen H., Yao D., “L_p Norm iterative sparse solution for EEG source localization,” IEEE Trans. Biomed. Eng. 54 (3), 400–409 (2007). doi: 10.1109/TBME.2006.886640. [PubMed] [Cross Ref]

23. He Z., Cichocki A., Zdunek R., Xie S., “Improved FOCCUS Method With conjugate gradient iterations,” IEEE Tras. Signal Process. 57 (1), 399–404 (2009). doi: 10.1109/TSP.2008.2007605. [Cross Ref]

24. Schweiger M., Arridge S. R., Delpy D. T., “Application of the finite-element method for the forward and inverse model in optical tomography,” J. Math. Imaging Vis. 3, 263–283 (1993). doi: 10.1007/BF01248356. [Cross Ref]

25. Vogel C. R., Computational Methods for Inverse Problems, Frontiers in Applied Mathematics (SIAM, Philadelphia, 2002). doi: 10.1137/1.9780898717570. [Cross Ref]

26. Arridge S. R., “A gradient-based optimization scheme for optical tomography,” Opt. Express 12(6), 213–226 (1998). doi: 10.1364/OE.2.000213. [PubMed] [Cross Ref]

27. Hansen P. C., O’Leary D. P., “The use of the L-curve in the regularization of discrete ill-posed problems,” SIAM J. Sci. Compt. 14(6), 1487–1503 (1993). doi: 10.1137/0914086. [Cross Ref]

28. Holder S., Electrical Impedance Tomography: Methods, History and Applications (Institute of Physics Publishing, Bristol, 2005).

29. Eda H., Oda I., Ito Y., Wada Y., Oikawa Y., Tsunazawa Y., Tuchiya Y., Yamashita Y., Oda M., Sassaroli A., Yamada Y., Tamura N., “Multichannel time-resolved optical tomographic imaging system,” Rev. Sci. Instrum. 70(9), 3595–3602 (1999). doi: 10.1063/1.1149965. [Cross Ref]

Articles from Biomedical Optics Express are provided here courtesy of
Optical Society of America