WO2005008195A2

WO2005008195A2 - Quantification of optical properties in scattering media using fractal analysis of photon distribution measurements

Info

Publication number: WO2005008195A2
Application number: PCT/IB2004/002303
Authority: WO
Inventors: David H. Burns; Claudia E. W. Gributs
Original assignee: Mcgill University
Priority date: 2003-07-16
Filing date: 2004-07-15
Publication date: 2005-01-27
Also published as: WO2005008195A3

Abstract

An optical property of a turbid medium to be tested is determined. A light measurement device (15) is used to obtain photon distribution data regarding a parameter from light exiting from the medium (one example is the time-of-flight intensity distribution). A fractal element analyzer (16) receives the distribution data and produces a count of a number of points of the distribution in a partition element for each data point in the distribution, for each of the partition element sizes. An optical property estimator (18) receives the count and compares characteristics of a relationship of the optical property of the medium to be tested to characteristics of a relationship of the optical property obtained for a plurality of calibration media (having known optical properties) to output a value for the optical property of the medium.

Description

QUANTIFICATION OF OPTICAL PROPERTIES IN SCATTERING MEDIA USING FRACTAL ANALYSIS OF PHOTON DISTRIBUTION MEASUREMENTS

FIELD OF THE INVENTION The invention relates to a method of calculating one or more optical properties of a turbid medium from time resolved light measurements, or polarization resolved light measurements, or spatially resolved light measurements.

BACKGROUND OF THE INVENTION Quantification in highly scattering media is not straightforward. Light paths of photons through scattering samples are described by statistical distributions, and the logarithmic Beer-Lambert relationship does not apply to such cases. Accurate optical characterization of samples such as biological tissues, pharmaceutical powders or bioreactors is therefore difficult. To optically characterize a turbid medium, there are several optical parameters which must be determined. In particular, there is the absorption and scattering coefficients, scattering anisotropy, media and scattering index of refractions, and optical chirality. The absorption coefficient, μ_a, represents the frequency of absorption events in a 1 millimeter path length. Similarly, the scattering coefficient, μ_s, is a measure of the number of scattering events per millimeter. The scattering anisotropy, g, is the average cosine of the angle that light is scattered at each scattering event. Both the media and the scattering centers in a sample will have an index of refraction, n_m and n_s. Also, a sample can contain optically active, chiral molecules, which can rotate the polarization and can be described by a degree of polarization, α. Both μ_a and μ_s can be obtained from time-resolved measurement of light traversing a sample, since the time spent by an individual photons in a medium before exiting it or be absorbed by it is a function of the number of scattering events as well as the absorption properties of said medium. They can also be determined by a space-resolved or angle-resolved measurement of light exiting a medium since the distance traveled by a photon in a medium is also a function of the number of scattering events and of the absorption properties of said medium. A polarization resolved measurement of light traversing a sample can as well be used to determine the scattering and the absorption coefficients of a medium. With the scattering and absorption determined, the other optical properties may also be determined using the above method as a model of light propagation. There has been considerable effort made to understand the utility of time distributions in measuring absorption and scattering in a sample. Understanding these profiles is important to obtain quantitative measurements in scattering media. One approach is to use the diffusion approximation to the radiative transport equation which is based on the two optical coefficients to model the time-resolved measurements. One of the limiting factors of this model is the description of the short photon paths through the sample, i.e. those photons that experience little scattering. Such photons have not undergone large multiple angle scattering and their behavior is inadequately described by the diffusion equation. This approach can also be mathematically quite involved. There is a need for a technique that could accurately estimate the absorption and scattering coefficients of a turbid medium. In this invention, a new and robust mathematical method is proposed that makes use of the fact that the photon distributions intrinsically contain absorption and scattering information.

SUMMARY OF THE INVENTION According to the invention, there is provided a data analysis method involving one element of fractal analysis of a photon distribution. The invention also provides an optical property measurement apparatus, operating as an imaging device or spectrometer, that measures the photon distribution and analyses it using an element of fractal analysis to obtain optical parameter data of a turbid medium. According to one broad aspect of the invention, there is provided a method of determining an optical property of a turbid medium to be tested, comprising obtaining a photon distribution in respect to a parameter from the measurements of exiting light from the medium, selecting a series of partition element sizes wherein each partition element corresponds to an area of the distribution, counting a number of points of the distribution in the partition element for each data point in the distribution, for each of the partition element sizes, and establishing a relationship between the number for each partition element size and the partition element size. The preceding steps are to be repeated for a plurality of calibration media having a known optical property. This then allows for the optical property of the medium to be tested to be determined by comparing characteristics of the relationship of the medium to be tested to characteristics of the relationship obtained for the plurality of calibration media. According to another broad aspect of the invention, there is provided an apparatus for determining an optical property of a turbid medium to be tested. A light measurement device is used to obtain photon distribution data regarding a parameter from light exiting from the medium. A fractal element analyzer receives the distribution data and produces a count of a number of points of the distribution in a partition element for each data point in the distribution, for each of the partition element sizes. An optical property estimator receives the count and compares characteristics of a relationship of the optical property of the medium to be tested to characteristics of a relationship of the optical property obtained for a plurality of calibration media to output a value for the optical property of the medium. One example of the photon distribution data is a time-of-flight intensity distribution.

