Multiple Projection Optical Diffusion Tomography with Plane Wave Illumination

发布于:2021-10-14 08:36:44

arXiv:physics/0505196v1 [physics.med-ph] 27 May 2005

Multiple Projection Optical Di?usion Tomography with Plane Wave Illumination
Departments of Radiology and Bioengineering, University of Pennsylvania, Philadelphia, PA 19104 Abstract. We describe a new data collection scheme for optical di?usion tomography in which plane wave illumination is combined with multiple projections in the slab imaging geometry. Multiple projection measurements are performed by rotating the slab around the sample. The advantage of the proposed method is that the measured data can be much more easily ?tted into the dynamic range of most commonly used detectors. At the same time, multiple projections improve image quality by mutually interchanging the depth and transverse directions, and the scanned (detection) and integrated (illumination) surfaces. Inversion methods are derived for image reconstructions with extremely large data sets. Numerical simulations are performed for ?xed and rotated slabs.

Vadim A. Markel? and John C. Schotland?

PACS numbers: 87.57.Gg,42.30.Wb

Submitted to: Phys. Med. Biol.

? vmarkel@mail.med.upenn.edu ? schotland@seas.upenn.edu

Multiple Projection Optical Tomography 1. Introduction

2

Tomographic imaging with di?use light, often referred to as optical di?usion tomography (ODT), is a novel biomedical imaging modality [1, 2]. Although ODT was introduced more than a decade ago, e?orts to bring it into the clinical environment are hampered by relatively low quality and spatial resolution of images. Therefore, optimization of image reconstruction algorithms for high-resolution ODT is of fundamental importance. In this paper we study the image reconstruction problem of ODT by combining three novel approaches. First, we employ analytic image reconstruction methods which allows the utilization of extremely large data sets [3, 4]. Second, we make use of multiple projections [5]. Here by multiple projections we mean multiple orientations of the measurement apparatus with respect to the medium. Finally, we utilize the recently proposed plane wave illumination scheme [6]. Each of these methods provides an advantage which is not lost when the techniques are combined. We begin by brie?y reviewing the approaches to ODT imaging mentioned above. Note that throughout this paper we consider the slab imaging geometry which is often used in mammography and small-animal imaging [7,8]. In order to obtain multiple projection measurements, a pair of parallel plates are rotated around the medium to be imaged which is assumed to be stationary and unperturbed. There is a direct relationship between the spatial resolution of images and the number of data points used for reconstruction [3]. Indeed, the reconstruction of an image with N voxels, in principle, requires at least N measurements. In practice, the ill-posedness of the image reconstruction problem and the presence of noise require that this number be larger than N . Measurements with up to 1010 data points are feasible with CCD camera-based instruments. However, many previous studies of the image reconstruction problem in ODT have been limited to relatively small data sets (e.g., 256 data points in Ref. [9], 900 data points in Ref. [10]). This can be explained by the high computational complexity of algebraic image reconstruction algorithms which scales as O (N 3 ). To ameliorate this di?culty, we have recently introduced a family of analytic image reconstruction algorithms that can utilize extremely large data sets [11–15]. These methods allow a dramatic reduction in computational complexity which, in turn, leads to a signi?cant improvement of spatial resolution of images. However, these methods have certain limitations. First, the data collection method described in Ref. [14] requires that measurements are taken for source-detector pairs separated by a distance which is much larger than the slab thickness. In practice, such measurements are technically di?cult to perform. Reduction of the required dynamic range of the detectors can be achieved by using plane wave illumination [6]. Note that due to the general theoretical reciprocity of sources and detectors, plane wave illumination and scanned detection is equivalent to integrated detection and scanned narrow beam illumination. However, in a practical situation, the di?erent nature of illuminating and detecting devices must be taken into account. For the sake of de?nitiveness, we consider below plane wave illumination and combine it

Multiple Projection Optical Tomography

3

