Diffuse optical tomography in thepresence of a chest wall
Han Y. Ban
David R. Busch
Saurav PathakFrank A. Moscatelli
Manabu MachidaJohn C. Schotland
Vadim A. Markel
Arjun G. Yodh

Downloaded From: http://biomedicaloptics.spiedigitallibrary.org/ on 02/11/2013 Terms of Use: http://spiedl.org/terms

Journal of Biomedical Optics 18(2), 026016 (February 2013)

Diffuse optical tomography in the presence of a chest wall
Han Y. Ban,a David R. Busch,a Saurav Pathak,a Frank A. Moscatelli,b Manabu Machida,c John C. Schotland,cVadim A. Markel,d and Arjun G. Yodhaa

University of Pennsylvania, Department of Physics and Astronomy, Philadelphia, Pennsylvania 19104-6396
Swarthmore College, Department of Physics and Astronomy, Swarthmore, Pennsylvania 19081
cUniversity of Michigan, Department of Mathematics, Ann Arbor, Michigan 48109-1043
dUniversity of Pennsylvania, Departments of Radiology, Bioengineering and Graduate Group in Applied Mathematics and Computational Science,
Philadelphia, Pennsylvania 19104b

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 suchas 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 thechest wall. Experiments demonstrate that the chest wall phantom can introduce severe image artifacts. We thenshow how these artifacts can be mitigated by exclusion of data affected by the chest wall. We also introduce anddemonstrate a linear algebraic reconstruction method well suited for very large data sets in the presence of a chest
wall. © 2013 Society of Photo-Optical Instrumentation Engineers (SPIE) [DOI: 10.1117/1.JBO.18.2.026016]
Keywords: image reconstruction; imaging systems; tomography.
Paper 12507 received Aug. 8, 2012; revised manuscript received Dec. 7, 2012; accepted for publication Jan. 14, 2013; published online
Feb. 7, 2013.

1

Introduction

Diffuse optical tomography (DOT) employs near-infrared lightto probe the optical properties of highly scattering biologicaltissues.1–3 DOT has shown promise in breast cancer imaging4–15 because optical methods are sensitive to changes inphysiological parameters such as blood volume and tissue oxygenation, alterations of which are biomarkers for cancer andother 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 havebeen used, though transmission measurements provide moresensitivity 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 availabilityof 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 morethan ∼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 potentialof these methods, but the simulations also indicate that large
Address all correspondence to: Han Y. Ban, University of Pennsylvania,
Department of Physics and Astronomy, Philadelphia, Pennsylvania 19104-6396.
Tel: +215-898-8422; E-mail: hban@physics.upenn.edu.

Journal of Biomedical Optics

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 andother 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 arepoorly 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 theutility 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 analyticaland algebraic data-intensive linearized image reconstruction
methods can produce reasonable results, provided the data
points are appropriately restricted to exclude measurementsthat are strongly influenced by the chest wall. Under these conditions, an absorbing target with subcentimeter features can beclearly reconstructed in the middle of a 6 cm slab, even when thechest wall is only 2 cm from the target. We obtain good imagesof the target even in the presence of a large chest wall phantom
0091-3286/2013/$25.00 © 2013 SPIE

026016-1

Downloaded From: http://biomedicaloptics.spiedigitallibrary.org/ on 02/11/2013 Terms of Use: http://spiedl.org/terms

