Imaging transverse flow velocity using spectral bandwidth of the doppler frequency shift in phase-resolved optical doppler tomography

[0001] The present application is related to U.S. Provisional Patent Application Ser. No. 60/363,580, filed on Mar. 12, 2002, which is incorporated herein by reference and to which priority is claimed pursuant to 35 USC 119.

[0002] This invention was made with Government support under Grant No. HL84218 awarded by the National Institute of Health and under F49620-00-1-0371 awarded by the U.S. Air Force. The U.S. Government has certain rights in this invention.

[0003] 1. Field of the Invention

[0004] The invention relates to an improvement in imaging transverse flow velocities using spectral bandwidth of the Doppler frequency shift in phase-resolved optical Doppler tomography.

[0005] 2. Description of the Prior Art

[0006] Optical Doppler tomography (ODT) combines Doppler velocimetry with optical coherence tomography (OCT) for noninvasive location and measurement of particle flow velocity in highly scattering media with micrometer-scale spatial resolution. The principle employed in ODT is very similar to that used in radar, sonar and medical ultrasound. ODT uses a low coherence source and optical interferometer to obtain high spatial resolution gating with a high speed scanning device such as a conventional rapid scanning optical delay line (RSOD) to perform fast ranging of microstructure and particle motion detection in biological tissues or other turbid media.

[0007] To detect the Doppler frequency shift signal induced by the moving particles, several algorithms and hardware schemes have been developed for ODT. The most straightforward method to determine the frequency shift involves the use of a small time fast Fourier transform (STFFT). However, the sensitivity of this method is mainly dependent on the FFT time window, which limits axial scanning speed and spatial resolution when measuring slowly moving blood flow in small vessels that requires high velocity sensitivity. However, a phase-resolved technique can decouple the Doppler sensitivity and spatial resolution while maintaining high axial scanning speed.

[0008] The limitation of phase-resolved ODT is that the dynamic range of the Doppler frequency shift measured is small because of an aliasing phenomenon caused by 2π ambiguity in the arctangent function. In Y. Zhao et. al., Opt. Lett. 25, 1358 (2000) on phase-resolved ODT, we reported on the use of the standard deviation of the Doppler spectrum to locate the microvasculature.

[0009] The invention is an improvement in an optical Doppler interferometer such as described in Z. Chen, et. al. U.S. patent application Ser. No. 09/574,560 filed May 19, 2000 for “Phase-Resolved Optical Coherence Tomography and Optical Doppler Tomography for Imaging Fluid Flow in Tissue with Fast Scanning and High Velocity Sensitivity” issued as U.S. Pat. No. 6,549,801 (2003),which is incorporated herein by reference. The Doppler bandwidth extracted from the standard deviation of the frequency shift in phase-resolved functional optical coherence tomography (F-OCT) is used to image the velocity component transverse to the optical probing beam. It is found that above certain threshold level the Doppler bandwidth is a linear function of transverse flow velocity and the effective numerical aperture of the optical objective in the sample arm determines the slope of this dependence.

[0010] The Doppler bandwidth allows accurate measurement of flow velocity without the precise determination of flow direction when the Doppler flow angle is within ±15 degrees perpendicular to the probing beam. Such an approach extends the dynamic range of flow velocity measurements obtained using phase-resolved F-OCT.

[0011] More particularly the invention is defined as an improvement in a method for tomographic imaging of a fluid flow using phase-resolved functional optical coherence tomography (F-OCT) comprising the steps of extracting a Doppler bandwidth above a certain threshold level from the standard deviation of the frequency shift as a linear function of flow velocity by use of a optical probing beam, and imaging the flow velocity component transverse to the optical probing beam.

[0012] An interferometer with a sample arm having an optical objective with an effective numerical aperture is used to extract the Doppler bandwidth. The step of extracting the Doppler bandwidth as a linear function of flow velocity extracts the Doppler bandwidth as a linear function with a slope determined by the effective numerical aperture of the optical objective in the sample arm.

[0013] The optical probing beam extracts the Doppler bandwidth from the standard deviation and the step of extracting the Doppler bandwidth above a certain threshold level from the standard deviation of the frequency shift is performed without the precise determination of flow direction when the Doppler flow angle is within ±15 degrees perpendicular to the probing beam.

[0014] The step of extracting a Doppler bandwidth as a linear function of flow velocity is performed as defined by $B_d=\frac{4V\sin \theta {\mathrm{NA}}_{\mathrm{eff}}}{\lambda }$
where B_d is the Doppler bandwidth, V is the flow velocity, λ is the wavelength of the probing beam, θ is the Doppler angle and, NA_eff is the effective numerical aperture and the Doppler bandwidth.

