2.1.

Instrumentation

We used the MI instrument platform first developed by Cuccia, described in detail in Refs. 15, 16, 17, 22, 23 and shown in Fig. 1a . The system consists of an intensity-stabilized $250\text{-}\mathrm{W}$ quartz-tungsten-halogen lamp (Newport, Irvine, California) focused on a spatial light modulator (DMD Discovery 1100, Vialux, Chemnitz, Sachsen, Germany), which is then projected through an objective lens (Thorlabs, Newton, New Jersey) whose magnification and focal length ensure complete illumination of the target by spatially modulated light. A Nuance multispectral imaging system (CRI Inc., Woburn, Massachusetts), composed of a liquid-crystal tunable filter ( $\lambda =650–1100\mathrm{nm}$ , $\Delta \lambda =10\mathrm{nm}$ ) and a $12\text{-}\text{bit}$ CCD camera, was used to acquire reflected images. Projection and collection are cross-polarized to reduce the contribution of specularly reflected light and to isolate diffuse reflectance.²⁴ We placed the samples on a custom-made, $z$ -translation, $x$ - and $y$ -rotation stage that allows positioning at a range of heights from $0\text{to}5\mathrm{cm}$ and tilt angles from $-40\text{to}+40\mathrm{deg}$ . MI data acquisition was performed using code written in LabVIEW (National Instruments, Austin, Texas), and image processing was performed using code written in MATLAB (MathWorks, Natick, Massachusetts).

Fig. 1

Instrumentation schematics: Intensity sine waves (grayscale stripes) are projected at a $20\text{-}\mathrm{deg}$ angle onto a sample using a spatial light modulator, and collected with a $12\text{-}\text{bit}$ CCD camera. Collection is cross-polarized with respect to projection, and wavelengths are filtered using the tunable filter on the camera. Orientation of the sine waves is shown for optical properties measurement (a) and profilometry measurement (b).

2.2.

Tissue-like Phantoms

Tissue-simulating phantoms with known optical properties were used to assess the qualitative and quantitative performance of the imaging technique. These phantoms, $96\times 96\times 10\mathrm{mm}$ , are made of polydimethylsiloxane and incorporate homogeneously distributed India ink as an absorber, and $\mathrm{Ti}{\mathrm{O}}_2$ as a scattering agent.²⁵ Spectral absorption and reduced scattering coefficients of the calibration phantom were verified using a two-distance, multifrequency FDPM measurement, which is a self-calibrating measurement.²⁶

2.3.

Phase Profilometry

Phase profilometry is a well-known measurement technique that is used for surface profile characterization.^{27, 28, 29} As shown in Fig. 1b, a spatially modulated intensity sine wave is projected onto the sample while the camera acquires at a fixed angle with respect to the projection axis. Maximum sensitivity to surface variations is obtained when the fringes of the intensity sine wave are orthogonal to the plane formed by the camera and projector optical axes. In Fig. 1b, the optical axes of the camera and projector are lying within the $\mathbf{x}\text{-}\mathbf{z}$ plane while the fringes are aligned along the $\mathbf{y}$ vector. The observed fringe phase will then be dependent on the geometry of the optical system and the sample surface height. Consequently, the sample height can be extracted from measurement of the fringes’ phase and either knowledge of the optical system geometry or calibration of the phase-height relationship. We chose to extract the fringes’ phase using a profilometry method that consists of projecting an intensity sine wave with three different phase shifts.²⁹ The fringe phase at each pixel is extracted from the three sine waves and unwrapped. Another alternative would be to perform a discrete Fourier transform (DFT) using a single intensity sine wave,²⁸ but in our case, the trade-off between phase resolution and spatial resolution renders the DFT approach inadequate. As the optical system geometry may vary in time, we chose to perform a multiple-height calibration scheme that allows the extraction of a linear relationship between the inverse of sample height variation and the inverse of phase variation.³⁰ This technique consists of vertically translating a horizontal flat surface and acquiring, for known heights, the phase at each pixel location.

It is important to note that the choice of spatial frequency used when performing phase profilometry is a trade-off between sensitivity to large changes of height, which is highest at low spatial frequencies, and height resolution, which is highest at high spatial frequencies. When measuring relatively smooth surfaces that do not present sudden height changes, such as flat or hemispheric surfaces, we performed phase unwrapping (i.e., phase was no longer constrained to $0–180\mathrm{deg}$ ). This permitted use of relatively high frequencies $(>0.1{\mathrm{mm}}^{-1})$ to minimize phase noise, while maintaining sensitivity to large height changes. Of course, this approach assumes that the phase does not vary by $>180\mathrm{deg}$ from one pixel to another, which is the case for smooth surfaces.

