Single snapshot of optical properties image quality improvement using anisotropic two-dimensional windows filtering

Abstract. Imaging methods permitting real-time, wide-field, and quantitative optical mapping of biological tissue properties offer an unprecedented range of applications for clinical use. Following the development of spatial frequency domain imaging, we introduce a real-time demodulation method called single snapshot of optical properties (SSOPs). However, since this method uses only a single image to generate absorption and reduced scattering maps, it was limited by a degraded image quality resulting in artifacts that diminished its potential for clinical use. We present filtering strategies for improving the image quality of optical properties maps obtained using SSOPs. We investigate the effect of anisotropic two-dimensional filtering strategies for spatial frequencies ranging from 0.1 to 0.4  mm−1 directly onto N=10 hands. Both accuracy and image quality of the optical properties are quantified in comparison with standard, multiple image acquisitions in the spatial frequency domain. Overall, using optimized filters, mean errors in predicting optical properties using SSOP remain under 8.8% in absorption and 7.5% in reduced scattering, while significantly improving image quality. Overall this work contributes to advance real-time, wide-field, and quantitative diffuse optical imaging toward clinical evaluation.


Introduction
Optical imaging methods capable of providing real-time information related to functional and structural conditions of living tissues are becoming increasingly popular. [1][2][3][4][5][6] Some medical fields in particular, such as surgery, are driving this need for realtime interpretable information content to aid decision-making in a time-constrained environment. [7][8][9][10][11][12] Several solutions have been investigated to fulfill this need, each with their pros and cons, ranging from raster scanning of microscopic information to wide-field acquisitions of diffused information. 13,14 For instance, fluorescence imaging to help surgeons visualizing structures of interest, such as lymph nodes, ureters, or micrometastasis, stands as an example where providing tissue-related information in real time plays a key role in clinical adoption. [15][16][17][18] Among all existing solutions, spatial frequency domain imaging (SFDI) has been recently developed with the unique capacity to provide rapidly quantitative diffuse information (absorption and reduced scattering) over large fields of view (typically, >10 × 10 cm 2 ). 14,19 Even more recently, the SFDI method has been improved to image optical properties more rapidly using two images 20 or even a single image. 21,22 The latter method, termed single snapshot of optical properties (SSOP), has the unique capacity to image fully in real time with the only limitation being the reflectance signal-to-noise ratio at high camera frame rates. However, improving speed of acquisition comes at the cost of image quality, a factor that cannot be neglected when aiming at translating an imaging method to the clinic. It is, therefore, necessary to provide at the same time speed and good quality images so they can be fully utilized by clinicians.
In this paper, we present means to improve image quality when using SSOPs. We investigate in particular two parameters: the filtration strategy in the spatial frequency domain, and the choice of spatial frequency. We present the effect of these parameters on both recovered optical properties and image quality, in comparison with standard, non-real-time SFDI. Finally, we propose an optimal set of parameters allowing to image optical properties accurately while maintaining good image quality.

