Light reflection imaging method for acquiring optical parameters and microstructures of tissues in a large area

The present application is a National Phase entry of PCT Application No. PCT/CN2015/093184, filed Oct. 29, 2015, which claims the benefit of Chinese Patent Application No. 201510666945.1, filed Oct. 15, 2015, which are incorporated herein by reference.

The present disclosure provides a light reflection imaging method for acquiring optical parameters and microstructures of tissues in a large area, which can be used in the fields of biomedical optics, remote sensing and the like.

Elastic scattering of light has long been used for analyzing random media. Reflectance spectroscopy and imaging are widely used noninvasive methods for measuring the optical properties of random media (e.g., atmosphere, oceans and tissues), including an absorption coefficient (μ a ) and a reduced scattering coefficient (μ s ′). These parameters can provide valuable information about the microstructures and the biochemical components of the media, and have been applied in the fields of cloud remote sensing, monitoring of cell apoptosis, skin characterization, cancer detection and the like. Since radiative transfer (RT) describes the propagation of light in random media, the reflectance of scattered light is essentially a difficult problem. Particularly in the case of a short light source-detector distance, the diffusion approximation usually adopted for the RT cannot work. Therefore, it is still difficult to quantify a phase function of a medium from reflectance measurement, particularly to acquire a medium phase function containing basic information of a microenvironment of the relevant medium. In the case of a random light source-detector distance, an accurate analysis model for reflectance is highly desirable. This model will be applied to rapid quantitative evaluation on optical properties, especially on a phase function of a random medium. Those reflectance empirical models in the case of a short light source-detector distance proposed recently have respective limitations. Since the phase function of the scattering medium has a significant influence on the sub-diffuse reflection in the case of a short light source-detector distance, an explicit analysis model that is related to the relation between the sub-diffuse reflectance and the phase function and can deduce the optical parameters (including the phase function) of the random medium from the reflectance distribution is extremely desirable.

The present disclosure relates to a light reflection imaging method for acquiring optical parameters and microstructures of tissues in a large area. The core of the method is to combine small-angle approximation (SAA) of radiative transfer to give a quantitative analysis relation between a near distance sub-diffuse scattering light reflectance and a scattering medium phase function.

The present disclosure provides a light reflection imaging method for acquiring optical parameters and microstructures of tissues in a large area, according to

respectively obtaining a low-frequency reflectance formula

I snake ( q )+ I diffuse ( q )+μ b /2μ t

and a high-frequency reflectance formula

I SAA ( q ),

wherein q c ˜2πβ,

obtaining an absorption coefficient μ a and a reduced scattering coefficient μ s ′ from the obtained reflectance through the low-frequency reflectance formula and an inversion method in sequence,

and obtaining an anisotropic factor g, a propagation length scale l θ and an isotropic scattering background p b of a random phase function of a scattering medium with the inversion method through the high-frequency reflectance formula.

The inversion method is a table lookup method or a formula fitting method.

The high-frequency reflectance formula

is obtained from an SAA diffusion function

S ( q,z;s ⊥0 )=exp(− iq·s ⊥0 z )exp[−∫ 0 z μ t ( z−ζ ) dζ+∫ 0 z μ s ( z−ζ )χ( qζ,z −ζ) d ζ]

Assuming that p Forward (θ) satisfies Gaussian distribution, the random phase function p(θ) of the scattering medium is modified into:

p SAA ⁡ ( θ ) = 1 - 2 ⁢ ⁢ p b π ⁢ ⁢ Θ 2 ⁢ exp ⁡ ( - θ 2 Θ 2 ) + 1 2 ⁢ ⁢ π ⁢ p b ,

from which

are obtained,

through p iso =2p b , the SAA diffusion function is transformed into

When ql t >>1 and p<<l t ,

the high-frequency reflectance can be simplified into

wherein μ t ′≡μ a +2p b μ s and the propagation length scale l Θ ≡Θ/μ s .

Parameters characterizing the microstructure of the system or a combination of these parameters are/is directly obtained using the high-frequency formulas I SAA (q) and I SAA (ρ) from high-frequency reflection maps of two or more media.

The low-frequency reflectance formulas

are obtained from the expressions of snake photons and diffuse photons