[0015] In the embodiment where the optical probing beam is a Gaussian beam, the standard deviation, σ, of the frequency shift is performed as defined by $\sigma =\frac{\pi V\sin \theta N{A}_{e\mathrm{ff}}}{8\lambda }.$
The standard deviation, σ, of the frequency shift is then performed as defined by $\sigma =\frac{\pi V\sin \theta N{A}_{e\mathrm{ff}}}{8\lambda }+b$
where b accounts for bandwidth spectrum broadening from contributions from Brownian motion and other factors that are independent of the macroscopic flow velocity.

[0016] The invention is further defined as an improvement in a method of using phase-resolved functional optical coherence tomography (F-OCT) comprising the steps of: measuring in an interferometer complex interference fringe signals {tilde over (Γ)}_ODT(i,j) for each pixel; generating for each pixel a variance matrix σ_i,j² based on normalized variance values; generating an intensity map I_i,j reflecting the axial distribution of backscattered intensities based on the amplitude, A(T), obtained by using the complex fringe signals; averaging the I_i,j over a predetermined set of scanning lines to obtain an axial intensity profile I_i; averaging temporally over the set, N, of axial scanning lines to find an axial variance profile σ_i², in which the variance values are weighted with respect to the intensity map I_i,j ${\sigma }_i^2=\frac{\sum _{j=1}^N{\sigma }_{i,j}^2{I}_{i,j}}{\sum _{j=1}^N{I}_{i,j}};$
generating a normalized standard deviation profile σ_i as an expression of the standard deviation of the Doppler power spectrum as a function of depth; and generating axial intensity profile I_i to provide information on the backscattered intensity as a function of depth.

[0017] In another embodiment, the invention is an improvement in a method of using phase-resolved functional optical coherence tomography (F-OCT) comprising the steps of measuring in an interferometer complex interference fringe signals {tilde over (Γ)}_ODT(i,j) for each pixel; generating for each pixel a variance matrix σ_i,j² based on normalized variance values; and generating an intensity map I_i,j reflecting the axial distribution of backscattered intensities based on the amplitude, A(T), obtained by using the complex fringe signals, where the generated tomographic intensity map I_i,j and the standard deviation map σ_i,j as maps of different particle sizes and different scattering parameters.

[0018] The invention also includes within its scope an apparatus for performing the foregoing improved methods. In this instance the optical apparatus used for performing the improved method is the same as or similar to that used in conventional phase-resolved functional optical coherence tomography (F-OCT), albeit combined with signal processing circuitry which performs the claimed data processing steps. Such signal processing circuitry includes analog or digital data processing circuits and in particular, all digital computers and/or digital signal processors of any type now known or later devised, which are capable of performing the claimed functions.

[0019] The present invention further includes the ability to measure a Doppler frequency shift proportional to a flow velocity component axial to the probing beam. Therefore, the flow direction can be determined from the flow velocity component axial to the probing beam and the flow velocity component transverse to the probing beam.

[0020] While the apparatus and method has or will be described for the sake of grammatical fluidity with functional explanations, it is to be expressly understood that the claims, unless expressly formulated under 35 USC 112, are not to be construed as necessarily limited in any way by the construction of “means” or “steps” limitations, but are to be accorded the full scope of the meaning and equivalents of the definition provided by the claims under the judicial doctrine of equivalents, and in the case where the claims are expressly formulated under 35 USC 112 are to be accorded full statutory equivalents under 35 USC 112. The invention can be better visualized by turning now to the following drawings wherein like elements are referenced by like numerals.

[0035] The invention and its various embodiments can now be better understood by turning to the following detailed description of the preferred embodiments which are presented as illustrated examples of the invention defined in the claims. It is expressly understood that the invention as defined by the claims may be broader than the illustrated embodiments described below.

[0036] The standard deviation of the Doppler spectrum provides quantitative information on flow velocity. Using a simple geometric optics model, the linear dependence of the Doppler bandwidth on the flow velocity is theoretically derived and it is found that the effective numerical aperture (NA) of the optical objective in the sample arm in the functional optical coherence tomography system depicted in FIG. 14 and described in U.S. Pat. No. 6,549,801 determines the slope of this dependence. When flow velocity is much higher than that of the Brownian motion of the moving particle, the linear dependence of the flow velocity on the Doppler bandwidth allows accurate measurement of flow velocity transverse to the optical axis of the probing light beam.

Principle of Transverse Velocity Measurement

