Improving image quality of diffuse optical tomography with a projection-error-based adaptive regularization method

Diffuse optical tomography (DOT) reconstructs the images of internal optical parameter distribution using noninvasive boundary measurements. The image reconstruction procedure is known to be an ill-posed problem. In order to solve such a problem, a regularization technique is needed to constrain the solution space. In this study, a projection-error-based adaptive regularization (PAR) technique is proposed to improve the reconstructed image quality. Simulations are performed using a diffusion approximation model and the simulated results demonstrate that the PAR technique can improve reconstruction precision of object more effectively. The method is demonstrated to have low sensitivity to noise at various noise levels. Moreover, with the PAR method, the detectability of an object located both at the center and near the peripheral regions has been increased largely. ©2008 Optical Society of America OCIS codes: (170.6960) Tomography; (170.3010) Image reconstruction techniques; (170.3880) Medical and biological imaging


Introduction
Interest in the use of near-infrared diffuse optical tomography (DOT) imaging for early detection of breast cancer has increased in recent years [1][2][3].DOT can be used to quantify hemoglobin concentration and blood oxygen saturation in tissue by imaging internal optical absorption and scattering [2,4], which allows non-invasive detection and diagnosis.
Near-infrared DOT imaging provides a number of advantages, including portability, real-time imaging and low instrumental cost [5,6] but is generally known to have low image resolution, which limits its further clinical application.Various methods, such as optimizing optode arrangements in the test fields [7,8], using approximated analytical approaches for solving the inverse problem [9,10], improving iteration strategies [11][12][13] etc., have been developed to improve DOT imaging.Among the numerous methods for enhancing image quality in DOT, the regularization approach has been shown to be effective [14,15] because it can decrease the ill-posed characteristics of the inverse matrix, consequently improving image quality.
Pogue et al. [14] showed that using spatially variant regularization parameters can diminish typical artifacts at source locations and improve the reconstruction precision of the entire image.In addition, the single-step Tikhonov regularization (STR) method [16] and the modified Levenberg-Marquardt regularization (LMR) approach [17] are frequently utilized to improve the image quality of near-infrared (NIR) tomography.These two methods have also been shown to be effective for enhancing the quantification accuracy of small objects [18].Recently, Davis et al. compared both methods under several conditions and concluded that the STR method is better than the LMR method since the STR images are smoother with fewer artifacts compared with the images using LMR [19].However, a common problem for these methods is that the choice of appropriate regularization is time-consuming and usually requires trial-and-error.
In this study, a projection-error-based adaptive regularization technique (PAR) is proposed for improving the image quality of NIR reconstructions.To illustrate the performance of this method, we reconstructed images using the proposed PAR technique under different imaging conditions.This new technique is also compared with the widely-used STR and LMR methods under same conditions.
The organization of this paper is as follows: in section 2, we review the theory of forward and inverse problems in DOT as well as introduce the PAR approach.In section 3, several simulations are demonstrated and their results are displayed in section 4. Finally, a discussion and conclusions are provided with respect to the PAR method in last section.

Theory
Photon propagation in tissue can be described by the diffusion equation, which is given for the frequency domain [20] by where Φ(r,ω) is the complex radiance at position r and modulation frequency ω (ω=2πf, in this work f =100 MHz).q 0 (r,ω) is an isotropic source; c is the wave speed in the medium; and μ a is an absorption coefficient.The optical diffusion coefficient κ(r) can be written as follows: ( ) where μ s ′ is a reduced scattering coefficient.The Robin-type (type III) boundary condition is used which best describes the light interaction from a scattering medium to the external air boundary.
Forward calculations with the diffusion equation are accomplished with a finite element method (FEM) [20], which is considered a flexible and accurate approach to modeling heterogeneous domains with arbitrary boundaries.The details of the FEM formulation of the forward model are given in [20][21][22][23][24].A 2-D imaging region which was 86 mm in diameter was used in the simulation studies.The region was separated into 1,785 nodes and 3,418 triangular elements.A symmetric circular configuration of 16 sources and 16 detectors was arranged around the edge of the imaging fields.This configuration produced 480 amplitude and phase measurements.The absorption and reduced scattering coefficients are μ a =0.01 mm -1 and μ s ′ =1.0 mm -1 , respectively.
For the inverse problem, an iterative solution is commonly used.We utilized the Newton-Raphson iterative method to find the optical parameter solutions by solving the minimum of the objective function: In Eq. (3), Φ c and Φ m are the calculated and measured radiance at the detectors, respectively.λ is the regularization parameter.μ and μ 0 are the current and initial estimates of the optical properties at each node, respectively.Based on a homogeneous estimate of the initial values for absorption and scattering properties, the updated optical distribution can be obtained with Eq. ( 4) by applying the Tikhonov regularization: ( ) where Δμ is the optical parameter update vector.H max is the maximum main diagonal element value of the matrix JJ T .Τhe regularization factor λ was selected by the adaptive technique given in the following section.J is the Jacobian matrix for the inverse problem and it maps the changes in log amplitude and phase to both absorption and diffusion changes at each node [25].Random Gaussian noise of 1% amplitude and 1 degree phase was added to the simulated data during image reconstruction (the 1% noise has been established to be the level of shot noise in experimental systems [16,26]).The iterative procedure is terminated when the difference of the objective function values between two successive iterations is less than 2%, which in the limit can ensure stable solution to be obtained [18].