February 2013 • Vol. 18(2)Ban et al.: Diffuse optical tomography in the presence of a chest wallthat introduces significant nonlinearities into the inverse problem due to its larger absorption coefficient compared to thebackground as well as due to its size. Moreover, we discovereda data restriction condition such that the presence of the chestwall phantom imposes minimal artifacts or distortions in theimage. 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 rolefor 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 detailsof the image reconstruction methods. Section 4 explains ourapproach to data restriction. Section 5 presents the results andSec. 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 sizeof 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 positionis illuminated, data is collected from the opposite side of thetank 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 thetank with the spacing p ¼ 0.416 mm.
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.2 cm−1 and the reduced scattering coefficient
μs0 ¼ 7.5 cm−1 . The tank is filled with a scattering fluid (μa ¼0.05 cm−1 and μs0 ¼ 7.5 cm−1 ); these background optical properties are similar to those used in previous in vitro and clinicalresearch. 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μs0 ¼ 5.0 cm−1 , dimensions 40 × 20 × 5.8 cm3 ) is suspended atvarious distances d from the top edge of the bar target (d ¼ 2, 5,

Fig. 2 Phantoms used in the experiment. (a) 6 mm thick bar target withμa ¼ 0.2 cm−1 and μs0 ¼ 7.5 cm−1 has slots 48 mm tall and 9 mm wide.
The outer dimensions are 60 × 50 mm2 . (b) The chest wall phantom
with μa ¼ 0.1 cm−1 and μs0 ¼ 5.0 cm−1 .

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 thechest wall phantom are shown in Fig. 2. Note that the chestwall 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.

3

Image Reconstruction

3.1

Linearized Integral Equations

This research employs linear image reconstruction methods
and CW data. In principle, one could also resort totime-38–43 or frequency-resolved44–47 measurements and nonlinear reconstruction methods1,31 to obtain a reconstructionof 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 independentmeasurements of the transmitted intensity are taken, one ina homogeneous (reference) slab and the other in a slab withthe target and the chest wall phantom present. We denotethese measurements by I 0 ð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 therespective surfaces of the tank. In this work, I 0 and I are
expressed in CCD counts. These measurements can be related
to the medium optical properties through the relations:Iðrd ; rs Þ ¼ Aðrd ÞBðrs ÞGðrd ; rs Þ;
I 0 ðrd ; rs Þ ¼ Aðrd ÞBðrs ÞG0 ðrd ; rs Þ:

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 eachsource position.
Journal of Biomedical Optics(1)

Here Aðrd Þ and Bðrs Þ are unknown coupling coefficientsassociated with the source and detection system which canbe 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 isknown 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

026016-2

Downloaded From: http://biomedicaloptics.spiedigitallibrary.org/ on 02/11/2013 Terms of Use: http://spiedl.org/terms

February 2013 • Vol. 18(2)Ban et al.: Diffuse optical tomography in the presence of a chest wallZGðrd ; rs Þ ¼ G0 ðrd ; rs Þ −VG0 ðrd ; rÞδαðrÞGðr; rs Þd3 r:

reconstruct the target at all. From the practical pointof 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)

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 , wherec is the average speed of light in the medium.
The procedure of linearization consists of making anapproximation 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

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 presenceof the chest wall, with zeroes. This is equivalent toreplacing a subset of optical measurements obtainedin the heterogeneous medium with the respective measurements obtained in the homogeneous (reference)
medium. Approach 2 is numerically efficient anddirectly 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 inthis approach results in the appearance of image artifacts that are poorly controlled and may be viewed as
undesirable.

Gðrd ; rs Þ Z
G0 ðrd ; rÞδαðrÞGðr; rs Þd3 r¼ G0 ðrd ; rs Þ exp −:G0 ðrd ; rs ÞV(3)Consequently, we define the data function ϕðrd ; rs Þ according toIðrd ; rs Þϕðrd ; rs Þ ¼ −G0 ðrd ; rs Þ ln
I 0 ðrd ; rs Þ(4)

and use Eq. (3) to obtain an integral equation of the formZV

G0 ðrd ; rÞδαðrÞG0 ðr; rs Þd3 r ¼ ϕðrd ; rs Þ:(5)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 0 Þ, 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 tothe slab geometry with experimental data have been reported inRefs. 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 separationd between the chest wall and the top of the bar target. At the two
smallest separations (d ¼ 2 cm and 5 cm), contamination of thedata function by the chest wall is significant. Under these circumstances, two approaches for data analysis are natural toconsider:1. Use all data points available (i.e., in a wide window onboth sides of the slab). This scheme will include, insome cases, data points strongly affected by thepresence of the chest wall. In this case, the analyticalimage reconstruction algorithm is applied as intended
by mathematical design, i.e., without additional
approximations. However, the strong nonlinearityof the inverse problem due to the presence of thechest wall phantom renders the underlying first
Rytov approximation invalid. This, in turn, resultsin poor-quality reconstructions or inability to
Journal of Biomedical Optics

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 totalof 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 Method

