OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 16, Iss. 7 — Mar. 31, 2008
  • pp: 5048–5060
« Show journal navigation

Imaging complex structures with diffuse light

Soren D. Konecky, George Y. Panasyuk, Kijoon Lee, Vadim Markel, Arjun G. Yodh, and John C. Schotland  »View Author Affiliations


Optics Express, Vol. 16, Issue 7, pp. 5048-5060 (2008)
http://dx.doi.org/10.1364/OE.16.005048


View Full Text Article

Acrobat PDF (387 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

We use diffuse optical tomography to quantitatively reconstruct images of complex phantoms with millimeter sized features located centimeters deep within a highly-scattering medium. A non-contact instrument was employed to collect large data sets consisting of greater than 107 source-detector pairs. Images were reconstructed using a fast image reconstruction algorithm based on an analytic solution to the inverse scattering problem for diffuse light.

© 2008 Optical Society of America

1. Introduction

The study of light propagation in highly scattering random media such as clouds, paint and tissue is of fundamental interest and considerable applied importance [1

1. M. C.W. van Rossum and T. M. Nieuwenhuizen, “Multiple scattering of classical waves: microscopy, mesoscopy and diffusion,” Rev. Mod. Phys. 71, 313–371 (1999). [CrossRef]

]. Novel physical effects such as the dynamic correlations of speckle patterns and weak localization arise from coherent effects in multiple scattering. Incoherent effects are also of interest and have lead to new experimental tools to probe the structure of highly-scattering media. One such technique, known as diffuse optical tomography (DOT), uses highly scattered near-infrared light to image biological tissue and provide functional information about physiological parameters such as blood volume and oxygenation [2–4

2. A. Yodh and B. Chance, “Spectroscopy and imaging with diffuse light,” Phys. Today pp. 34–40 (1995). [CrossRef]

]. Current clinical applications of DOT include breast imaging [5–12

5. S. B. Colak, M. B. van der Mark, G. W. Hooft, J. H. Hoogenraad, E. S. van der Linden, and F. A. Kuijpers, “Clinical optical tomography and NIR spectroscopy for breast cancer detection,” IEEE J. Sel. Top. Quantum Electron. 5, 1143–1158 (1999). [CrossRef]

] and functional brain mapping [13–21

13. D. A. Benaron, J. van Houten, D. C. Ho, S. D. Spilman, and D. K. Stevenson, “Imaging neonatal brain injury using light-based optical tomography,” Pediatr. Res. 35, A378–A378 (1994).

].

In an optically-thick medium such as the human breast, multiple scattering of light creates a fundamental obstruction to the direct formation of images. DOT overcomes this problem to some extent by solving an appropriate inverse problem, usually based on the diffusion equation, wherein the optical properties of a highly-scattering medium are reconstructed from boundary measurements. Since such inverse scattering problems are severely ill-posed [3

3. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, R41–R93 (1999). [CrossRef]

], the resultant image quality in DOT is expected to be poor; typical images resemble structureless ‘blobs’. As a result many accept that anatomically accurate DOT images cannot be reconstructed. Accordingly, the emphasis in DOT has been on functional imaging, and on multi-modality imaging, in which simultaneously acquired MRI or CT images are used to provide anatomical detail [4

4. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005). [CrossRef] [PubMed]

].

It has been suggested that the relatively low quality of images in DOT can be improved by the use of large data sets [22–26

22. J. C. Schotland, “Continuous wave diffusion imaging,” J. Opt. Soc. Am. A 14, 275–279 (1997). [CrossRef]

]. Moreover, data sets of approximately 108 measurements may be readily acquired with non-contact DOT systems [27–29

27. G. Turner, G. Zacharakis, A. Soubret, J. Ripoll, and V. Ntziachristos, “Complete-angle projection diffuse optical tomography by use of early photons,” Opt. Lett. 30, 409–411 (2005). [CrossRef] [PubMed]

]. For example, Turner et al. have shown that simple shapes can be imaged in an experiment using an optically-thin sample by employing time gating of early-arriving photons [27

27. G. Turner, G. Zacharakis, A. Soubret, J. Ripoll, and V. Ntziachristos, “Complete-angle projection diffuse optical tomography by use of early photons,” Opt. Lett. 30, 409–411 (2005). [CrossRef] [PubMed]

]. However, this approach uses photons that travel ballistically through the sample and is therefore not useful in optically thick tissues such as breast and brain.

In this contribution, we demonstrate for the first time that three-dimensional, quantitative images of the interior of such an optically thick medium, with millimeter scale complex structure, can be acquired using DOT. The reconstructed images exhibit spatially-resolved features that are much smaller than the overall size of the sample and are obtained by utilizing large data sets and fast image reconstruction algorithms developed in earlier work [23–26

23. V. A. Markel and J. C. Schotland, “Inverse scattering for the diffusion equation with general boundary conditions,” Phys. Rev. E 64, R035,601 (2001). [CrossRef]

, 30

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, 890–902 (2003). [CrossRef]

, 31

31. V. A. Markel and J. C. Schotland, “Symmetries, inversion formulas and image reconstruction for optical tomography,” Phys. Rev. E 70, 056,616 (2004). [CrossRef]

]. Our findings have important implications for biomedical applications, since we achieve image quality and spatial resolution sufficient for visualization of anatomical detail.

2. Methods

2.1. Instrumentation