PAR technique
Generally, the choice of an appropriate regularization can be evaluated or found empirically [14].However, theoretical analyses for selecting the regularization parameters exist [27].These analyses show that the regularization parameter is related to objective function.In our study, projection error is utilized as the objective function; hence the regularization parameter would be expected to be related to projection error.Furthermore, the value of regularization is required to vary with projection error.That is to say, a greater regularization value is needed for a larger projection error during iterations.If the projection error is very low, only a small regularization value is needed to regulate the ill-posed process.
Based on these considerations, we define the adaptive regularization parameter λ as follows: where ΔΦ is the projection error, which can be considered as the sum of the squared elements.The projection error also determines the accuracy of the match between the measured data Φ m and the calculated data Φ c .Moreover, projection error can reflect changes in object contrast, size, and position; thus different regularization values can be guaranteed and chosen for different imaging conditions.
In formula (5), the regularization parameter λ varies in the range from 1/3 to 1/2.It includes two primary considerations: (1) Due to the ill-posed characteristic of the inverse problem, the regularization should not be reduced with iterations to too small of a value.(2) The regularization value should be quantitatively similar to the diagonal of the matrix JJ T .In contrast to our previous study [28], expression (5) for the regularization parameter λ does not includes any empirical coefficients, so formula ( 5) is more convenient for users.

Evaluation criteria of reconstructed images
A contrast to noise ratio [29] (CNR) was used in this work to evaluate the reconstructed images.This ratio is defined as the difference between the averaged optical coefficient within the region of interest (ROI) and the region of background (ROB) divided by the weighted average of the standard deviations in the ROI and ROB.The CNR was obtained by the equation: where ω ROI is the fraction of the area of the region of interest with respect to the area of the whole image; ω ROB is defined as ω ROB =1-ω ROI ; μ ROI and μ ROB are the mean values of the object and the background regions in the reconstructed images; and σ ROI and σ ROB are the corresponding standard deviations.CNR indicates whether an object can be clearly seen in the reconstructed images.Greater CNR values are considered to imply a better image quality.A threshold of CNR=3 was chosen empirically as the lowest acceptable reconstruction for this paper.

