Quantitation and mapping of tissue optical properties using modulated imaging.
Journal: 2009/July - Journal of Biomedical Optics
ISSN: 1083-3668
Abstract:
We describe the development of a rapid, noncontact imaging method, modulated imaging (MI), for quantitative, wide-field characterization of optical absorption and scattering properties of turbid media. MI utilizes principles of frequency-domain sampling and model-based analysis of the spatial modulation transfer function (s-MTF). We present and compare analytic diffusion and probabilistic Monte Carlo models of diffuse reflectance in the spatial frequency domain. Next, we perform MI measurements on tissue-simulating phantoms exhibiting a wide range of l values (0.5 mm to 3 mm) and (micro(s) (')micro(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.
Relations:
Content
Citations
(136)
References
(24)
Grants
(1K+)
Organisms
(1)
Processes
(1)
Similar articles
Articles by the same authors
Discussion board
J Biomed Opt 14(2): 024012

Quantitation and mapping of tissue optical properties using modulated imaging

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 (t-PSF), or spreading of a propagating pulse in time,12 and the latter measuring the temporal modulation transfer function (t-MTF), or the attenuation and phase delay of a periodically varying photon density wave.35 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.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f1.jpg

Four measurement domains of turbid media characterization: time domain (top left), time frequency domain (bottom left), real spatial domain (top right), and spatial frequency domain (bottom right).

In diffuse optics, spatially resolved measurements have been generally limited to the real spatial domain. Here, the spatial point-spread function (s-PSF) is typically characterized by “multidistance” measurements,67 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. 911

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 (s-MTF) 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, λ=660 nm, Δλ=10 nm FWHM). To create the illumination patterns, 8-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-bit frame-transfer CCD camera (Roper Cascade 512F) capable of imaging up to 30 frames per second at full 512×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.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f2.jpg

Modulated imaging instrument platform.

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-bit CCD system, but the required dynamic range for many applications can be as low as 8 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.3461215 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

2φμeff2φ=3μtrq,
(1)

where φ is the fluence rate, q is the source, μtr=(μa+μs)is the transport coefficient, μeff=(3μaμr), μa is the absorption coefficient, μs=μs(1g)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:

q = q0(z)cos(kxx + α)cos(kyy + β),
(2)

with spatial frequencies (or repetencies) fx=(kx/2π) and fy=(ky/2π), and spatial phases α and β, extending infinitely in the tangential spatial dimensions, x and y, with some arbitrary dependence on depth, z.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f3.jpg

Schematic of modulated illumination source (in the x direction only) and the resulting modulated internal fluence rate with the same frequency and phase.

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 shift16 for normally incident light onto a homogeneous medium.)

φ = φ0(z)cos(kxx + α)cos(kyx + β).
(3)

Insertion of Eqs. (2) and (3) into Eq. (1) yields a 1-D second-order Helmholtz equation for the fluence rate as a function of depth, z:

d2dz2φ0(z)μeff2φ0(z)=3μtrq0(z),
(4)

where

μeff=(μeff2+kx2+ky2)1/21δeff.
(5)

Here, a plane wave with both x and y modulation gives rise to a wave propagating with a scalar attenuation coefficient μeff. Although spatial anisotropy may exist in real tissues, we focus on the characteristics of a 1-D projection to understand simple scalar photon density wave attenuation in multiply scattering media. Consequently, the subsequent discussion considers a single nonzero spatial frequency along the x dimension only, k=kx, with constant illumination along y (ky=0).

At zero spatial frequency (k=0), the effective penetration depth δeffis equivalent to that of a planar (constant) illumination source, δeff=(1/μeff). In general, however, μeff(and δeff) 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, φ0(z), 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 μeff with our new μeffterm.

As in the derivation of Svaasand et al.5 for planar photon density wave reflectance, we model an extended source as

q0(z)=P0μsexp(μtrz),
(6)

where P0 is the incident optical power. Conceptually, this represents a spatially distributed but angularly isotropic source introduced via scattering of a collimated, forward-directed beam. The solution for the resulting fluence rate is

φ0(z)=3P0aμeff2/μtr21exp(μtrz)+Cexp(μeffz),

a=μsμtr,
(7)

where a′ is the reduced albedo, and C is a constant determined by the choice of a boundary condition. Using the partial current boundary condition,15 the flux, j, is set proportional to the fluence at the interface z=0:

j|z0+φ|z0+3μtr=Aφ|z0+,

A=1Reff2(1+Reff);Reff0.0636n+0.668+0.710n1.440n2.
(8)

Here, A is the proportionality constant, and Reff is the effective reflection coefficient. C then becomes

