_{s}

^{′}/μ

_{a}) ratios (8 to 500), reporting an overall accuracy of approximately 6% and 3% in absorption and reduced scattering parameters, respectively. Sampling of only two spatial frequencies, achieved with only three camera images, is found to be sufficient for accurate determination of the optical properties. We then perform MI measurements in an in vivo tissue system, demonstrating spatial mapping of the absorption and scattering optical contrast in a human forearm and dynamic measurements of a forearm during venous occlusion. Last, metrics of spatial resolution are assessed through both simulations and measurements of spatially heterogeneous phantoms.

## 1.

## Introduction

Light transport in tissues is a complex process due to multiple scattering and absorption. Thus, at the core of every optical technique for quantitative tissue characterization is the ability to separate optical absorption from optical scattering effects by the detection of a remitted or transmitted light field. This remission (or transmission) is a function of time and space, yielding two general classes of quantitative techniques: time-resolved and spatially resolved measurements, respectively (see Fig. 1
). Time-resolved measurements are further broken down into time-domain and frequency-domain techniques: the first measuring the temporal point-spread function
$\left(t\text{-}\mathrm{PSF}\right)$
, or spreading of a propagating pulse in time,^{1, 2} and the latter measuring the temporal modulation transfer function
$\left(t\text{-}\mathrm{MTF}\right)$
, or the attenuation and phase delay of a periodically varying photon density wave.^{3, 4, 5} The time domain and frequency domain share an exact Fourier transform equivalency, although each has its trade-offs when considering real-life hardware and model-fitting constraints.

In diffuse optics, spatially resolved measurements have been generally limited to the real spatial domain. Here, the spatial point-spread function
$\left(s\text{-}\mathrm{PSF}\right)$
is typically characterized by “multidistance” measurements,^{6, 7} tracking the spatial dependence of a reflected or transmitted light field generated from a point-like illumination. The Fourier transform equivalent to the real spatial domain is the spatial frequency domain (SFD). While recent work has shown the use of spatially structured illumination techniques for manipulating diffractive optical systems,^{8} little has been reported for its use in characterization of *diffusive* systems.^{9, 10, 11}

In this paper, we describe a new imaging method, modulated imaging (MI), for quantitation and wide-field mapping of turbid media in the SFD. The spatial modulation transfer function
$\left(s\text{-}\mathrm{MTF}\right)$
of a turbid medium encodes both depth and optical property information, enabling both quantitation and tomographic imaging of the spatially varying medium optical properties.^{10} In this work, we present a detailed exposition and validation of the ability of MI to quantitatively recover homogeneous tissue optical properties. We present two homogeneous forward models of diffuse reflectance in the spatial frequency domain—the first, an analytic diffusion-based approach, and the second, a transport-based approach using Monte Carlo simulations. Next, we present reflectance measurements of tissue-simulating liquid phantoms exhibiting a wide range of absorption and scattering values. The optical properties of these samples are recovered by analysis with our analytic diffusion model using two inversion methods—the first, a least-squares multifrequency fitting algorithm, and the second, a rapid two-frequency lookup table approach. We then apply the technique to an *in vivo* tissue system, producing 2-D spatial maps of the absorption and reduced scattering contrast of a human forearm. Dynamic measurements are also acquired, demonstrating changes in forearm optical properties during venous occlusion. Last, we investigate metrics of spatial resolution and optical property contrast through both simulations and measurements of spatially heterogeneous phantoms.

## 2.

## MI Instrumentation

The MI apparatus is shown in Fig. 2 . Grayscale illumination patterns are generated using a light source in combination with a spatial light modulator (SLM). In this study, we used a simple digital projector (NEC HT1000), based on a digital micromirror-based digital light processing (DLP) light engine (Texas Instruments) and an ultra high performance (UHP) mercury lamp. The projector’s color filter wheel was removed, producing a broadband “white light” illumination of the sample, allowing us to use interference filters for detection of a narrow wavelength band (Andover Corporation, $\lambda =660\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ , $\Delta \lambda =10\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ FWHM). To create the illumination patterns, $8\text{-}\text{bit}$ grayscale bitmap images are generated using MATLAB (Mathworks, Inc.). They are then placed in a PowerPoint (Microsoft, Inc.) presentation file and automatically sequenced using the Microsoft Office ActiveX controls through an external LabVIEW (National Instruments, Inc.) program. The diffusely reflected light is captured by a $16\text{-}\text{bit}$ frame-transfer CCD camera (Roper Cascade 512F) capable of imaging up to 30 frames per second at full $512\times 512$ resolution. Specular reflection is avoided by illuminating at a small angle to normal incidence. Additionally, crossed linear polarizers can be added to further select the diffuse reflectance, useful for rough surfaces (such as skin), where specular light can be reflected at all angles.

The modularity of this system makes it very flexible. First, the field of view is limited only by the magnification of the illumination and collection optics (with fundamental resolution limits set by the physics of light transport). Second, the spectral range can be chosen by appropriate selection of light source, SLM, and imaging sensor. Last, for many applications, an MI system has the potential to be very low cost, capitalizing on the widespread availability of consumer-grade digital cameras and projection systems. Here, we use a research-grade $16\text{-}\text{bit}$ CCD system, but the required dynamic range for many applications can be as low as $8\phantom{\rule{0.3em}{0ex}}\text{bits}$ , depending on the required reflectance intensity (and thus optical property) resolution.

## 3.

## Theory and Measurement in the SFD

## 3.1.

### Diffusion Approximation

The concept of a temporally modulated scalar photon density wave in turbid media is well established—its dispersive, diffractive, and interference properties have been widely studied and used for both quantitation and image formation.^{3, 4, 6, 12, 13, 14, 15} The notion of *spatially* modulated photon density “standing” waves, however, has mainly been considered as a theoretical construct (i.e., as the Fourier transform representation of spatial point sources and perturbations), as opposed to a practical measurement modality employing periodic illumination. Our goal here is to provide a simple conceptual framework to understand the fluence rate and reflectance properties of spatially modulated photon density plane waves in the SFD. We formulate this within a diffusion context and then later extend the discussion to transport-based Monte Carlo simulations in order to extend the applicability of MI to low albedo and high spatial frequency regimes.

The time-independent form of the diffusion equation for a homogeneous medium is given by

where $\phi $ is the fluence rate, $q$ is the source, ${\mu}_{\mathit{tr}}=({\mu}_{a}+{\mu}_{s}^{\prime})$ is the transport coefficient, ${\mu}_{\mathit{eff}}={\left(3{\mu}_{a}{\mu}_{\mathit{tr}}\right)}^{1\u22152}$ , ${\mu}_{a}$ is the absorption coefficient, ${\mu}_{s}^{\prime}={\mu}_{s}(1-g)$ is the reduced scattering coefficient, and $g$ is the cosine of the average scattering angle. Imposing a semi-infinite geometry, as depicted in Fig. 3, we introduce a normally incident, periodically varying plane wave source:with spatial frequencies (or repetencies) ${f}_{x}=({k}_{x}\u22152\pi )$ and ${f}_{y}=({k}_{y}\u22152\pi )$ , and spatial phases $\alpha $ and $\beta $ , extending infinitely in the tangential spatial dimensions, $x$ and $y$ , with some arbitrary dependence on depth, $z$ .Assuming a linear medium (i.e., a response proportional to the input intensity), this sinusoidal source will give rise to a diffuse fluence rate with the same frequency and phase. (From symmetry considerations, there should be no lateral phase shift^{16} for normally incident light onto a homogeneous medium.)