BRIEF DESCRIPTION OF THE DRAWINGS For a more complete understanding of the present invention and the advantages thereof, reference is now made to the following descriptions taken in connection with the accompanying drawings in which: FIG. 1 is a schematic drawing illustrating the instrumentation set-up of a preferred embodiment of the present invention. BS: beam splitter, PD: photodiode, CFD: constant fraction discriminator, DL: delay line, ND: neutral density filter, S: sample, MCP-PMT: microchannel plate photomultiplier tube, TAC: time-to- amplitude converter, PC: personal computer; FIG 2A is a graph illustrating several photon time of flight distributions (TOF) that were obtained for samples having different absorption coefficients but same scattering coefficient; FIG. 2B is a graph illustrating several photon time of flight distributions (TOF) that were obtained for samples having different scattering coefficients but same absorption coefficient; FIG. 3 is a flowchart illustrating a method to determine at least one of /.a and //_s. Fig. 4 is a flowchart illustrating a method to calculate some fractal characteristics of photon TOF profiles ; Fig. 5 is a schematic drawing illustrating a method to compute a correlation sum; Fig. 6 is a graph illustrating the logarithm of a correlation sum with respect to the logarithm of a partition element size; Fig. 7 is a schematic drawing illustrating a matrix of calibration samples; Fig. 8 is a flowchart illustrating a method to determine relations between some regressions parameters and μ_s and μ_a; Fig. 9 is a schematic drawing illustrating the instrumentation set-up of a second embodiment fo the present invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT In an embodiment of the present invention, a method is provided for estimating optical parameters such as scattering and absorption coefficients (μ_a and μ_s) in turbid media from photon time of flight distributions using fractal dimension analysis of the signal. This new method greatly simplifies the analysis of the data The present invention and advantages thereof are best understood by referring to FIGS. 1 through 8, like numerals being used for corresponding parts of the various drawings. Measurements of the photon time of flight distributions are made using a single photon counting instrument. With reference to FIG.1, an ultrafast laser 1 produces repeated light pulses 2 that are split via a beam splitter 3 into two portions: approximately 4% of the light pulse intensity is sent to a PIN photodiode 5 and the remaining light is injected in a turbid medium 7 after going through a neutral density filter 4. Light injected in the turbid medium -referred therein as sample - is scattered and partly absorbed, and a cooled micro-channel plate photomultiplier tube (MCP-PMT) 6 detects photons individually coming out from the sample. In this particular embodiment the said PMT detects photons that are "reflected" by the sample and this detection geometry is therefore called a reflectance measurement. For each laser impulse, the time spent in the sample by each detected photon is determined by measuring the time delay between the PIN photodiode and the PMT signals. The signals detected are sent respectively to two constant fraction discriminators 8 that send two logic signals (a start signal and a stop signal) to a time-to amplitude converter 9. The time difference between the start and stop pulses, which is related to the time spent by a photon in the sample, is converted to a corresponding dc voltage sent to the measurement controller 11. A photon time of flight distribution can be constructed this way. In this preferred embodiment, light pulses of approximately 50 kW and having a time width of approximately 170 femtoseconds (at full width at half- maximum) were produced at a repetition rate of 75 MHz with a mode-locked TkSapphire laser. Typically, a few minutes of integration are necessary to obtain time-of-flight profiles with a total number of counts between 10⁶ to 10⁷ FIG.2A is a graph illustrating several photon time of flight distributions