The experimental apparatus consisted of a continuous-wave diode laser (Model TC40, SDL Inc., San Jose, CA) operating at 785 nm and coupled via a 100 µm multi-mode fiber to a pair of galvanometer-controlled mirrors (Innovations in Optics, Woburn, MA). The mirrors raster scanned the beam (on a 35×35 rectangular grid with 4 mm step size) over a 13.6×13.6 cm2 square on one side of an imaging tank containing a target. For each laser beam position, the light transmitted through the tank was collected by a f=25 mm F/0.95 lens and focused on a front-illuminated thermoelectric-cooled 16-bit CCD array (DV887ECS-UV, Andor Technology, Belfast, Ireland). A 20×20 cm2 square area of the opposite surface was mapped onto a square grid of 512×512 CCD pixels; this corresponds to a grid of detectors with 0.4 mm step. The signal recorded by each CCD pixel for a given source beam position defined an independent measurement. The total size of the data set recorded in a single experiment was (35×512)2≈3×108. Somewhat smaller subsets of the data were used for image reconstruction.

The imaging tank had an inner thickness L=6 cm, and was filled with a mixture of water, a highly scattering fat emulsion (Liposyn III, 30%, Abbott Laboratories, Chicago, IL), and India ink (Black India 4415, Sanford, Bellwood, IL). The absorption and the reduced scattering coefficients of the mixture at 786 nm were µa=0.05 cm-1 and µs=7.5 cm-1, respectively. The transport mean free path was *=1/(µa+µs) ≈ 1.3 mm. The diffuse wave number (i.e. the inverse characteristic length over which the diffuse waves decay exponentially) was kd3μa*1.1cm1 . The tank thickness was therefore much larger than *, and thus the experiments were carried out in the diffusion regime. In addition kdL≈6.6, and thus the transmitted light was substantially attenuated. These parameters are typical of human breast tissues in the NIR spectral region.

2.2. Targets

Targets were made of a mixture of silicone rubber (RTV-12, General Electric, Waterford, NY), titanium dioxide (T-8141,Sigma-Aldrich, St. Louis, MO), and carbon black (Raven 5000 Ultra Powder II, Columbian Chemicals Co., Marietta, GA) [32

32. R. Choe, “Diffuse Optical Tomography and Spectroscopy of Breast Cancer and Fetal Brain,” Ph.D. thesis, University of Pennsylvania (2005).

]. The ingredients were mixed in a proportion such that the reduced scattering coefficient of the targets was approximately matched to that in the surrounding fluid while the absorption coefficient µa was four times larger. Targets were made in the shape of letters (3 cm tall, 2 cm wide, 5 mm thick with individual components 3 mm in width), and bars (6 cm tall, 5 mm thick, with widths of 7–9 mm).

Titration experiments were performed using an 18 mm tall clear plastic cylinder with a diameter of 17 mm and 1 mm thick walls. It was positioned in the center of the tank. Laboratory tubing (outer diameter 2 mm) was used to flow fluid through the cylinder. During the first scan, the cylinder contained matching fluid identical to the fluid in the tank. Twelve titrations were then performed in which the ratio of ink concentration in the cylinder to that in the tank gave an expected absorption contrast ranging from 2:1 to 64:1.

2.3. Theoretical background

We use the diffusion approximation to the radiative transport equation [1

1. M. C.W. van Rossum and T. M. Nieuwenhuizen, “Multiple scattering of classical waves: microscopy, mesoscopy and diffusion,” Rev. Mod. Phys. 71, 313–371 (1999). [CrossRef]

, 33

33. A. Ishimaru, Wave Propagation and Scattering in Random Media (IEEE Press, Piscataway, NJ, 1997).

] to model the propagation of multiply-scattered light. The approximation is valid if µaµs. This condition was satisfied everywhere in the experimental medium. In the diffusion approximation, the electromagnetic energy density u(r) inside the medium obeys the diffusion equation

[·D(r)+cμa(r)]u(r)=S(r).
(1)

Here D=cℓ*/3=c/3(µa+µs) the diffusion coefficient, c is the average speed of light in the medium and S is the power density of the source. The diffusion equation (1) is supplemented by the boundary condition

u+n̂·u=0,
(2)

where is the extrapolation distance [33

33. A. Ishimaru, Wave Propagation and Scattering in Random Media (IEEE Press, Piscataway, NJ, 1997).

, 34

34. R. Aronson, “Boundary conditions for diffuse light,” J. Opt. Soc. Am. A 12, 2532–2539 (1995). [CrossRef]

]. In the same limit that is required for validity of the diffusion approximation, Dc/3µs. Thus, if the reduced scattering coefficient µs is the same inside the target and the surrounding medium, we have DD (0)=const. In this case, the goal of DOT is to reconstruct µa(r) from a set of boundary measurements of u(r), assuming D (0) is known. We, however, do not imply that reconstruction of µs (r) is not possible or not important in general. As was discussed in references [30

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, 890–902 (2003). [CrossRef]

, 31

31. V. A. Markel and J. C. Schotland, “Symmetries, inversion formulas and image reconstruction for optical tomography,” Phys. Rev. E 70, 056,616 (2004). [CrossRef]

], simultaneous reconstruction of µa and µs is possible with the use of either time-resolved or frequency resolved measurements. The fast image reconstruction method described in Section 2.4 can easily be generalized to the case of heterogeneous µs [30

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, 890–902 (2003). [CrossRef]

, 31

31. V. A. Markel and J. C. Schotland, “Symmetries, inversion formulas and image reconstruction for optical tomography,” Phys. Rev. E 70, 056,616 (2004). [CrossRef]

].

[D(0)2+cμa(r)]G(r,r)=δ(rr)
(3)

G(rd,rs)=G0(rd,rs)VG0(rd,r)cδμa(r)G(r,rs)d3r,
(4)

where the integral is taken over the volume of the sample and G0 is the Green’s function in the reference medium with µa=µ (0) a. G 0 satisfies Eq. (3) with the substitution µa(r)→µ (0) a and can be found analytically in an infinite slab [30

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, 890–902 (2003). [CrossRef]

]. The left-hand side of Eq. (4) is directly measurable but G(r,r s), which appears in the right-hand side, can not be measured since the point r lies inside the medium. Since G depends on δµa, the right-hand side of (4) is, generally, a nonlinear functional of the latter. Several approximate linearization methods can be used to transform Eq. (4) to an equation which is linear in δµa. We have used the first Rytov approximation, according to which, the right-hand side in (4) is replaced by the expression

