Determination of the scattering anisotropy with optical coherence tomography

In this work we demonstrate measurements with optical coherence tomography (OCT) of the scattering phase function in the backward direction and the scattering anisotropy parameter g. Measurements of the OCT attenuation coefficient and the backscattering amplitude are performed on calibrated polystyrene microspheres with a time-domain OCT system. From these measurements the phase function in the backward direction is determined. The measurements are described by the single scattering model and match Mie calculations very well. Measurements on Intralipid demonstrate the ability to determine the g of polydisperse samples and, for Intralipid, g = 0.35 ± 0.03 is measured, which is well in agreement with g from literature. These measurements are validated using the Intralipid particle size distribution determined from TEM measurements. Measurements of g and the scattering phase function in the backward direction can be used to monitor changes in backscattering, which can indicate morphological changes of the sample or act as contrast enhancement mechanism. ©2011 Optical Society of America OCIS codes: (170.4500) Optical coherence tomography; (170.3880) Medical and biological imaging; (290.7050) Turbid media. References and Links 1. D. J. Faber, F. J. van der Meer, M. C. G. Aalders, and T. G. van Leeuwen, “Quantitative measurement of attenuation coefficients of weakly scattering media using optical coherence tomography,” Opt. Express 12(19), 4353–4365 (2004). 2. F. J. van der Meer, D. J. Faber, D. M. Baraznji Sassoon, M. C. Aalders, G. Pasterkamp, and T. G. van Leeuwen, “Localized measurement of optical attenuation coefficients of atherosclerotic plaque constituents by quantitative optical coherence tomography,” IEEE Trans. Med. Imaging 24(10), 1369–1376 (2005). 3. F. J. van der Meer, D. J. Faber, J. Perree, G. Pasterkamp, D. B. Sassoon, and T. G. van Leeuwen, "Quantitative optical coherence tomography of arterial wall components," Lasers in Medical Science 20, 45-51 (2005). 4. V. M. Kodach, J. Kalkman, D. J. Faber, and T. G. van Leeuwen, “Quantitative comparison of the OCT imaging depth at 1300 nm and 1600 nm,” Biomed. Opt. Express 1(1), 176–185 (2010). 5. D. Levitz, L. Thrane, M. H. Frosz, P. E. Andersen, C. B. Andersen, S. Andersson-Engels, J. Valanciunaite, J. Swartling, and P. R. Hansen, “Determination of optical scattering properties of highly-scattering media in optical coherence tomography images,” Opt. Express 12(2), 249–259 (2004). 6. C. Xu, J. M. Schmitt, S. G. Carlier, and R. Virmani, “Characterization of atherosclerosis plaques by measuring both backscattering and attenuation coefficients in optical coherence tomography,” J. Biomed. Opt. 13(3), 034003 (2008). 7. Y. Wang, and R. Wang, “Autocorrelation optical coherence tomography for mapping transverse particle-flow velocity,” Opt. Lett. 35(21), 3538–3540 (2010). 8. J. Kalkman, R. Sprik, and T. G. van Leeuwen, “Path-length-resolved diffusive particle dynamics in spectraldomain optical coherence tomography,” Phys. Rev. Lett. 105(19), 198302 (2010). #141315 $15.00 USD Received 19 Jan 2011; revised 1 Mar 2011; accepted 3 Mar 2011; published 17 Mar 2011 (C) 2011 OSA 28 March 2011 / Vol. 19, No. 7 / OPTICS EXPRESS 6131 9. J. Pyhtila, R. Graf, and A. Wax, “Determining nuclear morphology using an improved angle-resolved low coherence interferometry system,” Opt. Express 11(25), 3473–3484 (2003). 10. A. E. Desjardins, B. J. Vakoc, G. J. Tearney, and B. E. Bouma, “Backscattering spectroscopic contrast with angle-resolved optical coherence tomography,” Opt. Lett. 32(21), 3158–3160 (2007). 11. J. Kalkman, A. V. Bykov, D. J. Faber, and T. G. van Leeuwen, “Multiple and dependent scattering effects in Doppler optical coherence tomography,” Opt. Express 18(4), 3883–3892 (2010). 12. L. Thrane, H. T. Yura, and P. E. Andersen, “Analysis of optical coherence tomography systems based on the extended Huygens-Fresnel principle,” J. Opt. Soc. Am. A 17(3), 484–490 (2000). 13. J. M. Schmitt, A. Knuttel, and R. F. Bonner, “Measurement of optical properties of biological tissues by lowcoherence reflectometry,” Appl. Opt. 32(30), 6032–6042 (11993). 14. D. J. Faber, and T. G. van Leeuwen, “Are quantitative attenuation measurements of blood by optical coherence tomography feasible?” Opt. Lett. 34(9), 1435–1437 (2009). 15. J. W. Goodman, Statistical optics (Wiley, New York, 1985). 16. T. G. van Leeuwen, D. J. Faber, and M. C. Aalders, “Measurement of the axial point spread function in scattering media using single-mode fiber-based optical coherence tomography,” IEEE J. Sel. Top. Quantum Electron. 9(2), 227–233 (2003). 17. C. Akcay, P. Parrein, and J. P. Rolland, “Estimation of Longitudinal Resolution in Optical Coherence Imaging,” Appl. Opt. 41(25), 5256–5262 (2002). 18. A. Unterhuber, B. Povazay, B. Hermann, H. Sattmann, W. Drexler, V. Yakovlev, G. Tempea, C. Schubert, E. M. Anger, P. K. Ahnelt, M. Stur, J. E. Morgan, A. Cowey, G. Jung, T. Le, and A. Stingl, “Compact, low-cost Ti:AlO laser for in vivo ultrahigh-resolution optical coherence tomography,” Opt. Lett. 28(11), 905–907 (2003). 19. H. J. van Staveren, C. J. M. Moes, J. van Marie, S. A. Prahl, and M. J. C. van Gemert, “Light scattering in Intralipid-10% in the wavelength range of 400-1100 nm,” Appl. Opt. 30(31), 4507–4514 (1991). 20. C. Chen, J. Q. Lu, H. F. Ding, K. M. Jacobs, Y. Du, and X. H. Hu, “A primary method for determination of optical parameters of turbid samples and application to intralipid between 550 and 1630 nm,” Opt. Express 14(16), 7420–7435 (2006). 21. J. R. Mourant, J. P. Freyer, A. H. Hielscher, A. A. Eick, D. Shen, and T. M. Johnson, “Mechanisms of light scattering from biological cells relevant to noninvasive optical-tissue diagnostics,” Appl. Opt. 37(16), 3586– 3593 (1998). 22. N. Bosschaart, D. J. Faber, T. G. van Leeuwen, M. C. G. Aalders, “Measurements of wavelength dependent scattering and backscattering coefficients by low-coherence spectroscopy,” J. Biomed.Opt. 16, 030503 (2011). 23. R. Michels, F. Foschum, and A. Kienle, “Optical properties of fat emulsions,” Opt. Express 16(8), 5907–5925 (2008). 24. B. Karamata, M. Laubscher, M. Leutenegger, S. Bourquin, T. Lasser, and P. Lambelet, “Multiple scattering in optical coherence tomography. I. Investigation and modeling,” J. Opt. Soc. Am. A 22(7), 1369–1379 (2005).


