## 1.

## Introduction

Many tissue types, including bone, muscle, skin, and white matter in the brain, have orientated internal structures such that the degree of optical scattering depends on the direction of light propagation. The spatial-dependence of the bulk (e.g., multiple) scattering intensity in these tissues is determined by their microscopic structure. Imaging the directionality of multiple light scattering may provide a noninvasive way to interrogate tissue microscopic structure, and may aid in the early detection and diagnosis of disease.

The orientational dependence of multiple light scattering from ordered structures has been measured in biological tissues including bovine tendon,^{1} chicken breast,^{2} skin,^{3, 4} bone,^{4} muscle,^{5} and dentin.^{6} It has also been tested in tissue phantoms.^{7, 8, 9} Researchers have modeled multiple scattering from ordered structures using Monte Carlo simulations of radiative transport,^{10} as well as random walk^{11, 12} and finite element implementations^{13} in the diffusive regime. In all of these studies, the emphasis has been to determine the orientation of large homogeneous scattering regions by looking for elliptical equal intensity profiles of the backscattered or transmitted light from a localized (e.g., point) source. Greater eccentricity of the ellipse indicates a higher degree of scattering structure, and the orientation of the ellipse gives the direction of the underlying structure. Such measurements of scattering orientation give the average scattering properties over an area of a few centimeters. This concept is similar to the classical scattering anisotropy, *g*, which describes the average cosine of the scattering angle over a single scattering length (ℓ), while scattering orientation applies to multiple light scattering based on length scales that are on the order of the transport scattering length, ℓ* ∼ ℓ/(1 − *g*).

Spatial frequency domain imaging (SFDI) is a wide-field method which measures the absorption and scattering properties of turbid media such as biological tissue.^{14, 15} It combines CCD camera-based imaging with diffuse optical methods in order to simultaneously determine the spatially varying optical properties throughout a sample. Instead of scanning a collimated source beam across the sample, optical properties at each detection location are determined simultaneously by measuring the attenuation of sinusoidal patterns of light which are projected onto the sample at varying spatial frequencies. SFDI has been used to quantitatively image stroke,^{16} brain injury,^{17} cortical spreading depression,^{18} layered structures in skin,^{19} and depth resolved fluorescent signals.^{20} It has also been used to perform tomographic imaging,^{21, 22} detect inhomogeneities,^{23} and determine optical properties over a broad range.^{24}

In this manuscript, we use SFDI to image spatially varying scattering orientation in highly scattering media. We begin by deriving a solution to an anisotropic diffusion equation in the spatial Fourier domain. The solution predicts that in turbid media, the attenuation of the projected sinusoidal patterns depends on the orientation of the projected patterns with respect to the subsurface scattering structures. When orientated along the same direction as the underlying scattering structures, the sinusoidal patterns attenuate more quickly with increasing frequencies. If the scattering structures are not parallel or orthogonal to the surface, there is a phase shift of the projected patterns. By rotating the projection patterns and measuring their attenuation and phase shift, we image the spatially varying orientation over a large field of view, and determine the relative direction of the underlying light scattering structures. We also introduce a contrast function, the scattering orientation index (SOI), that describes the degree of scattering orientation independent of other factors such as the overall amount of light scattering and absorption, both of which may also vary spatially. We validate our method using tissue simulating phantoms with fibrous structures, and with *ex vivo* samples of muscle and brain.

## 2.

## Theory

The propagation of light in turbid media is well described by the radiative transport equation (RTE) which appears as:

## Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} && {\bf \hat s} \cdot \nabla I({\bf r},{\bf \hat s}) + (\mu _a + \mu _s)I({\bf r},{\bf \hat s})\nonumber\\ && \quad -\, \mu _s \int {d^2 s'P({\bf \hat s},{\bf \hat s}')I({\bf r},{\bf \hat s}') = S({\bf r},{\bf \hat s})}. \end{eqnarray}\end{document} $$\begin{array}{ccc}& & \widehat{\mathbf{s}}\xb7\nabla I(\mathbf{r},\widehat{\mathbf{s}})+({\mu}_{a}+{\mu}_{s})I(\mathbf{r},\widehat{\mathbf{s}})\hfill \\ & & \phantom{\rule{1em}{0ex}}-\phantom{\rule{0.16em}{0ex}}{\mu}_{s}\int {d}^{2}{s}^{\prime}P(\widehat{\mathbf{s}},{\widehat{\mathbf{s}}}^{\prime})I(\mathbf{r},{\widehat{\mathbf{s}}}^{\prime})=S(\mathbf{r},\widehat{\mathbf{s}}).\hfill \end{array}$$_{a}is the absorption coefficient, μ