## Eq. 4

$\frac{{\mathrm{d}}^{2}}{\mathrm{d}{z}^{2}}{\phi}_{0}\left(z\right)-{\mu}_{\mathit{eff}}^{\prime 2}{\phi}_{0}\left(z\right)=-3{\mu}_{\mathit{tr}}{q}_{0}\left(z\right),$## Eq. 5

${\mu}_{\mathit{eff}}^{\prime}={({\mu}_{\mathit{eff}}^{2}+{k}_{x}^{2}+{k}_{y}^{2})}^{1\u22152}\equiv \frac{1}{{\delta}_{\mathit{eff}}^{\prime}}.$At zero spatial frequency
$(k=0)$
, the effective penetration depth
${\delta}_{\mathit{eff}}^{\prime}$
is equivalent to that of a planar (constant) illumination source,
${\delta}_{\mathit{eff}}=(1\u2215{\mu}_{\mathit{eff}})$
. In general, however,
${\mu}_{\mathit{eff}}^{\prime}$
(and
${\delta}_{\mathit{eff}}^{\prime}$
) are functions of *both* optical properties and spatial frequency of illumination. The 1-D form of Eq. 4 implies that the *amplitude* of the periodic wave,
${\phi}_{0}\left(z\right)$
, is independent of the tangential spatial dimensions
$x$
and
$y$
. As Eq. 4 is identical to the diffusion equation for a planar illumination, we can use existing planar geometry solutions by simply substituting
${\mu}_{\mathit{eff}}$
with our new
${\mu}_{\mathit{eff}}^{\prime}$
term.

As in the derivation of Svaasand
^{5} for planar photon density wave reflectance, we model an extended source as

## Eq. 6

${q}_{0}\left(z\right)={P}_{0}{\mu}_{s}^{\prime}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-{\mu}_{\mathit{tr}}z),$## 7.

^{15}the flux, $j$ , is set proportional to the fluence at the interface $z=0$ :

## 8.

## Eq. 9

$C=\frac{-{P}_{0}3{a}^{\prime}(1+3A)}{({\mu}_{\mathit{eff}}^{\prime 2}\u2215{\mu}_{\mathit{tr}}^{2}-1)({\mu}_{\mathit{eff}}^{\prime}\u2215{\mu}_{\mathit{tr}}+3A)},$## Eq. 10

${R}_{d}\left(k\right)=\frac{-{\phantom{\mid}j\mid}_{z\to {0}^{+}}}{{P}_{0}}=\frac{3A{a}^{\prime}}{({\mu}_{\mathit{eff}}^{\prime}\u2215{\mu}_{\mathit{tr}}+1)({\mu}_{\mathit{eff}}^{\prime}\u2215{\mu}_{\mathit{tr}}+3A)}.$For a given set of optical properties, the function ${R}_{d}\left(k\right)$ specifies the diffuse spatial modulation transfer function (MTF) of the medium. The simplicity of Eq. 10 allows some physical intuition of its properties. First, the frequency dependence of ${R}_{d}$ in the SFD is an inverse polynomial function of a single, positive-valued ratio, $({\mu}_{\mathit{eff}}^{\prime}\u2215{\mu}_{\mathit{tr}})$ , which fully describes the low-pass spatial filtering properties of homogeneous turbid samples within a steady-state diffusion context:

## Eq. 11

$\frac{{\mu}_{\mathit{eff}}^{\prime}}{{\mu}_{\mathit{tr}}}={\left(\frac{{\mu}_{\mathit{eff}}^{2}+{k}^{2}}{{\mu}_{\mathit{tr}}^{2}}\right)}^{1\u22152}={(3\frac{{\mu}_{a}}{{\mu}_{\mathit{tr}}}+\frac{{k}^{2}}{{\mu}_{\mathit{tr}}^{2}})}^{1\u22152}=\{\begin{array}{ll}\sqrt{3}{(1-{a}^{\prime})}^{1\u22152};& k=0\\ \frac{k}{{\mu}_{\mathit{tr}}};& k\u2aa2{\mu}_{\mathit{eff}}\end{array}\phantom{\}}.$The diffusion approximation to the radiative transport equation is valid when

and due to the anisotropic nature of light scattering, has been observed to be accurate approximately whenwhere $\rho $ describes the distance from collimated sources. Depending on the measurement technique (modality, geometry, calibration method, etc.) and metric of accuracy chosen, the practical*minimum*limit of $\rho $ is approximately in the range of $3{l}^{*}$ to $4{l}^{*}$ (Refs. 17, 18). The spatial frequency analog of the transport length ${l}^{*}$ is the transport spatial frequency, exactly equal to the transport spatial frequency (or transport coefficient), ${\mu}_{\mathit{tr}}={f}_{x,\mathit{tr}}=({k}_{\mathit{tr}}\u22152\pi )$ . If we relate the inverse of $\rho $ as a metric of spatial frequency, then Eq. 13 can be restated asWe therefore expect the

*maximum*spatial frequency accurately predicted by diffusion to be in the range of $1\u2215\left(3{l}^{*}\right)$ to $1\u2215\left(4{l}^{*}\right)$ , or $0.25{\mu}_{\mathit{tr}}$ to $0.33{\mu}_{\mathit{tr}}$ . Both albedo and source-distance requirements of the diffusion approximation limit the ratio $({\mu}_{\mathit{eff}}^{\prime}\u2215{\mu}_{\mathit{tr}})$ to be much less than one. In the following sections, we will highlight the quantitative power of this SFD diffusion model through (1) comparison to Monte Carlo simulations and (2) empirical demonstration of measurement accuracy in turbid phantom systems.

## 3.2.

### Monte Carlo Simulation in the SFD

Although more computationally intensive, it is desirable to have a forward model in the SFD that is valid for a greater range of both albedo and spatial frequency. A few transport-based approaches are available, including both direct numerical solution of the radiative transport equation^{19, 20} and Monte Carlo simulation.^{21, 22} For this paper, we have used “White” Monte Carlo (WMC) simulations^{23, 24} of a collimated point-source illumination to generate predictions of the spatially resolved, steady-state diffuse reflectance,
${R}_{d}\left(\rho \right)$
, for a given set of optical properties
${\mu}_{a}$
,
${\mu}_{s}$
, and
$g$
. This spatial point-spread function, provides an *impulse response*, and spatial frequency domain predictions of diffuse reflectance,
${R}_{d}\left(k\right)$
, are found by Fourier transformation of
${R}_{d}\left(\rho \right)$
. For a radially symmetric function such as
${R}_{d}\left(\rho \right)$
, the 2-D Fourier transform in the
$x\text{-}y$
plane reduces to a 1-D Hankel transform of order zero:

## Eq. 15

${R}_{d}\left(k\right)=2\pi \int \rho {J}_{0}\left(k\rho \right){R}_{d}\left(\rho \right)d\rho ,$## Eq. 16

${R}_{d}\left(k\right)=2\pi \sum _{i=1}^{n}{\rho}_{i}{J}_{0}\left(k{\rho}_{i}\right){R}_{d}\left({\rho}_{i}\right)\Delta {\rho}_{i}.$## 3.3.

### Simulations

Diffuse reflectance versus spatial frequency
$\left({\mathrm{mm}}^{-1}\right)$
, predicted by both the standard diffusion approximation (dashed lines) and Monte Carlo simulations (solid lines), is plotted in Fig. 4 (top) for varying values of
${l}^{*}$
, at a constant
$({\mu}_{s}^{\prime}\u2215{\mu}_{a})=100$
ratio (constant
${a}^{\prime}=0.99$
). Observe that as
${l}^{*}$
increases (or as
${\mu}_{\mathit{tr}}$
decreases), the diffuse MTF is rescaled toward lower spatial frequencies, indicating that less high-frequency content is preserved. This scaling with
${l}^{*}$
is consistent with our experience that high-scattering samples can retain very sharp (high-frequency) reflectance features. For example, reflectance from a point illumination is more localized in a high-scattering medium (like Spectralon, for instance), compared to a lower scattering medium (like *in vivo* tissue). Moreover, the frequency scaling of
${R}_{d}\left({f}_{x}\right)$
varies directly with
${\mu}_{\mathit{tr}}$
, or inversely with
${l}^{*}$
. This scaling of
${\mu}_{\mathit{tr}}$
and
${f}_{x}$
is directly evident in the
$({\mu}_{\mathit{eff}}^{\prime}\u2215{\mu}_{\mathit{tr}})$
ratio of Eq. 11 (diffusion approximation), and thus all five diffusion curves will coincide perfectly if plotted versus normalized spatial frequency,
$({f}_{x}\u2215{f}_{x,\mathit{tr}})=({f}_{x}\u2215{\mu}_{\mathit{tr}})=({f}_{x}\cdot {l}^{*})$
. This behavior is also retained in our Monte Carlo predictions to a high degree of accuracy. For instance, when plotted versus
$({f}_{x}\u2215{f}_{x,\mathit{tr}})$
(not shown), all five transport-based MTF curves fall within approximately 1% of each other, and this difference *decreases* further as the albedo is lowered. For all following simulations, therefore, we will plot the reflectance versus normalized spatial frequency,
$({f}_{x}\u2215{f}_{x,\mathit{tr}})$
. Conveniently,
${\mu}_{\mathit{tr}}=1\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$
$({l}^{*}=1\phantom{\rule{0.3em}{0ex}}\mathrm{mm})$
is a good approximate transport coefficient for many biological tissues, so for the high-albedo curves,
${f}_{x,\mathit{tr}}$
can be interpreted as
$\sim 1\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$
spatial frequency (
$1\text{-}\mathrm{mm}$
-spaced sinusoids).

Visual comparison of diffusion and Monte Carlo reflectance curves reveals that the diffusion solution slightly overestimates low-frequency components and underestimates the high-frequency components of the reflectance. This is partially due to our choice of a simple, mono-exponential extended source function [Eq. 6]. Analytical solutions that preserve higher order spatial moments of the source are available^{25, 26, 27} and will be examined in future work. Defining diffusion error as the percent difference of the diffusion prediction from Monte Carlo, we observe (for the case of
${\mu}_{s}^{\prime}\u2215{\mu}_{a}=100$
) a diffusion error of
$\pm 12\%$
for
${f}_{x}\u2a7d{\left(2{l}^{*}\right)}^{-1}=0.5{\mu}_{\mathit{tr}}$
.

In Fig. 4b, we further examine our transport and diffusion models when the albedo is decreased, ranging the $({\mu}_{s}^{\prime}\u2215{\mu}_{a})$ ratio from 30 (light gray lines) to 3 (black lines). Here, we plot the diffuse reflectance (top) and the diffusion error (middle) versus normalized spatial frequency, ${f}_{x,\mathit{tr}}$ . Again, diffusion overestimates reflectance at low spatial frequencies and underestimates reflectance at high frequencies. Furthermore, all diffusion lines seem to converge at high frequency $({f}_{x}\u2215{f}_{x,\mathit{tr}}\approx 1)$ , while an absolute offset between curves remains in the transport prediction.

For low frequencies (below $\sim 0.5\cdot {\mu}_{\mathit{tr}}$ ), the diffusion error remains less than $\pm 16\%$ at $({\mu}_{s}^{\prime}\u2215{\mu}_{a})=30$ $({a}^{\prime}=0.97)$ . For the lower albedo curves, we see that the absolute diffusion predictions are inaccurate, positively and negatively biased at low and high frequencies, respectively. In a real measurement, however, we always use a reference calibration sample (with known optical properties) to correct for these types of offsets. (See Sec. 3.4 for a detailed description of our calibration method.) We simulated this calibration procedure by generating forward Monte Carlo measurement data for both sample and calibration. The resulting diffusion error (not shown) using calibrated measurements of samples within $\pm 25\%$ of the reference phantom $({\mu}_{s}^{\prime}\u2215{\mu}_{a})$ exhibits a dramatic improvement in the shape and accuracy of measured data, reducing and flattening the diffusion error compared to the original “absolute” offsets. Specifically, we observe $<10\%$ error down to $({\mu}_{s}^{\prime}\u2215{\mu}_{a})=10$ for all frequencies. This result suggests that one can still achieve quantitatively accurate results through measurement calibration with a reference phantom of similar albedo.

## 3.4.

### Illumination, Imaging, and Calibration

The diffuse MTF of a turbid system can be measured in a transmission or reflection geometry. In practice, the illumination must be a superposition of AC (spatially modulated) and DC (planar) reflectance terms (i.e., we cannot illuminate with a negative scalar intensity). We therefore illuminate the sample with a spatial pattern of the form:

## Eq. 17

$S=\frac{{S}_{0}}{2}[1+{M}_{0}\phantom{\rule{0.2em}{0ex}}\mathrm{cos}(2\pi {f}_{x}x+\alpha )],$A host of signal processing schemes can be used to obtain
${M}_{\mathit{AC}}(x,{f}_{x})$
. Here, we employ a simple time-domain amplitude demodulation method,^{8, 28} illuminating a sinusoid pattern three times at the same spatial frequency, with phase offsets
$a=0$
,
$2\u22153\pi $
, and
$4\u22153\pi $
radians.
${M}_{\mathit{AC}}(x,{f}_{x})$
can then be calculated algebraically at each spatial location,
${x}_{i}$
, by:

## Eq. 20