2.4.

Modulated Imaging

The spatial dependence of light reflectance from turbid media has been characterized extensively.³¹ Two general approaches can be used to quantify turbid media optical properties using continuous (i.e., time-independent) light sources. The first “real domain” method measures the spatial point-spread function of a point source incident on the sample and is generally well suited for local, millimeter-scale interrogation using contact fiber probes.^{24, 32} This approach can be extended to noncontact geometries and large areas by raster scanning the point illumination source across the region of interest. In the “spatial frequency domain,” the reflectance modulation transfer function of the sample is evaluated in a wide-field imaging geometry.³³ MI, developed by Cuccia and described in detail in Refs. 15, 16, 17, 22, 23, is a spatial frequency domain technique that relies on the analysis of reflectance from a spatially modulated sinusoidal light source. Because the analysis is performed within the spatial frequency domain, fast acquisition of large FOVs (typically, $5\mathrm{cm}$ or greater) is possible. MI is a calibration-dependent, model-based method that allows the extraction of absorption and reduced scattering coefficients ( ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ ; i.e., optical properties) from a turbid medium at any wavelength.

The theory and methods of MI are described extensively in Refs. 16, 17. To place the present study in context, we provide the following brief overview of the method and its fundamental principles. The ac component, $I_{\mathrm{ac}}$ , of an intensity sine wave of spatial frequency $f_x$ and phase $\alpha$ , measured at a single wavelength and at the location $(x,y)$ on the sample has the following form:

During a typical measurement, the modulation amplitude $\left(M_{\mathrm{ac},\mathrm{ref}}\right)$ of a calibration phantom with known optical properties is measured at multiple illumination spatial frequencies $\left(f_x\right)$ , and the calibration phantom diffuse reflectance $\left(R_{\mathrm{d},\mathrm{ref}}\right)$ is predicted from the known optical properties using a forward light propagation model (which is analytical when using the diffusion approximation to the transport equation^{17, 34} or statistical when using Monte Carlo computations³⁵). Next, the $M_{\mathrm{ac}}$ of the unknown sample is measured and its diffuse reflectance extracted using the reference values ( $M_{\mathrm{ac},\mathrm{ref}}$ and $R_{\mathrm{d},\mathrm{ref}}$ ) from the calibration phantom, at each location $(x,y)$

Cuccia¹⁶ and Cuccia ¹⁷ showed that, by properly choosing the spatial frequencies, optical properties can be extracted using a quick, two-dimensional lookup table generated from either “white” Monte Carlo simulations or the diffusion approximation, and calibrated using a phantom with well-known optical properties. The premise of that technique is that dc, continuous wave illumination, is sensitive to both absorption and scattering whereas a certain range of higher spatial frequencies (typically $>0.1{\mathrm{mm}}^{-1}$ in the case of tissues) is primarily sensitive to scattering. This permits separation of the optical properties (i.e., absorption and reduced scattering) of an object. In the work presented here, we used the diffusion approximation to the transport equation as our forward light propagation model and optical properties were extracted using the rapid two-dimensional lookup table generated from the diffusion approximation to the transport equation. During our experiments, the desired high spatial frequency for optical properties measurement using the two-dimensional lookup table was in the range of $0.1–0.2{\mathrm{mm}}^{-1}$ . Because accurate knowledge of this spatial frequency is necessary for the light propagation model, a postacquisition calibration of the exact high-frequency value is measured using a ruler.

As shown in Fig. 1a, when measuring optical properties of phantoms, we purposely chose to project fringes that were not sensitive to the sample surface variations. This was simply achieved by orienting the fringes parallel to the plane formed by the camera-projector optical axes (see Sec. 2.3). Typically, MI is performed by first acquiring data from a single reference calibration phantom with well-known optical properties at a single specific height. Problems arise when the reference and the sample are not at the same height, which occurs commonly during in vivo measurements. In this case, the sample and the reference do not share the same incident light intensities [ $I_0$ in Eq. 3], and direct comparison of modulation amplitudes is no longer appropriate for accurate diffuse reflectance extraction and subsequent optical properties determination. The previous MI measurement protocol, without profile intensity correction, requires that the sample and reference heights be as close as possible to minimize errors.