with analytic image reconstruction methods. Note that plane wave illumination requires time- or frequency-resolved measurements. However, it can be seen that the number of degrees of freedom in the data is still insu?cient for unique, simultaneous reconstruction of the absorption and di?usion (or reduced scattering) coe?cients. This situation is similar to the nonuniqueness demonstrated in Ref. [16]. Therefore, we focus here on the reconstruction of absorbing inhomogeneities assuming that the di?usion coe?cient of the medium is constant. Reconstruction of purely absorbing inhomogeneities have been employed, for example, in breast imaging [17–20] or blood oxygenation level imaging [21, 22]. Second, it was shown in Ref. [3] that in the slab imaging geometry the depth resolution (in the direction perpendicular to the slab) is fundamentally di?erent from the transverse resolution (in the direction parallel to the slab surface). The depth resolution is much more sensitive to noise and the point-spread functions (PSFs) in the depth direction strongly depend on the location of the inhomogeneity. This results in image artifacts. In general, the non-uniformity of the PSF can be a serious problem if more than one inhomogeneity is present. To correct this situation, we have recently proposed multi-projection image reconstruction methods [5, 15]. Multiple projections render the depth and transverse directions mutually interchangeable. As a result, the PSF becomes more uniform and less position-dependent, and also more sharply peaked. Note that multiple projections have been used in X-ray imaging for some time. However, an important di?erence between ODT and X-ray computed tomography is that, in the ?rst case, tomographic imaging is possible in principle with a single projection while in the second case it is not. Perhaps, due to this fact, multiple projections in optical tomography have not been investigated until recently, except for the case of ballistic propagation without scattering (e.g. [23]), or in conjunction with a modi?ed version of X-ray backprojection tomography with phenomenological corrections introduced to compensate for scattering [24, 25]. In Ref. [15] we have developed a general theoretical formalism for inverting measurements obtained from multiple projections. In Ref. [5] image reconstruction with two orthogonal projections was numerically implemented. In this paper we implement the more general image reconstruction algorithm of Ref. [15] for treatment of more than two projections in conjunction with plane wave illumination. Note that the plane wave illumination is advantageous when measurement are limited by the dynamic range of detectors. If the dynamic range is not an important experimental factor, the traditional measurement scheme with point sources and point detectors is expected to provide superior image quality. We combine the advantageous features of these two approaches with the computational e?ciency of the analytic image reconstruction methods.

Multiple Projection Optical Tomography 2. Theory 2.1. Single projection

4

We assume that propagation of multiply-scattered light in tissue is described by the di?usion equation. In addition, we will work in the frequency domain with the sources harmonically modulated at the frequency ω and detectors which yield the oscillatory part of transmitted intensity. Then the density of electromagnetic energy in the medium u(r ) obeys the di?usion equation ? D0 ?2 u(r) + [α(r) ? iω ]u(r) = S (r) , (1)

where α(r ) is the position dependent absorption coe?cients, S (r ) is the source function and the D0 is the di?usion coe?cient. Consider a slab of thickness L with the plane of incidence located at x = ?L/2 and the detection plane at x = L/2. The medium is located in the region ?L/2 < x < L/2. If point-like sources and detectors are used (typically, thin optical ?bers), the data can be expressed as a function φ(ω, ρs , ρd ), where ρs and ρd are two-dimensional vectors specifying the location of the sources and detectors, respectively, on the slab surfaces. Using the ?rst Born approximation, we linearize the forward model by decomposing the absorption function α(r ) into a constant background and a small ?uctuating part, α(r ) = α0 + δα(r ). We seek to reconstruct the values of δα(r ) from the data φ(ω, ρs, ρd ). The usual mathematical formulation of the ODT inverse problem is based on the integral equation [26] φ(ω, ρs, ρd ) = where Γ(ω, ρs, ρd ; r ) = d2 qs d2 qd κ(ω, qs , qd ; x) (2π )4 × exp [iqs · (ρ ? ρs ) + iqd · (ρd ? ρ)] , Γ(ω, ρs , ρd ; r )δα(r )d3r , (2)

(3)

ρ is the transverse part of the vector r (r = (x, ρ)) and the form of κ(ω, qs, qd ; x) is determined from the boundary conditions on the surfaces of the slab and the expression which relates the measurable intensity to the energy density u(r ). The derivation of (2),(3) and explicit expressions for κ are given in Ref. [3]. Note that the general form of (2),(3) follows from the symmetry of the problem and is independent of the di?usion approximation. Next, we introduce the plane wave illumination scheme. Instead of using point sources located at points ρs , we illuminate the slab with a normally incident wide homogeneous beam of su?ciently large diameter (compared to transverse dimensions of the slab). At the same time we utilize point detectors. This ensures that the new data function ψ (ω, ρd) de?ned by ψ (ω, ρd) = φ(ω, ρs, ρd )d2 ρs (4)