Algebraic image reconstruction methods are obtained by discretization of the volume integral in Eq. (5) and computation of thepseudo-inverse for the obtained system of linearized equations.49 This method does not require large windows and canbe used with any data restriction. In terms of the two approaches
1 and 2 described in Sec. 3.2, approach 2 does not involve, inthe case of the algebraic method, any assumptions about thediscarded data points, while, in the analytical reconstruction
method, it is assumed that these data points are zero. On theother 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, thereconstructed 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 parallelepipedwith the dimensions 13.6 × 13.6 × 4.3 cm3 consisting of
N ¼ 21853 voxels. This parallelepiped was positioned from

026016-3

Downloaded From: http://biomedicaloptics.spiedigitallibrary.org/ on 02/11/2013 Terms of Use: http://spiedl.org/terms

February 2013 • Vol. 18(2)Ban et al.: Diffuse optical tomography in the presence of a chest wall

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 ofalgebraic equations Ax ¼ b where Amn is the M × N weightmatrix, 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 dimensionlessform by defining a dimensionless Green’s functions according
to G̃0 ðr; r 0 Þ ¼ D0 hG0 ðr; r 0 Þ, where D0 is the diffusion coefficient in the Intralipid solution. Then we haveAmn ¼ ðkd hÞ2 G̃0 ðrdm ; rn ÞG̃0 ðrn ; rsm Þ;(6)Iðrdm ; rsm Þ
bm ¼ −G̃0 ðrdm ; rsm Þ ln:
I 0 ðrdm ; rsm Þpﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
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þ
λ2 IÞ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 numericalapproach used by us) is the computation of the matrix productA 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, oncecomputed, 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 costis insignificant.
The computation of A A can be greatly accelerated by the useof 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 thecomputation 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 hwith very minimal degradation of image quality. However,the questions of computational efficiency are largely outsideof 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 bediagonalized; 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 þ λ2 IÞx ¼ A b directly by the conjugate-gradient descent method.
An additional feature of the algebraic method described hereis that it does not require that the set of detectors used be
Journal of Biomedical Optics

independent of the position of the source. We have taken advantage of this feature and have excluded the data points that arevery far off-axis. Specifically, for each source, we have used
only such detectors that are situated no further from the axisof 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 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
I 0 ð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 thenumerical value of D0 is not needed for the reconstructions
as long as we are interested in the dimensionless relative contrastxðrÞ ¼ δαðrÞ∕α0 . An additional parameter that enters theexpression for G0 ðr; r 0 Þ is the extrapolation distance l associated with the diffuse-nondiffuse boundary condition. The latter
was found to be equal to 0.83 mm.

4Data Restriction

The rejection of the datapoints deemed unreliable or too noisy isan established practice in DOT.51–54 Here this data restriction
methodology is used in a somewhat different context andwith a somewhat different goal compared to previous work.
In the usual approach, data points are rejected based on an apriori 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 theirlocation. Our purpose is not to suppress noise, but rather toinvestigate the reconstruction effects of the imaging window
restriction and the proximity of the chest wall phantom to thetarget. In particular, we are posing the following question:
how close the sources and detectors can be placed to a chestwall 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 topof 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 thelower 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.
The maximum available number of independent sourcedetector 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 thealgebraic image reconstruction), other data points were
eliminated by sampling of the detectors (in the algebraic reconstructions, every second detector was used and), and yet otherdata 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

026016-4

Downloaded From: http://biomedicaloptics.spiedigitallibrary.org/ on 02/11/2013 Terms of Use: http://spiedl.org/terms

