In vivo two-dimensional quantitative imaging of skin and cutaneous microcirculation with perturbative spatial frequency domain imaging (p-SFDI).

We introduce perturbative spatial frequency domain imaging (p-SFDI) for fast two-dimensional (2D) mapping of the optical properties and physiological characteristics of skin and cutaneous microcirculation using spatially modulated visible light. Compared to the traditional methods for recovering 2D maps through a pixel-by-pixel inversion, p-SFDI significantly shortens parameter retrieval time, largely avoids the random fitting errors caused by measurement noise, and enhances the image reconstruction quality. The efficacy of p-SFDI is demonstrated by in vivo imaging forearm of one healthy subject, recovering the 2D spatial distribution of cutaneous hemoglobin concentration, oxygen saturation, scattering properties, the melanin content, and the epidermal thickness over a large field of view. Furthermore, the temporal and spatial variations in physiological parameters under the forearm reactive hyperemia protocol are revealed, showing its applications in monitoring temporal and spatial dynamics.


Introduction
Biomedical optical imaging and diagnostic methods have been actively pursued in past decades due to some significant advantages, including non-invasiveness, high sensitivity to the chromophore content, and low cost. For skin imaging, optical coherence tomography (OCT) [1] and confocal microscopy (CM) [2] have emerged as high-resolution (e.g., 1-10µm) imaging methods. Their main drawback is the limitation to morphology and imaging depth. Spectroscopic techniques such as diffuse reflectance spectroscopy (DRS) [3], near-infrared spectroscopy (NIRS) [4], and Raman spectroscopy (RS) [5] have also been used to acquire information about the structure and biochemical composition of the tissue. Raster scanning is required for wide-field two-dimensional (2D) imaging and often not affordable for real-time diagnosis.
Recently, Spatial Frequency Domain Imaging (SFDI) [6][7][8][9][10][11][12][13] has attracted significant attention as one noncontact and wide-field quantitative modality for rapid mapping of optical properties (absorption and scattering coefficients) and physiological parameters (hemoglobin concentration, oxygen saturation, and the melanin content) of tissue. The optical and physiological properties are obtained by fitting the tissue modulation transfer function (MTF) to appropriate light migration models [14][15][16] or Monte Carlo simulations [17]. A look-up table has been utilized to approximate the relationship between the optical properties and the measured MTF values and speed up the two-dimensional imaging using SFDI [18]. In more complex situations such as imaging layered media, SFDI is limited in speed for a video-rate two-dimensional large field of view mapping with the traditional methods through a pixel-by-pixel inversion. As one potential solution, deep neural networks were proposed for optical parameter estimation with spatially resolved reflectance [19][20][21], reaching 300−100,000 times improved inversion speed. However, to the best of our knowledge, there have not been reports of a deep neural network inverse model that solves more than two parameters (absorption and reduced scattering coefficients) in SFDI. Therefore, new fast processing and interpretation of the two-dimensional SFDI measurement for the simultaneous mapping of the optical properties and physiological characteristics are desired, particularly for in vivo clinical applications.
In this article, we, for the first time, introduce perturbative spatial frequency domain imaging (p-SFDI) for fast two-dimensional mapping of the optical properties and physiological characteristics of skin and cutaneous microcirculation using spatially modulated visible light over a large fieldof-view (FOV). Our approach significantly shortens parameter retrieval time with perturbation analysis based on an analytical model of structured light reflection by layered media. More importantly, it substantially enhances the image reconstruction quality and largely avoids the random fitting errors caused by measurement noise inherent to the traditional methods through a pixel-by-pixel inversion. The efficacy of p-SFDI is then demonstrated by in vivo imaging the forearm of one healthy subject and monitoring the temporal and spatial cutaneous hemodynamics under the forearm reactive hyperemia protocol.