G0(rd,rs)exp[VG0(rd,r)cδμa(r)G0(r,rs)d3rG0(rd,rs)]. .

Note that this expression coincides with the right-hand side of (4) to leading order in δµa. We then define the data function as

ϕ(rd,rs)=G0(rd,rs)1n[T(rd,rs)],
(5)

where T(r d,r s)≡I(r d,r s)/I 0(r d,r s) is the experimentally measurable transmission function. Here I 0(r d,r s) is the result of a measurement in a reference homogeneous medium with δµa=0. We then obtain the following linear integral equation:

ϕ(rd,rs)=VG0(rd,r)cδμa(r)G0(r,rs)d3r.
(6)

The kernel Γ(r d,r s;r)=G 0(r d,r)G 0(r,r s) is known analytically and the right-hand side ϕ (r d,r s) is measured experimentally. Note that we have used the fact that I(r d,r s)/I 0(r d,r s)=G(r d,rs)/G 0(r d,r s) to eliminate the coupling constants Cd and Cs. Equations (5) and (6) are the main equations of linearized DOT. The method we have used for inverting (6) is briefly described below.

2.4. Inversion formula

The image reconstruction method used by us is based on the translational symmetry of the slab geometry. The latter allows one to write G 0 as

G0(ρ,z;ρ,z)=g(q;z,z)exp[iq·(ρρ)]d2q(2π)2,
(7)

where the explicit expressions for g(q; z, z′) are given in Ref. [30

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, 890–902 (2003). [CrossRef]

]. Note, both g(q; z, z′) and G 0(ρ, z;ρ′, z′) are real for continuous-wave measurements. We substitute this plane-wave decomposition of G 0 into (6) and take a four-dimensional Fourier transform with respect to the source and detector transverse coordinates to obtain

ϕ(qd,qs)=czdzsg(qs;zsz)g(qd;z,zd)δμ~a(qs+qd,z)dz,
(8)

where the Fourier-transforms are defined as

ϕ~(qd,qs)=ρd,ρsϕ(ρd,zd;ρs,zs)exp[i(ρd·qd+ρs·qs)],
(9)
δμ~a(q,z)=δμa(ρ,z)exp(iρ·q)d2ρ.
(10)

With the change of variables q d=q+p, q s=-p, we can now re-write (8) as

ψ(q,p)=zdzsκ(q,p;z)δμ~a(q,z)dz,
(11)

where ψ(q,p)=ϕ˜(q+p,-p) and κ(q,p; z)=cg(q+p; zs, z)g(-p; z, zd).

Thus, the problem has been reduced to inverting a set of linear one-dimensional integral equations with known kernels κ(q,p; z). We do so for a discrete set of samples of the Fourier variable q and then use the inverse Fourier transform to obtain δµa(ρ, z). Note that Eq. (11) is inverted for each sample of q by utilizing multiple discrete values of p. To this end, we construct analytically the pseudoinverse to the linear integral operator of the left-hand side of (11). The mathematical details of this construction are given in Ref. [30

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, 890–902 (2003). [CrossRef]

]. Computing δµa(ρ, z) by the inverse Fourier transform accounting for sampling and truncation of real-space data is discussed in Refs. [25

25. V. A. Markel and J. C. Schotland, “Effects of sampling and limited data in optical tomography,” Appl. Phys. Lett. 81, 1180–1182 (2002). [CrossRef]

,31

31. V. A. Markel and J. C. Schotland, “Symmetries, inversion formulas and image reconstruction for optical tomography,” Phys. Rev. E 70, 056,616 (2004). [CrossRef]

]. In Ref. [31

31. V. A. Markel and J. C. Schotland, “Symmetries, inversion formulas and image reconstruction for optical tomography,” Phys. Rev. E 70, 056,616 (2004). [CrossRef]

], we have considered symmetry-based image reconstruction algorithms from a more general point of view and a number of special cases have been derived, including the specific modality used in this paper (Section 3.C.1, “Fourier tomography”).

The computational complexity of the image reconstruction algorithm described in this section is O(NqN 3 p), where Nq is the number of discrete values of the vector q, while Np is the number of discrete values of p used to invert each one-dimensional equation (11). This is the computational cost of inverting the integral operator only. To this one should add the cost of Fourier-transforming the real-space data function and computing ψ(q,p). With the use of the fast Fourier transform, the latter scales as O(N 2 qlnNq). If the same amount of data and the same grid of voxels is used with a purely algebraic image reconstruction method, the computational cost of matrix inversion scales as O((NqNp)3). In the above estimate, we have assumed that the number of measurements is equal to the number of voxels. We thus see that the fast image reconstruction methods exploit the block structure of the linear operator that couples the data to the unknown function. Instead of inverting a large matrix of size NqNp, these methods require inversion of Nq matrices of the size Np each, thus gaining a factor of N 2 p in computation time. We note that reconstructions which utilized data sets of up to 107 source-detector pairs required less than one minute of CPU time on a 1.3GHz workstation.

2.5. Simulated data

Simulated data functions for a point absorber were generated within the linear approximation by representing δµa(r) as a delta-function. Data for infinite lattices of sources and detectors was computed directly in the Fourier domain using (8). For finite lattices we computed G 0 numerically using (7) and substituted the result into (9); in the latter formula, summation over real space variables was truncated.

We then added noise to the data function generated for the finite lattice as follows.We expect the standard deviation due to shot noise for repeated measurements of a single source detector pair to be equal to the square root of the mean number of detected photo-electrons ( I(rs,rd)) ). Propogating the error to the data function defined in equation (5) we get to first order σdata(rs,rd)=2I0(rs,rd)ϕ(rs,rd) . We scaled both the simulated data function and the simulated reference intensity I 0 to have the same maximum values as were experimentally measured. The resulting data function with shot noise was then

