Open Access
Add to BibSonomyAbstract
In this paper a new approach to improve contrast in optical subsurface imaging is presented. The method is based on time-resolved reflectance and selection of ballistic photons using polarization gating. Numerical studies with a statistical Monte Carlo method also reveal that weakly scattered diffuse photons can be eliminated by employing a small aperture and that the contrast improvement strongly depends on the single-scattering phase function. A possible experimental setup is discussed in the conclusions.
© 2010 Optical Society of America
1. Introduction
In diffuse optical imaging typically an unpolarized light source is used for time-resolved measurements. Moreover, for technical reasons the source is separated from the detector by at least a distance ρ = 2 – 4cm. Recently, it has been reported [1] that a contrast improvement can be achieved if ρ = 0 (null source-detector separation). The drawback of such an experimental setup is the presence of early photons, which in the case of time-correlated single-photon counting (the most common technique for diffuse imaging in the time domain) limits the counting rate of the employed electronics below saturation threshold [1]. Finally, a compromise was proposed, i.e. to consider small but not null source-detector separation in order to avoid blinding due to early photons, but still improve contrast and localization. On the other hand, polarization gating was considered [2–5]. Linearly polarized light was used in [2] to improve contrast in subsurface imaging by subtracting the image acquired with cross-polarized light from the one obtained with co-polarized light. Note that the photons which penetrate deeper into the tissue and thus experience more scattering events tend to have perpendicular polarization states. Spectral and polarization discrimination of backscattered light were employed in [3] to perform subsurface imaging. In [4, 5], both circularly and linearly polarized light was considered to detect a mirror embedded in scattering media composed of either small or large scattering particles. Circularly polarized light is preferred, if a scattering medium is composed of large scattering particles (g = 0.911 was used). Of course a mirror is an ideal target for subsurface imaging, since the helicity of the reflected light flips. This results in a nice signature compared to light backscattered from the background medium. In [6], reflecting, scattering and absorbing targets were probed using linearly and circularly polarized light in order to examine optical target properties in subsurface imaging. Monte Carlo simulations were used in [7] to investigate effects of scattering and absorption coefficients (μs and μa), anisotropy factor g and radius of the detector on the penetration depth of polarization gating. It was shown that with polarization gating one can detect photons, which penetrated several mean free path lengths into the medium. A polarization memory of circularly polarized light pulses results from the successive near-forward scattering events enhanced by a large anisotropy factor as shown in [8]. These near-forward scattering events maintain the same circular polarization. In [9] an analytic solution of the time-dependent vector radiative transfer equation is used to compute backscattering of circularly polarized light and it can be expected that circularly polarized light with reversed helicity dominates only when the detector is located close to the source. In [10] a time-independent Fokker-Planck approximation of the radiative transfer equation was used to compute light backscattered from forward scattering media. It was shown that the backscattered ring for circularly polarized incident light is mostly composed of light with the same helicity as the incident beam. Furthermore, in [11] it was shown that the observed ring structure of backscattered light is time dependent.
In this paper, we used polarized light to distinguish between diffuse and ballistic photons (photons with straight trajectories) in order to improve the signal-to-noise ratio. In all numerical experiments the light is detected in the time domain, which in the case of ballistic photons would require the use of a streak camera (high temporal resolution, e.g. 2ps [4]). Note that ballistic photons essentially experience one backward and many forward scattering events and therefore, unlike diffuse photons, their helicity is flipped. Furthermore, ballistic photons leaving the scattering medium tend to propagate in opposite direction than the incident photons (with little variance). Thus, a small acceptance angle of the detector effectively excludes most diffuse photons. It seems likely that such an experimental setup, which requires more expensive equipment (streak camera), does not experience the saturation problem discussed in [1]. For the studies discussed in this paper we considered circularly polarized laser light, since the considered single scattering phase functions are highly anisotropic (g > 0.7). For each test case the intensity, the V-component of the Stokes vector, the reached depth and the contrast are shown as functions of time-of-flight.
In section 2, the governing equations and some basic theory are presented, in section 3 the simulation setup and materials are described, computational studies are presented in section 4 and conclusions are given in section 5.
2. Governing equations
For all studies presented in this paper, the transport theory is considered, where the vector Boltzmann equation
for the Stokes vector I(x, s, λ, t) is solved with the Monte Carlo method implemented in our code Scatter3D [?, 12–15]. In Eq. (1) the position in physical space is denoted by x, s is the propagation direction, λ the wavelength and t the time. The vectors s and s′ form a so-called scattering plane, for which the single-scattering Mueller matrix M(s, s′) is defined. Extinction and scattering coefficients are denoted by γt and γs, respectively. A single scattering event requires rotation of the Stokes vector I = [I Q U V ]T into the scattering plane and its multiplication with the 4 × 4 single-scattering Mueller matrix (see Fig. 1). The components of the Stokes vector are defined as where El and Er are the complex parts of the phasors El = alexp(−iωt – iδ0) and Er = arexp(−iωt – iδ0 – iδ). Here, El and Er are the parallel and perpendicular components of the electric field vector with respect to the scattering plane. The phase difference between Er and El is denoted by δ, al = |al| and ar = |ar| are the amplitudes of El and Er, ω = 2πc/λ is the frequency and δ0 is the phase of El at t=0 (c is the speed of light).The change of the electric field vector has to be computed at every scattering event. A complete description of the method can been found in [16], but for completeness the most important part is shown below. At every scattering event, El = alexp(−iωt – iδ0) and Er = arexp(−iωt – iδ0 – iδ) are changed according to
where ϕ (azimuth angle) is the clockwise rotation angle around the pre-scattering propagation trajectory and θ is the polar scattering angle defined as shown in Fig. 1. Note that is the electric field vector before scattering. The scattering amplitude functions S1 and S2 in Eq. (3) are computed by the Lorenz-Mie theory and the detailed description of their computation can be found in [14]. Note that the light intensity I has to be preserved (without absorption), i.e. . Normalization factor F(θ, ϕ) in Eq. (3) is defined in [16]. Note that the Stokes vector can be computed from El and Er at any time for any scattering plane.3. Simulation setup and materials
As show in Fig. 2, the circular detector with radius R = 2.5mm is coaxial with the laser beam and source-detector separation is ρ = 0. Outgoing photons with all deflection angles from the surface normal were considered and in all test cases the center of the considered spherical inhomogeneity was placed on the line along the laser beam at different depths (z = 3mm and z = 5mm).
Fig. 2 Simulation setup with semi-infinite background medium and spherical inhomogeneity with diameter D = l2 – l1 placed at depth z and the center on the line of the incident laser beam. Source and detector are coaxial (ρ = 0). The trajectories t1 and t2 represent two possible paths of detected photons using polarization filtering.
Two different scattering media labeled as phantom 1 and phantom 2 were employed. Both phantoms can be considered as semi-infinite and are composed of scattering particles with diameter d = 1000nm embedded in silicon with nsi = 1.404. The particles in phantom 1 are TiO spheres with npa = 2.609 designed to mimic human tissue and in phantom 2 polystyrene spheres with npa = 1.59 are considered. The illumination is a near-infrared laser light with λ = 780nm perpendicularly incident on the scattering medium (propagation along the positive z-direction). The Stokes vector [see Eq. (2)] of incident light is I = [1 0 0 1]T. The background medium has the same scattering coefficient as the spherical inhomogeneity (μs = 10cm−1), but different absorption coefficients, i.e. μa,bg = 0.2cm−1 for the background medium and μa,target = 0.5cm−1 for the inhomogeneity. Matched refractive indices of phantoms and surrounding medium were used. To compute the single-scattering phase functions, which have anisotropy factors of g ≈ 0.72 and g ≈ 0.92 for the phantoms 1 and 2, respectively, the Lorenz-Mie theory was employed. For an easier overview Table 1 presents all considered test cases.
Table 1. Overview of all considered test cases. Incident light is circularly polarized I = [1 0 0 1]T with λ = 780nm
4. Results
For the test cases 1 and 2 the phantoms 1 and 2 are considered, respectively, and the spherical inhomogeneity with diameter D = 1mm is located at z = 3mm. The results are shown in Fig. 3.
Fig. 3 Results for the test cases 1 and 2 computed without polarization filtering using phantoms 1 and 2, respectively; the spherical inhomogeneity of diameter D = 1mm is placed at z = 3mm. Lines marked with circles and squares represent mean values for the test cases 1 and 2, respectively, and the lines with diamonds and asterisks show the corresponding standard deviations. Sub-plots (a), (b), (c) and (d) show detected intensity, V-component of the Stokes vector, reached depth and contrast as functions of time-of-flight, respectively.
The intensity difference between the two test cases shown in Fig. 3(a) is negligible, since the absorbed energy only depends on the absorption coefficient and time spent inside the medium. On the other hand, the detected V-components of the Stokes vector shown in Fig. 3(b) differ significantly and the standard deviations are large in both cases indicating the potential for contrast improvement using polarization filtering. In Fig. 3(c) the reached depth is shown as a function of time-of-flight and it can be observed that it is larger for test case 2 compared to test case 1, which is due to the more anisotropic single-scattering phase function resulting in less corrugated photon trajectories. As a consequence, the contrast
where I and I0 are the detected intensities with and without spherical inhomogeneity, respectively, is larger for test case 2 than for test case 1 [Fig. 3(d)]. Note that the distance l1 in Fig. 2 is approximately represented by the speed of light times half the time-of-flight for which C [see Eq. (4)] in Fig. 3(d) starts to rise, while l2 is approximately equal to the speed of light times half the time-of-flight corresponding to the contrast peak. Obviously it is possible to estimate depth and size of the target from the contrast curve, while zero source-detector separation leads to the best possible localization in the xy-plane. Figure 4 shows the standard deviations of the photon trajectories (distances between scattering positions and the axis of the incident laser beam) conditional on time-of-flight and V-component of the Stokes vector.Fig. 4 (Color online) Standard deviations (in cm) of photon trajectories for test case 2 conditional on the V-component of the Stokes vector and time-of-flight.
As expected, for the same time-of-flight trajectories of photons with larger V-value tend have a larger standard deviation. Employing polarization filtering, photons with corrugated trajectories can be eliminated, which improves the signal-to-noise ratio and consequentially the contrast. Similar results were obtained for phantom 1 as well. The 3D histogram shown in Fig. 5 gives a closer look at the amount of detected photons within a certain time-of-flight and V range.
Fig. 5 (Color online) 3D histogram of detected photons in test case 2 as a function of time-of-flight and V-component.
One can observe a large photon density for a time-of-flight below 30ps and V-values above −0.5, which is consistent with the ”bump” in Fig. 3(b) at a times-of-flight around 20ps. Compared to diffuse photons, the deflection angles α from the surface normal tend to be smaller for detected photons with more or less straight trajectories. In order to see that, the outgoing angles of the detected photons (for phantom 2) are shown in Fig. 6 as a function of time-of-flight and V-value. Furthermore, from the histogram in Fig. 5 and the results in Fig. 6 it becomes obvious that a large portion of weakly scattered diffuse photons can be removed by only accepting photons with α ≤ 30°.
Fig. 6 (Color online) Angles (in degrees) between photon propagation directions and outer surface normal as a function of V-component of the Stokes vector and time-of-flight. Results are shown for test case 2.
For the following test cases 3 and 4 phantom 2 and a spherical inhomogeneity with diameter D = 1mm and center at z = 5mm is considered, while polarization filtering was only applied for test case 4, i.e. photons with V ∉ [−1, −0.7] were rejected. Results are shown in Fig. 7.
Fig. 7 Results for the test cases 3 and 4. Lines marked with circles and squares represent mean values for the cases 3 and 4, respectively, while the lines with diamonds and asterisks represent the corresponding standard deviations. The sub-plots (a), (b), (c) and (d) show the detected intensity, V-component of the Stokes vector, reached depth and contrast as functions of time-of-flight, respectively.
Again, no difference between the detected intensities can be observed in Fig. 7(a). As expected, the standard deviation of V is much smaller for test case 4 due to polarization filtering [see Fig. 7(b)]. Since the average trajectories of the polarization filtered photons (test case 4) are less corrugated, the reached depths are slightly larger for the same time-of-flight compared to test case 3 [see Fig. 7(c)]. Consequently, a better contrast can be detected for test case 4 due to improved signal-to-noise ratio [see Fig. 7(d)]. Note that the contrast curve for test case 4 declines for larger times-of-flight and eventually reaches zero. This is due to the coaxial setup of laser and detector and the V-shaped photon trajectories. For large times-of-flight the same reflection angle as that of trajectory t2 in Fig. 2 leads to elimination for deeper penetrations. The same simulation setup as for the test cases 3 and 4 was used for the test cases 5 and 6, but phantom 1 instead of phantom 2 was considered. The results are plotted in Fig. 8.
Fig. 8 Results for the test cases 5 and 6. Lines marked with circles and squares represent mean values for the test cases 5 and 6, respectively, while the lines with diamonds and asterisks show the corresponding standard deviations. The sub-plots (a), (b), (c) and (d) show intensity, V-component of the Stokes vector, reached depth and contrast as functions of time-of-flight, respectively.
For the sub-plots (a), (b) and (c) the same conclusions can be made as for the results shown in the corresponding sub-plots of Fig. 7. Comparing Figs. 7(d) and 8(d) shows that the possible improvement due to polarization filtering is larger, if phantom 1 rather than phantom 2 is employed. This is due to the different anisotropy factor g, i.e. a larger g-value leads to less corrugated trajectories and consequentially in larger contrast. On the other hand, smaller g-values result in more diffuse photon scattering and thus polarization filtering has a higher impact. Figures 9(a) and 9(b) depict the reconstructed inclusions based on the contrast curves for test cases 5 and 6, respectively. The abscissa represents the scanning direction (x-axis) and the ordinate the depth (z-axis; the dimensions of both axes are in mm). The grayscale value represents the slope of the contrast curve showing that polarization gating leads to significant contrast improvements. A slope of the contrast curves is calculated as a difference between contrast values in two neighbor time bins divided by its width which is equal to 1; i.e. time bins of 1ps were used in all presented results. Thus the grayscale values in Fig. 9 are one order of magnitude lower than the contrast values in Fig. 8(d). For larger depths (z = 10mm and D = 2mm), in phantom 1 the polarization information gets lost and no contrast improvement can be achieved, while in phantom 2 a small contrast improvement is achieved and the peak contrast value is of the same order of magnitude as in the test cases 3 and 4. Next, the contrast curves with and without polarization gating were computed using phantom 2 and inclusion with D=2mm placed at depth of 15mm. The contrast values are one order of magnitude smaller compared to those in Fig. 7 and no contrast improvement was achieved using polarization gating. However, in this simulation a low signal-to-noise ratio is reported.
Fig. 9 Sub-plots (a) and (b) depict the reconstructed inclusion of test cases 5 and 6, respectively.
5. Conclusions
In conclusion we report that a significant portion of weakly scattered diffuse photons can be removed with a detector, which has a small enough acceptance angle, i.e. diffuse photons get eliminated much more likely than ballistic ones resulting in better contrast. For larger times-of-flight (> 30ps), the contrast can be improved by employing polarization filtering, if one is able to select only that light with flipped helicity within a certain ellipticity range (which depends on the inclusion depth). One can think about a potential experimental setup similar to the one used in [4], but with a variable wave-plate in front of the detector (after beam splitter) and additional quarter wave-plate on the illumination side (before beam splitter). The quarter wave-plate creates circularly polarized incident light. Several pulses, for which the variable wave-plate is adequately adapted, could be emitted and since the V-component of the Stokes vector monotonically increases with time-of-flight, the detected signal obtained by converting elliptically polarized backscattered light into linearly polarized light can be conditioned on the penetration depth. Obviously, such an approach is slower and more expensive than simpler ones, where all light with flipped helicity is detected [4] (using a quarterwave-plate and Wollaston prism), but there may exist situations where this is justified by the higher contrast. If one knew in advance the change of the V-component as a function of time-of-flight (for a given phantom and illumination), the contrast curve could be obtained without measurements in the time domain avoiding the use of an expensive streak camera. Therefore it is sufficient to measure the change only once for each commonly used phantom and illumination wavelength and it is possible to use the obtained curves for future imaging experiments, since in both phantoms the V-value of detected ballistic-photons (photons detected using polarization gating) monotonically changes with the time-of-flight and thereby uniquely determines their penetration depth. In all our simulation studies the incident laser light was circularly polarized and the anisotropy factor g was greater than 0.7 (which is the case for most biological tissues). Note that larger anisotropy factors lead to higher contrast values, but the improvement due to polarization filtering is greater for smaller g values.
Acknowledgments
Finally we would like to acknowledge the support by the Laboratory for Media Technology at the Swiss Federal Laboratories for Materials Science and Technology (EMPA) with regard to the provided computer power, which was required for the numerical studies presented in this paper.
References and links
1. A. Torricelli, A. Pifferi, L. Spinelli, R. Cubeddu, F. Martelli, S. Del Bianco, and G. Zaccanti, “Time-Resolved Reflectance at Null Source-Detector Separation: Improving Contrast and Resolution in Diffuse Optical Imaging,” Phys. Rev. Lett. 95, 078101 (2005). [CrossRef] [PubMed]
2. S. G. Demos and R. R. Alfano, “Optical polarization imaging,” Appl. Opt. 36, 150–155 (1997). [CrossRef] [PubMed]
3. S. G. Demos, H. B. Radousky, and R. R. Alfano, “Deep subsurface imaging in tissues using spectral and polarization filtering,” Opt. Express 7, 23–28 (2000). [CrossRef] [PubMed]
4. X. Ni and R. R. Alfano, “Time-resolved backscattering of circularly and linearly polarized light in a turbid medium,” Opt. Lett. 29, 2773–2775 (2004). [CrossRef] [PubMed]
5. S. A. Kartazayeva, X. Ni, and R. R. Alfano, “Backscattering target detection in a turbid medium by use of circularly and linearly polarized light,” Opt. Lett. 30, 1168–1170 (2005). [CrossRef] [PubMed]
6. R. Nothdurft and G. Yao, “Expression of target optical properties in subsurface polarization-gated imaging,” Opt. Express 13, 4185–4195 (2005). [CrossRef] [PubMed]
7. Y. Liu, Y.L. Kim, X. Li, and V. Backman, “Investigation of depth selectivity of polarization gating for tissue characterization,” Opt. Express 13, 601–611 (2005). [CrossRef] [PubMed]
8. A.D. Kim and M. Moscoso, “Backscattering of circularly polarized pulses,” Opt. Lett. 27, 1589–1591 (2002). [CrossRef]
9. W. Cai, X. Ni, S.K. Gayen, and R. R. Alfano, “Analytical cumulant solution of the vector radiative transfer equation investigates backscattering of circularly polarized light from turbid media,” Phys. Rev. E 74, 056605 (2006). [CrossRef]
10. A. D. Kim and M. Moscoso, “Backscattering of beams by forward-peaked scattering media,” Opt. Lett. 29, 74–76 (2004). [CrossRef] [PubMed]
11. K. G. Phillips, M. Xu, S. K. Gayen, and R. R. Alfano, “Time-resolved ring structure of circularly polarized beams backscattered from forward scattering media,” Opt. Express 13, 7954–7969 (2005). [CrossRef] [PubMed]
12. P. Jenny, S. Mourad, T. Stamm, M. Vöge, and K. Simon, “Computing light statistics in heterogeneous media based on a mass weighted probability density function method,” J. Opt. Soc. Am. A 24, 2206–2219 (2007). [CrossRef]
13. M. Šormaz, T. Stamm, S. Mourad, and P. Jenny, “Stochastic modeling of light scattering with fluorescence using a Monte Carlo-based multiscale approach,” J. Opt. Soc. Am. A 26, 1403–1413 (2009). [CrossRef]
14. M. Šormaz, T. Stamm, and P. Jenny, “Stochastic modeling of polarized light scattering using a Monte Carlo-based stencil method,” J. Opt. Soc. Am. A 27, 1100–1110 (2010). [CrossRef]
15. M. Šormaz, T. Stamm, and P. Jenny, “Influence of linear birefringence in the computation of scattering phase functions,” J. Biomed. Opt. 15, 055010 (2010). [CrossRef] [PubMed]
16. M. Xu, “Electric field Monte Carlo simulation of polarized light propagating in turbid media,” Opt. Express 12, 6530–6539 (2004). [CrossRef] [PubMed]