Simulations
In our experiments, the 2-D imaging field was a circle with a radius of 43 mm and the center at (0, 0).A circular object located in the center of this imaging field was moved along the y-axis from 0 to 40 mm in 2 mm steps.For each object position, the object size was changed from 1 mm to 20 mm with 1 mm intervals.The absorption and reduced scattering coefficients for the object were μ a =0.02 mm -1 and μ s ′=2.0 mm -1 , respectively.For each position and size pair, an image was reconstructed (with 1% noise) and its corresponding CNR value was calculated.For comparison, the STR and LMR [16,18] methods were also used under the same imaging conditions.The specific regularization technique used in this study for these two methods can be found in [19].The initial regularization value λ for STR and LMR in each object imaging depended on the object's position y on the y-axis, and a linear relationship between λ and y was adopted.For both these methods, a value of 1 was used as the initial regularization for the centrally (i.e.y=0) located objects, and for the peripherally (i.e.y=40) located objects, 10 and 100 were used respectively.It was noted that in LMR, the regularization parameter was decreased by a factor of 10 0.5 in each successive iteration.
To demonstrate the low sensitivity of the PAR method to various noise levels, we examined two different object positions: one was with an object embedded in the center of the imaging field; the other was with an object positioned near the edge of the field.The center of the imaging field was considered to be less sensitive to inclusions, and hence would have a lower resolution and poorer contrast recovery.Accordingly, imaging near the edge of the test fields could be expected to have a better quality.In both cases the object sizes were set to be 12 mm in diameter, and the interaction coefficients were μ a =0.02 mm -1 and μ s ′ =2.0 mm -1 , which provided a 2:1 contrast level relative to the background.Moreover, we added random Gaussian noise of 1% in amplitude and 1 degree phase (this is referred to as a 1% noise level hereafter), 5% amplitude and 5 degree phase (i.e. 5% noise level), and 10% amplitude and 10 degree phase (i.e.10% noise level) to the simulated measurement data.Absorption and scattering images were reconstructed using STR, LMR and PAR at these three noise levels, respectively.The STR regularization values that were used were 1 for the central object and 10 for the peripheral object.The LMR regularization values that were used were 1 and 100 respectively.
In order to illustrate the improvement of spatial resolution and contrast resolution in object detection by the PAR approach, we located an object in three positions: at the center, at 23 mm and 3 mm away from the edge of the imaging field.The object diameter was increased from 1 mm to 20 mm.The contrasts varied from 1.1 to 4.0 and 4.0 were regarded as the maximum contrast between normal tissue and heterogeneous tissue in the breast.The intervals for the diameters and contrast levels were 1.0 mm and 0.1, respectively.For each object size and contrast, the image was reconstructed with 1% added noise and the corresponding CNR values were calculated.

Improvement in image quality with the PAR method: CNR analysis
The CNR graphs obtained using STR and LMR and our proposed PAR are shown in Figs.1(a), 1(b), and 1(c), respectively.The fractional areas with CNR values more than 3 in Figs.1(a), 1(b), and 1(c) are 296, 290 and 372, respectively.For objects with diameters more than 3 mm, the CNR values from the PAR method are larger than 3 at every position.In addition, for a specified position and size, the CNR value using the PAR method is greater than that of the STR and LMR approaches.These results indicate that a better image quality of the reconstructed images has been achieved using the PAR technique.
The objects less than 3 mm in diameter can also be recognized by using the PAR method as long as they are embedded near the edge of the imaging fields.However, none of the objects can be recognized using the STR and LMR techniques when the object sizes are less than 3 mm.This demonstrates that the PAR method also perform well in the detection of small objects.

Improvement in reconstruction precision using the PAR method: projection error analysis
Figure 2 shows the plots of the projection error as a function of iteration number, of which the green, blue and red lines represent the STR, LMR and PAR methods, respectively.The plotted curves in Figs.2(a) and 2(c) illustrate objects located at the center of the test domains, with object diameters of 4 mm and 10 mm, respectively.The curves in Figs.2(b) and 2(d) illustrate situations in which objects of 4 mm and 10 mm were located near the peripheral region, respectively.The regularization value was chosen empirically as 1 for the center and 10 for the edge in the STR method.In the LMR approach 1 and 100 was chosen empirically as their initial regularization values for the centered and edged objects.
In Figs.2(a)-2(d), the final projection errors for the STR method were 0.14, 0.13, 0.26 and 0.42, respectively, and the values for the LMR approach were 0.18, 0.17, 0.24 and 0.35, respectively, In contrast, the final projection errors using the PAR technique were 0.08, 0.09, 0.09 and 0.10, respectively, which are much smaller than those obtained with either STR or LMR in each case (see Fig. 2).Similarly, the CNR values were higher when applying PAR in all four cases.These indicate that the PAR method can achieve much better reconstruction precision compared to the STR and LMR methods.