ϕnoise(rs,rd)=ϕ(rs,rd)+σdata(rs,rd)R.
(12)

Here R is a random variable with a Gaussian distribution and a variance of one. We note that since I 0ϕ, the shot noise depends primarily on the reference intensity, not on the data function itself. In order to simulate the effect of read noise and dark current, we determined the stardard deviation σbg of a background image taken with the light source blocked. This error was propagated as σdata=2σbg2 , and the resulting data function was calculated according to (12).

3. Experimental results

For the reconstructions shown in Fig. 1, the target was constructed of silicone rubber and shaped in the form of the letters “DOT” and “PENN”. In the first experiment, we placed the letters “DOT” one centimeter from the source plane and the letters “PENN” one centimeter from the detector plane, directly behind the “DOT” letters. The reconstruction is shown in Fig. 1(b). The letters are clearly visible. Note that the central slice from the middle of the tank is empty, as expected. In a second experiment we placed the letters “DOT” in the center of the tank (the letters “PENN” were not present) three centimeters from source and detector planes. The reconstruction is shown in Fig. 1(c). The letters are clearly reconstructed.

In Fig. 2 we show the data used to reconstruct the images. Only data corresponding to a single source beam position is shown. For each reconstruction we show, from left to right, the reference intensity I 0 when the target is not present (scattering fluid only), the intensity I when the target is present, and the Rytov data, ϕ=-log(I/I 0), which is used in the reconstruction algorithm (see Methods). Note that the letters cannot be identified by simply inspecting images of the transmitted light. Structure is visible in the Rytov data when the letters “PENN” are close to the detector plane (Fig. 2(a)). However, when the letters “DOT” are in the center of the tank (Fig. 2(b)) their shape is completely blurred.

In order to quantify the transverse resolution of the reconstructed images, we prepared several bar targets from the same material described above. The bars were 7 mm to 9 mm thick and placed consecutively (one at a time) in the center of the slab. Fig. 3(a) shows the corresponding reconstructions. As the bar widths decrease, the modulation depth between bars decreases. As can be seen, all but the 7 mm bar target are well resolved. Fig. 3(b) shows reconstructions for two experiments in which the 7mm bar target was positioned 1 cm from the source and detector planes, respectively. The bars in this figure are well resolved and the images are smoother and have fewer artifacts. As can be expected, the image is better resolved when the target is closer to the detector plane. This is because the detectors are sampled on a finer grid than the sources. Experiments were also performed with bar targets having an absorption contrast of 2:1. The resulting images were very similar to those acquired with a 4:1 absorption contrast bar targets, but they contained less contrast and more noise.

Fig. 1. Slices from three dimensional image reconstructions of the relative absorption coefficient (δµa/µ (0) a) for targets suspended in a 6 cmthick slab filled with highly scattering fluid. The three slices shown for each reconstruction correspond to depths of 1 cm (left), 3 cm (middle), and 5 cm (right) from the source plane. The field of view in each slice is 16 cm×16 cm. The quantity plotted is δµa/µ (0) a (a) Schematics of the positions of the letters during the experiments. Left: The target consists of letters “DOT” and “PENN”, suspended 1 cm and 5 cm from the source plane, respectively. Right: The target consists only of the letters “DOT” suspended 3 cm from the source plane, i.e., in the center of the slab. (b) Reconstructed image of the letters “DOT” and “PENN” (c) Reconstructed image of the letters “DOT”.

The images in Fig. 3 were reconstructed using approximately 107 measurements. The effect of changing the size of the data set was investigated by sampling the detectors on a grid with a step size of 2 mm, five times larger than the minimum experimentally available detector spacing. We found that increasing the number of detectors up to the experimentally available maximum did not improve image quality. That is, reconstructions performed using 2mmsource separations or with a denser sampling of CCD pixels were visually indistinguishable from those in Fig. 3. However, decreasing the number of data points does result in poorer image quality, as is illustrated in Fig. 4. From left to right, three separate reconstructions of the 8 mm bar target are shown. The data for these reconstructions were taken from a single experiment with the target positioned in the center of the tank. All reconstruction parameters are kept constant, except for the number of measurements used. The reconstruction on the left uses 8×106 measurements with sources and detectors sampled with 4 mm and 2 mm steps, respectively. In the center reconstruction, we use 4 mm spacing for both the sources and detectors, which corresponds to 2×106 data points. In the reconstruction on the right, the sampling is 8 mm for sources and 4 mm for detectors, or, approximately, 5×105 data points. It can be seen that as the number of measurements decreases, image artifacts become more prominent and resolution is lost. We note that the optimal number measurements for a given experiment will vary depending on factors such as experimental geometry and noise level. For example, smaller/larger data sets would be optimal if the reconstruction field of view was smaller/larger.

Fig. 2. Representative CCD data for the image reconstructions shown in Fig. 1. Each image corresponds to the measured light intensity for a single source beam position. The left column shows the reference intensity I 0 when the target is not present (scattering fluid only). The middle column shows the intensity I when the target is present. The right column shows the Rytov data, ϕ=-log(I/I 0), which is used in the reconstruction algorithm. (a) The target consists of the letters “DOT” and “PENN” (b) The target consists of the letters “DOT” only
Fig. 3. Reconstructed images of bar targets. Only slices drawn at the depth of the actual target are shown. (a) 7 mm to 9 mm bar targets located in the center of the tank. Here d denotes the width of the individual bars in the targets. (b) The 7 mm bar target located 1 cm from the source (left) and detector (right) planes.
Fig. 4. Images of the 8 mm bar target from a single experiment. All reconstruction parameters are held fixed except for the number of measurements used. From left to right, correspond to N=8×106, 2×106, and 5×105 measurements were used for the reconstruction.
Fig. 5. Reconstructed contrast of the absorption coefficient µa/µ (0) a between the cylinder and the tank vs. the expected contrast.