[0037] Phase-resolved ODT solves the compromise between spatial resolution and velocity sensitivity in the STFFT algorithm by mapping the phase change between sequential A-line signals for velocity detection. Phase-resolved ODT increases the imaging frame rate and it also increases the sensitivity of flow velocity detection by more than two orders of magnitude. The complex analytical function, {tilde over (Γ)}_j(t), of each A-line ODT digitized fringe signal, Γ_j(t), is calculated by the Hilbert transform
{tilde over (Γ)}_j(t)=Γ_j(t)+iH{Γ_j(t)}=A(t)e^iψ(t)  (1)
where H, A, ψ, j and t denotes the Hilbert transform, the amplitude and phase of the complex function {tilde over (Γ)}_j(t), A-line number and RSOD group delay time, respectively. The average Doppler frequency shift and standard deviation of the Doppler spectrum for each pixel can be obtained simultaneously by the following equation, respectively, $\begin{array}{cc}{f}_D\left(k\right)=\frac{1}{2\pi T}{\tan }^{-1}\left(\frac{\mathrm{Im}\left[\sum _{j=1}^n{\tilde{\Gamma }}_j\left(k\right)·{\tilde{\Gamma }}_{j+1}^{*}\left(k\right)\right]}{\mathrm{Re}\left[\sum _{j=1}^n{\tilde{\Gamma }}_j\left(k\right)·{\tilde{\Gamma }}_{j+1}^{*}\left(k\right)\right]}\right)& \left(2\right)\\ \sigma \left(k\right)=\sqrt{\frac{{\int }_{-\infty }^{+\infty }{\left(f-f_D\right)}^{2}P\left(f\right)df}{{\int }_{-\infty }^{+\infty }P\left(f\right)df}}=\frac{1}{2\pi }\sqrt{\frac{2}{T^2}\left(1-\frac{\sum _{j=1}^n{\tilde{\Gamma }}_j\left(k\right)·{\tilde{\Gamma }}_{j+1}^{*}\left(k\right)}{\sum _{j=1}^n{\tilde{\Gamma }}_j\left(k\right)·{\tilde{\Gamma }}_j^{*}\left(k\right)}\right)}& \left(3\right)\end{array}$
where k is the pixel index, T is the A-line scanning period, Im denotes imaginary part and Re denotes real part of the complex value, f is the Doppler frequency shift, P(f) is Doppler power spectrum. From equation (2), the velocity limit known as 2π ambiguity velocity range is obtained as, $\begin{array}{cc}{V}_{\max }^{f_d}=\frac{\lambda }{4\cos \theta }{f}_s& \left(4\right)\end{array}$
Equation (4) indicates that the maximum unambiguous velocity V_max^fd detected by the phase-resolved method is determined by the RSOD scanning speed f_s, the center wavelength of the optical source λ and cosine of the Doppler angle θ.

[0038] There are a number of factors that contribute to the broadening of the Doppler spectrum including Brownian motion, velocity gradient and turbulence and focusing beam geometry of the optical probing beam, which is not plane wave. When velocity is high, probing beam geometry dominates the broadening of the Doppler spectrum. The contribution from the focusing beam geometry can be derived from the optical geometry as show in FIG. 1. When the Doppler angle larger than tan⁻¹(2w/l_c), where w is waist radius of the optical Gaussian beam at the focus point and l_c is the optical coherence length, the Doppler bandwidth is determined by the difference between two extreme average Doppler frequency shift caused by the two optical rays at the two boundaries of the probing optical beam. The average Doppler frequency shift corresponding to these two optical rays can be written as $\begin{array}{cc}\begin{array}{cc}{f}_a=\frac{2V\cos \left(\theta -\varphi \right)}{\lambda }& \end{array}& \left(5\right)\\ f_b=\frac{2V\cos \left(\theta +\varphi \right)}{\lambda }& \left(6\right)\end{array}$
The estimated bandwidth of the Doppler signal derived from the above two equations is: $\begin{array}{cc}{B}_d=f_a-f_b=\frac{4V}{\lambda }\sin \mathrm{\varphi sin}\theta =\frac{4V\sin \theta N{A}_{e\mathrm{ff}}}{\lambda }& \left(7\right)\end{array}$
In this equation, B_d is the Doppler bandwidth, V is the flow velocity, φ is the optical aperture angle, θ is the Doppler angle and, NA_eff is the effective numerical aperture. The Doppler bandwidth B_d is defined by the above geometrical optics method.

[0039] For a Gaussian optical beam, the Doppler bandwidth B_1/e (full width at 1/e of maximum spectrum amplitude) is the inverse of the transit time spent by particles passing through the focus zone of the Gaussian optical beam. The relationship between B_1/e and B_d is: $\begin{array}{cc}{B}_{1/e}=\frac{\pi }{8}{B}_d& \left(8\right)\end{array}$
Considering the relationship between the standard deviation and the Doppler bandwidth for a Gaussian optical beam,
B_1/e=4σ  (9)
the relationship between the standard deviation and NA_eff can be derived as, $\begin{array}{cc}\sigma =\frac{\pi V\sin \theta N{A}_{e\mathrm{ff}}}{8\lambda }& \left(10\right)\end{array}$
If we include the contributions from Brownian motion and other factors that are independent of the macroscopic flow velocity, equation (10) can be modified as, $\begin{array}{cc}\sigma =\frac{\pi V\sin \theta N{A}_{e\mathrm{ff}}}{8\lambda }+b& \left(11\right)\end{array}$
where b accounts for bandwidth spectrum broadening from all of these sources. From the above theoretical result, the standard deviation is linearly dependent on velocity above a background value, and the effective numerical aperture NA_eff of the optical objective in the sample arm determines the slope of this velocity dependence.