${M}_{\mathit{AC}}({x}_{i},{f}_{x})=\frac{{2}^{1\u22152}}{3}{\{{[{I}_{1}\left({x}_{i}\right)-{I}_{2}\left({x}_{i}\right)]}^{2}+{[{I}_{2}\left({x}_{i}\right)-{I}_{3}\left({x}_{i}\right)]}^{2}+{[{I}_{3}\left({x}_{i}\right)-{I}_{1}\left({x}_{i}\right)]}^{2}\}}^{1/2},$## Eq. 21

${M}_{\mathit{DC}}\left({x}_{i}\right)=\frac{1}{3}[{I}_{1}\left({x}_{i}\right)+{I}_{2}\left({x}_{i}\right)+{I}_{3}\left({x}_{i}\right)].$In the frequency domain, a measurement
${M}_{\mathit{AC}}\left({f}_{x}\right)$
is the *product* of (1) the source intensity,
${I}_{0}$
; (2) the MTF of the illumination and imaging optical system,
${\mathit{MTF}}_{\mathit{\text{system}}}$
; and (3) the true turbid system MTF,
${R}_{d}$
:

## Eq. 22

${M}_{\mathit{AC}}({x}_{i},{f}_{x})={I}_{0}\cdot {\mathit{MTF}}_{\mathit{\text{system}}}({x}_{i},{f}_{x})\cdot {R}_{d}({x}_{i},{f}_{x}).$## Eq. 23

${R}_{d}({x}_{i},{f}_{x})=\frac{{M}_{\mathit{AC}}({x}_{i},{f}_{x})}{{M}_{AC,\mathit{ref}}({x}_{i},{f}_{x})}\cdot {R}_{d,\mathit{ref},\mathit{pred}}\left({f}_{x}\right).$This direct division-based correction for the system frequency response is an advantage of SFD measurement over other spatially resolved measurements, avoiding system PSF deconvolution in the real spatial domain, which can amplify measurement noise and uncertainties. Ideally, the surface contours of the sample and the phantom should be identical or be compensated numerically using surface profilometry.^{29}

Last, for a given modulation frequency, there are two unknowns in Eq. 10: ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ . We show here how measurements at as few as two spatial frequencies can be used to separate absorption and scattering. This is best visualized in a lookup table such as the one shown in Fig. 6, where ${R}_{d}\left(\mathit{DC}\right)$ and ${R}_{d}\left(\mathit{AC}\right)$ correspond to diffuse reflectance measurements at zero and nonzero spatial frequencies ${f}_{1}$ and ${f}_{2}$ , respectively. Gray and black contours correspond to constant absorption and reduced scattering, respectively. As a visual example, the dotted lines in Fig. 6 show that if ${R}_{d}\left(0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}\right)=0.55$ and ${R}_{d}\left(0.5\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}\right)=0.06$ , then ${\mu}_{a}\approx 0.03$ and ${\mu}_{s}^{\prime}\approx 1.4\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ , respectively. Notice the strongly orthogonal relationship between the absorption and scattering contour lines, indicating the ability to separate absorption and scattering values with maximal sensitivity. This is due to the large frequency range spanned by $0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ and $0.5\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ ( $\mathit{DC}$ and $\mathit{AC}$ ) frequencies. Correspondingly, as the $x$ - and $y$ -axis frequencies become closer to one another, these lines will become less orthogonal, and inversion coupling between absorption and scattering will increase. Last, both $\mathit{AC}$ and $\mathit{DC}$ measurements can be easily obtained with only three phase projections of a single illumination frequency [through Eqs. 20, 21], allowing rapid, high-resolution imaging of absorption and scattering contrast.

## 3.5.

### Inversion Methods

We use two inversion methods to calculate the absorption and reduced scattering from measurements of diffuse reflectance. First we use a “sweep” in spatial frequency space, analogous to the broadband frequency domain photon migration (FDPM) approach,^{30} producing an overdetermined set of measurements that can be fit to Eq. 10 via least-squares minimization. Second, we use a rapid two-frequency lookup table method, introduced in the previous paragraph, which uses cubic spline interpolation (the “griddata” method in MATLAB) of forward-model data at two spatial frequencies. On typical personal computers, this method is capable of millions of lookup calculations per second. In this initial work, we invert MTF measurements at each spatial location independently with our spatially homogeneous planar model. We acknowledge that while accurate for homogeneous media, this approach ignores any transverse or depth-resolved transport phenomena. We therefore expect significant partial volume effects in the recovered data in regions where the optical properties are spatially varying. For the remainder of this paper, we discuss optical property *maps* in terms of qualitative optical property contrast, in comparison to quantitative optical properties when measuring large spatially homogenous regions. In Sec. 5.4, we investigate further the resolution limits and partial volume effects of absorption and scattering contrast in the axial and transverse directions. For lateral step-function changes in optical properties, we find that measurements transition spatially between two quantitatively accurate values with a sigmoidal-like transition region.

Figure 7 visually depicts the entire data-mining process using the *in vivo* forearm data described in Sec. 4.2. Intensity data at each frequency (three phase images per frequency) are demodulated, calibrated, and fit using Eqs. 20, 23, 10, respectively. Data are processed separately for each pixel, generating spatial maps of absorption and scattering optical properties.

Compared to other spatially resolved methods, MI acquires coincident, axial “projection” measurements of optical contrast to quantify the optical properties at each $x\text{-}y$ spatial position, allowing a robust measurement of the average properties. Compared to “point” illumination measurements, MI samples only the low spatial frequency moments of the transfer function. These low frequencies $(<1\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1})$ are sufficient for separating the absorption from the scattering optical properties, reducing sensitivity to uncertainty inherent in measuring high-frequency spatial moments (i.e., reflectance close to the source).

## 4.

## Methods

## 4.1.

### Phantom Experiments

We performed a set of experiments to characterize the precision and accuracy of MI for measuring the homogeneous absorption and reduced scattering optical properties. Sixteen turbid phantoms were constructed using a single batch of Liposyn lipid emulsion and water-soluble nigrosin dye stock solutions for the scattering and absorbing properties, respectively. In the first eight phantoms, we varied the absorption coefficient,
${\mu}_{a}$
, over two orders of magnitude (logarithmically spaced between
$0.002\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}\u2a7d{\mu}_{a}\u2a7d0.12\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$
), with a constant scattering coefficient constant at
${\mu}_{s}^{\prime}=0.97\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$
. In the second set, we linearly varied
${\mu}_{s}^{\prime}$
$(0.32\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}\u2a7d{\mu}_{s}^{\prime}\u2a7d1.8\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1})$
, while holding the absorption coefficient constant at
${\mu}_{a}=0.0046\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$
. These values were calculated based on infinite-geometry, multifrequency
$\left(50\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}500\phantom{\rule{0.3em}{0ex}}\mathrm{MHz}\right)$
, multidistance
$\left(10\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}25\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\right)$
frequency-domain photon migration measurements^{15} of one of the Liposyn/nigrosin phantoms.