The absorption coefficient and the reduced scattering coefficient are directly obtained using the formula I snake, diffuse (q) or I snake, diffuse (ρ) from a low-frequency reflection map of one or more media.

In addition, the sub-diffuse scattering light reflectance can be first fitted to obtain the values of μ b /μ s , μ a /μ s and θ/μ s ; after the μ b /μ s is determined, it is supposed that

g =(1−2 p b )(1−Θ 2 /2)

and then the sub-diffuse light and diffuse light reflectance distribution are fitted using least squares fitting to obtain all optical parameters.

The present disclosure includes a turbid medium reflectance calculation method applicable at a random spatial distance and in an entire spatial frequency domain and a method of measuring the reflectance of a turbid medium at high and low spatial frequencies and inverting the obtained light reflectance to obtain optical parameters of the medium. The inversion method may be a table lookup method or a formula fitting method, etc. The measurement of sub-diffuse and diffuse light reflectance of the turbid medium can be used for measuring the optical properties of the turbid medium and microstructures including a phase function and the like in a large area. The phase function of the scattering medium carries basic information about the morphology and optical properties of a single scatterer. In addition, the analysis on the phase function can be used for predicting light propagation and detecting changes or inhomogeneity of microstructures in a random medium. The analysis model proposed by the present disclosure will be extensively and importantly applied in rapid quantification of all optical properties of a scattering medium, including a phase function and the like, especially in the fields of biomedical optics, remote sensing and the like.

FIG. 1 is a scattering schematic illustration of SAA photons according to the present disclosure.

FIG. 2 is a contrast schematic illustration of light reflectance of a suspension of polystyrene spheres having a diameter of 1.5 mm in water based on Monte Carlo simulation.

FIG. 3 is a schematic illustration 1 of approximating an exact Mie phase function using an SAA phase function.

FIG. 4 is a schematic illustration 2 of approximating the exact Mie phase function using the SAA phase function.

FIG. 5 is a schematic illustration 3 that approximates the exact Mie phase function using the SAA phase function.

FIG. 6 is a schematic illustration 1 of spatial distribution of reflectance of parallel beams perpendicularly incident at an origin at a radial position r away from the origin.

FIG. 7 is a schematic illustration 2 of spatial distribution of reflectance of parallel beams perpendicularly incident at the origin at the radial position r away from the origin.

FIG. 8 is a schematic illustration 3 of spatial distribution of reflectance of parallel beams perpendicularly incident at the origin at the radial position r away from the origin.

FIG. 9 is a schematic illustration 1 of reflectance of normal-incidence spatial modulated plane waves at a spatial frequency q.

FIG. 10 is a schematic illustration 2 of reflectance of normal-incidence spatial modulated plane waves at the spatial frequency q.

FIG. 11 is a schematic illustration 3 of reflectance of normal-incidence spatial modulated plane waves at the spatial frequency q.

FIG. 12 is a schematic illustration of a result of fitting the polystyrene suspension shown in FIG. 2 through the formula I SAA (q).

The embodiments of the present disclosure will be further illustrated below in conjunction with the accompanying drawings.

A collimated beam is incident on an interface of z=0 along the direction s in ≃ , and passes through a forward peak scattering medium, with a light reflectance of s in ≃ in the backscattering direction. In such medium, non-diffuse photons only undergo a few times of large-angle scattering, and are divided into n-order non-diffuse photons according to the times n of large-angle scattering. Since the first-order non-diffuse photons undergo multiple small-angle scattering and one large-angle scattering, the main contribution can be described as a near distance form using small-angle scattering approximation to radiative transfer. The first-order non-diffuse photons are referred to as SAA photons. The second-order non-diffuse photons are referred to as snake photons.

Backscattering of the SAA photons is mainly determined by the diffusion of a forward scattering angle and the backscattering efficiency. It is supposed that the phase function p(θ) (normalized into 2π∫ 0 π p(θ)sin θdθ=1) of the scattering medium is divided into a forward peak scattering component and an isotropic component, e.g.,

P SAA (θ)=(1−2 p b ) P Forward (θ)+(2π) −1 p b ( p b <<1).