Additionally, we developed a new, multiple spatial frequency processing algorithm for MI. The commonly employed algorithm assumes that there is only one single spatial frequency of interaction between the light and the sample. Because of surface profile variations, though, this intensity sine wave gets projected onto locally angled surfaces, and multiple frequencies of interaction are present. We overcame this problem by measuring the calibration phantom at multiple single frequencies, extracting the sample optical properties for all those frequencies, and interpolating the sample optical properties according to the measured local spatial frequency (obtained using the surface profile data).

2.5.

Profile-Based Correction Method

Because light intensity locally reemitted from a diffusive surface can be described using a Lambertian model varying in $\cos \left(\theta \right)∕r^2$ , where $r$ is the distance from the sample to the camera and $\theta$ the angle between the normal of the surface and the collection optical axis, we developed a surface profile correction method that consists of these two correction components. The first part addresses the intensity variation associated with surface height changes while the second part addresses the intensity variation due to surface angle. These corrections were applied to MI data. In order to validate this approach, we first assessed its performance in correcting height-dependent intensity variations only. This was made possible by vertically translating a horizontal flat tissue-simulating phantom perpendicular to the collection optical axis. In doing so, there were no angle-dependent contributions to the sample intensity variation. With the knowledge that the height-dependent intensity variations were corrected, we then tilted a flat tissue-simulating phantom, from $-40\text{to}+40\mathrm{deg}$ in both axes of rotation and studied the combination of both height- and angle-dependent corrections, because, when the sample is angled, contributions from height and angle cannot be studied independently.

2.6.

In Vitro Hemispheric Phantom Measurements

In addition to flat phantoms, we acquired data from a complex-shaped phantom that consisted of a flat homogenous phantom having on its surface a $2\text{-}\mathrm{cm}$ radius hemisphere of the same optical properties. One flat phantom with known optical properties was measured at five different heights, from $0\text{to}4\mathrm{cm}$ , every centimeter, for calibration purposes. The sample phantom data were measured at a height of $1\mathrm{cm}$ . Both sets of phantom data were processed for phase profilometry and optical properties as described above. The spatial frequency used for profilometry was $0.17{\mathrm{mm}}^{-1}$ and spatial frequencies used for optical property determination were dc $\left(0{\mathrm{mm}}^{-1}\right)$ and $0.12{\mathrm{mm}}^{-1}$ . FOV was $4\times 5\mathrm{cm}$ . Data were acquired at $670\mathrm{nm}$ and processed for height and optical property reconstructions, with and without profilometry correction (for height and angle). The first phantom served as the calibration reference for both height and optical properties.

2.7.

In Vivo Measurements

We acquired data from a human right hand, with the index finger having a rubber band-induced constriction of blood flow. One flat phantom with known optical properties was measured at six different heights, from $0\text{to}5\mathrm{cm}$ , every centimeter, for calibration purposes. The spatial frequency used for profilometry was $0.08{\mathrm{mm}}^{-1}$ and spatial frequencies used for optical properties were dc $\left(0{\mathrm{mm}}^{-1}\right)$ and $0.13{\mathrm{mm}}^{-1}$ . FOV was $7\times 9\mathrm{cm}$ . Data were acquired at $670\mathrm{nm}$ , and the sample surface smoothed using a Wiener filter ( $10\times 10\text{pixel}$ window size) to reduce the influence of surface angle errors originating from noise in the profilometry data. Data were processed for height and optical property extractions, with and without profilometry correction (for height and angle), and using the flat phantom as the calibration reference. Finally, optical property maps were smoothed using a Gaussian filter having a FWHM of $4\text{pixels}$ .

3.1.

Phase Profilometry

A flat phantom with homogenous, tissue-like optical properties was translated vertically (Figure 2a ), front tilted [Fig. 2b] and side tilted [Fig. 2c], and its surface reconstructed. Over $3\mathrm{cm}$ of vertical translation, measured height was always reconstructed within $1\mathrm{mm}$ or less of its theoretical height [Fig. 2d]. Over $\pm 20\mathrm{deg}$ of side tilt or front tilt, the surface of the phantom was reconstructed within $0.7\mathrm{deg}$ or less of its theoretical angle [Fig. 2e]. However, for front-tilt angles from $\pm 20\text{to}\pm 40\mathrm{deg}$ , there was proportionally increasing deviation due to systematic errors associated with phase extraction [Fig. 2e]. For front tilts, the phase of the sine wave is the result of a phase addition process along the $\mathbf{x}$ vector, causing errors to accumulate during phase extraction. This does not occur during side-tilts, where phase is constant along the $\mathbf{x}$ vector.