[0040] Knowing NA_eff, the measurement of the standard deviation can be used to determine the flow velocity. The flow velocity detected by the Doppler bandwidth method will be limited by the aliasing phenomenon being described by σ≦f_s. However, because the Doppler bandwidth and standard deviation are much smaller than the average frequency shift in the flow conditions encountered in most cases, the Doppler bandwidth can extend the detection range of flow velocity. The range of velocity detection using the Doppler bandwidth method is determined by the following: $\begin{array}{cc}{V}_{\max }^{\sigma }\sin \theta =\frac{8\lambda \left(f_s-b\right)}{\pi N{A}_{\mathrm{eff}}}& \left(12\right)\end{array}$
From equations (4) and (12), we can determine the ratio of the range of velocity detection for these two methods when they are directly applied with no phase unwrapping methods used, $\begin{array}{cc}M=\frac{V_{\max }^{\sigma }}{V_{\max }^{f_d}}=\frac{32}{\pi {\mathrm{NA}}_{\mathrm{eff}}\tan \theta }\left(1-\frac{b}{f_s}\right)& \left(13\right)\end{array}$

[0041] The experimental system used for the Doppler bandwidth measurement is conventional and has been previously described in Y. Zhao et. al., Opt Left. 25, 114 (2000); and Y. Zhao et. al., Opt. Lett. 25, 1358 (2000). The optical source for the fiber interferometer is a low coherence ASE broadband source, whose output power, center wavelength, and bandwidth is 5 mW, 1300 nm and 65 nm, respectively. The source is coupled to the source arm and split into reference and sample arms by a 3-dB 2×2 coupler. A rapid-scanning optical delay line (RSOD) in a group delay scanning mode is used for 500 Hz A-line scanning, and an E-O phase modulator at 500 kHz with a ramp signal is used as a carrier frequency producer to implement a heterodyne detection scheme.

[0042] In the sample arm shown in FIG. 1, an exchangeable optical collimator 10 and an optical objective 12 (20×, N.A., 0.35) are used to focus the optical beam to the center of a capillary tube whose inner diameter is 900 μm for M-mode imaging. A 0.1% intralipid solution composed of particles (0.356 micrometer diameter) is used as the turbid media, and a syringe pump driven by a high-resolution translation stage controls its flow through the capillary tube 14. In the detection arm, the signal from the photodetector is amplified with a bandpass preamplifier and then sent to an analog-to-digital conversion and data acquisition board sampling at 5 MHz, the number of data points for each A-line data acquisition is 4096. The A-line scanning of RSOD, phase modulation of the E-O phase modulator and data acquisition of the data acquisition board are synchronized in order to get precise alignment of the original phase of the sequential A-line scans.

Results and Analysis

[0043] The standard deviation images of the intralipid flow in the center section of the capillary tube for different NA_eff, Doppler flow angle θ are captured for a velocity from 0 to 960 μm/s. The position of the focus point of the optical probing beam is fixed to the center of the tube 14. The profiles of the flow standard deviation images for three velocities are shown in FIG. 2. We extract the center part of every flow standard deviation profile for a velocity from 0 to 960 μm/s to study the relationship between Doppler bandwidth and velocity. The extracted result from the center part of every standard deviation profile corresponds to the top of the parabolic flow profile in the tube. The experimental standard deviation data and theoretical curves from equation (11) are shown in FIG. 3 when NA_eff equals 0.09 and 0.05 at Doppler angle 77°. The solid and dotted lines are the theoretical curves and experiment results respectively. The theoretical and experimental curves are in good agreement for velocities higher than 300 μm/s. For NA_eff equals 0.09 and 0.05, the background b is 66 and 64.5 Hz, respectively. At low flow velocity, the standard deviation is constant. This is because at low velocity Brownian motion, which is independent of the macroscopic flow velocity, dominates the Doppler bandwidth. The regime for velocities higher than 300 μm/s is dominated by Doppler bandwidth broadening originated from the optical beam geometry. FIG. 3 indicates that the slope of the Doppler bandwidth dependence on velocity is proportional to NA_eff. The standard deviation experimental data and theoretical calculations for different angles and two NA_eff values 0.09 and 0.05 at a velocity of 698 μm/s are shown in FIG. 4. FIG. 4 indicates that the dependence of Doppler bandwidth is insensitive to Doppler angles around 90°. The background b when the effective NA equals 0.09 and 0.05 is 65.6 and 64.8 Hz, respectively. Therefore the Doppler bandwidth component in the standard deviation can be used in this situation to measure transverse velocity. It can be used in intravascular imaging blood flow using endoscopic optical coherence tomography where the Doppler angle is around 90°, and most of the blood vessels are parallel to the surface of the skin, so this method has its merit for this reason. From FIG. 3, it was determined that there is no aliasing phenomenon even though the velocity has reached 960 μm/s, which is much higher than the largest measurable velocity with the average Doppler frequency shift mapping method described by equation (4). This is because the standard deviation is only a small fraction of the average Doppler frequency shift as explained by equation (13). The ratio M theoretical curves are plotted in FIG. 5 for two different NA_eff. FIG. 5 indicates that the unambiguous detection range of velocities using the Doppler bandwidth mapping method can be many (10˜20) times larger than the corresponding average Doppler frequency shift mapping method when the Doppler angle is in the range from 70° to 85°. Therefore the Doppler bandwidth mapping method is a promising tool for overcoming the 2π ambiguity of velocity detection in the average Doppler frequency shift mapping method.