FIG. 1( a ) shows backscattering of first-order non-diffuse photons undergoing multiple small-angle scattering (“red” scatterer) and single large-angle scattering (“blue” scatterer). At the first or last scattering event, the photons can have isotropic scattering (dotted line) and forward peak scattering (solid line) using the probabilities P iso and 1−P iso respectively. (b) The introduction of S eff ′ improves the precision of SAA photons, and shows the contrast of light reflectance of a suspension of polystyrene spheres having a diameter of 1.5 mm in water based on Monte Carlo simulation.

In a layered medium, the SAA diffusion function for a collimated beam incident at the origin (r=0) in the direction s 0 to the depth z is given by the following formula:

S ( q,z;s ⊥0 )=exp(− iq·s ⊥0 z )exp[−∫ 0 z μt ( z −ζ) dζ+∫ 0 z μs ( z −ζ)χ( qζ,z −ζ) d ζ]

wherein q is a spatial frequency on the xy plane, s ⊥0 is the projection of S0 on the interface, μ t =μ s +μ a , in which μ s is a scattering coefficient, and χ(v, z) is two-dimensional Fourier transformation of (1−2p b )p Forward (θ,z). The second index item is expanded into

exp(−∫ 0 z μ t (ζ) d ζ)[1+∫ 0 z μ s ( z −ζ)χ( qζ,z −ζ) d ζ+ . . . ]

with a 0-order ballistic item, so that S coming from each scattering orders 0, 1, . . . can be identified.

Then, the reflectance of the SAA photons is obtained, which can be expressed as

In the above formula, the backscattering coefficient μ b ≡μ s (z)p(π,z), s ⊥ =s ⊥in +s ⊥out , S eff is a diffusion function for a valid medium with the same phase function and double absorption and scattering, S eff ′ is a diffusion function for a second valid medium with the same phase function and scattering and a absorption coefficient modified to 2μ a +2p b μ s . The second item of formula (2) above indicates that when the photons enter or escape from the medium, they pass through an isotropic scattering path instead of a forward peak scattering path, and the probabilities of the two paths are P iso and 1−P iso (P iso <<1) respectively (see FIG. 1 ). The introduction of the second item improves the accuracy of SAA (see FIG. 2 ). It should be noted that the ballistic item should be removed from the S eff ′ to avoid recalculation.

Assuming that p Forward (θ) satisfies Gaussian distribution, the random phase function p(θ) of the scattering medium is modified into

p SAA ⁡ ( θ ) = 1 - 2 ⁢ p b π ⁢ ⁢ Θ 2 ⁢ exp ⁡ ( - θ 2 Θ 2 ) + 1 2 ⁢ ⁢ π ⁢ p b , ( 3 )

from which

p b = 2 ⁢ ⁢ π ⁢ ∫ π / 2 π ⁢ p ⁡ ( θ ) ⁢ sin ⁢ ⁢ θ ⁢ ⁢ d ⁢ ⁢ θ ⁢ ⁢ and Θ 2 = 2 ⁢ ⁢ π 1 - 2 ⁢ ⁢ p b ⁢ ∫ 0 π ⁢ p ⁡ ( θ ) ⁢ θ 2 ⁢ sin ⁢ ⁢ θ ⁢ ⁢ d ⁢ ⁢ θ

are obtained,

The n-order moment of P SAA is obtained from

( 1 - 2 ⁢ ⁢ p b ) ⁡ [ 1 - n ⁡ ( n + 1 ) 4 ⁢ Θ 2 ] .

Specifically, the anisotropic factor (n=1) is

g =(1−2 p b )[1−Θ 2 /2].  (4)

Through p iso =2p b in formula (2), the SAA diffusion function, namely formula (1), is simplified into

Wherein erf is an error function. Compared with the order of other approximate solution truncated phase function moments of radiative transfer, the reflectance (2) of the SAA photons contains the contributions coming from all moments of the phase function.

The expression of P SAA is consistent with a unified Mie and fractal model of light scattering caused by tissues and cells. The Gaussian item captures the contribution from Mie scattering, and the isotropic scattering item is correlated with the refractive index fluctuation of the background. The root mean square scattering angle θ decreases with the size of the Mie scatterer (large structure) in the tissues and cells.

When ql t >>1 and ρ<<l t ,

the high-frequency reflectance can be simplified into