Spatial Frequency Domain Imaging
The theory of SFDI is described in the literature. 14,19 Briefly, it consists of analyzing the spatial modulation transfer function of a turbid medium for every pixel of an image at once. SFDI makes measurements relative to a tissue-mimicking calibration phantom of known optical properties and uses a model-based approach, typically diffusion theory or Monte Carlo simulations, to extract the optical properties. In the case of subsurface imaging, where optical properties are considered to be independently measured at every location in the image, a fast, precomputed lookup table (LUT) can be used to recover the optical properties from only two spatial frequency images. Typically, one DC image (e.g., 0 mm −1 ), sensitive to changes in both reduced scattering and absorption, and one AC image (e.g., 0.2 mm −1 ), mainly sensitive to changes in reduced scattering, are used. 14,23,24 The most common approach is to use a projector and a camera to acquire for each spatial frequency (e.g., f 0 ¼ 0 and f x ¼ 0.2 mm −1 ) three images phase shifted by 120 deg (φ1, φ2, and φ3), hence a total of six images (Fig. 1). When using the three-phase demodulation technique, an analytical expression can then be used to extract the AC and DC components for each pixel x i of the sample (M AC and M DC ): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 6 3 ; 5 0 0 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 6 3 ; with E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 7 4 1 Extracting the DC modulation could be performed in a different way by calculating the mean of the images obtained at three phases at a high spatial frequency [Eq. (21) in Ref. 14]. However, when performing this operation, dark noise and other constant noise [I offset in Eq. (3)] are not eliminated and lead to incorrect estimation of diffuse reflectance. In this particular case, as well as when performing SSOP, the contribution of dark noise and other constant noise is eliminated by subtracting a dark image from the acquired image. When using SFDI, the acquisition of multiple images along with Eq. (2) allows to eliminate constant noise from the image through image subtraction. We chose the latter method for our reference SFDI data processing.
Following demodulation, the medium diffuse reflectance R d;AC and R d;DC are obtained using a calibration reference with known optical properties: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 5 3 8 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 4 9 0 Using the medium diffuse reflectance, a precomputed LUT can be used then to recover reduced scattering (μ 0 s ) and absorption (μ a ). This process is summarized in Fig. 2(a). Fig. 1 Schematics of the SFDI and SSOP imaging system. A laser diode source is coupled to a digital micromirror device (DMD) using a 1-mm-diameter optical fiber. Sinusoidally modulated patterns are projected onto the field of view and collected using a sCMOS camera.

Single Snapshot of Optical Properties
The principle of SSOP 21,22 consists of acquiring a single image at a high spatial frequency (e.g., f x ¼ 0.2 mm −1 ) to reduce the acquisition time and extract the M AC and M DC images by filtering the acquired image directly in the spatial frequency domain [Fig. 2(b)]. A one-dimensional (1-D) or a two-dimensional (2-D) Fourier transform (FT) is performed on the image, the modulation frequency is detected in Fourier domain in the same manner as in FT profilometry, 25,26 and an ideal rectangular filter is used to separate low spatial frequency components from the high spatial frequency components. The M DC image is then obtained using simple inverse FT and the M AC image using Hilbert transform. Following this demodulation step, data are processed in an identical fashion as with SFDI: a calibration reference is used, the medium diffuse reflectance is obtained, and finally, reduced scattering and absorption maps are extracted. There are, however, major issues of this method that impact accuracy and image quality: first, the Gibbs effect due to the ideal rectangular filtering and second, the difficulty to properly extract the AC component with just one phase. In the following sections, we study the different filtering windows to obtain the M AC and M DC images and we investigate the means to reduce Gibbs effect, reduce fault value due to edge discontinuity, keep optical properties close to the reference multiphase SFDI values, and improve in the same time the visual aspect.

Design of the Single Snapshot of Optical Properties Demodulation Filters
From the literature, the classic rectangular filtering is not causal, and, even if it gives a most perfect response, is subject to ripple. 27 To reduce the ripple effect and improve image quality, other well-known windows were investigated, such as Sine, Hann, Blackman, and Gaussian. 27 These windows have known properties in 1-D (e.g., maximum side lobe level, side lobe roll-off rate, −3 dB main lobe width, and scalloping loss), but because images are two-dimensional, 2-D versions of these windows have to be generated, inherently altering their characteristics. There exists different ways to create 2-D windows from 1-D windows, such as simple product or circular rotation, 28 and for simplicity of implementation in Matlab, a simple product of the x and y directions of the filter characteristic was employed. The design rationale for 2-D anisotropic demodulation filters is as follows. The filter should capture as much information as possible in the frequency domain, information that should be related to the modulation frequency of interest (here, DC or AC). The cutoff frequencies should be chosen so information content in the spatial frequency domain is well separated between modulation frequencies (DC and AC). Because there is no known method to theoretically determine the best filter shape or cutoff frequencies, we developed our own method.
In the x direction, defined as the direction that contains spatially modulated information, filters were designed around the zero frequency for the DC filter and around the projected spatial frequency for the AC filter. Filters sizes in this direction were incremented by steps of 0.0067 mm −1 to study the effect of the cutoff frequency onto demodulated images. In the y direction, because there is no spatially modulated information, the filter shape was maximized to contain as much information as possible. Filters sizes in this direction were then kept constant and equal to the maximum spatial frequency (here, 6.83 mm −1 ). Two-dimensional filters were then obtained by multiplying the filters in the x and y directions, resulting in elliptical-shape anisotropic filters.
Finally, for each 2-D filter shape (with varying cutoff frequency in the x direction), a single-phase SSOP image of a complex in vivo sample (a human hand) was demodulated. The demodulated images were compared with those obtained using standard SFDI demodulation as a reference. The best filter parameters were chosen when demodulated images with SSOP best matched the ones obtained with SFDI. These operations were repeated for different projected spatial frequencies. Using this method, standard filter parameters were defined: center and cutoff frequencies for both DC and AC.

