In the last two decades, many theoretical and numerical methods have been developed to model photon migration in layered tissue structures. 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 One well-established method is Monte Carlo (MC) simulation, which can be used to characterize the light propagation in complex tissue structures.7, 10, 11 However, the MC method is very computationally intensive. Analytical methods have been developed by many research groups to model light propagation in layered tissue structures. The publications in this regard have mostly focused on estimation of the dependence of the surface reflectance measurements on the optical properties of the two tissue layers.1, 2, 13 The widely referenced paper by Kienle 3 provided solutions of time- and frequency-domain diffusion equations for two-layer tissue structures. These solutions could be used to estimate tissue optical properties and in some cases the layer thickness by using nonlinear fitting methods to the measured reflectance data.
To the best of our knowledge, there has been little effort to develop tomographic reconstruction algorithms suitable for layered tissue structures. In breast imaging, a semi-infinite medium is an accurate approximation when the patient is scanned in a supine position and the breast thickness is or more. However, the semi-infinite model is not accurate enough when breast tissue thickness is less than with the chest wall underneath. Under this condition, the chest wall contributes to the total reflectance measurements at the skin surface. The optical properties of breast tissue and the chest wall can be estimated using the two-layer analytical solution with the layer thickness provided by coregistered ultrasound. The accuracy of reconstructed target optical properties can be improved by implementing a layered structure in the imaging process as opposed to the semi-infinite tissue structure.15
In a recent publication, we discussed the development of the analytical solution for a two-layer tissue structure with a tilted interface to improve the estimation accuracy of the optical properties of breast tissue and the chest wall.12 In this paper, we extend the analytical solution of a two-layer tissue structure developed by Kienle 3 to tomographic imaging. The Born approximation and the dual-zone mesh scheme we developed earlier are used in the image reconstruction.16, 17 Simulations and phantom experiments are performed to validate the image reconstruction. Further, clinical examples are given to demonstrate the utility of this layer-structure-based tomography algorithm.
We briefly summarize the forward model of light propagation used in this paper in Sec. 2.1, and discuss our newly developed image reconstruction by combining the Born approximation with a dual-zone mesh scheme for the layered structure.
The diffusion equations were solved by Kienle 3 using the Fourier transform approach for a two-layer turbid medium with a semi-infinite second layer. We directly adopted this method for tomographic imaging and a brief description of the analytical solution is given here for the sake of completeness. A two-layer turbid medium is considered with upper layer thickness and a semi-infinite second layer. An infinitely thin beam of frequency-modulated light is incident perpendicular on the top surface of a two-layer turbid medium. The beam is scattered isotropically at the depth of , where the and are the absorption and reduced scattering coefficients of the first layer, and is in the range of for typical background breast-tissue absorption and scattering properties. The diffusion equation for the photon density wave at two layer structure is given as, are the diffusion coefficients, and are the absorption and reduced scattering coefficients, and are the fluence rates for the layers and 2, respectively. The modulation frequency of the source is , is the speed of light in the medium, and is the imaginary unit.
The diffusion equations were first transformed to ordinary differential equations with the use of a 2-D Fourier transformand are spatial frequencies in the and directions. The derived equations are then solved by using the following boundary conditions,1, 2. The results for at different tissue layers are , ; is the extrapolated boundary18 given by ; and depends on the refractive indices of the tissue layer and represents the fraction of photons internally diffused at the boundary.
The fluence rate in the spatial domain can be obtained by spatial Fourier inversion:is the Bessel function of zero order, and , where and are spatial coordinates of the source position. Thus, Eq. 11 combined with either Eq. 8 (for depth ), Eq. 9 (for depth ), or Eq. 10 (for ) provides the incident fluence rate at any depth inside the two-layer medium.
Image Reconstruction Using Two-Layer Tissue Model
In the Born approximation, the photon density wave originating from a source at and measured inside the tissue at , is a linear superposition of its incident (homogeneous) and the scattered (heterogeneous) parts given as . By using the first-order solution to the heterogeneous equation under the approximation that , we obtain , where represents detector position, represents the scattered field at due to the point scattered at , is the change of absorption coefficient at voxel at location , and is the diffusion coefficient of the homogeneous medium. Figure 1 illustrates the source position , field point , detector position , and used for computing the incident field.
When using the two-layer model, is calculated based on the depth of the voxel to which the fluence is being computed. For voxels in the three regions the fluence rates are represented by , and , and is replaced by or , depending on the voxel location in the layered structure. The Green function of the homogeneous equation at the detector is given as , which could be calculated again using spatial Fourier inversion from any one of Eq. 8 (if voxel depth ), Eq. 9 (if voxel depth ), and Eq. 10 (if ), depending on the voxel depth from the surface.
The dual-zone mesh scheme introduced by us earlier16 was used to segment the imaging volume into lesion and background regions designated as and , respectively. The weight matrix is calculated using the two-layer model as , where and are weight matrices for lesion and background volumes, respectively; and are the total absorption distributions of lesion and background volumes, respectively; and equals or , depending on the voxel location in the layered structure. The dual-zone mesh scheme significantly reduces the total number of voxels with unknown optical properties and dramatically improves the convergence of inverse mapping of target optical properties.16, 17
Materials and Methods
The imaging reconstruction algorithm developed in Sec. 2 was validated using forward data generated by the finite-element method (FEMLAB 6.2). The source-detector geometry is shown in Fig. 2 . There are total of nine sources and 10 detectors. The middle slot was used for a commercial ultrasound transducer. The geometry used in simulation was the same as that used in phantom and clinical measurements.
Two sets of simulations were performed. In the first set of simulations, a spherical target of diameter was embedded in the layered structure at a depth from the surface. The target center was placed at , and its absorption coefficient was varied from . Forward data was generated from two-layer structures whose optical properties were in the range of typical breast tissue and chest wall. For each set of target data, a reference data set was also generated using the same layered structure with no target in the medium. Perturbation data were used for imaging reconstruction of the target. Two sets of imaging reconstructions were performed. One set of images was reconstructed from known optical properties of the two-layer structure, while another set was obtained by using fitted optical properties from the reflectance data. The reconstruction results for the first set provided the ideal condition when no fitting error for the two-layer optical properties was considered. The second reconstruction set is relevant to realistic clinical experiments, where the first- and second-layer optical properties must be estimated from the reflectance measurements. For this set of simulations, image reconstruction was also performed by using a semi-infinite model and a two-layer model with true background optical properties and fitted background properties.
In the second set of simulations, the absorption coefficient of the -diam target was fixed at and the second-layer background optical absorption was varied from at intervals to simulate a range of chest-wall optical absorptions. For each set of background optical properties, forward data were generated with and without target. Image reconstruction was performed by using known background optical properties and fitted background optical properties. This set of simulations was intended to evaluate the effect of fitting error on reconstructed target absorption.
Phantom and Clinical Experiments
The experimental validation of the imaging algorithm was performed using an existing frequency-domain near-IR (NIR) system.19 It consisted of laser diodes of 780- and wavelengths. The outputs of the laser diodes were amplitude modulated at and optically switched to nine positions on a hand-held probe. Ten optical fibers of diameter were coupled to 10 parallel detection channels. Solid silicone phantoms were used to emulate the second layer over which intralipid solution was used for the top layer. All the phantom experiments were performed using a flat interface between the layers. The measurements were first made on simple two-layer structures with various layer depths and used as the reference data. The measurements were repeated for the same structures with solid target phantoms immersed in the intralipid layer at known depths. Two types of target phantoms were used; one with low absorption to emulate a benign lesion and one with a high absorption to emulate a cancerous lesion in the breast.20, 21 The layer thickness, target position, and target size were known in these phantom measurements as opposed to the clinical cases were we would estimate these parameters using a coregistered ultrasound probe.
The imaging algorithm was demonstrated with clinical examples. The study protocol was approved by a local Institutional Review Board (IRB) committee. Signed informed consent was obtained from all patients who agreed to participate in the study. In clinical experiments, NIR optical source fibers, detector fibers, and a commercial ultrasound transducer were integrated on a handheld probe.19 Patients were scanned in a supine position, while multiple sets of optical reflectance measurements were made with coregistered ultrasound images at lesion locations and symmetric locations of the contralateral normal breast. The measurements obtained at normal locations were used to estimate background optical properties, which were used to compute a weight matrix for imaging.
Estimation of Two-Layer Optical Properties and Image Reconstruction
Two-layer optical properties were estimated from fitting reflectance measurements made on the layered structure to analytical solution using the Nelder-Mead simplex algorithm as described in our previous paper.12 This algorithm minimizes a scalar-valued nonlinear function of multiple independent parameters by using only the function values, and it is chosen due to its simplicity and good convergence in spite of the slow convergence rate. We chose various initial values and found that initial values are not critical as long as the function converges. For each case, the initial values were chosen as 0.05, 5.0, 0.05, and for , , , and , respectively. In the phantom and clinical experiments, the thickness of the upper layer was estimated by using the coregistered ultrasound image. The nonlinear regression routine used here, fminsearch, is available in MATLAB.
In simulations, the voxel size was chosen as in the and dimensions for fine mesh and for coarse mesh. The corresponding voxel sizes were 0.25 and in phantom studies, and 0.5 and for clinical data. In all experiments, the voxel size in depth was in the fine-mesh zone and in the coarse-mesh zone. The choice of fine-mesh size is based on SNR and also spatial sampling considerations. The fine-mesh volume for all experiments was chosen at least 2.5 times larger than the target size and centered at the target position. Clinical experiments were performed by fitting the two-layer optical properties of the contralateral breast using the layer depth provided by coregistered ultrasound images. In all experiments, the dual-mesh algorithm was used for inversion with either known target location or size in simulation and phantom experiments or a priori information provided by coregistered ultrasound in clinical experiments.
In the first set of simulations, targets of various absorption properties were embedded in a layered structure of upper layer thickness. The first and second layer known optical properties were and and and , respectively. The target absorption coefficients were chosen as 0.07, 0.10, 0.15, 0.20, 0.25, and , and the reduced scattering coefficient was kept at . Reference data were also generated for the two-layer structure in the absence of the target. Figure 3 shows an example of reconstructed target absorption maps of an absorber ( and ) embedded in the layered-medium when the known optical properties were used in the reconstruction. Figure 3a shows the reconstruction result of using a semi-infinite model for the two-layer structure and known first-layer background values of and . Figure 3b shows the two-layer-model reconstruction using known two-layer background values of , , , and . The reconstructed maximum absorption values were (62%) and (73%), respectively. About 11% improvement in reconstruction accuracy was achieved using the layer-based reconstruction under this ideal scenario. When using a semi-infinite model to fit background optical properties of the two-layer structure, and were obtained. Figure 4a shows the result of semi-infinite-model-based reconstruction and the reconstructed maximum absorption coefficient was [Fig. 4a]. When fitted two-layer optical properties of , , , and were used in two-layer-model-based reconstruction, the reconstructed maximum absorption coefficient was [Fig. 4b]. About 7% improvements in reconstruction accuracy were achieved even with fitting errors of the two-layer optical properties.
Figure 5 shows the reconstructed maximum absorption coefficient versus true target absorption value with known and fitted optical properties of two-layer model (black star and purple cross lines) and the semi-infinite approximation for the layered structure (green triangle and red square lines), respectively. It is interesting to note that the two-layer model improves the reconstruction accuracy only at higher contrast end when known optical properties of two-layers are used. This is due to the use of a more accurate weight matrix for the two-layer structure. Because both semi-infinite and two-layer models use the linear Born approximation, the reconstructed target absorption is reduced at higher contrast. When the fitted two-layer optical properties are used, the improvement is distributed at both high and low contrast. The improvement at low contrast is due to the use of more accurate first-layer background properties. Because reconstructed target , a higher fitted background absorption using the semi-infinite model yields a higher target compared with using the two-layer model. This problem is less pronounced when the target to background ratio is higher.
In the second set of simulations, the first layer optical properties and the second layer reduced scattering coefficient were fixed at , and , and . However, the second layer background was varied from to simulate a range of chest-wall absorption. Table 1 lists the true and fitted two-layer optical properties and the corresponding reconstructed maximum target absorption coefficient. Target reconstruction results using true and fitted background properties for each set of background values are within 13%. The small difference results from the fitting error of first-layer .
Simulation results using the semi-infinite model and the two-layer model with fitted background values. Target center depth C=0.9 cm. Target radius R=0.5cm. Second layer depth D=1.5 cm. Refractive index (n 1=1.33, n 2=1.4).
|Two-layer Model with True Background Value||Two-layer Model with Fitted Background Value|
|True Background Value [μa1μs1′μa2μs2′] (cm−1)||Reconstructed Maximum μa (cm−1)||Fitted Background Value [μa1μs1′μa2μs2′] (cm−1)||Reconstructed Maximum μa (cm−1)|
|[0.02 7 0.02 7]||0.1525 (76%)||[0.0335 6.9 0.0174 7.05]||0.1369 (68%)|
|[0.02 7 0.06 7]||0.1398 (70%)||[0.0339 6.92 0.053 6.88]||0.1331 (66%)|
|[0.02 7 0.09 7]||0.1486 (74%)||[0.0326 6.98 0.083 6.51]||0.1738 (87%)|
|[0.02 7 0.12 7]||0.1725 (86%)||[0.0353 7.01 0.1041 6.1]||0.1823 (91%)|
|[0.02 7 0.15 7]||0.1814 (91%)||[0.0355 7.04 0.1286 6.51]||0.1934 (97%)|
Phantom experiments were performed with various upper layer thicknesses and two diameter spherical targets of high absorption (phantom 1) and low absorption (phantom 2) coefficients located at various depths. The calibrated absorption and reduced scattering coefficients obtained from the same phantom materials were and for phantom 1 and and for phantom 2 at . The upper layer was intralipid solution of measured optical properties and , and the bottom layer comprised silicon phantoms of measured optical properties and (silicon phantom 1) and and (silicon phantom 2).
The reconstruction results were obtained by using fitted two-layer optical properties from the reflectance measurements. The results of six different experiments of phantoms 1 and 2 are tabulated in Tables 2, 3 . The reconstruction accuracy calculated from the calibrated target values using the semi-infinite and two-layer models are also given in the tables. For the low-contrast phantom 2, because the two-layer model provides more accurate background value, the reconstructed absorption coefficient is greater than 80% more accurate than that obtained from the semi-infinite model at a shallow depth, such as . For the high-contrast phantom 1, the improvement of 10 to 20% is obtained when the two-layer analytic solution is used instead of the semi-infinite model. As expected, the difference in reconstructed values between the two models reduces as the layer depth increases. Because second-layer optical properties mainly affect distance measurements acquired at longer source-detector pairs and these measurements have lower SNRs compared with shorter pairs, the fitted optical properties of the second layer are less accurate than those of the first layer. However, the reconstructed target contrast is not very sensitive to the second layer properties.
Reconstruction results for phantom 1 ( μa=0.23cm−1 and μs′=5.45cm−1 ) using the semi-infinite model and the two-layer model with fitted background values.
|Target Center Depth (C) / Second Layer Depth (D)||Semi-Infinite Model||Two-Layer Model|
|Fitted Background Value [μaμs′] (cm−1)||Reconstructed Maximum μa (cm−1)||Fitted Background Value [μa1μs1′μa2μs2′] (cm−1)||Reconstructed Maximum μa (cm−1)|
|Phantom 1and siliconphantom 1||,||[0.048 5.33]||0.158 (68%)||[0.036 7.89 0.11 3.25]||0.201 (88%)|
|,||[0.039 4.76]||0.169 (73%)||[0.037 5.92 0.079 5.17]||0.206 (90%)|
|,||[0.037 5.40]||0.191 (83%)||[0.027 9.90 0.041 18.53]||0.198 (86%)|
|Phantom 1and siliconphantom 2||,||[0.054 4.79]||0.165 (72%)||[0.020 5.62 0.083 20.21]||0.205 (89%)|
|,||[0.037 5.42]||0.160 (70%)||[0.027 5.53 0.051 5.59]||0.173 (75%)|
|,||[0.033 5.57]||0.186 (81%)||[0.018 7.71 0.051 8.07]||0.187 (81%)|
Reconstruction results for phantom 2 ( μa=0.07cm−1 and μs′=5.50cm−1 ) using the semi-infinite model and the two-layer model with fitted background values.
|Target Center Depth (C) / Second Layer Depth (D)||Semi-Infinite Model||Two-Layer Model|
|Fitted Background Value [μaμs′] (cm−1)||Reconstructed Maximum μa (cm−1)||Fitted Background Value [μa1μs1′μa2μs2′] (cm−1)||Reconstructed Maximum μa (cm−1)|
|Phantom 2and siliconphantom 1||,||[0.065 4.82]||0.129 (184%)||[0.014 5.48 0.16 9.56]||0.0710 (101%)|
|,||[0.045 5.02]||0.104 (148%)||[0.019 7.87 0.099 41.74]||0.0889 (127%)|
|,||[0.042 4.88]||0.105 (150%)||[0.023 9.59 0.048 9.01]||0.0896 (128%)|
|Phantom 2 and silicon phantom 2||,||[0.054 5.04]||0.121 (173%)||[0.021 5.26 0.057 15.08]||0.0868 (124%)|
|,||[0.035 5.11]||0.0944 (135%)||[0.023 5.83 0.053 5.557]||0.0899 (128%)|
|,||[0.031 5.33]||0.103 (147%)||[0.031 5.93 0.058 6.26]||0.1000 (143%)|
An example of imaging results of phantom 1 using silicon phantom 1 as the second layer is given in Figs. 6 and 7 . Figure 6 shows the fitted background optical properties of , , , and versus iterations when the first layer depth is (third line in Table 1 two-layer model). The fitting procedure converges and the fitted background values at iteration of 272 are , , , and .
Figure 7 shows the absorption maps of high-contrast phantom 1 placed in Intralipid layers of 0.9-, 1.5-, and depth. The second layer (silicon phantom 1) was located at 1.5-, 2.1-, and depths, respectively. Figure 7a, 7c, 7e show reconstructed absorption maps obtained at using a semi-infinite model for the layer structure. Reconstructed target maximum absorption coefficients were 0.158 (68%), 0.169 (73%), and 0.191 (83%) , respectively. Figures 7b, 7d, 7f show reconstructed absorption maps obtained at the same wavelength using the two-layer model. The reconstructed target maximum absorption coefficients were 0.201 (88%), 0.209 (91%), and (86%), respectively. About 3 to 20% improvement in reconstruction accuracy was obtained in this set of experiments. The reconstruction results obtained at are very similar to the results at . Therefore, images of are not shown here.
The two-layer-model-based reconstruction was demonstrated using clinical data. Two examples of benign and malignant cases are given here to demonstrate the potential of this new reconstruction technique. Figure 8 shows the ultrasound image from the normal contralateral site [Fig. 8a] and lesion site [Fig. 8b] of a -old patient. The lesion was approximately deep with approximately radius in the direction (depth) and the lateral radius or shape of the lesion was approximately , as seen in the ultrasound image. We can see from Fig. 8a that the upper layer thickness is about with an almost flat interface between the chest wall and the breast tissue. The fitted optical properties of the breast and chest-wall layer from the contralateral site with no lesion were , , , and at , and , , , and at . The fitted two-layer optical properties at converge after 226 iterations, as shown in Figs. 9a, 9b, 9c, 9d . The fitted optical properties by using the semi-infinite method directly were , and at and , and at . Figures 10a, 10b, 10c show the absorption maps at 780 and and the computed total hemoglobin map, respectively, obtained by using a semi-infinite model. Figures 10d, 10e, 10f show the corresponding absorption and total hemoglobin maps obtained using the two-layer model. Both images show low absorption and therefore lower total hemoglobin concentration with layer-based reconstruction providing slightly lower values. Ultrasound-guided core biopsy confirmed that the lesion was a benign fibroadenoma.
Figure 11 shows the ultrasound images obtained from the contralateral and lesion sites of a -old patient. It can be seen from Fig. 11a that the lesion is about deep and the breast tissue and chest-wall layer interface is almost flat with a breast layer thickness of about . The fitted background optical properties of the breast and chest-wall tissue obtained from the contralateral normal breast were , , , and at , and , , , and at . The fitted optical properties by using the semi-infinite method directly were , and at and and at . The absorption maps and the computed total hemoglobin map obtained from with the semi-infinite approximation are shown in Figs. 12a, 12b, 12c and those from the two-layer model are shown in Figs. 12d, 12e, 12f. Both images show high absorption and therefore high total hemoglobin concentration with layer-based reconstruction providing higher values. Ultrasound-guided core biopsy confirmed that the lesion was an early stage invasive carcinoma with ductal and lobular features. These examples demonstrate that the layer-based imaging reconstruction has a great potential to improve the contrast between malignant and benign lesions and therefore provide more accurate diagnosis.
A new optical tomography method based on a two-layer-model analytical solution was developed and compared with the semi-infinite model using simulated data and phantom experiments. Clinical examples were given to demonstrate the potential of this method. It was found that the reconstructed target optical properties were more accurate when the second layer was present in less than depth and its effect was taken into account. This improvement is mainly attributed to the more accurate estimation of background optical properties and more accurate estimation of weight matrix for imaging reconstruction by considering the light propagation effect in the second layer. Our initial results indicate the potential of the layer-based imaging reconstruction method to improve discrimination between malignant and benign lesions.
The fitted first-layer optical properties were more accurate than those of the second layer. This is due to the lower SNR of measurements acquired at longer source-detector pairs. Simulations and phantom experiments showed that the reconstructed target absorption coefficient was not very sensitive to fitting errors of second-layer optical properties. Therefore, the fitting algorithm and the reported two-layer imaging reconstruction are more robust for use on clinical data where a range of optical properties and chest-wall depths may be present. The chest wall may potentially distort the distance measurements.
Optical properties of breast tissue and chest wall vary from patient to patient. From literature data and our own experience, the absorption coefficient is in the range of , and the reduced scattering coefficient is in the range of . The optical properties of the chest wall are not available in the literature. The only available data are optical properties of muscles and bones. In Refs. 22, 23, 24, the muscle was reported to have a large absorption and medium scattering coefficients with in the range of and in the range of . The bone is a highly absorbing and scattering medium. From the literature 25, 26, the bone was reported to have large absorption and scattering coefficients with in the range of and in the range of . We choose optical properties within this range for simulations and phantom experiments.
The two-layer analytical-solution-based imaging reconstruction is limited to a flat interface between first and second layers because the analytic solution is derived in the cylindrical coordinates with depth as a parameter. The breast tissue and chest-wall interface may not be flat in many clinical cases. These factors may introduce errors in estimated two-layer background optical properties, and therefore reconstructed lesion optical properties. More clinical cases will be analyzed to evaluate the robustness of this layer-based tomography method. Finite element methods (FEMs) have been used by many research groups for optical tomography applications27, 28 and are suitable for complex boundary conditions between the breast tissue and chest-wall layers. However, an FEM is computationally intensive, while the proposed analytical-solution-based reconstruction is computationally efficient.
In summary, a new tomographic method was introduced that has been shown to improve light quantification accuracy of lesions imbedded in a two-layer tissue structure. The algorithm was demonstrated by using ultrasound prior to it; however, it can be used alone for optical imaging reconstruction. The robustness of this algorithm will be validated in a larger patient cohort of different breast tissue and chest wall optical properties.
The authors are thankful for the funding support of this work from the National Institutes of Health (R01EB002136), The Patrick & Catherine Weldon Donaghue Medical Research Foundation, and the Army Medical Research and Materiel Command Predoctoral Training award to Chen Xu (81XWH-05-1-0299). The authors thank Dr. John Gamelin for reviewing the manuscript.