|
1.IntroductionDiffuse optical tomography (DOT) employs near-infrared light to probe the optical properties of highly scattering biological tissues.1–3 DOT has shown promise in breast cancer imaging4–15 because optical methods are sensitive to changes in physiological parameters such as blood volume and tissue oxygenation, alterations of which are biomarkers for cancer and other disease processes. Typical source-detector arrangements in DOT include the ring16 and slab14,17,18 geometries. In the latter case, both transmission14 and reflection19 measurements have been used, though transmission measurements provide more sensitivity for deep tissues since the detected light in this geometry is more likely to travel through all areas of interest. Recently, it has been demonstrated that the image quality in DOT,20–22 and in other imaging modalities such as inverse diffraction,23 can be significantly improved by utilization of large data sets in the reconstruction. The plane-parallel transmission geometry is particularly well suited for utilization of larger data sets, because it offers the possibility of noncontact scanning methods20,24–27 wherein a computer-controlled beam scanner on one side of the sample is used for illumination and a megapixel charge coupled device (CCD) camera is used for detection. A characteristic feature of noncontact scanning is the availability of very large data sets with up to ∼109 independent measurements, e.g., with ∼103 source positions and ∼106 CCD pixels per source. Naturally, utilization of data sets consisting of more than ∼105 independent measurements presents a serious computational challenge. For this reason, we developed fast algorithms capable of reconstructing very large data sets in simple imaging geometries (including the slab).28–32 Numerical simulations29 have demonstrated the full potential of these methods, but the simulations also indicate that large imaging windows are required. To obtain optimum resolution, for example, the dimensions on both sides of the slab where sources and detectors are scanned should be larger by a factor of about 3 in both transverse directions compared to the slab width, but in practice, a somewhat smaller ratio is expected to be sufficient due to the presence of experimental noise and other imaging system imperfections.20,21 Unfortunately, in clinical breast imaging, the large windows described above are not achievable due to physical limitations imposed by the chest wall, and its consequences in DOT are poorly understood. To address this problem, we employ engineered phantoms to study the effects of the chest wall on image reconstruction. These tissue phantoms are similar to many that have been employed with success in the DOT community for preclinical, in vitro investigations to ascertain the utility and limitations of various image reconstruction schemes.7,33,34 One advantage of phantom experiments is that we can compare reconstructions obtained under similar conditions but with different imaging windows, some of which would be physically unavailable in vivo. Based on these experimental results, we discuss and compare two approaches for reconstructing large data sets under the condition that the large imaging windows, which are required by the methods of Refs. 28 through 32, are unavailable; we then explore methods to ameliorate the chest wall effects. The main conclusion of this paper is that both the analytical and algebraic data-intensive linearized image reconstruction methods can produce reasonable results, provided the data points are appropriately restricted to exclude measurements that are strongly influenced by the chest wall. Under these conditions, an absorbing target with subcentimeter features can be clearly reconstructed in the middle of a 6 cm slab, even when the chest wall is only 2 cm from the target. We obtain good images of the target even in the presence of a large chest wall phantom that introduces significant nonlinearities into the inverse problem due to its larger absorption coefficient compared to the background as well as due to its size. Moreover, we discovered a data restriction condition such that the presence of the chest wall phantom imposes minimal artifacts or distortions in the image. The performance of both algebraic and analytic image reconstruction methods were compared under this condition and, while neither method is perfect, we believe that a role for both methods in DOT exists, the choice depending upon the particular clinical application. The reminder of this paper is organized as follows. In Sec. 2, we describe the experimental set up. Section 3 contains details of the image reconstruction methods. Section 4 explains our approach to data restriction. Section 5 presents the results and Sec. 6 contains a brief summary. 2.ExperimentThe experimental apparatus is shown schematically in Fig. 1. Briefly, a collimated continuous-wave (CW) 785 nm laser diode is coupled to a two-dimensional galvanometer scanner (Thorlabs GVS012). The laser beam with a focused spot size of 0.5 mm is raster scanned on a 35×35 square grid with a 4 mm spacing covering a 13.6×13.6 cm2 area on one side of the imaging tank (whose overall dimensions are 44×44×6 cm3). The thickness of the tank was chosen to be close to the average compression used in our previous clinical studies based on diffuse optical tomography.7,14 As each source position is illuminated, data is collected from the opposite side of the tank over a 21.2×21.2 cm2 field of view (FOV) area with a CCD camera (Andor, DV887ECS-UV, lens 25 mm F/0.95). The FOV was mapped to the grid of 512×512 CCD pixels. This corresponds to a rectangular grid on the surface of the tank with the spacing p=0.416 mm. Fig. 1DownloadA bar target is suspended in the mid-plane of the tank (3 cm from either surface) using monofilament fishing line. The target is made of silicon rubber (RTV-12, General Electric), titanium oxide (T-8141, Sigma-Aldrich) and carbon black (Raven 5000 Ultra Powder II), with the absorption coefficient μa=0.2 cm−1 and the reduced scattering coefficient μ′s=7.5 cm−1. The tank is filled with a scattering fluid (μa=0.05 cm−1 and μ′s=7.5 cm−1); these background optical properties are similar to those used in previous in vitro and clinical research. The contrast between the target and the surrounding fluid is purely absorptive with the ratio of about 4. A chest wall phantom (Biomimic, INO μa=0.1 cm−1 and μ′s=5.0 cm−1, dimensions 40×20×5.8 cm3) is suspended at various distances d from the top edge of the bar target (d=2, 5, 8, 11, 14, 17 cm). The optical properties of the chest wall phantom were chosen to mimic muscle tissue.35–37 Thus both absorptive and scattering contrast exists between the chest wall phantom and the background fluid. The bar target and the chest wall phantom are shown in Fig. 2. Note that the chest wall phantom almost entirely fills the imaging tank; the clearance between the chest wall phantom and the inner surfaces of the tank is 1 mm on both sides. Fig. 2Download3.Image Reconstruction3.1.Linearized Integral EquationsThis research employs linear image reconstruction methods and CW data. In principle, one could also resort to time-38–43 or frequency-resolved44–47 measurements and nonlinear reconstruction methods1,31 to obtain a reconstruction of the target and the chest wall phantom simultaneously; however, these approaches require more expensive and complex instrumentation, as well as more time-consuming computational schemes. For our linear approach, two independent measurements of the transmitted intensity are taken, one in a homogeneous (reference) slab and the other in a slab with the target and the chest wall phantom present. We denote these measurements by I0(rd,rs) and I(rd,rs), respectively, where rd and rs are the two-dimensional vectors specifying the lateral positions of the detector and the source on the respective surfaces of the tank. In this work, I0 and I are expressed in CCD counts. These measurements can be related to the medium optical properties through the relations: Here A(rd) and B(rs) are unknown coupling coefficients associated with the source and detection system which can be excluded from consideration as described below. G(rd,rs) and G0(rd,rs) are the Green’s functions for the diffusion equation in the slab with the inhomogeneities present and in the homogeneous (reference) slab, respectively. The latter is known analytically.32 The two Green’s functions are mathematically related by the Dyson equation. Under the assumption that the scattering properties of the medium are spatially uniform, the Dyson equation has the form Here the integral is taken over the volume of the slab and δα(r)=c[μa(r)−μa0] is the deviation of the absorption coefficient at a given point from the background value cμa0, where c is the average speed of light in the medium.The procedure of linearization consists of making an approximation to Eq. (2) so as to exclude the unknown function G(r,rs) from the right-hand side. In the first Rytov approximation,48 Eq. (2) is rewritten as Consequently, we define the data function ϕ(rd,rs) according to and use Eq. (3) to obtain an integral equation of the form Here the right-hand side contains the measurable data function, and the left-hand side is an integral transform of the contrast whose kernel is known analytically. This formulation of the linearized inverse problem of DOT is standard.48 The Green’s function G0(r,r′), which enters Eq. (5), is defined in Ref. 32.3.2.Analytical Reconstruction MethodThis method is described in detail in Ref. 32. Its applications to the slab geometry with experimental data have been reported in Refs. 20 and 21. The method requires the Fourier transform of the data function ϕ(rd,rs) with respect to both variables. To obtain this Fourier transform, the data function must be measured over sufficiently large areas so that the integrals involved (approximated by sums in practical implementations) have converged. This requirement translates into the need for sufficiently large imaging windows that was discussed above. In our experiments, some data points are strongly affected by the presence of the chest wall. The actual source and detector positions for the affected data points depend on the separation d between the chest wall and the top of the bar target. At the two smallest separations (d=2 cm and 5 cm), contamination of the data function by the chest wall is significant. Under these circumstances, two approaches for data analysis are natural to consider:
Note that in the implementation of the analytical method, we have followed closely Ref. 21 with some minor modifications. Thus, as in Ref. 21, we have used the change of variables ψ(q,p)=˜ϕ(q+p,−p), where ˜ϕ(qd,qs) is the Fourier transform of the data function defined in Eq. (4). Here we use a “symmetric version” in which ψ(q,p)=˜ϕ(q/2+p,q/2−p). The variables p and q were sampled as follows: qij=Δ(iˆx+jˆy), pkl=Δ(kˆx+lˆy), where Δ=2π/(331×p), −28≤i,j≤28 and 0≤k, l≤6. This corresponds to using 572=3249 “external” degrees of freedom and 72=49 “internal” degrees of freedom, where the terminology of Ref. 32 has been used, or a total of 49×3249≃1.6×105 independent Fourier-space data points. Using these parameters, the necessary computations require 7×109 floating point operations. On an Intel Core 2 Duo processor this translates into 7 min of computation time. 3.3.Algebraic Reconstruction MethodAlgebraic image reconstruction methods are obtained by discretization of the volume integral in Eq. (5) and computation of the pseudo-inverse for the obtained system of linearized equations.49 This method does not require large windows and can be used with any data restriction. In terms of the two approaches 1 and 2 described in Sec. 3.2, approach 2 does not involve, in the case of the algebraic method, any assumptions about the discarded data points, while, in the analytical reconstruction method, it is assumed that these data points are zero. On the other hand, the algebraic method requires explicit volume discretization in terms of voxels, while the analytical method allows one to reconstruct the target on any grid without much additional effort and is not dependent in any way on volume discretization. For the purpose of obtaining algebraic reconstructions, the reconstructed volume was divided into cubic voxels. The voxel size h was taken to be equal to 8 CCD pixels p. Thus, h=8×0.416 mm≃3.3 mm. The grid consisted of 41×41 voxels in the lateral direction and 17 voxels in the depth direction. Therefore, the discretized volume was a parallelepiped with the dimensions 13.6×13.6×4.3 cm3 consisting of N=21853 voxels. This parallelepiped was positioned from each surface of the slab. The target was situated approximately in the middle of the discretized volume. The discretization described above, and the approximation of the integral in Eq. (5) with a Riemann sum, results in a system of algebraic equations Ax=b where Amn is the M×N weight matrix, where M is the number of distinct source-detector pairs, N is the number of voxels, xn=δα(rn)/α0 is the vector of dimensionless contrast (n=1,…,N) and bm is the m-th data point (m=1,…,M). The equations are cast in dimensionless form by defining a dimensionless Green’s functions according to ˜G0(r,r′)=D0hG0(r,r′), where D0 is the diffusion coefficient in the Intralipid solution. Then we have Here kd=√α0/D0, rn is the position of the center of the n-th voxel while rdm, rsm are the detector and the source positions of the m-th data point used in the reconstruction.The pseudoinverse solution to the above system of equations is defined as the unique solution to the system (A*A+λ2I)x=A*b, where λ is the regularization parameter and I is the identity matrix. In our experiments, the number of measurements M is much larger than the number of voxels N (e.g., M∼107 and N∼104). Correspondingly, the most time consuming part of finding the pseudoinverse (at least, in the numerical approach used by us) is the computation of the matrix product A*A. In the most challenging case considered, we have used M=2×107 and N=2×104, which requires 8×1015 floating point operations. On an 8-core Xeon workstation with the peak performance of 56 GFlops, this translates into 40 h of computation time. However, this time-consuming procedure must be repeated only once for a given source-detector arrangement and given optical properties of the background medium (Intralipid). The latter did not change in the experiments reported in this paper. Thus, the resultant matrix A*A, once computed, was stored on a hard drive and re-used for image reconstruction with each new data set obtained, e.g., for various positions of the chest wall phantom. Note that computation of the projection, that is, of the N-component vector A*b involves only one matrix-vector multiplication and its computational cost is insignificant. The computation of A*A can be greatly accelerated by the use of the method proposed in Ref. 50. In this approach, the detectors are sampled for the purpose of computing the product A*A, but not for computing the projection A*b. In this way, the number of data points and the voxels used is not reduced but the computation time is shortened dramatically with no or minimal effects on the image quality. Indeed, we have verified that A*A can be computed by using the method in Ref. 50 in about 2 h with very minimal degradation of image quality. However, the questions of computational efficiency are largely outside of the scope of this paper, and we will adduce, therefore, only the images obtained by computing A*A without additional approximations. Note that the matrix A*A is small enough to be diagonalized; its eigenvectors and eigenvalues can be also stored on a hard drive for future use. However, in this study, we have solved the equation (A*A+λ2I)x=A*b directly by the conjugate-gradient descent method. An additional feature of the algebraic method described here is that it does not require that the set of detectors used be independent of the position of the source. We have taken advantage of this feature and have excluded the data points that are very far off-axis. Specifically, for each source, we have used only such detectors that are situated no further from the axis of the source than a given radius R. We have used R=6.25 cm, so that R is slightly larger than the width of the slab (6 cm). The justification for discarding the strongly off-axis data points is that these measurements contain predominantly noise. 3.4.Medium Parameters Used in the ReconstructionsTo perform quantitative reconstruction of the contrast, a few additional parameters must be specified. These parameters have been obtained by fitting the transmission function I0(rd,rs) in the homogeneous (reference) medium to the analytical prediction of the diffusion theory. The diffuse wave number kd was found to be equal to 0.72 cm−1. We note that the numerical value of D0 is not needed for the reconstructions as long as we are interested in the dimensionless relative contrast x(r)=δα(r)/α0. An additional parameter that enters the expression for G0(r,r′) is the extrapolation distance ℓ associated with the diffuse-nondiffuse boundary condition. The latter was found to be equal to 0.83 mm. 4.Data RestrictionThe rejection of the datapoints deemed unreliable or too noisy is an established practice in DOT.51–54 Here this data restriction methodology is used in a somewhat different context and with a somewhat different goal compared to previous work. In the usual approach, data points are rejected based on an a priori knowledge or a statistical model for the target without specific regard for the location of the rejected source-detector pairs. Here we reject certain data points based solely on their location. Our purpose is not to suppress noise, but rather to investigate the reconstruction effects of the imaging window restriction and the proximity of the chest wall phantom to the target. In particular, we are posing the following question: how close the sources and detectors can be placed to a chest wall to guarantee that the reconstruction of the target is not significantly distorted or contaminated by the latter. As discussed above, we denote the distance between the top of the bar target and the bottom of the chest wall phantom by d. We have collected the data for d=2, 5, 8, 11, 14, and 17 cm. The various data sets are graphically illustrated in Fig. 3. In particular, in Fig. 3(c), we show with red lines the positions of the lower edge of the chest wall phantom that correspond to different values of d and fall within the CCD FOV; these include d=8, 5, and 3 cm. In this figure, the positions of sources are shown by red dots. The drawing is to scale and a sample reconstruction is superimposed with the drawing to indicate the target shape and position. The larger, dark blue square corresponds to the FOV of the CCD camera. Fig. 3DownloadThe maximum available number of independent source-detector pairs is (512×35)2≃3.2×108. As discussed below, only a fraction of this data set was used in the reconstructions. Some data points have been eliminated by “windowing” (in the algebraic image reconstruction), other data points were eliminated by sampling of the detectors (in the algebraic reconstructions, every second detector was used and), and yet other data points have been eliminated by numerical data restriction, which is described below in this section. The maximum data set used in algebraic reconstruction consisted of ≃2×107 measurements; see more detailed information for various data restrictions in Table 1 below. Table 1Data restriction sizes.
We now discuss the data restriction. The latter was accomplished by removing all sources and detectors situated above one of the green lines shown in Fig. 3(c). In the case of algebraic reconstruction, these sources and detectors were simply not used, which did not amount to any additional approximation. In the case of the analytical reconstruction, it was assumed that the corresponding source-detector pairs produced zero data points bm (not zero intensity). The different data restrictions used are quantified as follows. There are 35 rows of sources. We define the quantity NR (Numerical data Restriction) as the number of the source row (counting from top to bottom) above which no sources and/or detectors are included in the reconstruction. Thus, if NR=5, the data are restricted to sources and detectors lying at the level of the fifth source row and below it. For example, NR=5 excludes the four topmost lines of sources and all detectors that lie above the fifth line of sources. In the reconstructions, we have used NR=5, 8, 9, 10 and, for each value of NR, we have performed image reconstruction with all available values of d. Table 1 summarizes the subsets of data used for each value of NR. 5.ResultsThe quantity plotted in all figures of this section is x(r)+1=α(r)/α0. From physical considerations, this function is nonnegative, since the medium is not amplifying. However, the image reconstruction reported here utilizes various approximations. This can result in reconstructing unphysical negative values of absorption, which are shown in the figures by the color black (the same color scale is used throughout). We note that the occurrence of negative absorption can be avoided by making use of a positivity constraint in the algebraic reconstruction method. The positivity constraint can be incorporated in the conjugate-gradient descent algorithm, which was used by us to invert the matrix A*A. However, we have found that the areas of negative absorption appear mostly in the case of the analytic (fast) image reconstruction method, which cannot incorporate the positivity constraint. On the other hand, the algebraic reconstructions have produced either no areas of negative absorption, or artifacts so severe (e.g., when d=2 cm and no numerical data restriction) that the use of the positivity constraint was not useful. In other words, we did not encounter a situation in which the positivity constraint was simultaneously numerically feasible and useful; therefore, it has not been used for producing the images shown in this section. Reconstructions of the central slice of the medium (3 cm from either of the slab surfaces) obtained with varying values of d and various numerical data restriction (NR) are shown in Fig. 4. It can be seen in Fig. 4(a) that the analytical inversion with no data restriction (the topmost row of images) produces severe image artifacts when chest wall is d=2 cm and d=5 cm away from the target. Data restriction with NR=5 results in a reasonable, yet suboptimal, image quality when d=5 cm, but not when d=2 cm. To remove the artifacts associated with the chest wall completely, NR=10 is required. Fig. 4DownloadHowever, NR=8, 9, 10 used in conjunction with the analytical reconstruction yields an additional image artifact, which is unrelated to the chest wall phantom. To see that this is true, consider the images for d=17 cm, which are not affected at all by the chest wall phantom, yet exhibit the additional artifact just mentioned. This artifact is shown as a black area where the reconstructed absorption coefficient is negative and, therefore, outside of the physically allowable range. We thus conclude that reconstructing the target by the analytical reconstruction method is feasible with the use of the appropriate data restriction, yet it results in an additional image artifact where the absorption in underestimated. The appearance of this artifact can be understood. As mentioned above, the data restriction used with the analytical reconstruction amounts to assuming that the truncated data points are zero. In other words, we assume that, in the presence of the target, the truncated source-detector pairs would have measured the same intensity as in the homogeneous slab, so that I(rd,rs)=I0(rd,rs) for the truncated source-detector pairs [see Eq. (4)]. The reconstruction algorithm seeks a contrast function δα(r), which is compatible with this assumption. For a purely absorbing target, however, the actual intensity I(rd,rs) is smaller than I0(rd,rs) when at least one of the points rd, rs is located not too far from the target (in the lateral direction) due to increased optical absorption. Whenever such data points are discarded, an artifact with negative δα is produced by the reconstruction algorithm to compensate for the absorption in the target. It can be seen that this artifact is located between the target and the region of source-detector pairs, which have been discarded. Of course, this analysis applies to the case when the position and optical contrast of the target is known. In general, it may be difficult to predict the position of this artifact or to distinguish it from a true occurrence of negative δα. There may also be a spatial overlap of the artifact and a true inhomogeneity. We now turn to the algebraic reconstructions [Fig. 4(b)]. For the unrestricted data set, the image quality is still poor. However, when the data restriction is gradually introduced, the artifacts disappear. Thus, in the case NR=10 and d=2 cm (the image in the bottom right corner), the target is clearly visible, and the image quality is about the same as with the use of the unrestricted data set and d=17 cm. Thus, introduction of the data restriction does not result in additional image artifacts or quality degradation when the algebraic method is used. In Figs. 5 and 6, we show slices drawn through the medium at different depths. Figure 5 displays the results of the analytical image reconstruction for d=5 cm and d=2 cm and Fig. 6 displays analogous data obtained by the algebraic reconstruction. In addition, we show in the right-most column of images the reconstruction averaged over the depth of the sample (that is, over different slices). Note that all reconstructed slices were used for the purpose of averaging, not only those shown in the figures. Note that, in all cases, we have reconstructed 13 slices separated by the distance of 8h≈3.328 mm with the central slice located exactly in the mid-plane of the slab. The “average” reconstruction was obtained by computing the arithmetic average of the reconstructions in all 13 slices. Fig. 5DownloadFig. 6DownloadThese averaged (“projection”) images correspond to the usual radiological projections obtained with a parallel beam of X-rays. The qualitative conclusions that can be drawn from Figs. 5 and 6 are the same as above. The analytical reconstruction produces reasonable image quality for the smallest chest wall-target separation d=2 cm and NR=10 but at the cost of an additional image artifact. The algebraic reconstruction is free from this artifact, but underestimates the image contrast relative to the analytic method (see below). The depth resolution is slightly better in algebraic reconstructions but, overall, much worse than the lateral resolution. This is typical for DOT images. One interesting feature observed in both types of image reconstruction is the following. The projection images discussed above are, generally, more stable and exhibit reasonable quality even when the individual slices contain severe artifacts. For example, consider the d=2 cm algebraic reconstructions without data restriction (Fig. 6). Even though all slices drawn through the medium are badly corrupted by the artifacts associated with proximity of the chest wall phantom, the projection image shows the target clearly. Moreover, the edge of the chest wall phantom is also clearly visible at the correct location. This result is somewhat unexpected and can be useful in the situations when the depth resolution is not of essence. We emphasize that obtaining the projections still requires knowledge of the three-dimensional distribution of the absorption coefficient; the projections cannot be computed or measured directly without such knowledge. We note that, in both types of image reconstructions, we see an underestimation of the contrast for the target phantom compared to the expected value. This underestimation can be attributed to the poor transverse (depth) resolution of the three-dimensional reconstruction which results in the “spreading” of the contrast in that direction. Indeed, consider the depth-integrated contrast, H(x,y)=∫[α(x,y,z)/α0−1]dz, where x, y are the coordinates in the plane of the slab and z is the transverse (depth) coordinate. Inside the target, we have α(x,y,z)/α0≃4 and the target thickness in the transverse direction is Δz=0.6 cm. Therefore, the actual value of H for a line passing through the target and perpendicularly to the slab surface is H≃1.8 cm. In the reconstructed images, the transverse thickness of the target is overestimated and is equal, approximately, to 2 cm while the quantity α(x,y,z)/α0 is underestimated and is equal, approximately, to 2. By using the reconstructed values to estimate the integrated contrast, we obtain H≃2 cm, which is reasonably close to the actual value. 6.Summary and DiscussionWe have used phantom experiments to investigate systematically the effects of the chest wall on diffusion optical tomography (DOT) of the breast. The results lead us to several promising conclusions. First, we have found that, when absorption contrast is of interest, simple CW instrumentation with linearized inversion can suffice. This finding was obtained in spite of the presence of the chest wall phantom in close proximity to the target, which first renders the inverse problem nonlinear and, second, differs from the background Intralipid and the target not only in absorption but also in scattering properties. Generally, under the conditions stated above, time- or frequency-resolved measurements and nonlinear image-reconstruction methods are required. We, however, have been able to bypass these complications by appropriately restricting the data points used in the reconstruction. We note that in clinical applications of DOT, the location of the chest wall relative to the sources and detectors is usually known; therefore, the approach of this paper to data restriction can be applied in vivo. Work remains, however, to optimize these approaches. This may require more sophisticated and/or data-driven algorithms for data rejection, as well as experimentation in vivo. In this paper, we have demonstrated that the rather severe effects of the chest wall can, in principle, be rectified by appropriate data restriction in conjunction with a linear image reconstruction. The paper shows that, by means of properly restricting the data points used in image reconstruction, it is possible to resolve a small absorptive target in the vicinity of a spatially and optically large inhomogeneity and that the quality of the reconstruction is almost unaffected by the chest wall. Interestingly, we have also found (see Figs. 5 and 6) that the image contrast, when averaged over the depth of a plane-parallel sample (we refer to this quantity as the projection), is not as sensitive to systematic errors encountered in image reconstruction as the individual slices drawn through the medium. Thus, under certain conditions, the nonlinearity of the inverse problem and the presence of a scattering contrast render our image reconstruction methods inadequate. Reconstructed slices show severe image artifacts in this case. The projection, however, is free from these artifacts and displays the target clearly. This finding may be significant since a projection of the type just discussed is similar to the usual radiographic projection, yet it displays the contrast specific to the near-infrared spectral range. Finally, we have developed and verified with experimental data an algebraic image reconstruction method, which is well suited for the use with the data sets restricted by the presence of the chest wall and capable of handling data sets as large as 2×107 independent measurements. AcknowledgmentsWe are grateful to Soren D. Konecky for valuable advice on the experimental design. This project was partially supported by the National Center for Research Resources and the National Institute of Biomedical Imaging and Bioengineering of the National Institutes of Health through grants R01 EB002109, P41 RR002305 and P41-EB015893, and by the National Science Foundation through grants DMS-1115616, DMS-1108969 and DMS-1115574. ReferencesS. R. Arridge,
“Optical tomography in medical imaging,”
Inverse Probl., 15
(2), R41
–R93
(1999). http://dx.doi.org/10.1088/0266-5611/15/2/022 INPEEY 0266-5611 Google Scholar
D. A. Boaset al.,
“Imaging the body with diffuse optical tomography,”
IEEE Signal Proc. Mag., 18
(6), 57
–75
(2001). http://dx.doi.org/10.1109/79.962278 ISPRE6 1053-5888 Google Scholar
S. R. ArridgeJ. C. Schotland,
“Optical tomography: forward and inverse problems,”
Inverse Probl., 25
(12), 123010
(2009). http://dx.doi.org/10.1088/0266-5611/25/12/123010 INPEEY 0266-5611 Google Scholar
S. B. Colaket al.,
“Clinical optical tomography and NIR spectroscopy for breast cancer detection,”
IEEE J. Sel. Top. Quantum Electron., 5
(4), 1143
–1158
(1999). http://dx.doi.org/10.1109/2944.796341 IJSQEN 1077-260X Google Scholar
D. J. HawryszE. M. Sevick-Muraca,
“Developments toward diagnostic breast cancer imaging using near-infrared optical measurements and fluorescent contrast agants,”
Neoplasia, 2
(5), 388
–417
(2000). http://dx.doi.org/10.1038/sj.neo.7900118 1522-8002 Google Scholar
T. McBrideet al.,
“A parallel-detection frequency-domain near-infrared tomography system for hemoglobin imaging of the breast in vivo,”
Rev. Sci. Instrum., 72
(3), 1817
–1825
(2001). http://dx.doi.org/10.1063/1.1344180 RSINAK 0034-6748 Google Scholar
J. P. Culveret al.,
“Three-dimensional diffuse optical tomography in the parallel plane transmission geometry: evaluation of a hybrid frequency domain/continuous wave clinical system for breast imaging,”
Med. Phys., 30
(2), 235
–247
(2003). http://dx.doi.org/10.1118/1.1534109 MPHYA6 0094-2405 Google Scholar
X. Inteset al.,
“In vivo continuous-wave optical breast imaging enhanced with indocyanine green,”
Med. Phys., 30
(6), 1039
–1047
(2003). http://dx.doi.org/10.1118/1.1573791 MPHYA6 0094-2405 Google Scholar
A. Liet al.,
“Tomographic optical breast imaging guided by three-dimensional mammography,”
Appl. Opt., 42
(25), 5181
–5190
(2003). http://dx.doi.org/10.1364/AO.42.005181 APOPAI 0003-6935 Google Scholar
R. Choeet al.,
“Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,”
Med. Phys., 32
(4), 1128
–1139
(2005). http://dx.doi.org/10.1118/1.1869612 MPHYA6 0094-2405 Google Scholar
T. Yateset al.,
“Time-resolved optical mammography using a liquid coupled interface,”
J. Biomed. Opt., 10
(5), 054011
(2005). http://dx.doi.org/10.1117/1.2063327 JBOPFO 1083-3668 Google Scholar
A. Corluet al.,
“Three-dimensional in vivo fluorescence diffuse optical tomography of breast cancer in humans,”
Opt. Express, 15
(11), 6696
–6716
(2007). http://dx.doi.org/10.1364/OE.15.006696 OPEXFF 1094-4087 Google Scholar
A. Cerussiet al.,
“Predicting response to breast cancer neoadjuvant chemotherapy using diffuse optical spectroscopy,”
Proc. Natl. Acad. Sci. U. S. A., 104
(10), 4014
–4019
(2007). http://dx.doi.org/10.1073/pnas.0611058104 Google Scholar
R. Choeet al.,
“Differentiation of benign and malignant breast tumors by in-vivo three-dimensional parallel-plate diffuse optical tomography,”
J. Biomed. Opt., 14
(2), 024020
(2009). http://dx.doi.org/10.1117/1.3103325 JBOPFO 1083-3668 Google Scholar
T. Durduranet al.,
“Diffuse optics for tissue monitoring and tomography,”
Rep. Prog. Phys., 73
(7), 076701
(2010). http://dx.doi.org/10.1088/0034-4885/73/7/076701 RPPHAG 0034-4885 Google Scholar
B. Pogueet al.,
“Initial assessment of a simple system for frequency domain diffuse optical tomography,”
Phys. Med. Biol., 40
(10), 1709
–1729
(1995). http://dx.doi.org/10.1088/0031-9155/40/10/011 PHMBA7 0031-9155 Google Scholar
D. Grosenicket al.,
“Time-domain scanning optical mammography: I. recording and assessment of mammograms of 154 patients,”
Phys. Med. Biol., 50
(11), 2429
–2449
(2005). http://dx.doi.org/10.1088/0031-9155/50/11/001 PHMBA7 0031-9155 Google Scholar
A. Pifferiet al.,
“Four-wavelength time-resolved optical mammography in the 680 980-nm range,”
Opt. Lett., 28
(13), 1138
–1140
(2003). http://dx.doi.org/10.1364/OL.28.001138 OPLEDP 0146-9592 Google Scholar
J. GeB. ZhuS. RegaladoA. Godavarty,
“Three-dimensional fluorescence-enhanced optical tomography using a hand-held probe based imaging system,”
Med. Phys., 35
(7), 3354
–3363
(2008). http://dx.doi.org/10.1118/1.2940603 MPHYA6 0094-2405 Google Scholar
Z.-M. Wanget al.,
“Experimental demonstration of an analytic method for image reconstruction in optical tomography with large data sets,”
Opt. Lett., 30
(24), 3338
–3340
(2005). http://dx.doi.org/10.1364/OL.30.003338 OPLEDP 0146-9592 Google Scholar
S. D. Koneckyet al.,
“Imaging complex structures with diffuse light,”
Opt. Express, 16
(7), 5048
–5060
(2008). http://dx.doi.org/10.1364/OE.16.005048 OPEXFF 1094-4087 Google Scholar
P. Bonfert-Tayloret al.,
“Information loss and reconstruction in diffuse fluorescence tomography,”
J. Opt. Soc. Am. A, 29
(3), 321
–330
(2012). http://dx.doi.org/10.1364/JOSAA.29.000321 JOAOD6 0740-3232 Google Scholar
S. ChaillatG. Biros,
“FaIMS: a fast algorithm for the inverse medium problem with multiple frequencies and multiple sources for the scalar Helmholtz equation,”
J. Comput. Phys., 231
(12), 4403
–4421
(2012). http://dx.doi.org/10.1016/j.jcp.2012.02.006 JCTPAH 0021-9991 Google Scholar
R. SchulzJ. RipollV. Ntziachristos,
“Noncontact optical tomography of turbid media,”
Opt. Lett., 28
(18), 1701
–1703
(2003). http://dx.doi.org/10.1364/OL.28.001701 OPLEDP 0146-9592 Google Scholar
J. RipollR. B. SchulzV. Ntziachristos,
“Free-space propagation of diffuse light: theory and experiments,”
Phys. Rev. Lett., 91
(10), 103901
(2003). http://dx.doi.org/10.1103/PhysRevLett.91.103901 PRLTAO 0031-9007 Google Scholar
J. RipollV. Ntziachristos,
“Imaging scattering media from a distance: theory and applications of noncontact optical tomography,”
Mod. Phys. Lett., 18
(28–29), 1403
–1431
(2004). http://dx.doi.org/10.1142/S0217984904007864 MPLAEQ 0217-7323 Google Scholar
G. Turneret al.,
“Complete-angle projection diffuse optical tomography by use of early photons,”
Opt. Lett., 30
(4), 409
–411
(2005). http://dx.doi.org/10.1364/OL.30.000409 OPLEDP 0146-9592 Google Scholar
V. A. MarkelJ. C. Schotland,
“Inverse scattering for the diffusion equation with general boundary conditions,”
Phys. Rev. E, 64
(3), R035601
(2001). http://dx.doi.org/10.1103/PhysRevE.64.035601 PLEEE8 1063-651X Google Scholar
V. A. MarkelJ. C. Schotland,
“Effects of sampling and limited data in optical tomography,”
Appl. Phys. Lett., 81
(7), 1180
–1182
(2002). http://dx.doi.org/10.1063/1.1495543 APPLAB 0003-6951 Google Scholar
V. A. MarkelV. MitalJ. C. Schotland,
“The inverse problem in optical diffusion tomography. III. Inversion formulas and singular value decomposition,”
J. Opt. Soc. Am. A, 20
(5), 890
–902
(2003). http://dx.doi.org/10.1364/JOSAA.20.000890 JOAOD6 0740-3232 Google Scholar
V. A. MarkelJ. A. O’SullivanJ. C. Schotland,
“Inverse problem in optical diffusion tomography. IV. Nonlinear inversion formulas,”
J. Opt. Soc. Am. A, 20
(5), 903
–912
(2003). http://dx.doi.org/10.1364/JOSAA.20.000903 JOAOD6 0740-3232 Google Scholar
V. A. MarkelJ. C. Schotland,
“Symmetries, inversion formulas, and image reconstruction for optical tomography,”
Phys. Rev. E, 70
(5), 056616
(2004). http://dx.doi.org/10.1103/PhysRevE.70.056616 PLEEE8 1063-651X Google Scholar
B. W. PogueM. S. Patterson,
“Review of tissue simulating phantoms for optical spectroscopy, imaging and dosimetry,”
J. Biomed. Opt., 11
(4), 041102
(2006). http://dx.doi.org/10.1117/1.2335429 JBOPFO 1083-3668 Google Scholar
A. E. Cerussiet al.,
“Tissue phantoms in multicenter clinical trials for diffuse optical technologies,”
Biomed. Opt. Express, 3
(5), 966
–971
(2012). http://dx.doi.org/10.1364/BOE.3.000966 BOEICL 2156-7085 Google Scholar
Y. ArdeshirpourQ. Zhu,
“Optical tomography method that accounts for tilted chest wall in breast imaging,”
J. Biomed. Opt., 15
(4), 041515
(2010). http://dx.doi.org/10.1117/1.3449570 JBOPFO 1083-3668 Google Scholar
A. KienleT. Glanzmann,
“In vivo determination of the optical properties of muscle with time-resolved reflectance using a layered model,”
Phil. Mag. B, 44
(11), 2689
–2702
(1999). http://dx.doi.org/10.1088/0031-9155/44/11/301 0141-8637 Google Scholar
P. Taroniet al.,
“In vivo absorption and scattering spectroscopy of biological tissues,”
Photochem. Photobiol. Sci., 2
(2), 124
–129
(2003). http://dx.doi.org/10.1039/b209651j PPSHCB 1474-905X Google Scholar
M. S. PattersonB. ChanceB. C. Wilson,
“Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,”
Appl. Opt., 28
(12), 2331
–2336
(1989). http://dx.doi.org/10.1364/AO.28.002331 APOPAI 0003-6935 Google Scholar
D. A. BenaronD. K. Stevenson,
“Optical time-of-flight and absorbance imaging of biologic media,”
Science, 259
(5100), 1463
–1466
(1993). http://dx.doi.org/10.1126/science.8451643 SCIEAS 0036-8075 Google Scholar
S. Andersson-Engelset al.,
“Time-resolved transillumination for medical diagnostics,”
Opt. Lett., 15
(21), 1179
–1181
(1990). http://dx.doi.org/10.1364/OL.15.001179 OPLEDP 0146-9592 Google Scholar
S. L. Jacques,
“Time-resolved reflectance spectroscopy in turbid tissues,”
IEEE Trans. Biomed. Eng., 36
(12), 1155
–1161
(1989). http://dx.doi.org/10.1109/10.42109 IEBEAX 0018-9294 Google Scholar
F. E. Schmidtet al.,
“Multiple-slice imaging of a tissue-equivalent phantom by use of time-resolved optical tomography,”
Appl. Opt., 39
(19), 3380
–3387
(2000). http://dx.doi.org/10.1364/AO.39.003380 APOPAI 0003-6935 Google Scholar
V. NtziachristosX. H. MaB. Chance,
“Time-correlated single photon counting imager for simultaneous magnetic resonance and near-infrared mammography,”
Rev. Sci. Instrum., 69
(12), 4221
–4233
(1998). http://dx.doi.org/10.1063/1.1149235 RSINAK 0034-6748 Google Scholar
E. Grattonet al.,
“The possibility of a near-infrared optical imaging system using frequency-domain methods,”
in Proc. 3rd Int. Conf. on Peace through Mind@ Brain Science,
183
–189
(1990). Google Scholar
J. B. FishkinE. Gratton,
“Propagation of photon-density waves in strongly scattering media containing an absorbing semiinfinite plane bounded by a straight edge,”
J. Opt. Soc. Am. A, 10
(1), 127
–140
(1993). http://dx.doi.org/10.1364/JOSAA.10.000127 JOAOD6 0740-3232 Google Scholar
B. Chanceet al.,
“Phase measurement of light absorption and scattering in human tissues,”
Rev. Sci. Instrum., 69
(10), 3457
–3481
(1998). http://dx.doi.org/10.1063/1.1149123 RSINAK 0034-6748 Google Scholar
B. W. PogueM. S. Patterson,
“Frequency-domain optical absorption spectroscopy of finite tissue volumes using diffusion theory,”
Phys. Med. Biol., 39
(7), 1157
–1180
(1994). http://dx.doi.org/10.1088/0031-9155/39/7/008 PHMBA7 0031-9155 Google Scholar
J. C. Schotland,
“Continuous wave diffusion imaging,”
J. Opt. Soc. Am. A, 14
(1), 275
–279
(1997). http://dx.doi.org/10.1364/JOSAA.14.000275 JOAOD6 0740-3232 Google Scholar
C. P. Gonataset al.,
“Optical diffusion imaging using a direct inversion method,”
Phys. Rev. E, 52
(4), 4361
–4365
(1995). http://dx.doi.org/10.1103/PhysRevE.52.4361 PLEEE8 1063-651X Google Scholar
V. A. MarkelZ.-M. WangJ. C. Schotland,
“Optical diffusion tomography with large data sets,”
Proc. SPIE, 5969 59691B
(2005). http://dx.doi.org/10.1117/12.628623 PSISDG 0277-786X Google Scholar
A. Blasiet al.,
“Investigation of depth dependent changes in cerebral haemodynamics during face preception in infants,”
Phys. Med. Biol., 52
(23), 6849
–6864
(2007). http://dx.doi.org/10.1088/0031-9155/52/23/005 0141-8637 Google Scholar
M. A. Franceschiniet al.,
“Assesment of infant brain development with frequency-domain near-infrared spectroscopy,”
Proc. Natl. Acad. Sci. U. S. A., 94
(12), 6468
–6473
(1997). http://dx.doi.org/10.1073/pnas.94.12.6468 PNASA6 0027-8424 Google Scholar
N. Roche-Labarbeet al.,
“Noninvasive optical measures of CBV, StO2, C.BF index, and rCMRO2 in human premature neonates' brains in the first six weeks of life,”
Hum. Brain Map, 31
(3), 341
–352
(2010). http://dx.doi.org/10.1002/hbm.v31:3 HBRME7 1065-9471 Google Scholar
F. Orihuela-Espinaet al.,
“Quality control and assurance in functional near infrared spectroscopy (fNIRS) experimentation,”
Phys. Med. Biol., 55
(13), 3701
–3724
(2010). http://dx.doi.org/10.1088/0031-9155/55/13/009 0141-8637 Google Scholar
|