Fig. 2

Profilometry measurements: (a) Surface maps of flat phantoms translated vertically, b) surface maps of flat phantoms rotated along the $y$ -axis (i.e., front tilt), (c) surface maps of flat phantoms rotated along the $x$ -axis (i.e., side tilt), (d) measured surface heights (squares) as a function of the expected theoretical surface heights (dashed line), and (e) measured surface angles for front tilt (squares) and side tilt (circles) as a function of the expected theoretical angles (dashed line).

3.2.

Optical Property Measurements of Flat Tissue-like Phantoms

Extraction of absorption and reduced scattering coefficients was impacted significantly by profilometry correction on a flat phantom with homogenous, tissue-like optical properties (Fig. 3 ). As described in Sec. 2, results from uncorrected data were processed using a calibration phantom measured at a single $0\text{-}\mathrm{cm}$ reference height and $0\mathrm{deg}$ in both axes of rotation. Results from corrected data have been processed using a calibration phantom measured at various heights and $0\mathrm{deg}$ in both axes of rotation.

Fig. 3

Optical property measurements of tissue-like phantoms: Recovered absorption coefficients [(a–c), $\text{mean}\pm \text{standard}$ deviation] and reduced scattering coefficients [(d–f), $\text{mean}\pm \text{standard}$ deviation] at $650\mathrm{nm}$ of (a,d) vertically translated phantoms, (b,e) front-tilted phantoms, (c,f) side-tilted phantoms having homogeneous optical properties. Profilometry-corrected data (squares) and uncorrected data (circles) are shown along with the expected theoretical values (dashed lines).

For vertical translation [Figs. 3a and 3d], the correction algorithm uses only height-dependent intensity correction. Profilometry-based correction reduced absorption coefficient error from 10% per centimeter to $<1\%$ over the entire range from $0\text{to}3\mathrm{cm}$ . Similarly, reduced scattering coefficient error was reduced from $\approx 10\%$ to $<1\%$ at every surface height. For both absorption and reduced scattering coefficients, the standard deviation of all pixels in an image was reduced from $⩾7\%$ to $⩽3\%$ of the mean pixel value.

For tilted phantoms [Figs. 3b, 3e, 3c, 3f], the correction algorithm uses a combination of height- and angle-dependent intensity correction. For angles within $\pm 40\mathrm{deg}$ , the correction reduced absorption coefficient error from 86% to 7% of the expected pixel value, and reduced scattering coefficient error from 10% to 4% of the expected pixel value. For absorption coefficients, the standard deviation of all pixels in an image was reduced from 13% to 6% of the mean pixel value, whereas for reduced scattering coefficients, standard deviation was reduced from 6% to 2% of the mean pixel value. It should be noted that the camera optical axis was not normal to the $0\mathrm{deg}$ horizontal phantom surface, but was angularly displaced ( $-20\mathrm{deg}$ front, $0\mathrm{deg}$ side). Consequently, measured values found their local extrema at this angle [see, for example, Fig. 3b and 3e]. We also note a tendency for overcorrection of tilted phantoms on the reduced-scattering maps (see Sec. 4).

3.3.

Optical Property Measurements of Hemispheric Tissue-like Phantoms

Although appearing to work well with flat surfaces, we next explored the use of profilometry-based correction on complex surfaces. Shown in Fig. 4 is a hemispheric phantom having homogeneous, tissue-like optical properties. Figure 4a shows optical properties maps (absorption and reduced scattered coefficients at $670\mathrm{nm}$ ) plotted on the sample surface, with and without correction. Figure 4b shows cross-sectional plots from Figure 4a made through the center of the hemisphere. Profilometry-corrected data (solid curve) and uncorrected data (thin dashed curve) are shown along with the expected theoretical values (thick dashed curves). In the absence of profilometry-based correction, reconstruction of absorption and reduced scattering coefficients results in large errors (28% of the pixels on average are within 10% of their expected values). These errors are most significant at points presenting a high surface angle, with deviations of $>100\%$ at angles of $>40\mathrm{deg}$ for absorption coefficients. However, when corrected for surface height and angle, the accuracy of optical property measurements is improved both qualitatively [Fig. 4a] and quantitatively {50% of the pixels on average are within 10% of their expected values [Fig. 4b]}.

Fig. 4