To investigate the capability for quantitative reconstruction of the absorption coefficient, we have performed a titration experiment. A clear plastic cylinder was positioned in the center of the tank, and laboratory tubing was used to flow fluid through the cylinder. During the first scan, the cylinder contained scattering fluid identical to the fluid in the tank. Twelve titrations were then performed in which the ratio of ink concentration in the cylinder, to that in the tank, gave an expected absorption contrast ranging from 2:1 to 64:1. The absorption contrast, i.e., the ratio of the absorption coefficient in the cylinder to that in the surrounding fluid, was taken from the corresponding reconstructed value for a single voxel located inside the cylinder. This voxel was chosen to be the voxel with the maximum reconstructed value of the absorption coefficient for the tenth titration. In Fig. 5, the reconstructed contrast µa/µ (0) a is plotted against the expected contrast. It can be seen that the absorption is quantitatively reconstructed with a linear dependence on ink concentration over nearly a decade in absorption contrast. Deviation from linearity occurs at higher concentrations, as expected.

4. Discussion

4.1. Transverse resolution and Fourier analysis of data

Factors that limit the ability to achieve even greater spatial resolution include the finite size of the source and detector grids, the shot noise in the measured light, the noise in the CCD array (a combination of dark current and read noise), and systematic errors. The latter are errors of the model. In particular, the light transmitted through the tank is described by the diffusion equation only approximately. There are also systematic errors associated with non-ideal optics (diffraction, multiple reflections on the lens and tank surfaces). Finally, there are systematic errors associated with the nonlinearity of the inverse problem of DOT.

To better understand the dependence of transverse spatial resolution on the factors listed above, we performed the following analysis. Let all inhomogeneities be confined in a thin slice of the medium parallel to the slab, z 0-h/2<z<z 0+h/2, where h is small. We seek to reconstruct the deviation of the absorption coefficient in this slice, δµa(ρ, z 0), from the respective value in the surrounding fluid, as a function of the transverse variable ρ=(x,y). The formula (11) derived in the Methods Section becomes in this case

ψ(q,p)=hκ(q,p)δμ˜a(q,z0).

Here ψ(q,p) is the measurable data function, κ(q,p) is known analytically (see Methods for precise definitions of ψ and κ), and δµa(q, z 0) is the transverse Fourier transform of δµa(ρ, z 0) with respect to ρ. The real-space function δµa(ρ, z 0) is obtained by inverse Fourier transformation, namely, δμa(ρ,z0)=1hexp(iq·ρ)ψ(q,p)κ(q,p)d2q(2π)2. The choice of p in the above formula is arbitrary since the problem is overdetermined (we use four-dimensional data to reconstruct a two-dimensional function); it is sufficient to choose p=0. We then use the data function ψ(q,0) for reconstruction. The latter is taken from experiment and contains noise. It is important to note that the “ideal” data function is rapidly (exponentially) decreasing with |q|. However, the experimental noise is approximately white, i.e., the amplitude of its Fourier transform is approximately constant in the range of |q| which is of interest. In image reconstruction, integration in the above formula is over a disc |q| < q max such that outside of this disc, the signal-to-noise ratio in ψ(q,0) becomes smaller than unity. It then follows that the minimum spatial feature that can be resolved has a characteristic transverse dimension of Δx~π/q max.

A few comments on the above analysis are necessary. First, if the target is not confined to a thin slice, the transverse resolution can only be lower than the estimate π/q max. Second, one can attempt to use additional degrees of freedom (i.e., N distinct values of p) to improve the resolution. This approach, however, is not expected to yield a significant improvement since the noise amplitude decreases no faster than 1/√N, while the function ψ(q,p) decreases exponentially with |q|.

To illustrate the above resolution estimate, we plot in Fig. 6 the power spectrum of the data for a point absorber located in the center of the tank. The data function ψ(q,0) is in this case cylindrically symmetric and we can write ψ(q,0)=Ψ(q), where q=|q|. The five curves shown in Fig. 6 contain (i) simulated ideal data corresponding to infinitely large and dense grids of noiseless sources and detectors, (ii) simulated data (without noise) for the finite grids of sources and detectors that were used in the experiments, (iii) simulated data for the finite grids with background noise added (i.e. Gaussian distributed noise with a variance equal to the variance of the detected signal when the laser is off), (iv) simulated data for the finite grid with shot noise added (i.e. Gaussian distributed noise with a variance equal to the number of photoelectrons detected experimentally), and (v) data from an experiment with a small (4 mm) highly absorbing target in the tank (see Methods for details of the simulations). The early onset of noise in the simulated data with shot noise, suggests the latter is the limiting factor in our experiments. We thus conclude that the other errors discussed above (e.g. the finite numbers of sources and detectors, CCD noise, errors in the diffusion model, non-ideal optics, and non-linearity) do not play a significant role in limiting image resolution. We see that noise begins to dominate the data at q max/π≈0.1 mm-1. This corresponds to a spatial resolution of ℓ1 cm. For objects closer to the surface, the power spectrum decays more slowly. For example, for a point absorber located 1 cm from the surface, shot noise begins to dominate the data at q max/π≈0.25 mm-1, which corresponds to a spatial resolution of ℓ4 mm.

Fig. 6. Power spectra of the data function |Ψ(q)|2 (defined in the text) for a small absorber located in the center of the tank. The different curves correspond to simulated ideal data, simulated data from finite grids of sources and detectors, simulated data with background noise, simulated data with shot noise, and to experimental data, as indicated.

4.2. Titration experiments

