Open Access
7 February 2013 Diffuse optical tomography in the presence of a chest wall
Author Affiliations +
Abstract
Diffuse optical tomography (DOT) has been employed to derive spatial maps of physiologically important chromophores in the human breast, but the fidelity of these images is often compromised by boundary effects such as those due to the chest wall. We explore the image quality in fast, data-intensive analytic and algebraic linear DOT reconstructions of phantoms with subcentimeter target features and large absorptive regions mimicking the chest wall. Experiments demonstrate that the chest wall phantom can introduce severe image artifacts. We then show how these artifacts can be mitigated by exclusion of data affected by the chest wall. We also introduce and demonstrate a linear algebraic reconstruction method well suited for very large data sets in the presence of a chest wall.

1.

Introduction

Diffuse optical tomography (DOT) employs near-infrared light to probe the optical properties of highly scattering biological tissues.13 DOT has shown promise in breast cancer imaging415 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,2022 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,2427 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).2832 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.

Experiment

The 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.6cm2 area on one side of the imaging tank (whose overall dimensions are 44×44×6cm3). 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.2cm2 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.416mm.

Fig. 1

Schematic of the experimental setup. A collimated CW 785 nm laser source at is raster scanned on one side of the imaging tank. The transmitted light on the detection plane is collected by a CCD for each source position.

JBO_18_2_026016_f001.png

A 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.2cm1 and the reduced scattering coefficient μs=7.5cm1. The tank is filled with a scattering fluid (μa=0.05cm1 and μs=7.5cm1); 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.1cm1 and μs=5.0cm1, dimensions 40×20×5.8cm3) 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.3537 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. 2

Phantoms used in the experiment. (a) 6 mm thick bar target with μa=0.2cm1 and μs=7.5cm1 has slots 48 mm tall and 9 mm wide. The outer dimensions are 60×50mm2. (b) The chest wall phantom with μa=0.1cm1 and μs=5.0cm1.

JBO_18_2_026016_f002.png

3.

Image Reconstruction

3.1.

Linearized Integral Equations

This research employs linear image reconstruction methods and CW data. In principle, one could also resort to time-3843 or frequency-resolved4447 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:

Eq. (1)

I(rd,rs)=A(rd)B(rs)G(rd,rs),I0(rd,rs)=A(rd)B(rs)G0(rd,rs).
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

Eq. (2)

G(rd,rs)=G0(rd,rs)VG0(rd,r)δα(r)G(r,rs)d3r.
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

Eq. (3)

G(rd,rs)=G0(rd,rs)exp[VG0(rd,r)δα(r)G(r,rs)d3rG0(rd,rs)].
Consequently, we define the data function ϕ(rd,rs) according to

Eq. (4)

ϕ(rd,rs)=G0(rd,rs)ln[I(rd,rs)I0(rd,rs)]
and use Eq. (3) to obtain an integral equation of the form

Eq. (5)

VG0(rd,r)δα(r)G0(r,rs)d3r=ϕ(rd,rs).
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 Method

This 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=2cm 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:

  • 1. Use all data points available (i.e., in a wide window on both sides of the slab). This scheme will include, in some cases, data points strongly affected by the presence of the chest wall. In this case, the analytical image reconstruction algorithm is applied as intended by mathematical design, i.e., without additional approximations. However, the strong nonlinearity of the inverse problem due to the presence of the chest wall phantom renders the underlying first Rytov approximation invalid. This, in turn, results in poor-quality reconstructions or inability to reconstruct the target at all. From the practical point of view, this approach is simply not feasible in vivo because the data, which are collected in our experiment over large windows, are physically unavailable in the case of a real chest wall.

  • 2. Truncate the data function by replacing data points, which are either deemed as physically unavailable in vivo or as too badly contaminated by the presence of the chest wall, with zeroes. This is equivalent to replacing a subset of optical measurements obtained in the heterogeneous medium with the respective measurements obtained in the homogeneous (reference) medium. Approach 2 is numerically efficient and directly applicable in vivo. In fact, we will demonstrate below that it produces images of reasonable quality, i.e., of much better quality than approach 1. However, the additional approximation involved in this approach results in the appearance of image artifacts that are poorly controlled and may be viewed as undesirable.

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/2p). The variables p and q were sampled as follows: qij=Δ(ix^+jy^), pkl=Δ(kx^+ly^), where Δ=2π/(331×p), 28i,j28 and 0k, l6. 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×32491.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 Method

Algebraic 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.416mm3.3mm. 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.3cm3 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 G˜0(r,r)=D0hG0(r,r), where D0 is the diffusion coefficient in the Intralipid solution. Then we have

