|
|
1.IntroductionSince Jobsis showed in 1977 that the relatively good transparency of biological materials to the near-infrared part of the spectrum can be used to help monitor hemoglobin and oxygen saturation information,1 the field of noninvasive optical imaging has experienced sustained growth and development toward clinically viable imaging systems. Studies have been performed using this modality to interrogate different tissues such as brain,2 3 4 breast,5 6 7 muscle,8 forearm,9 and skin,10 or for monitoring light-sensitive treatments such as photodynamic therapy.11 12 In the case of breast tissue, the technique provides new information that may be used to diagnose tumors based on their metabolic and functional status represented through vascularization, oxygenation, and water content. This information is not available from conventional imaging modalities such as mammography, and there is potential for combining these with the added information provided by near infrared (NIR) optical tomography to increase the specificity and sensitivity of clinical screening.13 In optical imaging, the dominant contrast results from increased hemoglobin at the location of the tumor due to neovascularization,14 which causes increased optical absorption,6 15 16 although scattering mechanisms also occur that may hold diagnostic promise as well. Tissue hypoxia typically found in malignant tumors14 is recorded in the oxygen saturation images of the breast, and water content images are also anticipated to provide additional information in cases such as fiberadenoma15 or fibrocystic disease. Recent studies have suggested that structural information related to density of the breast and risk of cancer may also be available through NIR optical tomography.17 18 19 20 Obtaining hemoglobin, oxygen saturation, and water images of the breast using diffuse optical tomography is a three stage process: 1. obtaining measurements of light reflectance from the breast; 2. applying a model of light propagation to recover the bulk optical properties of the breast such as absorption and reduced scattering coefficients; and 3. estimating concentrations of the underlying molecular chromophores in the tissue using their known spectral signatures. Hence, the quantitative reliability of chromophore concentration estimates from the last stage depends significantly on the results from stages 1 and 2. Several experimental techniques involving time domain,21 22 frequency domain,23 24 and continuous wave25 26 signals have been developed for quantification of tissue optical properties. Advancements in instrumentation are improving the quality of the measured data, but the process of obtaining the optical properties from the data is itself a complex problem. Several models have been developed, and of particular interest is the diffusion approximation to the radiation transport equation. This approximation is widely used because it can be computationally fast and robust when implemented using the finite element approach, providing a flexible way to model arbitrary tissue volumes.27 Image reconstruction results have shown accurate spatial recovery of heterogeneities, however, the common consensus seems to be that the spatial resolution is poor, which results physically from the dominance of scattering within tissue.28 29 There is also some limitation on the ability to uniquely separate absorption and scattering coefficients when both variations exist within tissues, although empirical methods have been devised to minimize this problem.30 Even when absorption and scattering images can be recovered (as in most cases), the quantitative characterization of heterogeneities in terms of their absorption coefficient has been found to decrease as the size of the heterogeneity decreases,31 which is an unfortunate side effect of having images that are essentially low-pass filtered in spatial frequency. Studies have shown that optical properties are recovered to within 15 of true value in cases of anomalies of size 25 mm,32 and 25 for heterogeneities of size 17 mm;30 however, for an object of 10 mm, the error limit can be much higher. In 3-D modeling of diffusion theory, Dehghani et al.; have shown that only 15 of the true value could be reached for an object of diameter 8 mm.33 Since this inaccuracy affects the resultant chromophore concentrations, the ability to quantify hemoglobin, oxygen saturation, and water suffers, and poses a limit on the reliability and diagnostic capability of optical tomography. Several studies have focused on developing innovative approaches to solving the problem of maximizing the accuracy of small objects. For example, Paulsen and Jiang implemented a technique to modify the objective function to minimize both the least-square error and the total variation of the field, which provided reconstruction of objects with sharp edges, thereby minimizing the total variation between pixels within the field.34 Other alternative objective functions are possible, and several have been attempted with limited success.35 36 Eppstein et al.; introduced the idea of using successive zonation approaches in image reconstruction, similar to the multigrid methods used in larger scale relaxation methods.37 Another concept along these lines is to use a priori information about the size and position of the tumor to enhance the quantitative accuracy of absorption and scattering of tumors, as demonstrated by Schweiger and Arridge,38 Pogue and Paulsen,39 Dehghani et al.;, and Brooksby et al.;40 This kind of high accuracy characterization would be especially useful to have in small objects. However, the a priori information is typically obtained from MRI and requires a combined MRI-NIR system, which is currently unavailable in most cases. While all of the techniques noted have their strengths and weaknesses, one overriding issue is the need for a problem-dependent calibration factor in the objective function, or the inclusion of a priori information in the process. This work is directed toward developing an algorithm for increased accuracy in quantification of heterogeneities without a priori information and without a strong need for empirically derived calibration factors. It applies a variant of the zonation technique proposed by Eppstein et al.;,37 and follows a systematic investigation of how to implement the approach. Specifically, the existing reconstruction algorithm has been modified to first start with the standard Newton minimization approach, and follow this up with zonation methods, which reduce the problem to one that is uniquely defined and does not require regularization. The modification examined is a three-step reconstruction, where step 1 consists of the existing algorithm that gives the approximate size and position of the anomaly, while step 2 repeats this reconstruction with a two-parameter regularization scheme where the anomaly (detected as a region of interest from step 1) is given a lower value compared to the background, yielding more freedom to update in areas with the largest change in optical properties. This type of reconstruction allows the heterogeneity to recover its true size, and the knowledge of size and position of the anomaly are extracted from this second step. Step 3 is a region-based reconstruction that uses the second step results as a priori information, to produce quantitatively accurate values for the optical parameters. This method substantially improves the image accuracy and has been validated using simulated and experimental data for anomalies between the sizes of 10 and 20 mm with differing contrasts and positions. Finally, the algorithm is applied to clinical subject data to examine how it changes the observed contrast in a tumor relative to the background breast properties. 2.MethodsThe image reconstruction is based on the frequency domain diffusion equation used to model light propagation in highly scattering media,21 41 which is given by where Φ(r,ω) is the isotropic fluence, D(r) is the diffusion coefficient, μa(r) is the absorption coefficient, c is the speed of light in the medium, ω is the modulation frequency, and q0(r,ω) is an isotropic source. The diffusion coefficient can be written as where μs ′ is the reduced scattering coefficient. Equation (1) has been modeled using finite element theory and the boundary condition applied is a type 3 Robin-type boundary condition relating Φ(r,ω) to the internal reflection at the boundary. Complete details of the implementation of this algorithm can be found elsewhere.27 28 42 Briefly, the core of the reconstruction scheme is a Newton-Raphson minimization method for iteratively updating the optical property parameters based on minimization of the standard sum of squared differences between the measured and calculated optical radiance at specific detector locations. This type of error is a measure of the fit of the model to the measured data, often referred to as the projection error, and is commonly used as a surrogate measure of the convergence of the algorithm.2.1.Step 1: Levenberg-Marquardt Reconstruction AlgorithmUsing a truncated Taylor series expansion for Φ to relate it to the update of the optical properties at each iteration gives42 where I is the Jacobian matrix consisting of derivatives of Φ with respect to μa and D for each source-detector location pair, Δμ is the vector containing the perturbations of μa and D, and Φ0 and Φc are the measured and calculated fluence detected at the boundary. Detailed derivations are available in Refs. 29 and 42. Equation (3) is a matrix equation of type Ax=b, where x is the solution vector to be found. The matrix A in this equation is generally ill-conditioned such that the solution x is known to be overwhelmed by data and rounding errors. Many regularization techniques have been studied for the purpose of obtaining stable solutions from this class of ill-posed problems.43 Regularization helps to keep the contributions from different errors within limits by constraining the problem with a priori estimates that bound the solution.Multiplying Eq. (3) by to make the matrix to be inverted square in the number of parameter estimates, we get which is regularized to yield where I is the identity matrix and λ is the regularizing factor that stabilizes the solution (i.e., matrix inversion). The Levenberg-Marquart technique is a widely used regularization technique (introduced by Ref. 44) that normalizes Eq. (5) and makes the Hessian diagonally dominant by using a diagonal matrix G, whose terms are comprised of inverse square roots of the corresponding diagonal terms of to produce45 The sequencing approach for λ as introduced by Marquardt46 47 depends on the projection error and can be written as where c is a constant ( in this work) and χ2 is the projection error. The value of λ in Eq. (6) obviously determines the nature of the relationship. If λ is too high, it dominates Eq. (5) and successive iterations do not yield much change in the update vector. However, if this damping factor is too small, then the problem will be dominated by noise in the data. In our studies, we have found empirically that a starting value of λ in the range of 1000 to 1 reasonably recovers the optical images,48 and noise begins to dominate for starting regularization values below 1. The stopping criterion is met when the projection error is within 2 of the previous iteration’s error. This typically occurs at 10 to 12 iterations, and all the reconstructions in step 1 in this study have begun with λ=10. We use a dual mesh scheme,34 where the Jacobian is calculated on a fine mesh (1785 nodes) and interpolated onto a coarse mesh (425 nodes) to minimize the number of unknowns being estimated in Eq. (6). This algorithm, when implemented, recovers spatial heterogeneities in approximately the correct position.30 Because of the ill-posed nature of the problem, the size of the anomaly is typically blurred with an approximately Gaussian-shaped profile,49 and the average quantification of optical properties is underestimated, with this underestimation depending on the heterogeneity size.30 49 However, since this reconstruction recovers the heterogeneity in its approximate position, this is used as step 1 of the three-step procedure.2.2.Step 2: Two Region RegularizationSince the regularization parameter governs the freedom of the pixels to change more from one iteration to the next, if a lower regularization parameter is specified in a region of interest (ROI) compared to the background, the pixels in the ROI have more opportunity to update so that the properties ultimately reach the true value, and the ROI reaches the true size. In step 2 of our reconstruction, λ in Eq. (6) is rewritten as where λ1 and λ2 are two selected regularization parameters that follow the relationship expressed in Eq. (7) with λ1<λ2.The ROI to be given lower regularization is determined by full width half maximum (FWHM) criterion on the final absorption or scattering coefficient images (when the algorithm has iterated to stopping criterion) from step 1. ROI computed this way is the region containing all the reconstructed nodes having optical parameters greater than or equal to a threshold, defined as the difference between maximum and half of the difference between the maximum and mean of the parameter throughout the image. It is assumed that FWHM gives a reasonable measure of the size of the recovered heterogeneity. Due to the blurring effect, the ROI defined from step 1 typically encompasses the actual heterogeneity as well as a transition zone. In comparison, when the final images from the second step are regionized in the same way using FWHM, a more accurate size estimate of the heterogeneity is obtained. The quantitative accuracy of the property value recovery from step 2 may still deviate from the true value, but it is expected to be closer to the true values relative to step 1. The best numerical values to use for λ1 and λ2 will depend on the object size as well as its contrast and position. In our algorithm, we have examined a series of values for λ1 and λ2 and chosen the optimal pair based on the projection error behavior—the pair giving the lowest projection error at the last iteration (when error change is less that 2) is used. 2.3.Step 3: Region-Based Reconstruction Using Spatial Information from Step 2Step 3 of the process is a region-based reconstruction that has been applied previously in the context of using a priori information.33 38 39 The implementation of this algorithm has been described in detail in Ref. 33. In short, it consists of homogenizing regions in the mesh, by updating selected zones uniformly so that the final image consists of a significantly reduced set of parameter values that represent the property estimates for the individual regions. This type of reconstruction, where the number of unknowns has been reduced to the number of regions, has shown quantitative accuracy of 97 and 93.4 in absorption and scatter for an 8-mm object as shown by Ref. 33 using 3-D reconstruction. Our simulations in 2-D have indicated similar high accuracy, given precise information about the size and position of the heterogeneity, as reported in the results section. In the third step, the geometric constraints are derived from step 2. The mesh is assigned regions through the zonation of areas based on FWHM peaks in the absorption or scatter images from step 2. This algorithm has proved to be sensitive to accurate position and size estimates for the anomaly, and hence, the approximate ROI size obtained from step 1 has proved to be inadequate for accurate quantification. Step 3 can be given an initial guess obtained by averaging the optical coefficients in the step 2 ROI and a regularization factor of 100. It can also be given the same initial guess as the two previous steps, and both converge to the same result with nearly identical projection error, the only difference in the final result being the number of iterations in which it converges. 3.ResultsThe results presented in the following sections focus on the need for quantitative accuracy in NIR tomography (Sec. 1) and demonstrate the application of the three-step reconstruction to simulated and experimental data and the sensitivity of the algorithm to both contrast and position of the anomaly (Sec. 2). Section 3 shows promising results when applied to clinical patient data with an abnormality. While these results are specific to breast imaging, the concept can be extended to all 2-D or 3-D FEM applications of the diffusion approximation for improved quantification in an ROI. 3.1.Section 1: Advantage of A Priori InformationThe boundary data at modulation frequency 100 MHz were generated using the forward solver42 for a circular disk of diameter 86 mm (mesh contains 1785 nodes) with a single absorbing heterogeneity. The model had a background μa=0.005 mm −1 and μs ′=1 mm −1; and data for Secs. 1 and 2 were generated for differing absorption coefficient in the anomaly ranging from 1.2 to 4.0 times the background in increments of 0.2. This type of data was generated for a circular object of diameter 10, 15, and 20 mm. Random Gaussian-distributed noise (1) was added, as in all simulated data used in this study (1 noise has been established to be the level of shot noise in the experimental system31). The object was positioned 10 mm from the edge in all cases. The initial guess for the reconstructions from all three steps was obtained from calibration of the data50 with a homogenous set of measurements generated on the same mesh with the same 1 noise. As noted, the initial guess for steps 2 and 3 can also be obtained from the average reconstructed values from the previous stage of the reconstruction. There is no difference in the two approaches, except in the number of iterations required for the solution to converge. Figure 1(a) shows the average value for absorption coefficient recovered in the region of interest for different diameters of the anomaly (10, 15, and 20 mm), for varying contrast, based on step 1. The reconstruction was stopped when the projection error was less than 2 of the previous iteration, and a starting regularization parameter of 10 was used in all cases. The average value was computed as the mean of the absorption coefficient within the FWHM of the peak absorption coefficient. The high frequency noise near the boundary typically found in our images was removed from consideration during the regionization. Figure 1(a) shows data that indicate quantitative accuracy suffers as the diameter of the anomaly decreases. Specifically, the mean percent error with respect to the true value used to generate the data has increased from 28 for a 20-mm object, to 38 for a 15-mm object, to 47.3 for a 10-mm object. A simple analysis has been pursued to evaluate how this percent error in absorption coefficient estimates propagates into the estimation of chromophore concentrations. The analysis begins with a known set of chromophore concentrations (Hb T , S t O 2, and water) and obtains the absorption coefficients corresponding to these concentrations using matrix multiplication [(μa)=[ε](C) where ε is the molar absorption spectra at the six wavelengths available in our system, as evaluated in Ref. 51, and C is the concentration of hemoglobin, oxy-hemoglobin, and water]. Random noise was added to these sets of absorption coefficients at the six wavelengths, and the concentrations were recovered using a constrained linear least-squares fit. This analysis was carried out for a starting concentration of total hemoglobin=30 μM, oxygen saturation=60, and water=60 (typical concentration found in the breast tissue5), and the percent error estimate used is the mean of 1000 such repetitions. The results show that, using the current six wavelengths used by our system, an error of nearly 50 in the reconstructed absorption coefficient for an anomaly of 10 mm diam can result in as much as 70 error in our estimation of water content in the breast. Since the accuracy in these parameters is likely a key determinant of the diagnostic utility of NIR information, this analysis shows the importance of quantitative accuracy in recovering the optical absorption coefficient, especially in small objects. Figure 1(a) Average reconstructed absorption coefficients plotted as a function of contrast for heterogeneities of different diameters compared to the true values, which illustrate the bias error as the size of the heterogeneity decreases. (b) Same as (a) for 10- and 20-mm heterogeneities using region-based reconstruction with a priori information on the position and size of the inclusion.  If the size and position of the anomaly is known a priori without steps 1 and 2, this information can directly be incorporated in our region-based reconstruction (step 3) to recover the average optical properties in the ROI and background. Results from such a reconstruction for an anomaly of size 10 and 20 mm for varying contrasts is shown in Fig. 1(b), where data were generated in a similar fashion as before, with 1 noise. Accurate results with less than 6 error are obtained for objects of both 10 mm as well as 20 mm diam (5.35 for 10 mm and 2 for 20 mm). Hence, knowledge of anomaly position and size is extremely useful and important in cases of objects of 10 mm in size. The next part of the study is targeted at extracting this information from the existing reconstruction. 3.2.Section 2: Application of Three Stage ReconstructionFigure 2 shows the application of the three-stage reconstruction to the data generated in Sec. 1, where the quantification of the average property estimate has been improved by the final stage of the reconstruction. As noted earlier, the regularization parameter for step 1 is 10, whereas for step 2, the optimal pair of the starting parameters has been chosen by searching over all sets of possible parameters, with λ1 in the range 10 to 25 (varied in steps of 5), and λ2 in the range 1 to 10 (in steps of 1). This range of parameters was chosen empirically after testing the algorithm for different ratios of λ. The pair giving the lowest projection error at the last iteration has been chosen for step 2, and the regularization pair thus obtained varies with contrast and size of the anomaly. The results from this stage are carried on to step 3. The computation time in a multiprocessing environment, for reconstructions on a fine mesh (∼2000 nodes), is approximately 3 min for convergence (∼15 iterations) for first and last step, and varies for the second step, depending on the number of permutations of regularization parameters. There is improvement in accuracy from a 47.3 mean error in step 1 to a 27.3 mean error in the final step of this three-stage reconstruction. There are certain spurious oscillations in the curve for step 3, and these oscillations probably stem from very small scale fluctuations present in the results from step 1, since step 1 provides the starting estimate for the ROI. There are certain artifacts in the image arising from the small size of the anomaly and the underdetermined nature of the reconstruction problem, which are also a possible reason for the oscillations. The reduced scattering coefficient is found to stay constant in the ROI, with a standard deviation of 0.06 mm−1, for varying contrast in absorption coefficient, and there is a 12 increase in the mean of the reduced scattering coefficient resulting from step 3. Figure 2(a) The average reconstructed absorption coefficients in the heterogeneity are shown for the 10-mm-diam anomaly, obtained using the original reconstruction (stage 1) and using a final step of the three stage reconstruction. These are compared against the true values for differing contrast levels. (b) Same as (a) but using a 15-mm-diam anomaly, and (c) same as (a) but using a 20-mm-diam anomaly. (d) Average reconstructed absorption coefficient in the background as a function of contrast, with a 10-mm-diam anomaly from the final step of three-stage reconstruction, which illustrates that the background is not strongly affected by the three stage algorithm (true value=0.005 mm−1).  The same analysis was carried out for data from an anomaly of 15 mm diam, where the results obtained from the final stage of the three-step reconstruction have been compared to the theoretical values and the original results in Fig. 2(b). As in the case of the 10-mm anomaly, Fig. 2(b) is encouraging and shows increased quantitative accuracy. Specifically, the error has dropped to 13 in the final reconstruction step compared to 38 in the original algorithm. Figure 2(c) shows analogous results for a 20-mm-diam anomaly. As expected, the algorithm achieves better quantification of larger objects (mean error of 5.5 in the final step, compared to 28.2 in the original reconstruction). The background absorption is found to stay constant for all three anomaly sizes, and results from a 10-mm case are shown in Fig. 2(d). To illustrate the importance of the sequence of steps in the three-stage reconstruction, a comparison of the results obtained by the three stage reconstruction, using just steps 1 and 3, is shown in Fig. 3(a) for the 15-mm anomaly case. The mean percent error improvement by using the three-stage reconstruction is 13 compared to 24 observed by using just steps 1 and 3. The plot of projection error as a function of iterations for different reconstruction steps (15-mm anomaly, contrast=2.2 case) is shown in Fig. 3(b). The projection error in step 3, which generates the most accurate average optical parameters, is higher than steps 1 and 2, which is expected since this is a homogeneous estimate in the anomaly and background, and a two-unknown problem is not sufficient to compensate for high frequency variations near the boundary. This observation is consistent in all the studies done, including the experimental and patient data. Figure 3(c) shows the difference in the projection error between starting the third step with same initial guess as the previous steps, relative to using an average of parameters from the previous steps of reconstruction. Both converge to the same projection error in step 3. The quantitative values of the recovered optical properties were also found to be almost identical. Figure 3(a) Comparison of average reconstructed absorption coefficient in the heterogeneity for the 15-mm anomaly, obtained using just step 1, using steps 1 and 3, without step 2, and using the sequence of steps 1, 2, and 3, which shows the advantage of the sequence of the three-stage algorithm. (b) Projection error versus iteration number for the three steps of the reconstruction algorithm when recovering an image from simulated data for a 15-mm-diam anomaly with a contrast of 2.2 with the background. (c) Projection error versus iteration number for the final step of reconstruction, using an initial guess from calibration compared to an initial guess from the previous step of reconstruction.  The results shown in Fig. 2 indicate that the three-stage reconstruction yields improved quantification for all degrees of contrast studies; however, all of the results have been based on a single location of the anomaly. Application of the algorithm for different anomaly positions at a fixed contrast (2:1, ROI relative to background) has been explored. Data were generated as before for an anomaly with the same background properties as for Figs. 1 and 2 (μa=0.005 mm −1 and μs ′=1 mm −1) with 1 random noise. This time the position of the anomaly was varied starting 10 mm from the edge to center in steps of 3 mm for both 10 and 20 mm diam. After calibration, the three-step reconstruction was carried out and Fig. 4(a) shows the reconstructed average absorption coefficient for step 1 (original reconstruction) and step 3 (final stage), along with the theoretical value. As seen from the figure, the final step provides better quantification compared to step 1. Overall, the reconstruction shows variation with the position, but no clear trend is visible. The mean error was decreased from 40 in step 1 to 19.47 in the final step. The variation in the reconstructed values in the final stage for different positions was 9.32 of their mean, suggesting that the algorithm is largely insensitive to position of the anomaly. To show the improvement in quantification and insensitivity to position using experimental data, data with similar parameters were collected. Specifically, an 84-mm cup phantom was made of epoxy resin, with India ink for absorption and titanium dioxide for scattering,31 (optical properties measured by the system) with a wall thickness of 5 mm. This cup was filled with Intralipid solution of the same optical properties as the phantom (μa=0.0045 mm −1 and μs ′=1.12 mm −1, measured independently). A cylindrical rod of diameter 10 mm (made similar to the cup phantom), with optical properties of μa=0.009 mm −1 and μs ′=1.5 mm −1, was moved starting 10 mm from the edge to the center in 4-mm increments, where measurements were acquired using a frequency domain NIR system.52 The three-step reconstruction was carried out on the data and the results are shown in Fig. 4(b). The graph shows that there is improvement in accuracy using the three-step reconstruction, and again similarly to Fig. 4(a), the variation in the absorption coefficient offers no clear trend. The error has been reduced from 27.56 in step 1 to 16.26 in step 3. There is an overshoot in the absorption coefficient values in some positions, which may be due to experimental artifacts. Figure 4(a) Average reconstructed absorption coefficient in the inclusion for different positions of a 10-mm region obtained from the original reconstruction (stage 1 only) and the final step of the three-stage algorithm. Results are also compared with true values. The measurement data were simulated and the position of the inclusion was varied from a location 10 mm from edge to the center, in increments of 3 mm. (b) Same as (a), except that experimental data is used from a 10-mm object. The position was varied from a location 10 mm from edge to the center in increments of 4 mm. (c) Same as (a) for a 20-mm heterogeneity with simulated data. (d) Same as (b) for a 20-mm object. The position of the inclusion was varied from 10 mm from the edge to the center in increments of 3 mm.  Figures 4(c) and 4(d) show the analyses for a 20-mm anomaly. The simulations were performed in the same way for different positions of the anomaly, and the trends are similar to Fig. 4(a): the simulated data show improvement in accuracy from 23.2 in the original reconstruction to 9.8 in the final stage. The same experiments were also conducted for a rod of diameter 20 mm, and quantification has been improved from 19.64 to 6.8 in this case. 3.3.Section 3: Clinical Abnormality—Absorption and Scattering Heterogeneity in Different PositionsAll the results presented in the previous sections had either an absorption-only heterogeneity or both absorption and scattering heterogeneity located in the same position. However, in the clinical environment, breast images with abnormalities may contain absorption and scattering heterogeneities in different positions. The algorithm can still be applied in this case, where the regionization from step 1 is based on both absorption and scattering images, and these regions are assigned lower regularization compared to the background for step 2. This ensures that both absorption and scattering heterogeneities are underdamped, and hence can update to the true size and quantitative property value. The regionization from step 2 is again absorption and scatter based, and the mesh is divided into different zones so that the region-based reconstruction for step 3 can generate accurate average optical parameters in the assigned areas. This algorithm was implemented and explored with simulated data using background properties μa=0.01 mm −1 and μs ′=1.0 mm −1, and a 15-mm absorption anomaly having twice the background (μa=0.02 mm −1) located with center coordinates at (0,−20) and a 15-mm scattering object of μs ′=1.5 mm −1 with center coordinates close to (0,20) and 1 random Gaussian noise. It was found that the scattering object always appears in reduced size compared to its true extent, while the absorption object is exaggerated. Empirically, it was found, increasing the scattering region obtained by FWHM of the scattering image by 20 of its radius and decreasing the absorption region by 20 of its radius in stages works with simulated data to bring the regionized size to its original expanse. The recovered size as it emerges in step 3 was found to be almost identical to the original 15-mm size of the anomalies, and there was no observable shift in position. The images from the three steps are shown in Figs. 5(a), 5(b), and 5(c). An accurate quantification with 13 error in μa and 5.8 error in μs ′ was found in the final step, and the values for the background and the anomaly for the initial and final steps are reported in Table 1. Figure 5(a), (b), and (c) Reconstructed images from the three stages of reconstruction, respectively, using simulated measurement data with 1 Gaussian noise generated from a phantom with background properties of: μa=0.01 mm −1 and μs ′=1.0 mm −1, and a 15-mm absorption object near 6 o’clock with μa=0.020 mm −1 and 15-mm scattering object with μs ′=1.5 mm −1 near 12 o’clock.  Table 1
This algorithm was automated and implemented on measurements from the left breast of a Caucasian subject of age 38 years, diagnosed with an intraductal carcinoma from biopsy. Mammograms showed that the cancer was in the 6:30 clock-face position and 3 cm from the nipple. The NIR breast exam was completed within an imaging study that has approval by the Dartmouth Institutional Review Board, where all subjects participate under informed consent. Figure 6 shows the application of the 3-step algorithm on this clinical data measured at 785 nm. From the reconstructed images in step 1, we can see that the increased absorption appears at the position of the tumor and increased scattering contrast presents in a different position. These two regions were given lower regularization, and step 2 was carried out, which yields comparable images to step 1. Step 3 shows the three-region reconstruction where the background, absorption anomaly, and scattering anomaly, whose size and position were obtained from step 2, were updated uniformly. Table 2 shows the average absorption and scattering values from the three steps, and as seen, the absorption coefficient is much higher in step 3 compared to steps 1 and 2, and is likely more accurate relative to previous results. The contrast has been enhanced in the final stage through increased absorption at the site of the tumor and lower background values. Figure 6(a), (b), and (c) Reconstructed images from the three stages of the reconstruction, using clinical data at 785-nm wavelength from the left breast of a subject having an intraductal carcinoma in the lower region of the breast. The measurements were recorded 30 mm from the nipple.  Table 2
4.Discussion and ConclusionsThe primary benefit of NIR diffuse optical tomography arises from the ability to image tissue volumes that characterize the breast. Early iterations of the reconstruction can yield images with the correct heterogeneity locations, but at these early iterations, the quantitative values of optical parameters are significantly underestimated (by nearly 50 or more). Later iterations show more noise and artifacts in the image, but the quantitative values of optical properties recovered are closer to the true values. The focus of the current study is quantitative property accuracy rather than overall image quality, and hence, all reconstructions have been allowed to proceed until the projection error reduction is less than 2. The simple analysis in Sec. 1 has shown that a 50 error in quantifying the absorption coefficient of an object translates into nearly 70 error in estimation of its water content, which indicates the importance of quantification in NIR tomography. The three-stage reconstruction has been implemented on simulated data with 1 random Gaussian noise for a range of typical contrasts in tumors of up to a maximum, absorption of 4 times the background. The results show improved quantification in 10-, 15-, and 20-mm-diam region data by bringing the reconstructed optical property values closer to tolerable error limits. The background absorption was found to stay constant, and hence improved contrast has also been observed using this method. The graphs in Fig. 2 show that the algorithm yields results that are insensitive to contrast compared to the original single-step algorithm where the higher contrasts were more difficult to recover. The algorithm has been tested for variation in position of the anomaly (from edge to center) and insensitive to object position because changes expected in positioning are treated in the first two steps. This trend was also evident in experimental data for anomaly sizes of 10 and 20 mm, and the experimental results show improved quantification as well, with error reductions from 28 to 16 in the 10-mm case and from 20 to 7 in the 20-mm case. The results in Sec. 3 indicate that the algorithm can accommodate both absorption and scattering heterogeneities, where accurate values within error limits of 15 for absorption and 6 for scattering have been obtained in simulated data with 1 random Gaussian noise. The algorithm has been successfully applied to clinical data, and the tabular column in Table 2 shows a much higher contrast of 3.8 in the absorption coefficient between tumor and background in the final step, compared to 1.23 in the original reconstruction. The main drawback with patient data is that there is no knowledge of the true optical properties; however, studies performed in phantoms and with simulated data suggest that the values from the final step of reconstruction are much closer to the actual properties relative to the original reconstruction. Tromberg et al.;16 showed measurements with a frequency domain photon migration (FDPM) system, where a contrast of approximately 3 fold was observed between tumor versus normal sites when the probe was placed just 5 mm lateral of the tumor center. Our results also show that quantification of contrast in the focal region has been increased substantially when the tumor is zoned through region-based reconstruction. The images of patient data in Fig. 6 also demonstrate the tradeoff between image quality and quantitative accuracy. In the three-stage reconstruction, the images from the final step contain homogeneous regions in the breast, and characterize the main anomalies. These may not be the best images, since different kinds of tissue, such as glandular, fatty, and fibrous, are present in the breast. However, the regionization improves the reliability of the modality in quantifying chromophore concentrations in the tumor, at the cost of characterization of other heterogeneities in the breast. Since this characterization is available from step 1, it is not completely lost, and the information from the final step is more reliable in quantifying the optical properties of the tumor. Finally, the method of defining the two regularization parameters in step 2 by the projection error approach is very robust, since this includes the change in these parameters due to changes in contrast and size of the anomalies. The algorithm is automated and can easily be executed in a multiprocessor computing environment. Future studies will aim to incorporate the three-step reconstruction for all different wavelength measurements of phantom solutions and clinical data, so that the chromophore concentrations can be calculated; and oxygenation status of the tumor will be studied in depth based on the more accurate quantification of oxygen saturation and other NIR parameters. AcknowledgmentsThis work has been supported through NIH grants PO1CA80139 and RO1CA69544 and the Norris Cotton Cancer Center through its shared services. The authors would like to acknowledge collaboration from Steven P. Poplack, Sandra Soho, and Christine Kogel in carrying out the clinical examinations. REFERENCES
F. F. Jobsis
,
“Non-invasive, infra-red monitoring of cerebral and myochardial oxygen sufficiency and circulatory parameters,”
Science , 198 1264
–1267
(1977). Google Scholar
J. C. Hebden
,
A. Gibson
,
R. M. Yusof
,
N. Everdell
,
E. M. Hillman
,
D. T. Delpy
,
S. R. Arridge
,
T. Austin
,
J. H. Meek
, and
J. S. Wyatt
,
“Three-dimensional optical tomography of the premature infant brain,”
Phys. Med. Biol. , 47
(23), 4155
–66
(2002). Google Scholar
G. Strangman
,
D. A. Boas
, and
J. P. Sutton
,
“Non-invasive neuroimaging using near-infrared light,”
Biol. Psychiatry , 52 679
–693
(2002). Google Scholar
J. S. Wyatt
,
“Cerebral oxygenation and haemodynamics in the foetus and newborn infant,”
Philos. Trans. R. Soc. London, Ser. B , 352 697
–700
(1997). Google Scholar
T. Durduran
,
R. Choe
,
J. P. Culver
,
L. Zubkov
,
M. J. Holboke
,
J. Giammarco
,
B. Chance
, and
A. G. Yodh
,
“Bulk optical properties of healthy female breast tissue,”
Phys. Med. Biol. , 47 2847
–2861
(2002). Google Scholar
B. W. Pogue
,
S. P. Poplack
,
T. O. McBride
,
W. A. Wells
,
O. K. S.
,
U. L. Osterberg
, and
K. D. Paulsen
,
“Quantitative hemoglobin tomography with diffuse near-infrared spectroscopy: Pilot results in the breast,”
Radiology , 218
(1), 261
–266
(2001). Google Scholar
N. Shah
,
A. Cerussi
,
C. Eker
,
J. Espinoza
,
J. Butler
,
J. Fishkin
,
R. Hornung
, and
B. Tromberg
,
“Noninvasive functional optical spectroscopy of human breast tissue,”
Proc. Natl. Acad. Sci. U.S.A. , 98
(8), 4420
–4425
(2001). Google Scholar
E. Gratton
,
S. Fantini
,
M. A. Franceschini
,
G. Gratton
, and
M. Fabiani
,
“Measurements of scattering and absorption changes in muscle and brain,”
Philos. Trans. R. Soc. London, Ser. B , 352 727
–735
(1997). Google Scholar
P. Van der Zee
,
M. Cope
,
S. R. Arridge
,
M. Essenpreis
,
L. A. Potter
,
A. D. Edwards
,
J. S. Wyatt
,
D. C. McCormick
,
S. C. Roth
,
E. O. Reynolds
et al.;,
“Experimentally measured optical pathlengths for the adult head, calf and forearm and the head of the newborn infant as a function of inter-optode spacing,”
Adv. Expt. Med. Biol. , 316 143
–153
(1992). Google Scholar
I. V. Meglinski
and
S. J. Matcher
,
“Computer simulation of the skin reflectance spectra,”
Comput. Methods Programs Biomed. , 70 179
–186
(2003). Google Scholar
B. C. Wilson
and
M. S. Patterson
,
“The physics of photodynamic therapy,”
Phys. Med. Biol. , 31 327
–360
(1986). Google Scholar
T. Johansson
,
M. S. Thompson
,
M. Stenberg
,
C. A. Klinteberg
,
S. Andersson-Engels
,
S. Svanberg
, and
K. Svanberg
,
“Feasibility study of a system for combined light dosimetry and interstitial photodynamic treatment of massive tumors,”
Appl. Opt. , 41
(7), 1462
–1468
(2002). Google Scholar
J. G. Elmore
,
M. B. Barton
,
V. M. Moceri
,
S. Polk
,
P. J. Arena
, and
S. W. Fletcher
,
“Ten-year risk of false positive screening mammograms and clinical breast examinations,”
N. Engl. J. Med. , 338
(16), 1089
–1096
(1998). Google Scholar
P. Vaupel
,
F. Kallinowski
, and
P. Okuneiff
,
“Blood flow, oxygen and nutrient supply, and metabolic microenvironment of human tumors: A review,”
Cancer Res. , 49 6449
–6465
(1989). Google Scholar
A. Cerussi
,
D. Jakubowski
,
N. Shah
,
F. Bevilacqua
,
R. Lanning
,
A. J. Berger
,
D. Hsiang
,
J. Butler
,
R. F. Holcombe
, and
B. J. Tromberg
,
“Spectroscopy enhances the information content of optical mammography,”
J. Biomed. Opt. , 7
(1), 60
–71
(2002). Google Scholar
B. J. Tromberg
,
N. Shah
,
R. Lanning
,
A. Cerussi
,
J. Espinoza
,
T. Pham
,
L. Svaasand
, and
J. Butler
,
“Non-invasive in vivo characterization of breast tumors using photon migration spectroscopy,”
Neoplasia , 2
(1,2), 26
–40
(2000). Google Scholar
B. W. Pogue
,
S. Jiang
,
H. Dehghani
,
C. Kogel
,
S. Soho
,
S. Srinivasan
,
X. Song
,
S. P. Poplack
, and
K. D. Paulsen
,
“Characterization of hemoglobin, water and NIR scattering in breast tissue: analysis of inter-subject variability and menstrual cycle changes relative to lesions,”
J. Biomed. Opt. , 9
(3), 541
–552
(2004). Google Scholar
S. Srinivasan
,
B. W. Pogue
,
S. Jiang
,
H. Dehghani
,
C. Kogel
,
S. Soho
,
J. J. Gibson
,
T. D. Tosteson
,
S. P. Poplack
, and
K. D. Paulsen
,
“Interpreting hemoglobin and water concentration, oxygen saturation and scattering measured in vivo by near-infrared breast tomography,”
PNAS , 100
(21), 12349
–12354
(2003). Google Scholar
M. S. Patterson
,
C. B.
, and
B. C. Wilson
,
“Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,”
Appl. Opt. , 28 2331
–2336
(1989). Google Scholar
B. Chance
,
S. Nioka
,
J. Kent
,
K. McCully
,
M. Fountain
,
R. Greenfield
, and
G. Holtom
,
“Time-resolved spectroscopy of hemoglobin and myoglobin in resting and ischemic muscle,”
Anal. Biochem. , 174 698
–707
(1988). Google Scholar
B. W. Pogue
,
M. S. Patterson
,
H. Jiang
, and
K. D. Paulsen
,
“Initial assessment of a simple system for frequency domain diffuse optical tomography,”
Phys. Med. Biol. , 40 1709
–1729
(1995). Google Scholar
S. Fantini
,
M. A. Franceschini-Fantini
,
J. S. Maier
,
S. A. Walker
,
B. Barbieri
, and
E. Gratton
,
“Frequency-domain multichannel optical detector for noninvasive tissue spectroscopy and oximetry,”
Opt. Eng. , 34 32
–42
(1995). Google Scholar
E. L. Hull
,
M. G. Nichols
, and
T. H. Foster
,
“Quantitative broadband near-infrared spectroscopy of tissue-simulating phantoms containing erythrocytes,”
Phys. Med. Biol. , 43
(11), 3381
–3404
(1998). Google Scholar
A. Kienle
,
L. Lilge
,
M. S. Patterson
,
R. Hibst
,
R. Steiner
, and
B. C. Wilson
,
“Spatially resolved absolute diffuse reflectance measurements for noninvasive determination of the optical scattering and absorption coefficients of biological tissue,”
Appl. Opt. , 35 2304
–2314
(1996). Google Scholar
S. R. Arridge
,
M. Schweiger
,
M. Hiraoka
, and
D. T. Delpy
,
“A finite element approach for modeling photon transport in tissue,”
Med. Phys. , 20
(2), 299
–309
(1993). Google Scholar
K. D. Paulsen
and
H. Jiang
,
“Spatially varying optical property reconstruction using a finite element diffusion equation approximation,”
Med. Phys. , 22
(6), 691
–701
(1995). Google Scholar
S. R. Arridge
and
M. Schweiger
,
“Image reconstruction in optical tomography,”
Philos. Trans. R. Soc. London, Ser. B , 352 717
–726
(1997). Google Scholar
T. O. McBride
,
B. W. Pogue
,
S. Jiang
,
U. L. Osterberg
,
K. D. Paulsen
, and
S. P. Poplack
,
“Initial studies of in vivo absorbing and scattering heterogeneity in near-infrared tomographic breast imaging,”
Opt. Lett. , 26
(11), 822
–824
(2001). Google Scholar
B. W. Pogue
,
C. Willscher
,
T. O. McBride
,
U. L. Osterberg
, and
K. D. Paulsen
,
“Contrast-detail analysis for detection and characterization with near-infrared diffuse tomography,”
Med. Phys. , 27
(12), 2693
–2700
(2000). Google Scholar
T. O. McBride
,
B. W. Pogue
,
E. Gerety
,
S. Poplack
,
U. L. Osterberg
, and
K. D. Paulsen
,
“Spectroscopic diffuse optical tomography for quantitatively assessing hemoglobin concentration and oxygenation in tissue,”
Appl. Opt. , 38
(25), 5480
–5490
(1999). Google Scholar
H. Dehghani
,
B. W. Pogue
,
J. Shudong
,
B. Brooksby
, and
K. D. Paulsen
,
“Three-dimensional optical-tomography: resolution in small-object imaging,”
Appl. Opt. , 42
(16), 3117
–3128
(2003). Google Scholar
K. D. Paulsen
and
H. Jiang
,
“Enhanced frequency-domain optical image reconstruction in tissues through total-variation minimization,”
Appl. Opt. , 35
(19), 3447
–3458
(1996). Google Scholar
H. W. Engl
and
W. Grever
,
“Using the L-curve for determining optimal regularization parameters,”
Numer. Math. , 69 25
–31
(1994). Google Scholar
A. Neumaier
,
“Solving ill-conditioned and singular linear systems: a tutorial on regularization,”
SIAM Rev. , 40
(3), 636
–666
(1998). Google Scholar
M. J. Eppstein
,
D. E. Dougherty
,
T. L. Troy
, and
E. M. Sevick-Muraca
,
“Biomedical optical tomography using dynamic parameterization and Bayesian conditioning on photon migration measurements,”
Appl. Opt. , 38
(10), 2138
–2150
(1999). Google Scholar
M. Schweiger
and
S. R. Arridge
,
“Optical tomographic reconstruction in a complex head model using apriori region boundary information,”
Phys. Med. Biol. , 44 2703
–2721
(1999). Google Scholar
B. W. Pogue
and
K. D. Paulsen
,
“High resolution near infrared tomographic imaging simulations of rat cranium using apriori MRI structural information,”
Opt. Lett. , 23
(21), 1716
–1718
(1998). Google Scholar
M. S. Patterson
,
J. D. Moulton
,
B. C. Wilson
,
K. W. Berndt
, and
J. R. Lakowicz
,
“Frequency-domain reflectance for the determination of the scattering and absorption properties of tissue,”
Appl. Opt. , 30
(24), 4474
–4476
(1991). Google Scholar
H. Jiang
,
K. D. Paulsen
,
U. L. Osterberg
,
B. W. Pogue
, and
M. S. Patterson
,
“Optical image reconstruction using frequency-domain data: simulations and experiments,”
J. Opt. Soc. Am. A , 13
(2), 253
–266
(1996). Google Scholar
K. Levenberg
,
“A method for the solution of certain nonlinear problems in least squares,”
Q. Appl. Math. , 2 164
–168
(1944). Google Scholar
P. M. Meaney
,
E. Demidenko
,
N. K. Yagnamurthy
,
D. Li
,
M. W. Fanning
, and
K. D. Paulsen
,
“A two-stage microwave image reconstruction procedure for improved internal feature extraction,”
Med. Phys. , 28
(11), 2358
–2369
(2001). Google Scholar
D. W. Marquardt
,
“An algorithm for least squares estimation of nonlinear parameters,”
J. Soc. Ind. Appl. Math. , 11 431
–441
(1963). Google Scholar
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
(13), 2950
–2961
(1999). Google Scholar
X. Song
,
B. W. Pogue
,
T. D. Tosteson
,
T. O. McBride
,
S. Jiang
, and
K. D. Paulsen
,
“Statistical analysis of nonlinearly reconstructed near-infrared tomographic images: Part II—Experimental interpretation,”
IEEE Trans. Med. Imaging , 21
(7), 764
–772
(2002). Google Scholar
S. Srinivasan
,
B. W. Pogue
,
S. Jiang
,
H. Dehghani
, and
K. D. Paulsen
,
“Validation of hemoglobin and water molar absorption spectra in near-infrared diffuse optical tomography,”
Proc. SPIE , 4955 407
–415
(2003). Google Scholar
T. O. McBride
,
B. W. Pogue
,
S. Jiang
,
U. L. Osterberg
, and
K. D. Paulsen
,
“Development and calibration of a parallel modulated near-infrared tomography system for hemoglobin imaging in vivo,”
Rev. Sci. Instrum. , 72
(3), 1817
–1824
(2001). Google Scholar
|
||||||||||||||||||||||||||||||||||||||||||||||||||