The saturation effect seen in Fig. 5 can be explained as follows. When the absorption of the target is very high, most of the light is absorbed as soon as it enters the target and almost no light reaches its interior; the absorption coefficient of the interior can be arbitrarily changed with no significant effect on the measured signal. This is a manifestation of the nonlinearity of the inverse problem of DOT. Note that the onset of nonlinearity in Fig. 5 occurs at absorber concentrations which are well beyond the range of contrast encountered in biological tissues. We finally note that the nonlinearity becomes stronger when the absorbing target becomes larger, as less and less light can penetrate the interior.

5. Conclusion

We have obtained high quality, quantitatively accurate, reconstructed images of complex structures deeply embedded in highly-scattering media. This result was achieved by using both a non-contact scanner capable of collecting large amounts of data and a fast image reconstruction algorithm capable of utilizing this data. We expect that these techniques will greatly improve the utility of diffuse optical methods in biomedical imaging.

Acknowledgments

This research was supported, in part, by the NIH grants R01-CA75124, R01-EB-002109, NTROI 1U54CA105480, P41RR02305, R21EB004524, R01-EB-004832 and by the NSF under grants DMS-0554100 and EEC-0615857.

References and links

1.

M. C.W. van Rossum and T. M. Nieuwenhuizen, “Multiple scattering of classical waves: microscopy, mesoscopy and diffusion,” Rev. Mod. Phys. 71, 313–371 (1999). [CrossRef]

2.

A. Yodh and B. Chance, “Spectroscopy and imaging with diffuse light,” Phys. Today pp. 34–40 (1995). [CrossRef]

3.

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, R41–R93 (1999). [CrossRef]

4.

A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005). [CrossRef] [PubMed]

5.

S. B. Colak, M. B. van der Mark, G. W. Hooft, J. H. Hoogenraad, E. S. van der Linden, and F. A. Kuijpers, “Clinical optical tomography and NIR spectroscopy for breast cancer detection,” IEEE J. Sel. Top. Quantum Electron. 5, 1143–1158 (1999). [CrossRef]

6.

D. J. Hawrysz and E. M. Sevick-Muraca, “Developments toward diagnostic breast cancer imaging using near-infrared optical measurements and fluorescent contrast agants,” Neoplasia 2, 388–417 (2000). [CrossRef]

7.

Y. Xu, X. J. Gu, L. L. Fajardo, and H. B. Jiang, “In vivo breast imaging with diffuse optical tomography based on higher-order diffusion equations,” Appl. Opt. 42, 3163–3169 (2003). [CrossRef] [PubMed]

8.

J. P. Culver, R. Choe, M. J. Holboke, L. Zubkov, T. Durduran, A. Slemp, V. Ntziachristos, B. Chance, and A. G. Yodh, “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, 235–247 (2003). [CrossRef] [PubMed]

9.

Q. I. Zhu, M. M. Huang, N. G. Chen, K. Zarfos, B. Jagjivan, M. Kane, P. Hedge, and S. H. Kurtzman, “Ultrasound-guided optical tomographic imaging of malignant and benign breast lesions: Initial clinical results of 19 cases,” Neoplasia 5, 379–388 (2003). [PubMed]

10.

X. Intes, J. Ripoll, Y. Chen, S. Nioka, A. G. Yodh, and B. Chance, “In vivo continuous-wave optical breast imaging enhanced with Indocyanine Green,” Med. Phys. 30, 1039–1047 (2003). [CrossRef] [PubMed]

11.

A. Li, E. L. Miller, M. E. Kilmer, T. J. Brukilacchio, T. Chaves, J. Stott, Q. Zhang, T. Wu, M. Chorlton, R. H. Moore, D. B. Kopans, and D. A. Boas, “Tomographic optical breast imaging guided by three-dimensional mammography,” Appl. Opt. 42, 5181–5190 (2003). [CrossRef] [PubMed]

12.

S. D. Jiang, B. W. Pogue, T. O. McBride, M. M. Doyley, S. P. Poplack, and K. D. Paulsen, “Near-infrared breast tomography calibration with optoelastic tissue simulating phantoms,” J. Electron. Imaging 12, 613–620 (2003). [CrossRef]

13.

D. A. Benaron, J. van Houten, D. C. Ho, S. D. Spilman, and D. K. Stevenson, “Imaging neonatal brain injury using light-based optical tomography,” Pediatr. Res. 35, A378–A378 (1994).

14.

J. P. van Houten, W. F. Cheong, E. L. Kermit, T. R. Machold, D. K. Stevenson, and D. A. Benaron, “Clinical measurement of brain oxygenation and function using light-based optical tomography,” Pediatr. Res. 39, 2273–2273 (1996). [CrossRef]

15.

A. Y. Bluestone, G. Abdoulaev, C. H. Schmitz, R. L. Barbour, and A. H. Hielscher, “Three-dimensional optical tomography of hemodynamics in the human head,” Opt. Express 9, 272–286 (2001). [CrossRef] [PubMed]

16.

G. Strangman, D. A. Boas, and J. P. Sutton, “Non-invasive neuroimaging using near-infrared light,” Biol. Psychiatry 52, 679–693 (2002). [CrossRef] [PubMed]

17.

J. C. Hebden, A. Gibson, R. M. Yusof, N. Everdell, E. M. C. Hillman, D. T. Delpy, S. R. Arridge, T. Austin, J. H. Meek, and J. S. Wyatt, “Three-dimensional optical tomography of the premature infant brain,” Phys. Med. Biol. 47, 4155–4166 (2002). [CrossRef] [PubMed]

18.

J. C. Hebden, “Advances in optical imaging of the newborn infant brain,” Psychophysiology 40, 501–510 (2003). [CrossRef] [PubMed]

19.

J. P. Culver, A. M. Siegel, J. J. Stott, and D. A. Boas, “Volumetric diffuse optical tomography of brain activity,” Opt. Lett. 28, 2061–2063 (2003). [CrossRef] [PubMed]