I SAA ⁡ ( q ) → μ b 2 ⁢ μ t ⁡ [ 1 + 2 · 6 1 / 3 ⁢ Γ ⁡ ( 4 3 ) ⁢ μ t μ s ⁢ ( l Θ ⁢ q ) - 2 / 3 - 2 3 ⁢ 6 2 / 3 ⁢ Γ ⁡ ( 2 3 ) ⁢ μ t ⁢ μ t ′ μ s 2 ⁢ ( l Θ ⁢ q ) - 4 / 3 ] ⁢

⁢ and ( 6 ) I SAA ⁡ ( ρ ) → μ b 2 ⁢ μ t [ ⁢ δ ⁡ ( ρ ) + 1 π ⁢ ( 2 3 ) 2 / 3 ⁢ Γ ⁡ ( 2 3 ) ⁢ μ t μ s ⁢ ( l Θ ⁢ q ) - 2 / 3 - 2 π ⁢ ( 2 3 ) 2 / 3 ⁢ Γ ⁡ ( 4 3 ) ⁢ μ t ⁢ μ t ′ μ s 2 ⁢ ( l Θ ⁢ q ) - 2 / 3 ] , ( 7 )

wherein μ t ′≡μ a +2p b μ s and the propagation length scale l Θ ≡Θ/μ s .

Parameters characterizing the microstructure of the system or a combination of these parameters are/is directly obtained by inversion using the high-frequency formulas I SAA (q) and I SAA (ρ) from high-frequency reflection maps of two or more media.

The low-frequency reflectance formulas

are obtained from the expressions of snake photons and diffuse photons

In the expression, β≡μ a +μ s ′, μ s ′≡μ s (1−g), g is an anisotropic factor, G( snake, diffuse ) is a Green's function of snake photons and diffuse photons respectively, G (snake) (r,r′)=exp(−β|r−r′|)/4π|r−r′| 2 is a Green's function of snake photons, and the snake photons are propagated along an isotropic source trajectory in an isotropic scattering turbid medium. In the Fourier domain, the reflectance of the snake photons and the diffuse photons can be simplified into:

wherein Q≡√{square root over (q 2 +3μ a μ s ′)}, Z e is an extrapolation length dependent on refractive index mismatch at the interface.

The absorption coefficient and the reduced scattering coefficient can be directly inverted using the formula I snake, diffuse (q) or I snake, diffuse (ρ) from a low-frequency reflection map of one or more media.

FIG. 3 , FIG. 4 and FIG. 5 compare the results of simulating the suspension of polystyrene spheres having the diameters of 1.5 μm (not absorbed and absorbed, μ a /μ s ′=0.16) and 0.49 μm (not absorbed) in water in a spatial domain and a Fourier domain based on Monte Carlo simulation as well as the light reflectance of SAA photons, snake photons and diffuse photons.

In the Monte Carlo simulation, the total number of incident photons is set to 106, assuming that the semi-infinite medium is matched with the surrounding refractive index. Absorption of the scattering medium is achieved by adding a non-zero imaginary part to the refractive index of the polystyrene particles. The parameters corresponding to the SAA photons are p b =0.0177 and Θ=0.451, p b =0.0171 and Θ=0.447, p b =0.0194 and Θ=0.587, respectively.

FIG. 6 , FIG. 7 and FIG. 8 compare the reflectance of SAA photons, and snake and diffuse photons and the results of simulating the suspension of polystyrene spheres having the diameters of 1.5 μm ( FIG. 6 : not absorbed, FIG. 7 : μ a /μ s ′=0.16; g=0.92) and 0.49 μm ( FIG. 8 : not absorbed, g=0.86) in water based on Monte Carlo simulation. The wavelength of incident light is 0.515 μm.

FIG. 3 to FIG. 5 show using the SAA phase function to approximate the exact Mie phase function. FIG. 6 to FIG. 8 show the spatial distribution of reflectance of parallel beams perpendicularly incident at an origin at a radial position r away from the origin. FIG. 9 , FIG. 10 and FIG. 11 show the reflectance of normal-incidence spatial modulated plane waves at a spatial frequency q.

Table 2, Comparison of parameters fitted by Monte Carlo simulation on the light reflectance of polystyrene sphere suspension (d=1.5 μm, μ a /μ s ′=0.16) and theoretical values.

The SAA photons and the combined snake and diffuse scattering photons well describe the reflectance of a light source and a detector in a real space over short and long distances ( FIGS. 6-8 ), and the reflectance in a frequency domain at high and low spatial frequencies ( FIG. 9-11 ).