[0044] The Doppler bandwidth can be used to measure flow velocity transverse to the probing beam direction. Using a geometric optics model, the linear dependence of the Doppler bandwidth on flow velocity is theoretically derived and it is found that the slope of this dependence is determined by the NA_eff of the optical objective. Above a certain threshold flow velocity, this linear relationship is in good agreement with experimental data. In the case where the angle between the probing optical beam and flow direction is within ±15 degrees to the perpendicular, the average Doppler frequency shift is very sensitive to angle position while the Doppler bandwidth is insensitive to flow direction. Linear dependence of the flow velocity on the Doppler bandwidth allows accurate measurement of flow velocity without precise determination of flow direction. Finally, this technique also extends the dynamic range of the average Doppler frequency shift mapping method used in the phase-resolved ODT.

[0045] Now referring to FIG. 6, the magnitude and direction of flow velocity vector {overscore (V)} is determined. As set forth above, the velocity transverse to the probing beam direction V₁, is proportional to the Doppler bandwidth. Also, the velocity axial to the probing beam direction V_A is proportional to the Doppler frequency shift. Therefore, velocity {overscore (V)} is determined, as shown, as well as the angle θ between probing optical beam and the flow direction.

Light Scattering in the Single and the Multiple Scattering Regimes

[0046] From known light scattering theory, the expected spectrum obtained from a monodisperse dilute suspension undergoing Browninan Motion has a Lorentzian shape. In the following analysis, we will differentiate between the single and the multiple scattering regimes.

[0047] In the single scattering regime, we apply the dynamic light scattering theory (DLS). Assuming weak scattering as well as non-interacting particles, the linewidth Ω of the power Doppler spectrum (in this case defined as the half width at half maximum) is proportional to the particle self-diffusion coefficient D_B as follows: $\begin{array}{cc}\Omega =q^2{D}_B;D_B=\frac{k_{B}T}{3\mathrm{\pi \eta }a}& \left(14\right)\end{array}$
where k_B is the Boltzman constant, T is the temperature, η is the viscosity of the suspending liquid and α is the hydrodynamic diameter of the scattering particle. Further, q denotes the photon momentum transfer, $q=2k\sin \frac{\theta }{2},$
where k is the wave vector in the scattering medium and θ is the scattering angle.

[0048] In an optical dense medium where light scatters multiple times before detection, the scattering angle, the polarization of the scattered wave, and consequently the momentum transfer q change with each scattering event. In order to describe this situation, we need to apply the theory of diffusing wave spectroscopy (DWS). When considering light diffusion through a medium where light scatters multiple times, the power Doppler spectrum is still a Lorentzian, but the linewidth Ω will be dependent upon the scattering properties of the medium and the photon pathlength within it. $\begin{array}{cc}\Omega =2k^2{D}_B\frac{s}{l^{*}};D_B=\frac{k_{B}T}{3\mathrm{\pi \eta }a}& \left(15\right)\end{array}$
where s is the pathlength traveled by the light in the medium, l*=l_s/(1−g) is the photon random walk length, l_s=1/μ_s is the photon scattering length, and g is the scattering anisotropy of the medium (defined as an average of the cosine of the scattering angle).
Standard Deviation Mapping of Dynamic Particle Properties

[0049] The linewidth Ω is related to the standard deviation/variance of the power Doppler spectrum. The ODT variance imaging technique offers a means for mapping the variance of the power Doppler spectrum with high scanning rates and high spatial resolution. The variance imaging technique is used to map particle dynamics.