20.

A. M. Siegel, J. P. Culver, J. B. Mandeville, and D. A. Boas, “Temporal comparison of functional brain imaging with diffuse optical tomography and fMRI during rat forepaw stimulation,” Phys. Med. Biol. 48, 1391–1403 (2003). [CrossRef] [PubMed]

21.

G. Q. Yu, T. Durduran, D. Furuya, J. H. Greenberg, and A. G. Yodh, “Frequency-domain multiplexing system for in vivo diffuse light measurements of rapid cerebral hemodynamics,” Appl. Opt. 42, 2931–2939 (2003). [CrossRef] [PubMed]

22.

J. C. Schotland, “Continuous wave diffusion imaging,” J. Opt. Soc. Am. A 14, 275–279 (1997). [CrossRef]

23.

V. A. Markel and J. C. Schotland, “Inverse scattering for the diffusion equation with general boundary conditions,” Phys. Rev. E 64, R035,601 (2001). [CrossRef]

24.

V. A. Markel and J. C. Schotland, “The inverse problem in optical diffusion tomography. II. Inversion with boundary conditions,” J. Opt. Soc. Am. A 19, 558–566 (2002). [CrossRef]

25.

V. A. Markel and J. C. Schotland, “Effects of sampling and limited data in optical tomography,” Appl. Phys. Lett. 81, 1180–1182 (2002). [CrossRef]

26.

J. C. Schotland and V. A. Markel, “Inverse scattering with diffusing waves,” J. Opt. Soc. Am. A 18, 2767–2777 (2001). [CrossRef]

27.

G. Turner, G. Zacharakis, A. Soubret, J. Ripoll, and V. Ntziachristos, “Complete-angle projection diffuse optical tomography by use of early photons,” Opt. Lett. 30, 409–411 (2005). [CrossRef] [PubMed]

28.

R. Schulz, J. Ripoll, and V. Ntziachristos, “Noncontact optical tomography of turbid media,” Opt. Lett. 28, 1701–1703 (2003). [CrossRef] [PubMed]

29.

Z.-M. Wang, G. Y. Panasyuk, V. A. Markel, and J. C. Schotland, “Experimental demonstration of an analytic method for image reconstruction in optical tomography with large data sets,” Opt. Lett. 30, 3338–3340 (2005). [CrossRef]

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, 890–902 (2003). [CrossRef]

31.

V. A. Markel and J. C. Schotland, “Symmetries, inversion formulas and image reconstruction for optical tomography,” Phys. Rev. E 70, 056,616 (2004). [CrossRef]

32.

R. Choe, “Diffuse Optical Tomography and Spectroscopy of Breast Cancer and Fetal Brain,” Ph.D. thesis, University of Pennsylvania (2005).

33.

A. Ishimaru, Wave Propagation and Scattering in Random Media (IEEE Press, Piscataway, NJ, 1997).

34.

R. Aronson, “Boundary conditions for diffuse light,” J. Opt. Soc. Am. A 12, 2532–2539 (1995). [CrossRef]

OCIS Codes
(170.3880) Medical optics and biotechnology : Medical and biological imaging
(110.0113) Imaging systems : Imaging through turbid media

ToC Category:
Imaging Systems

History
Original Manuscript: January 28, 2008
Revised Manuscript: March 24, 2008
Manuscript Accepted: March 24, 2008
Published: March 28, 2008

Virtual Issues
Vol. 3, Iss. 4 Virtual Journal for Biomedical Optics