_{s}is the scattering coefficient, [TeX:] $S({\bf r},{\bf \hat s})$ $S(\mathbf{r},\widehat{\mathbf{s}})$ is the source, [TeX:] ${\bf \hat s}$ $\widehat{\mathbf{s}}$ is a unit vector which specifies the propagation direction, and [TeX:] $P({\bf \hat s},{\bf \hat s}^\prime)$ $P(\widehat{\mathbf{s}},{\widehat{\mathbf{s}}}^{\prime})$ is the phase function. The phase function gives the probability for a photon traveling in the direction [TeX:] ${\bf \hat s}$ $\widehat{\mathbf{s}}$ to be scattered into the direction [TeX:] ${\bf \hat s}^\prime $ ${\widehat{\mathbf{s}}}^{\prime}$ . In order to model the propagation of multiply scattered light in turbid media with oriented structures, we use the anisotropic diffusion approximation to the RTE presented by Heino et. al.

^{13}In this model, light scattering is taken to have an isotropic component, which comes from randomly orientated or spherical structures, and an anisotropic component due to aligned cylindrical structures. Accordingly, the phase function appearing in Eq. 1 is the sum of an isotropic, and an anisotropic term. It appears as:

## Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} P({\bf \hat s},{\bf \hat s'}){\rm } &=& {\rm }(1 - f){\rm }H_3 ({\bf \hat s} \cdot {\bf \hat s'},g_0){\rm }\nonumber\\ && +\, {\rm }f{\rm }\delta (\cos \theta - \cos \theta '){\rm }H_2 (\varphi - \varphi ',g_ \bot). \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill P(\widehat{\mathbf{s}},{\widehat{\mathbf{s}}}^{\prime})& =& (1-f){H}_{3}(\widehat{\mathbf{s}}\xb7{\widehat{\mathbf{s}}}^{\prime},{g}_{0})\hfill \\ & & +\phantom{\rule{0.16em}{0ex}}f\delta (\mathrm{cos}\theta -\mathrm{cos}{\theta}^{\prime}){H}_{2}(\varphi -{\varphi}^{\prime},{g}_{\perp}).\hfill \end{array}$$*H*

_{2}and

*H*

_{3}are the two- and three-dimensional Henyey-Greenstein functions,

^{25}

*g*

_{⊥}and

*g*

_{0}are the two- and three-dimensional anisotropy factors, and

*f*represents the fraction of scatterers that are anisotropic. By following the standard P1 approximation to the RTE (Ref. 26) using the phase function of Eq. 2, one arrives at an anisotropic diffusion equation. For continuous-wave light the Greens function for the anisotropic diffusion equation obeys

## Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} - \nabla \cdot {\bf D}({\bf r}) \cdot \nabla G({\bf r},{\bf r'}) + c\mu _a ({\bf r})G({\bf r},{\bf r'}) = \delta ({\bf r} - {\bf r}'), \end{equation}\end{document} $$-\nabla \xb7\mathbf{D}\left(\mathbf{r}\right)\xb7\nabla G(\mathbf{r},{\mathbf{r}}^{\prime})+c{\mu}_{a}\left(\mathbf{r}\right)G(\mathbf{r},{\mathbf{r}}^{\prime})=\delta (\mathbf{r}-{\mathbf{r}}^{\prime}),$$*c*is the speed of light and

**D**is the diffusion coefficient. In contrast to the isotropic diffusion approximation in which the diffusion coefficient is a scalar, in the anisotropic case the diffusion coefficient

**D**is a tensor. If the scattering cylinders are aligned with one of the coordinate axes,

**D**can be represented as a diagonal matrix. For arbitrary cylinder directions the matrix will not be diagonal, but a change of coordinates can bring

**D**to a diagonal form (i.e.,

**D**

_{rot}=

**RDR**

^{T}where

**R**is a rotation matrix). For cylinders aligned along the first coordinate axis, the diffusion tensor is written as

## Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} D = \left({\begin{array}{c@{\quad}c@{\quad}c} B & 0 & 0 \\ 0 & A & 0 \\ 0 & 0 & A \\ \end{array}} \right), \end{equation}\end{document} $$D=\left(\begin{array}{ccc}\hfill B& \hfill 0& \hfill 0\\ \hfill 0& \hfill A& \hfill 0\\ \hfill 0& \hfill 0& \hfill A\end{array}\right),$$## Eq. 5

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} A = \frac{{c/3}}{{\mu _a + \mu _s [1 - fg_ \bot - (1 - f)g]}}{\rm,} \end{equation}\end{document} $$A=\frac{c/3}{{\mu}_{a}+{\mu}_{s}[1-f{g}_{\perp}-(1-f)g]},$$## Eq. 6

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} B = \frac{{c/3}}{{\mu _a + \mu _s (1 - g)(1 - f)}}{\rm.} \end{equation}\end{document} $$B=\frac{c/3}{{\mu}_{a}+{\mu}_{s}(1-g)(1-f)}.$$*B*is what the diffusion coefficient would be in the absence of the cylindrical scatters (i.e., the isotropic case), whereas

*A*contains the additional scattering due to the cylinders.

Heino et. al. solved Eq. 3 numerically using finite element methods. Here we show that for simple geometries and homogeneous optical properties, Eq. 3 can be solved analytically by decomposing the Greens function into plane waves as has been done for the isotropic diffusion equation.^{21, 27, 28} In the semi-infinite and slab geometries, there is translational invariance along the medium's surface. Let the surface of the diffusive medium be in the *x*-*y* plane, and let *z* represent the depth from the surface. Accordingly we can write Greens function as