Reduction in sensitivity to noise by the PAR method
Figures 3 and 4 show the images of the centered and edged objects reconstructed with STR, LMR and PAR at 1%，5% and 10% noise levels, respectively.The CNR values for these reconstructed images are listed in Table 1.As seen from the first rows of Figs. 3 and 4, STR, LMR and PAR techniques can all exhibit good image quality for objects located at the center and near the peripheral regions at a 1% noise level, but the reconstructed images for the PAR method have fewer artifacts and a higher contrast than the images obtained using the STR and LMR techniques.In response, the larger CNR values for the PAR method can be seen in Table 1.We also noticed that the image quality for the STR and LMR techniques dramatically decrease as noise is increased.When the added noise is 5% in amplitude and 5 degrees in phase, the centered objects can hardly be recognized from the background.The lower CNR values shown in Table 1 reveal the difficulty that these both methods had in reconstructing the targets.In comparison, the reconstructed images using the PAR approach still have a satisfactory image quality, even if the noise is increased to a 10% level.Clear targets and high object contrast can be observed in the reconstructed images in Figs.3(h) and 3(i) and Figs.4(h) and 4(i).This indicates that the PAR method has lower sensitivity to noise compared with the STR and LMR methods.
Furthermore, the quality of reconstructed images shown in Fig. 3 for both optical parameters is somewhat decreased compared to the image quality in Fig. 4.This is due to the fact that the objects in Fig. 3 are far away from the surface of the imaging fields.However, it is worth noting that the objects at the center can still be clearly recognized using our PAR method even at several noise levels.

Detectability using the PAR method
The CNR graphs obtained using PAR for three object positions (at the center, at 23 mm and 3 mm from the edge of the imaging fields) are shown in Figs.1(a), 1(b) and 1(c), respectively.
In the case of objects embedded at the centers of the imaging fields (i.e.Fig. 5(a)), it can be seen that the minimum detectable sizes and contrast values are 4.0 mm and 1.1, respectively.As the object moves from the center to the peripheral region, the fractional areas with CNR more than 3 increase from 501 to 530 to 581, which means that the range of detectable objects for the PAR is increased in terms of both object size and contrast, as observed in Figs.5(b) and 5(c).The results reported here show that smaller objects can be detected, compared with those reported by Pogue and Song in [16] and [29], so this method exhibits better detectability for small objects.This also indicates that the PAR method is effective in improving the spatial resolution and contrast resolution (i.e., the contrast required for detection at low resolution).The adaptive regularization method was applied in the reconstruction.

Discussion
In this paper, we have made a comparison among the STR, LMR and PAR methods.The results demonstrate that the PAR method significantly outweighs the STR and LMR approaches in improving image quality.The results also suggest that the PAR method exhibits a better performance in enhancing the accuracy of the reconstruction.In contrast, when applying both the STR and LMR techniques under various imaging conditions, the choice of the regularization values requires a trial-and-error process, which limits the reconstruction efficiency.In contrast, the PAR technique can adaptively choose an appropriate regularization value for different imaging conditions; therefore, PAR can obtain reconstructed objects more efficiently.
In our simulations, the same levels of random noise were added to all measurements, as is commonly done [22].However, it needs to be pointed out that in an actual imaging situation, the signal-to-noise ratio (SNR), not the noise level itself, is lower for detectors that are distant from the source.
A more realistic imaging system there might include systematic noise in the measured data besides the random noise.Pogue et al. noted that the total noise in their imaging system caused a measurement error approximately equal to a shot noise level of 3%~4% [16].Based on their study, we reconstructed the object images in the simulations for sensitivity to noise (i.e.Section 4.1.3)using the three noise levels of 1%, 5% and 10%.The data with 5% and 10% noise can be hypothesized as approximating the measurements obtained from medium and high noise-distributed imaging systems, respectively.The reconstruction for the 1% noise level should be viewed as a best case scenario and used as a comparison for the other reconstructed images (the effect of the 1% noise level on the measured data basically corresponds to the effect of a 12 mm-diameter inclusion in homogenous media).Three regularization methods for STR, LMR, and PAR have been applied in this study to perform reconstructions for these three noise levels and their reconstruction results have been presented in Figs. 3 and 4 for centered and edged objects, respectively.It was surprising that STR and LMR performed poor object reconstructions at the 5% noise level.In contrast, the objects reconstructed using PAR can be recognized clearly even at the 10% noise level.In particular, the reconstruction of a peripheral object displays a much lower sensitivity to noise.
Comparisons of the results obtained from the three regularization methods in three noise levels indicate that noise levels significantly affect the quality of reconstructed images, but the PAR method is the least influenced.It can also be observed that the noise level significantly impacts the projection error.When the noise level was increased from 1% to 10%, the projection error also increased.Hence, in the PAR approach, the large projection error leads to little change in the regularization value, which thus remains high throughout the whole reconstruction process.Nevertheless, in a high noise situation, constraining the minimum regularization to obtain a higher value can achieve optimal reconstructed images [15,26].Again, compared to the reconstructed images obtained using STR and LMR, the PAR images (Figs.3(i) and Figs.4(i)) demonstrate that a high regularization value is crucial for obtaining optimal images.In addition, we found that in this high noise situation, using a constant regularization value of 1/2 will produce almost the same quality of reconstruction as the images obtained by PAR.In fact, STR would not have performed so poorly under this high noise condition as the initial results indicated (Figs. 3 and 4) except for the fact that different regularization parameters for amplitude and phase were placed in the Hessian matrix (this can be seen in Ref. 19).As we show in Fig. 6, a significantly better reconstructed image can be obtained using STR when a regularization value of 1/2 is used for the terms associated with amplitude and phase in the Hessian matrix.
Simulations with non-randomly distributed (systematic) noise were also performed in this study.The simulations were designed so that part of the measured photon density was directed to 3 different detectors to simulate situations in which the contact between the detector and an object were either poor or ideal.Ten (10), 30 and 50% of the original photon density (OPD) from these three detectors was used as the measured data to simulate situations in which contact is poor.A measurement using 90% of the OPD was chosen to represent an ideal contact between the object and the detector.Note that the measured data for additional 13 detectors were respectively equal to their own OPD.A 12 mm diameter object was located at the edge of imaging field and the object image was reconstructed (with 1% noise) using STR, LMR and PAR under the same conditions.Figure 7 shows several reconstructed images for poor (30% of the OPD) and ideal (90% of the OPD) detector contact cases.As seen from the reconstructed images in Fig. 7, the degree of detector contact with the object significantly affects the image quality.When the contact was poor, the reconstructed images were all of very poor quality regardless of which of the three methods was used.Also, in these situations more artifacts were found near the sources.However, when the detector was in good contact with the object, the image quality for these three methods was greatly improved.Nevertheless, we found that our PAR method provided a better reconstruction than either STR or LMR.

Fig. 1 .
Fig. 1.Graphs of the CNR distribution of the reconstructed images using the (a) STR method, (b) LMR method, and (c) PAR method.A circular object located in the center of this imaging field was moved along the y-axis from 0 to 40 mm in 2 mm steps.At each position, the object size was increased from 1mm to 20 mm with 1mm intervals.

Fig. 2 .
Fig. 2. Comparisons of projection error and iteration number among the STR, LMR and PAR techniques.(a) and (c) illustrate objects 4mm and 10mm in diameter located at the center.The object sizes in (b) and (d) are the same as for (a) and (c) but are located near the edge.The CNR values indicate the image qualities of reconstructed objects with the three regularization methods.

Fig. 3 .
Fig. 3. Comparisons of reconstructed μ a and μ s ′ images with three regularization methods at different noise levels.The 12-mm-diameter objects are located in the center of the imaging fields.(a)-(c) for STR, (d)-(f) for LMR, and (g)-(i) for PAR.The reconstructions in the first row are with 1% noise, the second row with 5% noise, and 10% noise for the last row.

Fig. 4 .
Fig.4.The imaging conditions are the same as in Fig.3, except that the objects are located at the edge of the imaging fields.

Fig. 5 .
Fig. 5. Graphs of the CNR distribution of reconstructed images for three object positions using PAR: (a) central region and (b) 23 mm from the boundary and (c) 3 mm from the boundary.The adaptive regularization method was applied in the reconstruction.

Fig. 6 .
Fig. 6.STR reconstruction comparison for two regularization schemes: different regularization parameters for the amplitude and phase were placed in the Hessian matrix ((a) and (c)) and a constant regularization value of 1/2 was placed in the Hessian matrix ((b) and (d).The centrally located object appears in (a) and (b), and the peripheral one in (c) and (d).The object size is 12 mm in diameter.

Fig. 7 .
Fig. 7. Simulation comparisons of reconstructed μ a and μ s ′ images with three regularization methods in poor (i.e.30% OPD, (OPD represents the original photon density.))and ideal (90% OPD) detector-object contact cases.Three measurements were chosen, and 30 and 90% of the OPD of each measurement was respectively considered as the measured photon density.A 12-mm-diameter object was located at the edge of the imaging fields.(a) and (b) for STR, (c) and (d) for LMR, and (e) and (f) for PAR.

Table . 1
. The CNR values for the reconstructed μ a and μ s ′ images shown in Figs.3 and 4.