Imaging Spatially Varying Dynamic Media With Diffusing Correlation Waves

This invention relates generally to imaging. More specifically, this invention relates to using diffusing light fluctuations to image media with spatially varying dynamic and optical intensity fluctuations.

Imaging of the human body using visible light has long been a desirable goal. In the late 1970s, dynamic light scattering theory was applied to living tissue to measure blood flow. Multiple scattering from the blood occurred, resulting in a Doppler broadening of the laser line width as revealed by the intensity fluctuations of the remitted speckles. This type of imaging worked fairly well, since the body behaved as a static matrix, and there was only occasional scattering from moving particles within the matrix. In about 1981, Bonner and Nossal at the National Institute of Health, wrote down the theory for imaging of blood flow using light, thereby formalizing the theory of "diffusing wave spectroscopy".

By the mid-1980s, diffusing wave spectroscopy was used to image optically dense systems in which all particles are moving, but the medium is generally viewed as homogeneous; that is, there are no spatial variations in the dynamic or optical properties. The aforementioned techniques for imaging with diffuse light are therefore limited since they cannot characterize media with spatially varying, dynamic properties, but can only characterize average properties of the media with well-defined boundaries.

Therefore, the imaging art has not provided a technique for imaging heterogeneous media. The aforementioned prior imaging techniques have only been able to determine dynamic and optical characteristics of a homogeneous medium, such as a dense colloid, wherein the magnitude of fluctuation of the intensities is obtained from the temporal intensity, or the correlation from an emerging speckle.

The current theory and devices for analyzing homogeneous media do not suffice for analyzing heterogeneous fluctuating media which are analogous to tumors, burns, and other real world structures found in the human body. There thus exists a long-felt but unsatisfied need in the art for a method of imaging spatially varying dynamic media which has not heretofore been satisfied.

The aforementioned problems are solved and long-felt needs met by the present invention which in a preferred embodiment provides methods of imaging a medium having spatially varying dynamic properties or spatially varying optical properties. The methods preferably comprise the steps of irradiating the medium with a source of light which will diffuse through the medium, measuring at a detector position the temporal intensity fluctuations of photon streams that have been scattered from objects within the medium, determining the properties of the medium using measured power spectrums or temporal correlation functions, and reconstructing an image of the medium from the intensity fluctuations.

Systems for imaging a medium having spatially varying dynamic properties or spatially varying optical properties also solve the above-referenced long-felt needs. Preferably, the systems comprise a radiation source that irradiates the medium with light and for setting up a diffuse wave of photons which will scatter through the medium, a detector placed at a known position with respect to the radiation source for detecting light scattered through the medium, an analyzer interfaced to the detector for quantifying the intensity fluctuations received by the detector, thereby measuring a diffuse correlation wave, and a processor interfaced with the analyzing means for determining the properties of the medium and reconstructing an image of the medium.

The methods and systems provided in accordance with the present invention provide a significant advance in the imaging art. With the methods and systems described and claimed herein, it is possible to provide low resolution images of spatially varying, dynamic turbid media such as the human body. Such images of turbid media with spatially varying dynamic properties or spatially varying optical properties have also not been achieved in the imaging art prior to the invention of the systems and methods described herein.

The invention will be best understood by those with skill in the art by reading the foregoing detailed description of preferred embodiments of the invention in conjunction with the drawings which are first described briefly below.

FIG. 1 is a block diagram of a system for imaging turbid media with spatially varying dynamic properties or spatially varying optical properties provided in accordance with the present invention.

FIG. 2 is a graph of the experimental measurements of a temporal field or autocorrelation function verifying that a diffusing autocorrelation occurs.

FIGS. 3A and 3B are, respectively, an example of source detector orientation for imaging a medium having spatially varying dynamic properties or spatially varying optical properties, and a model system showing the image of a spatially dynamic medium provided in accordance with the present invention.

Referring now to the drawings wherein like reference numerals refer to like elements, FIG. 1 is a block diagram of a system provided in accordance with the present invention for imaging media having spatially varying dynamic properties or spatially varying optical properties.

A turbid medium is shown generally at 10, and comprises a substantially heterogeneous matrix 20 wherein there is embedded an object 30. The object 30 contains small particles which are moving in a non-ordered fashion. "Non-ordered fashion" means that the particles exhibit Brownian motion, turbulent motion, shear flow, random motion or any other motion which results in a change in the relative distance between the particles. This is what is meant by the term "dynamic" throughout. In an experimental setup which has been used to verify the systems of the present invention, the matrix 20 is a solid slab of titanium dioxide suspended in resin, and the object 30 is a small cavity with a 0.2 percent suspension of 0.296 micrometer diameter polystyrene spheres. The slab has dimensions of 15.times.15.times.8 centimeters.