Introduction
Optical coherence tomography (OCT) is a non-invasive imaging modality that is used to image the morphology of turbid media. An important application of OCT is to quantitatively determine the optical properties of tissue. The attenuation coefficient of the OCT signal can be used to characterize tissue [1][2][3] and to determine the OCT imaging depth [4]. In contrast to frequently reported measurements of the OCT attenuation coefficient, quantitative measurements of the OCT backscattering are reported less often. The OCT backscattering coefficient can be quantified [5] and used for tissue characterization, visualization, and/or contrast enhancement [6]. In addition, the variation of the OCT magnitude in time can provide quantitative information about flow [7] and diffusion [8]. Quantitative measurement of the OCT (or low coherence interferometry) magnitude as a function of angle and/or wavelength [9,10] can provide information on the scattering phase function and the size of the scatterer.
The scattering phase function, which describes the angular probability distribution of the scattered light, is parameterized by the scattering anisotropy parameter g, which is the average of the cosine of the scattering angle. The scattering anisotropy plays an important role in diffusive light transport, for example in diffuse reflectance spectroscopy. Also in OCT, g plays an important role as it determines the amount of backscattering and in the quantification of multiple scattering effects in (Doppler) OCT [11]. The possibility to extract g from OCT measurements is very attractive, because, in comparison with commonly used goniometric or integrating sphere methods, the non-invasive nature of OCT allows for in-vivo application.
Two models are available for a quantitative description of the OCT signal: the single scattering model and a comprehensive model based on the Extended Huygens-Fresnel principle [12]. The first model gives good results for weakly scattering, non-absorbing media [13] and describes the OCT signal by a single exponential decay function. The latter model can be used to estimate g from a fit of the model to the OCT signal in depth, but is valid only for samples with small-angle forward scattering (e.g. high g). The obtained root-mean-square scattering angle θ rms can be related to g if the phase function is known. In this model the scattering coefficient μ s and θ rms (or g) only can be separated when multiple scattering effects are significant. However, an increase of μ s can counteract a decrease of θ rms (i.e. increase of g) in the fit model and the fit parameters are not (statistically) independent [14].
Here we present a simple method based on the single scattering model to determine g from OCT measurements. We demonstrate experimentally, for the first time to our knowledge, that this model can be used for measurements of the scattering phase function in the backward direction and in specific cases for measuring g.