Spatial frequency domain imaging for the two-layer tissue
Spatial frequency domain imaging (SFDI) is a wide-field and noncontact method that uses spatially modulated structured illumination patterns: where f x, φ are the spatial modulation frequency and the phase angle. The reflectance takes the form of The ratios I DC /I (0) DC and I AC /I (0) AC represent the modulation transfer function (MTF) of reflected light at the spatial frequency 0 and f x . Three-phase demodulation from three measurements taken at φ, φ+2π/3, and φ+4π/3 phase delays are commonly used to compute MTF for each spatial frequency. The MTF values are then fitted to a light reflectance model to obtain absorption (µ a ) and reduced scattering (µ s ′ ) coefficients pixel-by-pixel in conventional SFDI [8,22]. Light reflectance from forward-peaked scattering media such as biological tissue at an arbitrary separation and any spatial frequency has been described in detail [15]. At low spatial frequencies, light reflection is diffusive. Diffuse light reflectance is dominated by the snake and diffuse photons and is given by and from spatially modulated incident light of unit intensity where q ≡ 2πf x , β ≡ µ a + µ ′ s , Q ≡ √︁ q 2 + 3µ a β, and l ≡ α/µ ′ s . Here q, β, l and α are the angular frequency, reduced attenuation coefficient, extrapolation length, and surface roughness [23], respectively. Tissue MTF measurements at a low spatial frequency can be used to determine optical parameters (µ a , µ s ′ , and α) by fitting to their theoretical prediction of the enhanced diffuse light reflectance I = I snake +I diffuse .
Biological tissue such as skin is highly structured and layered. Chromophores like melanin are typically confined within the epidermis, whereas hemoglobin is confined within the dermis beneath. Accounting for the two-layer structure is critical for accurately recovering hemoglobin concentration and oxygen saturation [7,24,25]. We have developed a novel two-layer model for SFDI, mapping the two-layer structure to an equivalent homogeneous medium [11]. The absorption coefficients within the epidermis and dermis are given by µ a,epidermis (λ) = ε melanin (λ)c melanin (4) and where ε Hbo2 , ε Hb , ε melanin are the molar extinction coefficients of oxygenated hemoglobin, deoxygenated hemoglobin, and melanin, respectively; and c THb , StO 2 , c melanin are the total hemoglobin concentration, oxygen saturation, and the melanin concentration, respectively. The equivalent absorption coefficient µ a (q, λ) is determined by requiring total light absorption to be equivalent within the two systems, i.e., Here h is the epidermal thickness, and L is the mean penetration depth for the modulated light at the spatial frequency f x and the wavelength λ [11,26], given by As light scattering by the epidermis and dermis is dominated by the fractal refractive index fluctuation [16,27,28], the reduced scattering coefficient can be written by a power law ignoring the difference in the scattering properties between the epidermis and dermis [7,11,24,25] where b is the scattering power. The reference wavelength λ 0 in Eq. (8) is set to the central wavelength in our experimental setup (540 nm). The pixel-by-pixel inversion for layered media [11,13] fits c THb , StO 2 , c melanin , h, µ s ′ (λ=540nm) , b, and α using the function "fmincon" in Matlab by minimizing the least squared error: at each point where i=1, 2, 3 represent the three wavelengths (Red (623 nm), Green (540 nm), and Blue (460 nm)) and the theoretical values of the modulation transfer functions mtf DC and mtf AC are computed with the enhanced diffuse light reflectance I = I snake +I diffuse at the spatial frequencies 0 and f x . Validation on the standard optical phantom (Biomimic, Canada) shows an error of less than 2% for the measured absorption and scattering coefficients comparing to their known values [11].

Fast two-dimensional mapping of parameters based on perturbative SFDI
In contrast to the traditional methods for the recovery of 2D maps through a pixel-by-pixel inversion, the perturbative SFDI comprises two steps: first grouping pixels into multiple clusters according to their measured values and recovering the optical and physiological properties for an "average" case for each cluster and then solving all pixels in that cluster through perturbation to the "average" one. We used k-means clustering [29][30][31] to group the pixels over the whole field of view into different groups according to their measured DC and AC MTF values.
The theoretical prediction of diffuse light reflectance I = I snake +I diffuse is a multivariable function of c THb , StO 2 , c melanin , h, µ s ′ (λ=540nm) , b, and α. It can be rewritten as to the first order where X is the vector of variable parameters (c THb , StO 2 , c melanin , h, µ s ′ , b, and α), and X 0 is the recovered optical and physiological properties for the "average" case whose MTFs are the average over one whole cluster. In our setup, DC and AC light reflectance at three wavelengths λ 1 =460 nm, λ 2 =540 nm, and λ 3 =623 nm are measured. Equation (10) yields (11) and can be simply expressed as ∆I = S n ∆X n (12) in a matrix notation. With the analytical light reflectance expressions, the partial derivative S n can be computed straightforwardly. We have found (see Sec. 3.1) that light reflectance depends close to linearly on each variable over the typical range for biological tissues, leading to a stable S n . The flowchart for fast two-dimensional mapping of parameters is outlined in Fig. 1. First, we extract tissue DC and AC MTFs at the multiple wavelengths with three-phase demodulation and group the pixels over the whole field of view into different clusters according to their MTF values using k-means clustering. The average of MTF values over all pixels in each cluster is computed to yield the "average" case <I>. The "average" case is fitted to obtain the "average" optical and physiological properties X 0 and the partial derivative S n at X 0 for each cluster. The perturbation to each pixel's optical and physiological properties is then given by the product of corresponding partial derivative inverse S n −1 and light reflectance perturbation ∆I, yielding their correct optical and physiological properties X 0 + ∆X n . Compared to the traditional pixel-by-pixel inversion, the speed of two-dimensional mapping of optical properties and physiological parameters based on perturbation SFDI increases approximately by a factor equal to the total number of pixels in one cluster. More importantly, this approach suppresses the sporadic fitting errors caused by measurement noise and substantially enhances the image reconstruction quality.