C=P03a(1+3A)(μeff2/μtr21)(μeff/μtr+3A,
(9)

yielding the diffuse reflectance, Rd(k):

Rd(k)=j|z0+P0=3Aa(μeff/μtr+1)(μeff/μtr+3A).
(10)

While the formulation shown is for a pure 1-D sinusoidal illumination pattern, an arbitrary illumination function can be modeled through linear superposition of sinusoids in both the x and y directions.

For a given set of optical properties, the function Rd(k) 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 Rd in the SFD is an inverse polynomial function of a single, positive-valued ratio, (μeff/μtr), which fully describes the low-pass spatial filtering properties of homogeneous turbid samples within a steady-state diffusion context:

μeffμtr=(μeff2+k2μtr2)1/2=(3μaμtr+k2μtr2)1/2={3(1a)1/2;k=0kμtr;kμeff.
(11)

The low- and high-frequency regimes have differential sensitivity to absorption and scattering properties, respectively. At k=0, Rd is a function only of the reduced albedo. In the additional limit of zero absorption, (μeff/μtr)0and Rd→1, implying that all incident photons are reflected back out of the turbid medium. At low spatial frequencies (k≪μeff), the absorption has a maximal effect on the reflectance. Approaching the high-frequency regime (k≫μeff), the denominator μtr (μsin the diffusion limit), becomes the only source of optical contrast. Both limits involve a ratio with respect to the transport coefficient, highlighting the natural length scale of light transport, l*=(1/μtr), the transport mean free path. In fact, dimensionless forms of the preceding fluence and reflectance solutions can be written in units of the transport spatial frequency, μtr, by the two substitutions μ̂a=(μa·l*), k̂=(k·l*), and =(z/l*), where μ^eff=(3μ^a+k^2)1/2and â′ = (1 − μ̂a).

The diffusion approximation to the radiative transport equation is valid when

μsμa,
(12)

and due to the anisotropic nature of light scattering, has been observed to be accurate approximately when

ρ ⪢ l,
(13)

where ρ 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 ρ is approximately in the range of 3l* to 4l* (Refs. 17 and 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), μtr=fx,tr=(ktr/2π). If we relate the inverse of ρ as a metric of spatial frequency, then Eq. (13) can be restated as

fxμtr1l*.
(14)

We therefore expect the maximum spatial frequency accurately predicted by diffusion to be in the range of 1/(3l*) to 1/(4l*), or 0.25μtr to 0.33μtr. Both albedo and source-distance requirements of the diffusion approximation limit the ratio (μeff/μ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 equation1920 and Monte Carlo simulation.2122 For this paper, we have used “White” Monte Carlo (WMC) simulations2324 of a collimated point-source illumination to generate predictions of the spatially resolved, steady-state diffuse reflectance, Rd(ρ), for a given set of optical properties μa, μs, and g. This spatial point-spread function, provides an impulse response, and spatial frequency domain predictions of diffuse reflectance, Rd(k), are found by Fourier transformation of Rd(ρ). For a radially symmetric function such as Rd(ρ), the 2-D Fourier transform in the x-y plane reduces to a 1-D Hankel transform of order zero:

Rd(k) = 2π∫ρJ0(kρ)Rd(ρ)dρ,
(15)

where J0(kr) is the zeroth-order Bessel function of the first kind. In our simulations, we bin ρ in n finite intervals Δρi. We can then calculate Rd(k) as

Rd(k)=2πi=1nρiJ0(kρi)Rd(ρi)Δρi.
(16)

In this paper, we generated the WMC data using 10 photons, a detector numerical aperture of 0.22. For all simulations, the index of refraction n and anisotropy factor g were set to 1.33 and 0.71, respectively, for direct comparison with our Liposyn phantom experiments. All radial bins had a spacing of Δρ=0.09 mm, making the maximum spatial frequency greater than 5 mm.

3.3 Simulations

Diffuse reflectance versus spatial frequency (mm), 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 (μs/μa)=100ratio (constant a′=0.99). Observe that as l* increases (or as μ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 Rd(fx) varies directly with μtr, or inversely with l*. This scaling of μtr and fx is directly evident in the (μeff/μtr)ratio of Eq. (11) (diffusion approximation), and thus all five diffusion curves will coincide perfectly if plotted versus normalized spatial frequency, (fx/fx,tr)=(fxtr)=(fx·l*). This behavior is also retained in our Monte Carlo predictions to a high degree of accuracy. For instance, when plotted versus (fx/fx,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, (fx/fx,tr). Conveniently, μtr=1 mm (l*=1 mm) is a good approximate transport coefficient for many biological tissues, so for the high-albedo curves, fx,tr can be interpreted as ~1 mm spatial frequency (1-mm-spaced sinusoids).

An external file that holds a picture, illustration, etc.
Object name is nihms165228f4.jpg

(a) Diffuse reflectance versus spatial frequency (mm) for varying values of l*, using both Monte Carlo simulations (solid lines) and the diffusion approximation (dashed lines), demonstrating accuracy of the diffusion approximation within 12% for fx≤(2l*)=0.5μtr. (b) Diffuse reflectance versus normalized spatial frequency (fx/fx,tr=fx·l*) for varying albedo (μs/μaratio), using both Monte Carlo simulations (solid lines) and the diffusion approximation (dashed lines), demonstrating degrading accuracy of diffusion with decreasing albedo. (Color online only.)

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 available2527 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 μs/μa=100) a diffusion error of ±12% for fx≤(2l*)=0.5μtr.

In Fig. 4(b), we further examine our transport and diffusion models when the albedo is decreased, ranging the (μs/μ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, fx,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 (fx/fx,tr≈1), while an absolute offset between curves remains in the transport prediction.

For low frequencies (below ~0.5·μtr), the diffusion error remains less than ±16% at (μs/μa)=30(a′=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 ±25% of the reference phantom (μs/μ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 (μs/μa)=10for 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:

S=S02[1+M0cos(2πfxx+α)],
(17)

where S0, M0, fx, and a are the illumination source intensity, modulation depth, spatial frequency, and spatial phase, respectively. In this simple case, the pattern is constant in the orthogonal y direction. In reflection mode, the diffusely reflected intensity, I, is a sum of AC and DC components:

I = IAC + IDC,
(18)

where the measured AC component of the reflected intensity, IAC, can be modeled as:

IAC = MAC(x, fx) · cos(2πfxx + α).
(19)

Here, MAC(x, fx) represents the amplitude envelope of the reflected photon density standing wave at frequency fx. Note first that an MAC can be a function of position, x, as shown in Fig. 5 (top). Additionally, multiple MAC(x, fx) curves can be sampled in parallel at each y pixel row using a 2-D camera, allowing spatial sampling of millions of reflectance values simultaneously.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f5.jpg

(a) Schematic of modulated reflectance and (b) demodulated AC and DC amplitudes.

A host of signal processing schemes can be used to obtain MAC(x, fx). Here, we employ a simple time-domain amplitude demodulation method,828 illuminating a sinusoid pattern three times at the same spatial frequency, with phase offsets a=0, 2/3π, and 4/3π radians. MAC(x, fx) can then be calculated algebraically at each spatial location, xi, by:

MAC(xi,fx)=21/23{[I1(xi)I2(xi)]2+[I2(xi)I3(xi)]2+[I3(xi)I1(xi)]2}1/2,
(20)

where I1, I2, and I3 represent the IAC image values at each location with shifted spatial phases. This differencing approach is convenient, as (1) it automatically removes features common to all three images, including the average image noise and digitization offset, and (2) it does not require knowledge of the spatial frequency, removing potential spatial calibration errors. Built in to this operation is an automatic subtraction of any constant ambient light present in each acquired image. The spatially varying DC amplitude, MDC(x), can be calculated as earlier with fx=0, or at any frequency of illumination using:

MDC(xi) = ⅓[I1(xi) + I2(xi) + I3(xi)].
(21)

In Fig. 5, we show a schematic of a spatially varying modulated reflectance (top) and its demodulated AC and DC amplitude (bottom) components.

In the frequency domain, a measurement MAC(fx) is the product of (1) the source intensity, I0; (2) the MTF of the illumination and imaging optical system, MTFsystem; and (3) the true turbid system MTF, Rd:

MAC(xi, fx) = I0 · MTFsystem(xi, fx) · Rd(xi, fx).
(22)

Therefore, we can simultaneously calibrate for the absolute intensity of the source and the MTF of the imaging system by performing a reference measurement, MAC,ref(x, fx), on a turbid phantom of known optical properties. Using a model prediction for the phantom diffuse reflectance, Rd,ref,pred(fx), we can write the diffuse reflectance at each spatial location as:

Rd(xi,fx)=MAC(xi,fx)MAC,ref(xi,fx)·Rd,ref,pred(fx).
(23)

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): μa and μs. 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 Rd(DC) and Rd(AC) correspond to diffuse reflectance measurements at zero and nonzero spatial frequencies f1 and f2, 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 Rd(0 mm)=0.55 and Rd(0.5 mm)=0.06, then μa≈0.03 and μs1.4mm1, 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 mm and 0.5 mm (DC and 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 AC and DC measurements can be easily obtained with only three phase projections of a single illumination frequency [through Eqs. (20) and (21)], allowing rapid, high-resolution imaging of absorption and scattering contrast.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f6.jpg

Two-frequency (DC versus AC) lookup table for rapid calculation of optical properties, generated from diffusion model forward predictions. Gray contours indicate constant absorption, black contours indicate constant reduced scattering. Dotted lines demonstrate the lookup method, translating DC and AC values into μa and μsparameters.

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), and (10), respectively. Data are processed separately for each pixel, generating spatial maps of absorption and scattering optical properties.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f7.jpg

Flow chart of MI data processing. Intensity data at each frequency (three phase images per frequency) are demodulated, calibrated, and fit using Eqs. (20), (23), and (10), respectively. Data are processed separately for each pixel, generating spatial maps of optical properties. Images are plotted within three standard deviations of the individual image mean in order to make them visually comparable. Notice the differential contrast in diffuse reflectance (Rd) versus spatial frequency (fx), the basis for separation absorption and scattering.

Compared to other spatially resolved methods, MI acquires coincident, axial “projection” measurements of optical contrast to quantify the optical properties at each xy 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 mm) 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).

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.3461215 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

2φμeff2φ=3μtrq,
(1)

where φ is the fluence rate, q is the source, μtr=(μa+μs)is the transport coefficient, μeff=(3μaμr), μa is the absorption coefficient, μs=μs(1g)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:

q = q0(z)cos(kxx + α)cos(kyy + β),
(2)

with spatial frequencies (or repetencies) fx=(kx/2π) and fy=(ky/2π), and spatial phases α and β, extending infinitely in the tangential spatial dimensions, x and y, with some arbitrary dependence on depth, z.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f3.jpg

Schematic of modulated illumination source (in the x direction only) and the resulting modulated internal fluence rate with the same frequency and phase.

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 shift16 for normally incident light onto a homogeneous medium.)

φ = φ0(z)cos(kxx + α)cos(kyx + β).
(3)

Insertion of Eqs. (2) and (3) into Eq. (1) yields a 1-D second-order Helmholtz equation for the fluence rate as a function of depth, z:

d2dz2φ0(z)μeff2φ0(z)=3μtrq0(z),
(4)

where

μeff=(μeff2+kx2+ky2)1/21δeff.
(5)

Here, a plane wave with both x and y modulation gives rise to a wave propagating with a scalar attenuation coefficient μeff. Although spatial anisotropy may exist in real tissues, we focus on the characteristics of a 1-D projection to understand simple scalar photon density wave attenuation in multiply scattering media. Consequently, the subsequent discussion considers a single nonzero spatial frequency along the x dimension only, k=kx, with constant illumination along y (ky=0).

At zero spatial frequency (k=0), the effective penetration depth δeffis equivalent to that of a planar (constant) illumination source, δeff=(1/μeff). In general, however, μeff(and δeff) 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, φ0(z), 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 μeff with our new μeffterm.

As in the derivation of Svaasand et al.5 for planar photon density wave reflectance, we model an extended source as

q0(z)=P0μsexp(μtrz),
(6)

where P0 is the incident optical power. Conceptually, this represents a spatially distributed but angularly isotropic source introduced via scattering of a collimated, forward-directed beam. The solution for the resulting fluence rate is

φ0(z)=3P0aμeff2/μtr21exp(μtrz)+Cexp(μeffz),

a=μsμtr,
(7)

where a′ is the reduced albedo, and C is a constant determined by the choice of a boundary condition. Using the partial current boundary condition,15 the flux, j, is set proportional to the fluence at the interface z=0:

j|z0+φ|z0+3μtr=Aφ|z0+,

A=1Reff2(1+Reff);Reff0.0636n+0.668+0.710n1.440n2.
(8)

Here, A is the proportionality constant, and Reff is the effective reflection coefficient. C then becomes

C=P03a(1+3A)(μeff2/μtr21)(μeff/μtr+3A,
(9)

yielding the diffuse reflectance, Rd(k):

Rd(k)=j|z0+P0=3Aa(μeff/μtr+1)(μeff/μtr+3A).
(10)

While the formulation shown is for a pure 1-D sinusoidal illumination pattern, an arbitrary illumination function can be modeled through linear superposition of sinusoids in both the x and y directions.

For a given set of optical properties, the function Rd(k) 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 Rd in the SFD is an inverse polynomial function of a single, positive-valued ratio, (μeff/μtr), which fully describes the low-pass spatial filtering properties of homogeneous turbid samples within a steady-state diffusion context:

μeffμtr=(μeff2+k2μtr2)1/2=(3μaμtr+k2μtr2)1/2={3(1a)1/2;k=0kμtr;kμeff.
(11)

The low- and high-frequency regimes have differential sensitivity to absorption and scattering properties, respectively. At k=0, Rd is a function only of the reduced albedo. In the additional limit of zero absorption, (μeff/μtr)0and Rd→1, implying that all incident photons are reflected back out of the turbid medium. At low spatial frequencies (k≪μeff), the absorption has a maximal effect on the reflectance. Approaching the high-frequency regime (k≫μeff), the denominator μtr (μsin the diffusion limit), becomes the only source of optical contrast. Both limits involve a ratio with respect to the transport coefficient, highlighting the natural length scale of light transport, l*=(1/μtr), the transport mean free path. In fact, dimensionless forms of the preceding fluence and reflectance solutions can be written in units of the transport spatial frequency, μtr, by the two substitutions μ̂a=(μa·l*), k̂=(k·l*), and =(z/l*), where μ^eff=(3μ^a+k^2)1/2and â′ = (1 − μ̂a).

The diffusion approximation to the radiative transport equation is valid when

μsμa,
(12)

and due to the anisotropic nature of light scattering, has been observed to be accurate approximately when

ρ ⪢ l,
(13)

where ρ 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 ρ is approximately in the range of 3l* to 4l* (Refs. 17 and 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), μtr=fx,tr=(ktr/2π). If we relate the inverse of ρ as a metric of spatial frequency, then Eq. (13) can be restated as

fxμtr1l*.
(14)

We therefore expect the maximum spatial frequency accurately predicted by diffusion to be in the range of 1/(3l*) to 1/(4l*), or 0.25μtr to 0.33μtr. Both albedo and source-distance requirements of the diffusion approximation limit the ratio (μeff/μ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 equation1920 and Monte Carlo simulation.2122 For this paper, we have used “White” Monte Carlo (WMC) simulations2324 of a collimated point-source illumination to generate predictions of the spatially resolved, steady-state diffuse reflectance, Rd(ρ), for a given set of optical properties μa, μs, and g. This spatial point-spread function, provides an impulse response, and spatial frequency domain predictions of diffuse reflectance, Rd(k), are found by Fourier transformation of Rd(ρ). For a radially symmetric function such as Rd(ρ), the 2-D Fourier transform in the x-y plane reduces to a 1-D Hankel transform of order zero:

Rd(k) = 2π∫ρJ0(kρ)Rd(ρ)dρ,
(15)

where J0(kr) is the zeroth-order Bessel function of the first kind. In our simulations, we bin ρ in n finite intervals Δρi. We can then calculate Rd(k) as

Rd(k)=2πi=1nρiJ0(kρi)Rd(ρi)Δρi.
(16)

In this paper, we generated the WMC data using 10 photons, a detector numerical aperture of 0.22. For all simulations, the index of refraction n and anisotropy factor g were set to 1.33 and 0.71, respectively, for direct comparison with our Liposyn phantom experiments. All radial bins had a spacing of Δρ=0.09 mm, making the maximum spatial frequency greater than 5 mm.

3.3 Simulations

Diffuse reflectance versus spatial frequency (mm), 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 (μs/μa)=100ratio (constant a′=0.99). Observe that as l* increases (or as μ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 Rd(fx) varies directly with μtr, or inversely with l*. This scaling of μtr and fx is directly evident in the (μeff/μtr)ratio of Eq. (11) (diffusion approximation), and thus all five diffusion curves will coincide perfectly if plotted versus normalized spatial frequency, (fx/fx,tr)=(fxtr)=(fx·l*). This behavior is also retained in our Monte Carlo predictions to a high degree of accuracy. For instance, when plotted versus (fx/fx,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, (fx/fx,tr). Conveniently, μtr=1 mm (l*=1 mm) is a good approximate transport coefficient for many biological tissues, so for the high-albedo curves, fx,tr can be interpreted as ~1 mm spatial frequency (1-mm-spaced sinusoids).

An external file that holds a picture, illustration, etc.
Object name is nihms165228f4.jpg

(a) Diffuse reflectance versus spatial frequency (mm) for varying values of l*, using both Monte Carlo simulations (solid lines) and the diffusion approximation (dashed lines), demonstrating accuracy of the diffusion approximation within 12% for fx≤(2l*)=0.5μtr. (b) Diffuse reflectance versus normalized spatial frequency (fx/fx,tr=fx·l*) for varying albedo (μs/μaratio), using both Monte Carlo simulations (solid lines) and the diffusion approximation (dashed lines), demonstrating degrading accuracy of diffusion with decreasing albedo. (Color online only.)

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 available2527 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 μs/μa=100) a diffusion error of ±12% for fx≤(2l*)=0.5μtr.

In Fig. 4(b), we further examine our transport and diffusion models when the albedo is decreased, ranging the (μs/μ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, fx,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 (fx/fx,tr≈1), while an absolute offset between curves remains in the transport prediction.

For low frequencies (below ~0.5·μtr), the diffusion error remains less than ±16% at (μs/μa)=30(a′=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 ±25% of the reference phantom (μs/μ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 (μs/μa)=10for 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:

S=S02[1+M0cos(2πfxx+α)],
(17)

where S0, M0, fx, and a are the illumination source intensity, modulation depth, spatial frequency, and spatial phase, respectively. In this simple case, the pattern is constant in the orthogonal y direction. In reflection mode, the diffusely reflected intensity, I, is a sum of AC and DC components:

I = IAC + IDC,
(18)

where the measured AC component of the reflected intensity, IAC, can be modeled as:

IAC = MAC(x, fx) · cos(2πfxx + α).
(19)

Here, MAC(x, fx) represents the amplitude envelope of the reflected photon density standing wave at frequency fx. Note first that an MAC can be a function of position, x, as shown in Fig. 5 (top). Additionally, multiple MAC(x, fx) curves can be sampled in parallel at each y pixel row using a 2-D camera, allowing spatial sampling of millions of reflectance values simultaneously.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f5.jpg

(a) Schematic of modulated reflectance and (b) demodulated AC and DC amplitudes.

A host of signal processing schemes can be used to obtain MAC(x, fx). Here, we employ a simple time-domain amplitude demodulation method,828 illuminating a sinusoid pattern three times at the same spatial frequency, with phase offsets a=0, 2/3π, and 4/3π radians. MAC(x, fx) can then be calculated algebraically at each spatial location, xi, by:

MAC(xi,fx)=21/23{[I1(xi)I2(xi)]2+[I2(xi)I3(xi)]2+[I3(xi)I1(xi)]2}1/2,
(20)

where I1, I2, and I3 represent the IAC image values at each location with shifted spatial phases. This differencing approach is convenient, as (1) it automatically removes features common to all three images, including the average image noise and digitization offset, and (2) it does not require knowledge of the spatial frequency, removing potential spatial calibration errors. Built in to this operation is an automatic subtraction of any constant ambient light present in each acquired image. The spatially varying DC amplitude, MDC(x), can be calculated as earlier with fx=0, or at any frequency of illumination using:

MDC(xi) = ⅓[I1(xi) + I2(xi) + I3(xi)].
(21)

In Fig. 5, we show a schematic of a spatially varying modulated reflectance (top) and its demodulated AC and DC amplitude (bottom) components.

In the frequency domain, a measurement MAC(fx) is the product of (1) the source intensity, I0; (2) the MTF of the illumination and imaging optical system, MTFsystem; and (3) the true turbid system MTF, Rd:

MAC(xi, fx) = I0 · MTFsystem(xi, fx) · Rd(xi, fx).
(22)

Therefore, we can simultaneously calibrate for the absolute intensity of the source and the MTF of the imaging system by performing a reference measurement, MAC,ref(x, fx), on a turbid phantom of known optical properties. Using a model prediction for the phantom diffuse reflectance, Rd,ref,pred(fx), we can write the diffuse reflectance at each spatial location as:

Rd(xi,fx)=MAC(xi,fx)MAC,ref(xi,fx)·Rd,ref,pred(fx).
(23)

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): μa and μs. 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 Rd(DC) and Rd(AC) correspond to diffuse reflectance measurements at zero and nonzero spatial frequencies f1 and f2, 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 Rd(0 mm)=0.55 and Rd(0.5 mm)=0.06, then μa≈0.03 and μs1.4mm1, 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 mm and 0.5 mm (DC and 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 AC and DC measurements can be easily obtained with only three phase projections of a single illumination frequency [through Eqs. (20) and (21)], allowing rapid, high-resolution imaging of absorption and scattering contrast.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f6.jpg

Two-frequency (DC versus AC) lookup table for rapid calculation of optical properties, generated from diffusion model forward predictions. Gray contours indicate constant absorption, black contours indicate constant reduced scattering. Dotted lines demonstrate the lookup method, translating DC and AC values into μa and μsparameters.

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), and (10), respectively. Data are processed separately for each pixel, generating spatial maps of absorption and scattering optical properties.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f7.jpg

Flow chart of MI data processing. Intensity data at each frequency (three phase images per frequency) are demodulated, calibrated, and fit using Eqs. (20), (23), and (10), respectively. Data are processed separately for each pixel, generating spatial maps of optical properties. Images are plotted within three standard deviations of the individual image mean in order to make them visually comparable. Notice the differential contrast in diffuse reflectance (Rd) versus spatial frequency (fx), the basis for separation absorption and scattering.

Compared to other spatially resolved methods, MI acquires coincident, axial “projection” measurements of optical contrast to quantify the optical properties at each xy 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 mm) 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, μa, over two orders of magnitude (logarithmically spaced between 0.002 mm≤μa≤0.12 mm), with a constant scattering coefficient constant at μs=0.97mm1. In the second set, we linearly varied μs(0.32mm1μs1.8mm1), while holding the absorption coefficient constant at μa=0.0046 mm. These values were calculated based on infinite-geometry, multifrequency (50 to 500 MHz), multidistance (10 to 25 mm) frequency-domain photon migration measurements15 of one of the Liposyn/nigrosin phantoms.

MI measurements were performed on each sample. Thirty spatial frequencies of illumination were chosen between 0 mm and 0.13 mm, 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 nm with an approximate 75×75 mm illumination area, a 50×50 mm camera field of view, and an integration time of 100 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×20 pixel binning was performed on each image to speed computation, resulting in low-resolution, 15×15 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, Rd,ref,pred(fx). Then, for each spatial frequency and each spatial location, we use Eq. (23) to calculate Rd(fx) 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 μa and μsusing both least-squares minimization by a simplex search algorithm (in “fminsearch” MATLAB) and via the two-frequency lookup table approach using the lowest (0 mm) and highest (0.13 mm) 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×48 mm field of view. Four evenly spaced spatial frequencies between 0 and 0.15 mm were collected and analyzed. The imaging system was identical to that described earlier, except for the inclusion of a 640±10 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×40 mm region of the volar forearm. Measurements were performed at a wavelength of 800±10 nm, near the hemoglobin isosbestic point of 805 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 mm were acquired every 4 s for a period of 13 min. After 2.5 min of baseline acquisition, an arm cuff was pressurized to 100 mm Hg for 6.5 min and subsequently released at minute 9.

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, μa, over two orders of magnitude (logarithmically spaced between 0.002 mm≤μa≤0.12 mm), with a constant scattering coefficient constant at μs=0.97mm1. In the second set, we linearly varied μs(0.32mm1μs1.8mm1), while holding the absorption coefficient constant at μa=0.0046 mm. These values were calculated based on infinite-geometry, multifrequency (50 to 500 MHz), multidistance (10 to 25 mm) frequency-domain photon migration measurements15 of one of the Liposyn/nigrosin phantoms.

MI measurements were performed on each sample. Thirty spatial frequencies of illumination were chosen between 0 mm and 0.13 mm, 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 nm with an approximate 75×75 mm illumination area, a 50×50 mm camera field of view, and an integration time of 100 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×20 pixel binning was performed on each image to speed computation, resulting in low-resolution, 15×15 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, Rd,ref,pred(fx). Then, for each spatial frequency and each spatial location, we use Eq. (23) to calculate Rd(fx) 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 μa and μsusing both least-squares minimization by a simplex search algorithm (in “fminsearch” MATLAB) and via the two-frequency lookup table approach using the lowest (0 mm) and highest (0.13 mm) 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×48 mm field of view. Four evenly spaced spatial frequencies between 0 and 0.15 mm were collected and analyzed. The imaging system was identical to that described earlier, except for the inclusion of a 640±10 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×40 mm region of the volar forearm. Measurements were performed at a wavelength of 800±10 nm, near the hemoglobin isosbestic point of 805 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 mm were acquired every 4 s for a period of 13 min. After 2.5 min of baseline acquisition, an arm cuff was pressurized to 100 mm Hg for 6.5 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. 8(a) and 8(b), respectively. In solid black lines, we show the corresponding fits using the diffusion-based reflectance model [Eq. (10)].

An external file that holds a picture, illustration, etc.
Object name is nihms165228f8.jpg

Liposyn experiment data (black circles) fit to the SFD diffusion model of Eq. (10) (gray lines). (a) As absorption increases, reflectance at low spatial frequencies decreases dramatically, while high-frequency data shows little sensitivity to absorption contrast. (b) Increasing scattering gives rise to an increase in reflectance, with scattering contrast observed at all frequencies. All results show excellent visual agreement between measurement data and diffusion fits.

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 res-caling 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 (μs=0.32mm1,l*3mm). This seems consistent with the l* plots of Fig. 4, where we would expect model breakdown at or before fx=1/(2l*), or 0.16 mm.

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×15 pixel optical property maps (not shown). Thin black lines are drawn to indicate the expected values in each experiment.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f9.jpg

(a) Recovered versus expected optical properties for (a) absorption variation and (b) scattering variation experiments, using both multifrequency (black lines) and 2-frequency lookup table (gray lines) fitting methods. (c) and (d) The corresponding recovered “constant” values for reduced scattering and absorption are shown respectively. Error bars indicate the corresponding standard deviations of 15×15 pixel optical property maps (generally <1%), indicating both high precision and spatial uniformity over the field of view. Thin black lines indicate the expected values in each experiment. Circles indicate the phantom used for calibration.

In the absorption variation experiment, recovered versus expected absorption shows excellent linearity over two orders of magnitude, ranging from high to low albedo (μs/μa=500toμs/μ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×15 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 experimentScattering experiment

Average error
(%)
30-frequency fit2-frequency fitAverage error
(%)
30-frequency fit2-frequency fit
μa error4.744.85μa error7.5111.4
μserror2.982.29μserror3.0510.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 (0.002 mm). 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 fx=0) correspond to spatial periods larger than the projector’s illumination area (fx<0.013 mm). 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 Rd, 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 (Rd) versus spatial frequency (fx) 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. 10(a), we show the calibrated diffuse reflectance at fx=0 mm, and in Fig. 10(d), the average multifrequency diffuse reflectance data (black squares) and fit (gray line). In Figs. 10(b) and 10(c), we show spatial maps of the absorption and reduced scattering data, respectively, and in Figs. 10(e) and 10(f), 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.3132 Based on a diffusion-based sensitivity analysis, we predict for this experiment 1/e sampling depths ranging from 2 mm to 3.3 mm for low (0 mm) and high (0.15 mm) spatial frequencies, respectively.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f10.jpg

(a) Planar diffuse reflectance, (b) optical absorption, and (c) reduced scattering maps. (d) Average multifrequency diffuse reflectance data and least-squares fit to Eq. (10) and pixel histograms of (e) absorption and (f) scattering maps, 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 (0.46 mm) compared to the background (0.23 mm); in the same region, optical scattering showed little contrast, with ~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 nm. Figure 11(a) 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. 11(b). 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 min, absorption decreases to the baseline over the course of approximately 2 min. Maximal increases in absorption were observed to be approximately 0.012 mm globally and 0.017 mm 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.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f11.jpg

Venous occlusion data of the volar forearm measurement at 800 nm. (a) Diffuse reflectance map (left) and optical absorption map (right) measured at the baseline. Dotted lines in the reflectance map indicate the regions of interest for time-course analysis. (b) Region-wise average changes in optical absorption (top) and reduced scattering (bottom) are shown for the whole image field (gray lines) and a region absent of any obvious large vessels (black lines). 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. (Color online only.)

5.4 Sensitivity, Contrast, and Resolution

Detecting changes in μa and μsrequires 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 (ΔRd) versus normalized spatial frequency resulting from a 1% change in absorption or scattering, for sample (μs/μ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.)

An external file that holds a picture, illustration, etc.
Object name is nihms165228f12.jpg

Sensitivity of reflectance to 1% change in optical properties for absorption (solid lines) and reduced scattering (dashed lines) optical properties versus normalized spatial frequency. Results indicate that scattering contrast is retained at higher frequencies than absorption contrast. Dotted lines indicate detection thresholds for 12- and 14-bit dynamic range measurements.

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 Rd. We have added dotted lines to Fig. 12 to show the approximate corresponding camera detection limits with 12- and 14-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 (μs/μa)=100at fx/fx,tr=0.1, a 12-bit measurement can resolve (in intensity) a 1% change in scattering but not absorption, while a 14-bit camera can resolve both. For (μs/μa)=3, where optical properties are closer in magnitude, this 12-bit absorption detection criterion occurs at a higher frequency of fx/fx,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-s integration times. Both intensity stability (~10% decrease over 4h) 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 (μs/μ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. 13(a), 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. 13(a), left; μa,left=0.01 mm, μa,right=0.03 mm; [μs,both=1.0mm1] or scattering [Fig. 13(a), right; μa,both=0.02 mm; μs,left=0.5mm1, μs,right=1.5mm1] were measured at 660 nm both directly [Fig. 13(a), top] and through a 2-mm homogeneous layer [Fig. 13(a), bottom] with optical properties μa=0.01 mm, μs=0.5mm1. Both were calibrated by a homogeneous phantom with optical properties μa=0.02 mm, μs=1.0mm1. Nine spatial frequencies between 0 mm and 0.11 mm were measured and used to calculate optical property maps using least-squares regression to our diffusion reflectance model.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f13.jpg

Preliminary contrast-resolution study. (a) Gelatin phantoms with a 3× step in either absorption (left) or scattering (right) were measured directly (top) and through a 2-mm homogeneous layer (bottom). (b) Edge-response profiles for absorption (top) and scattering (bottom) media reveal loss of both resolution and quantitative contrast through the homogeneous layer (gray lines) compared to that at the surface (black lines).

Recovered optical property spatial profiles were averaged over the vertical direction, shown in Fig. 13(b) 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 mm), 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 mm and 0.05 mm, respectively, for surface perturbations, and 0.5 mm and 0.25 mm, respectively, for perturbations 2 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-y resolution to scale with the number of spatial sampling points (# of detectors).33

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. 8(a) and 8(b), respectively. In solid black lines, we show the corresponding fits using the diffusion-based reflectance model [Eq. (10)].

An external file that holds a picture, illustration, etc.
Object name is nihms165228f8.jpg

Liposyn experiment data (black circles) fit to the SFD diffusion model of Eq. (10) (gray lines). (a) As absorption increases, reflectance at low spatial frequencies decreases dramatically, while high-frequency data shows little sensitivity to absorption contrast. (b) Increasing scattering gives rise to an increase in reflectance, with scattering contrast observed at all frequencies. All results show excellent visual agreement between measurement data and diffusion fits.

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 res-caling 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 (μs=0.32mm1,l*3mm). This seems consistent with the l* plots of Fig. 4, where we would expect model breakdown at or before fx=1/(2l*), or 0.16 mm.

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×15 pixel optical property maps (not shown). Thin black lines are drawn to indicate the expected values in each experiment.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f9.jpg

(a) Recovered versus expected optical properties for (a) absorption variation and (b) scattering variation experiments, using both multifrequency (black lines) and 2-frequency lookup table (gray lines) fitting methods. (c) and (d) The corresponding recovered “constant” values for reduced scattering and absorption are shown respectively. Error bars indicate the corresponding standard deviations of 15×15 pixel optical property maps (generally <1%), indicating both high precision and spatial uniformity over the field of view. Thin black lines indicate the expected values in each experiment. Circles indicate the phantom used for calibration.

In the absorption variation experiment, recovered versus expected absorption shows excellent linearity over two orders of magnitude, ranging from high to low albedo (μs/μa=500toμs/μ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×15 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 experimentScattering experiment

Average error
(%)
30-frequency fit2-frequency fitAverage error
(%)
30-frequency fit2-frequency fit
μa error4.744.85μa error7.5111.4
μserror2.982.29μserror3.0510.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 (0.002 mm). 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 fx=0) correspond to spatial periods larger than the projector’s illumination area (fx<0.013 mm). 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 Rd, 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 (Rd) versus spatial frequency (fx) 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. 10(a), we show the calibrated diffuse reflectance at fx=0 mm, and in Fig. 10(d), the average multifrequency diffuse reflectance data (black squares) and fit (gray line). In Figs. 10(b) and 10(c), we show spatial maps of the absorption and reduced scattering data, respectively, and in Figs. 10(e) and 10(f), 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.3132 Based on a diffusion-based sensitivity analysis, we predict for this experiment 1/e sampling depths ranging from 2 mm to 3.3 mm for low (0 mm) and high (0.15 mm) spatial frequencies, respectively.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f10.jpg

(a) Planar diffuse reflectance, (b) optical absorption, and (c) reduced scattering maps. (d) Average multifrequency diffuse reflectance data and least-squares fit to Eq. (10) and pixel histograms of (e) absorption and (f) scattering maps, 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 (0.46 mm) compared to the background (0.23 mm); in the same region, optical scattering showed little contrast, with ~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 nm. Figure 11(a) 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. 11(b). 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 min, absorption decreases to the baseline over the course of approximately 2 min. Maximal increases in absorption were observed to be approximately 0.012 mm globally and 0.017 mm 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.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f11.jpg

Venous occlusion data of the volar forearm measurement at 800 nm. (a) Diffuse reflectance map (left) and optical absorption map (right) measured at the baseline. Dotted lines in the reflectance map indicate the regions of interest for time-course analysis. (b) Region-wise average changes in optical absorption (top) and reduced scattering (bottom) are shown for the whole image field (gray lines) and a region absent of any obvious large vessels (black lines). 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. (Color online only.)

5.4 Sensitivity, Contrast, and Resolution

Detecting changes in μa and μsrequires 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 (ΔRd) versus normalized spatial frequency resulting from a 1% change in absorption or scattering, for sample (μs/μ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.)

An external file that holds a picture, illustration, etc.
Object name is nihms165228f12.jpg

Sensitivity of reflectance to 1% change in optical properties for absorption (solid lines) and reduced scattering (dashed lines) optical properties versus normalized spatial frequency. Results indicate that scattering contrast is retained at higher frequencies than absorption contrast. Dotted lines indicate detection thresholds for 12- and 14-bit dynamic range measurements.

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 Rd. We have added dotted lines to Fig. 12 to show the approximate corresponding camera detection limits with 12- and 14-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 (μs/μa)=100at fx/fx,tr=0.1, a 12-bit measurement can resolve (in intensity) a 1% change in scattering but not absorption, while a 14-bit camera can resolve both. For (μs/μa)=3, where optical properties are closer in magnitude, this 12-bit absorption detection criterion occurs at a higher frequency of fx/fx,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-s integration times. Both intensity stability (~10% decrease over 4h) 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 (μs/μ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. 13(a), 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. 13(a), left; μa,left=0.01 mm, μa,right=0.03 mm; [μs,both=1.0mm1] or scattering [Fig. 13(a), right; μa,both=0.02 mm; μs,left=0.5mm1, μs,right=1.5mm1] were measured at 660 nm both directly [Fig. 13(a), top] and through a 2-mm homogeneous layer [Fig. 13(a), bottom] with optical properties μa=0.01 mm, μs=0.5mm1. Both were calibrated by a homogeneous phantom with optical properties μa=0.02 mm, μs=1.0mm1. Nine spatial frequencies between 0 mm and 0.11 mm were measured and used to calculate optical property maps using least-squares regression to our diffusion reflectance model.

An external file that holds a picture, illustration, etc.
Object name is nihms165228f13.jpg

Preliminary contrast-resolution study. (a) Gelatin phantoms with a 3× step in either absorption (left) or scattering (right) were measured directly (top) and through a 2-mm homogeneous layer (bottom). (b) Edge-response profiles for absorption (top) and scattering (bottom) media reveal loss of both resolution and quantitative contrast through the homogeneous layer (gray lines) compared to that at the surface (black lines).

Recovered optical property spatial profiles were averaged over the vertical direction, shown in Fig. 13(b) 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 mm), 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 mm and 0.05 mm, respectively, for surface perturbations, and 0.5 mm and 0.25 mm, respectively, for perturbations 2 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-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 (0.5 mm to 3 mm) and (μs/μ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.

David J. Cuccia, Modulated Imaging, Inc., 1002 Health Sciences Road, Irvine, California 92612;
Contributor Information.
Address all correspondence to: David J. Cuccia, Modulated Imaging, Inc., 1002 Health Sciences Rd., Irvine, CA 92612. Tel: 949-824-8367; Fax: 949-824-6969; moc.gnigamidetaludom@aiccuc.divad
Current address: Ircam, 1 Place Igor Stravinsky, Paris, 75001, France.

Abstract

We describe the development of a rapid, noncontact imaging method, modulated imaging (MI), for quantitative, wide-field characterization of optical absorption and scattering properties of turbid media. MI utilizes principles of frequency-domain sampling and model-based analysis of the spatial modulation transfer function (s-MTF). We present and compare analytic diffusion and probabilistic Monte Carlo models of diffuse reflectance in the spatial frequency domain. Next, we perform MI measurements on tissue-simulating phantoms exhibiting a wide range of l* values (0.5 mm to 3 mm) and (μ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.

Keywords: structural illumination, turbid media, optical properties, diffuse optical spectroscopy, wide-field imaging, spatial frequency domain
Abstract

Contributor Information

David J. Cuccia, Modulated Imaging, Inc., 1002 Health Sciences Road, Irvine, California 92612.

Frederic Bevilacqua, University of California, Irvine, Beckman Laser Institute, 1002 Health Sciences Road, Irvine, California 92612.

Anthony J. Durkin, University of California, Irvine, Beckman Laser Institute, 1002 Health Sciences Road, Irvine, California 92612.

Frederick R. Ayers, University of California, Irvine, Beckman Laser Institute, 1002 Health Sciences Road, Irvine, California 92612.

Bruce J. Tromberg, University of California, Irvine, Beckman Laser Institute, 1002 Health Sciences Road, Irvine, California 92612.

Contributor Information

References

  • 1. Welch AJ, Yoon G, van Gemert MJPractical models for light distribution in laser-irradiated tissue. Lasers Surg. Med. 1987;6(6):488–493.[PubMed][Google Scholar]
  • 2. Yoo K, Alfano RDetermination of the scattering and absorption lengths from the temporal profile of a backscattered pulse. Opt. Lett. 1990;15(5):276–278.[PubMed][Google Scholar]
  • 3. O’Leary MA, Boas DA, Chance B, Yodh AGRefraction of diffuse photon density waves. Phys. Rev. Lett. 1992;69(18):2658–2661.[PubMed][Google Scholar]
  • 4. Fishkin JB, Gratton EPropagation of photon-density waves in strongly scattering media containing an absorbing semi-infinite plane bounded by a straight edge. J. Opt. Soc. Am. A. 1993;10(1):127–140.[PubMed][Google Scholar]
  • 5. Svaasand LO, Spott T, Fishkin JB, Pham T, Tromberg BJ, Berns MWReflectance measurements of layered media with diffuse photon-density waves: a potential tool for evaluating deep burns and subcutaneous lesions. Phys. Med. Biol. 1999;44(3):801–813.[PubMed][Google Scholar]
  • 6. Bonner RF, Nossal N, Havlin S, Weiss GHModel for photon migration in turbid biological media. J. Opt. Soc. Am. A. 1987;4(3):423–432.[PubMed][Google Scholar]
  • 7. Farrell TJ, Patterson MS, Wilson BA diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo. Med. Phys. 1992;19(4):879–888.[PubMed][Google Scholar]
  • 8. Neil MAA, Juskaitis R, Wilson TMethod of obtaining optical sectioning by using structured light in a conventional microscope. Opt. Lett. 1997;22(24):1905–1907.[PubMed][Google Scholar]
  • 9. Dognitz N, Wagnieres GDetermination of tissue optical properties by steady-state spatial frequency-domain reflectometry. Lasers Med. Sci. 1998;13:55–65.[PubMed][Google Scholar]
  • 10. Cuccia DJ, Bevilacqua F, Durkin AJ, Tromberg BJModulated imaging: quantitative analysis and tomography of turbid media in the spatial-frequency domain. Opt. Lett. 2005;30(11):1354–1356.[PubMed][Google Scholar]
  • 11. Joshi A, Bangerth W, Sevick-Muraca ENoncontact fluorescence optical tomography with scanning patterned illumination. Opt. Express. 2006;14:6516–6534.[PubMed][Google Scholar]
  • 12. Patterson M, Chance B, Wilson BTime-resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties. Appl. Opt. 1989;28:2331–2336.[PubMed][Google Scholar]
  • 13. Boas DA, O’Leary MA, Chance B, Yodh AGScattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytic solution and applications. Proc. Natl. Acad. Sci. U.S.A. 1994;91(11):4887–4891.[Google Scholar]
  • 14. Li X, et al Near-field diffraction tomography with diffuse photon density waves. Phys. Rev. E. 2000;61(4 Pt B):4295–4309.[PubMed][Google Scholar]
  • 15. Haskell RC, Svaasand LO, Tsay TT, Feng TC, McAd-ams MS, Tromberg BJBoundary conditions for the diffusion equation in radiative transfer. J. Opt. Soc. Am. A. 1994;11(10):2727–2741.[PubMed][Google Scholar]
  • 16. Bassi A, Cuccia DJ, Durkin AJ, Tromberg BJSpatial shift of spatially modulated light projected on turbid media. J. Opt. Soc. Am. A. 2008;25(11):2833–2839.[Google Scholar]
  • 17. Yoo KM, Liu F, Alfano RRWhen does the diffusion approximation fail to describe photon transport in random media? Phys. Rev. Lett. 1990;64(22):2647–2650.[PubMed][Google Scholar]
  • 18. Fishkin JB, Fantini S, vandeVen MJ, Gratton EGigahertz photon density waves in a turbid medium: theory and experiments. Phys. Rev. E. 1996;53(3):2307–2319.[PubMed][Google Scholar]
  • 19. Ishimaru A Wave Propagation and Scattering in Random Media. New York: IEEE Press; 1978. [PubMed][Google Scholar]
  • 20. Kim ADTransport theory for light propagation in biological tissue. J. Opt. Soc. Am. A. 2004;21(5):820–827.[PubMed][Google Scholar]
  • 21. Jacques SL, Wang L. Monte Carlo modeling of light transport in tissues. In: Welch AJ, van Gemert MJC, editors. Optical-Thermal Response of Laser-Irradiated Tissue. New York: Plenum; 1995. pp. 73–99. [PubMed]
  • 22. Marquet P, Bevilacqua F, Depeursinge C, de Haller EB. Determination of reduced scattering and absorption coefficients by a single charge-coupled-device array measurement. Part I: Comparison between experiments and simulations. Opt. Eng. 1995;34:2055–2063.[PubMed]
  • 23. Kienle A, Patterson MADetermination of the optical properties of turbid media from a single Monte Carlo simulation. Phys. Med. Biol. 1996;41:2221–2227.[PubMed][Google Scholar]
  • 24. Swartling J, Pifferi A, Enejder AM, Andersson-Engels AAccelerated Monte Carlo models to simulate fluorescence spectra from layered tissues. J. Opt. Soc. Am. A. 2003;20(4):714–727.[PubMed][Google Scholar]
  • 25. Joseph JH, Wiscombe WJ, Weinman JADelta-Eddington approximation for radiative flux-transfer. J. Atmos. Sci. 1976;33(12):2452–2459.[PubMed][Google Scholar]
  • 26. Spott T, Svaasand LOCollimated light sources in the diffusion approximation. Appl. Opt. 2000;39(34):6453–6465.[PubMed][Google Scholar]
  • 27. Carp SA, Prahl SA, Venugopalan VRadiative transport in the delta-P1 approximation: accuracy of fluence rate and optical penetration depth predictions in turbid semi-infinite media. J. Biomed. Opt. 2004;9(3):632–647.[PubMed][Google Scholar]
  • 28. Carlson AB Communication Systems. New York: McGraw-Hill; 1988. [PubMed][Google Scholar]
  • 29. Gioux S, Mazhar A, Cuccia D, Durkin A, Tromberg BJ, Frangioni JV Frontiers in Optics 2008/Laser Science (LS) XXIV. Rochester, NY: OSA Technical Digest (CD), paper FWP2, Optical Society of America; 2008. Fluorescence lifetime and spatially modulated light for image-guided surgery (invited) [PubMed][Google Scholar]
  • 30. Pham T, Coquoz O, Fishkin JB, Anderson E, Tromberg BJA broad bandwidth frequency domain instrument for quantitative tissue optical spectroscopy. Rev. Sci. Instrum. 2000;7:2500–2513.[PubMed][Google Scholar]
  • 31. Vidal BCThe part played by the mucopolysaccharides in the form birefringence of the collogen. Protoplasma. 1964;59:472–479.[PubMed][Google Scholar]
  • 32. Roth S, Freund IOptical second-harmonic scattering in rat-tail tendon. Biopolymers. 1981;20(6):1271–1290.[PubMed][Google Scholar]
  • 33. Markel VA, Schotland JCSymmetries, inversion formulas, and image reconstruction for optical tomography. Phys. Rev. E. 2004;70(5 Pt 2):056616.[PubMed][Google Scholar]
Collaboration tool especially designed for Life Science professionals.Drag-and-drop any entity to your messages.