Single scattering model of the OCT signal
In a loss less time-domain OCT system without focus tracking the OCT detector current i det (z) as a function of depth z is equal to the backscatter profile of the sample as function of z convoluted with the complex coherence function γ(2z/c) [15]. For a perfect mirror in air positioned in the sample arm located at z=0 (a single reflector), the OCT detector current signal is where η is the detector conversion factor from the incident light power to the electric current, Re{} is the real part of the complex coherence function, c is the speed of light, r mirror is the field reflection coefficient of the mirror, h(z) is the confocal point spread function [16], δ(z)=1 for z=0 and δ(z)=0 for all other z.  (2) where it is assumed that h(0)=1, i.e. the mirror is in the focus. For a scattering medium the situation is more complicated. A one-dimensional single scattering model is assumed where homogenously distributed scatterers all add coherently to the OCT signal. Assuming that the OCT signal for a homogenous scattering medium is the sum of all scattering contributions (see section 4.3 for a discussion), the detector current is with n med the group refractive index of the medium, µ b,NA the effective backscattering coefficient (quantifying the part of the light that is backscattered into the detection NA of the OCT system); µ s the scattering coefficient. The factor 2 in the exponent of (3) accounts for the round-trip attenuation to and from depth z. In a scattering medium with attenuation µ s the amplitude of the OCT signal can be found by extrapolating the attenuated OCT signal to z=0. The square of the OCT signal at the interface is where coherence length l c is defined in single pass and according to Schmitt et al. [13] and Goodman [15]. Note that this coherence length definition is a factor 8ln 2 0.75   smaller than the commonly used definition related to the axial resolution in OCT, which is defined as the full width at half maximum of the Gaussian-shaped coherence point spread function of the OCT amplitude [17]. The constant Q describes the heterodyne intensity back-coupling efficiency from a scattering medium compared to that of a mirror and ranges from 0 to 1 (see section 4.3). From this analysis it can be observed that µ b can be determined by dividing (2) by (4). The backscattering coefficient is with the backscattering coefficient µ b proportional to the scattering coefficient times the phase function integrated over the NA in the backscattering direction, p NA . Note that by taking the ratio of two OCT measurements additional loss factors in the OCT system do not influence the determination of µ b .
In the absence of absorption, the scattering coefficient µ s can be determined from the slope of the OCT signal. For media with absorption and described by the single scattering approximation, the light travels in a ballistic way and Beer's law can be applied to calculate the total OCT attenuation coefficient µ t , which equals µ t =µ s +µ a . Consequently, µ s can be obtained by subtracting the absorption coefficient from the total attenuation coefficient obtained from the slope of the OCT signal. The scattering phase function in the backscattering direction can be obtained by integrating the scattering phase function p(θ) over angles from π -NA to π. The phase function integrated over the NA, p NA describes the fraction of scattered photons which are detected by OCT system, i.e. p NA =µ b,NA /µ s . Consequently, p NA can be determined using (5) and µ s Since p NA is related to the phase function, the scattering anisotropy g can be determined from a determination of p NA if the shape of the phase function is known a priori.