Effect of Spatial Frequency
As our interest is to separate DC and AC content in the frequency domain, the value of the projected spatial frequency is of paramount importance. The higher the spatial frequency, the lower the cross talk between the DC and AC components and the easier it is to separate contributions from the M DC and M AC . Within the limits of our projection system, a range of four spatial frequencies were tested: f x ¼ 0.1, 0.2, 0.3, and 0.4 mm −1 . An additional acquisition at f x ¼ 0 mm −1 was performed for processing with SFDI. For each spatial frequency, SSOP processing was performed using the AC and DC filters, and the resulting images were compared both in optical properties value and in image quality.

Instrumental Setup, Phantoms, and Samples
The instrumental setup was custom-built using a DMD (Vialux, Germany) for the projection of custom patterns, fiber-coupled to a 665-nm laser diode (LDX Optronics, Maryville, Tennessee). The projection system projects a sine wave pattern over a 150 × 150 mm 2 field of view at 47-cm working distance. About 1024 × 1024 images are acquired using a scientific monochrome 16 bits CMOS camera (PCO AG, pco.edge 5.5, Kelheim, Germany). Polarizers (PPL05C; Moxtek, Orem, Utah), arranged in a crossed configuration, are used to minimize the contribution from specular reflections at the surface of the sample.
Silicone-based optical phantoms were built using titanium dioxide (TiO 2 ) as a scattering agent and India ink as an absorbing agent. 29,30 One large calibration phantom was made (210 mm × 210 mm × 20 mm in size) with reduced scattering μ 0 s ¼ 1.08 mm −1 and absorption μ a ¼ 0.012 mm −1 at 665 nm. Anatomically relevant in vivo samples were used to best evaluate the potential of the proposed method for improved accuracy and image quality in clinically relevant conditions.
Note: d is the field of view width in mm.
Journal of Biomedical Optics 071611-4 July 2019 • Vol. 24 (7) Because hands are having complex variations in height, shape, and heterogeneous composition, they are ideal samples to evaluate the processing method proposed in this work.

Data Processing and Results Analysis
Data were processed using both SSOP and SFDI methods as described previously. 21,22 Custom processing code in Matlab (MathWorks, Natick, Massachusetts) using an LUT extraction method was generated from white Monte Carlo and, therefore, valid for high spatial frequencies beyond the diffusion limit. 14 For the evaluation of the results accuracy, optical properties maps of N ¼ 10 hands, complex heterogeneous samples, obtained using the SSOP method for each filter set and each spatial frequency were compared with optical properties maps obtained from the SFDI method. Variations of optical properties between the two methods were calculated in terms of quality and fidelity with a custom percentage error formula over the entire image: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 6 3 ; 5 4 6 Error% ¼ 100 × ðSFDI n;m − SSOP n;m Þ SFDI n;m : (6) For the evaluation of the visual quality of the images, we measured the full width at half-maximum (FWHM) of a vein of N ¼ 10 hands onto the optical properties maps obtained from the SSOP method using the AC and DC filters in comparison with the SFDI method. Images were also analyzed visually to take into account degradations, such as ripples and edge artifacts. Mean and standard deviations were calculated for the N ¼ 10 samples.