February 2013 • Vol. 18(2)Ban et al.: Diffuse optical tomography in the presence of a chest wallFig. 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 theinner light blue square indicates the reconstruction region. The red dots indicate the source positions. A sample reconstruction is superimposed with thedrawing to illustrate the target shape and position. The red lines indicate the three lowest position of the chest wall phantom (other positions are outsideof 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.

measurements; see more detailed information for various data
restrictions in Table 1 below.
We now discuss the data restriction. The latter was accomplished by removing all sources and detectors situated above oneof 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 zerodata 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 thereconstruction. Thus, if NR ¼ 5, the data are restricted tosources 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 lineof 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 isxð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
Table 1 Data restriction sizes.

Numerical data
restriction (NR)

Number ofsourcesNumber of distinct
source-detector pairs

No restriction

35 × 35 ¼ 1225

20591492

535 × 31 ¼ 1085

164791528

35 × 28 ¼ 980

14661145

935 × 27 ¼ 945

14074728

10

35 × 26 ¼ 910

13457507

Journal of Biomedical Optics

approximations. This can result in reconstructing unphysicalnegative values of absorption, which are shown in the figuresby 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 wasused by us to invert the matrix A A. However, we havefound that the areas of negative absorption appear mostly inthe case of the analytic (fast) image reconstruction method,which cannot incorporate the positivity constraint. On theother hand, the algebraic reconstructions have produced either
no areas of negative absorption, or artifacts so severe (e.g., whend ¼ 2 cm and no numerical data restriction) that the use of thepositivity constraint was not useful. In other words, we did notencounter 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 thissection.
Reconstructions of the central slice of the medium (3 cmfrom either of the slab surfaces) obtained with varying valuesof d and various numerical data restriction (NR) are shownin 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 andd ¼ 5 cm away from the target. Data restriction with NR ¼ 5
results in a reasonable, yet suboptimal, image quality whend ¼ 5 cm, but not when d ¼ 2 cm. To remove the artifacts associated with the chest wall completely, NR ¼ 10 is required.
However, NR ¼ 8, 9, 10 used in conjunction with the analytical reconstruction yields an additional image artifact, whichis unrelated to the chest wall phantom. To see that this is true,
consider the images for d ¼ 17 cm, which are not affected at allby the chest wall phantom, yet exhibit the additional artifact justmentioned. This artifact is shown as a black area where thereconstructed 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 theabsorption in underestimated.
The appearance of this artifact can be understood. As mentioned above, the data restriction used with the analytical

026016-5

Downloaded From: http://biomedicaloptics.spiedigitallibrary.org/ on 02/11/2013 Terms of Use: http://spiedl.org/terms

February 2013 • Vol. 18(2)Ban et al.: Diffuse optical tomography in the presence of a chest wallFig. 4 Images of the central slice obtained by analytical (a) and algebraic (b) reconstruction methods. Different columns show data obtained withthe 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.

Journal of Biomedical Optics

026016-6

Downloaded From: http://biomedicaloptics.spiedigitallibrary.org/ on 02/11/2013 Terms of Use: http://spiedl.org/terms

February 2013 • Vol. 18(2)Ban et al.: Diffuse optical tomography in the presence of a chest wallFig. 5 Slices through the medium drawn at different depths (from the plane of sources) as indicated. Analytical image reconstruction method withd ¼ 5 cm and d ¼ 2 cm.

Journal of Biomedical Optics

026016-7

Downloaded From: http://biomedicaloptics.spiedigitallibrary.org/ on 02/11/2013 Terms of Use: http://spiedl.org/terms

February 2013 • Vol. 18(2)Ban et al.: Diffuse optical tomography in the presence of a chest wallFig. 6 Same as in Fig. 5 but obtained by algebraic reconstruction.

Journal of Biomedical Optics

026016-8

Downloaded From: http://biomedicaloptics.spiedigitallibrary.org/ on 02/11/2013 Terms of Use: http://spiedl.org/terms

February 2013 • Vol. 18(2)Ban et al.: Diffuse optical tomography in the presence of a chest wallreconstruction amounts to assuming that the truncated data
points are zero. In other words, we assume that, in the presenceof the target, the truncated source-detector pairs would havemeasured the same intensity as in the homogeneous slab, sothat Iðrd ; rs Þ ¼ I 0 ð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 apurely absorbing target, however, the actual intensity Iðrd ; rs Þ issmaller than I 0 ð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 thereconstruction 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 havebeen discarded. Of course, this analysis applies to the casewhen 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 isclearly visible, and the image quality is about the same aswith 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 andFig. 6 displays analogous data obtained by the algebraic
reconstruction. In addition, we show in the right-most columnof images the reconstruction averaged over the depth of the sample (that is, over different slices). Note that all reconstructedslices were used for the purpose of averaging, not only thoseshown 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 thearithmetic average of the reconstructions in all 13 slices.
These averaged (“projection”) images correspond to theusual radiological projections obtained with a parallel beamof 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 atthe 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 typicalfor 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
Journal of Biomedical Optics

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 chestwall phantom is also clearly visible at the correct location.
This result is somewhat unexpected and can be useful in thesituations 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 −1dz, 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 isequal, 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 theactual 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 promisingconclusions.
First, we have found that, when absorption contrast isof interest, simple CW instrumentation with linearized inversion
can suffice. This finding was obtained in spite of the presence ofthe chest wall phantom in close proximity to the target, which
first renders the inverse problem nonlinear and, second, differsfrom the background Intralipid and the target not only in absorption but also in scattering properties. Generally, under theconditions 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, tooptimize these approaches. This may require more sophisticatedand/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 alinear image reconstruction. The paper shows that, by means ofproperly restricting the data points used in image reconstruction,
it is possible to resolve a small absorptive target in the vicinity ofa spatially and optically large inhomogeneity and that the qualityof 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

026016-9

Downloaded From: http://biomedicaloptics.spiedigitallibrary.org/ on 02/11/2013 Terms of Use: http://spiedl.org/terms

February 2013 • Vol. 18(2)Ban et al.: Diffuse optical tomography in the presence of a chest wall

plane-parallel sample (we refer to this quantity as the projection), is not as sensitive to systematic errors encountered inimage reconstruction as the individual slices drawn throughthe medium. Thus, under certain conditions, the nonlinearityof 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 thenear-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 presenceof 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 theexperimental 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, DMS1108969 and DMS-1115574.

References
1. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl.
15(2), R41–R93 (1999).
2. D. A. Boas et al., “Imaging the body with diffuse optical tomography,”IEEE Signal Proc. Mag. 18(6), 57–75 (2001).
3. S. R. Arridge and J. C. Schotland, “Optical tomography: forward andinverse problems,” Inverse Probl. 25(12), 123010 (2009).
4. S. B. Colak et al., “Clinical optical tomography and NIR spectroscopy
for breast cancer detection,” IEEE J. Sel. Top. Quantum Electron.
5(4), 1143–1158 (1999).
5. D. J. Hawrysz and E. M. Sevick-Muraca, “Developments towarddiagnostic breast cancer imaging using near-infrared optical measurements and fluorescent contrast agants,” Neoplasia 2(5), 388–417(2000).
6. T. McBride et 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).
7. J. P. Culver et al., “Three-dimensional diffuse optical tomography in theparallel plane transmission geometry: evaluation of a hybrid frequency
domain/continuous wave clinical system for breast imaging,” Med.
Phys. 30(2), 235–247 (2003).
8. X. Intes et al., “In vivo continuous-wave optical breast imaging
enhanced with indocyanine green,” Med. Phys. 30(6), 1039–1047(2003).
9. A. Li et al., “Tomographic optical breast imaging guided bythree-dimensional mammography,” Appl. Opt. 42(25), 5181–5190(2003).
10. R. Choe et al., “Diffuse optical tomography of breast cancer during neoadjuvant chemotherapy: a case study with comparison to MRI,” Med.
Phys. 32(4), 1128–1139 (2005).
11. T. Yates et al., “Time-resolved optical mammography using a liquid
coupled interface,” J. Biomed. Opt. 10(5), 054011 (2005).
12. A. Corlu et al., “Three-dimensional in vivo fluorescence diffuse optical
tomography of breast cancer in humans,” Opt. Express 15(11), 6696–
6716 (2007).
13. A. Cerussi et al., “Predicting response to breast cancer neoadjuvant
chemotherapy using diffuse optical spectroscopy,” Proc. Natl. Acad.
Sci. U. S. A. 104(10), 4014–4019 (2007).

