Detecting structural information of scatterers using spatial frequency domain imaging.

We demonstrate optical phantom experiments on the phase function parameter γ using spatial frequency domain imaging. The incorporation of two different types of scattering particles allows for control of the optical phantoms’ microscopic scattering properties. By laterally structuring areas with either TiO2 or Al2O3 scattering particles, we were able to obtain almost pure subdiffusive scattering contrast in a single optical phantom. Optical parameter mapping was then achieved using an analytical radiative transfer model revealing the microscopic structural contrast on a macroscopic field of view. As part of our study, we explain several correction and referencing techniques for high spatial frequency analysis and experimentally study the sampling depth of the subdiffusive parameter γ.


Introduction
The analysis of backscattered light from biological tissue can reveal extensive and valuable information on tissue metabolism and structure. More precisely, diffuse optical techniques enable determination of tissue oxygenation and perfusion states [1][2][3] and are used for tissue differentiation with respect to the physical state, 4 structural alignment, 5 and malignancy. 6 All of this information is accessible by measurement of the absorption coefficient μ a and the reduced scattering coefficient μ 0 s . With the aim to also quantify microscopic structural tissue properties, there is a growing attempt for measurement of subdiffusive light scattering, which is detected at very short distances to the point of incidence. [7][8][9][10][11][12][13][14][15][16][17][18] Based on the small number of scattering interactions, subdiffusive scattering is strongly influenced by the scattering phase function of the scattering particles. The phase function states the probability distribution for light to be redirected into the different scattering angles at each scattering interaction. Generally, this function is heavily dependent on the size, shape, and relative refractive index of the scattering particles. Therefore, subdiffusive reflectance contains additional information on the micro-and ultrastructure of refractive index fluctuations in tissue and may provide insights into the underlying size distribution of scattering particles. 10,19 As microstructural properties are often linked to the dysplastic state of tissue, 20 there is the aspiration to also use diffuse optical techniques for noninvasive and possibly inter-operative histopathologic differentiation and margin assessment.
In our subsequent phantom study, we apply spatial frequency domain imaging (SFDI) to well-defined optical scattering phantoms to retrieve their microscopic scattering contrast. [7][8][9] Similar to a previous study by Kanick et al., 7 we make use of high spatial frequency analysis to image this microscopic and subdiffusive scattering contrast on a macroscopic scale through quantification of the phase function parameter γ. 19,21,22 In contrast to this previous study, we demonstrate subdiffusive imaging on a single heterogeneous yet well-defined tissue phantom and also reveal microscopic contrast for the case of closely matched macroscopic absorption and scattering properties.
After a description of our experimental setup in Sec. 2, we present our theoretical evaluation approach in Sec. 3. In this context, we also provide useful instructions for experimental reference and correction mechanisms on spatial frequency domain (SFD) reflectance data. In Sec. 4, we present our subdiffusive imaging experiments and demonstrate both mapping and depth sensitivity of the phase function parameter γ. Section 5 presents our conclusions

