Open Access
Add to BibSonomy
Abstract
We have proposed the profile-based intensity and frequency corrections for single-snapshot spatial frequency domain (SFD) imaging to mitigate surface profile effects on the measured intensity and spatial frequency in extracting the optical properties. In the scheme, the spatially modulated frequency of the projected sinusoidal pattern is adaptively adjusted according to the sample surface profile, reducing distortions of the modulation amplitude in the single-snapshot demodulation and errors in the optical property extraction. The profile effects on both the measured intensities of light incident onto and reflected from the sample are then compensated using Minnaert’s correction to obtain the true diffuse reflectance of the sample. We have validated the method by phantom experiments using a highly sensitive SFD imaging system based on the single-pixel photon-counting detection and assessed error reductions in extracting the absorption and reduced scattering coefficients by an average of 40% and 10%, respectively. Further, an in vivo topography experiment of the opisthenar vessels has demonstrated its clinical feasibility.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Spatial frequency domain (SFD) imaging is a noncontact wide field imaging modality that quantitively maps the optical properties of biological tissues, including the absorption coefficient (${\mu _a}$) and reduced scattering coefficient (${\mu ^{\prime}_s}$), at multiple wavelengths, from which the biochemical component distributions are then extracted over the sample surface [1,2]. In recent years, development of the single snapshot imaging has propelled a considerable speed-up of the amplitude demodulation that further improves the real-time capability of the technique [3–5]. The real-time implementation will render SFD imaging a potential tool in many clinical applications, e.g., blood oxygen monitoring, burn degree assessment, and fluorescence-guided surgery [1–13].
As a specific wide field imaging modality of structural illumination that employs space harmonic modulation of incident light intensity, SFD imaging has been challenged in realistic scenarios by how to mitigate the surface profile effects on both the spatial frequency and measured intensity. This issue arises from the fact that the measured intensity and spatial frequency of the reflected light from an uneven surface can be deviated from those from a planar one. Thus, direct application of the established methods in the optical property extraction without correcting the above deviation definitely exerts adverse influences on the performances of the methods. Several works have been devoted to elimination of the surface profile effects. Nguyen et. al. has proposed a scheme of phantom calibration where a reference phantom having the identical shape with the subject was fabricated using three-dimensional (3D) printing technique, and the light distributions on the same-shaped surfaces of the subject and reference phantom were measured to eliminate the profile effects on the calculation of the true diffuse reflectance [14,15]. Gioux et. al. has introduced an intensity correction method by compensating the measured reflected intensity using Lambertian model, and described a multi-frequency processing algorithm by interpolating the optical properties of the sample extracted from pre-measured multiple frequencies to solve the spatial frequency deviations due to the curvature of the sample [16]. On the basis of Gioux’s work, Zhao et. al. has introduced a modified angle correction method using Minnaert’s model to prevent over correction of a steep surface [17].
Despite of their promising results, there are some drawbacks with these approaches. Although the phantom calibration is universal for the demodulation modes, it is apparently impractical due to difficulties in fabricating a profile-adapted phantom for each subject. The model-based corrections, despite of requiring no additional hardware modification and taking less pre-experiment time, have only compensated deviations of the intensity reflected from the sample to the camera, more involved correction is necessary for an imaging system with high angle between the projection and collection optical axes. The multi-frequency processing algorithm has disadvantages in size and storage, which makes it difficult to implement in a field-programmable gate array. Moreover, current approaches pay less attention to amplitude distortion of the single-snapshot demodulation caused by spatial frequency deviations. Presently, there are two kinds of single-snapshot demodulation methods; one isolates the planar (DC) and spatially modulated (AC) components using spectrum filtering to extract the DC and AC modulation amplitudes from Fourier inverse transformation and Hilbert demodulation [3,4], and the other uses the orthogonality of harmonic functions to extract DC and AC amplitudes by spatial filtering the image over a sliding window [5]. As a window-wise method, the single-snapshot demodulation is dependent of the adjacent pixels in space that requires the sinewave frequency to remain constant. The frequency broadening increases the difficulty of filtering DC and AC components from the spectrum and the errors of amplitude extraction from the envelope, thus distorting the results. Figure 1 compares the performance of the three-phase demodulation and the single-snapshot demodulation of a frequency deformed sinewave. The three-phase demodulation, as a pixel-wise demodulation method, extracts the modulation amplitude independently for every single pixel of the image [12]. Although the surface curvature distorts the frequency projected, the modulation amplitude is extracted accurately since the phase difference of each pixel between the sequentially projected frames remains unchanged. But the three-phase demodulation requires multiple projections that degrades the real-time capability of SFD imaging.
Fig. 1. Performances of the three-phase demodulation and single-snapshot demodulation based on spectrum filtering. The signal to be demodulated is a sinewave that varies continuously in frequency from 0.1 mm−1 to 0.15 mm−1.
To solve the problems above, we have proposed the profile-based intensity and frequency corrections for single-snapshot SFD imaging. In the scheme, spatial frequency deviations of the sinusoidal light on the sample is corrected by projecting a pre-distorted sinusoidal pattern, in which the spatially modulated frequency of each pixel is adjusted adaptively based on the surface profile of the sample. When the pre-distorted projection pattern is decoded on the sample, the sinusoidal light of uniformly corrected spatial frequency can be generated on the profiled surface. The scheme directly compensates the spatial frequency deviations of the sinusoidal light and thus reduces distortions of single-snapshot demodulation. Meanwhile, the optical properties of the sample can be simply extracted from a single measurement, requiring no additional calibration measurements and interpolation steps of multiple frequencies. After the frequency correction, the intensity deviations incident onto and reflected from the sample are simultaneously compensated using a calibration-based height correction as well as a Minnaert’s model-based incident and reflection angle correction.
The proposed scheme has been validated by phantom experiments on a highly sensitive SFD imaging system based on single-pixel photon-counting detection, incorporating the phase profilometry technique for 3D profile measurements. The results showed that the profile-adaptive projection corrected the spatial frequency deviations of the sinusoidal modulated light on the uneven surface with errors of <2% at each pixel, and consequently reduced distortions of the modulation amplitude in the single-snapshot demodulation by an average of 9%. The complete frequency and intensity corrections subsequently reduced the errors in recovering the absorption and reduced scattering coefficients by an average of 40% and 10%, respectively. Then an in vivo topography experiment of the opisthenar vessels further demonstrated the clinical feasibility of the method.
2. Methods
The optical properties of a sample (diffusive medium) are measured from single-snapshot SFD imaging, using the spectrum filtering method for amplitude demodulation. Heights of the sample are measured using a phase profilometry technique. The SFD imaging system based on single-pixel photon-counting detection is used for both the single-snapshot SFD imaging and the phase profilometry. Based on 3D surface maps, a pre-distorted projection pattern is designed for frequency correction of the sinusoidal light projected on the sample. The measured reflectance images are then compensated by the intensity correction.
2.1 Single-snapshot SFD imaging
SFD imaging projects spatially modulated sinusoidal light on a sample with a certain spatial frequency, f, in the x-direction as shown in Fig. 2, and extracts the absorption and reduced scattering coefficients of the sample based on the analysis of the measured diffuse reflectances. The intensity of the reflected sinusoidal light, I(x,f), at a single wavelength is formulated as
where MDC(x, f) is the DC component and MAC(x, f) is the AC modulation amplitude, and φ0 is the phase offset. The diffuse reflectance, Rd(x, f), is calculated from the following equation, where I0 is the intensity of the light source and ${F_{system}}({x,f} )$ is the system response (the modulation transfer function of the illumination and imaging optical system). Rd(x, f) is used to recover the absorption and reduced scattering coefficients of the sample through a two-dimensional look-up table of Rd(x, f) versus the optical properties, which is generated using a forward photon transport model [2,18].Fig. 2. Illustration of the SFD imaging system showing the arrangements of the projector, sample and camera with the x-y-z coordinate.
To improve the standard three-phase demodulation method, several fast demodulation methods based on spectrum filtering and spatial filtering have been proposed to extract the modulation amplitudes from a single SFD image. The spectrum filtering method first transforms the image to the Fourier domain and separates the spectrums into the DC and AC components using the low-pass and high-pass filters. The DC components are then obtained through the inverse Fourier transform, while the AC modulation amplitudes are extracted using both the inverse Fourier transform and the Hilbert demodulation [3,4]. The spatial filtering method uses a sinusoidal and cosine filtering kernel with the same modulation frequency to recover the AC modulation amplitude [5]. The size of the filtering window is determined by the spatial modulation period, and the DC component is extracted through a direct averaging over the same filtering window. However, a problem arises for both the spectrum and spatial filtering methods when targeting the tissue with the sinusoidal light deformed by the surface curvature.
Besides the demodulation distortions, the inconsistencies between the spatial frequency and intensity on the sample and those on the calibration phantom also produce errors in the optical property extraction. Typically, when a planar calibration phantom with known optical properties is used, the diffuse reflectances of the sample is normalized by Eq. (3) to eliminate the system response, ${F_{system}}({x,f} )$ in Eq. (2),
where ${M_{AC,ref}}(x,f)$ and ${R_d}_{,ref}(x,f)$ are the AC modulation amplitude and the diffuse reflectance of the calibration phantom, respectively. Equation (3) is accurate provided that the surface profiles of the sample and the calibration phantom are identical. For in vivo measurements of profiled tissues, the spatial frequency and intensity deviations have to be compensated.2.2 Phase profilometry for 3D profile measurements
The phase profilometry technique is applied to 3D profiling of the sample, using the SFD imaging system to project a set of sinusoidal fringes with a selected spatial frequency, ${f_p}$, onto the sample in the y-direction to calculate the surface heights, $h(x,y)$, from the measured phases, $\Phi (x,y)$ [19,20]. In general, a phase profilometer reconstructs the height profiles of the samples using the triangle similarity of the geometric structure of the system,
where H is the vertical height of the projector or camera, d is the horizontal distance between the projector and camera, and $\overline {\textrm{AB}}$ is one side of similar triangles. $\overline {\textrm{AB}}$ is obtained by $\Delta \Phi /2\pi {f_p}$, where $\Delta \Phi $ is the difference between the phase measured from sample surface and measured from a reference plane. Equation (4) can be rewritten as a height-phase function of $h = H/(1 + 2\pi {f_p}d/\Delta \Phi )$. This height-phase function requires the projector and camera to be at the same height. For more general situations, though, a height difference exists in the system. To solve the height difference, we described a modified height-phase function.As shown in Fig. 3, for point P on the sample, the side $\overline {SE} $ of the triangle $\varDelta SPE$ used to derive the height-phase function changes from d to $d - \Delta l$ as a result of the height difference between the projector and camera. The value of $\Delta l$ must be compensated for each pixel in the height-phase function to ensure the accuracy of the reconstructed height. Using the similar triangles $\varDelta \textrm{DEF}$ and $\varDelta \textrm{DAC}$, $\Delta l$ is calculated by $L \cdot ({H_2} - {H_1})/{H_2}$, where L is the distance between the pixel and acquisition center in the reference plane to be determined for each pixel, and H1 and H2 are the heights of the projection lens and the collection lens, respectively. The modified height-phase function is then formulated as
2.3 Spatial frequency correction based on the profile-adaptive projection
Generally, SFD imaging uses the same sinusoidal projection pattern for both the calibration phantom and the samples, regardless of the morphological differences between them. For the sample with an uneven surface, though, the projection pattern is deformed due to the surface profile variation and eventually generates the sinusoidal light on the surface with locally varied spatial frequencies [16].
Figure 4 illustrates the mechanism of how the projection pattern is deformed on an uneven surface when using the standard sinusoidal projection with the spatial frequency of f. When the SFD imaging system is setup, the modulated frequency of the standard sinusoidal pattern is determined by $f^{\prime} = f \cdot (d{x_{ref}}/dx^{\prime})$, where $d{x^\prime }$ and $d{x_{ref}}$ are the physical sizes of the pixels in the projection surface and the reference plane, respectively. For the calibration measurement, the sinusoidal pattern is linearly projected onto the flat phantom to generate the sinusoidal light of the desired spatial frequency, f, at each position. However, when the same pattern is projected onto an uneven surface, the physical size of the ith pixel on the surface, $d{x_i}$, varies due to the height variation, and the sinusoidal pattern is therefore nonlinearly projected on the surface. The actual spatial frequency at the ith pixel, fi, is given by the derivative of the phase with respect to the pixel position, $2\pi {f_i} = d{\varphi _i}/d{x_i}$. Since the phase at ith pixel, φ′, remains the same as φ, the spatial frequency at the ith pixel is calculated by
The frequency deviations could be corrected by the proposed profile-adaptive projection method, in which the sinusoidal pattern is pre-distorted according to the surface profile of each sample. To encode the pre-distortion and generate the sinusoidal light with the spatial frequency of f at each pixel on the sample, the modulation frequency of the sinusoidal pattern should be reversely adjusted at each pixel according to the distribution of $d{x_i}$. At the pixels with $d{x_i}$ larger than $d{x_{ref}}$ where the spatial frequency decreases, the modulation frequency of the projection pattern increases. Conversely, the modulation frequency decreases at the pixels with $d{x_i}$ smaller than $d{x_{ref}}$. For each line of the pre-distorted sinusoidal pattern, the modulation frequency at each pixel, ${f^{\prime}_{i,cor}}$, is calculated by
The physical size of ith pixel, $d{x_i}$, is determined from the measured 3D surface data of the sample. The coordinate system is given as shown in Fig. 3, where point O is the projection center on the reference plane, and S the optical center of the projector. The coordinate of each projection position on the reference plane can be determined within the coordinate system. The line from point S to the projection position constitutes the projection optical axis (the dashed line in the figure) of each pixel. Given the 3D coordinates of the surface, the intersection of each projection optical axis at the sample can be calculated. And the $d{x_i}$ is then obtained from the coordinates of adjacent intersections.
2.4 Incident and reflected intensity corrections
The surface profile of the sample mainly affects the light intensity in two aspects: 1) the intensity incident on the sample is affected by the surface heights and angles relative to the projector, 2) the measured intensity reflected from the sample is affected by the surface heights and angles relative to the camera. The intensity, following Lambertian model, is proportional to the cosine of the angle, and inversely proportional to the height. In addition to the surface profile of the sample, the effects are influenced by the geometry of the imaging system, including the position of the projector relative to the camera and their divergence angles. For instance, the projection optical axis is dependent of the projection divergence thus affecting the value of incident angle, and the angle between the projection and collection optical axes (${\beta _z}$) determines the difference in the incident and reflection angles. For small divergence and ${\beta _z}$ angles, the values of incident and reflection angles are close at each position on the sample, thus correcting one of the angles can be sufficient. For large divergence and ${\beta _z}$ angles, however, the intensities of light incident onto and reflected from the sample are required to be corrected simultaneously.
The height-dependent intensity deviation is compensated using the calibration-based height correction described in [16]. The calibration phantom is pre-measured at multiple heights, and the AC modulation amplitude for the calibration phantom, ${M_{ref}}(x,f)$, is extracted and interpolated to establish the coefficient matrix which described the modulation amplitude versus height at each pixel. For the sample of given height maps, ${M_{ref}}(x,f)$ is calibrated pixel-by-pixel to the same height as the sample.
The angle-dependent intensity deviation is compensated by the correction of two parts: incident intensity and reflection intensity. Reflection intensity correction using Minnaert’s model has been applied in prior work by Zhao et. al. [17], as shown in Eq. (8),
where $\theta$ is the reflection angle between the collection optical axis and the normal to the sample surface, as shown in Fig. 5. The collection optical axis is collimated (the divergence angle of the camera <10°) and perpendicular to the reference plane. $\theta$ is calculated for each pixel using the 3D surface data of the sample. ${k _1}$ refers to the Minnaert’s constant that prevents overcorrection of $\theta$. The value of ${k _1}$ varies from 0 to 1, depending on the nature of the sample, spatial frequencies and spectral wavelengths [21].Fig. 5. Normal to the sample surface (vn), incident angle ($\beta$) and reflection angle ($\theta$).
Similarly, the incident intensity on the sample is compensated using the Minnaert’s correction of the incident angle,
where ${I_0}(x,f)$ refers to the intensity incident on the calibration phantom and is obtained by ${M_{ref}}(x,f)\textrm{ }/{R_d}_{,ref}(x,f)$. $\beta$ is the incident angle between the projection optical axis and the normal to the sample surface, and the projection optical axis is inclined to the reference plane with the angle of ${\beta _z}$. Since the projection light is uncollimated (the divergence angle of the projector >40°), ${\beta _z}$ is calculated for each pixel.The diffuse reflectances of the sample are calculated from the corrected data using Eq. (3), formulated as
2.5 Experimental validation
The proposed scheme was validated using the previously developed SFD imaging system based on the single-pixel lock-in photon counting detection [22]. The lock-in photon-counting detection features high sensitivity and large dynamic range that enables the larger modulation depth of excitation light to effectively increase the signal-to-noise ratio of the AC components, and it enables the simultaneous measurement of multiple wavelengths and the rejection of ambient light. Sinusoidal fringe patterns for the measurements of the optical properties and of the 3D profile were projected onto the sample along and perpendicular to the acquisition plane, respectively. The spatial frequencies used for both of the measurements were 0.1 mm−1. The SFD reflectance images and phase images were measured using single-pixel imaging based on two-dimensional discrete cosine transform (DCT), which spatially compressed the pixel-array images using the sampling patterns composed of the DCT kernel matrices and recovered the images by applying an inverse DCT to the DCT coefficients acquired from the single-pixel detector [23,24]. In order to reduce the number of sampling patterns, the images are recovered by a selective acquisition of DCT coefficients concentrated in DC component and AC modulation frequency that fully utilizes the sparsity of sine harmonic in cosine base. The spatial frequency distribution of the selected DCT coefficients for recovering the SFD reflectance and phase images are shown in Fig. 6. Each DCT coefficient is obtained using the corresponding DCT kernel matrix. The number of the pixels of the recovered image is 128×128, with the imaging field of view of 40×40 mm2. The phases and the modulation amplitudes were extracted using 4-step phase shifting and the single-snapshot demodulation based on spectrum filtering, respectively.
Fig. 6. Spatial frequency distribution of the selected DCT coefficients for recovering SFD reflectance images (red box) and phase images (green box).
Tissue-like phantoms were fabricated for validation, using India ink as absorber, titanium dioxide (TiO2) as scatterer, and epoxy resin as substrate solvent [25]. The phantoms were then cured and made into different morphologies, including a homogenous flat phantom, a cylindrical phantom and a hemispherical phantom. By calculating the amount of scatterers and absorbers, the optical properties of the phantoms were determined. The flat phantom was vertically moved and tilted to simulate the changes in the surface height and angle, respectively. Based on the measured surface profile, the performance of the profile-adaptive projection to correct the spatial frequency and intensity deviations for reducing the errors in the extracted optical properties was assessed. The homogenous cylindrical and hemispherical phantoms were then used to verify the performance of the frequency and intensity corrections for more complex samples. In addition to the phantom experiments, an in vivo SFD topography of opisthenar vessels was obtained to further demonstrate the clinical feasibility of the proposed method.
3. Results
3.1 Profile corrections of the flat phantom
At first, the flat phantom was vertically translated from 0 to 30 mm in the z-direction relative to the reference plane with a step of 10 mm. The 3D surface map was measured at each height, as shown in Fig. 7(a). The reconstructed surface heights matched well with the expected values with the error of <1 mm for all the pixels, demonstrating the effectiveness of the modified height-phase function in the surface reconstruction. The height error originates from the measurement noise, and systematic errors associated with geometry and phase extraction. For subsequent corrections, this error mainly reduces the accuracy of calculated incident and reflection angles, but the resulting error of angle cosine is < 5% for all the pixels. The frequency and height corrections are less affected by this height error since it is much smaller than the heights of projector and camera (< 1%). Based on the measured height, the profile-adaptive projection pattern was generated. The projection patterns with and without the frequency correction were projected onto the surface, respectively, and the reflected sinusoidal light patterns were measured and compared with that measured for the reference plane. The results are shown in Fig. 7(a). Before the frequency correction, the spatial frequency of the sinusoidal light increased with the height by more than 11% per centimeter. The profile-adaptive projection corrected the frequency deviations effectively at each height, with the mean frequency error of the sinusoidal light < 2%.
Fig. 7. Frequency correction of the flat phantom. (a) The surface maps (left), the intensity profiles at height of 30 mm (center) and mean spatial frequency versus height (right) of the vertically translated phantom. (b) The surface map (left), intensity profiles at frequency of 0.1 mm−1 (center) and spatial frequency at each pixel (right) of the tilted phantom.
The flat phantom was then tilted up 30 deg relative to the reference plane along the x-axis. The sinusoidal light with the spatial frequencies of 0.05 mm−1, 0.1 mm−1, 0.15 mm−1 and 0.2 mm−1 were projected on the tilted phantom using the standard and profile-adaptive projection modes, respectively. The measured spatial frequency on the phantom at each pixel was calculated and compared between the two projection modes. When using the standard sinusoidal projection, the frequency of the sinusoidal light broadened over the inclined surface, with the maximum frequency increased by more than 20%. The profile-adaptive projection well corrected the frequency broadening of the sinusoidal light over the surface. The spatial frequency and phase of the sinusoidal light on the tilted phantom after correction remained the same as that on the reference plane, with the mean error of < 3% at each pixel, as shown in Fig. 7(b).
The optical properties of the flat phantom at each vertical height were then extracted from the intensity corrected and uncorrected data, respectively. The results are shown in Fig. 8. Without the profile corrections, the absorption coefficients were underestimated as the height increased, while the reduced scattering coefficients were overestimated. The mean errors of ${\mu _a}$ and ${\mu ^{\prime}_s}$ increased by more than 22% and 5% per centimeter, respectively. The frequency and intensity corrections mitigated the profile effects on the optical property extractions effectively, improving the results of both the absorption and reduced scattering coefficients at all heights. After correction, the mean errors of ${\mu _a}$ and ${\mu ^{\prime}_s}$ were reduced to <2% and <1% at each height, respectively. We also compared the results when applying the frequency and intensity corrections separately. The frequency correction alone reduced the mean error of ${\mu _a}$ by <10%, while the intensity correction alone can reduce the mean error by 60%. Relatively, the intensity correction was more effective than the frequency correction to improve the results of the absorption coefficient because the intensity correction compensated the DC and AC components simultaneously, while the frequency correction worked only on the AC component. However, since the reduced scattering coefficient was more sensitive to the AC diffuse reflectance which increased with the frequency correction while it decreased with the intensity correction as the surface height increased, the results of the reduced scattering coefficient failed to be improved by either correction alone.
Fig. 8. The optical properties of the flat phantom at each height extracted from the corrected and uncorrected data, at the wavelength of 635 nm. Mean of (a) the absorption coefficient and (b) the reduced scattering coefficient.
3.2 Profile corrections of the cylindrical and hemispherical phantoms
The cylindrical and hemispherical phantoms were then the samples to verify the performance of the profile corrections for uneven surfaces. The variations in the height and angle of the surface were 15 mm and 75 degrees at the maximum, respectively.
At first, the cylindrical phantom was used to validate the frequency correction for reducing the AC modulation amplitude distortions in single-snapshot demodulation. Both the standard sinusoidal and the profile-adaptive projections were used for measurements. The modulation amplitudes were extracted using the single-snapshot demodulation based on spectrum filtering. The AC reflectances of each projection mode were calculated and compared with that measured using the three-phase demodulation, which was not affected by frequency deviations. The reflectance measured by the three-phase demodulation in each projection mode was treated as the “gold standard” of the demodulation result.
The results are shown in Fig. 9. Before the frequency correction, the sinusoidal light was deformed over the surface of the cylindrical phantom, with the spatial frequency increasing in the center region with large surface heights and decreasing in the region with large surface angles on both sides. The AC reflectances extracted from single-snapshot demodulation were therefore distorted when compared with that measured by the three-phase demodulation. The profile-adaptive projection corrected the deformations of the sinusoidal light over the whole surface and consequently reduced the AC reflectance distortion. The mean error of the extracted AC reflectance was reduced from 11% to 2% after the frequency correction.
Fig. 9. The performance of the profile-adaptive projection for reducing distortions of the AC modulation amplitudes in single snapshot demodulation. (Left) The measured SFD images using the two projection modes, and (right) the calculated AC reflectances within the red dotted boxes in the SFD images.
The optical properties of the homogeneous cylindrical phantom were then extracted from the profile corrected and uncorrected data, respectively. The results are shown in Fig. 10. Before the profile corrections, the extracted optical properties were affected significantly by the surface profile. The mean errors of ${\mu _a}$ and ${\mu ^{\prime}_s}$ were more than 130% and 25%, respectively. The frequency and intensity corrections mitigated the profile effects of the surface with large heights and angles by compensating the deviations of the measured frequencies and reflectances, and consequently improved the results of both ${\mu _a}$ and ${\mu ^{\prime}_s}$. Overall, the mean errors of ${\mu _a}$ and ${\mu ^{\prime}_s}$ were reduced by more than 80% and by about 20%, respectively.
Fig. 10. The optical properties of the cylindrical phantom extracted from the uncorrected and corrected data, at the wavelength of 450 nm. (a) The 3D maps and (b) the profile curves at y = 20 mm of ${\mu _a}$ and ${\mu ^{\prime}_s}$. The black arrows in the color bars indicate the true values (${\mu _a}$ = 0.006 mm−1, ${\mu ^{\prime}_s}$ = 1 mm−1).
The homogeneous hemispherical phantom with a more complex surface was then used as a sample. The results of the absorption and reduced scattering coefficients extracted from the profile corrected and uncorrected data are shown in Fig. 11. Before the profile corrections, the ${\mu _a}$ values of more than half of the whole pixels were extracted with the mean error of > 40%. The errors occurred mainly on the side of the hemisphere away from the projector where the incident intensity attenuated significantly due to the increase in the incident angle. Although the reflection angle also increased over the surface, it varied less than the incident angle. Therefore, the results of absorption coefficient were more affected by the incident angles. This further proved the significance of the incident intensity correction. After the profile corrections, the ${\mu _a}$ values of more than 90% of the whole pixels were extracted close to the true value, with the mean error of < 9%. Similar results were obtained in extraction of ${\mu ^{\prime}_s}$. Overall, the mean errors of ${\mu _a}$ and ${\mu ^{\prime}_s}$ were reduced from > 40% to < 10% and > 14% to < 5%, respectively.
Fig. 11. The optical properties of the hemispherical phantom extracted from the uncorrected and corrected data, at the wavelength of 635 nm. (a) The 3D maps and (b) the profile curves at y = 20 mm and x = 0 mm of ${\mu _a}$ and ${\mu ^{\prime}_s}$. The black arrows in the color bars indicate the true values (${\mu _a}$ = 0.05 mm−1, ${\mu ^{\prime}_s}$ = 1 mm−1).
3.4 In vivo topography of opisthenar vessel
Figure 12 shows the results of in vivo topography of an opisthenar vessel. The measured 3D surface map is shown in Fig. 12(b). Heights of the opisthenar decreased slightly along the x-direction with the height differences of < 5 mm, and the region of the blood vessel was raised in height about 2 mm from the opisthenar tissues. The surface height and angle variations caused the errors in the optical property extractions. Before the profile corrections, the absorption coefficients of the blood vessel region and opisthenar tissues with large surface angles were overestimated, and the diameter of the blood vessel was reconstructed thicker than its true value. After the corrections, the absorption and reduced scattering coefficients of the opisthenar tissues and blood vessel were more physiologically reasonable, as shown in Fig. 12(c).
Fig. 12. In vivo topography of an opisthenar vessel. (a) White light image. The black box indicates the imaging field of view. (b) The measured 3D surface map. (c) The extracted absorption and reduced scattering coefficients from the uncorrected and corrected data, measured at the wavelength of 520 nm.
4. Discussion and conclusion
Several features of the proposed method are further discussed. Although the SFD imaging typically projects sinusoidal fringes along the acquisition plane which is less sensitive to the surface variations, the fringe deformation is still unavoidable on samples with large surface variations. Phantom experiments have validated the performance of the frequency correction to improve the AC amplitudes in single-snapshot demodulation based on spectrum filtering, though, the DC amplitudes may also be improved more or less since the correction eliminates the frequency broadening which otherwise may lead to the overlapped spectrum of the AC and DC components. And when multiple frequencies are demodulated from a single image, the frequency broadening will further increase the difficulty of filter design in both spectrum and spatial filtering methods. The profile-adaptive projection method can correct sinusoidal light distortions at different modulation frequencies and orientations on sample surface, so it is expected to effectively reduce the single-snapshot demodulation errors mentioned above.
In addition to improving the demodulation accuracy, the frequency correction also reduces the error of optical property extraction caused by spatial frequency inconsistency. Based on a Monte Carlo simulation, the effect of the spatial frequency deviation on extraction of the optical properties is further investigated [26]. Figure 13 shows the diffuse reflectance versus the spatial frequency of SFD imaging for samples with the varied optical properties. Diffuse reflectance decreases nonlinearly with the spatial frequency and shows a varied trend of decline with the varied optical properties. Relatively, the spatial frequency deviations have greater impact on the optical property extraction within the frequency range of 0.1-0.2 mm−1, since the diffuse reflectance decreases most rapidly in this range. And the diffuse reflectance is more sensitive to the changes in the absorption property at low frequencies, and more sensitive to the changes in the scattering property at high frequencies. Therefore, the frequency correction is more useful for improving accuracy of the extracted reduced scattering coefficient. In the phantom experiments, the frequency correction improved the extracted absorption and reduced coefficients by an average of 3% and 8%, respectively. Compared to the multi-frequency processing algorithm, the proposed profile-adaptive projection method is model-included to directly correct the frequency deviation on the sample, thus avoiding errors from interpolations when using multiple look-up tables. Moreover, the proposed frequency correction method is also applied to the SFD diffuse optical tomography, in which the multi-frequency processing algorithm is not applicable.
Fig. 13. The diffuse reflectance versus the spatial frequency of SFD imaging for samples with the varied optical properties, calculated using a Monte Carlo simulation.
The value of Minnaert’s constant, k, in the intensity correction reflects the profile effect on the light intensity. As discussed in [17], larger k values are required to correct the DC reflectances, which means that the light intensity at low frequency is more affected by morphology. Besides, we have also observed that the k values in the incident angle correction are larger than that in the reflection angle correction, which further proves the necessity of the incident intensity correction. Previous works only corrected the reflection angle and still obtained good results. It is likely because that the projector in their system is spatially close to the camera with small divergence angle (1.8 degree) and small ${\beta _z}$ angles (0 degree), which synchronize the incident and reflected angles at each position on the sample [16,17]. Thus, the incident intensity deviation is compensated to some extent in the correction of the reflection angle. Similar to the profile-adaptive projection in the frequency correction, the incident intensity deviation is also expected to be compensated using a pre-distorted projection pattern, by adjusting the modulation amplitudes based on the sample surface profile. But it may face the insufficient quantization gray levels of the spatial light modulator and the difficulty in determining the Minnaert’s constant.
There are limitations of our method in some respects. For example, the intensity correction includes both the incident and reflection angles, which makes the reflectance of the sample with large morphological changes more likely be overcorrected. The method has presented good correction performance for morphology changes of the sample within 3 cm in height and 70° in angle. However, in the case of surface angle >70°, slight overcorrections occurred. Additionally, since the value of ${R_d}_{,cor}$ in Eq. (11) refers to the constant term of the linear model and is assumed to be consistent throughout the image, the Minnaert’s constant calculated by a linear regression in the intensity correction may be invalid when the sample has large heterogeneity leading to greatly differed diffuse reflectances. The performance of the frequency and intensity corrections are also affected by the accuracy of profilometry. We simply modified the height-phase function to overcome the height difference of the system, however, the accuracy relies on the measured system parameters. A direct mapping algorithm can avoid system parameter measurements, but the error of a multi-height calibration scheme is also inevitable [27]. Moreover, the error also arises from the insufficient sampling of the high-frequency information that represents sharp changes in sample morphology and heterogeneity by single-pixel imaging. This could be improved by increasing the spatial frequency of the selected DCT coefficients.
Although being validated using the system based on single-pixel camera, our method is also applicable to the array camera-based SFD imaging system. Comparatively, single-pixel camera has limitations in the resolution and imaging speed, but it brings the superiority in cost, sensitivity and multi-wavelength measurements [28,29]. The current work takes about 10 seconds to acquire each SFD image or phase image. 5 images are required for each experiment (4-step phase images for 3D reconstruction and a single-snapshot image for optical property measurement), hence a total of 50 seconds is needed. However, for in vivo imaging of moving organs and real time monitoring of rapidly changing physiological parameters, the imaging speed has to be further improved. 3D profilometry technique of dynamic objects could be an alternative to phase profilometry [30,31]. Recently, artificial neural networks have provided new approaches for both single-pixel imaging and SFD imaging [32,33], improving both the imaging speed and quality. And the optical properties and 3D profile are expected to be simultaneously extracted from a single image [34,35]. These approaches are expected to be applied in our method. Since the frequency correction in this study requires a change in the projection pattern to adapt for the surface profile of the sample, 2 sequential projections of 3D profiling and SFD imaging are required, which can be a challenge to achieve real-time capability of the method. However, in some static imaging applications where the tissue is generally morphologically unchanged, the surface profile of the tissue can be measured in advance for subsequent corrections.
In conclusion, we have for the first time proposed a profile-adaptive projection scheme to correct the spatial frequency deviations for single-snapshot SFD imaging and have described the simultaneous corrections of the intensities incident onto and reflected from the samples. The method has improved the accuracy of single-snapshot SFD imaging in extracting the optical properties of the sample with surface curvatures. Further work will focus on the application of the frequency and intensity corrections to more complex samples and to SFD diffuse optical tomography for reconstructing the depth information of tissue heterogeneities.
Funding
National Natural Science Foundation of China (61575140, 62075156, 81871393).
Acknowledgments
The authors acknowledge the insightful suggestions and revisions given by Professor Yukio Yamada from University of Electro-Communications, Japan.
Disclosures
The authors declare that there are no conflicts of interest related to this article.
References
1. D. J. Cuccia, F. Bevilacqua, A. J. Durkin, and B. J. Tromberg, “Modulated imaging: quantitative analysis and tomography of turbid media in the spatial-frequency domain,” Opt. Lett. 30(11), 1354–1356 (2005). [CrossRef]
2. D. J. Cuccia, F. Bevilacqua, A. J. Durkin, F. Ayers, and B. J. Tromberg, “Quantitation and mapping of tissue optical properties using modulated imaging,” J. Biomed. Opt. 14(2), 024012 (2009). [CrossRef]
3. J. Vervandier and S. Gioux, “Single snapshot imaging of optical properties,” Biomed. Opt. Express 4(12), 2938–2944 (2013). [CrossRef]
4. K. P. Nadeau, A. J. Durkin, and B. J. Tromberg, “Advanced demodulation technique for the extraction of tissue optical properties and structural orientation contrast in the spatial frequency domain,” J. Biomed. Opt. 19(5), 056013 (2014). [CrossRef]
5. M. Xu, Z. L. Cao, W. H. Lin, X. L. Chen, L. F. Zheng, and B. X. Zeng, “Single snapshot multiple frequency modulated imaging of subsurface optical properties of turbid media with structured light,” AIP Adv. 6(12), 125208 (2016). [CrossRef]
6. S. D. Konecky, C. M. Owen, T. Rice, P. A. Valdés, K. Kolste, B. C. Wilson, F. Leblond, D. W. Roberts, K. D. Paulsen, and B. J. Tromberg, “Spatial frequency domain tomography of protoporphyrin IX fluorescence in preclinical glioma models,” J. Biomed. Opt. 17(5), 056008 (2012). [CrossRef]
7. J. Q. Nguyen, C. Crouzet, T. Mai, K. Riola, D. Uchitel, L. H. Liaw, N. Bernal, A. Ponticorvo, B. Choi, and A. J. Durkin, “Spatial frequency domain imaging of burn wounds in a preclinical model of graded burn severity,” J. Biomed. Opt. 18(6), 066010 (2013). [CrossRef]
8. S. Bélanger, M. Abran, X. Intes, C. Casanova, and F. Lesage, “Real-time diffuse optical tomography based on structured illumination,” J. Biomed. Opt. 15(1), 016006 (2010). [CrossRef]
9. M. Ghijsen, B. Choi, A. J. Durkin, S. Gioux, and B. J. Tromberg, “Real-time simultaneous single snapshot of optical properties and blood flow using coherent spatial frequency domain imaging (cSFDI),” Biomed. Opt. Express 7(3), 870–882 (2016). [CrossRef]
10. G. T. Kennedy, R. Stone II, A. C. Kowalczewski, R. Rowland, J. H. Chen, M. L. Baldado, A. Ponticorvo, N. Bernal, R. J. Christy, and A. J. Durkin, “Spatial frequency domain imaging: a quantitative, noninvasive tool for in vivo monitoring of burn wound and skin graft healing,” J. Biomed. Opt. 24(07), 1 (2019). [CrossRef]
11. R. Hachadorian, P. Bruza, M. Jermyn, A. Mazhar, D. Cuccia, L. Jarvis, D. Gladstone, and B. Pogue, “Correcting Cherenkov light attenuation in tissue using spatial frequency domain imaging for quantitative surface dosimetry during whole breast radiation therapy,” J. Biomed. Opt. 24(7), 071609 (2018). [CrossRef]
12. S. Gioux, A. Mazhar, and D. J. Cuccia, “Spatial Frequency Domain Imaging in 2019: Principles, Applications and Perspectives,” J. Biomed. Opt. 24(7), 071613 (2019). [CrossRef]
13. D. J. Wirth, M. Sibai, B. C. Wilson, D. W. Roberts, and K. Paulsen, “First experience with spatial frequency domain imaging and red-light excitation of protoporphyrin IX fluorescence during tumor resection,” Biomed. Opt. Express 11(8), 4306–4315 (2020). [CrossRef]
14. T. T. Nguyen, H. N. Le, M. Vo, Z. Y. Wang, L. Luu, and J. C. Ramella-Roman, “Three-dimensional phantoms for curvature correction in spatial frequency domain imaging,” Biomed. Opt. Express 3(6), 1200–1214 (2012). [CrossRef]
15. P. Diep, S. Pannem, J. Sweer, J. Lo, M. Snyder, G. Stueber, Y. Y. Zhao, S. Tabassum, R. Istfan, J. J. Wu, S. Erramilli, and D. Roblyer, “Three-dimensional printed optical phantoms with customized absorption and scattering properties,” Biomed. Opt. Express 6(11), 4212–4220 (2015). [CrossRef]
16. S. Gioux, A. Mazhar, D. J. Cuccia, A. J. Durkin, B. J. Tromberg, and J. V. Frangioni, “Three-dimensional surface profile intensity correction for spatially modulated imaging,” J. Biomed. Opt. 14(3), 034045 (2009). [CrossRef]
17. Y. Y. Zhao, S. Tabassum, S. Piracha, M. S. Nandhu, M. Viapiano, and D. Roblyer, “Angle correction for small animal tumor imaging with spatial frequency domain imaging (SFDI),” Biomed. Opt. Express 7(6), 2373–2384 (2016). [CrossRef]
18. T. A. Erickson, A. Mazhar, D. J. Cuccia, A. J. Durkin, and J. W. Tunnell, “Lookup-table method for imaging optical properties with structured illumination beyond the diffusion theory regime,” J. Biomed. Opt. 15(3), 036013 (2010). [CrossRef]
19. V. Srinivasan, H. C. Liu, and M. Halioua, “Automated phase-measuring profilometry: a phase mapping approach,” Appl. Opt. 24(2), 185–188 (1985). [CrossRef]
20. J. Geng, “Structured-light 3D surface imaging: a tutorial,” Adv. Opt. Photonics 3(2), 128–160 (2011). [CrossRef]
21. M. Minnaert, “The reciprocity principle in lunar photometry,” Astrophys. J. 93, 403–410 (1941). [CrossRef]
22. T. X. Li, Z. P. Qin, X. Hou, M. Dan, J. Li, L. Zhang, Z. X. Zhou, and F. Gao, “Multi-wavelength spatial frequency domain diffuse optical tomography using single-pixel imaging based on lock-in photon counting,” Opt. Express 27(16), 23138 (2019). [CrossRef]
23. B. L. Liu, Z. H. Yang, X. Liu, and L. A. Wu, “Coloured computational imaging with single-pixel detectors based on a 2D discrete cosine transform,” J. Mod. Opt. 64(3), 259–264 (2017). [CrossRef]
24. S. Dokouzyannis and T. Tziortzios, “High throughput and energy efficient two-dimensional inverse discrete cosine transform architecture,” IET Image Process. 7(5), 533–541 (2013). [CrossRef]
25. B. W. Pogue and M. S. Patterson, “Review of tissue simulating phantoms for optical spectroscopy, imaging and dosimetry,” J. Biomed. Opt. 11(4), 041102 (2006). [CrossRef]
26. A. R. Gardner and V. Venugopalan, “Accurate and efficient Monte Carlo solutions to the radiative transport equation in the spatial frequency domain,” Opt. Lett. 36(12), 2269–2271 (2011). [CrossRef]
27. W. S. Zhou and X. Y. Su, “A direct mapping algorithm for phase-measuring profilometry,” J. Mod. Opt. 41(1), 89–94 (1994). [CrossRef]
28. M. Torabzadeh, Y. Park, R. A. Bartels, A. J. Durkin, and B. J. Tromberg, “Compressed single pixel imaging in the spatial frequency domain,” J. Biomed. Opt. 22(3), 030501 (2017). [CrossRef]
29. E. Aguénounon, F. Dadouche, W. Uhring, N. Ducros, and S. Gioux, “Single snapshot imaging of optical properties using a single-pixel camera: a simulation study,” J. Biomed. Opt. 24(7), 071612 (2019). [CrossRef]
30. H. T. Yun, B. W. Li, and S. Zhang, “Pixel-by-pixel absolute three-dimensional shape measurement with modified Fourier transform profilometry,” Appl. Opt. 56(5), 1472–1480 (2017). [CrossRef]
31. J. S. Hyun and S. Zhang, “Superfast 3d absolute shape measurement using five binary patterns,” Opt. Lasers Eng. 90, 217–224 (2017). [CrossRef]
32. C. F. Higham, R. Murray-Smith, M. J. Padgett, and M. P. Edgar, “Deep learning for real-time single-pixel video,” Sci. Rep. 8(1), 2369 (2018). [CrossRef]
33. Y. Y. Zhao, Y. Deng, F. Bao, H. Peterson, R. Istfan, and D. Roblyer, “Deep learning model for ultrafast multifrequency optical property extractions for spatial frequency domain imaging,” Opt. Lett. 43(22), 5669–5672 (2018). [CrossRef]
34. M. van de Giessen, J. P. Angelo, and S. Gioux, “Real-time, profile corrected single snapshot imaging of optical properties,” Biomed. Opt. Express 6(10), 4051–4062 (2015). [CrossRef]
35. E. Aguenounon, J. T. Smith, M. Al-Taher, M. Diana, X. Intes, and S. Gioux, “Real-time, wide-field and high-quality single snapshot imaging of optical properties with profile correction using deep learning,” Biomed. Opt. Express 11(10), 5701–5716 (2020). [CrossRef]