Multiple Projection Optical Tomography

5

has the same number of degrees of freedom as the unknown δα(r ) (two spatial directions and the frequency ω ). Thus, the inverse problem is well determined. The integral equation (2) can now be transformed to ψ (ω, ρd) = d2 q κ(ω, 0, q ; x) exp[iq · (ρd ? ρ)]δα(r )d3r . (2π )2 (5)

If ψ is measured for N di?erent modulation frequencies and the sources are placed on a square lattice with step size h, Eq. (5) can be inverted using the methods described in [14]. The SVD pseudo-inverse solution is given by δα(r ) = h2
FBZ

d2 u ?(ω ′, u) . (6) P ? (ω, u; r ) ω |M ?1(u)|ω ′ ψ exp(?iu · ρ) (2π )2 ′ ω,ω

Here the vector u is in the ?rst Brillouin zone (FBZ) of the lattice of sources, namely, ?π/h < uy,z ≤ π/h and P (ω, u; r ) =
v

κ(ω, 0, u + v ; x) exp(iv · ρ) ,

(7)

? + nz z ?). The elements where v are reciprocal lattice vectors of the form v = (2π/h)(ny y of matrix the M (u) are given by ω | M (u )| ω ′ = where ω |M1(q )|ω ′ =
L/2 ?L/2

M1 (u + v ) ,
v

(8)

κ(ω, 0, q ; x)κ?(ω ′ , 0, q ; x)dx

(9)

(the inverse matrix M ?1 (u) must be appropriately regularized [27]) and the Fourier ?(ω, u) is de?ned as transformed data function ψ ?(ω, u) = ψ
ρd

ψ (ω, ρd ) exp(iu · ρd ) .

(10)

Note that, if δα is reconstructed only at points which are commensurate with the lattice of sources, the factor exp(iv ·ρ) is equal to unity and the function P becomes independent of ρ. Note also that κ and M1 can be calculated in terms of elementary functions [15]. 2.2. Multiple projections We now consider inclusion of multiple projections. Let the sources and detectors be rotated around the sample as illustrated in Fig. 1. We assume that the rotations do not √ disturb the medium inside the cylinder x2 + y 2 < L/2 and that the unknown function δα vanishes outside the same region. The space inside the slab but outside the above cylindrical region is assumed to have the background values of the coe?cients α0 and D0 . Experimentally, this can be implemented, for example, by rotating an imaging apparatus around a sample suspended in matching ?uid. We introduce cylindrical coordinates r = (R, z, ?) with the z -axis being the axis of rotation. If the data are measured for

Multiple Projection Optical Tomography

6

θ

Figure 1. A sketch of the experimental set up with rotating slab. The axis of rotations (Oz ) is perpendicular to the plane of the ?gure and coincides with the axis of the cylinder R < L/2 inside which reconstructions are performed. Locations of sources and detectors are given in a local reference frame which rotates together with the slab.

Nθ di?erent orientations, where the respective angles θn are equally spaced and given by θn = 2π (n ? 1)/Nθ , n = 1, . . . , Nθ , the reconstruction formula (6) can be generalized to [15]: 2πh2 Nθ
Nθ π/h

δα(r ) =

n=1 ?π/h

duz exp[?i(uz z + n?)] 2π ω,ω ′

π/h ?π/h

duy

π/h ?π/h

du′y P ? (ω, u, n; r ) (11)

?(ω ′ , u′ , uz , n) . × ω, uy |M ?1 (uz , n)|ω ′, u′y ψ y Here


P (ω, u, n; r ) = a(ω, q , m; R) =
k =?∞ v 2π 0

a(ω, u + v , n + Nθ k ; R) exp[i(Nθ k? + vz z )] , κ(ω, 0, q ; R cos ?) exp[i(qy R sin ? ? m?)]d? ,

(12) (13)

the elements of the matrix M (uz , n) are given by


ω, uy |M (uz , n)|ω ′, u′y =

′ vz k =?∞ vy ,vy

′ ω, uy + vy |M1 (uz + vz , n + Nθ k )|ω ′, u′y + vy ,

(14)
′ ω, qy |M1 (qz , m)|ω ′, qy = L/2 0 ′ a(ω, qy , qz , m; R)a? (ω ′ , qy , qz , m; R)RdR

(15)

and the Fourier-transformed data function is ?(ω, u, n) = ψ
ρd ,θ

ψ (ω, ρd, θ) exp[i(u · ρd + nθ)] .

(16)

Multiple Projection Optical Tomography

7

Note that in (16) we have explicitly included the dependence of the data function on the ′ angle of orientation θ. The functions a(ω, q , m; R) and ω, qy |M1 (qz , m)|ω ′ , qy can be, in general, expressed in terms of modi?ed Bessel functions. The corresponding integrals (13) and (15) are calculated in the Appendix for the case of purely absorbing boundaries. A few comments on the reconstruction formula (11) are necessary. First, there is an apparent di?erence between the variables uz , n and uy , ω . The ?rst set of variables correspond (after Fourier transformation of the data) to the variables z, θ. These are the variables with respect to which the unperturbed medium is translationally invariant, and they can be referred to as “external” variables. The variables ω, uy are “internal” variables: they do not correspond to any translational invariance of the system. Second, the reconstruction algorithm (16) involves integration over the continuous variables uy and u′y and inversion of the operator M (uz , n) whose matrix elements depend on continuous indices. However, if the variables uy , u′y are discretized and the corresponding integration in (16) is replaced by a summation, then M (uz , n) becomes a discrete matrix. The resulting reconstruction formula is no longer an SVD pseudo-inverse on the whole set of data ψ (ω, ρd , θ). However, it is a pseudo-inverse solution on the set of the Fourier?(ω, uy , uz , θ) where uy takes only discrete values. Third, it can transformed data ψ be veri?ed that in the case Nθ = 1, the reconstruction formula (16) reduces to (6). ? is four Fourth, we note that the number of degrees of freedom in the data-function ψ (ω, uy , uz and n). Thus, when the number of rotations is large, it is su?cient to use only one or a few values of the variable uy , in which case the inverse problem is still well determined. It can be argued that the reconstruction algorithm is then “numerical” in one dimension and “analytic” in two. § However, when only a small number of projections is taken, we must use a relatively large number of discrete values of uy . By doing so, we increase the size of the matrix M whose SVD must be found numerically. The inverse solution (11) is then “numerical” in two dimensions and “analytic” in one. A similar algorithm (numerical in two dimensions and analytic in one dimension) was proposed and implemented in [5], where the image reconstruction area was rectangular rather than cylindrical, but only two orthogonal projections were allowed. In contrast, the full potential of the image reconstruction algorithm proposed here is realized when Nθ is large. 3. Numerical Results 3.1. Single projection We have implemented the proposed reconstruction algorithm using computer-generated data and the following parameters: the slab thickness was chosen to be the same as the cw di?use wavelength, L = 2π D0 /α0 (for most biological tissues, this corresponds to
§ If the number of rotations and the number of discrete values of uy are both large, it should be possible to recover the absorption and scattering coe?cients uniquely and simultaneously. This theoretical possibility is not discussed in this paper.

Multiple Projection Optical Tomography
d d0 0.125L d d0 0.25L d d0 0.375L

8

d d0 0.5L

d d0 0.625L

d d0 0.75L

d d0 0.875L

a ?Α a.u.

b

x L

x L

-0.5

0

0.5 -0.5

0

0.5

Figure 2. Tomographic slices parallel to the slab surface drawn through the medium at di?erent depths d (from the plane of scanned detection) with the small absorber lying in the center of the ?eld of view at the same depth d0 = d, and the pointspread functions representing depth resolution (a,b). The curves are plotted on the same scale (a) and normalized to their own maxima (b). For curves (a,b), the point absorber depth is d0 = 0.25L (solid line), d0 = 0.5L (short dash) and d0 = 0.75L (long dash).

L ? 6cm); the lattice step was chosen to be h = L/40 and we have used N = 25 di?erent modulation frequencies which range from ω = 0 to ω = 10α0 (the maximum frequency corresponds to ? 1.6GHz); the ?eld of view was chosen to be L × L and, ?nally, we have generated forward data for a single point (delta-function) absorber which is located in the center of the ?eld of view but at di?erent depths. Absorbing boundary conditions were imposed on the surface of the slab. The corresponding expression for the function κ(ω, 0, q ; x) is given in the Appendix. The results of reconstructions are shown in Fig. 2. The density plots represent tomographic slices of the medium drawn at di?erent depths d (the distance from the plane of scanned detection) parallel to the slab surfaces. The depth of the absorbing inhomogeneity, d0 , was in each case equal to d; thus the slices represent the depthdependent y ? z PSFs. Each density plot has a linear color scale and is normalized to its own maximum. As expected, the PSFs become broader when the point absorber approaches the illuminated plane. The last two panels (a,b) show the PSFs in the depth direction (x) for point absorbers located at d0 = 0.25L, d0 = 0.5L and d0 = 0.75L.

Multiple Projection Optical Tomography

9

Note that the approximate half-widths of these curves are 0.06L, 0.09L and 0.09L, respectively. The analysis of Fig. 2 suggests that the PSFs are depth-dependent. Moreover, the PSFs have di?erent integral weights. Thus, the point absorbers which are closer to the plane of scanned detection result in higher peaks in the reconstructed images. The width of the PSFs also depends on depth of the point absorber. This potentially constitutes a serious problem for three dimensional tomographic imaging. 3.2. Multiple projections We have implemented numerically the multi-projection image reconstruction formula (11). Note that in the multi-projection case there are two choices for graphically representing the tomographic slices. In one case, the slices are perpendicular to the axis of rotation. The image then is reconstructed in a circle. This choice is convenient for studying the radial and angular resolutions. Another possibility is to construct cylindrical slices R = Rimage = const, and map them onto rectangles. The image is then reconstructed in the rectangular area 2πRimage × (zmax ? zmin ), where zmax and zmin are the maximum and minimum values of z , chosen arbitrarily. We start with the discussion of circular slices. The results of numerical implementation of the reconstruction formula (16) are shown in Fig. 3 for four di?erent orientations of the slab, namely θ = 0, π/2, π, 3π/2. We have used 23 equally spaced values of uy ranging from ?π/h to π/h and 15 equally spaced modulation frequencies ranging from 0 to 10α0 ; otherwise, the parameters are the same as in Fig. 2. The inhomogeneity was located as speci?ed in the ?gure legend. The white spots in the density plots illustrate the depth PSFs. The graphs (a,b) show the same PSFs in a more quantitative way by plotting δα along the diameter of the cylinder which intersects all three inhomogeneities. As expected, using four di?erent projections improves the image quality by interchanging the source and detector planes, and the depth and transverse directions. Moreover, using more projections than four does not change the results substantially, as is illustrated in Figs. 4 and 5. However, when a large number of projections is taken, the inverse problem becomes well determined even when a relatively small number of “internal” degrees of freedom uy is used. This makes the reconstruction formulae computationally e?cient. Thus, the computation time required for producing data for Fig. 5 is more than an order of magnitude less than that for Fig. 3, yet the image quality is similar. We have veri?ed that three discrete values of uy is also su?cient for Nθ = 20 (taking a single value uy = 0 results in a slight decrease in image quality; data not shown). Although the images shown in Figs. 3-5 are similar, the best image quality is, in fact, attained in Fig. 4. Here the approximate half-widths of the PSF in the R direction are 0.05L for the inhomogeneity located at R0 = 0, 0.04L for the inhomogeneity at R0 = 0.25L, ?0 = 0; and 0.03L for the inhomogeneity at R0 = 0.375L, ?0 = π . These

Multiple Projection Optical Tomography
z 0 R0 0 R0 0.25L,
0

10

0 R0 0.375L,

0

Π

z 0.05L R0 0 R0 0.25L,
0

0 R0 0.375L,

0

Π

z 0.1L R0 0 R0 0.25L,
0

0 R0 0.375L,

0

Π

a ?Α a.u.

b ?Α a.u.

x L

x L

-0.5

0

0.5 -0.5

0

0.5

Figure 3. Circular slices illustrating radial, angular and z resolution. All point absorbers are in the z = 0 plane, and the point-spread functions representing depth (R) resolution (a,b). The radial and angular coordinates of the point absorber, R0 and ?0 , are speci?ed in the ?gure legends. First row of images: slices at z = 0; second row: slices at z = 0.05L; third row: slices at z = 0.1L. Images (a-b): reconstruction along the diameter that crosses all three inhomogeneities. In (a,b) solid line corresponds to R0 = 0.375L and ?0 = π , short dash to R0 = 0 and long dash to R0 = 0.25L and ?0 = 0 Four projections, 15 modulation frequencies and 23 discrete values of uy are used.

Multiple Projection Optical Tomography
z 0 R0 0 R0 0.25L,
0

11

0 R0 0.375L,

0

Π

z 0.05L R0 0 R0 0.25L,
0

0 R0 0.375L,

0

Π

z 0.1L R0 0 R0 0.25L,
0

0 R0 0.375L,

0

Π

a ?Α a.u.

b

x L

x L

-0.5

0

0.5 -0.5

0

0.5

Figure 4. Same as in Fig. 3 but 20 projections are used.

Multiple Projection Optical Tomography
z 0 R0 0 R0 0.25L,
0

12

0 R0 0.375L,

0

Π

z 0.05L R0 0 R0 0.25L,
0

0 R0 0.375L,

0

Π

z 0.1L R0 0 R0 0.25L,
0

0 R0 0.375L,

0

Π

a ?Α a.u.

b

x L

x L

-0.5

0

0.5 -0.5

0

0.5

Figure 5. Same as in Fig. 3 but 40 projections and only three discrete values of uy are used.

Multiple Projection Optical Tomography

13

values should be compared to the respective values given in the discussion of Fig. 2. In particular, the inhomogeneity located at h0 = 0.5L in Fig. 2 corresponds to the inhomogeneity at R0 = 0 in Figs. 3-5 and is the most “di?cult” to reconstruct since it is located deep inside the medium. It can be seen that the PSF half-width in the image of this particular inhomogeneity is reduced by approximately the factor of 2 due to the use of multiple projections. In addition, the relative heights of the maxima of the PSFs in Fig. 3-5 do not di?er as much as in Fig. 2. This is expected to reduce image artifacts. Now we consider the cylindrical slices. From the computational point of view, the use of cylindrical slices is a more natural way to display reconstructed images. This is evident from the inversion formulae (11),(12). Indeed, it can be seen that when the reconstructed image is rasterized so that the variables z and ? are placed on lattices with steps h and 2π/Nθ , respectively, the function P (ω, u, n; R, z, ?) becomes independent of z and ?. Then the dependence of reconstructed images on these two variables is only due to the exponent in the integral (11) and the reconstruction formula, with respect to these two variables, is reduced to a Fourier transform. In Fig. 6 we have used three discrete values of uy with 40 di?erent projections and slices are drawn as described in the ?gure caption. Fig. 6(a) illustrates image reconstruction with noiseless data. It can be directly compared to slices shown in Fig. 2. To demonstrate the stability of image reconstruction, we have added random Gaussian noise to the data function at the level of 1% of the average absolute value of the data. The result is shown in Fig. 6(b). As is well known, inclusion of noise tends to decrease spatial resolution. It can be seen that this e?ect is stronger for inhomogeneities that are deeper inside the medium. We have demonstrated earlier that multi-projection imaging is more stable in the presence of noise than the single projection technique [5]. 4. Summary In summary, we have presented a new experimental modality and computationally e?cient image reconstruction algorithms for optical di?usion tomography employing plane wave illumination with multiple projections. Note that due to reciprocity, plane wave illumination and scanned detection is equivalent to illumination by a scanned narrow beam and integrated detection (e.g., with the use of time-resolved CCD camera). The following speci?c conclusions can be drawn ? Use of plane wave illumination may be simpler experimentally than the traditional approach in which point-like sources and detectors are scanned because measurements with a much smaller dynamic range are required. ? In a single projection experiment, the image quality is relatively good when the point absorber is close to the scanned surface and deteriorates as it approaches the plane of illumination. This situation should be contrasted with the traditional point source/point detector modality [14], where the image quality is low for inhomogeneities located in the center of a slab and improves when the

Multiple Projection Optical Tomography
R R0 0.9 L 2 R R0 0.9 L 2

14

R R0 0.8 L 2

R R0 0.8 L 2

R R0 0.7 L 2

R R0 0.7 L 2

R R0 0.6 L 2

R R0 0.6 L 2

R R0 0.5 L 2

R R0 0.5 L 2

R R0 0.4 L 2

R R0 0.4 L 2

R R0 0.3 L 2

R R0 0.3 L 2

(a)

(b)

Figure 6. Cylindrical slices illustrating z and ?-resolution for zero noise level (a) and for 1% noise-to-signal ratio (b). The point absorbers are located in the z = 0 plane at radial depths R0 as indicated. The cylindrical surfaces with radii R = R0 (directly intersecting the inhomogeneity) are shown as projections onto a plane; the length of the vertical side of each rectangle is equal to L and of the horizontal side to 2πR. Forty projections, 25 modulation frequencies and 9 discrete values of uy are used for reconstruction.

Multiple Projection Optical Tomography

15

inhomogeneity approaches either of the imaging surfaces. For a point inhomogeneity in the center of a slab, the image quality is slightly better for the traditional (point source/point detector) modality (cf. [14]). ? Rotating the slab around the sample removes many of the de?ciencies of the plane wave illumination scheme by interchanging the scanned and integrated detection surfaces and depth and transverse directions. A minimum of four projections is required for such an interchange. ? When only four rotations are used, a large number of discrete values of the “internal” variable uy must be utilized in the reconstruction. Alternatively, a large number of projections can be used with a small number of discrete values of uy . The second approach is much more computationally e?cient but requires more complicated measurements. The quality of images is similar in both cases. ? The plane wave illumination approach allows one to signi?cantly simplify reconstruction formulae, both in single- and multiple-projection imaging. ? If only small number of projections is used (two or four) an alternative approach may be used, which is purely numerical in two dimensions and analytic in one dimension [5]. For a large number of projections, the algorithm reported here is computationally more e?cient. This work was supported in part by the AFOSR under the grant F41624-02-1-7001 and by the NIH under grant P41RR0205. Appendix: Calculation of the functions a(ω, q , m; R) and M1 (qz , m). The function a(ω, q , m; R) is de?ned by (13). To evaluate the integral, we must specify the function κ(ω, 0, q ; x). Explicit expressions for κ are given in [3] for general boundary conditions. In this paper we consider absorbing boundaries for which κ is given by the expression κ(ω, 0, q ; x) = ?? D0
2

sinh[k (L/2 ? x)] sinh[Q(L/2 + x)] , sinh(kL) sinh(QL)

(A1)

where ?? = 3D0 /c is the transport mean free path, c is the average speed of light in the √ medium, k = (α0 ? iω )/D0 is the complex di?use wavenumber, Q = q 2 + k 2 and q = (qy , qz ). Generalization to mixed boundaries of Robin type is straightforward and is not discussed here. Then, the expression for a(ω, q , m; R) becomes ?? D0
2

a(ω, q , m; R) = ×
2π 0

1 sinh(kL) sinh(QL)

sinh[k (L/2 ? R cos ?)] sinh[Q(L/2 + R cos ?)] exp[i(qy R sin ? ? m?)]d? . (A2)

Multiple Projection Optical Tomography This can be equivalently rewritten as ?? D0
2

16

a(ω, q , m; R) =

1 4 sinh(kL) sinh(QL)

× exp (Q + k )L/2 Fm (Q ? k )R, iqy R ? exp (?Q + k )L/2 Fm (?Q ? k )R, iqy R ? exp (Q ? k )L/2 Fm (Q + k )R, iqy R + exp (?Q ? k )L/2 Fm (?Q + k )R, iqy R (A3) where
2π 0

Fm (u, v ) =

exp[u cos ? + v sin ? ? im?]d? = 2π



u2 + v 2 u + iv

m

√ Im ( u2 + v 2 ) , (A4)

and Im (x) is the modi?ed Bessel function of the ?rst kind. Note that (A4) is well de?ned, including the case v = iu. The expressions (A3) and (A4) de?ne a(q , m; R). Next, we need to calculate the matrix elements of M1 (qz , m). This integral contains sixteen terms of the form s1 s2 s3 s4 π 2 (?? /D0 )2 exp [(s1 k + s2 Q + s3 k ′ + s4 Q′ )L/2] 4 sinh(kL) sinh(QL) sinh(k ′ L) sinh(Q′ L) ×? ?
0

? ?

(?s1 k + s2
L/2

Q)2

?

2 qy

(?s3

k′

+ s4

Q′ )2 ?

?

′ )2 (qy

(?s1 k + s2 Q ? qy )(?s3

k′

+ s4

Q′

′) qy

?m ? ? ?

×

2 I ′ ′ 2 ′ 2 Im R (?s1 k + s2 Q)2 ? qy m R (?s3 k + s4 Q ) ? (qy ) RdR .

(A5)

where sk = ±1, the sixteen terms correspond to sixteen possible permutations of the signs of sk and the primed variables should be understood as follows: k ′ = ′ )2 + q 2 + (k ′ )2 . The integral in (A5) is evaluated with (α0 ? iω ′ )/D0 and Q′ = (qy z the use of
? ? ? ? ? ? ? ? ? ? ?

c 0

xIn (ax)In (bx)dx =

c [aI (ac)I (bc) ? bI (ac)I (bc)] n n n+1 a2 ? b2 n+1
2 2 1 c2 + n2 I 2 (ac) ′ ?c [ I ( ac )] + n 2 2 a2 n

, ,

a=b, (A6) a=b.

This completely de?nes all the functions necessary for implementation of the multiprojection reconstruction algorithm.

Multiple Projection Optical Tomography References
[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]

17

S. R. Arridge, Inverse Problems 15, R41 (1999). D. A. Boas et al., IEEE Signal Proc. Mag. 18, 57 (2001). V. A. Markel and J. C. Schotland, J. Opt. Soc. Am. A 19, 558 (2002). V. A. Markel and J. C. Schotland, J. Opt. Soc. Am. A 20, 890 (2003). V. A. Markel and J. C. Schotland, Opt. Lett. 29, 2019 (2004). M. Xu, M. Lax, and R. R. Alfano, J. Opt. Soc. Am. A 18, 1535 (2001). M. Franceschini et al., Proc. Natl. Acad. Sci. USA 94, 6468 (1997). V. Ntziachristos, A. Yodh, M. Schnall, and B. Chance, Proc. Natl. Acad. Sci. USA 97, 2767 (1999). B. W. Pogue, T. O. McBride, U. L. Ostererg, and K. D. Paulsen, Opt. Express 4, 270 (1999). J. P. Culver, V. Ntziachristos, M. J. Holboke, and A. G. Yodh, Opt. Lett. 26, 701 (2001). J. C. Schotland, J. Opt. Soc. Am. A 14, 275 (1997). J. C. Schotland and V. A. Markel, J. Opt. Soc. Am. A 18, 2767 (2001). V. A. Markel and J. C. Schotland, Phys. Rev. E 64, R035601 (2001). V. A. Markel and J. C. Schotland, Appl. Phys. Lett. 81, 1180 (2002). V. A. Markel and J. C. Schotland, Phys. Rev. E 70, 056616(19) (2004). S. R. Arridge and W. R. B. Lionhart, Opt. Lett. 23, 882 (1998). S. B. Colak et al., IEEE J. Selected Topics in Quantum Electronics 5, 1143 (1999). D. J. Hawrysz and E. M. Sevick-Muraca, Neoplasia 2, 388 (2000). J. P. Culver et al., Med. Phys. 30, 235 (2003). X. Intes et al., Med. Phys. 30, 1039 (2003). J. P. van Houten et al., Pediatric Research 39, 2273 (1996). J. P. Culver et al., J. of Cerebral Blood Flow and Metabolism 23, 911 (2003). C. S. Brown, D. H. Burns, F. A. Spelman, and A. C. Nelson, Appl. Opt. 31, 6247 (1992). S. B. Colak et al., Appl. Opt. 36, 180 (1997). C. L. Matson and H. L. Liu, J. Opt. Soc. Am. A 16, 1254 (1999). C. P. Gonatas, M. Ishii, J. S. Leigh, and J. C. Schotland, Phys. Rev. E 52, 4361 (1995). F. Natterer, The mathematics of computerized tomography (Wiley, New York, 1986).


相关推荐

最新更新

猜你喜欢