[0050] The probe configuration is illustrated in FIG. 7. The scanning parameters are listed in Table 1. A 1310 nm OCT/ODT system includes a 1310 nm SLD source with 10.2 mW output power and a heterodyne detection technique which gives a dynamic range of 120 dB. Applying confocal optics with low NA, and thus a large confocal parameter makes it possible to perform axial scans using a static focus. This gives the possibility for higher scan rates compared to what is the case when using dynamic focusing. Both the confocal configuration of the probe and the coherence gate effectively discriminates multiple scattered light. Adequate measurements can thus be obtained even, when the axial scanning range exceeds the length of the confocal parameter. However, we can not, assume that the lateral resolution is constant all over the scanning range, which was set to 900 μm.

[0051]

	TABLE 1
	Parameter	Symbol	Value	Unit
	Source bandwidth	Δλ	80	nm
	Center frequency of the source	λ0	1310	nm
	Coherence length	lc	10	μm
	Confocal parameter	b	120	μm
	Source power	P	10.2	mW
	Lateral resolution	Δx	10	μm
	Modulation frequency	fmod	500	kHz
	Numerical aperture	NA	0.13	—
	Sampling rate	fsamp	5	MHz
	Scan rate	fs	500	Hz

[0052] The measurements are extremely sensitive to environmental vibrations. As a consequence, all measurements were done on an optical table with laminar air flow isolation (Newport Corporation, Irvine Calif.).

[0053] In the experiments we applied eight different homogeneous suspensions containing monodisperse polystyrene microsphere particles (Polysciences, Inc., Warrington Penn.), with diameters ranging From 0.2 to 2.0 μm. The suspensions were kept at room temperature T=294 K for a certain time to assure equal temperature in all the suspensions before starting the measurements. The temperature was also measured after data acquisition. In some of the experiments, we wanted to immobilize the polystyrene particles in a gel before doing measurements. For this purpose the particle suspensions were mixed in a 50:50 volume ratio with a 4% gelatin suspension.

[0054] During data acquisition, a fixed position of the focus was used, both axially arid laterally. The scanning frequency was f_s=500 Hz. The focus spot was placed right below the end of the axial scanning range. A schematic overview of the post processing technique used to obtain the data is given in FIG. 7. In step 16 a variance matrix σ_i,j² is calculated based on the complex interference fringe signals {tilde over (Γ)}_ODT(i,j). In total, n=4000 sequential axial scans are performed, and 200 sequential axial scans are used to build up each axial scan line, thus j∈[1,n/200]. The variance matrix σ_i,j² contains the normalized variance values.

[0055] In step 18 an intensity map I_i,j is calculated, reflecting the axial distribution of backscattered intensities. The intensity values are based on the amplitude A(T), obtained by using the complex fringe signals. As in step 16, 200 axial scans were used to construct each intensity line.

[0056] In step 20 the intensity is averaged over n/200 scanning lines, to obtain an axial intensity profile I_i. This is equivalent to performing temporal averaging.

[0057] In step 22 when averaging temporally over n/200 axial scanning lines to find an axial variance profile σ_i², the variance values are weighted with respect to the intensity map I_i,j ${\sigma }_i^2=\frac{\sum _{j=1}^{n/200}{\sigma }_{i,j}^{{}_{}2}{I}_{i,j}}{\sum _{j=1}^{n/200}{I}_{i,j}}$
By weighting the variance with respect to the intensity, we assure that the points with the largest signal-to-noise ratio (SNR) ratio contribute most to the results. The signal-to-noise ratio is a crucial factor, which will have influence on the broadness of the Doppler spectrum.

[0058] In step 24 the normalized standard deviation profile σ_i expresses the standard deviation of the Doppler power spectrum as a function of depth, and may identify the transition region between the single and the multiple scattering regimes.

[0059] In step 26 the axial intensity profile I_i provides information on the backscattered intensity as a function of depth in the suspension. It is an important tool for identifying the spatial extent of the single scattering regime.

[0060] FIG. 8 lists the normalized standard deviation of the Doppler spectrum in the single scattering regime, measured in polystyrene suspensions. These values are obtained by spatial averaging of the axial standard deviation profile over the range from 0 to 30 μm below the liquid surface. This corresponds to an optical depth less than one scattering length (1/μ_s) for all suspensions, except for the suspension with d=2.092 μm. FIG. 8 also shows a graphical presentation of the results.

[0061] We are interested in the value of the normalized standard deviation obtained in the absence of Brownian motion. This value was experimentally determined by mixing the polystyrene suspensions with d=0.356 μm and d=0.548 μm with a 4% gelatin suspension. In such a mixture, the viscosity factor η in equation (14) and (15) comes close to infinity, and thus the effects originating from particle movement can be neglected. The horizontal line in FIG. 8 illustrates the intrinsic value of the standard deviation, found to have the value 0.4.