Journal of Biomedical Optics

14. R. Choe et 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).
15. T. Durduran et al., “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. 73(7), 076701 (2010).
16. B. Pogue et al., “Initial assessment of a simple system for frequency
domain diffuse optical tomography,” Phys. Med. Biol. 40(10),1709–1729 (1995).
17. D. Grosenick et al., “Time-domain scanning optical mammography: I.
recording and assessment of mammograms of 154 patients,” Phys. Med.
Biol. 50(11), 2429–2449 (2005).
18. A. Pifferi et al., “Four-wavelength time-resolved optical mammography
in the 680 980-nm range,” Opt. Lett. 28(13), 1138–1140 (2003).
19. J. Ge, B. Zhu, S. Regalado, and A. Godavarty, “Three-dimensional
fluorescence-enhanced optical tomography using a hand-held probe
based imaging system,” Med. Phys. 35(7), 3354–3363 (2008).
20. Z.-M. Wang et al., “Experimental demonstration of an analytic method
for image reconstruction in optical tomography with large data sets,”
Opt. Lett. 30(24), 3338–3340 (2005).
21. S. D. Konecky et al., “Imaging complex structures with diffuse light,”
Opt. Express 16(7), 5048–5060 (2008).
22. P. Bonfert-Taylor et al., “Information loss and reconstruction in diffuse
fluorescence tomography,” J. Opt. Soc. Am. A 29(3), 321–330 (2012).
23. S. Chaillat and G. Biros, “FaIMS: a fast algorithm for the inverse
medium problem with multiple frequencies and multiple sources forthe scalar Helmholtz equation,” J. Comput. Phys. 231(12), 4403–4421(2012).
24. R. Schulz, J. Ripoll, and V. Ntziachristos, “Noncontact optical tomography of turbid media,” Opt. Lett. 28(18), 1701–1703 (2003).
25. J. Ripoll, R. B. Schulz, and V. Ntziachristos, “Free-space propagation ofdiffuse light: theory and experiments,” Phys. Rev. Lett. 91(10), 103901(2003).
26. J. Ripoll and V. Ntziachristos, “Imaging scattering media from a distance: theory and applications of noncontact optical tomography,”Mod. Phys. Lett. 18(28–29), 1403–1431 (2004).
27. G. Turner et al., “Complete-angle projection diffuse optical tomography
by use of early photons,” Opt. Lett. 30(4), 409–411 (2005).
28. V. A. Markel and J. C. Schotland, “Inverse scattering for the diffusion
equation with general boundary conditions,” Phys. Rev. E 64(3),R035601 (2001).
29. V. A. Markel and J. C. Schotland, “Effects of sampling and limited
data in optical tomography,” Appl. Phys. Lett. 81(7), 1180–1182 (2002).
30. V. A. Markel, V. Mital, and J. 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).
31. V. A. Markel, J. A. O’Sullivan, and J. C. Schotland, “Inverse problem in
optical diffusion tomography. IV. Nonlinear inversion formulas,” J. Opt.
Soc. Am. A 20(5), 903–912 (2003).
32. V. A. Markel and J. C. Schotland, “Symmetries, inversion formulas, andimage reconstruction for optical tomography,” Phys. Rev. E 70(5),056616 (2004).
33. B. W. Pogue and M. S. Patterson, “Review of tissue simulating phantoms for optical spectroscopy, imaging and dosimetry,” J. Biomed. Opt.
11(4), 041102 (2006).
34. A. E. Cerussi et al., “Tissue phantoms in multicenter clinical trials for
diffuse optical technologies,” Biomed. Opt. Express 3(5), 966–971(2012).
35. Y. Ardeshirpour and Q. Zhu, “Optical tomography method that accounts
for tilted chest wall in breast imaging,” J. Biomed. Opt. 15(4), 041515(2010).
36. A. Kienle and T. 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).
37. P. Taroni et al., “In vivo absorption and scattering spectroscopy of biological tissues,” Photochem. Photobiol. Sci. 2(2), 124–129 (2003).
38. M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue
optical properties,” Appl. Opt. 28(12), 2331–2336 (1989).
39. D. A. Benaron and D. K. Stevenson, “Optical time-of-flight and absorbance imaging of biologic media,” Science 259(5100), 1463–1466 (1993).
40. S. Andersson-Engels et al., “Time-resolved transillumination for medical diagnostics,” Opt. Lett. 15(21), 1179–1181 (1990).