## Eq. 7

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} G({\bf r},{\bf r'}) = \int {\frac{{d^{{\rm }2} q}}{{\left({2\pi } \right)^2 }}g ({{\bf q},z,z'})} \exp [ {i{\bf q} \cdot ({\bm \rho }' - {\bm \rho })} ], \end{equation}\end{document} $$G(\mathbf{r},{\mathbf{r}}^{\prime})=\int \frac{{d}^{2}q}{{\left(2\pi \right)}^{2}}g\left(\mathbf{q},z,{z}^{\prime}\right)\mathrm{exp}\left[i\mathbf{q}\xb7({\mathit{\rho}}^{\prime}-\mathit{\rho})\right],$$*x*

**e**

_{x}+

*y*

**e**

_{y}, and

**q**is the two-dimensional wave number specifying the frequency and direction of the plane wave. In order for the Greens function of Eq. 7 to satisfy the anisotropic diffusion Eq. 3,

*g*(

**q**,

*z*,

*z*

^{′}) must satisfy the one-dimensional equation

## Eq. 8

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \left[ {\frac{{\partial ^2 }}{{\partial z^2 }} - 2i\beta \frac{\partial }{{\partial z}} - \gamma } \right]g\left({{\bf q},z,z'} \right) = \frac{{ - 1}}{{D_{zz} }}\delta \left({z - z'} \right), \end{equation}\end{document} $$\left[\frac{{\partial}^{2}}{\partial {z}^{2}}-2i\beta \frac{\partial}{\partial z}-\gamma \right]g\left(\mathbf{q},z,{z}^{\prime}\right)=\frac{-1}{{D}_{zz}}\delta \left(z-{z}^{\prime}\right),$$## Eq. 9

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \beta = \frac{{D_{xz} q_x + D_{yz} q_y }}{{D_{zz} }}, \end{equation}\end{document} $$\beta =\frac{{D}_{xz}{q}_{x}+{D}_{yz}{q}_{y}}{{D}_{zz}},$$## Eq. 10

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \gamma = \frac{{D_{xx} q_x^2 + D_{yy} q_y^2 + 2D_{xy} q_x q_y + c\mu _a }}{{D_{zz} }}. \end{equation}\end{document} $$\gamma =\frac{{D}_{xx}{q}_{x}^{2}+{D}_{yy}{q}_{y}^{2}+2{D}_{xy}{q}_{x}{q}_{y}+c{\mu}_{a}}{{D}_{zz}}.$$*g*(

**q**,

*z*,

*z*

^{′}) we must specify the boundary conditions.

^{29}Here we concentrate on the semi-infinite geometry since we measure the backscattered light in our experiments, and the thickness of the sample was large compared to other spatial scales in the problem. Transmission through a slab will be investigated in future work. Accordingly we have the boundary conditions

## Eq. 11

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} g({{\bf q},0,z'}) + \ell \frac{\partial }{{\partial z}}g({{\bf q},0,z'}) = 0, \end{equation}\end{document} $$g\left(\mathbf{q},0,{z}^{\prime}\right)+\ell \frac{\partial}{\partial z}g\left(\mathbf{q},0,{z}^{\prime}\right)=0,$$## Eq. 12

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} g({{\bf q},\infty,z'}) < \infty, \end{equation}\end{document} $$g\left(\mathbf{q},\infty ,{z}^{\prime}\right)<\infty ,$$*z = z*

^{′}, the function

*g*(

**q**,

*z*,

*z*

^{′}) is continuous, and its derivative is discontinuous.

^{30}The solution for the plane wave components of Greens function is then

## Eq. 13

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} g({\bf q},z,z') = \frac{{\exp [ {i\beta ({z - z'})}]}}{{2D_{zz} Q}}\left\{\vphantom{\frac{{1 - \ell ({i\beta + Q})}} {{1 - \ell ({i\beta - Q})}}} \exp [ { - Q{\rm }| {z - z'}|{\rm }} ] \right.\nonumber\\ \left. - {\rm }\frac{{1 - \ell ({i\beta + Q})}} {{1 - \ell ({i\beta - Q})}}\exp[ { - Q({z + z'})} ]\right\}, \end{eqnarray}\end{document} $$\begin{array}{c}\hfill g(\mathbf{q},z,{z}^{\prime})=\frac{\mathrm{exp}\left[i\beta \left(z-{z}^{\prime}\right)\right]}{2{D}_{zz}Q}\left\{\mathrm{exp}\left[-Q|z-{z}^{\prime}|\right]\right.\\ \hfill \left.-\frac{1-\ell \left(i\beta +Q\right)}{1-\ell \left(i\beta -Q\right)}\mathrm{exp}\left[-Q\left(z+{z}^{\prime}\right)\right]\right\},\end{array}$$## Eq. 14

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} Q = \sqrt {\gamma - \beta ^2 }. \end{equation}\end{document} $$Q=\sqrt{\gamma -{\beta}^{2}}.$$*z*= 0), then the amplitude and phase of

*g*(

**q**,

*z*,

*z*

^{′}) can be written in a simpler form where the amplitude is

## Eq. 15

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \left| {g\left({{\bm q},z} \right)} \right| = \frac{\ell }{{D_{zz} }}\exp \left({ - Q{\rm }z} \right){\rm } [ {\left({1 + \ell Q} \right)^2 + \left({\ell \beta } \right)^2 } ]^{ - 1/2}, \end{equation}\end{document} $$\left|g\left(\mathit{q},z\right)\right|=\frac{\ell}{{D}_{zz}}\mathrm{exp}\left(-Qz\right){\left[{\left(1+\ell Q\right)}^{2}+{\left(\ell \beta \right)}^{2}\right]}^{-1/2},$$## Eq. 16

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\rm Arg}(g) = - \beta {\rm }z + \tan ^{ - 1} \left({\frac{{\ell \beta }}{{1 + \ell Q}}} \right). \end{equation}\end{document} $$\mathrm{Arg}\left(g\right)=-\beta z+{\mathrm{tan}}^{-1}\left(\frac{\ell \beta}{1+\ell Q}\right).$$*Q*, the extrapolation length ℓ, and the depth

*z*. As the sinusoidal pattern propagates into the scattering medium it is attenuated exponentially with decay constant

*Q*, and shifted in space linearly at a rate determined by β. The functions

*Q*and β in turn depend on the scattering and absorption coefficients, the two- and three-dimensional anisotropy factors, the fraction of anisotropic scatterers, the frequency of the projected pattern, and the orientation of the anisotropic scatterers in relation to both the surface and the projected patterns. Let θ be the angle between the surface and the anisotropy direction, and let φ be the angle in the

*x*-

