Translator Disclaimer
1 December 2011 Imaging scattering orientation with spatial frequency domain imaging
Author Affiliations +
Optical imaging techniques based on multiple light scattering generally have poor sensitivity to the orientation and direction of microscopic light scattering structures. In order to address this limitation, we introduce a spatial frequency domain method for imaging contrast from oriented scattering structures by measuring the angular-dependence of structured light reflectance. The measurement is made by projecting sinusoidal patterns of light intensity on a sample, and measuring the degree to which the patterns are blurred as a function of the projection angle. We derive a spatial Fourier domain solution to an anisotropic diffusion model. This solution predicts the effects of bulk scattering orientation on the amplitude and phase of the projected patterns. We introduce a new contrast function based on a scattering orientation index (SOI) which is sensitive to the degree to which light scattering is directionally dependent. We validate the technique using tissue simulating phantoms, and ex vivo samples of muscle and brain. Our results show that SOI is independent of the overall amount of bulk light scattering and absorption, and that isotropic versus oriented scattering structures can be clearly distinguished. We determine the orientation of subsurface microscopic scattering structures located up to 600 μm beneath highly scattering (μ′s = 1.5 mm−1) material.



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 walk11, 12 and finite element implementations13 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.



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} ŝ·I(r,ŝ)+(μa+μs)I(r,ŝ)μsd2sP(ŝ,ŝ)I(r,ŝ)=S(r,ŝ).
Here [TeX:] $I({\bf r},{\bf \hat s})$ I(r,ŝ) is the specific intensity, μa is the absorption coefficient, μs is the scattering coefficient, [TeX:] $S({\bf r},{\bf \hat s})$ S(r,ŝ) is the source, [TeX:] ${\bf \hat s}$ ŝ is a unit vector which specifies the propagation direction, and [TeX:] $P({\bf \hat s},{\bf \hat s}^\prime)$ P(ŝ,ŝ) is the phase function. The phase function gives the probability for a photon traveling in the direction [TeX:] ${\bf \hat s}$ ŝ to be scattered into the direction [TeX:] ${\bf \hat 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} P(ŝ,ŝ)=(1f)H3(ŝ·ŝ,g0)+fδ(cosθcosθ)H2(ϕϕ,g).
Here, 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} ·D(r)·G(r,r)+cμa(r)G(r,r)=δ(rr),
where 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=B000A000A,

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=c/3μa+μs[1fg(1f)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=c/3μa+μs(1g)(1f).
Note that 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(r,r)=d2q2π2g(q,z,z)exp[iq·(ρρ)],
where [TeX:] ${\bm\rho}$ ρ = 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} 2z22iβzγgq,z,z=1Dzzδzz,

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} β=Dxzqx+DyzqyDzz,

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} γ=Dxxqx2+Dyyqy2+2Dxyqxqy+cμaDzz.
In order to solve for 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(q,0,z)+zg(q,0,z)=0,

Eq. 12

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} g({{\bf q},\infty,z'}) < \infty, \end{equation}\end{document} g(q,,z)<,
where ℓ is the extrapolation length. In addition we require that the at the point 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} g(q,z,z)=exp[iβ(zz)]2DzzQexp[Q|zz|]1(iβ+Q)1(iβQ)exp[Q(z+z)],

Eq. 14

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} Q = \sqrt {\gamma - \beta ^2 }. \end{equation}\end{document} Q=γβ2.
In our SFDI experiments, the sample is illuminated with sinusoidal patterns of light intensity. The amplitude and phase shift of the projected patterns inside and remitted from the sample are governed by Eq. 13. If the source is located at the surface (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} gq,z=DzzexpQz[1+Q2+β2]1/2,
and the phase shift is

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} Arg(g)=βz+tan1β1+Q.
The amplitude and phase shift depend on β, 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} β=(BA)cosϕcosθsinθAcos2θ+Bsin2θ,
where 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=Aq2+cμaB,
which is independent of the projection angle φ. However, for the parallel case we have

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||=(Asin2ϕ+Bcos2ϕ)q2+cμaA.
Since AB, 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.