026016-10

Downloaded From: http://biomedicaloptics.spiedigitallibrary.org/ on 02/11/2013 Terms of Use: http://spiedl.org/terms

February 2013 • Vol. 18(2)Ban et al.: Diffuse optical tomography in the presence of a chest wall
41. S. L. Jacques, “Time-resolved reflectance spectroscopy in turbid tissues,” IEEE Trans. Biomed. Eng. 36(12), 1155–1161 (1989).
42. F. E. Schmidt et al., “Multiple-slice imaging of a tissue-equivalent phantom by use of time-resolved optical tomography,” Appl. Opt. 39(19),3380–3387 (2000).
43. V. Ntziachristos, X. H. Ma, and B. Chance, “Time-correlated single
photon counting imager for simultaneous magnetic resonance and nearinfrared mammography,” Rev. Sci. Instrum. 69(12), 4221–4233 (1998).
44. E. Gratton et 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, Hamamatsu, Japan, pp. 183–189 (1990).
45. J. B. Fishkin and E. Gratton, “Propagation of photon-density waves instrongly scattering media containing an absorbing semiinfinite planebounded by a straight edge,” J. Opt. Soc. Am. A 10(1), 127–140 (1993).
46. B. Chance et al., “Phase measurement of light absorption and scatteringin human tissues,” Rev. Sci. Instrum. 69(10), 3457–3481 (1998).
47. B. W. Pogue and M. S. Patterson, “Frequency-domain optical absorption spectroscopy of finite tissue volumes using diffusion theory,” Phys.
Med. Biol. 39(7), 1157–1180 (1994).