*y*plane between the wave number and the anisotropy direction (see Fig. 1). Without loss of generality, we can assume the wave number points in the

*x*direction. Then we have

## Eq. 17

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \beta = \frac{{(B - A){\rm }\cos \varphi {\rm }\cos \theta {\rm }\sin \theta }}{{A\cos ^2 \theta + B\sin ^2 \theta }}, \end{equation}\end{document} $$\beta =\frac{(B-A)\mathrm{cos}\varphi \mathrm{cos}\theta \mathrm{sin}\theta}{A{\mathrm{cos}}^{2}\theta +B{\mathrm{sin}}^{2}\theta},$$*A*and

*B*were defined in Eqs. 5, 6. Notice that if the scatterers are either parallel or perpendicular to the surface of the medium, there will be no phase shift. The perpendicular and parallel cases also yield simple solutions for

*Q*. For the perpendicular case we have

## Eq. 18

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} Q_ \bot = \sqrt {\frac{{Aq^2 + c\mu _a }}{B}}, \end{equation}\end{document} $${Q}_{\perp}=\sqrt{\frac{A{q}^{2}+c{\mu}_{a}}{B}},$$## Eq. 19

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} Q_{||} = \sqrt {\frac{{(A\sin ^2 \varphi + B\cos ^2 \varphi)q^2 + c\mu _a }}{A}}. \end{equation}\end{document} $${Q}_{\left|\right|}=\sqrt{\frac{(A{\mathrm{sin}}^{2}\varphi +B{\mathrm{cos}}^{2}\varphi ){q}^{2}+c{\mu}_{a}}{A}}.$$*A*⩽

*B*,

*Q*

_{||}has its maximum when φ = 0. According to Eq. 15, this means the measured amplitude will be a minimum when the wave number of the projected sine wave is aligned with the scattering structures of the medium. This is the central theoretical result of this manuscript. It means that when β is small, one can determine the orientation of the underlying structure of the medium by simply looking for the angle at which the projected pattern blurs the most (i.e., its amplitude is smallest). In Fig. 2, the amplitude ∣

*g*∣ is plotted as a function of φ for five different projection frequencies using optical properties: μ

_{a}= 0.1 mm

^{−1}, μ

_{s}= 10 mm

^{−1},

*g*

_{0}= 0.9,

*g*

_{⊥}= 0.9,

*f*= 0.5,

*z*= 0, and θ = 0. Notice that as the frequency increases, the overall amplitude decreases, but the change in amplitude as a function of projection angle increases. Thus higher frequencies are more sensitive to scattering orientation, but also give less overall signal. Since the overall amplitude decreases exponentially with depth, higher frequencies are most beneficial when interrogating superficial volumes.

## 3.

## Methods

## 3.1.

### Data Acquisition

We used a custom SFDI device developed in collaboration with Modulated Imaging, Inc. (Irvine, California). Light from a light emitting diode centered at 850 nm was coupled via an integrating rod to a digital micromirror device (DMD). The DMD projected an 8 bit sinusoidal intensity pattern onto the sample with a spatial frequency of 0.4 mm^{−1}. The light remitted from the sample was detected with a 12 bit CCD array. To make a complete measurement, the sinusoidal pattern was projected serially at three different phases separated by 2π/3 radians, and a separate image was acquired by the CCD for each phase. Then the orientation of the sinusoidal pattern was rotated by 5 deg, and the remitted light was imaged again for each of the three projection phases. This process was repeated until the projection pattern was rotated a full 180 deg. Thus a complete measurement consisted of 36 discrete angles, each with 3 phases, for a total of 108 CCD exposures.

## 3.2.

### Data Analysis

For each projection angle, the amplitude of the remitted light at every detector location was calculated according to