Eq. (6)

Amn=(kdh)2G˜0(rdm,rn)G˜0(rn,rsm),bm=G˜0(rdm,rsm)ln[I(rdm,rsm)I0(rdm,rsm)].
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., M107 and N104). 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.25cm, 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 Reconstructions

To 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.72cm1. 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 Restriction

The rejection of the datapoints deemed unreliable or too noisy is an established practice in DOT.5154 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. 3

Models for data restriction. (a) Photograph of the drained imaging tank illustrating the position the target with respect to the chest wall phantom. (b) Schematic of the imaging tank. (c) Illustration of the various data sets used in the reconstructions. The dark blue square is the CCD FOV while the inner light blue square indicates the reconstruction region. The red dots indicate the source positions. A sample reconstruction is superimposed with the drawing to illustrate the target shape and position. The red lines indicate the three lowest position of the chest wall phantom (other positions are outside of the CCD FOV) while the green lines illustrate the restricted data sets where all the sources and detectors situated above a given green line have been discarded.

JBO_18_2_026016_f003.png

The maximum available number of independent source-detector pairs is (512×35)23.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 1

Data restriction sizes.

Numerical data restriction (NR)Number of sourcesNumber of distinct source-detector pairs
No restriction35×35=122520591492
535×31=108516479152
835×28=98014661145
935×27=94514074728
1035×26=91013457507

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.

Results

The 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=2cm 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=2cm and d=5cm away from the target. Data restriction with NR=5 results in a reasonable, yet suboptimal, image quality when d=5cm, but not when d=2cm. To remove the artifacts associated with the chest wall completely, NR=10 is required.

Fig. 4

Images of the central slice obtained by analytical (a) and algebraic (b) reconstruction methods. Different columns show data obtained with the chest wall phantoms at different distances d from the bar target. Different rows of images correspond to different data restrictions NR, as indicated. The color bar applies to all other images shown below.

JBO_18_2_026016_f004.png

However, 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=17cm, 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=2cm (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=17cm. 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=5cm and d=2cm 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 8h3.328mm 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. 5

Slices through the medium drawn at different depths (from the plane of sources) as indicated. Analytical image reconstruction method with d=5cm and d=2cm.

JBO_18_2_026016_f005.png

Fig. 6

Same as in Fig. 5 but obtained by algebraic reconstruction.

JBO_18_2_026016_f006.png

These 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=2cm 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=2cm 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)/α01]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)/α04 and the target thickness in the transverse direction is Δz=0.6cm. Therefore, the actual value of H for a line passing through the target and perpendicularly to the slab surface is H1.8cm. 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 H2cm, which is reasonably close to the actual value.

6.

Summary and Discussion

We 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.

Acknowledgments

We 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.

References

1. 

S. 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

2. 

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

3. 

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

4. 

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

5. 

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

6. 

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

7. 

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

8. 

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

9. 

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

10. 

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

11. 

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

12. 

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

13. 

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

14. 

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

15. 

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

16. 

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

17. 

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

18. 

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

19. 

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

20. 

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

21. 

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

22. 

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

23. 

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

24. 

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

25. 

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

26. 

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

27. 

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

28. 

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

29. 

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

30. 

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

31. 

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

32. 

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

33. 

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

34. 

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

35. 

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

36. 

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

37. 

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

38. 

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

39. 

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

40. 

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

41. 

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

42. 

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

43. 

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

44. 

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

45. 

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

46. 

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

47. 

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

48. 

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

49. 

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

50. 

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

51. 

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

52. 

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

53. 

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

54. 

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
© 2013 Society of Photo-Optical Instrumentation Engineers (SPIE) 0091-3286/2013/$25.00 © 2013 SPIE
Han Y. Ban, David R. Busch, Saurav Pathak, Frank A. Moscatelli, Manabu Machida, John C. Schotland, Vadim A. Markel, and Arjun G. Yodh "Diffuse optical tomography in the presence of a chest wall," Journal of Biomedical Optics 18(2), 026016 (7 February 2013). https://doi.org/10.1117/1.JBO.18.2.026016
Published: 7 February 2013
Lens.org Logo
CITATIONS
Cited by 6 scholarly publications and 1 patent.
Advertisement
Advertisement
RIGHTS & PERMISSIONS
Get copyright permission  Get copyright permission on Copyright Marketplace
KEYWORDS
Chest

Reconstruction algorithms

Image restoration

Diffuse optical tomography

Sensors

Absorption

Image quality


CHORUS Article. This article was made freely available starting 07 February 2014

Back to Top