Experiment procedure
The schematic diagram for the digital micro-mirror device (DMD)-based SFDI system is shown in Fig. 2. The illumination comprises three monochromatic light of three wavelengths: 623 nm (Red), 540 nm (Green), and 460 nm (Blue), respectively. The sinusoidal fringe pattern produced on a DMD (DLP LightCrafter TM 4500, Texas Instruments) is magnified by a factor of four and projected onto the specimen. The spatial modulation frequency in our reported experiments is 0.2mm −1 . A pellicle beamsplitter is used to integrate the modulated illumination system and imaging system without introducing ghosting. A CCD camera (Point Gray Grasshop3 GS3-U3-51S5C) was used to record backscattered light reflectance images. The exposure time was set at 20 ms. The SNR for the recorded intensity at one single pixel is approximately 28 dB. The field of view on the specimen surface is 27mm×18 mm, and one millimeter on the specimen surface corresponds to 60 pixels on the CCD in the reported experiments. As the DMD's output light intensity is not linear and uniform over the whole field of view, calibration was performed to correct the nonlinearity and nonuniformity based on the diffuse reflectance from a Lambertian reflectance standard (Hefei Xingyue Luminous Technology Applications Institute, China) [10]. The modulation transfer function of the optical system was also calibrated by imaging resolution test targets (R1L1S1P, Thorlabs, Inc.) over the spatial frequency range. The diffuse reflectance from the Lambertian reflectance standard was used to determine the illumination light intensity. Color correction was further performed to unmix the three channels of the camera [11]. The absolute MTF of the specimen is computed after incorporating all the above factors [10,12].

Simulation study for fast two-dimensional mapping of optical and physiological parameters
Two typical types of skin tissues (normal skin and melanotic nevus) were used to simulate and verify the robustness of perturbative SFDI. The typical physiological and optical parameters (<c THb >, <StO 2 >, <c melanin >, <µ s ′ >, <b>, <α>, <h>) of normal skin and melanotic nevus were set at (0.02 mM, 0.7, 4.5 mM, 1.6mm −1 , 1, 0.7, 0.1 mm) and (0.01 mM, 0.5, 8.0 mM, 1.0mm −1 , 0.5, 0.4, 0.2 mm), respectively, according to our earlier findings [11]. The reduced scattering coefficient (µ s ′ ) was for the central wavelength 540 nm. Figure 3 shows the dependence of DC and AC light reflectance at the three wavelengths on each parameter within the human tissue range for normal skin and melanotic nevus tissues. The dependence is approximately linear, resulting in a stable slope S n . The calculated sensitivity of light reflectance on seven parameters at the reference points for two different skin types is summarized in Table 1.  To evaluate the image quality improvement with p-SFDI against the pixel-by-pixel least-squares fitting method, we simulated the AC and DC MTFs at three wavelengths (460 nm, 540 nm, 623 nm) and spatial frequencies f = 0 and 0.2mm −1 for a layered medium of properties similar to a nevus based on the enhanced diffusion model. The accuracies of the recovered optical and physiological parameters by p-SFDI and pixel-by-pixel inversion were compared under different noise levels (see Fig. 4 and Table 2). P-SFDI outperforms the pixel-by-pixel inversion with less mean squared relative errors, particularly at a higher noise level. Even at a low noise level (1% noise), p-SFDI improves image reconstruction quality significantly, removing image artifacts and outlining the nevus morphology more accurately by suppressing the adverse effects of measurement noise (see Fig. 4).

In vivo two-dimensional fast imaging of forearm skin parameters of one healthy subject
Forearm skin of one healthy subject was demonstrated for the 2D fast mapping of the skin optical properties and physiological characteristics. The c THb , StO 2 , melanin area-density concentration (c melanin ×h), µ s ′ , b, α and h maps were recovered, respectively, by perturbative SFDI (see Fig. 5(a)) and pixel-by-pixel least-squares fitting (see Fig. 5(b)). The 100pixel×100pixel nevus region of interest (ROI) was selected in the image window, and the parameter k was set to 20 for   the k-means clustering for data analysis. The mean and standard deviation of the seven parameters retrieved by the two different methods was displayed in Fig. 6 and Table 3. The results show that the average optical and physiological parameters retrieved by the two methods agree. However, perturbative SFDI performs significantly better in recovering the nevus' morphology thanks to its superior measurement noise suppression (see Fig. 5). The appreciable noise interference at the upper right corner of the ROI (see Fig. 5(b)) largely disappears in Fig. 5(a). The nevus region is observed to have much higher melanin content and lower oxygen saturation than surrounding areas, agreeing with the findings by Ref. [32].