(TOF) 12 that were obtained using the above described method for samples having different absorption coefficients but the same scattering coefficient: the graph presents the total number of photons measured 13 at a specific time delay 14. In the case of FIG. 2B, the samples had the same absorption coefficient but different scattering coefficients. Clearly these measured photon time of flight distributions contain information on both the absorption and scattering properties of the sample, as it can be seen from these graphs that a change of μ_a orμ_s modified the shape of the TOF 12. The time of flight distribution 12 is then analyzed with in view to determine its fractal characteristics 16. FIG. 3 is a flowchart illustrating how this method is implemented. Once the fractal characteristics 16 of the medium have been determined they can then be compared with pre- established standard curves 17 relating fractal characteristics to known optical values, and at least one of μ_a orμ_s can thus be determined 18. Fig. 4 gives more details on the method 27 that is used in the preferred embodiment to calculate some TOF fractal characteristics. The first step 19 is to calculate a greatest distance, Γ_MAX_> between any two data points of the TOF. The second step 20 is to select a' series of different values of partition element size smaller than r_MAX- A partition element size is then selected from these values in third step 21. Then in step 4, a correlation sum of the TOF for this selected partition element size is computed 22. In a preferred embodiment of this invention, the correlation sum C(r) is defined as: c(_r) (1)

Where Ν is the number of points in the TOF distribution, H is the Heaviside function, η and η are the coordinates of two points in the data set and r is the partition element size. Practically, the number of points that lie within a distance r of each data point are counted, summed over the data set and normalized. In step 5 the logarithm of C(r), log (C(r)), is calculated as well as the logarithm of the partition element size, log(r), 23. The process from step 3 to step 5 is then repeated until the entire series of selected partition element sizes have been scanned 24. Fig. 5 illustrates schematically the above process. In step 6, a linear relation is calculated between the logarithm of above computed correlation sums and the logarithm of partition element sizes.

25. This linear relation is of the form log (C(r)) = m*log(r) + b and can be determined by performing a linear least square regression. In step 7, the slope, m, and the parameter, b, of the linear relation are calculated as well as a parameter X_b which is equal to -b/m. The slope, m, is equal to the fractal dimension of the TOF data set. Referring to Fig. 4, after having completed steps 1 through 5, the logarithm of the calculated correlation values can be plotted in respect to the logarithm of the partition element size as shown in Fig. 6. This plot may help to discard points of the curve 30 that are associated to experimental noise, since noisy data points have a fractal dimension different from the fractal dimension of the remaining points. These noisy data points 35 are associated to the smaller partition element sizes as illustrated in Fig. 6. The remaining points 34 are thus used to calculate the linear relation 24 as described in step 6 of Fig. 4. Values such as m 31, b 33, and X 32 are this way obtained. The pre-established relations 17 between the μ_s and μ_a coefficients are obtained by performing a calibration 48 of the instrument described in Fig. 1. The said calibration is carried out by using a series of calibration samples having known μ_s and μ_a and having the same thickness as the sample 7 being analysed . Fig. 7 illustrates a possible series of calibration samples that can be chosen for this purpose, a calibration sample 36 being represented as a box. In this case 2x56 samples were used, each sample 36 representing a unique combination of μ_s and μ_a coefficients among 7 different levels of μ_s 37 and 8 different levels of μ_a 38. For each of these calibrations samples a TOF is acquired by applying the same technique previously described. Then the data analysis previously described is carried out and the values of m 31, b33, and x 32 for each calibration sample are this way obtained. These said parameters will be therein referred to as regression parameters 47. Several data analysis steps 48 are then performed in order to find which combination of regression parameters gives the most parsimonious estimation of the scattering coefficient μ_s and the absorption coefficient μ_a 46. Fig. 8 illustrates how this method 48 can be implemented. As a first step a combination of the regression parameters 47 is chosen 41. Then as a second step 42 a series of samples having the same known μ_s coefficient are selected as a first validation set 39, the remaining samples constituting a . calibration set. Using the calibration set, the best linear relation between the chosen regression parameters 47 and the known μ_s values is calculated 43 (step 3). This linear relation is of the form: μ_s = ff.m + ff₂b + <7₃x + ₄ (2) where σi, o₂, ₃, and α₄ are calibration coefficients.

In step 4 estimated values of the scattering coefficient, /_SEST_> are obtained 44 by inserting the m.b.and xb values of the samples of the validation set into equation (2). Step 2 to step 4 are repeated for all possible validation sets so that at the end of the process we can compare the estimated scattering values / J_SESTJ obtained with (2) to their known values. This is done in step 5 45 where the coefficient of variation, C.V. 47, of the model represented by equation (2) is computed. C.V is given by equation (3): C.V. = ^SEE (3) mean(μ_s) where SEE is the standard error of estimation given by: ^is -isEst 'j ' (4) N^#param with μis : known scattering coefficient of the i validation set Pis _Estj : estimated scattering coefficient obtained with the j sample of the i validation set N: number of μ_s estimates #param: number of parameters (m, b or xb) used in the model (2)

