2.1.

Spatial Frequency Domain Imaging

SFDI instrumentation, data acquisition, and data analysis have previously been described in detail.⁵ Briefly, structured light is projected onto a sample using a spatial light modulator (SLM), and a camera detects the diffusely reflected light emitted from the boundary of the sample. In three-phase SFDI, three frames of data are acquired at relative modulation phases of 0, 120, and 240 deg for each AC spatial frequency. These phase-offset images are applied to a simple formula to extract the AC information content pixel-by-pixel, as shown in Eq. (1). A simulation of this demodulation approach using three phase-offset modulation patterns is shown in Fig. 1. All data processing and computation used to produce figures were performed using the MATLAB® software suite (MATLAB and Statistics Toolbox Release 2011b, The MathWorks, Inc., Natick, Massachusetts).

Fig. 1

(a) Simulation of three-phase demodulation technique on a turbid sample with uniform scattering and a circular absorbing lesion in the center of the field of view. Images are acquired at relative modulation phases of 0, 120, and 240 deg. These images are then applied to a three-phase demodulation formula [Eq. (1)], which extracts the AC information content from the sample pixel-by-pixel. (b) Plot of reflectance cross-sections taken from the center row (dashed line shown in a) of each AC + DC image, and the resulting AC + DC information envelope (left). A cross section of the demodulated DC (taken from additional planar image) and AC information content (right) is also shown. The DC cross section is more sensitive than AC to absorption, whereas the AC cross section is more sensitive than DC to scattering.

Next, the demodulated intensity data from the sample is calibrated to that of a tissue-simulating phantom having known optical properties. This calibrated diffuse reflectance data at each spatial frequency is then applied to a light transport model such as diffusion or Monte Carlo, from which ${\mu }_a$ and ${{\mu }_s}^{\prime }$ maps are determined. Finally, these maps are generated at several wavelengths and are fit to known chromophore spectra to quantitatively determine the concentration of relevant chromophores in the sample.

In order to decouple scattering from absorption, at least two spatial frequencies are required. In the case where the minimal number of frames are taken, 0 (planar illumination) and $0.2{\mathrm{mm}}^{-1}$ are typically used. Since there is no spatial variation in the spatial modulation pattern at $0{\mathrm{mm}}^{-1}$, only one phase is required, so a minimum of four frames (per wavelength) total are needed in three-phase SFDI, consisting of a single DC frame and three phase-offset AC frames. After demodulation, a fast lookup table is employed to determine ${\mu }_a$ and ${{\mu }_s}^{\prime }$. The primary limitation with respect to the implementation of real-time SFDI is data acquisition time. In particular, as we currently practice the technique, the need to acquire three frames of data for each AC spatial frequency limits the speed of SFDI data acquisition. We have developed a technique to address this bottleneck by reducing the number of frames required for each modulating frequency from three to one.

2.2.

Two-Dimensional Hilbert Demodulation Technique

The Hilbert transform is a ubiquitous tool in signal processing with a wide variety of applications in the communication field.¹² The general principle is that a modulating double-sideband signal such as a sine or a cosine contains redundant information; only one sideband is needed to extract the modulated information content. Using the Hilbert transform, one can derive a single-sideband expression for this modulated signal with no loss of information. This single-sideband expression allows for the extraction of the demodulated information content and phase map of the modulated signal. Recently, the concept of applying the Hilbert transform using spiral phase functions in 2-D Fourier space to demodulate 2-D curved patterns in space was developed by Larkin et al.^13,14 We have adopted the general concept of this approach and have applied it to the SFDI workflow.

The modulated reflectance images obtained in SFDI can be described by Eqs. (2) and (3), where $f_{x,y}$ is the modulating spatial frequency, and ${\varnothing }_{x,y}$ is the phase.

The purpose of demodulation is to extract the AC diffuse reflectance term $R_{\mathrm{AC}}(x,y)$ from the detected amplitude $I(x,y)$. Using Euler’s theorem, a cosine function can be expressed as the sum of two complex exponentials or sidebands in the frequency domain. As mentioned previously, the Hilbert transform is used to obtain a single-sideband expression for a double-sideband function such as a cosine. Since a single-sideband function can be expressed as a complex exponential, demodulation is straightforward. The magnitude of the single-sideband expression for SFDI modulation results in the diffuse reflectance we wish to obtain.