OCT measurements
For our study we utilize a home-built time-domain OCT system, which is described in detail in Ref [4]. In summary, light from a Fianium light source is filtered to obtain a spectrum centered at 1300 nm, which is coupled into a fixed focus OCT setup. The axial resolution of this system is 9.7 ± 0.1 µm as was determined from the full width at half maximum of the OCT magnitude point spread function. The signal to noise ratio is 90 dB. The Rayleigh length, measured with a mirror in the sample arm, is 292 ± 9 µm. The corresponding Gaussian beam waist is 9.6 ± 0.2 µm and the numerical aperture NA of the sample arm lens is 0.043 ± 0.001. Prior to the experiment, the OCT system is calibrated for quantitative measurements of the backscattered power. Due to the limited dynamic range of the OCT system, the power from the mirror is measured using different calibrated neutral density filters in the sample arm. From the dependence of the OCT signal on the optical attenuation the reflected power and the OCT magnitude can be directly compared to the signals for the scattering sample (no optical attenuation in the sample arm) in Eq. (6).
OCT measurements on suspensions of scatterers are performed in a 1 mm thick glass cuvette. The cuvette is placed in the sample arm at ~70° angle relative to the incident beam to avoid specular reflections in the OCT signal. The sample arm beam is focused at the first glass-medium interface. Measurements for every solution are performed independently for 5 times. For each measurement, the average of 100 A-scans is taken. After background subtraction the OCT signal magnitude is corrected for the confocal point spread function [16]. The OCT attenuation coefficient is determined with a two parameter single exponential fit of the measured OCT signal in depth (corrected for the refractive index of water). To reduce the effects of multiple scattering, the scatterers concentration is kept low to create samples with scattering coefficients below 5 mm 1 (as calculated with Mie theory). In addition, only the first 500 data points (190 µm) of the OCT signal (starting at ~60 µm depth after sample front surface) are used for fitting the single exponential decay. The scattering coefficient is calculated by subtracting the water absorption coefficient (µ a = 0.2 mm 1 ) from the fitted attenuation coefficient. Finally, the scattering cross section σ scat is calculated by dividing the scattering coefficient by the known particle concentration. The OCT signal amplitude at the front glass-sample boundary is determined from the exponential fit by extending the fitting line to zero depth (see Fig. 1). Zero depth was determined from the crossing of the OCT signal with the vertical drop line at half height. With this method the OCT signal magnitude is determined on an absolute scale in [mW 1/2 ] units. This is indicated in Fig. 1 on the right hand side scale for a measurement of a scattering medium for 400 nm diameter particles. The peak of the backscattered power in the heterodyne (OCT) signal is of the order of 90 picoWatts. Fig. 1. OCT measurement of backscattered power from a suspension of 400 nm diameter polystyrene microspheres in water. The vertical scale on the left is converted to absolute units of square root of power using the power calibration and is indicated on the right. The attenuation coefficient is determined from the single exponential decay fit (solid red line); the amplitude of the OCT signal from the sample surface is determined by extending the exponential fit to zero depth (dashed red line). Zero depth is indicated (dashed blue line).