MI measurements were performed on each sample. Thirty spatial frequencies of illumination were chosen between $0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ and $0.13\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ , corresponding to a total of 90 images per phantom (three spatial phases per frequency). The interfrequency spacing was chosen to accurately capture the MTF curvature of all phantoms, with finer spacing at low frequencies and coarser spacing at high frequencies, accordingly. All measurements were taken at $660\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ with an approximate $75\times 75\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ illumination area, a $50\times 50\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ camera field of view, and an integration time of $100\phantom{\rule{0.3em}{0ex}}\mathrm{ms}$ . The individual phantoms were measured in a randomized order, and measurements were repeated three times to allow for statistical averaging.

Modulation images of the AC reflectance were obtained at each frequency using Eq. 20. At full CCD resolution, the pixel-by-pixel demodulation approach results in approximately 250,000 separate measurements of reflectance per spatial frequency, highlighting the statistical power of the technique. As the lipid solutions were expected to be highly homogeneous, $20\times 20\phantom{\rule{0.3em}{0ex}}\text{pixel}$ binning was performed on each image to speed computation, resulting in low-resolution, $15\times 15\phantom{\rule{0.3em}{0ex}}\text{pixel}$ modulation images. The resulting 30 images provide a quantitative AC amplitude measurement at each of 100 spatial locations within the field of view. For calibration, a single phantom from the entire set of 16 was chosen as the reference (second-lowest absorption phantom). Using the reference’s known optical properties (determined from infinite-geometry FDPM measurements), we calculate a model-based prediction for the reflectance, ${R}_{d,\mathit{ref},\mathit{pred}}\left({f}_{x}\right)$ . Then, for each spatial frequency and each spatial location, we use Eq. 23 to calculate ${R}_{d}\left({f}_{x}\right)$ of the sample. Having retained some low-resolution spatial data, we can calculate a standard deviation of recovered values within an image as an indicator of measurement precision.

The diffusion model of Eq. 10 was used to solve for ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ using both least-squares minimization by a simplex search algorithm (in “fminsearch” MATLAB) and via the two-frequency lookup table approach using the lowest $\left(0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}\right)$ and highest $\left(0.13\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}\right)$ spatial frequencies. For each phantom, each spatial sampling point was separately analyzed, generating images of recovered absorption and scattering values. As these were homogeneous samples, a mean and a standard deviation were calculated to represent each optical property image result, characterizing the accuracy and precision of MI, respectively.

## 4.2.

###
*In Vivo* Human Forearm Experiments

MI data were collected on a normal human forearm over a $72\times 48\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ field of view. Four evenly spaced spatial frequencies between 0 and $0.15\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ were collected and analyzed. The imaging system was identical to that described earlier, except for the inclusion of a $640\pm 10\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ bandpass detection filter and crossed linear polarizers, which reject specular reflection from rough surfaces and maximize our sensitivity to the diffuse component of the light. In idealized liquid phantom experiments, we have performed measurements with and without crossed polarizers and found the difference in recovered optical properties to be typically less than 2 to 3%.

In order to demonstrate the sensitivity of our system to physiological perturbations, we performed a standard venous occlusion study on a $29\times 40\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ region of the volar forearm. Measurements were performed at a wavelength of $800\pm 10\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ , near the hemoglobin isosbestic point of $805\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ . Measured changes in absorption at this wavelength are insensitive to oxygenation and therefore reflect only that of total hemoglobin. Multifrequency reflectance data at 0 and $0.135\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ were acquired every $4\phantom{\rule{0.3em}{0ex}}\mathrm{s}$ for a period of $13\phantom{\rule{0.3em}{0ex}}\mathrm{min}$ . After $2.5\phantom{\rule{0.3em}{0ex}}\mathrm{min}$ of baseline acquisition, an arm cuff was pressurized to $100\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ Hg for $6.5\phantom{\rule{0.3em}{0ex}}\mathrm{min}$ and subsequently released at minute 9.

## 5.

## Results and Discussion

## 5.1.

### Phantom Experiment Results

The average measured diffuse reflectance versus spatial frequency is plotted in Fig. 8, showing the absorption variation and scattering variation measurement sets in Figs. 8a and 8b, respectively. In solid black lines, we show the corresponding fits using the diffusion-based reflectance model [Eq. 10].

The absorption experiment data demonstrate that increasing absorption causes a *decrease* in reflectance, with absorption contrast residing primarily in the low-frequency regime. (Notice that all absorption curves converge at high frequency.) Conversely, the scattering data indicate that increasing scattering causes an *increase* in reflectance amplitude and a rescaling to higher spatial frequencies (i.e., a decrease in
${l}^{*}$
), with contrast apparent at all spatial frequencies.

All model-based fits of Fig. 8 (solid lines) demonstrate excellent visual agreement with the data, with typical errors less than 0.02. This is particularly satisfying, as all measurements were calibrated with a single reference phantom (second-lowest absorption phantom). The largest model-data deviation appears in the high-frequency range of the lowest scattering phantom ( ${\mu}_{s}^{\prime}=0.32\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ , ${l}^{*}\approx 3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ ). This seems consistent with the ${l}^{*}$ plots of Fig. 4, where we would expect model breakdown at or before ${f}_{x}=1\u2215\left(2{l}^{*}\right)$ , or $0.16\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ .

In Fig. 9, we plot on the left and right the recovered optical properties for absorption and scattering variation measurements, respectively. Multifrequency and two-frequency lookup table interpolation results are shown in black and gray, respectively. For each set, the varied value is plotted versus the expected value on the horizontal axis, and the corresponding value held constant is shown below on a separate axis. Error bars indicate the corresponding standard deviations of the recovered $15\times 15\phantom{\rule{0.3em}{0ex}}\text{pixel}$ optical property maps (not shown). Thin black lines are drawn to indicate the expected values in each experiment.

In the absorption variation experiment, recovered versus expected absorption shows excellent linearity over two orders of magnitude, ranging from high to low albedo ( ${\mu}_{s}^{\prime}\u2215{\mu}_{a}=500$ to ${\mu}_{s}^{\prime}\u2215{\mu}_{a}=8$ ). When absorption is very small, a slight overestimation trend is observed independent of calibration choice. This trend is discussed further in Sec. 5.2. The experiment’s recovered scattering values show less than 10% deviation from the expected value in all cases. Similar linearity is observed in the scattering variation experiment, albeit with slightly more fluctuation. Absorption values in this case demonstrate less than 15% deviation from the expected value, except in the lowest scattering (largest ${l}^{*}$ ) case. Standard deviations of the recovered $15\times 15\phantom{\rule{0.3em}{0ex}}\text{pixel}$ optical property maps are mostly less than 1% (maps not shown), indicating both high precision and spatial uniformity over the field of view.

In Table 1, we summarize the accuracy of all recovered values by showing the average percent deviation for multifrequency and two-frequency lookup table inverse models. On average, we have 6% and 3% deviation in absorption and reduced scattering, respectively. The two-frequency lookup table errors are generally comparable to those of the multifrequency method. In real-world measurements of spatially heterogeneous systems, we expect that multifrequency fits will provide a more stable measure of the average optical properties. Nevertheless, in situations where speed and/or resolution is desired, this method shows promise to provide rapid feedback while retaining quantitative accuracy.

## Table 1