Filters Design
We evaluated a total of 39 common filters (13 DC filters and 26 AC filters) for a total number of filtering combinations of 338.
All filters are detailed in Sec. 6. For ease of interpretation, we present only the results obtained with the filters having the best performances: the standard DC and AC rectangular filters previously used with SSOP (DC1 and AC1, respectively), a DC Blackman filter (DC2), a DC Sine filter (DC3), an AC Blackman filter (AC2), and a half Blackman slope filter (AC3). All combinations of these filters are presented in the following sections leading to nine filter sets tested. A detailed summary of these filters illustrated in Fig. 3 is provided in Table 1.

Absorption images
Results quantifying accuracy for the selected filters on N ¼ 10 hands are shown in Fig. 4, color-coded for ease of interpretation.
Overall the mean error expressed in percentage value for absorption varies from 25.31% in the worst case to 8.84% in the best case. Low mean error percentage value indicates a better quality, and we added a colorbar along the numerical values for easy interpretation. Results for all 338 filters combinations are shown in Sec. 6. More precisely, the best combinations to extract absorption, as a function of spatial frequency, are as follows: Globally, these results indicate that DC3 and DC2 filters allow to extract absorption properties more accurately than the rectangular filter (DC1). The DC3 filter consistently shows best results for all spatial frequencies. When combined with AC filters, best results are obtained using the AC2 and AC3 filters.

Reduced scattering images
Results quantifying accuracy for the selected filters on N ¼ 10 hands are shown in Fig. 5, color-coded for ease of interpretation.
Overall the mean error expressed in percentage value for reduced scattering varies from 14.47% in the worst case to 7.46% in the best case. Results for all 338 filters combinations are shown in Sec. 6. More precisely, the best combinations to extract reduced scattering, as a function of spatial frequency, are as follows: Globally, these results indicate that most accurate reduced scattering properties are extracted with the AC3 filter at low spatial frequency and with the AC2 at high spatial frequency. When combined with DC filters, best results are obtained using either DC2 or DC3 filters.

Quantification of Image Visual Quality
Assessing the visual quality of an image is not straightforward. To overcome this issue, we chose to quantify the image quality by measuring the FWHM of the veins of N ¼ 10 hands on the extracted absorption images for all filters in comparison with SFDI. An example of a chosen vein is shown in Fig. 7 (white bar). FWHM % errors for all filter combinations are shown in Fig. 6, color-coded for ease of interpretation. Results for all 338 filters combinations are shown in Sec. 6.  Results can be summarized as follows: Overall, the visual quality of the image increases with spatial frequency. At 0.4 mm −1 , DC1 combinations and DC3 and AC3 combinations allow extraction of absorption properties with the best visual quality, compared to SFDI absorption images.

Qualitative Analysis
To appreciate the improvement in image quality, we provide here an example of the N ¼ 10 hands imaged in the study. Figures 7 and 8 show, respectively, absorption and reduced scattering images obtained with SFDI (first column), SSOP using the standard rectangular filters DC1 and AC1 (SSOP std), and two optimized SSOP filters: DC1 and AC3 (SSOP opt1) and DC3 and AC3 (SSOP opt2).
Qualitatively, one can first appreciate the image improvement with increasing spatial frequency and filter combination choice. Fine details such as the veins in absorption appear more clearly as spatial frequency is increased and ripples around the edges as well as on the background are strongly diminished as filtration quality is improved.
Finally, a movie of a hand was acquired making a quarter turn with tight fingers at the beginning and spread fingers at the end to demonstrate the filtering improvement onto data acquired in real time (Fig. 9). This movie is acquired at a spatial frequency of 0.4 mm −1 , a true frame rate of 50 images per second, and processed with the standard rectangular DC1 and AC1 filters combination [ Fig. 9(a)] and an improved DC3 and AC3 filters combination [ Fig. 9(b)].