Scattering samples
Polystyrene microspheres (Thermo Scientific, USA) are used as scatterers, certified by NIST traceable procedures, with the concentration of the sample calculated based on the used dilution (1 wt.% concentration). The microspheres have mean diameters of 203 ± 5, 400 ± 9, 596 ± 6, 799 ± 9, and 994 ± 10 nm and size distribution standard deviations of 4.7, 7.3, 7.7, 4.8 and 10 nm, respectively. Mie calculations are performed based on mean diameters to calculate the scattering cross section and phase function of the polystyrene particles. The refractive index of water (n water =1.32) and polystyrene are used (n polyst =1.57 [18]) as input. From Mie calculations, the scattering anisotropy, i.e. the g of these microspheres, is calculated to be: 0.07, 0.29, 0.62, 0.73, and 0.81 for increasing sphere diameter. Also, from the calculated phase functions, the scattering efficiency in the backscattering direction p NA is calculated by integrating the phase function over the NA of the sample arm focusing lens.
Measurements performed on Intralipid samples are used as an example of our technique to a polydisperse medium. Intralipid is an aqueous suspension of polydisperse lipid droplets, which is often used as a tissue phantom for optical measurements. Recently we showed that for high Intralipid concentrations multiple and dependent scattering effects play an important role and lead to nonlinear changes in the scattering coefficient with concentration [11]. To avoid these effects, measurements with low concentrations Intralipid (0.63, 1.25 and 2.5 wt.%) are performed. The samples are prepared by dilution of a single batch of 20 wt.% Intralipid (Fresenius-Kabi) with deionized water. The refractive indexes used for Mie calculations of Intralipid are: n water =1.32 for water; n lip =1.46 for lipid droplets [19].
To correlate our OCT measurements to the size of the scatterers in the polydisperse Intralipid suspension, transmission electron microscopy (TEM) measurements are performed to determine the size distribution of the scatterers in Intralipid. For TEM imaging the samples were cryoprotected with glycerin and frozen in liquid ethane at the temperature of liquid nitrogen. The samples were replicated in a vacuum better than 3×10 7 mbar at a temperature of 120°C. A platinum layer (2 nm) was evaporated at an angle of 45° and carbon (20 nm) was deposited at an angle of 90°. The replicas were cleaned with house hold bleach and collected on 300 mesh copper grids. The replicas were then imaged with a FEI T2 electron microscope at a magnification of 26500. Images were collected with a SIS Velata camera in a 2048 × 2048 format. QWin-pro (Leica) software was used to determine the diameters. In total, 2019 particle diameters were measured.

Polystyrene microspheres
From an OCT measurement, like the measurement shown in Fig. 1, the scattering coefficient for a solution of particles is determined. Figure 2(a) shows the scattering cross section of polystyrene microspheres for the different diameters obtained from the scattering coefficient. The experimental results are compared to Mie calculations and good agreement is observed (typical error is within 10%). Consequently, it can be assumed that multiple scattering effects are negligible and that the single scattering model is valid for a description of the OCT signal. From the data in Fig. 2(a) and the measured OCT magnitude at the interface, p NA is determined using Eq. (6). Measurements of p NA for all diameters are used to calculate the average heterodyne intensity back-coupling efficiency Q = 0.26 ± 0.04, which is used to compare p NA to Mie calculations in Fig. 2(b). The experimental points match the calculated values reasonably well and the oscillations in p NA due to the diameter dependent lobe structure of the backscattering efficiency are clearly observable.

Intralipid
The phase function in the backscattering direction p NA can be used to estimate the average size of the scatterers and, consequently, g. To determine the average diameter, the crossing point of the horizontal line through the experimental point of p NA and the calculated curve of p NA for each scatterer size has to be found. The obtained crossing point specifies the average particle diameter (Fig. 2(b)). Finally, the average particle diameter corresponds to a scattering anisotropy g (Fig. 2(c)).
To demonstrate our method on non-calibrated samples we apply it to Intralipid as a polydisperse suspension of scatterers. Following the same procedure as for polystyrene microspheres, p NA of Intralipid is determined. For low particle concentration, p NA is independent of the concentration of scatterers and the p NA values for all measured Intralipid concentrations (1.30·10 4 , 1.26·10 4 and 1.11·10 4 for 0.63, 1.25 and 2.5 wt.% Intralipid correspondingly) are averaged. The resulting value (p NA = (1.22 ± 0.21)·10 4 ) is plotted in Fig. 2(b) (although the refractive indexes of polystyrene spheres and Intralipid droplets are different, this has negligible influence on the calculated p NA ). From Fig. 2(b), the crossing point of the line of the calculated p NA of polystyrene microspheres and the measured p NA of Intralipid corresponds to a particle diameter of 438 ± 21 nm. Consequently, from Fig. 2(c) where we calculate g for varying particle diameter, this particle diameter corresponds to the scattering anisotropy g = 0.35 ± 0.03. This value agrees very well with previously reported measurements of g for Intralipid at 1300 nm g = 0.32 ± 0.07 [20]. Mie calculations of g versus particle diameter at 1300 nm for polystyrene microspheres. Arrows show the average particle diameter and the g of Intralipid. Dotted lines indicate the limits of our method for particle diameter and scattering anisotropy (grey zoneregion of applicability).
To validate our experimental results, Mie theory calculations are performed based on the Intralipid particle size distribution measured by TEM, see Fig. 3. The fraction of particles f vs. particle diameter d is given in the Table 1. The mean particle diameter is calculated to be 214 nm. Based on the particle size distribution in Table 1, the average scattering phase function p NA of Intralipid in the backscattering direction is calculated using [21]