Experimental Setup
A sketch of our SFDI setup is given in Fig. 1. The projection unit is composed of a halogen light source, which couples into a flexible light guide for subsequent wavelength selection using a filter wheel. The wheel has a selection of eight coated bandpass interference filters for peak transmission at 500, 550, 600, 650, 700, 750, 800, and 900 nm and FWHM of 10 nm (Edmund Optics). Sinusoidal intensity modulation at spatial frequencies (0 mm −1 ≤ f ≤ 0.85 mm −1 ) is achieved by a digital micromirror device (0.7 XGA VIS, Discovery 4100 Development Kit, Vialux, Germany) and obliquely projected onto the sample at an angle of 35 deg with digital precorrection for contortion effects. Diffuse reflectance is captured from a 38 mm × 38 mm imaging area by a cooled CCD camera (QSI 640i, Quantum Scientific Imaging) at a numerical detection aperture of 0.09.
Our tissue phantoms are made of transparent epoxy resin with suspended titanium dioxide (TiO 2 , Larit reinweiß, Lange & Ritter, Germany) or aluminum oxide (Al 2 O 3 , T530SG, Bassermann Minerals GmbH, Germany) scattering particles. A detailed description of the phantom fabrication process and the optical and physical properties of the epoxy material, as well as further analysis on the employed scattering particles, can be found in Ref. 23.
To probe the subdiffusive differences of the two particle types, two pure phantoms were prepared with either TiO 2 or Al 2 O 3 scattering particles. In order to closely match the reduced scattering value in both phantoms, the weight concentration of Al 2 O 3 particles in the first phantom was ∼40 times larger than that of the TiO 2 particles in the second phantom. The large difference in weight concentration is based on the much lower scattering efficiency and higher anisotropy of the Al 2 O 3 scatterers owing to their larger diameter d and lower refractive index n. According to data from the distributing company, the median diameter of the Al 2 O 3 particles is d 50;Al 2 O 3 ¼ 1.5 AE 0.4 μm, while we determined the mean diameter of the TiO 2 particles to bed TiO 2 ¼ 250 AE 30 nm using scanning electron microscopy. A similarly large difference exists for the refractive index values n with n TiO 2 ¼ 2.7 AE 0.2 24 and n Al 2 O 3 ≈ 1.77. 25 As is discussed further in Ref. 23, both particle types exhibit variations in shape and often come in angled and elongated nonspherical forms. Phase function modeling using the Mie theory is, therefore, inappropriate.
After suspension of the scattering particles in the liquid epoxy resin and subsequent solidification, all phantoms were polished using a series of increasingly fine sandpapers (grit P160-P2500). The smoothness of the phantoms' surfaces was required to reduce the effect of specular surface reflections. While correcting for this effect, as will be described in Sec. 3, we chose not to use crossed linear polarizers in order for the scalar light propagation model to fully agree with our measurement setting. 26