Discussion
In this study, we have investigated different 2-D filtering strategies allowing improvement of the SSOP method's visual aspect while preserving good optical properties measurement accuracy. Additionally, this work provided an opportunity to precisely design equations for 2-D anisotropic filters in the spatial frequency domain. First, we were able to define DC and AC filters shape and cutoff frequencies by comparing demodulation accuracy with SFDI as a reference. Then, the filters that give best results in terms of optical properties and image quality were combined, applied onto complex in vivo objects (N ¼ 10 hands), and results presented in this paper. The extracted optical properties at different spatial frequencies from these filters were compared to SFDI, and their visual aspects were assessed with measurement at mid-height of the veins.
In summary, using an acquisition spatial frequency of 0.4 mm −1 , combinations using the DC1 and AC3 or DC3 and AC3 filters allow extraction of absorption and reduced scattering properties with good quantitative value while providing an improved visual aspect at the same time. In addition, not quantified in the visual quality, ripples are significantly reduced from both the background phantom and the in vivo sample.
However, this work is not without limitations, in particular regarding the quantitative comparison metrics, the precision of the comparison reference using SFDI, and the image quality metrics.
Regarding the quantitative value metric, the formula used [Eq. (6)] gives a unique value for all the pixels of the image. This implies that the pixels at the background and at the edges of the hand have the same importance as the pixels at the center of the hand. For example in Fig. 7 (absorption images), at 0.4 mm −1 vein, structures of the hand using SSOP std are better visualized than using SSOP opt2. But because there are more errors on SSOP std edges and background, the mean error percentage of SSOP opt2 is lower. Therefore, when choosing a filtering combination, a compromise is necessary between the defaults at the edges of the image, the defaults in the background, and the accuracy of the values in the region of interest, depending of the imaging application.
Regarding the use of SFDI as a reference for the quantitative comparison, noise related to multiple-phases acquisition can result in stripe artifacts, as seen in reduced scattering in Fig. 8. Because amplitude modulation and thereby signal-to-noise significantly decrease when increasing spatial frequency, 14 SFDI demodulation is becoming more noisy when SSOP demodulation still gives good results. This increase in noise partially explains why the mean error percentage values for reduced scattering are not decreasing when increasing the spatial frequency from 0.3 to 0.4 mm −1 as one would expect.
It should also be noted that the percent error was computed over the entire hand giving a result that appears to be high. We purposely chose this approach to evaluate errors on both flat zones (back of the hand, for instance) and zones with height changes (edge of the fingers, for instance). The errors resulting from the edge artifacts that we attempted to reduce in this work contribute largely to a high error value for the entire hand. Clinically, the accuracy necessary for proper use of this technology remains unknown and dependent on the clinical scenario.
We chose to perform the image quality assessment by measuring the width of a vein on the back of N ¼ 10 hands. This choice was motivated by the fact that it reflects the image quality metric expected by the end-user of the method. We did investigate other known image quality metrics, such as BRISQUE, NIQUE, or SSIM. However, none of these metrics provided satisfactory results for quantifying image quality. This is partially due to the fact that these algorithms have been developed for a specific context of image quality that does contain features in the entire image. In our case, in vivo images contain both a flat, featureless zone (the reference phantom) and a zone containing features (the hand). This results in an inaccurate assessment of the impact of image artifacts, such as ripples, on flat zones, leading to a disagreement with simple visual evaluation. We, therefore, preferred to use a simpler metric reflecting the resolution of absorption within the image, by measuring the width a vein.
Overall, this quantitative assessment proves to be effective for absorption images, showing an image quality improvement as spatial frequency increases, as one would expect and subjectively confirms. However, this measurement does not take into account the background of the image or the edge ripples. Finding a way to better quantify image quality would certainly improve such study. In addition, because veins at the back of hands are common absorption features, we were able to perform this measurement on several hands. Unfortunately, scattering features that would be similar between individuals are not common. We, therefore, did not perform the image quality assessment onto scattering images. However, this is not a strong limitation since it is well known that scattering images do not suffer quality degradation as much as absorption images. This is due to the fact that higher spatial frequencies are preserved in scattering images and, as a consequence, image quality is higher.
Of interest as well, this study offers a better understanding of both quantification and image quality as a function of spatial frequency. From a pure image analysis standpoint, it may appear obvious that increasing spatial frequency would improve demodulation in the case of SSOP by providing better separation between DC and AC content. However, image quality suffers when increasing spatial frequency with an increase in noise due to lower signal-to-noise ratio that, depending on the instrumentation used, can lead to inaccuracies in extracting optical properties. 14 This topic is currently being investigated in particular regarding instrumental contributions to both bias and image quality. In addition, it is well known that increasing the spatial frequency results in more superficial measurement in the spatial frequency domain. 14 This is well demonstrated in the reduced scattering maps (Fig. 8) with a noticeable change in properties with spatial frequency. On the contrary, absorption maps remain nearly identical for all spatial frequencies (Fig. 7). This result is interesting in the context of surgical guidance where functional parameters, derived from absorption, are particularly useful. This result justifies the use of a high spatial frequency acquisition in the interest of image quality while preserving accurate functional measurements. More studies are necessary to understand the implication of this result onto measuring reduced scattering in the context in vivo applications. Finally, it is interesting to note that some filter combinations perform better for absorption accuracy while others perform better for reduced scattering accuracy. One may, therefore, consider the option of using different filter combinations for extracting absorption and reduced scattering to obtain as accurate results as possible. This, however, comes at the cost of an increase in processing time.