Summary of recovered average optical properties for absorption and scattering variation experiments. Overall accuracies of recovered absorption and reduced scattering parameters are approximately 6% and 3%, respectively.

Absorption experiment | Scattering experiment | ||||
---|---|---|---|---|---|

Average error (%) | 30-frequency fit | 2-frequency fit | Average error (%) | 30-frequency fit | 2-frequency fit |

${\mu}_{a}$ error | 4.74 | 4.85 | ${\mu}_{a}$ error | 7.51 | 11.4 |

${\mu}_{s}^{\prime}$ error | 2.98 | 2.29 | ${\mu}_{s}^{\prime}$ error | 3.05 | 10.2 |

## 5.2.

### Separation of Absorption and Scattering

In our inverse fitting results, the largest errors were observed for the lowest values of absorption and scattering. As a planar imaging modality, MI samples relatively superficial volumes and therefore short photon paths losing sensitivity to low absorption and low scattering contrast where the length scales of photon interaction are very long. Furthermore, the finite size of projection and illumination set a physical limit on the low-frequency MTF sampling. For instance, in the absorption experiment, we observed an overestimation of the absorption when absorption was very low
$\left(0.002\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}\right)$
. In this regime, both diffusion and Monte Carlo models predict the diffuse reflectance to have a steep, decreasing slope at low frequency. However, the lowest four experimental illumination frequencies (including
${f}_{x}=0$
) correspond to spatial periods *larger* than the projector’s illumination area
$({f}_{x}<0.013\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1})$
. As sampling theory in this domain enforces a frequency bandwidth greater than the fundamental illumination frequency, we expect inaccuracy in these lowest frequency measures. Specifically, based on the low-pass MTF shape of
${R}_{d}$
, we expect amplitude *underestimation*, and therefore absorption *overestimation*, due to the high sensitivity of absorption at low frequencies. One strategy to systematically account for this effect is to equivalently model an illumination with a finite spatial extent.

## 5.3.

###
*In Vivo* Human Forearm Results

In Fig. 7 (middle), we showed diffuse reflectance images
$\left({R}_{d}\right)$
versus spatial frequency
$\left({f}_{x}\right)$
for the *in vivo* human forearm experiment. Notice the differential contrast in diffuse reflectance as illumination frequency increases, forming the basis for separation of absorption and scattering. In addition, high frequencies will sample a more superficial region of the tissue, which is expected to have a lower contribution from deeper vascular features. In Fig. 10, we further show the recovered optical property maps after multifrequency MTF fitting at each pixel. In Fig. 10a, we show the calibrated diffuse reflectance at
${f}_{x}=0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$
, and in Fig. 10d, the average multifrequency diffuse reflectance data (black squares) and fit (gray line). In Figs. 10b and 10c, we show spatial maps of the absorption and reduced scattering data, respectively, and in Figs. 10e and 10f, we show the corresponding pixel histograms, respectively. Absorption contrast from the underlying veins is the dominant feature in the optical property maps. A vertical feature of lowered scattering appears in the middle of the image. This feature is coincident with a large superficial tendon, which may be acting effectively as a light guide due to its generally higher index compared to tissue matrix.^{31, 32} Based on a diffusion-based sensitivity analysis, we predict for this experiment 1/e sampling depths ranging from
$2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}3.3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
for low
$\left(0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}\right)$
and high
$\left(0.15\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}\right)$
spatial frequencies, respectively.

In the absorption map, we identify with dotted lines a region of interest containing a prominent vein. In this region, we observe 100% contrast in the recovered absorption values over the vein
$\left(0.46\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}\right)$
compared to the background
$\left(0.23\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}\right)$
; in the same region, optical scattering showed little contrast, with
$\sim 10\%$
spatial variation. While these data show clear separation of absorption and scattering optical contrast, the values derived from a homogeneous model exhibit partial volume effects that diminish the actual absorption contrast of the vein beneath the surface. Our current homogeneous model of reflectance prevents absolute quantitation in the presence of lateral and depth-dependent optical properties. In the next section, however, we attempt to empirically assess these partial volume effects through *edge-response* experiments in gelatin phantoms.

Venous occlusion measurements on the volar forearm are shown in Fig. 11 and demonstrate the accumulation and dissipation of blood volume via the absorption coefficient at $800\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ . Figure 11a shows the diffuse reflectance map (left) and optical absorption map (right) measured at the baseline. Two regions of interest were chosen to show both global changes over the entire image (gray lines) and changes in region absent of any obvious large vessels (black lines), presumably containing only microvasculature. Region-wise average changes in optical properties from the baseline are shown in Fig. 11b. Absorption (top) and reduced scattering (bottom) are respectively plotted within 40% and 10% of the corresponding measured baseline values. As expected, the absorption in either region begins to increase steadily after arm cuff pressurization at minute 2.5. After release of the cuff at $9\phantom{\rule{0.3em}{0ex}}\mathrm{min}$ , absorption decreases to the baseline over the course of approximately $2\phantom{\rule{0.3em}{0ex}}\mathrm{min}$ . Maximal increases in absorption were observed to be approximately $0.012\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ globally and $0.017\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ in the microvascular region, corresponding to approximate 15% and 28% increases from the baseline. The larger increase in absorption observed in the microvascular region may be explained by the fact that the microvasculature is more susceptible to pooling, while the larger vessels are less reactive. Small fluctuations in the measured reduced scattering were observed to be $<5\%$ globally and $<2\%$ for the microvascular region.

## 5.4.

### Sensitivity, Contrast, and Resolution

Detecting changes in ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ requires the corresponding change in reflectance to be above the measurement noise floor. We examined how a homogeneous perturbation in optical properties gives rise to reflectance contrast at each spatial frequency using our Monte Carlo forward model. This was done by taking a numerical derivative of reflectance with respect to a change in a given optical property, generated for normalized frequencies between 0 and 1. In Fig. 12, we show the change in diffuse reflectance $\left(\Delta {R}_{d}\right)$ versus normalized spatial frequency resulting from a 1% change in absorption or scattering, for sample $({\mu}_{s}^{\prime}\u2215{\mu}_{a})$ ratios of 100 (black lines) and 3 (gray lines). The solid and dashed lines denote the sensitivity to absorption and reduced scattering, respectively. (Note: We use a negative absorption perturbation and a positive scattering perturbation in order to produce reflectance contrast of the same polarity.)

The sensitivity profiles in Fig. 12 reveal how changes in absorption and scattering change the reflected light at each modulation frequency, indicating that absorption contrast is attenuated more rapidly with frequency compared to scattering. This difference is intuitive physically, as high-frequency illumination should sample only short path-length phenomena, losing sensitivity to long length-scale processes such as absorption. From the graph, we observe that 1% change in absorption or scattering produces at most an approximate 0.3% change in ${R}_{d}$ . We have added dotted lines to Fig. 12 to show the approximate corresponding camera detection limits with 12- and $14\text{-}\text{bit}$ intensity resolution (i.e., the detection limit for the physical contrast—actual detection limits will depend on the noise characteristics of the particular camera/imaging system). The figure illustrates, for example, that with $({\mu}_{s}^{\prime}\u2215{\mu}_{a})=100$ at ${f}_{x}\u2215{f}_{x,\mathit{tr}}=0.1$ , a $12\text{-}\text{bit}$ measurement can resolve (in intensity) a 1% change in scattering but not absorption, while a $14\text{-}\text{bit}$ camera can resolve both. For $({\mu}_{s}^{\prime}\u2215{\mu}_{a})=3$ , where optical properties are closer in magnitude, this $12\text{-}\text{bit}$ absorption detection criterion occurs at a higher frequency of ${f}_{x}\u2215{f}_{x,\mathit{tr}}=0.3$ .