Our new SFDI demodulation approach employs a 2-D Hilbert transform to SFDI frames by applying a spiral phase function to the image in 2-D Fourier space. One unique aspect of this approach is that it can demodulate frames whose modulation patterns are rotated or arbitrarily oriented. That is, the wavenumber of the modulating pattern can have arbitrary directionality with respect to the lateral imaging axes $(x,y)$. The spiral function is described in Eq. (4)

In order to implement the Hilbert demodulation technique, the following steps are performed. First, the DC component of the spatially modulated image which consists of both AC and DC components [$I(x,y)$ from Eq. (2)] is removed. A 2-D fast Fourier transform (FFT) is then applied to the resulting AC image [$M(x,y)$ from Eq. (2)]. In 2-D Fourier space, the transformed AC image is multiplied by a map generated using Eq. (4), having the same dimensions as the AC image. Next, an inverse FFT is applied to this product. The resulting image is similar to the original AC image, except that the modulating “cosine” is now a “sine,” i.e., the phase of the modulating wave is shifted by 90 deg. Then, the magnitude of this “sine” image is taken, which accounts for the complex contribution of the transformed map due to the orientation angle of the modulating wave, shown in Eq. (4). The resulting term, $H(x,y)$, represents the Hilbert transform of the original AC image [$M(x,y)$]. Finally, $H(x,y)$ is multiplied by the complex unit and added to the AC component of the original AC image. The resulting magnitude is the demodulated AC diffuse reflectance, denoted in Eq. (5) by $R(x,y)$

A walkthrough of the Hilbert technique is shown in Fig. 2 using a simulated DC and AC + DC image. In this simulation, the sample is highly reflective, such that the sinusoidal pattern is kept intact as the light reaches the boundary of the sample. In reality, we apply this technique to turbid tissue simulating samples or biological samples, which will be demonstrated in the experimental results, but this virtual sample was chosen to clearly illustrate the Hilbert demodulation concept. We expect that the Hilbert technique will yield results comparable in quality to those shown in Fig. 2 so long as the reflected modulation pattern is greater in amplitude than the camera noise. This constraint should typically be satisfied when looking at biological samples.

Fig. 2

Simulation of the two-dimensional (2-D) Hilbert demodulation method on a highly reflecting surface. First, the DC component of the modulated image is removed, and a fast Fourier transform (FFT) is performed on the AC + DC image. The resulting 2-D map in Fourier space is then multiplied by a spiral phase function, consisting of a continuous, radially-varying map ranging in value from $-1$ to $+1$ in real and imaginary space. An inverse FFT is performed on the map, resulting in an image whose magnitude is the original modulated image phase-shifted by 90 deg. This image is multiplied by the imaginary unit and added to the original image. The magnitude of this image results in the demodulated diffuse reflectance of the AC component from the original AC + DC image.

Here, we begin with a DC image at a uniform intensity of one, and an AC + DC image [Eq. (2)] with an intensity varying from zero to one, with a modulation pattern oriented diagonally. First, the DC component is removed from the AC + DC image, and an FFT is performed on the resulting AC image. Next, the spatial frequency map of the transformed AC image is multiplied with the complex spiral function map. The resulting map is then FFT inverted, and the magnitude is taken. This “magnitude” image is then multiplied by the imaginary unit and added to the initial AC image (before Hilbert transform). The magnitude of this sum results in the demodulated AC diffuse reflectance, which is uniform at an intensity of $\sim 0.5$.

It should be noted that the demodulated AC images obtained in Fig. 1 and subsequent figures using the Hilbert technique contain residual ringing artifacts, which are due to the fact that the spatially projected sinusoidal patterns are cutoff by the boundaries of the image and are, therefore, finite in length. We refer to the degradation in image quality as resulting from ringing artifacts. In theory, these artifacts should increase in severity as the number of periods in the modulation pattern in the image decreases. Therefore, one potential strategy to minimizing these artifacts, particularly at lower spatial frequencies, is to use an SFDI instrument having a large field-of-view. Alternatively, a window function such as a Gaussian or Hamming filter could be applied to the image in postprocessing,¹⁵ which will be the subject of future work.