Fig. 1

Schematic of the experimental geometry. Gray rectangles represent cylindrical light scattering structures in a turbid medium. Left: top view of the medium where φ is the angle in the x-y plane between the wave number q of the sinusoidal illumination pattern, and the projection of the cylindrical structures onto the x-y plane. Right: side view of the medium where θ is the angle between the surface and the cylinder direction.


Fig. 2

Plot of the calculated amplitude of remitted light as a function of the angle φ between the wave number of the projected sinusoidal pattern, and the direction of the cylindrical scattering structures. The amplitudes due to 5 projection frequencies ranging from 0.1 to 0.5 mm−1 are shown. The minimum occurs when the wave number is aligned with the scattering structures. Higher frequencies are more sensitive to anisotropy, but have smaller average amplitude.





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.


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} A(x,y)=21/23{I1(x,y)I2(x,y)2+I1(x,y)I3(x,y)2+I2(x,y)I3(x,y)2}1/2.

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} SOI=maxϕg(ϕ)g(ϕ+π/2)g(ϕ)+g(ϕ+π/2),
where φ is the projection angle. This function was calculated for every pixel on the detector array in order to create an image.

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.

Fig. 3

Plots showing the effect on the SOI when various optical properties are changed. The SOI is most sensitive to changes in the fraction f of light scattering structures which are anisotropic. Changing the values of μa and μs has almost no effect on the SOI. Likewise, when the anisotropy factors g 0 and g are both varied in the same manner, there is little effect on the SOI. However, when the ratio between g 0 and g changes, the relative amounts of isotropic and anisotropic scattering changes, leading to a change in the SOI. In the graph above, increasing g 0 while g is fixed leads to an increase in the SOI due to the decrease in isotropic scattering.



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.

Fig. 4

Results from an imaging experiment in which nine fibrous samples were attached to the surface of an isotropic light scattering phantom. The fibrous samples were rotated so that the orientations of the fibers for each sample were different. (a) A black and white photograph of the entire sample. Inset is a confocal reflectance image of one of the fibrous samples showing the alignment of the microscopic fibers. (b) Image of the SOI indicating that the fibrous samples have more anisotropic scattering than the underlying phantom. (c) Thresholded image of the SOI superimposed over the black and white photograph. (d) Color coded image of the orientation of the light scattering structures in which each color represents a different angle as denoted by the color wheel to the left. (e) Thresholded image of the structure orientation superimposed over the black and white photograph. (f) Plots of the measured amplitude |g(φ)| as a function of angle for three individual CCD pixels corresponding to the center of three fibrous samples in the bottom row. (g) Plots of the same CCD pixels in (f) but after the spatial and angular filtering described in the text.




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.

Fig. 5

Images of the (a) SOI and (b) structure orientation from an experiment in which a layer of highly scattering material was put on top of the fibrous samples. The top layer was varied in thickness from 0 to 600 μm.


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.

Fig. 6

Images of ex vivo samples of (a) chicken breast and (b) rat brain. From left to right: a black and white photograph, a SOI image, and a structure orientation image are shown for each sample. SOI and orientation images are superimposed over the black and white photograph after a threshold was applied.




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 models19 and tomographic imaging21 for SFDI in isotropic media, and we expect to be able to extend these approaches to anisotropic tissues.



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.


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.



A. Kienle, C. Wetzel, A. Bassi, D. Comelli, P. Taroni, and A. Pifferi, “Determination of the optical properties of anisotropic biological media using an isotropic diffusion model,” J. Biomed. Opt., 12 (1), 014026 (2007). Google Scholar


G. Marquez, L. V. Wang, S. P. Lin, J. A. Schwartz, and S. L. Thomsen, “Anisotropy in the absorption and scattering spectra of chicken breast tissue,” Appl. Op, 37 (4), 798 –804 (1998). Google Scholar


S. Nickell, M. Hermann, M. Essenpreis, T. J. Farrell, U. Kramer, and M. S. Patterson, “Anisotropy of light propagation in human skin,” Phys. Med. Biol., 45 (10), 2873 –2886 (2000). Google Scholar