## Eq. 20

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} A(x,y) &=& \frac{{2^{1/2} }}{3}\{ {\rm }\left[ {I_1 (x,y) - I_2 (x,y)} \right]^2 + \left[ {I_1 (x,y) - I_3 (x,y)} \right]^2\nonumber\\ && + \left[ {I_2 (x,y) - I_3 (x,y)} \right]^2 \}^{1/2}. \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill A(x,y)& =& \frac{{2}^{1/2}}{3}\{{\left[{I}_{1}(x,y)-{I}_{2}(x,y)\right]}^{2}+{\left[{I}_{1}(x,y)-{I}_{3}(x,y)\right]}^{2}\hfill \\ & & {+{\left[{I}_{2}(x,y)-{I}_{3}(x,y)\right]}^{2}\}}^{1/2}.\hfill \end{array}$$We have discussed this demodulation and the measurement of the amplitude previously.^{14} The phase shift of the projected wave can also be calculated,^{31} and has been used for tomographic imaging.^{21} However, in this preliminary study we only looked at the amplitude of the remitted light. The tissue simulating phantoms measured in this manuscript had structures primarily orientated parallel to the surface, and accordingly little phase shift was expected [see Eqs.
16, 17]. In contrast, our *ex vivo* samples may have had structures which were not parallel to the surface. Such structures would have been difficult to detect with an amplitude only measurement because they would not have caused a large change in the measured amplitude as the projection angle was rotated. After demodulation, the amplitude was low pass filtered in space using a 4 × 4 pixel median filter, and filtered in angle using a Gaussian filter with σ = 20 deg. In order to compensate for systematic errors, the amplitude for each projection angle was calibrated by dividing it by the amplitude for that angle measured using an isotropic tissue simulating phantom.

We then determined the orientation of the light scattering structures, as well as a SOI. In order to image the orientation, for each pixel we found the projection angle for which the measured amplitude was a minimum, and used these values to create an image of the angular orientation of the anisotropic light scatterers. As a quantitative measure we defined a contrast function, which we referred to as the SOI. It varied from 0 to 1 as:

## Eq. 21

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} {\rm SOI} = \mathop {\rm max}\limits_\varphi \left\{ {\frac{{\left| {g(\varphi)} \right| - \left| {g(\varphi + \pi /2)} \right|}}{{\left| {g(\varphi)} \right| + \left| {g(\varphi + \pi /2)} \right|}}} \right\}, \end{equation}\end{document} $$\mathrm{SOI}=\underset{\varphi}{\mathrm{max}}\left\{\frac{\left|g\left(\varphi \right)\right|-\left|g(\varphi +\pi /2)\right|}{\left|g\left(\varphi \right)\right|+\left|g(\varphi +\pi /2)\right|}\right\},$$The SOI is primarily sensitive to the fraction of scatters *f* which are anisotropic, and to differences between the two- and three-dimensional anisotropy factors *g*
_{0} and *g*
_{⊥}. To demonstrate this point we calculate the SOI for baseline optical properties of μ_{a} = 0.1 mm^{−1}, μ_{s} = 10 mm^{−1}, *g*
_{0} = 0.5, *g*
_{⊥} = 0 .5, *f* = 0.5, *z* = 0, and θ = 0. We then individually vary the optical properties from 10% to 200% of baseline values. Results are plotted in Fig. 3. Changes in the values of μ_{a} and μ_{s} have almost no effect on the SOI; whereas changes in the value of *f* have the largest effect. Note, the SOI goes to zero when *f* goes to zero. When both *g*
_{0} and *g*
_{⊥} are changed together, there is little effect on the SOI. However, when either *g*
_{0} and *g*
_{⊥} is changed alone, then the relative amount of isotropic and anisotropic scattering is changed, and the SOI is affected. For example, in Fig. (3), *g*
_{0} is increased while *g*
_{⊥} is held fixed. This results in a decrease in isotropic scattering while the anisotropic scattering does not change. As a result, the SOI increases.

## 3.3.

### Sample Preparation

In order to validate our imaging method, we scanned both tissue simulating phantoms and *ex vivo* tissue samples using the procedure outlined above. Tissue simulating phantoms contained both isotropic regions and anisotropic regions in which the orientation of the anisotropy was varied. The phantoms utilized a large silicon base containing titanium dioxide for scattering and India ink for absorption. We have used this type of phantom extensively.^{32} Sections of a pleated air filter were soaked in water and then pressed against the silicon base. The fibers composing the filters clearly had a preferred orientation and filter sections were rotated so that the phantom contained regions with differing fiber orientations [see Fig. 4a]. In order to evaluate our ability to image subsurface scattering orientation, thin sheets of highly scattering (μ^{′}
_{s} = 1.5 mm^{−1}) silicon phantoms were pressed over the filter sections. The sheets had thicknesses ranging from 200 to 600 μm. *Ex vivo* samples consisted of chicken breast (i.e., muscle) and a rat brain. The chicken breast was cut into slices (∼5 mm thick) and each slice was positioned on top of the silicon base with a different orientation. The rat brain was taken from an 8 week old male Sprague Dawley rat immediately follow euthanasia using pentobarbital. The brain was carefully removed and sectioned into two slices along the coronal plane with a thickness of ∼5 mm.