The intensity resolution of the measurement can be further improved through either spatial binning or averaging multiple acquisitions. Imaging with a high-resolution camera, therefore, allows flexibility between spatial resolution and intensity (optical property) resolution. Sources of noise that may limit measurement precision and accuracy include light-source fluctuation (jitter), long-time source intensity stability, and spatial nonuniformity of the projection. In our measurements, we approximate a source fluctuation of approximately 0.1% for $1\text{-}\mathrm{s}$ integration times. Both intensity stability ( $\sim 10\%$ decrease over $4\phantom{\rule{0.3em}{0ex}}\mathrm{h}$ ) and spatial nonuniformity (20%) were also present but were corrected to first order by using periodic reflectance standard measurements and phantom calibration (see Sec. 3.4). In the case of our human forearm measurements, the background $({\mu}_{s}^{\prime}\u2215{\mu}_{a})$ ratio was approximately 100. In the vein region, we observed absorption contrast that is clearly resolved by our CCD camera.

In order to further understand the lateral and depth-dependent partial volume effects, such as those observed in the arm mapping experiments, we performed an exploratory contrast-resolution study using heterogeneous optical phantoms. Eight homogeneous gelatin phantoms were fabricated using nigrosin as the absorber and Liposyn as the scattering agent. Heterogeneous phantoms, shown in Fig. 13a, were assembled by placing two gelatin slabs of differing optical properties adjacent to one another. Gelatin phantoms with a 300% step in either absorption [Fig. 13a, left; ${\mu}_{a,\mathit{\text{left}}}=0.01\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ , ${\mu}_{a,\mathit{\text{right}}}=0.03\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ ; ${\mu}_{s,\mathit{\text{both}}}^{\prime}=1.0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ ] or scattering [Fig. 13a, right; ${\mu}_{a,\mathit{\text{both}}}=0.02\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ ; ${\mu}_{s,\mathit{\text{left}}}^{\prime}=0.5\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ , ${\mu}_{s,\mathit{\text{right}}}^{\prime}=1.5\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ ] were measured at $660\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ both directly [Fig. 13a, top] and through a $2\text{-}\mathrm{mm}$ homogeneous layer [Fig. 13a, bottom] with optical properties ${\mu}_{a}=0.01\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ , ${\mu}_{s}^{\prime}=0.5\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ . Both were calibrated by a homogeneous phantom with optical properties ${\mu}_{a}=0.02\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ , ${\mu}_{s}^{\prime}=1.0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ . Nine spatial frequencies between $0\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ and $0.11\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$ were measured and used to calculate optical property maps using least-squares regression to our diffusion reflectance model.

Recovered optical property spatial profiles were averaged over the vertical direction, shown in Fig. 13b for absorption (top) and scattering (bottom) media. The results reveal a diffuse *edge-response* function that is both depth and optical property dependent. For both absorption and scattering experiments, we observe a degradation of spatial resolution and quantitative contrast through the homogeneous layer (gray lines) compared to that at the surface (black lines). Specifically, the measured contrast values through the homogeneous layer are approximately 15% and 5% of those at the surface, for absorption and scattering, respectively. This level of absorption contrast is comparable with the measured arm vein absorption in Fig. 10 (approximately
$0.06\phantom{\rule{0.3em}{0ex}}{\mathrm{mm}}^{-1}$
), which corresponds to approximately 18% of whole blood absorption (assuming 70% tissue oxygen saturation). We further observe that the spatial resolution of the recovered scattering maps is consistently better than that for absorption. Put differently, scattering profiles appear to preserve higher spatial frequencies than those of absorption. This property of the measured heterogeneous media is consistent with the homogeneous phantom MTF results shown in Fig. 8, where we noted that absorption contrast appeared mainly at low frequencies, while scattering contrast appeared at all measured frequencies.