We performed a side-by-side comparison of our advanced Hilbert demodulation technique to three-phase demodulation. To generate the data used to produce the images analyzed in this section, we employed a clinical SFDI system at a wavelength of 658 nm. This system has been described previously.¹⁶ In our first experiment, we compared diffuse reflectance maps obtained on a tissue-simulation phantom at multiple spatial frequencies. This phantom consists of a silicone foundation consisting of India ink as an absorbing agent and titanium dioxide as the scattering agent. The fabrication technique used to make these phantoms has been described previously.¹⁷ Next, we compared ${\mu }_a$ and ${{\mu }_s}^{\prime }$ maps extracted from an in vivo human forearm. Finally, we evaluated SOI maps taken on a structural orientation phantom consisting of a silicone-based bottom layer (described above), and a top layer composed of sections of pleated air filters at various orientation angles. The SOI of these phantoms has been previously evaluated.¹¹

3.1.

Tissue Phantom Reflectance Experiment

Demodulation in the SFDI workflow allows for the extraction of information content in the SFD, which is used to generate ${\mu }_a$ and ${{\mu }_s}^{\prime }$ maps. Figure 3 shows a comparison between three-phase SFDI and the Hilbert technique of demodulated diffuse reflectance maps of a homogeneous tissue-simulating phantom. Maps of demodulated reflectance at an AC spatial frequency of $0.2{\mathrm{mm}}^{-1}$ are shown in the figure. Here, only the first phase (0 deg) intensity image is applied to the Hilbert technique, whereas intensity images at three phases (0, 120, and 240 deg) are applied to three-phase demodulation [shown in Eq. (1)]. Average demodulated diffuse reflectance results are shown at five spatial frequencies evenly distributed from 0 to $0.2{\mathrm{mm}}^{-1}$ which is taken from the region of interest (ROI) shown in the black box. These spatial frequencies are typically employed in the SFDI workflow, and instrumentation and models have been shown to perform adequately in this range.^5,18 Here, we see good demodulation quality across the entire field-of-view, with pixel intensity differences between three-phase and Hilbert generally within 5%. We also show agreement in diffuse reflectance values between the two techniques, with mean reflectance values within 1% within the ROI for spatial frequencies of 0.05, 0.1, 0.15, and $0.2\mathrm{m}{\mathrm{m}}^{-1}$. It should be noted that, although the sample in this case is homogeneous, the reflectance intensity over the field-of-view is not. This is due to the inhomogeneity of the light source, which is accounted for during calibration.

Fig. 3

Demodulated reflectance images at $0.2\mathrm{m}{\mathrm{m}}^{-1}$ of a tissue-simulating phantom using (a) conventional, three-phase spatial frequency domain imaging (SFDI) (4 frames), and (b) the advanced, Hilbert-based technique (2 frames). (c) Map of percent difference in demodulated reflectance between three-phase SFDI and the Hilbert-based technique. (d) Plot of mean diffuse reflectance versus spatial frequency for the region of interest (ROI), shown in (a).

Accurate demodulation, and thus determination of AC information content in the SFD, is a necessary component of SFDI and is what allows for the generation of ${\mu }_a$ and ${{\mu }_s}^{\prime }$ maps. In the following section, we present ${\mu }_a$ and ${{\mu }_s}^{\prime }$ maps extracted from a volar forearm using the Hilbert demodulation technique and compare this data directly to ${\mu }_a$ and ${{\mu }_s}^{\prime }$ maps derived using three-phase demodulation.

3.2.

In Vivo Volar Forearm Experiment

Optical property maps of a human volar forearm were calculated using a fast lookup table method, employing spatial frequencies of 0 and $0.2\mathrm{m}{\mathrm{m}}^{-1}$. Figure 4 shows a comparison of ${\mu }_a$ and ${{\mu }_s}^{\prime }$ maps of the forearm using the three-phase (left column) and Hilbert (right column) demodulation techniques. Here, we see agreement in optical property values between the two techniques, with the difference in mean ${\mu }_a$ and ${{\mu }_s}^{\prime }$ values in the ROI (shown in the black box) being 0.2% and 0.15%, respectively.

Fig. 4

In vivo optical property results taken from a human volar forearm. (a) Absorption (${\mu }_a$) and (b) reduced scattering (${{\mu }_s}^{\prime }$) coefficient maps derived from 4-frame, three-phase SFDI (left) and the 2-frame, Hilbert (right) demodulation techniques. For ${\mu }_a$ and ${{\mu }_s}^{\prime }$, the difference in optical property calculations over the region of interest (ROI, shown in black box) is 0.2% and 0.15%, respectively.

