BACKGROUND OF THE DISCLOSURE

[0001]
1.) Field of the Invention

[0002]
This invention relates to tomography and, more particularly, to method and concomitant system wherein an image of an object is directly reconstructed from measurements of scattered radiation using a nonlinear reconstruction technique.

[0003]
2.) Description of the Background Art

[0004]
There has been considerable interest in the inverse scattering problem (ISP) for diffuse light. The basic physical problem consists of reconstructing the spatial distribution of the optical absorption and diffusion coefficients inside a highlyscattering medium from intensity measurements on the boundary of the medium.

[0005]
The equations describing scattering of diffuse light from fluctuations in the absorption and diffusion coefficients oz and D are, in general, nonlinear. Thus numerical methods for solving the nonlinear inverse problem have been widely studied and are typically based upon Newton's method. A limitation of this approach is its computational complexity which arises from the fact that the forward problem must be solved at each iteration of the algorithm.

[0006]
Approaches to the inverse problem based upon linearization of the forward problem have also been explored. In this method, the integral equations of diffuse light propagation are expanded and linearized in α and D. These equations can then be solved with the use of analytic inversion formulas. The use of inversion formulas is especially attractive due to computational efficiency. Representative of this technique are the disclosures of U.S. Pat. No. 5,905,261, the Background section of which is incorporated herein by reference.

[0007]
The art is devoid of a methodology, and concomitant system, wherein the nonlinear equations describing diffusion and absorption of an image are directly solved to thereby effect direct, but generalized, reconstruction of the image. That is, in the past, only explicit inversion formulas for the case of linearized ISP have been obtained; explicit inversion formulas for nonlinear ISP case have not been devised.
SUMMARY OF THE INVENTION

[0008]
These and other shortcomings are obviated in accordance with the present invention via a technique whereby an object is irradiated with a source of radiation and then waves diffusively scattered from the object are processed with a prescribed nonlinear mathematical algorithm to reconstruct the tomographic image.

[0009]
In accordance with one broad method aspect of the present invention, the method for generating a tomographic image of an object includes: (i) irradiating the object with a source of radiation; (ii) measuring a transmitted intensity due to diffusively scattered radiation wherein said transmitted intensity is related to at least one coefficient characterizing the image by a nonlinear integral operator; and (iii) directly reconstructing the image by executing a prescribed mathematical algorithm, determined with reference to said integral operator, on said transmitted intensity.

[0010]
In accordance with another broad method aspect of the present invention, the method for generating a tomographic image of an object ineludes: (i) irradiating the object with a source of radiation; (ii) measuring a transmitted intensity due to diffusively scattered radiation wherein said transmitted intensity is related to at least one coefficient characterizing the image by a nonlinear integral operator; and (iii) directly reconstructing the image by executing a prescribed mathematical algorithm, determined with reference to said integral operator, on said transmitted intensity, said algorithm further relating said coefficient to said transmitted intensity by another nonlinear integral operator.

[0011]
In accordance with yet another broad method aspect of the present invention, the method for generating a tomographic image of an object includes: (i) irradiating the object with a source of radiation, (ii) measuring a transmitted intensity due to diffusively scattered radiation wherein said transmitted intensity is related to at least one coefficient characterizing the image by a nonlinear integral operator, and (iii) directly reconstructing the image by executing a prescribed mathematical algorithm, determined with reference to said integral operator, on said transmitted intensity, said mathematical algorithm expressed as a functional series expansion for said coefficient in powers of said transmitted intensity.

[0012]
Broad system aspects of the present invention are commensurate with the broad method aspects.

[0013]
The features of the this approach are at least twofold: (i) the approach provides a formally exact solution to the ISP in diffusion tomography. The approach may be viewed as a nonlinear inversion formula whose first term coincides with the pseudoinverse solution to the linearized ISP. The higher order terms represent systematically improvable nonlinear corrections which, in principle, can be computed to arbitrarily high order. Thus, it is only necessary to solve the linear ISP in order to formally solve the nonlinear ISP; and (ii) the approach to the ISP differs from Newtontype iterative methods. This follows from the fact that such prior methods require the conventional forward problem to be solved for each iteration.
BRIEF DESCRIPTION OF THE DRAWING

[0014]
[0014]FIG. 1 is a highlevel block diagram of a system for directly reconstructing the tomographic image of the scatterer

[0015]
FIGS. 2A2D depict direct reconstruction results for an exemplary scatterer using different ratios of parameters;

[0016]
FIGS. 3A3D depict direct reconstruction results for an exemplary scatterer using the same parameters different from FIGS. 2A2D that exemplify increased resolution with higherorder terms;

[0017]
[0017]FIG. 4 depicts the convergence of the higher order corrections to the true profile;

[0018]
[0018]FIG. 5A is a highlevel flow diagram of the methodology for directly reconstructing the tomographic image of the scatterer; and

[0019]
[0019]FIG. 5B is another highlevel flow diagram of the methodology for directly reconstructing the tomographic image of the scatterer commensurate with the functional series technique.
BACKGROUND OF THE DISCLOSURE

[0020]
1. System

[0021]
As depicted in highlevel block diagram form in FIG. 1, system 100 is a tomography system for generating an image of an scatterer/object using measurements of scattered waves emanating from an object in response to waves illuminating the object. In particular, object 105 is shown as being under investigation. System 100 is composed of: source 120 for probing the object 105; data acquisition detector 130 for detecting the scattering data corresponding to the scattered waves from object 105 at one or more locations proximate to object 105; position controller 140 for controlling the locations of detectors 130 and sources 120; and computer processor 150, having associated input device 160 (e.g., a keyboard) and output device 170 (e.g., a graphical display terminal). Computer processor 150 has as its inputs positional information from controller 140 and the measured scattering data from detector 130.

[0022]
Computer 150 stores a computer program which implements the direct reconstruction algorithm; in particular, the stored program processes the measured scattering data to produce the image of the object under study using a prescribed mathematical algorithm. The algorithm is, generally, based upon a functional expansion of the absorption and diffusion coefficients of the object in terms of tensor products of the measured scattering data. The algorithm will be described in detail below.

[0023]
2. Overview of the Underlying Mathematical Formalism

[0024]
We begin by setting forth the relevant mathematical formalism which serves as a backdrop to the point of departure in accordance with the present invention. We assume that the energy density u(r, t) of diffuse light in an inhomogeneous medium obeys the diffusion equation
$\begin{array}{cc}\frac{\partial u\ue8a0\left(r,t\right)}{\partial t}=\nabla \xb7\left[D\ue8a0\left(r\right)\ue89e\nabla u\ue8a0\left(r,t\right)\right]\alpha \ue8a0\left(r\right)\ue89eu\ue8a0\left(r,t\right)+S\ue8a0\left(r,t\right),& \left(1\right)\end{array}$

[0025]
where α(r) and D(r) are the positiondependent absorption and diffusion coefficients, and S(r, t) is the power density of the source. We further assume that the source is harmonically modulated with angular frequency ω. In addition to (1), the energy density must satisfy boundary conditions on the surface of the medium (or at infinity in the case of free boundaries) of the general form

u+l{circumflex over (n)}·∇u=0, (2)

[0026]
where l is the socalled extrapolation length and {circumflex over (n)} is an outward pointing normal. Note that when l=0 we obtain purely absorbing boundaries and when l→∞ purely reflecting boundaries.

[0027]
In general, the socalled Green's function may be directly related to the intensity measured by a point detector when the medium is illuminated by a point source. The Green's function G(r_{1}, r_{2}) for the frequencydomain diffusion equation obeys the integral equation

(r′,r _{2})=G _{0}(r _{1} ,r _{2})−∫d ^{3} rG _{0}(r _{1} ,r)V(r)G(r,r _{2}), (3)

[0028]
where G_{0 }is the Green's function for a homogeneous medium with absorption α_{0 }and diffusion constant D_{0}. We have also introduced the notation

V(r)≡δα(r)−∇·δD(r)∇, (4)

[0029]
where δα(r)=α(r)−α
_{0 }and δD(r)=D(r)−D
_{0}. The unperturbed Green's function G
_{0}(r, r′) obeys the boundary condition (2) and satisfies
$\begin{array}{cc}\left({\nabla}^{2}\ue89e{k}^{2}\right)\ue89e{G}_{0}\ue8a0\left(r,{r}^{\prime}\right)=\frac{1}{{D}_{0}}\ue89e\delta \ue8a0\left(r{r}^{\prime}\right),& \left(5\right)\end{array}$

[0030]
where the diffuse wave number k is given by
$\begin{array}{cc}{k}^{2}=\frac{{\alpha}_{0}\uf74e\ue89e\text{\hspace{1em}}\ue89e\omega}{{D}_{0}}.& \left(6\right)\end{array}$

