3.1.

Diffusion Approximation

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

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

Fig. 3

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

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

As in the derivation of Svaasand ⁵ for planar photon density wave reflectance, we model an extended source as

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

The diffusion approximation to the radiative transport equation is valid when

3.2.

Monte Carlo Simulation in the SFD

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

3.3.

Simulations

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

Fig. 4

(a) Diffuse reflectance versus spatial frequency $\left({\mathrm{mm}}^{-1}\right)$ 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 $f_x⩽{\left(2l^{*}\right)}^{-1}=0.5{\mu }_{\mathit{tr}}$ . (b) Diffuse reflectance versus normalized spatial frequency $(f_x∕f_{x,\mathit{tr}}=f_x\cdot l^{*})$ for varying albedo ( ${\mu }_s^{\prime }∕{\mu }_a$ ratio), 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 available^{25, 26, 27} and will be examined in future work. Defining diffusion error as the percent difference of the diffusion prediction from Monte Carlo, we observe (for the case of ${\mu }_s^{\prime }∕{\mu }_a=100$ ) a diffusion error of $\pm 12\%$ for $f_x⩽{\left(2l^{*}\right)}^{-1}=0.5{\mu }_{\mathit{tr}}$ .

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

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

3.4.

Illumination, Imaging, and Calibration

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

Fig. 5

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

A host of signal processing schemes can be used to obtain $M_{\mathit{AC}}(x,f_x)$ . Here, we employ a simple time-domain amplitude demodulation method,^{8, 28} illuminating a sinusoid pattern three times at the same spatial frequency, with phase offsets $a=0$ , $2∕3\pi$ , and $4∕3\pi$ radians. $M_{\mathit{AC}}(x,f_x)$ can then be calculated algebraically at each spatial location, $x_i$ , by:

In the frequency domain, a measurement $M_{\mathit{AC}}\left(f_x\right)$ is the product of (1) the source intensity, $I_0$ ; (2) the MTF of the illumination and imaging optical system, ${\mathit{MTF}}_{\mathit{\text{system}}}$ ; and (3) the true turbid system MTF, $R_d$ :

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.²⁹

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

Fig. 6

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 ${\mu }_a$ and ${\mu }_s^{\prime }$ parameters.

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

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

Fig. 7

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, 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 $\left(R_d\right)$ versus spatial frequency $\left(f_x\right)$ , 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 $x\text{-}y$ spatial position, allowing a robust measurement of the average properties. Compared to “point” illumination measurements, MI samples only the low spatial frequency moments of the transfer function. These low frequencies $(<1{\mathrm{mm}}^{-1})$ are sufficient for separating the absorption from the scattering optical properties, reducing sensitivity to uncertainty inherent in measuring high-frequency spatial moments (i.e., reflectance close to the source).

4.1.

Phantom Experiments

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

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

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

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

4.2.

In Vivo Human Forearm Experiments

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

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

5.1.

Phantom Experiment Results

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

Fig. 8

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 rescaling to higher spatial frequencies (i.e., a decrease in $l^{*}$ ), with contrast apparent at all spatial frequencies.

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

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

Fig. 9