Citation
Soren D. Konecky, George Y. Panasyuk, Kijoon Lee, Vadim Markel, Arjun G. Yodh, and John C. Schotland, "Imaging complex structures with diffuse light," Opt. Express 16, 5048-5060 (2008)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-7-5048


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. M. C. W. van Rossum and T. M. Nieuwenhuizen, "Multiple scattering of classical waves: microscopy, mesoscopy and diffusion," Rev. Mod. Phys. 71, 313-371 (1999). [CrossRef]
  2. A. Yodh and B. Chance, "Spectroscopy and imaging with diffuse light," Phys. Today 34-40 (1995). [CrossRef]
  3. S. R. Arridge, "Optical tomography in medical imaging," Inverse Problems 15, R41-R93 (1999). [CrossRef]
  4. A. P. Gibson, J. C. Hebden, and S. R. Arridge, "Recent advances in diffuse optical imaging," Phys. Med. Biol. 50, R1-R43 (2005). [CrossRef] [PubMed]
  5. S. B. Colak, M. B. van der Mark, G. W. Hooft, J. H. Hoogenraad, E. S. van der Linden, and F. A. Kuijpers, "Clinical optical tomography and NIR spectroscopy for breast cancer detection," IEEE J. Sel. Top. Quantum Electron. 5, 1143-1158 (1999). [CrossRef]
  6. D. J. Hawrysz and E. M. Sevick-Muraca, "Developments toward diagnostic breast cancer imaging using nearinfrared optical measurements and fluorescent contrast agants," Neoplasia 2, 388-417 (2000). [CrossRef]
  7. Y. Xu, X. J. Gu, L. L. Fajardo, and H. B. Jiang, "In vivo breast imaging with diffuse optical tomography based on higher-order diffusion equations," Appl. Opt. 42, 3163-3169 (2003). [CrossRef] [PubMed]
  8. J. P. Culver, R. Choe, M. J. Holboke, L. Zubkov, T. Durduran, A. Slemp, V. Ntziachristos, B. Chance, and A. G. Yodh, "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, 235-247 (2003). [CrossRef] [PubMed]
  9. Q. I. Zhu, M. M. Huang, N. G. Chen, K. Zarfos, B. Jagjivan, M. Kane, P. Hedge, and S. H. Kurtzman, "Ultrasound-guided optical tomographic imaging of malignant and benign breast lesions: Initial clinical results of 19 cases," Neoplasia 5, 379-388 (2003). [PubMed]
  10. X. Intes, J. Ripoll, Y. Chen, S. Nioka, A. G. Yodh, and B. Chance, "In vivo continuous-wave optical breast imaging enhanced with Indocyanine Green," Med. Phys. 30, 1039-1047 (2003). [CrossRef] [PubMed]
  11. A. Li, E. L. Miller, M. E. Kilmer, T. J. Brukilacchio, T. Chaves, J. Stott, Q. Zhang, T. Wu, M. Chorlton, R. H. Moore, D. B. Kopans, and D. A. Boas, "Tomographic optical breast imaging guided by three-dimensional mammography," Appl. Opt. 42, 5181-5190 (2003). [CrossRef] [PubMed]
  12. S. D. Jiang, B. W. Pogue, T. O. McBride, M. M. Doyley, S. P. Poplack, and K. D. Paulsen, "Near-infrared breast tomography calibration with optoelastic tissue simulating phantoms," J. Electron. Imaging 12, 613-620 (2003). [CrossRef]
  13. D. A. Benaron, J. van Houten, D. C. Ho, S. D. Spilman, and D. K. Stevenson, "Imaging neonatal brain injury using light-based optical tomography," Pediatr. Res. 35, A378-A378 (1994).
  14. J. P. van Houten, W. F. Cheong, E. L. Kermit, T. R. Machold, D. K. Stevenson, and D. A. Benaron, "Clinical measurement of brain oxygenation and function using light-based optical tomography," Pediatr. Res. 39, 2273- 2273 (1996). [CrossRef]
  15. A. Y. Bluestone, G. Abdoulaev, C. H. Schmitz, R. L. Barbour, and A. H. Hielscher, "Three-dimensional optical tomography of hemodynamics in the human head," Opt. Express 9, 272-286 (2001). [CrossRef] [PubMed]
  16. G. Strangman, D. A. Boas, and J. P. Sutton, "Non-invasive neuroimaging using near-infrared light," Biol. Psychiatry 52, 679-693 (2002). [CrossRef] [PubMed]
  17. J. C. Hebden, A. Gibson, R. M. Yusof, N. Everdell, E. M. C. Hillman, D. T. Delpy, S. R. Arridge, T. Austin, J. H. Meek, and J. S. Wyatt, "Three-dimensional optical tomography of the premature infant brain," Phys. Med. Biol. 47, 4155-4166 (2002). [CrossRef] [PubMed]
  18. J. C. Hebden, "Advances in optical imaging of the newborn infant brain," Psychophysiology 40, 501-510 (2003). [CrossRef] [PubMed]
  19. J. P. Culver, A. M. Siegel, J. J. Stott, and D. A. Boas, "Volumetric diffuse optical tomography of brain activity," Opt. Lett. 28, 2061-2063 (2003). [CrossRef] [PubMed]
  20. A. M. Siegel, J. P. Culver, J. B. Mandeville, and D. A. Boas, "Temporal comparison of functional brain imaging with diffuse optical tomography and fMRI during rat forepaw stimulation," Phys. Med. Biol. 48, 1391-1403 (2003). [CrossRef] [PubMed]
  21. G. Q. Yu, T. Durduran, D. Furuya, J. H. Greenberg, and A. G. Yodh, "Frequency-domain multiplexing system for in vivo diffuse light measurements of rapid cerebral hemodynamics," Appl. Opt. 42, 2931-2939 (2003). [CrossRef] [PubMed]
  22. J. C. Schotland, "Continuous wave diffusion imaging," J. Opt. Soc. Am. A 14, 275-279 (1997). [CrossRef]
  23. V. A. Markel and J. C. Schotland, "Inverse scattering for the diffusion equation with general boundary conditions," Phys. Rev. E 64, R035,601 (2001). [CrossRef]
  24. V. A. Markel and J. C. Schotland, "The inverse problem in optical diffusion tomography. II. Inversion with boundary conditions," J. Opt. Soc. Am. A 19, 558-566 (2002). [CrossRef]
  25. V. A. Markel and J. C. Schotland, "Effects of sampling and limited data in optical tomography," Appl. Phys. Lett. 81, 1180-1182 (2002). [CrossRef]
  26. J. C. Schotland and V. A. Markel, "Inverse scattering with diffusing waves," J. Opt. Soc. Am. A 18, 2767-2777 (2001). [CrossRef]
  27. G. Turner, G. Zacharakis, A. Soubret, J. Ripoll, and V. Ntziachristos, "Complete-angle projection diffuse optical tomography by use of early photons," Opt. Lett. 30, 409-411 (2005). [CrossRef] [PubMed]
  28. R. Schulz, J. Ripoll, and V. Ntziachristos, "Noncontact optical tomography of turbid media," Opt. Lett. 28, 1701- 1703 (2003). [CrossRef] [PubMed]
  29. Z.-M. Wang, G. Y. Panasyuk, V. A. Markel, and J. C. Schotland, "Experimental demonstration of an analytic method for image reconstruction in optical tomography with large data sets," Opt. Lett. 30, 3338-3340 (2005). [CrossRef]
  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, 890-902 (2003). [CrossRef]
  31. V. A. Markel and J. C. Schotland, "Symmetries, inversion formulas and image reconstruction for optical tomography," Phys. Rev. E 70, 056,616 (2004). [CrossRef]
  32. R. Choe, "Diffuse Optical Tomography and Spectroscopy of Breast Cancer and Fetal Brain," Ph.D. thesis, University of Pennsylvania (2005).
  33. A. Ishimaru, Wave Propagation and Scattering in Random Media (IEEE Press, Piscataway, NJ, 1997).
  34. R. Aronson, "Boundary conditions for diffuse light," J. Opt. Soc. Am. A 12, 2532-2539 (1995). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited