We report a novel demodulation method that enables single snapshot wide field imaging of optical properties of turbid media in the Spatial Frequency Domain (SFD). This Single Snapshot Multiple frequency Demodulation (SSMD) method makes use of the orthogonality of harmonic functions to extract the modulation transfer function (MTF) at multiple modulation frequencies simultaneously from a single structured-illuminated image at once. The orientation, frequency, and amplitude of each modulation can be set arbitrarily subject to the limitation of the implementation device. We first validate and compare SSMD to the existing demodulation methods by numerical simulations. The performance of SSMD is then demonstrated with experiments on both tissue mimicking phantoms and *in vivo* for recovering optical properties by comparing to the standard three-phase demodulation approach. The results show that SSMD increases significantly the data acquisition speed and reduces motion artefacts. SSMD exhibits excellent noise suppression in imaging as well at the rate proportional to the square root of the number of pixels contained in its kernel. SSMD is ideal in the implementation of a real-time spatial frequency domain imaging platform and will open up SFDI for vast applications in imaging and monitoring dynamic turbid medium and processes.

## I. INTRODUCTION

Spatial Frequency Domain Imaging (SFDI) is one emerging non-contact and quantitative wide field modality for rapid mapping of the optical properties of turbid media such as biological tissue.^{1–3} It has attracted significant attention recently as SFDI has been proved effective in probing the micro-architecture of the sub-surface of the turbid medium over a wide field^{4} and has promising clinical applications in hemodynamics monitoring,^{5,6} brain,^{7} burns,^{8} pathology evaluations,^{9} and surgery guidance.^{10} SFDI has also been applied to reduce the absorption contamination in fluorescence planar imaging.^{11} One important issue in SFDI is the imaging speed and motion artefacts.^{12} In order to decouple scattering from absorption, the measurement of the modulation transfer function (MTF) of the turbid medium at no less than two modulation frequencies is required. The standard three phase demodulation method requires three images at different phase delay (0, 2π/3, and 4π/3) to compute MTF at one spatial frequency (one AC component). A minimum of 4 images (one DC and one AC) or 6 images (two AC components of different frequency) will be needed to recover the optical property map, which in turn hinders SFDI in real time applications such as when monitoring fast temporal dynamics.

Recently, several fast demodulation methods have been proposed to improve upon the standard 3-phase demodulation approach for SFDI. The single snapshot optical properties method (SSOP) developed by Vervandier and Gioux^{2} uses a line by line Fourier transform of an SFDI intensity image and processes in the frequency space to extract one DC and one AC image from a single frame of data to recover the optical property. Although this approach reduces notably the data acquisition time, one needs to pick a cutoff frequency to divide the spatial frequency spectrum into the DC and AC components. An alternative approach proposed by Nadeau et al^{13} uses a 2-D Hilbert transform by employing a complex-valued spiral phase function in the Fourier space. The Hilbert demodulation method requires one single frame of image to recover one AC component at the specific frequency. Although it increases the data acquisition speed by two or threefold versus the standard three-phase approach, a minimal of two frames of images are still needed to generate the optical property map. This Hilbert transform method was later extended to Multifrequency Synthesis and Extraction (MSE) which extracts one fundamental and its harmonic AC components using square wave projection patterns.^{14} The multiple AC components in MSE are the Fourier components of the square wave. The relative strength of the AC components is hence fixed and decays with the frequency. This will limit the sensitivity and accuracy of the recovery of the AC components of higher frequency in MSE as the MTF of a turbid medium drops rapidly with the modulation frequency.

In this work, we present a novel single snapshot multiple frequency demodulation (SSMD) method to extract one DC and multiple AC components simultaneously from a single structured-illuminated image at once using the orthogonality of harmonic functions. The multiple AC components can have arbitrary orientations, amplitudes, and modulation frequencies subject to the limitation of the implementation device, offering the design freedom to achieve the best signal to noise performance over the medium response at all the modulation frequencies of interest. After the validation and the comparison of SSMD to the existing demodulation methods by numerical simulations, the performance of SSMD is demonstrated with experiments on both tissue mimicking phantoms and *in vivo* for recovering optical properties by comparing with the standard three-phase approach. SSMD not only is much faster and reduces motion artefacts but also exhibits excellent noise suppression in imaging as well. The application of SSMD may open up SFDI for real time mapping of the optical properties of a dynamic turbid medium and monitoring of fast temporal evolutions.

## II. MATERIAL AND METHODS

### A. Spatial frequency domain imaging

Before introducing SSMD, we will briefly review the principle of SFDI.^{3} SFDI maps the optical properties (the absorption coefficient $\mu a$ and the reduced scattering coefficient $\mu s\u2032)$ of a turbid medium from its reflectance images under structured illuminations. Assuming an incident intensity pattern of $I0=IDC(0)+IAC(0)cos(2\pi fxx+2\pi fyy+\varphi )$ in SFDI, the backscattering light intensity can be described as:

where *f*_{x}, *f*_{y} are the modulating spatial frequencies along the *x* and *y* directions, respectively, and *ϕ* is the phase. The demodulation is to extract the AC amplitude $IAC(x,y)$ from $I(x,y)$. The ratio $IAC/IAC(0)$ represents the modulation transfer function (MTF) of the turbid medium at the spatial frequency $(fx,fy)$ after removing the imaging system response characteristics through calibration. In the standard three-phase demodulation, three images are acquired at three equally separated phases *ϕ*_{1,2,3}=(0, $2\pi /3$, $4\pi /3)$ for each AC spatial frequency. The AC amplitude $IAC(x,y)$ is solved pixel-by-pixel as

where *I*_{1,2,3} is the backscattered light intensity corresponding to the phases *ϕ*_{1,2,3}, respectively. Once the MTF of the turbid medium at two or more modulation frequencies is obtained, the optical properties of the medium are then determined by comparing to the light transport model such as radiative transfer or diffusion in the Fourier domain.^{3,15–17}

### B. Single snapshot multiple frequency demodulation technique

The main power of structured illumination with digital micromirror devices (DMDs) is that the pattern of the incident light on the surface of a specimen can be freely controlled and altered with precise knowledge. SSMD adopts a structured illumination consisting of multiple harmonics simultaneously. The light pattern on the surface of the specimen and the backscattered light intensity is given, respectively, by:

where *k*≥1 is the number of AC components, $IAC,i(0)$ and $IAC,i$ are, respectively, the incident and backscattering amplitudes of the *i*-th AC component modulated at the spatial frequency $(fx,i,fy,i)$ and of phase $\varphi i$. Note that the spatial modulation frequencies can be prescribed exactly by DMD (subject to the limitation of DMD imposed by its the pixelated nature and size) once the SFDI system has been setup and calibrated.

The key in SSMD is that the *k* different spatial frequencies are chosen to take the form of $fx,i=mi/T1$ and $fy,i=ni/T2$ where $\u2212T1<mi<T1$, $\u2212T2<ni<T2$, and $mi$, $ni$, $T1$, and $T2$ are all integers. In other words, $T1$ and $T2$ define the smallest common multiples of all the AC modulating spatial periods along the horizontal and vertical directions, respectively. As in the standard three-phase demodulation method, $IAC,i(x,y)$ is further assumed to be slow varying and can be treated as a constant over the sliding window of size *T*_{1} × *T*_{2} pixels. Each AC component $IAC,i$ is given by

where the integration is performed over the sliding windowσ. Here, in order to remove the requirement of the knowledge of the phase $\varphi i$, a sinusoid and cosine kernel with the same spatial modulation frequency $(fx,i,fy,i)$ and of size *T*_{1} × *T*_{2}, respectively, is applied to filter $I(x,y)$, and the two resulting images are then combined to recover the corresponding AC component.

A schematic diagram of the SSMD is illustrated in Fig.1. As one example, the superimposed image consists of one uniform DC background, a “中” AC image and a “华” AC image, where “中” is modulated horizontally at the spatial frequency *f*_{1} and “华” is modulated along the 45^{o} direction at the spatial frequency *f*_{2}. First, we determine the sliding window (kernel) size *T*_{1} × *T*_{2} that *T*_{1} equals the least common multiple of $1/f1$ and $2/f2$, and *T*_{2} equals $2/f2$. A sinusoid and cosine kernel of size *T*_{1} × *T*_{2}, respectively, is then applied to filter the superimposed image and the two resulting images are combined to recover the corresponding AC component according to Eq. (5). The DC image is extracted by a direct averaging over the $T1\xd7T2$ window (i.e., a unity kernel). Each individual image of “中” and “华” is then demodulated with high fidelity (the rightmost column of Fig. 1).

### C. Experimental setup

The experimental setup is shown in Fig. 2. The light source is Digital Micromirror Device (DMD, DLP®LightCrafter^{TM} 4500, Texas Instruments) and the red channel (623nm) was used in the reported experiments. The DMD is operated in the pattern mode and the trigger time for one displayed image set at 10ms. A Canon 5D Mark III camera is used to record the reflectance images with its exposure time set at 100ms. Our SFDI imaging system is designed to achieve a high resolution imaging over a field of view of centimeter squared (hence significantly different from the conventional SFDI system^{1,3}). A pellicle beam splitter is used to integrate the modulated illumination subsystem and imaging subsystem without ghosting. The sample is conjugate to both the DMD and the camera sensor chip. The magnification power of the illumination subsystem and the imaging subsystem is 2 and 0.81, respectively. The illuminated area on the surface of the specimen was 13 × 13mm in the reported experiments in which 1mm corresponds to 34 pixels in DMD and 130 pixels in the camera. The angle of view for the 13 × 13mm imaging area is 3.33 degrees. The sample is slightly tilted (∼5 degrees) to avoid the specular reflection to be imaged.

As the output light intensity from the DMD is not linear, calibration was first performed to correct this nonlinearity under planar illumination by comparing to the diffuse reflectance from a Lambertian reflectance standard. We then placed a camera exactly in the position of the sample and directly measured the structured-light field to obtain the MTF of the illumination subsystem over the spatial frequency range 0∼10 line pairs/mm. The MTF of the imaging subsystem was also measured by imaging resolution test targets (R1L1S1P, Thorlabs, Inc.) over the same spatial frequency range. The system MTF at each spatial frequency we used in experiments was computed by interpolation based on the above MTF calibration of the illumination and imaging subsystems.

## III. RESULTS

### A. Numerical simulations

In order to validate SSMD, numerical simulations are performed to compare the results obtained with the standard three-phase, the Hilbert demodulation, the SSOP, and the SSMD methods. We first generate an AC + DC image with various levels of Poisson noise (the case of 10% Poisson noise is displayed in Fig. 3(a)). The “华” AC component is modulated horizontally with a period of 5 pixels. Demodulation results obtained from the four different methods are exhibited in Fig. 3(b–e) for the DC component and Fig. 3(f–i) for the AC component. Table I summarizes the performance of recovering the AC component. It is evident the standard three-phase method performs best when the noise is low whereas SSMD demonstrates excellent noise suppression capability and achieves the best accuracy among the four methods with increasing noise (noise ≥ 5%). The recovered DC components by methods other than the standard three-phase approach contain the residue from the AC component with the recovery by SSOP showing the most amount of residue. We then increase the resolution of the image by 5 times and the modulation period of the “华” AC component also by 5 times to 25 pixels. The image contaminated by different levels of Poisson noise is again demodulated by the four different methods. The demodulated results show little difference from Fig. 3 and the performance comparison stays essentially unchanged (see Table II).

Relative Noise . | Standard three-phase method . | Hilbert demodulation method . | SSOP method . | SSMD method . |
---|---|---|---|---|

10% | 28.20% | 35.01% | 30.76% | 16.62% |

5% | 14.10% | 19.65% | 18.55% | 14.52% |

1% | 2.82% | 10.83% | 12.31% | 13.78% |

Relative Noise . | Standard three-phase method . | Hilbert demodulation method . | SSOP method . | SSMD method . |
---|---|---|---|---|

10% | 28.20% | 35.01% | 30.76% | 16.62% |

5% | 14.10% | 19.65% | 18.55% | 14.52% |

1% | 2.82% | 10.83% | 12.31% | 13.78% |

Relative Noise . | Standard three-phase method . | Hilbert demodulation method . | SSOP method . | SSMD method . |
---|---|---|---|---|

10% | 31.03% | 39.25% | 38.41% | 14.04% |

5% | 15.51% | 21.80% | 21.55% | 13.94% |

1% | 3.10% | 11.51% | 11.85% | 13.91% |

Relative Noise . | Standard three-phase method . | Hilbert demodulation method . | SSOP method . | SSMD method . |
---|---|---|---|---|

10% | 31.03% | 39.25% | 38.41% | 14.04% |

5% | 15.51% | 21.80% | 21.55% | 13.94% |

1% | 3.10% | 11.51% | 11.85% | 13.91% |

We further simulate a case consisting of two AC and one DC components as shown in Fig. 4(a) where the “中” AC image is modulated vertically with a period of *T*_{y} =5 pixels and the “华” AC image is modulated diagonally with a period of *T*_{x} = *T*_{y} =10 pixels. SSMD and Hilbert demodulation methods are used to extract the two AC components. SSMD successfully recovers all AC components yet Hilbert demodulation method fails as it is designed for only one single AC component (see Fig. 4(b–d)).

### B. Homogenous phantom experiment

The performance of SSMD is then compared with the standard three-phase demodulation method on their accuracy of retrieving the modulation transfer function of an Intralipid-2% suspension (Intralipid-10% diluted with de-ionized water with the volume ratio 1:4) in a quartz container. First, the standard three-phase demodulation is applied to measurements (a total 18 images) at the following single modulation frequencies 0.487, 0.689, 0.975, 1.706, 2.413, and 3.412 (mm^{-1}). Second, SSMD is applied to measurements of the same single modulation frequencies, those modulated at two frequencies simultaneously (0.487 and 0.689 mm^{-1}; 0.487 and 0.975 mm^{-1}; 0.689 and 0.975 mm^{-1}; 1.706 and 2.413 mm^{-1}; 1.706 and 3.412 mm^{-1}; and 2.413 and 3.412 mm^{-1}), and those modulated at three frequencies simultaneously (0.487, 0.689 and 0.975 mm^{-1}; and 1.706, 2.413 and 3.412 mm^{-1}). The recovered values and their standard deviations for MTF are displayed in Fig. 5. There is a slight elevation of the MTF at the highest modulation frequency, in particular, for the three-phase method due to the degrade of the signal to noise ratio (SNR) in measurements at higher frequencies. Fig. 5 also shows that SSMD exhibits smaller standard deviations than the standard three-phase method, similar to what being observed in numerical simulations (Tables I and II) under moderate to high noise levels.

The 2D maps of the optical properties of the Intralipid-2% suspension reconstructed by the three-phase method and SSMD are then compared in Fig. 6. The incident beam is spatially modulated at *f*_{1} = 0.487mm^{-1} vertically and/or *f*_{2} = 0.689mm^{-1} diagonally. The optical properties of the medium are reconstructed by a look-up table^{1–3} based on the enhanced diffusion approximation.^{18} The three-phase method is applied to measurements modulated at one single frequency alone and SSMD to measurements modulated at either single or two frequencies simultaneously. Both methods recover the optical properties of the Intralipid-2% well with the absorption and reduced scattering coefficients of 3.10 × 10^{-4}mm^{-1} and 2.63mm^{-1} from *f*_{1} demodulation and 3.69 × 10^{-4}mm^{-1} and 3.14mm^{-1} from *f*_{2} demodulation, respectively. The same suspension was measured by an oblique incidence diffuse reflectance method^{19} which yields $\mu a$=3.0 × 10^{-4} mm^{-1} and $\mu s\u2032$ = 2.60 mm^{-1}. The recovered $\mu a$ and $\mu s\u2032$ at *f*_{2} are slightly overestimated compared to the expected values partly attributed to the gradual loss of the accuracy of the diffusion approximation at higher modulation frequencies.^{20} The optical property maps produced with SSMD overall show much less fluctuations comparing to those with the standard three-phase method. The performance of SSMD also holds up well when modulating at either single or multiple frequencies simultaneously. Note the intensity of the reflectance inside the central region is stronger than the surrounding region owing to the boundary effects. The boundary effects also show up clearly in the recovered optical properties maps, in particular, nearby the left edge in the *f*_{1} case with SSMD.

To validate the experimental system and data analysis procedure, we have also simulated the SFDI measurement and recovered the AC component with the standard three-phase demodulation method for an incident structured-illumination identical to the experimental condition upon a semi-infinite medium ($\mu a$ = 3.0 × 10^{-4} mm^{-1} and $\mu s\u2032$= 2.60 mm^{-1}) for the case of *f*_{1}=0.487mm^{-1}. Fig. 7 compares the experimental and simulation results. Figure 7(a, c) are the measured and simulated light reflectance images, respectively. Their horizontal profiles (dash lines) are shown in (e). The two horizontal profiles agree very well with each other. Figure 7(b, d) are the recovered AC MTF from the experimental and simulation data. Their horizontal profiles (dash lines) are shown in (f). Due to the noise and imperfections in the intensity modulation of the incident light, the experimental three-phase method data exhibits more oscillation than simulation. Even so, the two profiles still coincide well with each other. The stripes observable in the recovered AC component by the standard three-phase method on the experimental data originate from the boundary effect and persist in the simulation result as well.

### C. In vivo back of hand experiment

Lastly, spatially modulated light was used to image the back of one hand of a healthy male. The spatial modulation frequency is set at 3.412mm^{-1} vertically. Three images with phases 0, $2\pi /3$, and $4\pi /3$ were taken in sequence. The standard three-phase method used all three images to recover the MTF whereas SSOP and SSMD used only the first image to compute MTFs. We used an analytical light reflectance model incorporating both diffuse and non-diffuse photons to reconstruct the optical properties from measured MTFs.^{20} The anisotropy factor was assumed to be 0.87.^{21,22} The reconstructed optical properties of the back of the hand over a square of size 10 × 10 mm is shown in Fig. 8(a). The magnified region is shown in Fig. 8(b) for the highlighted region of size 1 × 1 mm corresponding to one hair pore. The white dash lines highlight a vein which exhibits higher absorption and lower reduced scattering. The optical property map generated from the standard three-phase method suffers from motion artefacts, which renders an incorrect optical property map. In contrary, as both SSOP and SSMD only used one single frame of image to retrieve all optical properties, motion artefacts were significantly reduced in their reconstructed optical properties. Furthermore, the reconstructed optical property maps with SSMD were found to be more accurate than those with SSOP and were able to show hair pores on the back of hand (see Fig. 8). The average absorption and reduced scattering coefficients recovered by SSMD over the whole imaged window is 0.019mm^{-1} and 1.64mm^{-1}, respectively, within the range of the reported optical properties of skin.^{21,22}

## IV. DISCUSSION

In this work, we described and validated the Single Snapshot Multiple frequency Demodulation method applied to the Spatial Frequency Domain Imaging. SSMD can extract the DC and multiple AC components of different orientation, modulating frequencies and amplitudes from a single snapshot at once. Comparing to the conventional three-phase, Hilbert demodulation methods, or the SSOP method, in terms of noise suppression the standard three-phase method performs best at a low noise level whereas SSMD achieves the best performance at a moderate to high noise levels (see Fig. 3 and Tables I and II). SSMD effectively removes the oscillations in the reconstructed optical property maps when the three-phase demodulation is employed (see Fig. 6). This robustness of SSMD originates from its construction (Eqs. (3)–(6)). Even for the case of only one AC component of a specific spatial frequency *f*_{1}, there are always noise and imperfections in the modulation pattern of the incident beam that contributing to the error in the retrieval of the MTF at that frequency *f*_{1}. Noting that noise is wide-band and the imperfections in the modulation pattern contain higher order harmonics of the base frequency, SSMD’s retrieval formalism (Eq. (5)) locks in at the modulating frequency *f*_{1} alone and effectively suppresses noise and imperfections in the modulating pattern. This gives rise to the very desirable behavior of SSMD that the relative error in the recovered AC components stays relatively flat even with increasing noise, unlike other competing demodulation approaches (see Tables I and II). This also explains why SSMD obtains more accurate MTF at higher spatial modulating frequencies than the standard three-phase method (see Fig. 5).

To understand the noise suppression capability of SSMD comparing to the standard three-phase approach, consider a set of reflectance measurement given by Eq. (1) with iid noise of magnitude ε superimposed at each pixel. The extracted AC component by Eq. (2) yields a variance of $\u27e8\Delta IAC2\u27e9=23\epsilon 2$ with the standard three-phase approach whereas it is given by $\u27e8\Delta IAC2\u27e9=2T1T2\epsilon 2$ via Eq. (5) in SSMD of a kernel of size *T*_{1} × *T*_{2} when max(*T*_{1},*T*_{2})>2. The error in the AC component recovered by SSMD due to noise hence decreases proportional to the square root of the number of pixels contained in the kernel window.

We would like to point out the superior noise suppression of SSMD comes with a trade-off. The amplitude of the AC component is obtained by filtering the reflectance image with a kernel (see Eq. (5)). This requires that the optical properties of the medium vary slowly and can be regarded unchanged over the kernel window. Such an assumption is also made in SFDI in general by assuming the local medium to be homogeneous to extract the optical properties. When the kernel size is larger than the scale over which the optical properties of the medium vary, SSMD may lead to a sacrifice in the spatial resolution and introduce error. One notable example is Fig 3 and Table I where the error introduced by the insufficient resolving power of the sharp edges in “华” with SSMD degrades its performance, losing to the standard three-phase method at the low noise level. SSMD is thus particularly suited for un-mixing modulations at moderate or high spatial frequencies where the above assumption is valid. It is also worth pointing out that the kernel window is not always the smaller the better and SNR is also an essential factor in imaging. The relative error in the determination of the AC component with SSMD reduces to $3/T1T2$ times of that with the standard three-phase method. This noise reduction is especially relevant for dynamic systems where repeated measurements are impossible. SSMD hence expands the range over which MTF can be correctly measured than that with the current demodulation approaches for, in particular, real-time SFDI applications.

Another major advantage of SSMD is that it can simultaneously resolve multiple AC components of different frequencies from a single snapshot at once. The orientation, frequency and amplitude of the multiple AC components can be chosen freely to optimize the signal to noise performance over the medium response at all the modulation frequencies of interest. This reduces motion artefacts and is critical for real time imaging of dynamic turbid media. To eliminate entirely motion artefacts, the single image acquisition speed (exposure time) should still be much faster than the motion of the sample. For example, the standard three-phase demodulation method fails to generate a correct optical property map of the back of the hand due to motion artefacts as the standard three-phase demodulation method requires at least three images to be taken (see Fig. 8). On the other hand, by reducing motion artefacts, SSMD clearly resolves structures of the back of the hand. Moreover, as the probing depth of SFDI varies with the modulation frequency, the plurality of the modulation frequencies within a single snapshot in SSMD enables the characterization of the optical properties of the turbid medium at different depths simultaneously at once as well.

## V. CONCLUSION

We have presented a novel Single Snapshot Multiple frequency Demodulation (SSMD) method enabling single snapshot wide field imaging of optical properties of turbid media in the Spatial Frequency Domain. SSMD makes use of the orthogonality of harmonic functions and extracts the modulation transfer function (MTF) at multiple modulation frequencies and of arbitrary orientations and amplitudes subject to the implementation device simultaneously from a single structured-illuminated image at once. The excellent performance of SSMD has been demonstrated with experiment on both tissue mimicking phantoms and *in vivo* for recovering optical properties by comparing to the conventional three-phase approach. SSMD not only increases significantly the data acquisition speed and reduces motion artefacts but also exhibits excellent noise suppression in imaging as well. SSMD is ideal in the implementation of a real-time spatial frequency domain imaging platform, which will open up SFDI for vast applications in, for example, mapping the optical properties of a dynamic turbid medium and monitoring fast temporal evolutions.

## ACKNOWLEDGMENTS

This work is supported by National Nature Science Foundation of China Award# 81470081, Key Research Program of Zhejiang Natural Science Foundation of China Award# Z16H180005, and US DOD Award# W81XWH-10-1-0526. We also thank Wenlei Yu for the help in the oblique reflectometry measurements.