Open Access
Add to BibSonomyAbstract
Diffuse optical tomography (DOT) has much lower sensitivity in deep tissues than in superficial tissues, which leads to poor depth resolution. In this paper, a layer-based sigmoid adjustment (LSA) method is proposed for reducing sensitivity contrast in the depth dimension. Using this method, differences in image quality between layers can be effectively reduced. As a result, positioning errors of less than 3 mm can be obtained in the depth dimension for all depths from -1 cm to -3 cm.
©2007 Optical Society of America
1. Introduction
In 1977 Jobsis [1] showed that near-infrared light (wavelengths of 650nm to 950nm) could penetrate several centimeters into biological tissues, which made it possible to “see” activations non-invasively as they happened in deep tissues. Among the methods developed to reconstruct images of activations in biological tissues [2–5], diffuse optical tomography (DOT) has become important [6–9]. In DOT, near-infrared light is introduced into tissues through optical fibers. Near-infrared photons are then scattered many times within the biological tissues. These can be considered as random walks or as diffuse propagations. During propagation, some photons are absorbed by absorbers such as hemoglobins (including oxy-hemoglobin and deoxy-hemoglobin) and water. Afterwards, the photons that emerge are collected by detectors placed several centimeters away from the sources.
It has been well established that a detected signal from a source-detector pair can be influenced by changes in hemoglobin concentrations over a large area [10–12], which can be defined as the field of view (FOV) of this source-detector pair. In addition, detected signals are much more sensitive to activations arising near the sources and detectors than from those at other positions within the FOV. Reflectance diffuse optical tomography (rDOT) has been used for adult human brains, and its sensitivity has been shown to drop off exponentially in the depth dimension in tissues deeper than several millimeters [13]. Since reconstructed activations tend to be located in regions having higher sensitivity, rDOT has poor depth resolution, especially in deep tissues.
In practice, time-domain (TD) measurements can achieve better depth resolution than continuous wave (CW) ones [14–16]. However, time-domain measurements are limited by expensive instruments and low spatio-temporal resolution. This paper concentrates on improving the depth resolution of CW measurements. This would have the effect of making CW-DOT an appropriate tool for imaging adult brains by proving high spatial, temporal and depth resolution.
Using information previously obtained from anatomical magnetic resonance imaging (MRI) can achieve better depth localization of rDOT [13, 17]. In such methods, a 3-D model of the brain was obtained using high resolution MRI. Then voxels belonging to cortical tissues were separated and images were reconstructed based on these cortical-tissues-only voxels. Under the constraint obtained by using the cortical-tissues-only voxels, reconstructed activations can only be located on the superficial surfaces of cortical tissues, because within the cortical tissues, the sensitivity of rDOT dropped off exponentially.
Fig. 1. The multicentered geometry used in this paper. The sources and detectors are indicated by dots and open circles, respectively. The imaging region is indicated by the block.
Spatial variant regularization (SVR) is a family of methods in which different regularization parameters are used for different voxels [18, 19]. SVR methods can achieve constant spatial resolution in the imaging region, and thus benefit the image quality of DOT.
In this paper for the first time to our knowledge, a layer-based sigmoid adjustment (LSA) method is proposed to balance the sensitivity contrast in the depth dimension by a direct adjustment in the forward matrix. The LSA method is based on the fact that sensitivity drops off significantly in the depth dimension in tissues deeper than several millimeters, but at the same time sensitivity contrasts within same layer usually remain unapparent. To evaluate the proposed LSA method, simulations were performed using a semi-infinite model.
2. Methods
2.1. Multicentered geometry
The multicentered geometry [20] is shown in Fig. 1. The side of the equilateral hexagon is 4 cm. The multicentered geometry has 7 sources in the central part of the hexagon; and around these sources 24 detectors are placed. In Fig. 1, source 2 is located on the line defined by source 1 and detector 1, and is 1.15 cm away from source 1. Detector 13 is located 0.58 cm below the top edge and is equidistant from detectors 1 and 2. Sources 3 to 7 and detectors 14 to 24 can be placed in the same way as source 2 and detector 13, respectively.
In this multicentered geometry, a signal from any source can be received by all of the detectors, so a total of 168 measurements can be achieved. The imaging region in the x-y plane is 6.0×6.0 cm2, as indicated by the block in Fig. 1. In this paper we chose depths from -1.0 cm to -3.0 cm, based on the fact that the distance from the surface of the scalp to the surface of the cortical tissues of an adult human would be no less than 1.0 cm. The voxel size in the reconstructed images is 1×1×1 mm3.
Fig. 2. The overall sensitivity (a) within the x-y slice with z = -2.0 cm , and (c) within the x-z slice with y = 0 cm . The profiles (b) with y = 0 cm in the x-y slice and (d) with x = 0 cm in the x-z slice are shown in the first column. The corresponding (e) x-y slice and (g) x-z slice, together with the corresponding profiles in the (f) x-y slice and (h) x-z slice after adjustment with LSA parameter of 400.0 are shown in the second column.
2.2. Forward model
The propagation of near-infrared photons in biological tissues can be described by a diffusion equation [21]. In practice, the analytic solution of diffusion equations can be obtained under the semi-infinite assumption [21–23]. The scattering portion of the analytic solution is hard to determine using CW instruments. Fortunately, the scattering portion remains nearly constant during any single simulation, so we considered only changes in the absorption portion, which was defined as in Ref. [24]:
In which ΔOD is the change in optical density. Φ0 is the simulated photon fluency in a semi-infinite model, and Φpert is the simulated photon fluency with the absorbers included. The second part of Eq. (1) is known as the modified Beer-Lambert law, in which Δμa(r) is the change in absorption coefficient, and L(r) is the effective average path length of the photons. In this paper, Gaussian noise was added to Φ0 and Φpert to achieve an approximate signal to noise ratio (SNR) of 1000.
The discrete form of the modified Beer-Lambert law can be written as y = Ax , in which y contains the measured changes in optical density from all the measurements; x is the unknown changes in absorption coefficient in all of the voxels; and matrix A contains the mapping from the absorption space to the measurement space. The elements in matrix A can also be viewed as representations of the sensitivity of the measurements to changes in the absorption coefficients. The sum of matrix A by rows indicates the overall sensitivity in each voxel, as shown in Fig. 2.
2.3. Layer-based sigmoid adjustment method
As can be observed from the first column of Fig. 2, the sensitivity dropped off sharply as depth increased, while within the same layer the contrast in sensitivity was not significant. To smooth the sensitivity in the depth dimension, a layer-based sigmoid adjustment (LSA) was performed using the following steps.
First, the sigmoid function was defined by the following equation:
In Eq. (2), a is the LSA parameter and the plots of the sigmoid function for several a ’s are shown in Fig. 3. The range of β was chosen subjectively and remained unchanged throughout the simulations.
Second, assuming that the reconstructed images contained L layers, β was then sampled equally to get β 1,...,βL.
Last, γ(β 1) was the LSA coefficient for the bottom layer and γ(βi) was the LSA coefficient for the ith layer. The forward matrix A was then adjusted according to the following equation:
here D is a diagonal matrix. Equal LSA coefficients were assigned to every voxel within the same layer.
Fig. 3. The plots of the sigmoid function with LSA parameters of 200 and 400. β was sampled equally from -3.0 to 5.5, and the corresponding γ was calculated and used as the LSA coefficients.
2.4. Image reconstruction
We reconstructed images using the following three Tikhonov regularization methods [25, 26], the first one of which is the DOT method without the LSA adjustment; the second one is defined as the LSA method; and the third one is the SVR method taking D -2 as the spatial variant parameter.
Here A and à are the forward matrices before and after the LSA adjustment. x̂ is the reconstructed image; y is the simulated ΔOD . λa , λb and λc are the regularization parameters for the three equations.
Equation (5) can be rewritten as:
Equation (6) can be rewritten as:
Since Eq. (8) is equivalent to:
we can see that the LSA method and the SVR method are formally similar. In this paper, we compared the image quality of the LSA method to the SVR method.
In Eq. (4) λa = λ ∙ s max, parameter λ is generated by the L-curve method and matrix A [27, 28], s max is the maximal singular value of matrix A . In this paper, we compared the image quality of Eqs. (5) and (6) with various λb and λc to select the best regularization parameters:
Here is generated by the L-curve method and matrix Ã, s̃ max is the maximal singular value of matrix Ã.
2.5. Evaluation of image quality
To calculate the contrast to noise ratio (CNR), the reconstructed image was divided into two parts, a region of interest (ROI) and a region of background (ROB). The ROI was the same region as that of the located objects, while the ROB was defined as the residual region in the reconstructed image. After the mean values and standard deviations of the reconstructed absorption coefficients over the ROI and ROB were calculated, the CNR value was calculated using the following equation [29]:
where wROI is the division of the area of the ROI by the area of the whole image, and wROB is defined as wROB = 1-wROI ∙ μROI and μROB are the mean values of the objects and background regions in the reconstructed images. σROI and σROB are the corresponding standard deviations.
The detected objects are defined as voxels with absorption coefficients above the half-maxima in the reconstructed images. The centers of the detected objects were calculated by averaging the coordinates of the voxels within the objects. Then the positioning errors (PE) became the differences between the real positions of the objects and the centers of the detected objects. In this paper, only the PE in the z dimension is evaluated.
3. Simulations
The background absorption and reduced scattering coefficients used in the semi-infinite model are 0.1 cm-1 and 10 cm-1, respectively. The perturbations in absorption coefficients, which can also be viewed as objects, are simulated by spheres with radii of 3 mm and absorption coefficients of 0.3 cm-1.
To select the best regularization parameters for the LSA method and the SVR method, objects were located from -1.0 cm to -3.0 cm in depth, and LSA parameters were varied from 1.0 to 400.0. The intervals of the depths and the LSA parameters were 1 mm and 10.0, respectively. For each depth and LSA parameter an image was reconstructed and the corresponding CNR value was calculated. The CNR maps are shown in Fig. 4. The quotients obtained by dividing standard deviations of CNR values over all depths by corresponding mean values are shown in Fig. 5.
For an illustrational analysis of the improvement in image quality by the LSA method, the real images with objects located at depths of -1.0 cm, -1.3 cm, -1.6 cm, -1.9 cm, -2.1 cm, -2.4 cm, -2.7 cm and -3.0 cm are shown in the first column of Fig. 6. The reconstructed images without and with the LSA adjustment are shown in the second and the third columns of Fig. 6. The reconstructed images using the SVR method are shown in the last column of Fig. 6. The CNR and PE values for the reconstructed images are listed in Table 1. The LSA parameters used to generate matrices D for the LSA method and for the SVR method are 400.0 and 20.0, respectively.
Fig. 4. The CNR maps for the LSA methods with (a) λb = λ ∙ s̃ max, (b) λb = λ ∙ s̃ max, and (c) λb = ∙ s̃ max . The CNR maps for the SVR methods with (d) λc = λ ∙ s max, (e) λc = λ ∙ s̃ max, and (f) λc = ∙ s̃ max.
Fig. 5. The quotients obtained by dividing standard deviations of CNR values over all depths by the corresponding mean values for (a) the LSA method and (b) the SVR method. The red lines are for {λb, λc} = λ ∙ s max. The blue lines are for {λb, λc} = λ ∙ s̃ max. The black lines are for {λb, λc} = ∙ s̃ max.
4. Results
4.1. Reduction in sensitivity contrast by the LSA method
The proposed LSA method is effective in reducing the sensitivity contrast in the depth dimension. Note from Figs. 2(c) and 2(d) that the maximal level of sensitivity at a depth of -3.0 cm is about 1/154 of the top layer before adjustment. In DOT, the measured signals are influenced by changes in absorption coefficients both in the superficial layers and in the deep layers. The much smaller sensitivities in the deep layers have trivial weights in the reconstruction of images; so reconstructed objects tend to be found in superficial layers, and large PE can be found in the depth dimension.
However, the maximal sensitivity at a depth of -3.0 cm is about 1/4 of the maximal sensitivity in the whole imaging region with an LSA parameter of 400.0, as shown in Figs. 2(f) and 2(h). The low sensitivity contrast balances the importance of different layers in the reconstruction of images and helps to locate the activations at the right depth.
4.2. Selection of best regularization parameters
The LSA method can achieve the best image quality for λb = λ ∙ s max because the brighter area, which indicated larger CNR values, is larger in Fig. 4(a) than that in Figs. 4(b) and 4(c). The same conclusion can be reached through comparing the three lines in Fig. 5(a) since the LSA method taking λb = λ ∙ s max as the regularization parameter can achieve the most uniform image quality. Comparison of the three lines in Fig. 5(b) shows that the SVR method, using λc = λ ∙ s max as the regularization parameter, can achieve slightly better image quality than using the other two as the regularization parameter.
From the previous statement we can conclude that the best regularization parameters for the LSA and SVR methods are λb = λc = λa. In the rest of this paper, the regularization parameters used in the image reconstructions for the LSA and SVR methods are λb = λc = λa
We also chose LSA parameters of 400.0 for the LSA method and of 20.0 for the SVR method to reconstruct the images in Fig. 6 since the most uniform image quality can be reached using these numbers for the two methods, as can be seen in the red lines in Fig. 5.
Fig. 6. The (a) real images, reconstructed images (b) without and (c) with applying the LSA method, and (d) with the SVR method. The objects were located at depths of -1.0 cm, -1.3 cm, -1.6 cm, -1.9 cm, -2.1 cm, -2.4 cm, -2.7 cm and -3.0 cm (from top to bottom).
Table 1. The CNR and PE values for the reconstructed images shown in Fig. 6.
4.3. Improvement in image quality by the LSA method
The reconstructed images without any adjustment in the forward matrix are clearly biased in favor of the superficial layers. This can be observed from the second column of Fig. 6. As a result, high CNR values and small PE values can be seen in Table 1 for depths from -1.0 cm to -1.3 cm. At the same time much smaller CNR values and much larger PE values can be found for depths from -1.9 cm to -3.0 cm. This indicates that DOT without any adjustment has poor performance in the depth dimension.
In contrast, the images reconstructed by the SVR method from -1.3 cm to -2.1 cm are clearly biased in favor of the deep layers, as can be seen in the last column of Fig. 6. In particular, noise in the deep layers is emphasized so that the reconstructed images for -1.0 cm are too blurred to recognize the reconstructed objects.
Unlike the previous two methods, the LSA method can achieve a satisfactory image quality for both superficial objects and deep objects. As can be seen in Fig. 6, the LSA method can achieve much better image quality than the DOT without any adjustment for deep objects, and can achieve much better image quality than the SVR method for superficial objects. In particular, PE values of less than 3.0 mm can be obtained for all depths from -1.0 cm to -3.0 cm by the LSA method, as can be seen in the fifth column of Table 1. This indicates a significant improvement in depth resolution by the LSA method.
Since near-infrared photons can only penetrate a few centimeters into biological tissues, DOT cannot image tissues deeper than the maximal depth that near-infrared photons can reach. Nevertheless, the LSA method can reduce the deviations in image quality in the depth dimension within the visual field of DOT. This will benefit the detection of human brain activations.
4.4. Relationship between the LSA method and the SVR method
The poor image quality of the SVR method for superficial objects is mainly because the pre-multiplied D on the right-hand side of Eq. (9) overemphasizes deep layers. However, the LSA method can moderately balance the superficial layers and deep layers, so better image quality can be achieved by the LSA method. This indicates that the LSA method and the SVR method are different methods, even though they are formally similar.
The LSA method directly balances the sensitivity contrast in the depth dimension. The SVR method uses smaller regularization parameters for regions with lower sensitivity. The results in this paper indicate that the former is the better method for improving depth resolution of DOT.
5. Discussion
In this paper, the LSA parameters used are larger than 1.0, so the LSA method will multiply each element in the forward matrix A by a number larger than 1.0. The enlarged elements in the adjusted forward matrix à will underestimate changes in absorption coefficients.
The underestimation of changes in absorption coefficients is a common problem in DOT. In this paper, the underestimation in absorption coefficients has little influence on the recognition of reconstructed objects. However, this problem is still an important one and should be left to be resolved by future advances.
6. Conclusions
Large contrasts in sensitivity between deep and superficial tissues are believed to be the main reason for the poor depth resolution of DOT. In this paper, we demonstrated for the first time to our knowledge that the depth resolution of DOT could be significantly improved by direct adjustments in the sensitivity matrix. Rather than providing a more nearly accurate forward model for solving the inverse problem, the proposed layer-based sigmoid adjustment (LSA) method balances the sensitivity in the depth dimension. Simulations performed in a semi-infinite model indicated that the LSA method could effectively reduce the sensitivity contrast between the deep and superficial tissues, and positioning errors of less than 3 mm could be achieved for depths from -1 cm to -3 cm. Further analysis of the simulation results showed that the proposed LSA method had a better image quality than that obtained by the SVR method.
Acknowledgments
The authors thank the reviewers for their constructive suggestions. They also thank Drs. Rhoda E. and Edmund F. Perozzi for editing and English language assistance. This work was partially supported by the Natural Science Foundation of China, Grant Nos. 30425004 and 60121302, and the National Key Basic Research and Development Program (973), Grant No. 2004CB318107.
References and links
1. F. Jobsis, “Non-invasive infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory Parameters,” Science 198,1264–1267 (1977). [CrossRef] [PubMed]
2. M. A. Franceschini, V. Toronov, M. E. Filiaci, and E. Gratton, “On-line optical imaging of the human brain with 160-ms temporal resolution,” Opt. Express 6,49–57 (2000). [CrossRef] [PubMed]
3. S. Fantini, M. A. Fanceschini, and E. Gratton, “Non-invasive optical mapping of the piglet brain in real time,” Opt. Express 4,308–314 (1999). [CrossRef] [PubMed]
4. G. Taga, Y. Konishi, A. Maki, T. Tachibana, M. Fujiwara, and H. Koizumi, “Spontaneous oscillation of oxy- and deoxy- hemoglobin changes with a phase difference throughout the occipital cortex of newborn infants observed using non-invasive optical topography,” Neurosci. Lett. 282,101–104 (2000). [CrossRef] [PubMed]
5. V. Ntziachristos and B. Chance, “Probing physiology and molecular function using optical imaging: applications to breast cancer,” Breast Cancer Res. 3,41–46 (2001). [CrossRef] [PubMed]
6. A. M. Siegel, J. J. A. Marota, and D. A. Boas, “Design and evaluation of a continuous-wave diffuse optical tomography system,” Opt. Express 4,287–298 (1999). [CrossRef] [PubMed]
7. S. R. Hintz, D. A. Benaron, A. M. Siegel, A. Zourabian, D. K. Stevenson, and D. A. Boas, “Bedside functional imaging of the premature infant brain during passive motor activation,” J. Perinat. Med. 29,335–343 (2001). [PubMed]
8. A. Bluestone, G. Abdoulaev, C. Schmitz, R. Barbour, and A. Hielscher, “Three-dimensional optical tomography of hemodynamics in the human head,” Opt. Express 9,272–286 (2001). [CrossRef] [PubMed]
9. E. M. C. Hillman, J. C. Hebden, M. Schweiger, H. Dehghani, F. E. W. Schmidt, D. T. Delpy, and S. R. Arridge, “Time resolved optical tomography of the human forearm,” Phys. Med. Biol. 46,1117–1130 (2001). [CrossRef] [PubMed]
10. D. A. Boas, J. P. Culver, J. Stott, and A. K. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult head,” Opt. Express 10,159–170 (2002). [PubMed]
11. E. Okada and D. T. Delpy, “Near-infrared light propagation in an adult head model. I. Modeling of low-level scattering in the cerebrospinal fluid layer,” Appl. Opt. 42,2906–2914 (2003). [CrossRef] [PubMed]
12. E. Okada and D. T. Delpy, “Near-infrared light propagation in an adult head model. II. Modeling of low-level scattering in the cerebrospinal fluid layer,” Appl. Opt. 42,2915–2922 (2003). [CrossRef] [PubMed]
13. D. A. Boas and A. M. Dale, “Simulation study of magnetic resonance imaging-guided cortically constrained diffuse optical tomography of human brain function,” Appl. Opt. 44,1957–1968 (2005) [CrossRef] [PubMed]
14. A. Kienle, M. S. Patterson, N. Dognitz, R. Bays, G. Wagnieres, and H. van den Bergh, “Noninvasive determination of the optical properties of two-layered turbid media,” Appl. Opt. 37,779–791 (1998). [CrossRef]
15. J. Steinbrink, H. Wabnitz, H. Obrig, A. Villringer, and H. Rinneberg, “Determining changes in NIR absorption using a layered model of the human head,” Phys. Med. Biol. 46,879–896 (2001). [CrossRef] [PubMed]
16. M. Kohl-Bareis, H. Obrig, J. Steinbrink, J. Malak, K. Uludag, and A. Villringer, “Noninvasive monitoring of cerebral blood flow by a dye bolus method: separation of brain from skin and skull signals,” J. Biomed. Opt. 7,464–470 (2002). [CrossRef] [PubMed]
17. B. W. Pogue and K. D. Paulsen, “High-resolution near-infrared tomographic imaging simulations of the rat cranium by use of a priori magnetic resonance imaging structural information,” Opt. Lett. 23,1716–1718 (1998). [CrossRef]
18. B. W. Pogue, T. O. McBride, J. Prewitt, U. L. Osterberg, and K. D. Paulsen, “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt. 38,2950–2961 (1999). [CrossRef]
19. J. P. Culver, A. M. Siegel, J. J. Stott, and D. A. Boas, “Volumetric diffuse optical tomography of brain activity,” Opt. Lett. 28,2061–2063 (2003). [CrossRef] [PubMed]
20. Q. Zhao, L. Ji, and T. Z. Jiang, “Improving performance of reflectance diffuse optical imaging using a multicentered mode,” J. Biomed. Opt. 11, (2006) (in press). [CrossRef] [PubMed]
21. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. Eng. 15,R41–R93 (1999). [CrossRef]
22. X. Cheng and D. A. Boas, “Diffuse optical reflectance tomography with continuous-wave illumination,” Opt. Express 3,118–123 (1998). [CrossRef] [PubMed]
23. D. A. Boas, A. M. Dale, and M. A. Franceschini, “Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy,” NeuroImage 23,275–288 (2004). [CrossRef]
24. D. T. Delpy, M. Cope, P. van der Zee, S. Arridge, S. Wray, and J. Wyatt, “Estimation of optical pathlength through tissue from direct time of flight measurement,” Phys. Med. Biol. 33,1433–1442 (1988). [CrossRef] [PubMed]
25. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. Eng. 15,R41–R93 (1999). [CrossRef]
26. D. A. Boas, A. M. Dale, and M. A. Franceschini, “Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy,” NeuroImage 23,275–288 (2004). [CrossRef]
27. L. Wu, “A parameter choice method for Tikhonov regularization,” Electron. Trans. Numer. Anal. 16,107–128 (2003).
28. P. C. Hansen and D. O’Leary, “The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems,” SIAM J. Sci. Comput. 14,1487–1503 (1993). [CrossRef]
29. X. Song, B. W. Pogue, S. Jiang, M. M. Doyley, and H. Dehghani, “Automated region detection based on the contrast-to-noise ratio in near-infrared tomography,” Appl. Opt. 43,1053–1062 (2004). [CrossRef] [PubMed]