A source of light 40 is coupled by a multimode fiber 50 to the medium 10. Preferably, the source 40 is a coherent source of energy. Even more preferably, the source is a laser. Most preferably, the source is an argon ion laser which outputs energy in the 514 nanometer green wavelength range.

A fiber 60 is placed at a known position with respect to the matrix 20, and will pick up light which diffuses through the medium 10. A detector 70 is interfaced to the fiber, and by standard gain techniques creates a signal which is representative of the intensity fluctuations at the fiber resulting from the photons which have diffused through the medium 10 and which may have scattered from the particles in object 30. Any detector which can produce a gain, for example a photomultiplier tube, can be interfaced with the fiber 60, and it will be recognized by those with skill in the art that a fiber-detector combination will produce the results required in accordance with the invention. A preferred embodiment of fiber 60 is a single mode fiber.

In accordance with the present invention, a digital autocorrelator 90 is interfaced with the photomultiplier tube 70 in order to observe the intensity fluctuations of the signal speckle. The digital autocorrelation device 90, which is a well-known electronic system that is commercially available, measures the temporal intensity autocorrelation function of the photons received by the detector. As will be described in more detail below, this autocorrelation function can be used to obtain the scattering, absorption and dynamic characteristics of the medium in accordance with the methodology described herein.

In accordance with the invention, autocorrelation functions are measured with the source and collecting fibers individually positioned at different locations with respect to the object 30. A computer processor 100 with the appropriate software that implements the correlation diffusion theory described below determines the scattering, absorption and dynamic characteristics of the medium from the diffusion correlation wave and thereby reconstructs an image of the medium 10. The computer can be any known, standard processor which can utilize the correlation information output by the autocorrelator 90. In this manner, a reconstructed image of the medium 10 having the object 30 therein can be produced as a function of the scattering and absorption of the diffuse correlation wave as it propagates diffusing through the medium 10.

It has been demonstrated that transmissions and remission measurements of diffuse light intensity emerging from heterogeneous turbid media can provide adequate information for the construction of low resolution images of the spatially varying absorption and scattering coefficients. When the medium is dynamic, the time-dependent density fluctuations of the sample are impressed upon the temporal behavior of the diffusing light. For homogeneous fluctuating turbid media, the temporal intensity autocorrelatiLon function of an emerging speckle of diffuse light can be analyzed to provide information about the temporal fluctuations within the sample. This technique has been dubbed "diffusing wave spectroscopy". Just as diffuse photon density waves can be used to probe the spatial variations in absorption and scattering coefficients within heterogeneous media, the position-dependent measurements of the temporal autocorrelation function of diffusing light fields can be used to generate images of temporal fluctuations within heterogeneous turbid media. The techniques of the present invention provide a new contrast mechanism for imaging within turbid media. Specifically, these techniques can distinguish regions of strong and weak density fluctuations such as calcium deposits versus soft tissue or blood flow versus stasis.

The field and intensity temporal autocorrelation functions of coherent light emanating from a scattering medium provide information about the dynamic properties and/or the optical properties of the medium. Alternatively, it is well known that a power spectrum analysis of the light fluctuations can also be applied to characterize the dynamic properties and/or optical properties. Either embodiment may be applied in accordance with the invention.

The normalized temporal field autocorrelation function of light that has diffused through a homogeneous, turbid colloid is given by: ##EQU1## where E is the electric field of the emerging speckle. P(s) is the probability of a photon traveling a length s between the source and detector. D.sub.s is the diffusion coefficient of the scattering particles. k.sub.o is the wavenumber of the photon in the medium, and .mu..sub.s ', is the reduced scattering factor for the diffusing photon. The brackets denote an ensemble average which, for an ergodic system is the same as a temporal average. This form for the temporal field autocorrelation function is accurate for homogeneous, fluctuating turbid media, such as a colloid. Here, homogeneous refers to the absorption, scattering, and dynamic properties of the medium.

The steady state correlation transport equation is:

Here, G.sub.1 (r,.OMEGA.,.tau.) is the unnormalized temporal field-autocorrelation function (i.e. G.sub.1 (r,.OMEGA.,.tau.)=<E(0)E*(.tau.)>) which is a function of position r. direction .OMEGA., and correlation time .tau.. The scattering and absorption coefficients are respectively .mu..sub.s and .mu..sub.a and .mu..sub.t =.mu..sub.s +.mu..sub.a. Furthermore, g.sub.1.sup.s (.OMEGA.,.OMEGA.',.tau.) is the normalized temporal field-autocorrelation function for single scattering, given by:

and f(.OMEGA.,.OMEGA.') is the phase function. Finally, v is the speed of light in the medium and S(r..OMEGA.)is the distribution of light sources. Note that when .tau.=0 the steady state correlation transport equation is reduced to the steady state photon transport equation.

Using the P.sub.N approximation or spherical Harmonics method, a diffusion equation for correlation is derived from this more general transport equation. The correlation diffusion equation is as follows:

where ##EQU2## is the photon diffusion coefficient. Notice that in the correlation diffusion equation, G.sub.1 (r..tau.) is no longer a function of .OMEGA., and the reduced scattering coefficient .mu.'.sub.s appears instead of the scattering coefficient .mu..sub.s. For a homogeneous system, the solution of the correlation equation is identical to the solution of Eq. (1). The homogeneous solution for correlation diffusion formally has also been derived from the scalar wave equation for the electric field and the angular-averaged dynamic structure factor. The correlation diffusion equation presented here, however, is particularly useful when considering solutions for media having spatially varying dynamic properties or spatially varying optical properties.

Given a piecewise homogeneous system, a solution to the diffusion equation is found by satisfying the differential equation in each homogeneous region while matching the boundary conditions at the surface between regions. An analytic solution exists for a spherical inhomogeneity embedded in an otherwise homogeneous system.

To simplify the bookkeeping, Eq. (4) can be re-written as

where K.sup.2 (.tau.)=-3 (.mu.'.sub.3 +.mu..sub.a) (.mu..sub.a +2.mu..sub.s D.sub.s k.sub.o.sup.2 .tau.) and it has been assumed that a point source in space at position r.sub.s. Thus, 2.mu.'.sub.s D.sub.s k.sub.o.sup.2 .tau. represents the absorption of correlation due to dynamic processes. When .tau.=0 there is no dynamic absorption and eq. (5) is simply reduced to the photon diffusion equation.

The solution of eq. (5) outside of a spherical inhomogeneity is ##EQU3## where the sphere is centered on the origin and the source is placed on the z-axis. Here, h.sub.1 .sup.(1) are the Hankel functions of the first kind and Y.sub.1.sup.0 (.theta.,.phi.) are the spherical harmonics. The coefficient A.sub.1 is found by matching the appropriate boundary conditions on the surface of the sphere. These boundary conditions are the same as for photon diffusion namely G.sub.1.sup.out (r,.tau.)=G.sub.1.sup.in (r..tau.) and -D.sub.p.sup.out n..gradient.G.sub.1.sup.out (r,.tau.)=-D.sub.p.sup.in n..gradient.G.sub.1.sup.in (r..tau.) on the surface of the sphere, where G.sub.1.sup.in (r..tau.) is the correlation function inside the spherical object and n is the normal vector to the sphere. Applying these boundary conditions: ##EQU4## where j.sub.1 are the spherical Bessel functions of the first kind x=K.sup.out (.tau.)a. y=K.sup.in (.tau.)a, a is the radius of the sphere r.sub.s =(r=Z.sub.s .theta.=.pi.,.theta.=0), and j.sub.1 ' and h.sub.1.sup.(1)' are the first derivatives of the functions j.sub.t and h.sub.1.sup.(1) with respect to the argument.

For systems with spatially varying dynamic properties, the correlation diffusion equation is approximated as:

where K.sup.2 (r. .tau.)=-3 (.mu.'.sub.s +.mu..sub.a) (.mu..sub.a +2.mu.'.sub.s D.sub.s (r)k.sub.o.sup.2 .tau.) and D.sub.s (r)=D.sub.s.sup.avg +.delta.D.sub.s (r). The average particle self-diffusion coefficient of the systems is given by D.sub.s.sup.avg, and .delta. D.sub.s (r) are the variations with respect to this average. In this form, K.sup.2 accounts for dynamics governed by Brownian motion. K.sup.2 is easily modified to account for other types of dynamics as defined before.

Assuming a solution of the form G.sub.1 (r..tau.)=U.sub.o (r..tau.)exp (.PHI.(r. .tau.)), where U.sub.o (r. .tau.) is the solution of the homogeneous correlation diffusion equation, and apply the Rytov approximation, although other approximations may also be implemented, the solution of the heterogeneous diffusion equation is:

where G(r',r,.tau.) is the Green's function for the homogeneous diffusion equation. The perturbations due to any heterogeneity are thus contained in .PHI.(r,.tau.). Measurements of the perturbed temporal field autocorrelation function using several different source-detector pairs can be used along with knowledge of the optical properties and average particle self-diffusion coefficient of the medium to determine .delta. D.sub.s (r). The reconstruction is facilitated by inversion of eq. (9).

Referring now to FIG. 2, the results of experiments show the perturbation of diffusing correlation by fluctuating spherical objects. In particular, the intensity autocorrelation function for a semi-infinite highly-scattered solid containing a spherical cavity filled with a highly-scattered colloid was measured. This simulates the human body.

The experimental apparatus is shown in FIG. 1. The light source which was the 514 nanometer line of an argon ion laser (2.5 Watts and operated with an etalon to increase the source coherence). As discussed above, the beam is coupled into a multi-mode fiber for flexibility in positioning the source on the sample. Light is collected from the sample using a single mode fiber in order to observe the intensity fluctuation of single speckles/modes. The collected light is delivered to a photon-counting photo-multiplier tube, and finally, a digital autocorrelator is used to determine the temporal intensity autocorrelation function. Autocorrelation functions are measured with the source and collecting fibers individually positioned at different locations with respect to the sphere.

The autocorrelation functions are measured at different locations with respect to the sphere by individually positioning the source and collecting fibers at different locations within the sample. Note that the temporal intensity autocorrelation function is related to the temporal field autocorrelation function by the Seigert relation: ##EQU5##

In FIG. 2, the measured decay of the field autocorrelation function (symbols) compared to the theoretical prediction (lines) was plotted. The correlation function for the measurement made near the fluctuating sphere decays more rapidly than when measured further from the sphere. This is expected since in the former case a large fraction of the detected photons have sampled the liquid. The difference observed between the measurements at X=+1.0 cm indicates that the sphere is not located at X=0. The theoretical best fit to the experimental data using three free parameters produced good agreement between experiment and theory and reasonable values for the free parameters. The tree parameters and the fits were: X.sub.o =0.08 cm and the optical properties of the solid .mu..sub.s '=2.0 cm.sup.-1 and .mu..sub.a =0.08 cm.sup.-1. The sphere was not perfectly centered at X.sub.o =0.0 cm as observed in the experimental data and the fit.

In order to reconstruct an image with a spatially varying particle diffusion coefficient, it is necessary to measure the temporal intensity autocorrelation function with several source-detector configurations so as to sample all volume elements within the sample. To generate the necessary data for the image reconstruction, the analytic solution is used for a spherical inhomogeneity embedded in an otherwise homogeneous system. The system is an infinite medium with homogeneous optical properties, .mu..sub.s 40 =10.0 cm.sup.-1 .mu..sub.a =0.05 cm.sup.-1 and and D.sub.s =0. A 0.5 cm diameter sphere is centered at x=0.5 cm. y=0.0 cm. and z=0.0 cm. Two different values for the object D.sub.s were used 1.0.times.10.sup.-8 cm.sup.2 s.sup.-1 and 5.0.times.10.sup.-9 cm.sup.2 s.sup.-.

Referring to FIG. 3A, the system is illustrated to indicate the various source-detector positions relative to the spherical inhomogeneity. For each source-detector pair, the correlation function was calculated for correlation times ranging from 40 .mu.s to 200 .mu.s in steps of 40 .mu.s. In FIG. 3B, The reconstructed images for the two different systems are presented. The spherical object is clearly distinguished in each case, from the grey scale we see the expected factor of two difference in the self-diffusion coefficient of the scattering particles in the sphere.

Measurements were made with a single source and detector separated by 170.degree.. The source and detector were displaced 2.3 cm from the center of the circle and were rotated as a unit about the center of the circle, in intervals of 10.degree.. The lightly shaded circle indicates the region reconstructed in the images presented in (b)-(c) while the dark circle indicates the actual position of the inhomogeneity. The inhomogeneity is positioned 0.5 cm laterally from the center and has a diameter of 0.5 cm. Images reconstructed with the two different values of D.sub.s for the object are presented in (b)-(c). The object D.sub.s equals 1.0.times.10.sup.-8 cm.sup.2 s.sup.-1 and 5.0.times.10.sup.-9 cm.sup.2 s.sup.-1 in (b) and (c) respectively. The reconstructed values for D.sub.s are 1.0.times.10.sup.-8 cm.sup.2 s.sup.-1 and 5.0.times.10.sup.-9 cm.sup.2 s.sup.-1 respectively normalized to indicate that the relative values were accurately reconstructed.

Diffusion photons are sensitive to the random motion of the scattering particles. This information is accessible through correlation measurements of the intensity fluctuations of a single speckle of remitted or transmitted light. Thus, the correlation diffusion equation accurately predicts the decay of the temporal intensity autocorrelation function for heterogeneous turbid media. Furthermore, it has been demonstrated that the correlation diffusion equation can be used in an imaging algorithm to reconstruct the dynamic heterogeneities of a turbid media. Sensitivity to dynamic heterogeneities can be exploited to gain contrast in identifying calcified tumors and skeletal structures, as well as distinguishing regions of the body which have varying blood flow, or which are ischemic or in stasis.

There have thus been described certain preferred embodiments of systems and methods for imaging turbid media having spatially varying dynamic properties or spatially varying optical properties provided in accordance with the present invention. While preferred embodiments have been described and disclosed, it will be recognized by those with skill in the art that modifications are within the true spirit and scope of the invention. The appended claims are intended to cover all such modifications.