Power Spectra in the Single Scattering Regime

[0062] To demonstrate that fast scanning optical Doppler tomography systems can be used to image and quantify Brownian Motion within highly scattering media, we determined the power Doppler spectra for separate monodisperse suspension of polystyrene microspheres. The measurements were done in the single scattering regime, and the linewidth Ω of the spectra was used to estimate an experimental Brownian diffusion coefficient {circumflex over (D)}_B for each sample.

Recording and Processing of Data

[0063] The data was obtained by using an ODT technique with conventional post processing, with no lateral or axial movement of the focus during data acquisition. The focus spot was placed 20 μm below the liquid surface. The data was acquired using 4000 sequential axial scans, which corresponds to a total acquisition time of 8 seconds at a scanning frequency of 500 Hz.

[0064] For each sample, we apply a post processing technique as illustrated in FIG. 9. Sequential steps 91 through 95 are shown. This made it possible to perform a Lorentzian fit to the measured discrete power Doppler spectra for various particle sizes. Knowing the Lorentzian approximation of the power Doppler spectrum, the linewidth Ω as a function of particle size in the single scattering regime, could be obtained.

[0065] Referring to step 91, the raw data of the interference fringe was obtained using a 12 bit A/D converter with a sampling frequency of 5 MHz. The real valued interference fringe functions Γ_ODT(T_i) were received and stored in a personal computer, using software developed at Beckman Laser Institute (Irvine, Calif.).

[0066] As described, the focus was fixed during data acquisition, and the interference fringe data can be described by the two dimensional matrix Γ(i,j), where i denotes the axial coordinate defined by the coherence gate, and j denotes the scan number, i∈[0,900], j∈[1, n]; n=4000.

[0067] Referring to step 92, the matrix Γ(i,j) was exported to MATLAB, where the FFT of Γ(i,j) was calculated. We thus obtained a two-dimensional matrix containing discrete power Doppler spectra, described by β(i,j).

[0068] Referring to steps 93 and 94, power Doppler spectra in each axial scanning line j ranging in depth from 20 to 37 μm below the liquid surface were averaged (spatial averaging corresponding to the pixel range a≦i≧b). Using selected scattering parameters, we conclude that this depth interval theoretically lies within the single scattering regime. Thereafter, for the selected depth range, we averaged the power Doppler spectrum over 4000 subsequent axial-lines (temporal averaging).

[0069] Referring to step 95, the discrete power Doppler spectrum obtained after performing FFT and averaging were exported to Microcal Origin (Microsoft) where a Lorentzian curve fit was performed. In order to mathematically being able to able to perform the Lorentzian fits to the power Doppler spectra, a power offset value has to be introduced. This is due to the background noise.

[0070] The linewidth Ω need from this Lorentzian fit, was used to determine {circumflex over (D)}_B, the experimental Brownian diffusion coefficient. {circumflex over (D)}_B was compared with the theoretical values, given by equation (15) with the results summarized below.

Results

[0071] FIG. 10 shows the discrete power Doppler spectrum and the Lorentzian fit for the mono disperse suspension containing the particle size d=0.356 μm.

[0072] FIG. 11 shows the right half of the discrete power Doppler spectra obtained for the polystyrene particles with diameters d=0.548 μm and d=0.992 μm. Lorentzian fits are inserted. The intention is to visualize the change in linewidth Ω of power Doppler spectra for different polystyrene suspensions. For this purpose, showing the right half of the spectrum is sufficient.

[0073] FIG. 12 illustrates the resulting relationship between the inverse particle size and the measured linewidth Ω. Error bars show the standard deviation in linewidth Ω, measured from the Lorentzian fits for 5 independent measurements. The solid line indicates a linear regression fit to the data points.

[0074] One problem that arises when measuring the power Doppler spectrum, is the background signal originating from the EOM. The effect can be observed in FIG. 10. Notice the data points with higher power values than what should be expected for frequencies around f=0 Hz.

[0075] Further, from FIG. 12 we can conclude that the statistical uncertainty in the estimate of Ω is largest for d=2.092 μm and d=0.202 μm. For d=0.202 μm, which yields Ω=171 Hz, the linewidth Ω is relatively close to the maximum detectable linewidth of 250 Hz, defined by the Nyquist sampling theorem. For d=2.092 μm, and Ω=24 Hz, the small linewidth complicates the process of fitting the Lorentzian curve to the discrete power Doppler spectrum.

[0076] As can be observed in FIG. 10, a Lorentzian curve approximates the acquired power Doppler spectra in the single scattering regime. As further shown in FIGS. 11 and 12, the linewidth Ω appears to be proportional to the inverse of the particle diameter, which is in accordance with equation (14).