Conclusion
Filtering strategies for SSOP were presented, using 2-D anisotropic windows for filtering in the Fourier domain and tested onto N ¼ 10 hands. Using custom mean error percentage computation and FWHM vein measure, this study demonstrates that it is possible using Blackman and Sine windows at high spatial frequencies (typically, 0.4 mm −1 ) to significantly improve image quality while preserving accuracy in optical properties measurements. This work contributes to the development of real-time, quantitative, and wide-field diffuse optical imaging methods with the ultimate goal of providing real-time objective feedback during surgery. Tables 2 and 3 show the characteristics of all filters tested during this study. Figs. 10-20 show all results from all filters combinations tested during this study.

Appendix
Exponential bandpass

Triangular bandpass
Welch bandpass Lanczos bandpass Hann bandpass Hamming bandpass Blackman bandpass Nuttall bandpass Journal of Biomedical Optics 071611-10 July 2019 • Vol. 24(7) Gaussian bandpass Note: d is the field of view width in mm.
Journal of Biomedical Optics 071611-11 July 2019 • Vol. 24(7) Table 3 Bandpass and highpass AC filter equations with their cutoff frequency and size.

Harris highpass
Blackman-Nuttall highpass Note: d is the field of view width in mm.
Journal of Biomedical Optics 071611-14 July 2019 • Vol. 24(7) Fig. 10 Mean percentage error over 10 measured hands for absorption at 0.1-mm −1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right).

Fig. 11
Mean percentage error over 10 measured hands for absorption at 0.2-mm −1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right).
Journal of Biomedical Optics 071611-15 July 2019 • Vol. 24(7) Fig. 12 Mean percentage error over 10 measured hands for absorption at 0.3-mm −1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right).

Fig. 13
Mean percentage error over 10 measured hands for absorption at 0.4-mm −1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right).
Journal of Biomedical Optics 071611-16 July 2019 • Vol. 24 (7) Fig. 14 Mean percentage error over 10 measured hands for reduced scattering at 0.1-mm −1 spatial frequency using all filters combinations. Color-coding is used for ease of interpretation (scale on the right).

Disclosures
The authors declare that there are no conflicts of interest related to this paper.