Monitoring cutaneous hemodynamics
The cutaneous hemodynamics for one healthy subject under the forearm reactive hyperemia protocol was then imaged by perturbative SFDI. The measurement was taken for 4 minutes during occlusion (cuff pressure: 200 mmHg) and 3 minutes after cuff release at the rate of 3 images per second. The photo of the imaged area and 2D maps for optical and physiological properties of the human forearm skin at one minute before and after cuff release were shown in Fig. 7. Box and whisker plots were shown in Fig. 8 for hemoglobin concentration and oxygen saturation (StO 2 ) in the venous region. The average and standard deviation of c THb and StO 2 were 0.0165 ± 0.0009 mM, 0.578 ± 0.043 upon cuff occlusion, whereas they jumped to 0.0245 ± 0.0009 mM, 0.835 ± 0.071 upon the cuff release, respectively. The other parameters remained essentially unchanged.   More interestingly, the spatial distribution of the total hemoglobin also expands, suggesting the recruitment of idle capillaries upon the cuff release [33].

Discussion and summary
We have presented perturbative SFDI (p-SFDI) for fast two-dimensional mapping of the optical properties and physiological characteristics for skin and cutaneous microcirculation using spatially modulated visible light. Various methods have been proposed to obtain the DC and AC amplitudes from demodulating structured light reflectance [9,10,34,35]. One crucial trade-off in demodulation is between spatial resolution and noise suppression. The three-phase demodulation for SFDI enjoys the highest spatial resolution but is prone to measurement noise [10]. Other methods rely on different forms of filtering. The single snapshot multiple frequency demodulation (SSMD) extracts DC and AC components from a single snapshot at a √ N times higher signal-to-noise ratio (SNR), where N is the total number of pixels inside the SSMD kernel [10]. However, SSMD comes at the cost of reduced spatial resolution when the scale over which the medium's optical properties vary is smaller than the kernel size.
The proposed p-SFDI provides a unique approach to improve the robustness and SNR of the three-phase demodulation while retaining its superior spatial resolution. The key of p-SFDI is by grouping pixels into multiple clusters by k-means clustering according to their measured MTF values and recovering the optical and physiological parameters for an "average" case, followed by solving all pixels in one cluster through perturbation to the "average" one. The validity of the perturbation approach is ensured by maintaining sensitivities within each cluster (covering ∼10% of the AC and DC range; see Fig. 3) to be approximately constant. The DC and AC MTF values of the "average" case are the average of respective MTF values for all pixels within the cluster, whose SNR increases by roughly √ N ′ times with N' being the total number of pixels in one cluster. A well-posed matrix inversion then finds the perturbation to the "average" case and recovers the true values. As a consequence of this two-step approach, the proposed p-SFDI suppresses measurement noise and enhances the image reconstruction quality without losing three-phase demodulation's spatial resolution.
Compared to the traditional pixel-by-pixel inversion, p-SFDI not only largely avoids the image artifacts caused by measurement noise but also significantly shortens parameter retrieval time. The image retrieval speed is much faster by p-SFDI than by pixel-by-pixel least-squares fitting. We used a non-linear fitting routine "fmincon" in Matlab (MathWorks Inc., USA) to fit the optical and physiological parameters from the measurement. The fitting routine is iterative and time-consuming. The speedup of p-SFDI over the pixel-by-pixel inversion is approximately N' times. For the images presented in Fig. 4, the pixel-by-pixel fitting method took nearly over 5 hours, while p-SFDI took only 40s. Varying the number of clusters used in p-SFDI has negligible effects on the recovery of optical properties and physiological characteristics. The difference is about 1% when the number of clusters varies between 20 and 500. A cluster number of ∼20 works well for biological tissue of a limited range for their properties to achieve image reconstruction accuracy and speed requirements.
We note the p-SFDI processing time scales with the number of clusters but not the image size. For a full-sized image of 1620×1080 pixels, the inverse process remained approximately 40s with p-SFDI on a PC desktop with 8 GB of RAM and an AMD R7-1700 eight-core processor, excluding the processing time required by k-means clustering (k=20). When a database of the "average" case is pre-computed, p-SFDI will enable 2D imaging optical and physiological properties over a large field of view in real-time using the fast k-means clustering algorithms [36,37].
The efficacy of this perturbative SFDI has been demonstrated by in vivo imaging forearm of one healthy subject, yielding accurate determination of cutaneous hemoglobin concentration, oxygen saturation, the scattering properties, the melanin content, and the epidermal thickness over a large field of view in real-time. The temporal and spatial variations in physiological parameters under the forearm reactive hyperemia protocol have also been revealed. This work paved the way for 2D imaging and monitoring temporal and spatial dynamics in vivo over a large field of view in real-time with perturbative SFDI.