Journal of Biomedical Optics

48. J. C. Schotland, “Continuous wave diffusion imaging,” J. Opt. Soc.
Am. A 14(1), 275–279 (1997).
49. C. P. Gonatas et al., “Optical diffusion imaging using a direct inversion
method,” Phys. Rev. E 52(4), 4361–4365 (1995).
50. V. A. Markel, Z.-M. Wang, and J. C. Schotland, “Optical diffusion tomography with large data sets,” Proc. SPIE 5969, 59691B(2005).
51. A. Blasi et al., “Investigation of depth dependent changes in cerebral
haemodynamics during face preception in infants,” Phys. Med. Biol.
52(23), 6849–6864 (2007).
52. M. A. Franceschini et al., “Assesment of infant brain development withfrequency-domain near-infrared spectroscopy,” Proc. Natl. Acad. Sci.
U. S. A. 94(12), 6468–6473 (1997).
53. N. Roche-Labarbe et al., “Noninvasive optical measures of CBV, StO2,
C.BF index, and rCMRO2 in human premature neonates' brains in thefirst six weeks of life,” Hum. Brain Map 31(3), 341–352 (2010).
54. F. Orihuela-Espina et al., “Quality control and assurance in functional
near infrared spectroscopy (fNIRS) experimentation,” Phys. Med. Biol.
55(13), 3701–3724 (2010).

026016-11

Downloaded From: http://biomedicaloptics.spiedigitallibrary.org/ on 02/11/2013 Terms of Use: http://spiedl.org/terms

February 2013 • Vol. 18(2)