A key motivation for spatially modulating light in SFDI is to decouple scattering from absorption. Therefore, quantitative optical property mapping is an essential feature of the technique. Figure 4 demonstrates the ability of our new demodulation technique to produce ${\mu }_a$ and ${\mu }_{\mathrm{s}}\prime$ maps in biological tissue that agree with the conventional demodulation method, which suggests that we can reduce this new technique to practice. The payoff using the new Hilbert technique is the reduction in frames of data required to derive ${\mu }_{\mathrm{a}}$ and ${{\mu }_s}^{\prime }$, and thus an increase in data acquisition speed. In the case shown in Fig. 4, the number of frames reduced using the Hilbert technique over a three-phase demodulation is from four to two, resulting in a twofold increase in imaging speed. However, the payoff in speed increases further if more spatial frequencies are employed, asymptotically approaching a threefold increase.

Although the primary benefit of using the Hilbert technique over a three-phase demodulation is the increase in imaging speed, there are additional benefits. In particular, as shown in Fig. 2, the Hilbert technique can demodulate rotated or oriented sinusoidal patterns using only one AC phase. Acquiring reflectance data at multiple sinusoidal pattern orientation angles is used to characterize tissue structural orientation in the SFD. The following section shows orientation angle and contrast maps on tissue structural orientation phantoms using both the Hilbert and the three-phase demodulation techniques.

3.3.

Scattering Orientation Experiment

By rotating SFDI modulation patterns and acquiring diffuse reflectance maps at several projection angles, we can determine the orientation angle and magnitude of structures in biological tissue.¹¹ Therefore, if we wish to use the Hilbert technique to probe tissue orientation, we must verify that it can produce similar structural orientation contrast as three-phase SFDI. Figure 5 shows scattering orientation angle and contrast maps for the three-phase and Hilbert techniques at an AC spatial frequency of $0.2\mathrm{m}{\mathrm{m}}^{-1}$. The orientation angle of the structure being probed is determined by the angle at which minimum diffuse reflectance is detected. We see here that the average orientation angle for all three ROIs (shown in white boxes) is within 1 deg, which is well within our angular resolution of 5 deg.

Fig. 5

Structural orientation results on structural orientation phantoms consisting of air filters with known structural anisotropy. Orientation angle maps derived from demodulated reflectance images using the (a) three-phase and (b) the Hilbert technique. Scattering orientation index (SOI) maps using (c) three-phase SFDI and (d) the Hilbert technique. ROIs were analyzed in three filters (white boxes). The difference mean orientation angle determined by the left, top, and bottom ROIs is 0, 1, and 0.75 deg, respectively. The difference in mean SOI in the left, top, and bottom ROIs between the Hilbert technique and conventional SFDI is 7.8%, 1.7%, and 0.27%, respectively.

To assess the degree to which the underlying structures are oriented, we used a normalized quantity known as the SOI, shown in Eq. (7).¹¹ Here, the SOI is determined by maximizing, for all projection angles, the reflectance taken at a given angle subtracted from the reflectance taken at the orthogonal projection angle, divided by the sum. The tissue structure orientation phantoms used consist of rectangular-shaped, pleated air filters having significant structural orientation, placed on top of a tissue-simulating phantom having minimal structural orientation

Figure 5 shows SOI maps of tissue structural orientation phantoms using the three-phase (left column) and Hilbert (right column) demodulation techniques. In general, the SOI values obtained using the Hilbert technique are within 10% of those obtained using three-phase demodulation. In particular, the mean difference in the SOI is well within 10% for the three ROIs, and within 2% for two out of three of the ROIs. This demonstrates an overall agreement in SOI between results obtained using the Hilbert and three-phase techniques.

Since the characterization of structural orientation in SFDI uses multiple projection angles of sinusoidal patterns, several frames of data are required. In the case shown in Fig. 5, the angular resolution in orientation analysis is 5 deg. Since the orientation angle has a range of 0 to 180 deg, 36 projection angles were employed. Using the three-phase technique, this results in a total of 108 frames, whereas the Hilbert technique requires only 36 frames (one frame per projection angle). Thus, the Hilbert technique in this example increases imaging speed threefold over three-phase demodulation.