Defining spatial resolution as the distance at which the edge-response contrast is reduced by a 90% (all reflectance contrast remained well above our system’s minimum *intensity* resolution), we determined the resolution of absorption and scattering contrast to be
$0.3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
and
$0.05\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
, respectively, for surface perturbations, and
$0.5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
and
$0.25\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
, respectively, for perturbations
$2\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
beneath the surface. Although these changes exhibit characteristics that suggest the capability to resolve optical contrast on small spatial scales, actual performance will also be dependent on illumination spatial frequency and system noise. In future work, we aim to perform more rigorous contrast-resolution analyses using heterogeneous computational models and phantoms to determine quantitative resolution limits as a function of depth. In general, at a given depth, we expect the depth resolution to scale with the measurement precision and number of frequencies (# of *sources*), and
$x\text{-}y$
resolution to scale with the number of spatial sampling points (# of detectors).^{33}

## 6.

## Conclusion

We have presented a theoretical framework and instrumental platform for SFD measurement in turbid media, compared analytic diffusion and Monte Carlo–based transport forward models, and performed quantitative measurements of phantom systems and an *in vivo* human forearm. In the phantom validation measurements, we demonstrate excellent accuracy in optical properties (approximately 6% and 3% in absorption and reduced scattering, respectively) over a wide range of
${l}^{*}$
values
$\left(0.5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\right)$
and
$({\mu}_{s}^{\prime}\u2215{\mu}_{a})$
ratios (8 to 500). In the *in vivo* forearm spatial mapping, we report both imaging of optical absorption and scattering contrast with millimeter-scale resolution and dynamic imaging of physiological perturbations in blood volume due to venous occlusion.

A modulated imaging (MI) system can obtain quantitative optical properties in turbid homogeneous systems and maps of optical property contrast in tissues with a noncontact reflectance measurement. This combination of advantages make it particularly suited to imaging of static and dynamic processes in *in vivo* biological tissues, particularly for the field medical diagnostics. In ongoing studies, we are extending the method to multispectral imaging for quantitative functional mapping of both intrinsic and extrinsic tissue chromophores and evaluating depth-resolved imaging models, including both multilayer and tomographic approaches.

## Acknowledgments

We would like to express our sincere gratitude to all the people who provided advice and assistance toward this work, including Sophie Chung, Carole Hayakawa, Andrea Bassi, Katya Bhan, and Jae Gwan Kim and for support from the Laser Microbeam and Medical Program (LAMMP), an NIH Biomedical Technology Resource, Grant No. P41-RR01192; the U.S. Air Force Office of Scientific Research, Medical Free-Electron Laser Program (F49620-00-2-0371 and FA9550-04-1-0101); and the Beckman Foundation. Dr. Cuccia would also like to acknowledge support from the NSF Graduate Research Fellowship Program.

## References

**,” Lasers Surg. Med., 6 (6), 488 –493 (1987). https://doi.org/10.1002/lsm.1900060603 0196-8092 Google Scholar**

*Practical models for light distribution in laser-irradiated tissue***,” Opt. Lett., 15 (5), 276 –278 (1990). 0146-9592 Google Scholar**

*Determination of the scattering and absorption lengths from the temporal profile of a backscattered pulse***,” Phys. Rev. Lett., 69 (18), 2658 –2661 (1992). https://doi.org/10.1103/PhysRevLett.69.2658 0031-9007 Google Scholar**

*Refraction of diffuse photon density waves***,” J. Opt. Soc. Am. A, 10 (1), 127 –140 (1993). 0740-3232 Google Scholar**

*Propagation of photon-density waves in strongly scattering media containing an absorbing semi-infinite plane bounded by a straight edge***,” Phys. Med. Biol., 44 (3), 801 –813 (1999). https://doi.org/10.1088/0031-9155/44/3/020 0031-9155 Google Scholar**

*Reflectance measurements of layered media with diffuse photon-density waves: a potential tool for evaluating deep burns and subcutaneous lesions***,” J. Opt. Soc. Am. A, 4 (3), 423 –432 (1987). https://doi.org/10.1364/JOSAA.4.000423 0740-3232 Google Scholar**

*Model for photon migration in turbid biological media***,” Med. Phys., 19 (4), 879 –888 (1992). https://doi.org/10.1118/1.596777 0094-2405 Google Scholar**

*A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties**in vivo***,” Opt. Lett., 22 (24), 1905 –1907 (1997). https://doi.org/10.1364/OL.22.001905 0146-9592 Google Scholar**

*Method of obtaining optical sectioning by using structured light in a conventional microscope***,” Lasers Med. Sci., 13 55 –65 (1998). 0268-8921 Google Scholar**

*Determination of tissue optical properties by steady-state spatial frequency-domain reflectometry***,” Opt. Lett., 30 (11), 1354 –1356 (2005). https://doi.org/10.1364/OL.30.001354 0146-9592 Google Scholar**

*Modulated imaging: quantitative analysis and tomography of turbid media in the spatial-frequency domain***,” Opt. Express, 14 6516 –6534 (2006). https://doi.org/10.1364/OE.14.006516 1094-4087 Google Scholar**

*Noncontact fluorescence optical tomography with scanning patterned illumination***,” Appl. Opt., 28 2331 –2336 (1989). https://doi.org/10.1364/AO.28.002331 0003-6935 Google Scholar**

*Time-resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties***,” Proc. Natl. Acad. Sci. U.S.A., 91 (11), 4887 –4891 (1994). https://doi.org/10.1073/pnas.91.11.4887 0027-8424 Google Scholar**

*Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytic solution and applications***,” Phys. Rev. E, 61 (4), 4295 –4309 (2000). https://doi.org/10.1103/PhysRevE.61.4295 1063-651X Google Scholar**

*Near-field diffraction tomography with diffuse photon density waves***,” J. Opt. Soc. Am. A, 11 (10), 2727 –2741 (1994). https://doi.org/10.1364/JOSAA.11.002727 0740-3232 Google Scholar**

*Boundary conditions for the diffusion equation in radiative transfer***,” J. Opt. Soc. Am. A, 25 (11), 2833 –2839 (2008). https://doi.org/10.1364/JOSAA.25.002833 0740-3232 Google Scholar**

*Spatial shift of spatially modulated light projected on turbid media***,” Phys. Rev. Lett., 64 (22), 2647 –2650 (1990). https://doi.org/10.1103/PhysRevLett.64.2647 0031-9007 Google Scholar**

*When does the diffusion approximation fail to describe photon transport in random media***,” Phys. Rev. E, 53 (3), 2307 –2319 (1996). https://doi.org/10.1103/PhysRevE.53.2307 1063-651X Google Scholar**

*Gigahertz photon density waves in a turbid medium: theory and experiments***,” J. Opt. Soc. Am. A, 21 (5), 820 –827 (2004). https://doi.org/10.1364/JOSAA.21.000820 0740-3232 Google Scholar**

*Transport theory for light propagation in biological tissue***,” Optical-Thermal Response of Laser-Irradiated Tissue, 73 –99 Plenum, New York (1995). Google Scholar**

*Monte Carlo modeling of light transport in tissues***,” Opt. Eng., 34 2055 –2063 (1995). https://doi.org/10.1117/12.204798 0091-3286 Google Scholar**

*Determination of reduced scattering and absorption coefficients by a single charge-coupled-device array measurement. Part I: Comparison between experiments and simulations***,” Phys. Med. Biol., 41 2221 –2227 (1996). https://doi.org/10.1088/0031-9155/41/10/026 0031-9155 Google Scholar**

*Determination of the optical properties of turbid media from a single Monte Carlo simulation***,” J. Opt. Soc. Am. A, 20 (4), 714 –727 (2003). https://doi.org/10.1364/JOSAA.20.000714 0740-3232 Google Scholar**

*Accelerated Monte Carlo models to simulate fluorescence spectra from layered tissues***,” J. Atmos. Sci., 33 (12), 2452 –2459 (1976). https://doi.org/10.1175/1520-0469(1976)033<2452:TDEAFR>2.0.CO;2 0022-4928 Google Scholar**

*Delta-Eddington approximation for radiative flux-transfer***,” Appl. Opt., 39 (34), 6453 –6465 (2000). https://doi.org/10.1364/AO.39.006453 0003-6935 Google Scholar**

*Collimated light sources in the diffusion approximation***,” J. Biomed. Opt., 9 (3), 632 –647 (2004). https://doi.org/10.1117/1.1695412 1083-3668 Google Scholar**

*Radiative transport in the delta-P1 approximation: accuracy of fluence rate and optical penetration depth predictions in turbid semi-infinite media***,” Frontiers in Optics 2008/Laser Science (LS) XXIV, FWP2 (2008) Google Scholar**

*Fluorescence lifetime and spatially modulated light for image-guided surgery (invited)***,” Rev. Sci. Instrum., 7 2500 –2513 (2000). https://doi.org/10.1063/1.1150665 0034-6748 Google Scholar**

*A broad bandwidth frequency domain instrument for quantitative tissue optical spectroscopy***,” Protoplasma, 59 472 –479 (1964). 0033-183X Google Scholar**

*The part played by the mucopolysaccharides in the form birefringence of the collogen***,” Biopolymers, 20 (6), 1271 –1290 (1981). https://doi.org/10.1002/bip.1981.360200613 0006-3525 Google Scholar**

*Optical second-harmonic scattering in rat-tail tendon***,” Phys. Rev. E, 70 (5), 056616 (2004). https://doi.org/10.1103/PhysRevE.70.056616 1063-651X Google Scholar**

*Symmetries, inversion formulas, and image reconstruction for optical tomography*