[0077] Using equation (14), an experimental value of the Brownian diffusion coefficient, {circumflex over (D)}_B, could be calculated based on the measured linewidths Ω. Neglecting the measurement for the largest particles size, due to the difficulties performing the Lorentzian fit, the accuracy of {circumflex over (D)}_B compared to the theoretical value was better than 7.5%.

Power Spectra in the Multiple Scattering

[0078] Using the experimental setup described herein, the linewidth Ω of the power Doppler spectra may be studied as a function of depth in the polystyrene suspension. This is analogous to prior experiments where broadening of the spectrum with increasing depth was observed. This was interpreted as being caused by multiple scattering of the detected photons.

[0079] According to equation (15) the linewidth Ω in the multiple scattering regime should increase linearly with path length s traveled by the light in the medium. This is experimentally verified as described herein.

Recording and Processing of the Data

[0080] The focus spot was placed right below the end of the axial scanning range. Since the axial scanning range was 900 μm, we needed a suspension with a relatively short scattering length l=1/μ_s in order to detect photons originating from the multiple scattering regime. The polystyrene suspension with particle diameter d=0.548 μm with a scattering; length of 93 μm was selected.

[0081] The acquisition technique and the postprocessing methods are equivalent to the ones described above. However, in this case we selected non-overlapping depth intervals, [a,b], distributed all over the total scanning range. Power Doppler spectra for several depths could thus be obtained.

[0082] The results are illustrated in FIG. 13 where the measured linewidth Ω as a function of the pathlength s traveled by the light in the suspension of polystyrene particles with diameter d=0.548 μm. The horizontal solid line shows the approximation found in the single scattering regime. The broken line shows the approximation found for the multiple scattering regime, obtained by using equation (15).

[0083] As we detect photons having longer pathlengths s, the linewidth of the power Doppler spectrum increases rapidly. As s increases, according to the Nyquist limit, aliasing may represent a problem when using a scanning frequency f_s of 500 Hz. Constructing a scanning mechanism capable of increasing f_s may represent a solution.

[0084] Furthermore, the signal power decreases exponentially with depth, and as we detect power spectra deeper in the suspension, the poor signal-to-noise ratio (SNR) complicates the process of performing the Lorentzian fits.

[0085] The estimated value of {circumflex over (D)}_B in the multiple scattering regime, given in table 3.4, is approximately 2.3 times smaller than expected. The fast scanning OCT system used in our experiment seems not capable of resolving particle dynamics in the multiple scattering regime. Increasing f_s and SNR may improve our experimental results.

[0086] FIG. 13 indicates that for a photon pathlength in the sample between 300 μm and 400 μm, we observe that the linewidth starts to increase almost linearly. This is an indication identification of the transition region between the two scattering regimes may be possible.

[0087] Many alterations and modifications may be made by those having ordinary skill in the art without departing from the spirit and scope of the invention. Therefore, it must be understood that the illustrated embodiment has been set forth only for the purposes of example and that it should not be taken as limiting the invention as defined by the following claims. For example, notwithstanding the fact that the elements of a claim are set forth below in a certain combination, it must be expressly understood that the invention includes other combinations of fewer, more or different elements, which are disclosed in above even when not initially claimed in such combinations.

[0088] The words used in this specification to describe the invention and its various embodiments are to be understood not only in the sense of their commonly defined meanings, but to include by special definition in this specification structure, material or acts beyond the scope of the commonly defined meanings. Thus if an element can be understood in the context of this specification as including more than one meaning, then its use in a claim must be understood as being generic to all possible meanings supported by the specification and by the word itself.

[0089] The definitions of the words or elements of the following claims are, therefore, defined in this specification to include not only the combination of elements which are literally set forth, but all equivalent structure, material or acts for performing substantially the same function in substantially the same way to obtain substantially the same result. In this sense it is therefore contemplated that an equivalent substitution of two or more elements may be made for any one of the elements in the claims below or that a single element may be substituted for two or more elements in a claim. Although elements may be described above as acting in certain combinations and even initially claimed as such, it is to be expressly understood that one or more elements from a claimed combination can in some cases be excised from the combination and that the claimed combination may be directed to a subcombination or variation of a subcombination.

[0090] Insubstantial changes from the claimed subject matter as viewed by a person with ordinary skill in the art, now known or later devised, are expressly contemplated as being equivalently within the scope of the claims. Therefore, obvious substitutions now or later known to one with ordinary skill in the art are defined to be within the scope of the defined elements.

[0091] The claims are thus to be understood to include what is specifically illustrated and described above, what is conceptionally equivalent, what can be obviously substituted and also what essentially incorporates the essential idea of the invention.