Theory and Evaluation Procedure
The model we apply to our experimental data is a solution to the radiative transfer equation for oblique projection of sinusoidal intensity patterns. This solution was derived by Liemert and Kienle assuming a semi-infinite scattering medium and can be found in Ref. 27. A recent improvement of the solution approach is given in Ref. 28. On the basis of computation time, we limit the computational precision to the order N ¼ 17.
The presence of two very distinct scattering particle types (TiO 2 and Al 2 O 3 ) with very different sizes and refractive indices renders precise yet universal phase function modeling difficult. Using forward calculations with goniometrically measured phase function data of the two scattering particles for λ ¼ 750 nm, 23 we find that the Reynolds-McCormick scattering phase function may serve as an acceptable yet well-established model function to approximate the γ-reflectance characteristic of both particle types. 29 The employed Reynolds-McCormick model phase function takes the two input parameters g 0 and α and, allows for more variability than the related Henyey-Greenstein function, which is a special case of the former for α ¼ 0.5.
Owing to both the approximate nature of our phase function model and uncertainties in our experimental phase function data, our evaluation remains subject to a potential absolute error in γ of ∼0.1. This error could be determined by numerous forward calculations using both measured and model-based phase function data. The interested reader may also find a comprehensive analysis on the influence of the model phase function on the quantification of the subdiffusive parameter γ in our recent study in Ref. 22.
To best approximate the subdiffusive scattering properties of our scattering particles, we set the input parameter g 0 ¼ 0.9 in our Reynolds-McCormick model function and use the following empirical relationship for the input parameter α to allow for fitting of γ independent of any other phase function parameter.
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 4 7 7 This equation is based on a large set of forward calculations and a high accuracy fit (R 2 ¼ 0.9999) of the αðγÞ relationship.
To enable rapid derivation of absorption, reduced scattering, and the subdiffusive phase function parameter γ, a precomputed look-up table of SFD reflectance values is used by the trustregion-reflective fit algorithm (MATLAB ® , MathWorks). The error introduced by the use of the discretized look-up table is negligibly small as could be confirmed by selective evaluation without precomputed values. The four dimensional look-up table contains values for the 12 employed spatial frequencies in the range 0 mm −1 ≤ f ≤ 0.85 mm −1 for many possible combinations of μ a , μ 0 s , and γ. Each of these dimensions is discretized into 40 parameter steps with a logarithmic spacing of μ a in the range 5 × 10 −4 mm −1 ≤ μ a ≤ 5 × 10 −2 mm −1 and linear spacing of both μ 0 s and γ with 0.5 mm −1 ≤ μ 0 s ≤ 2.5 mm −1 and The derivation of the phantoms' optical properties is not based on raw experimental imaging data, but follows a series of computational postprocessing steps. Some of these steps are of special importance for high spatial frequency analysis. In the following, we will elaborate on three postprocessing steps that account for absolute referencing of the reflectance intensity and correction of the imaging data for system related defocusing and specular surface reflections.
The first step is the common and well-known demodulation principle. In this step, the three evenly phase-shifted images R 1 ðx; yÞ, R 2 ðx; yÞ, and R 3 ðx; yÞ are turned into the actual SFD reflectance image R AC;raw ðx; yÞ. E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 1 1 1 R AC;raw still represents a raw image with pixel values measured in counts. In the next step, these raw values need to be transferred into actual reflectance intensities, which are independent of both the intensity and spatial heterogeneity of the projection system and the spatial heterogeneity of the detection system. This step is often approached by measurement of a spatially homogeneous reference phantom with very well-known optical properties or reflectance intensities. In so doing, measurement inaccuracies from other imaging modalities may be transferred and the underlying scattering phase function of the reference phantom has to be known precisely.
In order not to rely on other imaging modalities and possible phase function uncertainties, we follow a slightly different approach by measuring the optical properties of a very homogeneous and intensely scattering epoxy phantom (μ 0 s ≈ 15 mm −1 ). The only assumptions we make on this phantom are the spatial homogeneity of its optical properties along with its predetermined refractive index values. 23 We then determine the phantom's optical properties (i.e., μ a , μ 0 s , and γ) at every imaging point using the product of light source intensity and detection sensitivity as an additional multiplicative fit variable a. In a next evaluation step, the derived optical properties are averaged and assumed constant for all imaging points, which leaves the factor aðx; yÞ as a single and robust fit variable in a second reevaluation of the phantom. This fit variable aðx; yÞ quantifies the ratio between R AC;raw and the actual reflectance value R AC as predicted by the theoretical model. E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 6 3 ; 4 5 5 R AC ðx; y; fÞ ¼ aðx; yÞR AC;raw ðx; y; fÞ: (3) All subsequent measurements are multiplied by the obtained map aðx; yÞ, serving as an intensity reference that quantifies both the light source intensity and the detection efficiency at the same time. Note that aðx; yÞ is independent of spatial frequency. One further postprocessing step corrects for optical defocusing in the spatial frequency projections, which is especially required for high spatial frequencies and oblique projection geometries. We achieve this by measuring the SFD reflectance from an anodized plate of aluminum. Due to the absence of volume scattering, the reflectance from the aluminum originates solely from the metal surface and the very thin anodization layer, which is typically <25 μm thick. Using this reflectance measurement, the amount of defocusing at a spatial frequency f can be quantified by the ratio bðx; y; fÞ, where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 6 3 ; 2 7 0 bðx; y; fÞ ¼ R DC;Al ðx; yÞ R AC;Al ðx; y; fÞ ; and R DC;Al corresponds to the aluminum reflectance at spatial frequency zero, which is commonly obtained by averaging R 1 , R 2 , and R 3 at any spatial frequency. For very sharp projections and low spatial frequencies, b ≈ 1 and no correction is required. For high spatial frequencies, b may become much larger based on a decrease in R DC;Al . We perform postcorrection for projection defocusing by multiplying R AC ðx; y; fÞ of all subsequent measurements with bðx; y; fÞ.
In a last correction step, we take into account specular surface reflection. Surface roughness related specular reflections have a potentially very strong impact on high spatial frequency analysis. The use of crossed linear polarizers can diminish this effect while complicating the theoretical modeling of the underlying light propagation. 26 Through the use of a nonscattering transparent epoxy phantom with closely matched surface roughness, we were able to quantify the specular reflections and to subtract them from our actual phase measurements. In addition to this correction step, we chose a rather large projection angle in our measurement setup to further reduce the influence of specular surface reflections.
At the end of the above postprocessing steps, every computed reflectance image with a pixel resolution of 2048 × 2048 pixels was binned with bin sizes of 8 × 8 pixels to reduce the overall data volume and to enhance the signal-to-noise ratio. Especially for high spatial frequency analysis, it is recommended to perform this binning step only after computation of the AC-reflectance value, as binning may otherwise add to the demodulation of the phase data.
Both our experimental setup as well as our evaluation approach have been repeatedly verified by comparative studies employing various other measurement modalities, such as spatially and time resolved reflectance spectroscopy on a large set of well-defined tissue phantoms.