Optical property measurements of a hemispheric, tissue-like phantom having homogeneous optical properties: (a) Recovered absorption coefficients (top) and reduced scattering coefficients (bottom) at $670\mathrm{nm}$ either without (left), or with (right), profilometry-based correction for height and angle. Arrows indicate the expected theoretical value. (b) Cross-sectional plots from (a), through the center of the hemisphere, of recovered absorption coefficients (left) and reduced scattering coefficients (right) at $670\mathrm{nm}$ . Profilometry-corrected data (solid curve) and uncorrected data (thin dashed curve) are shown along with the expected theoretical values (thick dashed curves).

When reconstructing only those pixels that fall within the range of the calibration phantom ( $0–3\mathrm{cm}$ in vertical translation and $\pm 40\mathrm{deg}$ in front tilt and side tilt), deviation of the recovered absorption coefficient from the expected value was reduced from 65% to 9% with profilometry-based correction. Deviation of the recovered reduced scattering coefficient from the expected value was actually slightly degraded after correction, increasing from 1.5% to 3% over the entire hemispheric surface. There was also a trend toward overcorrection with increasing surface angle [Fig. 4b]. The standard deviation of measurements at any given pixel was reduced from 51% to 25% of the mean pixel value for absorption coefficient, and from 20% to 7% of the mean pixel value for reduced scattering coefficient, with the incorporation of surface height maps.

When the entire hemispheric surface is analyzed, without restricting reconstruction to $\pm 40\mathrm{deg}$ [Fig. 4b], absorption maps are still improved by a similar factor but corrected reduced scattering results are slightly degraded, increasing to 8.4% average deviation from expected and showing no improvement in standard deviation over the uncorrected case. This is thought to be due to very steep angles at the junction of the sphere and the flat phantom, resulting in both insensitivity to optical properties and inappropriate Lambertian reflectance correction (see Sec. 4).

Note is made of a tiny dark spot at the top center of the hemisphere, seen in both the absorption and reduced scattering coefficient maps [Fig. 4a], which is caused by excessive specular reflection. A second tiny dark spot in the reduced scattering coefficient maps is due to a small speck of nondissolved $\mathrm{Ti}{\mathrm{O}}_2$ , visible to the naked eye, which arose during the phantom-making process. Note is also made of an area at the back of the hemisphere where the profile cannot be reconstructed due the low intensity of the illumination. This phenomenon is typical when performing phase profilometry using single-angled illumination.³⁶ We also note a tendency for overcorrection of both scattering and absorption maps, increasing with surface angle (see Sec. 4).

3.4.

In Vivo Optical Property Measurements

The optical properties of an even more complex object, a human hand with the index finger constricted with a rubber band, were reconstructed with and without phase profilometry-based correction. A reflectance image of the hand and its surface map is shown in Fig. 5a . Recovered absorption and reduced scattering coefficients, with and without surface map correction, are shown in Fig. 5b. Consistent with in vitro results, the standard deviation of measurements made on the constricted finger section was reduced from 44% to 30% of the mean pixel value for absorption coefficient, and from 29% to 23% of the mean pixel value for reduced scattering coefficient. Artifacts due to movement of the hand during data acquisition were visible, especially in scattering maps (see Sec. 4). Additionally, noise originating from phase profilometry acquired at low spatial frequency is clearly visible and degrades the results obtained from angle-dependent correction (see Sec. 4). Most important, profilometry-based correction greatly improved the uniformity of in vivo optical property measurements, as shown in the cross-sectional plot, along the $y$ -axis of the constricted finger, of the absorption data [Fig. 5c].

Fig. 5

In vivo imaging of a human hand: (a) Reflectance image of a right hand having an index finger constriction (rubber band), taken at $670\mathrm{nm}$ (left), and the 3-D plot of the hand obtained using phase profilometry (right). Dashed box indicates the region of interest for cross-sectional quantitative analysis. (b) Recovered absorption coefficients (top) and reduced scattering coefficients (bottom) at $670\mathrm{nm}$ either without (left), or with (right), profilometry-based correction for height and angle. Note increased absorption at $670\mathrm{nm}$ on the constricted finger due to an increase in deoxyhemoglobin and total hemoglobin. (c) Cross section along the $Y$ -axis through the constricted finger [ $X=\text{average}$ of pixels 275 to 300; see dashed box in (a)] of the hand shown in (b). Corrected (solid curve) and uncorrected (dashed curve) absorption coefficients are plotted.