[0031]
It can be shown that the change in the intensity of transmitted light (at the modulation frequency ω) due to spatial fluctuations in α(r) and D(r) is given by the integral equation

Φ(r _{1} ,r _{2})=β∫G _{0}(r _{1} ,r)V(r)G(r,r _{2})d ^{3} r. (7)

[0032]
Here the data function Φ(r
_{1}, r
_{2}) is proportional to the change in intensity relative to a reference medium with absorption α
_{0 }and diffusion constant D
_{0}, r
_{1 }and r
_{2 }denote the coordinates of the source and detector, and
$\begin{array}{cc}\beta =\{\begin{array}{cc}1& \mathrm{for}\ue89e\text{\hspace{1em}}\ue89e\mathrm{free}\ue89e\text{\hspace{1em}}\ue89e\mathrm{boundaries}\\ {\left(1+{l}^{*}/l\right)}^{2}& \mathrm{for}\ue89e\text{\hspace{1em}}\ue89e\mathrm{boundary}\ue89e\text{\hspace{1em}}\ue89e\mathrm{conditions}\ue89e\text{\hspace{1em}}\ue89e\mathrm{expressed}\ue89e\text{\hspace{1em}}\ue89e\mathrm{by}\ue89e\text{\hspace{1em}}\ue89e\left(2\right)\end{array}& \left(8\right)\end{array}$

[0033]
with l*=3D_{0}/c.

[0034]
The forward problem in diffusion tomography is defined as the problem of computing the data function Φ from the scattering potential η=(δα, δD). More precisely, the integral equation (7) may be regarded as defining a nonlinear operator K from the Hilbert space of scattering potentials H_{1 }into the Hilbert space of scattering data H_{2}. The fact that K is nonlinear may be understood by examining the perturbation expansion for Φ in powers of V. Using the series expansion for the Green's function G, which can be obtained by iterating the integral equation (3), and the definition of the data function (7) we obtain the required expansion:

Φ(r _{1} ,r _{2})=β∫d ^{3} rG _{0}(r _{1} ,r)V(r)G _{0}(r,r _{2})+β∫d ^{3} rd ^{3} r′G _{0}(r _{1} ,r)V(r)G _{0}(r,r′)V(r′)G _{0}(r′,r _{2})+ (9)

[0035]
If only the first term in the series is retained we refer to this as the weakscattering approximation.

[0036]
The inverse problem in diffusion tomography is defined as recovering η from measurements of Φ. The standard numerical approach to this nonlinear problem is to employ a functional Newton's method. This results in an iterative algorithm of the form

η_{n+1}=η_{n} +M _{n} ^{+}(Φ−K{η _{n}}), n=1,2, . . . , (10)

[0037]
where M
_{n} ^{+}denotes the pseudoinverse of the functional derivative
$\begin{array}{cc}{M}_{n}=\frac{\delta \ue89e\text{\hspace{1em}}\ue89eK}{\delta \ue89e\text{\hspace{1em}}\ue89e\eta}\ue89e{}_{\eta ={\eta}_{n}}.& \left(11\right)\end{array}$

[0038]
In accordance with the present inventive subject matter, we consider an alternative to the use of Newton's method. In particular, we construct a formally exact analytic solution to the nonlinear ISP. This solution, which we refer to as the inverse scattering series, has the form of a functional series expansion for 7 in powers of the data function 4). The first term in the expansion corresponds to the pseudoinverse solution to the linearized inverse problem. The higher order terms may be interpreted as nonlinear corrections to the singular value decomposition (SVD) inversion formulas of the linearized inverse problem.

[0039]
The remainder of this description is organized as follows. In Section 2 we derive the inverse scattering series for diffusion tomography in its most general form, independent of geometry and the type of boundary conditions.

[0040]
In Section 3 we consider the inverse problem in the biplanar geometry.

[0041]
Section 4 details exemplary numerical results for the nonlinear reconstruction of a spherical inhomogeneity in the biplanar geometry.

[0042]
Section 5 discusses a flow diagram of the process in accordance with the present invention.

[0043]
Section 6 presents mathematical properties of the inverse scattering series and the derivation of the data function for a spherical inhomogeneity in the biplanar geometry with free boundaries.

[0044]
2. INVERSE PROBLEM—Inverse Scattering Series

[0045]
In this section we present the construction of the inverse scattering series for diffusion tomography.

[0046]
The scattering series (9) can be rewritten in the form

Φ(r′,r _{2})=∫d ^{3} rK _{1} ^{i}(r _{1} ,r _{2} ,r)η_{i}(r)+∫d ^{3} rd ^{3} r′K _{2} ^{ij}(r _{1} ,r _{2} ,r,r′)η_{i}(r)η_{j}(r′)+. . . , (12)

[0047]
where
$\begin{array}{cc}\eta \ue8a0\left(r\right)=\left(\begin{array}{c}{\eta}_{1}\ue8a0\left(r\right)\\ {\eta}_{2}\ue8a0\left(r\right)\end{array}\right)=\left(\begin{array}{c}\delta \ue89e\text{\hspace{1em}}\ue89e\alpha \ue8a0\left(r\right)\\ \delta \ue89e\text{\hspace{1em}}\ue89eD\ue8a0\left(r\right)\end{array}\right),& \left(13\right)\end{array}$

[0048]
the action of the operator V has been taken into account and summation over repeated indices is implied with i, j=1, 2. The components of the operators K_{1 }and K_{2 }are given by

K _{1} ^{1}(r _{1} ,r _{2} ;r)=βG _{0}(r _{1} ,r)G _{0}(r,r _{2}) (14)

K _{1} ^{2}(r _{1} ,r _{2} ;r)=β∇_{r} G _{0}(r _{1} ,r)·∇_{r} G _{0}(r,r _{2}), (15)

K _{2} ^{11}(r _{1} ,r _{2} ;r,r′)=−βG _{0}(r′,r)G _{0}(r,r′)G _{0}(r′,r _{2}), (16)

K _{2} ^{12}(r _{1} ,r _{2} ;r,r′)=−βG _{0}(r′,r)∇_{r′} G _{0}(r,r′)·∇_{r′} G _{0}(r′,r _{2}), (17)

K _{2} ^{21}(r _{1} ,r _{2} ;r,r′)=−β∇_{r} G _{0}(r _{1} ,r)·∇_{r} G _{0}(r,r′)G _{0}(r′,r _{2}), (18)

K _{2} ^{22}(r _{1} ,r _{2} ;r,r′)=−β∇_{r} G _{0}(r _{1} ,r)·∇_{r}{∇_{r′} G _{0}(r,r′)·∇_{r′} G _{0}(r′,r _{2})} (19)

[0049]
The components of K
_{n }are given by
$\begin{array}{cc}\begin{array}{c}{K}^{{i}_{1\ue89en}\ue89e\dots \ue89e\text{\hspace{1em}}\ue89e{i}_{n}}\ue8a0\left({r}_{1},{r}_{2};{R}_{1},\dots \ue89e\text{\hspace{1em}},{R}_{n}\right)={\left(1\right)}^{n}\ue89e\sum _{{\alpha}_{1},\text{\hspace{1em}}\ue89e\dots \ue89e\text{\hspace{1em}},{\alpha}_{n}}^{\text{\hspace{1em}}}\ue89e\text{\hspace{1em}}\ue89e\frac{{\partial}^{{i}_{1}1}\ue89e{G}_{0}\ue8a0\left({r}_{1},{R}_{1}\right)}{\partial {R}_{1\ue89e{\alpha}_{1}}^{{i}_{1}1}}\ue89e\frac{{\partial}^{{i}_{1}\ue89e{i}_{2}2}\ue89e{G}_{0}\ue8a0\left({R}_{1},{R}_{2}\right)}{\partial {R}_{1\ue89e{\alpha}_{1}}^{{i}_{1}1}\ue89e\partial {R}_{2\ue89e{\alpha}_{2}}^{{i}_{2}1}}\times \cdots \\ \times \frac{{\partial}^{{i}_{n1}+{i}_{n}2}\ue89e{G}_{0}\ue8a0\left({R}_{n1},{R}_{n}\right)}{\partial {R}_{n1,{\alpha}_{n1}}^{{i}_{n1}1}\ue89e\partial {R}_{n\ue89e\text{\hspace{1em}}\ue89e{\alpha}_{n}}^{{i}_{n}1}}\ue89e\frac{{\partial}^{{i}_{n}1}\ue89e{G}_{0}\ue8a0\left({R}_{n},{r}_{2}\right)}{\partial {R}_{n\ue89e\text{\hspace{1em}}\ue89e{\alpha}_{n}}^{{i}_{n}1}},\end{array}& \left(20\right)\end{array}$

[0050]
where α_{1}, . . . , α_{n }label Cartesian components of the vectors R_{1}, . . . , R_{n }and no summation should be performed over indexes which are not explicitly present in the sum (i.e., for i_{k}=1).

[0051]
Observe that (12) is a functional power series expansion each term of which is multilinear in η. Thus we can expand Φ in tensor powers of

Φ=
K _{1} η+K _{2}η
η+. . . , (21)

[0052]
Here K
_{1 }is a linear operator which maps the Hilbert space H
_{1 }into the Hilbert space H
_{2 }and K
_{2 }may be interpreted as a tensor operator which maps H
_{1} H
_{1 }into H
_{2}.

[0053]
If the spatial fluctuations in α and D are sufficiently small, the series (21) may be truncated after its first term. This results in an effective linearization of the forward scattering problem with Φ=K_{1}η. The corresponding linear ISP has the solution η=K_{1} ^{+}Φ, where K_{1} ^{+}denotes the pseudoinverse of K_{1}. To construct the solution to the nonlinear ISP we act on (21) with the pseudoinverse operator K_{1} ^{+}and use the identity K_{1} ^{+}K_{1}=I_{H1}. We thus obtain

η=
K _{1} ^{+} Φ−K _{1} ^{+} K _{2}η
η+. . . . (22)

[0054]
Next, by iterating this result we find that

η=
K _{1} ^{+} Φ−K _{1} ^{+} K _{2} K _{1} ^{+} K _{1} ^{+}Φ
η+. . . , (23)

[0055]
which is a functional expansion for η in tensor powers of Φ. We will refer to (23) as the inverse scattering series for diffusion tomography.

[0056]
Several comments on the above result are necessary. First,

[0057]
(23) provides a formally exact solution to the inverse problem in diffusion tomography. It may be viewed as a nonlinear inversion formula whose first term coincides with the pseudoinverse solution to the linearized ISP. The higher order terms represent systematically improvable nonlinear corrections which, in principle, can be computed to arbitrarily high order. Thus, it is only necessary to solve the linear ISP in order to formally solve the nonlinear ISP. Second, (23) may also be obtained by formal inversion of the functional power series (9). This results in an explicit formula for the coefficient
_{n }of the nth term in (23):
$\begin{array}{cc}{\ue52a}_{n}=\left(\sum _{p=1}^{n1}\ue89e\text{\hspace{1em}}\ue89e{\ue52a}_{p}\ue89e\sum _{{i}_{1}+\text{\hspace{1em}}\ue89e\cdots \ue89e\text{\hspace{1em}}+{i}_{p}=n}^{\text{\hspace{1em}}}\ue89e\text{\hspace{1em}}\ue89e{K}_{{i}_{1}}\otimes \cdots \otimes {K}_{{i}_{p}}\right)\ue89e{\ue52a}_{1}\otimes \cdots \otimes {\ue52a}_{1},& \left(24\right)\end{array}$

[0058]
where
_{1}, =K
_{1} ^{+}. See Section 6.1 Third, it may be seen that the coefficients of all the higher order terms in (23) have K
_{1} ^{+} as a prefactor. As a result, to any finite order, the spatial resolution of images reconstructed using the nonlinear theory can never exceed the resolution of images reconstructed by linear means alone. Fourth, the shortrange propagation of diffusive waves implies that the forward scattering problem in diffusion tomography is weakly nonlinear. This is precisely the condition under which the inverse scattering series may be expected to exhibit fast convergence. Finally, the approach to the ISP based on (23) differs from Newtontype methods. This follows from the fact that such methods require the forward problem to be solved for each iteration.

[0059]
3. Nonlinear Inversion in the Plane Geometry

[0060]
The inverse scattering series was developed in a form which is independent of geometry. We now specialize to the case of the planar geometry. Other cases including the cylindrical and spherical geometries may also be considered.

[0061]
3.1 Inversion Formulas

[0062]
In the planar geometry measurements are taken on two parallel planes. Sources are taken to be located on the z=0 plane and detectors on the plane z=L. The object
105 to be imaged lies between the planes in the region
0<z<L. In this geometry, the unperturbed Green's function is given by the planewave decomposition
$\begin{array}{cc}{G}_{0}\ue8a0\left(r,{r}^{\prime}\right)=\int \frac{{d}^{2}\ue89eq}{{\left(2\ue89e\pi \right)}^{2}}\ue89eg\ue8a0\left(q;z,{z}^{\prime}\right)\ue89e\mathrm{exp}\ue8a0\left[\uf74e\ue89e\text{\hspace{1em}}\ue89eq\xb7\left(\rho {\rho}^{\prime}\right)\right],& \left(25\right)\end{array}$

[0063]
where we have used the notation r=(ρ, z). In the case of free boundaries, the function g(q; z, z′) is given by
$\begin{array}{cc}g\ue8a0\left(q;z,{z}^{\prime}\right)=\frac{\mathrm{exp}\ue8a0\left[Q\ue8a0\left(q\right)\ue89e\uf603z{z}^{\prime}\uf604\right]}{2\ue89eQ\ue8a0\left(q\right)\ue89e{D}_{0}}& \left(26\right)\end{array}$

[0064]
and in the case of boundary conditions of the type expressed by equation (2)
$\begin{array}{cc}g\ue89e\left(q;z,{z}^{\prime}\right)=\frac{l}{{D}_{0}}\ue89e\frac{\mathrm{sinh}\ue8a0\left[Q\ue8a0\left(q\right)\ue89e\left(L\uf603z{z}^{\prime}\uf604\right)\right]+Q\ue8a0\left(q\right)\ue89el\ue89e\text{\hspace{1em}}\ue89e\mathrm{cosh}\ue8a0\left[Q\ue8a0\left(q\right)\ue89e\left(L\uf603z{z}^{\prime}\uf604\right)\right]}{\mathrm{sinh}\ue8a0\left(Q\ue8a0\left(q\right)\ue89eL\right)+2\ue89eQ\ue8a0\left(q\right)\ue89el\ue89e\text{\hspace{1em}}\ue89e\mathrm{cosh}\ue8a0\left(Q\ue8a0\left(q\right)\ue89eL\right)+{\left(Q\ue8a0\left(q\right)\ue89el\right)}^{2}\ue89e\mathrm{sinh}\ue8a0\left(Q\ue89e\left(q\right)\ue89eL\right)},& \left(27\right)\end{array}$

[0065]
where

Q(q)≡(q ^{2} +k ^{2})^{1/2} (28)

[0066]
and we have assumed that either r or r′ lies on one of the measurement planes.

[0067]
We will find it advantageous to rewrite the inverse scattering series (23) in the form

η=η^{(1)}+η^{(2)}+ (29)

η^{(1)} =K _{1} ^{+}Φ (30)

η
_{(2)} =K _{1} ^{+} K _{2}η
^{(1)} η
^{(1)}, (31)

[0068]
where η^{(1) }is the solution to the linearized ISP and η^{(2) }is the first nonlinear correction. In the planar geometry, since the measurement planes have translational symmetry, it is natural to express (30) and (31) in the Fourier basis of twodimensional plane waves. In this representation (30) and (31) become

η^{(1)}(r)=∫d ^{2} q _{1} d ^{2} q _{2} K _{1} ^{+}(r,q _{1} ,q _{2})Φ(q _{1} ,q _{2}) (32)

η^{(2)}(r)=∫d ^{2} q _{1} d ^{2} q _{2} ∫d ^{3} r′d ^{3} r″K _{1} ^{+}(r;q _{1} ,q _{2})K _{2}(q _{1} ,q _{2} ;r′,r″)η^{(1)}(r′)η^{(1)}(r″) (33)

[0069]
Here

Φ(q _{1} ,q _{2})=∫d ^{2} p _{1} d ^{2} p _{2} exp{i(q _{1} ·p _{1} +q _{2} ·p _{2})}Φ(p _{1} ,z _{1} ,p _{2} ,z _{2}), (34)

K(q _{1} ,q _{2})=∫d ^{2} q _{1} d ^{2} q _{2} exp{i(q _{1} ·p _{1} +q _{2} ·p _{2})}K(p _{1} ,z _{1} ,p _{2} ,z _{2};·). (35)

[0070]
Note that according to (32) and (33), once K_{1} ^{+}(r; q_{1}, q_{2}) is determined the inverse problem is solved, in principle.

[0071]
3.2 Singular Value Decomposition of K_{1} ^{+}

[0072]
The SVD of the pseudoinverse operator K
_{1} ^{+}is given by
$\begin{array}{cc}{K}_{1}^{+}\ue8a0\left(r;{q}_{1},{q}_{2}\right)=\int \frac{1}{\sigma}\ue89e{f}_{\sigma}\ue8a0\left(r\right)\ue89e{g}_{\sigma}^{*}\ue8a0\left({q}_{1},{q}_{2}\right)\ue89e\uf74c\sigma ,& \left(36\right)\end{array}$

[0073]
where σ is the singular value associated with the singular functions f_{σ}and g_{σ}. The singular functions are eigenfunctions with eigenvalues σ_{2 }of the positive selfadjoint operators K_{1}*K_{1 }and K_{1}K*_{1}:

K _{1} ^{*} K _{1} f _{σ}=σ^{2} f _{σ}, (37)

K _{1} K _{1} ^{*} g _{σ}=σ^{2} g _{σ}. (38)

[0074]
In addition, the singular functions are related by

K _{1} f _{σ} =σg _{σ}, (39)

K _{1} ^{*} g _{σ} =σf _{σ}. (40)

[0075]
To proceed further, we require an explicit expression for K_{1 }(q_{1}, q_{2}; r). This is obtained by using the definitions of G_{0 }from Eq. (25) and K_{1 }from Eqs. (14,15) to carry out the Fourier transformation in Eq. (35):

K _{1}(q _{1} ,q _{2} ;r)=κ(q _{1} ,q _{2} ;z)exp{i(q _{1} +q _{2})·p}, (41)

[0076]
where the components of κ are given by
$\begin{array}{cc}\begin{array}{c}{\kappa}_{1}\ue8a0\left({q}_{1},{q}_{2};z\right)=\beta \ue89e\text{\hspace{1em}}\ue89eg\ue8a0\left({q}_{1};{z}_{1},z\right)\ue89eg\ue8a0\left({q}_{2};z,{z}_{2}\right),\\ \begin{array}{c}{\kappa}_{1}\ue8a0\left({q}_{1},{q}_{2};z\right)=\ue89e\beta (\frac{\partial g\ue8a0\left({q}_{1};{z}_{1},z\right)}{\partial z}\ue89e\frac{\partial g\ue8a0\left({q}_{2};z,{z}_{2}\right)}{\partial z}\\ \ue89e{q}_{1}^{\prime}\xb7{q}_{2}\ue89eg\ue8a0\left({q}_{1};{z}_{1},z\right)\ue89eg\ue8a0\left({q}_{2};z,{z}_{2}\right)).\end{array}\end{array}& \left(42\right)\end{array}$

[0077]
Using (41), we find that the matrix elements of the operator K
_{1},K*
_{1 }are given by
$\begin{array}{cc}\begin{array}{c}{K}_{1}\ue89e{K}_{1}^{*}\ue8a0\left({q}_{1},{q}_{2};{q}_{1}^{\prime},{q}_{2}^{\prime}\right)=\ue89e\int {\uf74c}^{2}\ue89eQ\ue89e\text{\hspace{1em}}\ue89e\delta \ue8a0\left(Q{q}_{1}{q}_{2}\right)\ue89e\delta \ue8a0\left(Q{q}_{1}^{\prime}{q}_{2}^{\prime}\right)\times \\ \ue89eM\ue8a0\left(\frac{1}{2}\ue89e\left({q}_{1}{q}_{2}\right),\frac{1}{2}\ue89e\left({q}_{1}^{\prime}{q}_{2}^{\prime}\right);Q\right),\end{array}& \left(44\right)\end{array}$

[0078]
where
$\begin{array}{cc}M\ue8a0\left(P,{P}^{\prime};Q\right)={\left(2\ue89e\pi \right)}^{2}\ue89e{\int}_{0}^{L}\ue89e\text{\hspace{1em}}\ue89e\uf74cz\ue89e\text{\hspace{1em}}\ue89ek\ue8a0\left(Q/2+P,Q/2P;z\right)\ue89e{k}^{*}\ue8a0\left(Q/2+{P}^{\prime},Q/2{P}^{\prime};z\right).& \left(45\right)\end{array}$

[0079]
To find the singular functions, we make the ansatz

g _{QQ′}(q _{1} , q _{2})∫d ^{2} PC _{Q′}(P; Q)δ(q _{1} −Q/2−P)δ(q _{2} −Q/2+P), (46)

[0080]
where Q and Q′ are twodimensional wavevectors. Eq. (38) now implies that

∫d ^{2} P′M(P, P′; Q)C _{Q′}(P′; Q)=σ_{QQ′} C _{Q′}(P; Q), (47)

[0081]
that is C
_{Q′}(P; Q) is an eigenfunction of M(P, P′; Q) labeled by Q′ with eigenvalue σ
_{QQ′} ^{2}. Note that since M(Q) is selfadjoint, the C
_{Q′}(P; Q) may be taken to be orthonormal. The singular functions f
_{QQ′} may be found from (40) by direct calculation:
$\begin{array}{cc}{f}_{{\mathrm{QQ}}^{\prime}}\ue8a0\left(r\right)=\frac{1}{{\sigma}_{{\mathrm{QQ}}^{\prime}}}\ue89e\int {\uf74c}^{2}\ue89eP\ue89e\text{\hspace{1em}}\ue89e\mathrm{exp}\ue8a0\left(\uf74e\ue89e\text{\hspace{1em}}\ue89eQ\xb7\rho \right)\ue89e{\kappa}^{*}\ue8a0\left(Q/2+P,Q/2P;z\right)\ue89e{C}_{{Q}^{\prime}}\ue8a0\left(P;Q\right).& \left(48\right)\end{array}$

[0082]
It follows that the SVD of K
_{1} ^{+}is given by the expression
$\begin{array}{cc}{K}_{1}^{+}\ue8a0\left(r;{q}_{1},{q}_{2}\right)=\int {\uf74c}^{2}\ue89eQ\ue89e\text{\hspace{1em}}\ue89e{\uf74c}^{2}\ue89e{Q}^{\prime}\ue89e\frac{1}{{\sigma}_{{\mathrm{QQ}}^{\prime}}}\ue89e{f}_{{\mathrm{QQ}}^{\prime}}\ue8a0\left(r\right)\ue89e{g}_{{\mathrm{QQ}}^{\prime}}^{*}\ue8a0\left({q}_{1},{q}_{2}\right).& \left(49\right)\end{array}$

[0083]
The above expression for the SVD of K
_{1} ^{+}may be simplified by using the spectral decomposition
$\begin{array}{cc}{M}^{1}\ue8a0\left(P,{P}^{\prime};Q\right)=\int {\uf74c}^{2}\ue89e{Q}^{\prime}\ue89e\frac{1}{{\sigma}_{{\mathrm{QQ}}^{\prime}}^{2}}\ue89e{C}_{{Q}^{\prime}}\ue8a0\left(P;Q\right)\ue89e{C}_{{Q}^{\prime}}^{*}\ue8a0\left({P}^{\prime};Q\right)& \left(50\right)\end{array}$

[0084]
and the explicit expressions for the singular functions. Eq. (49) thus becomes

K _{1} ^{+}(r;q _{1} ,q _{2})=∫d ^{2} Qd ^{2} P′exp(−iQ·p)M ^{−1}(P,P′;Q)×κ^{*}(Q/2+P,Q/2−P;z)δ(q _{1} −Q/2−P)δ(q _{2} −Q/2+P). (51)

[0085]
Using this result, along with (32), we obtain η^{(1)}(r), the solution to the linearized ISP:

η^{(1)}(r)=∫d ^{2} Qd ^{2} Pd ^{2} P′exp(−iQ·ρ)M ^{−1}(P,P′;Q)×κ^{*}(Q/2+P,Q/2−P;z)Φ(Q/2+P,Q/2−P). (52)

[0086]
Note that the above inversion formula is based on a direct calculation of the pseudoinverse solution rather than a construction of the SVD of the linearized forward scattering operator.

[0087]
4. Numerical Results

[0088]
We now illustrate the inversion formulas derived with numerical examples. We work in the planar geometry with free boundary conditions. In addition, we assume a priori that there are no inhomogeneities in the diffusion coefficient (δD=0). This allows the use of a single modulation frequency which we set to zero. In this case the linearized inversion formula (52) can be written in the form

δα^{(1)}(r)=∫d ^{2} Qd ^{2} Pd ^{2} P′exp(−iQ·ρ)M ^{−1}(P,P′;Q)×κ_{1} ^{*}(Q/2+P,Q/2−P;z)Φ(Q/2+P,Q/2−P), (53)

[0089]
where
$\begin{array}{cc}M\ue8a0\left(P,{P}^{\prime};Q\right)={\left(2\ue89e\text{\hspace{1em}}\ue89e\pi \right)}^{2}\ue89e{\int}_{0}^{L}\ue89e\uf74cz\ue89e\text{\hspace{1em}}\ue89e{\kappa}_{1}\ue8a0\left(Q/2+P,Q/2P;z\right)\ue89e{\kappa}_{1}^{*}\ue8a0\left(Q/2+{P}^{\prime},Q/2{P}^{\prime};z\right).& \left(54\right)\end{array}$

[0090]
In practice, this formula must be discretized. Namely, we chose the vectors Q to occupy a square lattice with unit step size Δq=k
_{1}=L
^{−1 }inside the circle Q≦40Δq. The vectors P, P′ were chosen on a onedimensional lattice coinciding with the xaxis; the spacing was Δp=5Δq=5ΔL
^{−1}, and 0≦P<40L
^{−1}. Thus a total of eight different vectors P were used (including P=0). Note that for numerical calculation of M
^{−1}, the operator M becomes a square matrix which can be diagonalized by the methods of linear algebra. In order to avoid numerical instability, the calculation of M
^{−1 }must be regularized. In particular, we replace 1/σ by R(σ) where R is a suitable regularizer. The effect of regularization is to limit the contribution of small singular values to the reconstruction. The simplest way to do this is to simply cutoff all σ below some cutoff σ
_{c}. That is, we set
$\begin{array}{cc}R\ue8a0\left(\sigma \right)=\frac{1}{\sigma}\ue89e\theta \ue8a0\left(\sigma {\sigma}_{c}\right),& \left(55\right)\end{array}$

[0091]
θdenoting the usual Heavyside step function.

[0092]
The forward data were calculated for a spherical absorbing inhomogeneity. The data function Φ(q_{1}, q_{2}) is given in the Section 6.2. The diffuse wave numbers are given by k_{1} ^{2}=α_{0}/D_{0 }outside the sphere and k_{2} ^{2}=(α_{0}+a)/D_{0 }inside the sphere. Thus, inside the sphere we have δα=a=D_{0}(k_{2} ^{2}−k_{1} ^{2}) and we have chosen D_{0}=1, k_{1}=L^{−1}. Reconstructions were carried out in the volume −L≦x, y≦L, 0<z<L with the center of absorbing sphere placed at the point (0, 0, L/2). The sphere's radius was taken to be R=0.4L and the corresponding size parameter of the sphere was k_{1}R=0.4.

[0093]
The results of linear reconstruction (δα^{(1)}) are shown in FIGS. 2A2D for k_{1}=0.99k_{2}, k_{1}=0.8k_{2}, k_{1}=0.8k_{2}, and k_{1}=0.5k_{2}, respectively. Generally, tomographic slices are drawn at different depths z for an absorbing sphere with different values of k_{2}. The solid black circles indicate the physical boundary of the absorbing sphere. The reconstructed images are normalized by the “true” value of δα, D_{0}(k_{2} ^{2}−k_{1} ^{2}). A linear gray scale is employed; white corresponds to 1 and black to 0. In the case k_{1}=0.99k_{2 }(FIG. 2A), the weak scattering approximation is quite accurate. As a result, the quality of reconstructed images is high even to lowest order in δα. As the mismatch of the absorption inside and outside the sphere becomes larger, the quality of the linearized inversion decreases. In particular, a false dark area in the center of the sphere develops. In the case k_{1}=0.5k_{2 }(δα/α_{0}=3) shown in FIG. 2D, only a thin outer shell is reconstructed, while the inside area of the sphere is almost completely black. This is explained by the fact that in the weak scattering approximation the inhomogeneities are probed by the unperturbed incident field. However, the field inside more absorbing areas, such as the absorbing sphere in this numerical example, differs from the incident field. In essence, the field does not penetrate into areas with very high absorption. The linearized reconstruction “interprets” this fact as the absence of inhomogeneity. Note that although linearized inversion is not accurate when δα/α_{0}>1, it still reconstructs the correct shape of an inhomogeneity. Thus, although the internal region of the absorbing sphere in the case k_{1}=0.5k_{2 }is not reconstructed, the overall spherical shape is reconstructed correctly. It is expected that this is a general property of the linearized reconstruction and is not limited to spherical shapes.

[0094]
We now consider the first nonlinear correction δα^{(2) }which is given by

δα^{(2)}(r)=∫d ^{2} Qd ^{2} Pd ^{2} P′exp(−i−Q·ρ)M ^{−1}(P,P′;Q)×κ*_{1}(Q/2+P,Q/2−P;z)Φ^{(1)}(Q/2+P,Q/2−P), (56)

[0095]
where

Φ^{(1)}(q _{1} ,q _{2})=∫d ^{3} rd ^{3} r′K _{2} ^{11}(q _{1} ,q _{2} ;r,r′)δα^{(1)}(r)δα^{(1)}(r′). (57)

[0096]
Here K
_{2} ^{11}(q
_{1}, q
_{2}; r
_{1}, r
_{2}) is obtained by Fourier transformation of (16):
$\begin{array}{cc}\begin{array}{c}{K}_{2}^{11}\ue8a0\left({q}_{1},{q}_{2};r,{r}^{\prime}\right)=\ue89e\frac{{G}_{0}\ue8a0\left(r,{r}^{\prime}\right)}{{\left(2\ue89e\text{\hspace{1em}}\ue89e\pi \ue89e\text{\hspace{1em}}\ue89e{D}_{0}\right)}^{2}\ue89eQ\ue8a0\left({q}_{1}\right)\ue89eQ\ue8a0\left({q}_{2}\right)}\ue89e\mathrm{exp}\ue8a0\left[\uf74e\ue8a0\left({q}_{1}\xb7{\rho}_{1}+{q}_{2}\xb7{\rho}_{2}\right)\right]\times \\ \ue89e\mathrm{exp}\ue8a0\left[Q\ue8a0\left({q}_{1}\right)\ue89e\uf603z{z}_{1}\uf604Q\ue8a0\left({q}_{2}\right)\ue89e\uf603z{z}_{2}\uf604\right].\end{array}& \left(58\right)\end{array}$

[0097]
The quantity Φ^{(1)}(q_{1}, q_{2}) is calculated by numerical integration using the precomputed function δα^{(1)}.

[0098]
The reconstructed absorption coefficient with the first nonlinear correction (δα=δα^{(1)}+δα^{(2)}) is shown in FIGS. 3A3D for the same k_{1 }values of FIGS. 2A2D As can be seen by comparing the panels with k_{1}=0.8k_{2 }and k_{1}=0.9k_{2 }in FIGS. 2B and 2C and FIGS. 3B and 3C, the effect of the first nonlinear correction is to fill in the voids that are seen in the linearized reconstruction. To illustrate this point more quantitatively, we plotted the reconstructed function δα, with and without the first nonlinear correction, as a function of the xcoordinate on the onedimensional line determined by z=const, y=0. The results are shown in FIG. 4. In particular, in this figure, the reconstructed function δα(r) is calculated along the line z=const (as indicated by the following legend for z) y=0, x ∈[−L, L]: (a) Long dash: δα=δα^{(1) }(linearized inversion); (b) solid line: δα=δα^{(1)}+δα^{(2) }(first nonlinear correction); and (c) short dash: the true profile of δα. The effect of filling in the void in the center of the sphere is especially well manifested in the case k_{1}=0.9k_{2}.

[0099]
As expected, the first nonlinear correction had no significant effect in the cases k_{1}=0.99k_{2 }and k_{1}=0.5k_{2 }(FIGS. 3A and 3C, respectively). In the first case, the linearized inversion already provides accurate results and all higherorder corrections are small. In the second case, the weak scattering approximation is strongly violated for the forward problem, and very high order corrections must be included to obtain convergence (provided the series converges at all).

[0100]
5. Flow Diagrams

[0101]
An embodiment illustrative of the methodology carried out by the subject matter of the present invention is set forth in highlevel flow diagram 500 of FIG. 5A in terms of the illustrative system embodiment shown in FIG. 1. With reference to FIG. 5A, the processing effected by block 510 enables source 120 and data acquisition detector 130 so as to measure the scattering data emanating from scatterer 105 due to illumination from source 120. These measurements are passed to computer processor 150 from data acquisition detector 130 via bus 131. Next, processing block 520 is invoked to formulate the nonlinear operator relating at least one coefficient characterizing an image of the scatterer/object to the measured scattering data (e.g., equation (9)). In turn, processing block 530 is operated to execute the nonlinear reconstruction algorithm using the nonlinear operator formulation. Finally, as depicted by processing block 540, the reconstructed tomographic image corresponding to a is provided to output device 170 in a form determined by the user; device 170 may be, for example, a display monitor or a more sophisticated threedimensional display device.

[0102]
Another embodiment illustrative of the methodology carried out by the subject matter of the present invention is set forth in highlevel flow diagram 550 of FIG. 5B in terms of the illustrative system embodiment shown in FIG. 1. With reference to FIG. 5B, the processing effected by block 560 enables source 120 and data acquisition detector 130 so as to measure the scattering data emanating from scatterer 105 due to illumination from source 120. These measurements are passed to computer processor 150 from data acquisition detector 130 via bus 131. Next, processing block 570 is invoked to compute the linear and tensor operators, that is, the operators given by equations (14) and (15) for the linear operator and equations (16)(19) for the tensor operator. In turn, processing block 580 is operated to execute the reconstruction algorithm set forth in equation (23), that is, the inverse scattering series, which is a functional expansion for η in tensor products of Φ, is computed. Finally, as depicted by processing block 590, the reconstructed tomographic image corresponding to η is provided to output device 170 in a form determined by the user; device 170 may be, for example, a display monitor or a more sophisticated threedimensional display device.

[0103]
6.1 Inversion of Series

[0104]
In this Section we show that the inverse scattering series (23) may be obtained by formal inversion of the forward scattering series

Φ=
K _{1} η+K _{2} ηη+K _{3}η
η
η
η+. . . . (59)

[0105]
To proceed, we assume that η may be expressed as a functional expansion in Φ:

η=
_{1}Φ+
_{2}Φ
Φ+
_{3}Φ
Φ
Φ+. . . , (60)

[0106]
where
_{1 }is a linear operator which maps the Hilbert space H
_{2 }into the Hilbert space H
_{1 }and
_{n }is a tensor operator which maps H
_{2 } . . .
H
_{2 }(n copies) into H
_{1 }for n≧2. To find the
's we substitute the expression (59) for Φ into (60) and equate terms with the same tensor power of η. We thus obtain the relations


_{2} K _{1} K _{1}+
_{1} K _{2}=0 (62)

_{3} K _{1} K _{1} K _{1}+
_{2} K _{1} K _{2}+
_{2} K _{2} K _{1}+
_{1} K _{3}=0 (63)

[0107]
[0107]
$\begin{array}{cc}\sum _{p=1}^{n1}\ue89e{p}_{}\ue89e\sum _{{i}_{1}+\cdots +{i}_{p}=n}\ue89e{K}_{{i}_{1}}\otimes \cdots \otimes {K}_{{i}_{p}}+{n}_{}\ue89e{K}_{1}\otimes \cdots \otimes {K}_{1}=0,& \left(64\right)\end{array}$

[0108]
which may be solved for the
's with the result


_{2}=−
_{1}K
_{2} _{1} _{1} (66)

_{3}=−(
_{2} K _{1} K _{2}+
_{2} K _{2} K _{1}+
_{1} K _{3})
_{1} _{1} _{1} (67)

[0109]
[0109]
$\begin{array}{cc}{n}_{}=\left(\sum _{p=1}^{n1}\ue89e{p}_{}\ue89e\sum _{{i}_{1}+\cdots +{i}_{p}=n}\ue89e{K}_{{i}_{1}}\otimes \cdots \otimes {K}_{{i}_{p}}\right)\ue89e{1}_{}\otimes \cdots \otimes {1}_{}.& \left(68\right)\end{array}$

[0110]
It can be seen that the above expressions for
_{1 }and
_{2 }agree with (23). In addition, (68) provides a general formula for the coefficients of the higher order terms in (
23). We note this formula implies that the nonlinear correction of order n involves all forward operators up to order n.

[0111]
6.2 The Data Function for a Spherical Inhomogeneity

[0112]
In this Section we calculate the data function for an absorbing spherical inhomogeneity in an infinite medium. Note that the scattering of diffusing waves from a sphere is analogous to Mie scattering in electromagnetic theory.

[0113]
Consider a spherical inclusion whose properties differ from the surrounding homogeneous background. We assume that δD=0 and δα=α=const inside a spherical region r−r
_{0}<R. We will work in a reference frame whose origin coincides with the center of the sphere, ro. In this case, the Green's function can be represented as
$\begin{array}{cc}G\ue8a0\left(r,{r}^{\prime}\right)=\sum _{l=0}^{\infty}\ue89e\sum _{m=l}^{l}\ue89e{g}_{l}\ue8a0\left(r,{r}^{\prime}\right)\ue89e{Y}_{l\ue89e\text{\hspace{1em}}\ue89em}^{*}\ue8a0\left({\hat{r}}^{\prime}\right)\ue89e{Y}_{l\ue89e\text{\hspace{1em}}\ue89em}\ue8a0\left(\hat{r}\right),& \left(69\right)\end{array}$

[0114]
where the Y
_{lm}'s are spherical harmonics. If both the sources and detectors are located outside the sphere (r, r′>R) then
$\begin{array}{cc}{g}_{l}\ue8a0\left(r,{r}^{\prime}\right)=\frac{2\ue89e{k}_{1}}{\pi \ue89e\text{\hspace{1em}}\ue89e{D}_{0}}\ue8a0\left[{i}_{l}\ue8a0\left({k}_{1}\ue89e{r}_{<}\right)\ue89e{k}_{l}\ue8a0\left({k}_{1}\ue89e{r}_{>}\right){F}_{l}\ue89e{k}_{l}\ue8a0\left({k}_{1}\ue89er\right)\ue89e{k}_{l}\ue8a0\left({k}_{1}\ue89e{r}^{\prime}\right)\right],& \left(70\right)\end{array}$

[0115]
where k
_{l }(x) and i
_{l}(x) are the modified spherical Bessel and Hankel functions of the first kind, r
_{<} and r
_{>} are the lesser and greater of r and r′, the Mie coefficient F
_{l }is given by
$\begin{array}{cc}{F}_{l}=\frac{{k}_{2}\ue89e{i}_{l}\ue8a0\left({k}_{1}\ue89eR\right)\ue89e{i}_{l}^{\prime}\ue8a0\left({k}_{2}\ue89eR\right){k}_{1}\ue89e{i}_{l}\ue8a0\left({k}_{2}\ue89eR\right)\ue89e{i}_{l}^{\prime}\ue8a0\left({k}_{1}\ue89eR\right)}{{k}_{2}\ue89e{i}_{l}^{\prime}\ue8a0\left({k}_{2}\ue89eR\right)\ue89e{k}_{l}^{\prime}\ue8a0\left({k}_{1}\ue89eR\right){k}_{1}\ue89e{i}_{l}\ue8a0\left({k}_{2}\ue89eR\right)\ue89e{k}_{l}^{\prime}\ue8a0\left({k}_{1}\ue89eR\right)},& \left(71\right)\end{array}$

[0116]
and k
_{1 }and k
_{2 }are the wavenumbers outside and inside the sphere:
$\begin{array}{cc}\begin{array}{cc}{k}_{1}^{2}=\frac{{\alpha}_{0}\uf74e\ue89e\text{\hspace{1em}}\ue89e\omega}{{D}_{0}},& {k}_{2}^{2}=\frac{{\alpha}_{0}+a\uf74e\ue89e\text{\hspace{1em}}\ue89e\omega}{{D}_{0}}.\end{array}& \left(72\right)\end{array}$

[0117]
By observing that in an infinite medium the unperturbed Green's function G
_{0}(r, r′) can be written as
$\begin{array}{cc}{G}_{0}\ue8a0\left(r,{r}^{\prime}\right)=\frac{2\ue89e{k}_{1}}{\pi \ue89e\text{\hspace{1em}}\ue89e{D}_{0}}\ue89e\sum _{l=0}^{\infty}\ue89e\sum _{m=l}^{l}\ue89e{i}_{l}\ue8a0\left({k}_{1}\ue89e{r}_{<}\right)\ue89e{k}_{l}\ue8a0\left({k}_{1}\ue89e{r}_{>}\right)\ue89e{Y}_{l\ue89e\text{\hspace{1em}}\ue89em}^{*}\ue8a0\left({\hat{r}}^{\prime}\right)\ue89e{Y}_{l\ue89e\text{\hspace{1em}}\ue89em}\ue8a0\left(\hat{r}\right),& \left(73\right)\end{array}$

[0118]
we see that the first term in (70) can be identified as the incident field, while the second term represents the scattered field. Consequently, the data function θ(r
_{1}, r
_{2}) is given by
$\begin{array}{cc}\phi \ue8a0\left({r}_{1},{r}_{2}\right)=\frac{2\ue89e{k}_{1}}{\pi \ue89e\text{\hspace{1em}}\ue89e{D}_{0}}\ue89e\sum _{l=0}^{\infty}\ue89e\sum _{m=l}^{l}\ue89e{F}_{l}\ue89e{k}_{l}\ue8a0\left({k}_{1}\ue89e{r}_{1}\right)\ue89e{k}_{l}\ue8a0\left({k}_{1}\ue89e{r}_{2}\right)\ue89e{Y}_{l\ue89e\text{\hspace{1em}}\ue89em}^{*}\ue8a0\left({\hat{r}}_{2}\right)\ue89e{Y}_{l\ue89e\text{\hspace{1em}}\ue89em}\ue8a0\left({\hat{r}}_{1}\right).& \left(74\right)\end{array}$

[0119]
The above expression is valid in a reference frame whose origin is at the center of the sphere. The corresponding expression in an arbitrary reference frame is obtained by making the transformation r_{1,2 }→r_{1,2}−r_{0}.

[0120]
We now calculate the Fourier transformed data function Φ(q_{1}, q_{2}) which is defined by

Φ(q _{1} ,q _{2})=∫d ^{2}ρ_{1} d ^{2}ρ_{2} exp[i(q _{1}·ρ_{1} +q _{2}·ρ_{2})]Φ(ρ_{1} ,z _{1},ρ_{2} ,z _{2}). (75)

[0121]
We will show that
$\begin{array}{cc}\phi \ue8a0\left({q}_{1},{q}_{2}\right)=\frac{{\pi}^{2}}{2\ue89e{D}_{0}\ue89e{k}_{1}\ue89eQ\ue8a0\left({q}_{1}\right)\ue89eQ\ue8a0\left({q}_{2}\right)}\ue89e\mathrm{exp}\ue8a0\left[i\ue8a0\left({q}_{2}{q}_{1}\right)\xb7{\rho}_{c}Q\ue8a0\left({q}_{1}\right)\ue89e\uf603{z}_{0}{z}_{1}\uf604Q\ue8a0\left({q}_{2}\right)\ue89e\uf603{z}_{2}{z}_{0}\uf604\right]\times \sum _{l=0}^{\infty}\ue89e\left(2\ue89el+1\right)\ue89e{F}_{l}\ue89e{P}_{l}\ue8a0\left(\frac{{q}_{1}\xb7{q}_{2}+\gamma \ue8a0\left({z}_{1},{z}_{2}\right)\ue89eQ\ue8a0\left({q}_{1}\right)\ue89eQ\ue8a0\left({q}_{2}\right)}{{k}_{1}^{2}}\right),& \left(76\right)\end{array}$

[0122]
where P
_{l}(x) are the Legendre polynomials and
$\begin{array}{cc}\gamma \ue8a0\left({z}_{1},{z}_{2}\right)=\{\begin{array}{cc}1,& \mathrm{if}\ue89e\text{\hspace{1em}}\ue89e{z}_{1}={z}_{2}\\ 1,& \mathrm{if}\ue89e\text{\hspace{1em}}\ue89e{z}_{1}\ne {z}_{2}.\end{array}& \left(77\right)\end{array}$

[0123]
Note that the arguments of the Legendre polynomials in (76) can be greater than unity; however the Mie coefficients F
_{l }decay with l faster than exponentially so that the convergence of the series is guaranteed. To derive (
76) we start by expanding Φ(r
_{1}, r
_{2}) in a threedimensional Fourier integral. To this end we define I
_{lm} ^{1}(p
_{1}) and I
_{lm} ^{2 }(p
_{2) }such that
$\begin{array}{cc}{k}_{l}\ue8a0\left({k}_{1}\ue89e{r}_{1}\right)\ue89e{Y}_{\mathrm{lm}}\ue8a0\left({\hat{r}}_{1}\right)=\int \frac{{d}^{3}\ue89e{p}_{1}}{{\left(2\ue89e\pi \right)}^{3}}\ue89e{I}_{\mathrm{lm}}^{1}\ue8a0\left({p}_{1}\right)\ue89e\mathrm{exp}\ue8a0\left({\mathrm{ip}}_{1}\xb7{r}_{1}\right),& \left(78\right)\\ {k}_{l}\ue8a0\left({k}_{1}\ue89e{r}_{2}\right)\ue89e{Y}_{\mathrm{lm}}^{*}\ue8a0\left({\hat{r}}_{2}\right)=\int \frac{{d}^{3}\ue89e{p}_{2}}{{\left(2\ue89e\pi \right)}^{3}}\ue89e{I}_{\mathrm{lm}}^{2}\ue8a0\left({p}_{2}\right)\ue89e\mathrm{exp}\ue8a0\left({\mathrm{ip}}_{2}\xb7{r}_{2}\right),& \left(79\right)\end{array}$

I _{lm} ^{1}(p _{1})=∫k _{l}(k _{1} r)Y _{lm}({circumflex over (r)})exp(ip _{1} ·r)d ^{3} r, (80)

I _{lm} ^{2}(p _{2})=∫k _{l}(k _{1} r)Y _{lm} ^{*}({circumflex over (r)})exp(ip _{2} ·r)d ^{3} r. (81)

[0124]
The integrals (80) and (81) are evaluated by using the identity
$\begin{array}{cc}\mathrm{exp}\ue8a0\left(\mathrm{ip}\xb7r\right)=4\ue89e\pi \ue89e\sum _{l=0}^{\infty}\ue89e\sum _{m=l}^{l}\ue89e{i}^{l}\ue89e{j}_{l}\ue8a0\left(\mathrm{pr}\right)\ue89e{Y}_{\mathrm{lm}}^{*}\ue8a0\left(\hat{p}\right)\ue89e{Y}_{\mathrm{lm}}\ue8a0\left(\hat{r}\right),& \left(82\right)\end{array}$

[0125]
where j
_{l}(x) are spherical Bessel functions of the first kind. After expanding the exponents in (80) and (81) according to (82) and integrating over the angular variables we obtain
$\begin{array}{cc}{I}_{\mathrm{lm}}^{1}\ue8a0\left({p}_{1}\right)=4\ue89e\pi \ue89e\text{\hspace{1em}}\ue89e{i}^{l}\ue89e{Y}_{\mathrm{lm}}\ue8a0\left({\hat{p}}_{1}\right)\ue89e{\int}_{0}^{\infty}\ue89e{r}^{2}\ue89e{k}_{l}\ue8a0\left({k}_{1}\ue89er\right)\ue89e{j}_{l}\ue8a0\left({p}_{1}\ue89er\right)\ue89e\uf74cr,& \left(83\right)\\ {I}_{\mathrm{lm}}^{2}\ue8a0\left({p}_{2}\right)=4\ue89e\pi \ue89e\text{\hspace{1em}}\ue89e{i}^{l}\ue89e{Y}_{\mathrm{lm}}^{*}\ue8a0\left({\hat{p}}_{2}\right)\ue89e{\int}_{0}^{\infty}\ue89e{r}^{2}\ue89e{k}_{l}\ue8a0\left({k}_{1}\ue89er\right)\ue89e{j}_{l}\ue8a0\left({p}_{2}\ue89er\right)\ue89e\uf74cr.& \left(84\right)\end{array}$

[0126]
The onedimensional integrals are easily calculated using the formula
$\begin{array}{cc}{\int}_{0}^{\infty}\ue89e{x}^{2}\ue89e{k}_{l}\ue8a0\left(\mathrm{ax}\right)\ue89e{j}_{l}\ue8a0\left({b}_{x}\right)\ue89e\uf74cx=\frac{\pi}{2\ue89ea}\ue89e\frac{{\left(b/a\right)}^{l}}{{a}^{2}+{b}^{2}}.& \left(85\right)\end{array}$

[0127]
Combining the above results, we can write the threedimensional Fourier expansion of Φ(r
_{1}, r
_{2}):
$\begin{array}{cc}\phi \ue8a0\left({r}_{1},{r}_{2}\right)=\frac{1}{{\left(2\ue89e\pi \right)}^{3}\ue89e{D}_{0}\ue89e{k}_{1}}\ue89e\sum _{\mathrm{lm}}^{\text{\hspace{1em}}}\ue89e{\left(1\right)}^{l}\ue89e{F}_{l}\ue89e\int {d}^{3}\ue89e{p}_{1}\ue89e{d}^{3}\ue89e{p}_{2}\ue89e\frac{{\left({p}_{1}\ue89e{p}_{2}/{k}_{1}^{2}\right)}^{l}}{\left({p}_{1}^{2}+{k}_{1}^{2}\right)\ue89e\left({p}_{2}^{2}+{k}_{1}^{2}\right)}\times {Y}_{\mathrm{lm}}\ue8a0\left({\hat{p}}_{1}\right)\ue89e{Y}_{\mathrm{lm}}^{*}\ue8a0\left({\hat{p}}_{2}\right)\ue89e\mathrm{exp}\ue8a0\left[i\ue8a0\left({p}_{1}\xb7\left({r}_{1}{r}_{0}\right)+{p}_{2}\xb7\left({r}_{2}{r}_{0}\right)\right)\right].& \left(86\right)\end{array}$

[0128]
where we have made the shift r_{1,2}→r_{1,2}−r_{0}. Thus the expression (86) is valid in an arbitrary reference frame

[0129]
Next we decompose the threedimensional vectors as P
_{1}=−q
_{1}+t
_{1}ê
_{z}, P
_{2}=q
_{2}+t
_{2}ê
_{z }and r
_{0}=ρ
_{0}+z
_{0}ê
_{z}. Taking into account that p
_{1,2} ^{2}=q
_{1,2} ^{2}+
_{1,2} ^{2}, this immediately leads to the following expression for the twodimensional Fourier transform of Φ(r
_{1}, r
_{2}):
$\begin{array}{cc}\phi \ue8a0\left({q}_{1},{q}_{2}\right)=\frac{2\ue89e\mathrm{\pi exp}\ue8a0\left[i\ue8a0\left({q}_{1}+{q}_{2}\right)\xb7{\rho}_{0}\right]}{{D}_{0}\ue89e{k}_{1}}\ue89e\sum _{\mathrm{lm}}^{\text{\hspace{1em}}}\ue89e{\left(1\right)}^{l}\ue89e{F}_{l}\ue89e{J}_{\mathrm{lm}}^{1}\ue8a0\left({q}_{1}\right)\ue89e{J}_{\mathrm{lm}}^{2}\ue8a0\left({q}_{2}\right),& \left(87\right)\end{array}$

[0130]
where
$\begin{array}{cc}{J}_{\mathrm{lm}}^{1}\ue8a0\left({q}_{1}\right)={\int}_{\infty}^{\infty}\ue89e\uf74ct\ue89e\frac{{\left(\sqrt{{q}_{1}^{2}+{t}^{2}}/{k}_{1}\right)}^{l}}{{q}_{1}^{2}+{t}^{2}+{k}_{1}^{2}}\ue89e{Y}_{\mathrm{lm}}\ue8a0\left(\frac{{q}_{1}+t\ue89e\text{\hspace{1em}}\ue89e{\hat{e}}_{z}}{\sqrt{{q}_{1}^{2}+{t}^{2}}}\right)\ue89e\mathrm{exp}\ue8a0\left[\mathrm{it}\ue8a0\left({z}_{0}{z}_{1}\right)\right]& \left(88\right)\\ {J}_{\mathrm{lm}}^{2}\ue8a0\left({q}_{2}\right)={\int}_{\infty}^{\infty}\ue89e\uf74ct\ue89e\frac{{\left(\sqrt{{q}_{2}^{2}+{t}^{2}}/{k}_{1}\right)}^{l}}{{q}_{2}^{2}+{t}^{2}+{k}_{1}^{2}}\ue89e{Y}_{\mathrm{lm}}^{*}\ue8a0\left(\frac{{q}_{2}+t\ue89e\text{\hspace{1em}}\ue89e{\hat{e}}_{z}}{\sqrt{{q}_{2}^{2}+{t}^{2}}}\right)\ue89e\mathrm{exp}\ue8a0\left[\mathrm{it}\ue8a0\left({z}_{0}{z}_{1}\right)\right].& \left(89\right)\end{array}$

[0131]
Although the integrands in (88) and (89) contain square roots, they are analytic functions of t. This can be seen by examining the explicit expressions for the spherical harmonics in terms of the associated Legendre polynomials and observing that the square roots in question are raised to an even power for any l and m. Therefore, (88) and (89) can be evaluated by analytic continuation of the integrands into the complex plane and contour integration. The result is
$\begin{array}{cc}{J}_{\mathrm{lm}}^{1}\ue8a0\left({q}_{1}\right)=\frac{\pi \ue89e\text{\hspace{1em}}\ue89e{i}^{l}}{Q\ue8a0\left({q}_{1}\right)}\ue89e{Y}_{\mathrm{lm}}\ue8a0\left(\frac{{q}_{1}+i\ue89e\text{\hspace{1em}}\ue89e\mathrm{sgn}\ue89e\text{\hspace{1em}}\ue89e\left({z}_{0}{z}_{1}\right)\ue89eQ\ue8a0\left({q}_{1}\right)\ue89e{\hat{e}}_{z}}{{\mathrm{ik}}_{1}}\right)\ue89e\mathrm{exp}\ue8a0\left[Q\ue8a0\left({q}_{1}\right)\ue89e\uf603{z}_{0}{z}_{1}\uf604\right],& \left(90\right)\\ {J}_{\mathrm{lm}}^{2}\ue8a0\left({q}_{2}\right)=\frac{\pi \ue89e\text{\hspace{1em}}\ue89e{i}^{l}}{Q\ue8a0\left({q}_{1}\right)}\ue89e{Y}_{\mathrm{lm}}^{*}\ue8a0\left(\frac{{q}_{2}+i\ue89e\text{\hspace{1em}}\ue89e\mathrm{sg}\ue89e\text{\hspace{1em}}\ue89en\ue89e\text{\hspace{1em}}\ue89e\left({z}_{0}{z}_{2}\right)\ue89eQ\ue8a0\left({q}_{2}\right)\ue89e{\hat{e}}_{z}}{{\mathrm{ik}}_{1}}\right)\ue89e\mathrm{exp}\ue8a0\left[Q\ue8a0\left({q}_{2}\right)\ue89e\uf603{z}_{0}{z}_{2}\uf604\right].& \left(91\right)\end{array}$

[0132]
Note that the spherical harmonic functions in the above expressions are analytically continued to complex angles; the arguments of Y_{lm }in (90) and (91) are complex unit vectors with the property a·a=1.

[0133]
Finally, we use the addition theorem to perform the summation over the index m in (87):
$\begin{array}{cc}\sum _{m=1}^{l}\ue89e{Y}_{\mathrm{lm}}\ue8a0\left(\frac{{q}_{1}+i\ue89e\text{\hspace{1em}}\ue89e\mathrm{sgn}\ue89e\text{\hspace{1em}}\ue89e\left({z}_{0}{z}_{1}\right)\ue89eQ\ue8a0\left({q}_{1}\right)\ue89e{\hat{e}}_{z}}{{\mathrm{ik}}_{1}}\right)\ue89e{Y}_{\mathrm{lm}}^{*}\ue8a0\left(\frac{{q}_{2}+i\ue89e\text{\hspace{1em}}\ue89e\mathrm{sg}\ue89e\text{\hspace{1em}}\ue89en\ue89e\text{\hspace{1em}}\ue89e\left({z}_{0}{z}_{2}\right)\ue89eQ\ue8a0\left({q}_{2}\right)\ue89e{\hat{e}}_{z}}{{\mathrm{ik}}_{1}}\right)=\frac{2\ue89el+1}{4\ue89e\pi}\ue89e{P}_{l}\ue8a0\left(\mathrm{cos}\ue89e\text{\hspace{1em}}\ue89e\theta \right),& \left(92\right)\end{array}$

[0134]
where θ is the angle between the two complex vector arguments of the spherical harmonic functions in (92). The cosine of this angle is obviously given by
$\begin{array}{cc}\mathrm{cos}\ue89e\text{\hspace{1em}}\ue89e\theta =\frac{{q}_{1}\xb7{q}_{2}+\mathrm{sgn}\ue8a0\left({z}_{0}{z}_{1}\right)\ue89e\text{\hspace{1em}}\ue89e\mathrm{sgn}\ue8a0\left({z}_{0}{z}_{2}\right)\ue89eQ\ue8a0\left({q}_{1}\right)\ue89eQ\ue8a0\left({q}_{2}\right)}{{k}_{1}^{2}}.& \left(93\right)\end{array}$

[0135]
Taking into account that sgn(z_{0}−z_{1})sgn(z_{0}−Z_{2})=γ(z_{1}, Z_{2}), we obtain the result (76).

[0136]
Although the present invention has been shown and described in detail herein, those skilled in the art can readily devise many other varied embodiments that still incorporate these teachings. Thus, the previous description merely illustrates the principles of the invention. It will thus be appreciated that those with ordinary skill in the art will be able to devise various arrangements which, although not explicitly described or shown herein, embody principles of the invention and are included within its spirit and scope. Furthermore, all examples and conditional language recited herein are principally intended expressly to be only for pedagogical purposes to aid the reader in understanding the principles of the invention and the concepts contributed by the inventor to furthering the art, and are to be construed as being without limitation to such specifically recited examples and conditions. Moreover, all statements herein reciting principles, aspects, and embodiments of the invention, as well as specific examples thereof, are intended to encompass both structural and functional equivalents thereof. Additionally, it is intended that such equivalents include both currently know equivalents as well as equivalents developed in the future, that is, any elements developed that perform the function, regardless of structure.

[0137]
In addition, it will be appreciated by those with ordinary skill in the art that the block diagrams herein represent conceptual views of illustrative circuitry embodying the principles of the invention.