(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\times 15\text{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 ( ${\mu }_s^{\prime }∕{\mu }_a=500$ to ${\mu }_s^{\prime }∕{\mu }_a=8$ ). When absorption is very small, a slight overestimation trend is observed independent of calibration choice. This trend is discussed further in Sec. 5.2. The experiment’s recovered scattering values show less than 10% deviation from the expected value in all cases. Similar linearity is observed in the scattering variation experiment, albeit with slightly more fluctuation. Absorption values in this case demonstrate less than 15% deviation from the expected value, except in the lowest scattering (largest $l^{*}$ ) case. Standard deviations of the recovered $15\times 15\text{pixel}$ optical property maps are mostly less than 1% (maps not shown), indicating both high precision and spatial uniformity over the field of view.

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

Table 1

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

5.2.

Separation of Absorption and Scattering

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

5.3.

In Vivo Human Forearm Results

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

Fig. 10

(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 $\left(0.46{\mathrm{mm}}^{-1}\right)$ compared to the background $\left(0.23{\mathrm{mm}}^{-1}\right)$ ; in the same region, optical scattering showed little contrast, with $\sim 10\%$ spatial variation. While these data show clear separation of absorption and scattering optical contrast, the values derived from a homogeneous model exhibit partial volume effects that diminish the actual absorption contrast of the vein beneath the surface. Our current homogeneous model of reflectance prevents absolute quantitation in the presence of lateral and depth-dependent optical properties. In the next section, however, we attempt to empirically assess these partial volume effects through edge-response experiments in gelatin phantoms.

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

Fig. 11

Venous occlusion data of the volar forearm measurement at $800\mathrm{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) Regionwise 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 ${\mu }_a$ and ${\mu }_s^{\prime }$ requires the corresponding change in reflectance to be above the measurement noise floor. We examined how a homogeneous perturbation in optical properties gives rise to reflectance contrast at each spatial frequency using our Monte Carlo forward model. This was done by taking a numerical derivative of reflectance with respect to a change in a given optical property, generated for normalized frequencies between 0 and 1. In Fig. 12, we show the change in diffuse reflectance $\left(\Delta R_d\right)$ versus normalized spatial frequency resulting from a 1% change in absorption or scattering, for sample $({\mu }_s^{\prime }∕{\mu }_a)$ ratios of 100 (black lines) and 3 (gray lines). The solid and dashed lines denote the sensitivity to absorption and reduced scattering, respectively. (Note: We use a negative absorption perturbation and a positive scattering perturbation in order to produce reflectance contrast of the same polarity.)

Fig. 12

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\text{-}\text{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 $R_d$ . We have added dotted lines to Fig. 12 to show the approximate corresponding camera detection limits with 12- and $14\text{-}\text{bit}$ intensity resolution (i.e., the detection limit for the physical contrast—actual detection limits will depend on the noise characteristics of the particular camera/imaging system). The figure illustrates, for example, that with $({\mu }_s^{\prime }∕{\mu }_a)=100$ at $f_x∕f_{x,\mathit{tr}}=0.1$ , a $12\text{-}\text{bit}$ measurement can resolve (in intensity) a 1% change in scattering but not absorption, while a $14\text{-}\text{bit}$ camera can resolve both. For $({\mu }_s^{\prime }∕{\mu }_a)=3$ , where optical properties are closer in magnitude, this $12\text{-}\text{bit}$ absorption detection criterion occurs at a higher frequency of $f_x∕f_{x,\mathit{tr}}=0.3$ .

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

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

Fig. 13

Preliminary contrast-resolution study. (a) Gelatin phantoms with a $3\times$ step in either absorption (left) or scattering (right) were measured directly (top) and through a $2\text{-}\mathrm{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. 13b for absorption (top) and scattering (bottom) media. The results reveal a diffuse edge-response function that is both depth and optical property dependent. For both absorption and scattering experiments, we observe a degradation of spatial resolution and quantitative contrast through the homogeneous layer (gray lines) compared to that at the surface (black lines). Specifically, the measured contrast values through the homogeneous layer are approximately 15% and 5% of those at the surface, for absorption and scattering, respectively. This level of absorption contrast is comparable with the measured arm vein absorption in Fig. 10 (approximately $0.06{\mathrm{mm}}^{-1}$ ), which corresponds to approximately 18% of whole blood absorption (assuming 70% tissue oxygen saturation). We further observe that the spatial resolution of the recovered scattering maps is consistently better than that for absorption. Put differently, scattering profiles appear to preserve higher spatial frequencies than those of absorption. This property of the measured heterogeneous media is consistent with the homogeneous phantom MTF results shown in Fig. 8, where we noted that absorption contrast appeared mainly at low frequencies, while scattering contrast appeared at all measured frequencies.

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