In order to find what would be the best relation (2) to use, steps 1 to 5 are repeated 49. The best relation is determined in step 6 46. The best relation is the one that uses the smallest number of regression parameters 47 and that has a corresponding C.V. that is not significantly larger than the smallest one found overall. Once determined, this best relation can be used as a pre-established relation 17 to determine μ_s of an unknown sample. The pre-established relation to determine μ_a of an unknown sample can be found by following the calibration method 48 just described with the following changes: at step 3 the regression parameters of each sample k of the validation set L are divided by the value of its best estimated scattering coefficient previously calculated. The validation sets are naturally also different and Fig. 7 illustrates a possible validation set 40 use to obtain a relation to determine μ_a. The following examples are given to illustrate the accuracy level this preferred embodiment can provide in the determination of the of μ_a and μ_s coefficients. Two independent sets of samples with a combination of 7 scattering and 8 absorption levels were analyzed at separate times and in random order. The scattering coefficients values of these samples were ranging from 100 mm^"1 to 250 mm^"1 and the absorption coefficients values were ranging form 0 to 0.035 mm^"1. After following the previously described data analysis, it was possible to determine the scattering coefficients with a C.V. value of about 3% by using only 2 regression parameters. The C.V. value increases slightly to about 4% when only one regression parameter is used. Equivalently, the absorption coefficients were determine with a C.V. value less than 10%. It is understood that the method that has been described in details can equally be used in the case of a transmittance measurement. It is also understood that the data analysis just described can be applied as well to photon distributions obtained from polarization-resolved measurement where linear or circular polarization light is injected in a medium and scattered light at a series of states of polarization is detected and measured. Since the properties of the medium have the ability to change the state of polarization of the incoming light, the time dimension in Fig.2A and 2B is replaced by a state of polarization dimension. Fig.9 illustrates schematically an experimental system that could be used to obtain a polarized-resolved photon distribution measurement. In this second embodiment, a He-Ne laser is launched onto the turbid medium after having been intensity modulated by a chopper and after polarization modulated by a photo-elastic modulator (PEM). The reflected light is collected via a beam splitter and a state of polarization θ (the electric field vector makes an angle θ with the horizontal) of the reflected light is selected via the rotating polarizer . This light intensity, \(θ), is detected by a photomultiplier and an electrical signal proportional to the light intensity is sent to a lock-in amplifier. The lock-in amplifier also receives a signal from the chopper. A synchronous detection of the light scattered by the medium is thus possible with increase sensitivity. It can be shown (as. in Studinski et.al, Journal of Biomedical Optics, July 2000, vol.5, No.3) that the scattered light intensity, l(0), is a function of the sample optical properties. Then by using the same data analysis steps that were described earlier, the optical parameters of the medium can be inferred after calibration of the instrumentation via a calibration set of media having known optical properties. It is understood that the method that has been described in detail can equally be used on photon distributions obtained from an optical coherent reflectometry or optical coherence tomography system. In these systems a phase resolved measurement of the backscattered light is measured using an interferometric technique comparing the phase of the backscattered light an unscattered beam. Since the properties of the medium change the relative distance of the incoming light, the time dimension in Fig.2A and 2B is replaced by a phase or distance state dimension. It is understood that the method that has been described in detail can equally be used on photon distributions obtained from an angle resolved optical system measuring the backscattered light at different angles in respect to the light illumination orientation. Since the properties of the medium change the relative angle distribution of the incoming light, the time dimension in Fig.2A and 2B is replaced by an angle dimension. Likewise, the same analysis steps caη be used for the analysis of a space-resolved photon distribution that is generated by measuring exiting light from different locations at various distances from the point of light injection in the medium. It is understood that the method that has been described in detail can equally be used where a multi-linear model is developed for estimation of the scattering and absorption coefficients from the correlation dimension vector. It will be understood that numerous modifications thereto will appear to those skilled in the art. Accordingly, the above description and accompanying drawings should be taken as illustrative of the invention and not in a limiting sense. It will further be understood that it is intended to cover any variations, uses, or adaptations of the invention following, in general, the principles of the invention and including such departures from the present disclosure as come within known or customary practice within the art to which the invention pertains and as may be applied to the essential features herein before set forth, and as follows in the scope of the appended claims.

Claims

WHAT IS CLAIMED IS:

1. A method of determining an optical property of a turbid medium to be tested, the method comprising: a) obtaining a photon distribution in respect to a parameter from the measurements of exiting light from said medium; b) selecting a series of partition element sizes wherein each partition element corresponds to an area of said distribution; c) counting a number of points of said distribution in the partition element for each data point in said distribution, for each of the partition element sizes; d) establishing a relationship between said number for each partition element size and said partition element size ; e) repeating steps a) to d) for a plurality of calibration media having known said optical property; f) determining said optical property of said medium to be tested by comparing characteristics of said relationship of said medium to be tested to characteristics of said relationship obtained for said plurality of calibration media.

2. The method as claimed in claim 1 , wherein said optical property comprises at least one of an absorption coefficient, μ_a, or of a scattering coefficient, μ_s, and wherein said photon distribution comprises a photon time-of-flight distribution (TOF)

3. The method as claimed in claim 2, wherein step (a) comprises: i) injecting a pulse of light into the turbid medium; and ii) detecting and counting photons exiting said medium as a function of time to- generate a TOF.

4. The method as claimed in claim 2, wherein generating said TOF comphses a fixed geometry between said light injection and said light detection.

5. The method as claimed in claim 1 , wherein step (a) comprises: i) injecting polarized light into the turbid medium; and ii) detecting light exiting said medium as a function of its state of polarization to generate a photon polarization distribution.

6. The method as claimed in claim 5, wherein generating said polarization distribution comprises a fixed geometry between said light injection and said light detection.

7. The method as claimed in claim 1 , wherein step (a) comprises: i) injecting light into the turbid medium; and ii) detecting light exiting said medium as a function of different medium positions relatively to the light injection location to generate a photon spatial distribution.

8. The method as claimed in claim 1 , wherein step (a) comprises: i) injecting light into the turbid medium; and ii) detecting light exiting said medium as a function of different exit angles relatively to the light injection plan to generate a photon angular distribution.

9. The method as claimed in claim 1 , wherein obtaining said photon distribution comprises a reflectance measurement.

10. The method as claimed in claim 1, wherein said partition element comprises a circle centered on a data point of the photon distribution.

11. The method as claimed in claim 1, wherein step (b) comprises calculating a greatest distance, Γ_MAX, between any two data points of said photon distribution and selecting a series of different values of partition element size smaller than ΓMAX.

12. The method as claimed in claim 1 , wherein step (c) comprises computing a correlation sum.

13. The method as claimed in claim 1 , wherein said turbid medium and calibration media comprise homogeneous media.

14. The method as claimed in claim 2, wherein said turbid medium is a homogeneous media and wherein step c) comprises computing a correlation sum and said relationship comprises a linear relationship between a logarithm of said correlation sum computed at said selected partition element sizes and a logarithm of said sizes.

15. The method as claimed in claim 14, wherein in step f) said characteristics comprise the parameters that defined the said linear relationship among which a fractal dimension of said photon distribution.

16. The method as claimed in claim 15, wherein determining at least one of μ_a and μ_s coefficients of said medium to be tested comprises: a) establishing an equation between a combination of said parameters, obtained from said plurality of calibration media and their known μ_a and μ_s coefficients; b) determining at least one of μ_a and μ_s coefficients of the medium to be tested by: i) obtaining the parameters, further referred to as test parameters, for said medium to be tested; ii) calculating said μ_a and μ_s of the medium to be tested using said equation established in step a) and said test parameters.

17. The method as claimed in claim 1 , wherein said turbid medium and calibration media comprise a stack of different homogeneous medium layers.

18. The method as claimed in claim 17, wherein in step c) comprises computing a correlation sum and wherein said relationship comprises a relationship between a logarithm of said correlation sum computed at said selected partition element sizes and a logarithm of said sizes.

19. The method as claimed in claim 18, wherein in step f) said characteristics comprise the parameters that defined the said relationship.

20. An apparatus for determining an optical property of a turbid medium to be tested, the apparatus comprising: a light measurement device for obtaining photon distribution data regarding a parameter from light exiting from said medium; a fractal element analyzer receiving said distribution data and producing a count of a number of points of said distribution in a partition element for each data point in said distribution, for each of the partition element sizes; an optical property estimator receiving said count and comparing characteristics of a relationship of the optical property of said medium to be tested to characteristics of a relationship of the optical property obtained for a plurality of calibration media to output a value for the optical property of said medium.

21. The apparatus as claimed in claim 20, wherein said apparatus performs the method of determining an optical property of a turbid medium to be tested as defined in any one of claims 1 to 19.