A. Sviridov, V. Chernomordik, M. Hassan, A. Russo, A. Eidsath, P. Smith, and A. H. Gandjbakhche, “Intensity profiles of linearly polarized light backscattered from skin and tissue-like phantoms,” J. Biomed. Opt., 10 (1), 014012 (2005). Google Scholar


T. Binzoni, C. Courvoisier R. Giust, G. Tribillon, T. Gharbi, J. C. Hebden, T. S. Leung, J. Roux, and D. T. Delpy, “Anisotropic photon migration in human skeletal muscle,” Phys. Med. Biol., 51 (5), N79 –N90 (2006). Google Scholar


A. Kienle and R. Hibst, “Light guiding in biological tissue due to scattering,” Phys. Rev. Lett., 97 (1), 018104 (2006). Google Scholar


J. C. Hebden, J. J. Guerrero, V. Chernomordik, and A. H. Gandjbakhche, “Experimental evaluation of an anisotropic scattering model of a slab geometry,” Opt. Lett., 29 (21), 2518 –2520 (2004). Google Scholar


P. M. Johnson, S. Faez, and A. Lagendijk, “Full characterization of anisotropic diffuse light,” Opt. Express, 16 (10), 7435 –7446 (2008). Google Scholar


P. M. Johnson and A. Lagendijk, J. Biomed. Opt., 14 (5), 054036 (2009). Google Scholar


A. Kienle, F. K. Forster, and R. Hibst, “Anisotropy of light propagation in biological tissue,” Opt. Lett., 29 (22), 2617 –2619 (2004). Google Scholar


L. Dagdug, G. H. Weiss, and A. H. Gandjbakhche, “Effects of anisotropic optical properties on photon migration in structured tissues,” Phys. Med. Biol., 48 (10), 1361 –1370 (2003). Google Scholar


O. K. Dudko, G. H. Weiss, V. Chernomordik, and A. H. Gandjbakhche, “Photon migration in turbid media with anisotropic optical properties,” Phys. Med. Biol., 49 (17), 3979 –3989 (2004). Google Scholar


J. Heino, S. Arridge, J. Sikora, and E. Somersalo, “Anisotropic effects in highly scattering media,” Phys. Rev. E, 68 (3 Pt 1), 031908 (2003). Google Scholar


D. J. Cuccia, F. Bevilacqua, A. J. Durkin, F. R. Ayers, and B. J. Tromberg, “Quantitation and mapping of tissue optical properties using modulated imaging,” J. Biomed. Opt., 14 (2), 024012 (2009). Google Scholar


D. J. Cuccia, F. Bevilacqua, A. J. Durkin, and B. J. Tromberg, “Modulated imaging: quantitative analysis and tomography of turbid media in the spatial-frequency domain,” Opt. Lett., 30 (11), 1354 –1356 (2005). Google Scholar


D. Abookasis, “Imaging cortical absorption, scattering, and hemodynamic response during ischemic stroke using spatially modulated near-infrared illumination,” J. Biomed. Opt., 14 (2), 024033 (2009). Google Scholar


J. R. Weber, D. J. Cuccia, W. R. Johnson, G. H. Bearman, A. J. Durkin, M. Hsu, A. Lin, D. K. Binder, D. Wilson, and B. J. Tromberg, “Multispectral imaging of tissue absorption and scattering using spatial frequency domain imaging and a computed-tomography imaging spectrometer,” J. Biomed. Opt., 16 (1), 011015 (2011). Google Scholar


D. J. Cuccia, D. Abookasis, R. D. Frostig, and B. J. Tromberg, “Quantitative in vivo imaging of tissue absorption, scattering, and hemoglobin concentration in rat cortex using spatially-modulated structured light,” In Vivo Optical Imaging of Brain Function, 2nd Ed.CRC, Boca Raton, FL (2009). Google Scholar