## 4.

## Results

Figure 4 demonstrates the ability of SFDI to image the orientation of scattering structures at the surface of a turbid medium. Nine fiber samples from a pleated air filter were fixed at various orientations on the surface of an isotropic tissue simulating phantom and imaged as described above. A black and white photograph [Fig. 4a] shows the position of the fibrous samples. The inset is a confocal reflectance image of one of the fibrous samples in which the orientation of the microscopic fibers is visible. The SOI image [Fig. 4b] clearly indicates that the scattering from the fiber samples is more anisotropic than that of the surrounding medium. We then used an arbitrary threshold to aid in the display of the SOI. In Fig. 4c, only values of the SOI above the threshold value indicated on the color bar appear, and these values are superimposed over the black and white photograph. Figure 4d is a color coded image of the orientation of the microscopic fibers composed of the fibrous samples. The color wheel to the left gives the correspondence between each color and the direction it denotes. Starting from the bottom left and moving upwards column by column, the orientation of the fibrous samples is rotated clockwise. This is clearly shown in the image by the overall color of each fibrous sample. The areas between the fiber samples also are assigned values. This occurs because the processing code simply picks the projection angle for which each pixel has its minimum value. As expected, in areas where the scattering is isotropic, the assigned values are random. This motivates using the threshold from the SOI image to choose which pixels to display in the orientation image. In Fig. 4e, only pixels in which the contrast was above its threshold value are shown superimposed over the black and white image. For these pixels, the direction of the underlying structures was clearly measured. In Fig. 4f we plot the amplitude as a function of the projection angle for three individual CCD pixels taken from the centers of the fibrous samples in the bottom row. The plots in Fig. 4g are for the same three CCD pixels after the spatial and angular filtering described in Sec. 3.2. Each pixel is assigned an angle corresponding to the minimum of its curve.

A potential advantage of using SFDI to image scattering orientation, as opposed to methods involving polarization, is that diffuse optical methods such as SFDI use multiply scattered light. Thus, they are able to probe underneath the tissue surface at depths where light loses its polarization. In order to test the ability of SFDI to image subsurface scattering orientation, we placed thin sheets of highly scattering material (μ^{′}
_{s} = 1.5 mm^{−1}) over the fibrous samples. In Figs. 5a, 5b, we show images of SOI and fiber orientation, as the thickness of the overlying medium is increased from 0 to 600 μm. As the thickness increases, the SOI becomes smaller, and the image quality deteriorates. Nonetheless, even at a depth of 600 μm, the method assigns the correct orientation to much of the image. We quantified the degradation of image quality with depth by picking a 50 × 50 pixel (11.5 × 11.5 mm) region of interest at the center of each fibrous sample and calculating the standard deviation of angular values measured for each depth within the region of interest. We obtained standard deviations of 3.8, 4.7, 12, and 33 deg for depths of 0, 200, 400, and 600 μm, respectively. We also attempted to image through a silicon sheet with a thickness of 1 mm. At this depth no scattering orientation was detected. While it may be possible to penetrate more deeply using lower spatial frequencies and/or longer wavelengths, it is clear from Fig. 2 that lower frequencies come with the price of decreased sensitivity to the scattering orientation.