Method
With the demonstrated procedure it is possible to determine g for monodisperse polystyrene particles with diameters smaller than ~500 nm (at an OCT operating wavelength of 1300 nm). With increasing particle diameter oscillations of the scattering phase function become present and cause the measured value of p NA to coincide with a multitude of diameters, which makes a unique determination of g impossible (Fig. 2(b)). However, using additional measurements (e.g. at different wavelengths [22] or different angles [9]) or using additional information (e.g. the particle concentration) the average size of the scatterers and the scattering anisotropy g also can be determined for particles larger than ~500 nm.
However, even without the possibility to perform quantitative measurements, relative changes in the scattering phase function can be used to detect changes in tissue morphology and/or act as a contrast enhancement mechanism in OCT images. Furthermore, in contrast with the backscattering coefficient, the scattering phase function does not depend on the concentration of scatterers and therefore can distinguish between changes in scattering due to refractive index contrast changes or due to changes in the scatterer size.
The accuracy in the determination of p NA and g depends on a large number of factors. As can be seen from Eq. (6), μ s , l c , NA and Q are parameters that need to be accurately determined to calculate p NA . The errors in all these parameters are summarized in the determination of the heterodyne intensity back-coupling efficiency Q, which has an error of 15%. Any systematic error in Q and its relation to the theoretical value needs to be investigated (see section 4.3) as well as methods to reduce the error in Q.

Results
The physical size distribution for Intralipid (Fig. 3) that we obtained shows a size distribution with a decreasing number of particles with increasing size. The average physical diameter is 214 nm, which is larger than 97 nm, reported by van Staveren et al. [19], which was measured for a different weight fraction solution and of Intralipid from a different manufacturer. Also our average size is larger than 123 nm, as was measured by Michels et al. [23], however, that value was obtained with an indirect method using integrating spheres measurements correlated with Monte Carlo simulations.
For polydisperse particle solutions the described method results in values for the size d and scattering anisotropy g that are averaged over the size distribution. However, since these properties are measured using an optical technique instead of a direct size measurement technique, the averages tend to underestimate the contribution of small particles since the scattering strength decreases strongly with decreasing size. For Intralipid, for example, the difference between the average physical diameter of the scatterers (214 nm) and the optically measured average scatterer diameter (438 nm) demonstrates the strong impact of the size determination method being used (TEM versus OCT). This is a well known phenomenon in, for example, dynamic light scattering where the different sizes are weighted with their respective scattering strength.
Like for the average diameter, also for g the described method is more sensitive to larger diameters [21]. However, since g is a parameter describing diffusive light transport, its optically weighted value is the only relevant parameter. The experimental value we obtain for Intralipid (g = 0.35 ± 0.03) is in good agreement with the value obtained using diffuse reflection spectroscopy (0.32 ± 0.07). In comparison to g determined from the Intralipid particle size distribution measured with TEM (g = 0.40), the value we obtained from our OCT measurements is lower. This difference can be explained by the underestimation of the small diameter particles in the TEM measurements.

Efficiency factor Q
The heterodyne detection efficiency factor Q for small scatterers is of fundamental interest in OCT [6,13,24] and we identified two major reasons that cause a reduction of the heterodyne detection efficiency Q for a homogenous distribution of scatterers compared to a mirror to values less than 1. First, the OCT magnitude for a sum of scatterers depends on the random axial position of the scatterers in the sample (causing speckle). As a result, the envelope of the OCT signal (e.g. the mean value of the speckle pattern) is smaller than the addition of all individual particle envelopes. Second, the OCT backscatter efficiency for small particles depends on lateral offset of the scatterer in the focused beam. The spherical wave originating from a particle is less efficiently coupled back to the single mode fiber compared to the planar wavefront reflected from the mirror. A full quantification of these effects is subject to future work.

Conclusions
In conclusion, we have presented experimental measurements of the scattering phase function in the backward direction p NA by OCT and demonstrated the feasibility of such measurements to determine the scattering anisotropy g. The proposed method is well applicable to scatterers with diameters below ~500 nm (for measurements at 1300 nm). We applied this method to a polydisperse Intralipid solution and obtain g = 0.35 ± 0.03, which is close to that found in literature.