When the distance is greater than β −1 , this division on the backscattering light is consistent with the general radial distribution independent from the specific form of a phase function. Formula (8) perfectly describes the light reflectance of a forward peak scattering medium with low and medium absorption (μ a /μ s ′<1). Table 1 shows the similarity calculated within 0≤q/μ s ′≤100 between the Fourier space model and the Monte Carlo simulation. The similarity between the model and the Monte Carlo simulation is defined as follows

R 2 = 1 - SSR SST ,

wherein SSR=Σ(log I modal −log I MC ) 2 is a mean square error of the model, SST=Σ(log I MC − log I MC ) 2 is a true value deviation, log I MC is the mean of log I MC . R 2 indicates perfect overlap matching of the model and the true value. The similarity value is calculated within the range 0≤q/μ s ′≤100 in the Fourier space. When the absorption increases (μ a /μ s ′>≈1), the photons are no longer completely randomized in the propagation direction thereof or the light is no longer scattered toward forward peaks in the medium, and the precision of formula (8) decreases.

The closed form of the light reflectance is:

The closed form (8) of the light reflectance and the limiting forms (6) and (7) show, in addition to well-known transfer mean free path l t =1/(1−g)μ s acting on light diffusion and reflectance at low spatial frequencies (long distance), there is a new propagation length scale l Θ ≡Θ/μ s characterizing the reflectance of sub-diffuse light at high spatial frequencies (short distance), which is determined only by the angle width of the forward peak portion of the phase function. For the forward peak scattering medium, an internal relation between the above two dimensions can be further established through formula (4), and approximately satisfies μ s 3 l t l θ 2 =2. The high-frequency (short distance) light reflectance can be conveniently measured through formulas (6, 7). FIG. 12 shows the results of fitting the polystyrene suspension shown in FIG. 3 to FIG. 11 through formula (6). The diffusion angle Θ of the scattering light is determined by the precision of 2.4%, 3.6% and 5.6%, respectively.

A complete group of optical properties of the turbid medium is determined by fitting the light reflectance at low and high spatial frequencies, and the group of optical properties include an SAA phase function defined entirely by the diffusion angle Θ of the scattering light, an isotropic scattering background p b of a random phase function of a scattering medium, a backscattering coefficient μ b , a scattering coefficient μ s , an absorption coefficient μ a and an anisotropic factor g.

For example, table 2 shows the parameter results of fitting the light reflectance of polystyrene suspension (d=1.5 μm, μ a /μ s ′=0.16) based on Monte Carlo simulation.

In addition, the sub-diffuse scattering light reflectance can be first fitted to obtain the values of μ b /μ s , μ a /μ s and Θ/μ s , g=(1−2p b )(1−Θ 2 /2) is assumed after the μ b /μ s is determined, and then the sub-diffuse light and the diffuse light reflectance distribution are fitted using least squares fitting to obtain all optical parameters.

The group of optical parameters obtained from the light reflectance agree well with the theoretical values. In particular, the precision of the extracted scattering characteristics and phase function parameters (μ, g, Θ) is very good, partly because of the constraints between g and Θ. g and Θ characterize the diffuse scattering and sub-diffuse scattering light reflectance caused by the forward peak scattering medium, respectively.

In summary, the present disclosure proposes a reflectance analysis model describing a forward peak scattering medium at a random light source and detector distance, and its application in analyzing the optical parameters and phase functions of the forward peak scattering medium. When the absorption is at medium to low levels, the analysis model shows excellent performance over the entire spatial scale, and it has been successfully proved that the model can be applied to accurately determine the optical properties and phase functions of the turbid medium from the sub-diffuse and diffuse light reflectance of the medium. The phase function of the scattering medium carries basic information about the morphology and the optical properties of a single scatterer. In addition, the analysis on the phase function can also be used for predicting light propagation and detecting the changes or inhomogeneities of microstructures in a random medium. With the development of spatial modulated illumination in recent years, it has become possible to quickly measure the reflectance of sub-diffuse and diffuse scattering light in a wide field of view. Therefore, the analysis model proposed by the present disclosure will be extensively and importantly applied in rapid quantification of all optical properties of a scattering medium, including a phase function and the like, especially in the fields of biomedical optics, remote sensing and the like.