Experimental Results
Our initial experiments relate to two pure scattering phantoms with either TiO 2 or Al 2 O 3 as the suspended scattering particles. These two homogeneous phantoms are hardly distinguishable to the naked eye yet have some differences in their reduced scattering and absorption properties. As previously mentioned, both particle types differ strongly in their scattering phase function based on the differences in particle diameter and refractive index. The two phantoms and all further phantoms were characterized with respect to μ a , μ 0 s , and γ using SFDI at 12 linearly spaced frequencies in the range 0 mm −1 ≤ f ≤ 0.85 mm −1 . Figures 2(a)-2(c) show the result of our measurements on the pure scattering phantoms averaged over five phantom positions versus wavelength λ. The given error bars correspond to the standard deviation of the positional average. Figure 2(a) reveals similar absorption values in the visible to near-infrared wavelength range for both tissue phantoms. The absorption is found to be very small in both phantoms and is only slightly increased for the Al 2 O 3 phantom based on the intrinsic absorption of Al 2 O 3 scatterers and its higher weight concentration. Figure 2(b) shows the different wavelength dependences of μ 0 s for both phantoms. The two dashed curves in Fig. 2(b) are a fit of the power law μ 0 s ðλÞ ∝ λ −α to the reduced scattering values yielding α Al 2 O 3 ¼ 0.58 and α TiO 2 ¼ 1. 15. The stronger decline of μ 0 s ðλÞ for the TiO 2 phantom relates to the smaller size of the TiO 2 scattering particles and their stronger refractive index wavelength dependence. Owing to the different wavelength dependence of the reduced scattering value of both pure scattering phantoms, an almost matched reduced scattering value can be found at the wavelength λ ¼ 750 nm in Fig. 2(b). In addition to the overall small differences in absorption, this wavelength is well suited for our next step, in which we aim to perform subdiffusive imaging in the near absence of macroscopic absorption and scattering contrast in a combined heterogeneous phantom.
Relating to the wavelength dependence of scattering, the microstructural differences of the two pure phantoms become especially apparent in Fig. 2(c), showing the obtained γ values for both phantoms. The larger size of Al 2 O 3 particles relates to both more forward peaked and reduced high angle scattering leading to less subdiffusive backscattering and a larger γ value. We compare our results for γ (dashed curves) to those obtained by goniometric measurements on the same scattering particles in a similar refractive index medium. These comparative values are given by the two solid curves in Fig. 2(c), with the shaded areas representing their measurement errors. 23,30 In spite of common features in the compared data curves, we find absolute deviations for both particle types. For TiO 2 , we attribute the systematic deviation of Δγ ≈ 0.1 to the approximate nature of our phase function model and residual specular surface reflections, which lower the measured γ values. The larger deviation for the Al 2 O 3 particles may have its cause in differences in the effective particle size distributions due to slow sedimentation of the bigger-sized Al 2 O 3 particles during phantom curing and the surface insensitivity of the goniometric measurement approach.
Using the derived optical properties of the two pure scattering phantoms from Figs. 2(a)-2(c), it is possible to create a laterally heterogeneous phantom composed of both phantom materials with only microscopic subdiffusive scattering contrast at λ ¼ 750 nm. We achieved this by engraving the name of our institute into the back side of the Al 2 O 3 phantom with a milling machine and a rounded milling head of 4 mm diameter. Using a syringe, the on average 2-mm-deep inscription was then filled with TiO 2 suspended epoxy resin with optical properties and TiO 2 particle concentration equal to that of the TiO 2 reference phantom [see Fig. 2(d)]. After hardening and repeated surface polishing, the lateral structuring of the so obtained phantom is hardly recognizable with the eye based on the small macroscopic scattering and absorption contrast in the visible range. Indeed, according to Fig. 2(b), the composed phantom has hardly any macroscopic scattering contrast at λ ¼ 750 nm. We used this wavelength to perform SFDI on the heterogeneous phantom in an almost complete absence of μ a and μ 0 s contrast. Figures 3(a)-3(e) give SFD reflectance images [i.e., grayscale maps of R AC ðfÞ] at various spatial frequencies for λ ¼ 750 nm and illustrate the increasing microscopic scattering contrast with spatial frequency. These reflectance images were scaled by our intensity reference and corrected for image blurring and specular surface effects, as described previously.
The increased subdiffusive backscattering of TiO 2 scatterers gives rise to a strong contrast at spatial frequency f ¼ 0.35 mm −1 . Radiative transfer model based analysis finds the largest absolute derivative of R AC with respect to γ near this spatial frequency and reveals that only five percent of the SFD reflectance signal is based on photons traveling deeper than 1 mm for the given phantom optical properties. 22 We denote this penetration depth by d 5 and provide all corresponding depth values in the caption of Fig. 3. To best compare the achieved microscopic imaging contrast of Figs. 3(a)-3(e), the gray-scale coloring was adjusted for every image such that the ratio of pixel values corresponding to the brightest and darkest color is constant among all reflectance images.
A very interesting effect can be observed in Figs. 3(c) and 3(d). In these images, the contrast of the engraved letters mostly stems from a bright contour at the inner border of the engraving. This contour is ∼0.5 mm wide and explicable by a slight bias for subdiffusive light in the boundary area to reach the phantom surface on the TiO 2 side.  By fitting of our radiative transfer model to the experimental data, optical parameter maps for μ a , μ 0 s , and γ were obtained and are given in the top row of Fig. 4. As expected from the optical properties of the pure reference phantoms [Figs. 2(a)-2(c)], the laterally structured phantom has very little absorption and reduced scattering contrast at λ ¼ 750 nm, as can be seen in the left and central parameter maps of Fig. 4. The apparent change in optical properties on the lower and upper left corners of the images is a boundary artifact based on the finite projection area, which extends the imaging field of view by only 10 to 20 mm to each side. 31 The third parameter map in Fig. 4 depicts the strong γ contrast relating to the different types of scattering particles and, thus, demonstrates the capability of high-frequency analysis for macroscopic imaging of quantified microscopic scattering properties.
Three line profiles in the bottom row of Fig. 4 further quantify the obtained contrast in the three parameter maps. The line profiles correspond to the single pixel lines indicated by the dashed lines in the parameter images above.
The quantification of γ is noticeably impaired by surface roughness as can be seen by the scratches in the γ map. The apparent increase in γ at the left end of the imaging area is an artifact caused by residual specular surface reflections in connection with the oblique projection geometry.
The optical properties displayed in Fig. 4 agree well with those measured for the pure reference phantoms (Fig. 2), even though our model based assumptions for a semi-infinite turbid media do not reflect the phantom's lateral and axial structuring. The absorption map together with its line profile in Fig. 4 depict an almost complete absence of absorption contrast. This corresponds to the relatively long optical path lengths required for such absorption intensities to have a meaningful impact on the reflectance signal.
The parameter sensitivity toward μ 0 s is mostly given by intermediate optical path lengths (i.e., intermediate spatial frequencies), which have dimensions similar to the depth and width of our engraving. Therefore, we observe the obtained μ 0 s contrast at wavelengths other than λ ¼ 750 nm to be diminished (not shown) as compared to the expected contrast based on the homogeneous phantoms. This also reflects the model's assumption of a semi-infinite and homogeneous sample geometry and its inability to correctly describe the remaining macroscopic scattering heterogeneity. 1 On the contrary, we always find the obtained γ maps also for wavelengths other than λ ¼ 750 nm to retain most of the contrast anticipated from our reference measurements. This is because the corresponding short photon paths are mostly limited to the inside or the outside of the engraved region.
When using γ maps for the study and assessment of microscopic tissue structure, it is important to know the corresponding sampling depth of the obtained γ values. To answer this question, we took an experimental approach by adding additional scattering layers on top of the engraved phantom [see Fig. 5(a)]. These layers are made of the same epoxy phantom material and incorporate TiO 2 scatterers in concentrations very similar to that of the phantom engraving. By iteratively adding layers of increasing thickness onto the phantom, the γ contrast obtained from further reflectance measurements is more and more reduced and eventually disappears. One thereby finds the depth regime that mostly accounts for the quantification of γ. Figure 5 shows the obtained γ maps at λ ¼ 700 nm for no added layer [ Fig. 5 layer are coated with a microscopy immersion index-matching fluid to reduce Fresnel reflections both on the surface and between the layer and phantom. The fluid greatly reduces the noise in the γ maps but gives rise to a small interference effect as seen by the bow pattern in Figs. 5(b)-5(e). As expected, the γ contrast of the inscription diminishes with increasing effective layer thickness and is even canceled at the layer thickness of dμ 0 s ¼ 1.18 as can be seen in Fig. 5(e). When comparing the coloring (i.e., the γ values) within the "ilm"inscription of Figs. 5(b) and 5(c), a slight decrease in γ value can be observed, which we attribute to the added interface between the cover layer and phantom. In spite of the used index-matching fluid, the added interface and the finite thickness of the oil layer have a small influence on the backscattering intensity. However, it can be observed in Figs. 5(c)-5(e) that γ remains almost constant in the engraved area independent of the surface layer thickness, as the surface layer is composed of the same scattering particles. By contrast, the lower subdiffusive backscattering of the outer Al 2 O 3 particles gets more and more masked by the TiO 2 cover layer and is even hidden at the effective thickness of dμ 0 s ¼ 1.18 (Fig. 5(e)). As a conclusion from this cover layer experiment, we find the sampling depth of γ to lie in the range of 0.78 < dμ 0 s < 1.18. This finding is consistent with the basic conception for ballistic light to lose its directionality after a propagation distance of roughly d ¼ 1∕μ 0 s . An even smaller sampling depth of γ images may be achieved by sole consideration of very high spatial frequencies.

Conclusions
Through our phantom study, we have not only demonstrated the feasibility to fabricate laterally structured microscopic scattering contrast in tissue simulating phantoms, but also described its measurement using high spatial frequency analysis in SFDI. The use of high spatial frequency analysis allows for mapping of the phase function parameter γ and can, thus, provide macroscopic images of the underlying microscopic structure. In addition to the control of the phantoms' macroscopic absorption and reduced scattering values, the use of two different scattering agents allowed us to design distinct and well-defined subdiffusive scattering properties. By measurement of γ, these subdiffusive scattering properties, which ultimately relate to the size and the refractive index of the scattering particles, could be quantified.
The presented measurement approach has potential application for in vivo tissue analysis, where the quantification of subdiffusive scattering may provide important insights into the tissue microstructure, which is otherwise only accessible through histology. Dysplastic tissue alterations often come along with major changes in the cellular structure, and it is highly desirable to be able to detect these changes without the need for tissue biopsies. Subdiffusive reflectance imaging represents one promising approach toward such in situ measurements.
In view of its potential biomedical application, the strong influence of specular surface reflections even for our almost ideally flat phantoms underscores the need to find strategies to correct for these imaging artifacts in nonideal imaging situations.
In a previous study, we reported on the importance of correct phase function modeling for precise γ quantification. 22 For this phantom study, the Reynolds-McCormick phase function was used to model the subdiffusive reflectance characteristics of the employed scattering particles. In view of similar studies on biological tissue, a fractal size distribution of Mie scatterers may provide a useful model for microscopic tissue scattering. In this case, the retrieved γ values could favorably relate to the fractal dimension and, thus, to the size distribution of the underlying scatterers. 19,[32][33][34] With respect to tissue histology, this information may allow for valuable insights into the microscopic tissue structure averaged over the tissue surface to a depth of about 1∕μ 0 s .