J. R. Weber, D. J. Cuccia, A. J. Durkin, and B. J. Tromberg, “Noncontact imaging of absorption and scattering in layered tissue using spatially modulated structured light,” J. Appl. Phys., 105 102028 (2009). Google Scholar


A. Mazhar, D. J. Cuccia, S. Gioux, A. J. Durkin, J. V. Frangioni, and B. J. Tromberg, “Structured illumination enhances resolution and contrast in thick tissue fluorescence imaging,” J. Biomed. Opt., 15 (1), 010506 (2010). Google Scholar


S. D. Konecky, A. Mazhar, D. Cuccia, A. J. Durkin, J. C. Schotland, and B. J. Tromberg, “Quantitative optical tomography of sub-surface heterogeneities using spatially modulated structured light,” Opt. Express, 17 (17), 14780 –14790 (2009). Google Scholar


S. Belanger, M. Abran, X. Intes, C. Casanova, and F. Lesage, “Real-time diffuse optical tomography based on structured illumination,” J. Biomed. Opt., 15 (1), 016006 (2010). Google Scholar


A. Bassi, C. D'Andrea, G. Valentini, R. Cubeddu, and S. Arridge, “Detection of inhomogeneities in diffusive media using spatially modulated light,” Opt. Lett., 34 (14), 2156 –2158 (2009). Google Scholar


T. A. Erickson, A. Mazhar, D. Cuccia, A. J. Durkin, and J. W. Tunnell, “Lookup-table method for imaging optical properties with structured illumination beyond the diffusion theory regime,” J. Biomed. Opt., 15 (3), 036013 (2010). Google Scholar


L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the Galaxy,” Astrophys. J., 93 70 –83 (1941). Google Scholar


L. V. Wang and H. Wu, Biomedical Optics: Principles and Imaging, John Wiley & Sons, Inc., Hoboken, NJ (2007). Google Scholar


V. Lukic, V. A. Markel, and J. C. Schotland, “Optical tomography with structured illumination,” Opt. Lett., 34 (7), 983 –985 (2009). Google Scholar


V. A. Markel and J. C. Schotland, “Symmetries, inversion formulas, and image reconstruction for optical tomography,” Phys. Rev. E, 70 (5 Pt 2), 056616 (2004). Google Scholar


V. A. Markel and J. C. Schotland, “Inverse problem in optical diffusion tomography. II. Role of boundary conditions,” J. Opt. Soc. Am. A, 19 (3), 558 –566 (2002). Google Scholar


G. B. Arfken and H. J. Weber, Mathematical Methods for Physicists, 6th Ed.Academic Press, Burlington, MA (2005). Google Scholar


A. Bassi, “Spatial shift of spatially modulated light projected on turbid media,” J. Opt. Soc. Am. A, 25 (11), 2833 –2839 (2008). Google Scholar


M. E. Moseley, Y. Cohen, J. Kucharczyk, J. Mintorovitch, H. S. Asgari, M. F. Wendland, J. Tsuruda, and D. Norman, “Diffusion-weighted MR imaging of anisotropic water diffusion in cat central nervous system,” Radiology, 176 (2), 439 –445 (1990). Google Scholar


L. O. Svaasand, T. Spott, J. B. Fishkin, T. Pham, B. J. Tromberg, and M. W. Berns, “Reflectance measurements of layered media with diffuse photon-density waves: a potential tool for evaluating deep burns and subcutaneous lesions,” Phys. Med. Biol., 44 (3), 801 –813 (1999). Google Scholar


A. Kienle, “Anisotropic light diffusion: an oxymoron,” Phys. Rev. Lett., 98 (21), 218104 (2007). Google Scholar
© 2011 Society of Photo-Optical Instrumentation Engineers (SPIE) 1083-3668/2011/16(12)/126001/8/$25.00
Soren D. Konecky, Tyler B. Rice, Anthony J. Durkin, and Bruce J. Tromberg "Imaging scattering orientation with spatial frequency domain imaging," Journal of Biomedical Optics 16(12), 126001 (1 December 2011).
Published: 1 December 2011

Back to Top