Finally we evaluated the ability of SFDI to image scattering orientation in *ex vivo* samples of muscle (in this case chicken breast) and brain. As shown in Fig. 6a, four slices of chicken breast were positioned over the tissue simulating phantom such that the orientation of their muscle fibers was varied in 45 deg increments. The chicken breast clearly exhibited scattering anisotropy, and the imaged orientation matched our expectations, albeit with some biological variation. Figure 6b displays the results from imaging coronal slices of the brain of a rat. White matter composing the corpus callosum is visible in the black and white image. We expect that the oriented myelinated axons predominately found in white matter, will result in anisotropic diffusion of diffuse light. The myelinated axons are already known to cause the anisotropy diffusion of water, and form the basis for diffusion tensor magnetic resonance imaging (DT-MRI).^{33} In the contrast image, areas of high contrast tend to correspond to white matter. Although we are unable to validate that the axon directions are correctly determined in this preliminary experiment, the orientation of white matter structures found with SFDI could be compared with DT-MRI in future work. These *ex vivo* experiments demonstrate the potential for using SFDI as a tool to image scattering orientation in biological tissues.

## 5.

## Discussion

In this work we use a model of light propagation based on the diffusion approximation to the RTE. This approximation is known to be most accurate when 1. the reduced scattering coefficient is much larger than the absorption coefficient, 2. the distance between sources and detectors is much greater than one over the reduced scattering coefficient, and 3. the volume of interest is far from boundaries compared to one over the reduced scattering coefficient. In our experiment, conditions 2 and 3 are not met since the spatially distributed source is co-localized with all of the detectors, and we are only sensitive to anisotropy within ∼1 mm from the surface. Nonetheless, we (and others) have shown that the diffuse approximation can be used to accurately predict the diffuse reflectance from isotropic media.^{34} However, we have also shown a slight increase in accuracy in our SFDI measurements of isotropic media when modeling light transport with Monte Carlo solutions to the RTE as opposed to using the diffusion approximation.^{14} We do expect the same qualitative results for SFDI using the anisotropic diffusion approximation as we would obtain by solving the full RTE. However, the extent to which the diffusion approximation can be used to model light propagation in anisotropic turbid media has been controversial,^{9, 35} and will be the subject of further SFDI research.

In this manuscript, we did not fit for the exact optical properties of the medium. Instead we determined the orientation and degree of anisotropy by finding the angle corresponding to the minimum measured amplitude, and defining a contrast function. Indeed, the anisotropic diffusion model contains many free parameters, and we have yet to determine to what extent we can solve for them using SFDI measurements.

Finally, it is important to point out that although we are using SFDI to image heterogeneous media, our solution to the anisotropic diffusion equation assumed the media was homogeneous. We analyzed the data from each CCD pixel independently assuming that the optical properties of the sample were sufficiently homogeneous such that the optical properties of neighboring locations did not affect the properties assigned to a given location. This was sufficient for demonstrating our sensitivity to scattering orientation, and for determining the orientation of the light scattering structures in two-dimensions. However, this approach has limitations. For example, the results of Fig. 5 clearly show that the measured contrast is determined not only by the properties of the fiber sample, but also by its depth. Thus, more sophisticated approaches involving layered models and/or tomographic imaging may be useful to facilitate quantitative measurements of scattering orientation in heterogeneous tissues. We have implemented both layered models^{19} and tomographic imaging^{21} for SFDI in isotropic media, and we expect to be able to extend these approaches to anisotropic tissues.

## 6.

## Conclusion

We used SFDI to image the orientation of light scattering structures in highly scattering media. This was accomplished by rotating the orientation of sinusoidal patterns of light intensity projected on the sample, and measuring the attenuation in amplitude of the patterns as a function of the projection angle. We derived an analytic solution in the Fourier domain to the anisotropic diffusion equation, and defined a contrast function, the scattering orientation index, which gives a measure of the degree to which multiple light scattering has a preferred direction. We validated our approach by imaging tissue simulating phantoms as well as *ex vivo* samples of muscle and brain. We were able to distinguish between regions with isotropic and directional scattering, image the orientation of scattering structures embedded at a depth of 600 μm in highly scattering isotropic media (μ^{′}
_{s} = 1.5 mm^{−1}), and determine the scattering orientation within a standard deviation of 3.8 deg using individual CCD pixels. This technique may aid in the diagnosis of tissues such as skin, muscle, bone, and brain; in which the loss of orientated structure may be associated with disease.

## Acknowledgments

Support for this work was provided by the NIH NCRR Laser Microbeam and Medical Program (LAMMP, P41-RR01192), the Military Photomedicine Program (AFOSR Grant No. FA9550-08-1-0384), and the Beckman Foundation. Soren D. Konecky was supported by a fellowship from the Hewitt Foundation for Medical Research.