Spatiotemporal processing in photoplethysmography for skin microcirculatory perfusion imaging

: Technological advances in the real-time visualization of cutaneous microcirculation aim to realize benefits including high-resolution imaging, suppressed noise, and robust temporal coherence. Photoplethysmography (PPG), a noninvasive technique that measures single or multiple points of relative blood volume changes in blood vessels under the skin, shows potential as a signal candidate for visualizing blood vessels and tracking blood flow. However, challenges still remain, such as extracting/image reconstruction of the blood vessel/flow signal in a precise frequency window ( < 0.2 Hz) from a noisy image that is caused by the loss of spatial coherence of the light source in a turbid biological tissue. We attempted to overcome this challenge by adopting a combination of direct-contact-type, lens-less, conformable imagers and singular value decomposition (SVD) in this study. We focused on the numerical analysis of SVD for discriminating the tissue and vein blood flow in PPG for reconstructing blood fluidic images, followed by a complete demonstration of skin microcirculation blood tracking in the vessel visualization process when applying our lens-less, conformable, wearable imagers.


Introduction
The progression of the real-time visualization of cutaneous microcirculation involves developing intravital microscopy, nailfold video capillaroscopy, orthogonal polarization spectral imaging techniques, sidestream dark field imaging, and laser doppler imaging (LDI), aiming to achieve high-resolution imaging, suppressed noise, robust temporal coherence, etc., for the study and diagnosis of vascular conditions such as trauma, aneurysms, and vascular malformations [1].
Recently, skin microcirculatory has been studied extensively with the development of methods such as laser Doppler flowmetry (LDF), LDI, and optical coherence tomography (OCT). [2][3][4][5]. Such sensing and imaging technologies evolve with the development of precision of blood perfusion rate measurement, expansion of target areas, and new visualization methods (e.g., holography) for under-skin blood fluidic activity. The fundamental principle of such optical techniques is detecting coherent light that is transmitted through or reflected by tissue containing flowing blood [6]. Following calibration, the extent of coherent light scattering is used to calculate the blood flow velocity or perfusion rate [7]. To extract Doppler-shifted signal components in LDI or OCT, singular value decomposition (SVD) has been studied and applied to diagnose microcirculatory activities based on its ability to separate signal components in different frequency bands effectively in the high-frequency domain [8,9]. Additionally, both Fourier and SVD applications in laser Doppler methods have been reported, demonstrating strong performance for filtering noise and revealing blood flows at different fluidic velocities in retinal motion [10,11]. SVD is a powerful mathematical model for LDI, and for improving ultrasound vessel imaging using advanced SVD-based algorithms, which have been developed by several labs [12][13][14].
Interestingly, photoplethysmography (PPG) signals have been demonstrated to be feasible for monitoring skin microcirculation during the early stages of analyzing its similarity to LDF signals in low-frequency bands (<2 Hz) [15,16]. Mizeva et al. [17] reported the correlation between the PPG and LDF of microvascular low-frequency oscillations. Their study demonstrated a strong correlation (approximately 0.9) between PPG and LDF for specific microvascular physiological bands in myogenic activity (0.052-0.145 Hz), neurogenic activity (0.021-0.052 Hz), and endothelial activity (0.0095-0.021 Hz). Therefore, they suggested that PPG has the potential to be a low-cost replacement for low-frequency-band vascular activity assessment.
Multiple point measurements in tissue circulation monitoring systems have been developed to monitor vascular activities based on the PPG method in combination with signal analysis techniques, such as the standard Fourier-based time-frequency analysis method in our previous work [18]. The sensing mechanism of PPG measures reflected or transmitted light through tissue or blood vessels, thereby capturing signal variations in blood volume changes corresponding to heartbeats. Based on this mechanism, the PPG method is powerful for sensing signals in the low-frequency domain (0-2 Hz), which makes it possible to design less complicated and costly signal acquisition systems, and raises the possibility of developing wearable PPG sensors [19,20].
However, although PPG has robust temporal coherence when captured through direct-contacttype sensors, clutter filtering is a challenge because noise stemming from clutter tissue is extremely close to the blood vessel/flow signals (frequency window <0.2 Hz). One of the reasons is the spatial coherence loss of an NIR LED light source as it penetrates turbid biological tissue [21]. Considering the outlook for high resolution, wearable or smartphone connected blood flow monitoring devices, it is necessary to find an alternative for complicated math models (e.g., high order Fourier filters) that can secure suppressed noise blood flow imaging efficiently while avoiding the need for high system computational ability. Additionally, such noise cannot be resolved via conventional low-order Fourier filters because it is difficult to be extracted if it resides in the middle of each close and precise frequency window of blood flow signal (i.e., blood flow signals are corrupted with tissue clutter noise signals) [22]. It directly affects the spatial coherence of direct-contact-captured PPG (i.e., highly noisy vein imaging). The other major challenge is that, although taking advantage of the Doppler effect and clutter filtering using SVD is promising for estimating relative blood flow velocity, it cannot track blood flow positions to estimate absolute blood flow velocity.
In this study, we first performed a) blood flow tracking/absolute velocity estimation and b) vessel visualization simultaneously by tracing the positions of relative volume changes in blood vessels (i.e., PPG) during blood transport. To achieve this, we combined two imaging technologies: a) direct-contact-type, lens-less, conformable imagers [23] to collect high-quality PPG signals; and b) SVD with numerical analysis to determine the singular value thresholds for clutter, noise filtering, and reconstructing vessel and blood flow signals from PPG signals.

Biometric signal capture
We used organic photodiodes coupled to an active matrix backplane with the low-temperature polycrystalline silicon thin-film transistor readout circuit setup, presented in [23]. Veins were observed using transmitted light. An optical lens (Canon Macro Lenz EF-S 60 mm) was placed on the imager, and its focus was adjusted. The distance between the finger and lens was 200 mm. The veins were imaged using all sensors (252 × 256). An LED with a wavelength of 850 nm (OSRAM, SFH4052) was placed on top of the finger. The light intensity of the LED was 51 µW/cm 2 at a distance of 330 mm and 520 µW/cm 2 at a distance of 100 mm. For the hand palm vein measurement, conformable imagers (128 × 168) directly contacted the palm, and an LED light source (NISSIN ELECTRONIC CO. LTD, CR-30S-L) with a wavelength of 850 nm (40 mW/cm 2 ) was placed on the back of the hand. Figure 1(a) presents consecutive images generated by transmitted light traveling through the finger to the image sensor array. The collected signal information includes static biometric signals (blue pixels), dynamic biometric signal fluctuations caused a by more vascular activity inside the finger veins (yellow pixels), and null signals from the non-illuminated pixel area (red pixels). In our definitions for dynamic/static biometric signals shown in Fig. 1(b), static biometric signals (i.e., low-frequency fluctuating signals) in the non-vein generate stronger digital signals compared to a dynamic biometric signal. However, dynamic biometric signals in vein areas yield lower DC component digital signals, but result in more intense fluctuation because of a more tense vascular activity in the veins. In this study, SVD aimed to extract dynamic biometric signals to visualize blood fluidic behavior, which is difficult to achieved using the standard Fourier method, as will be discussed in the section 3.2.
To process data in the SVD method, consecutive stacks of images from a 3D array (n x , n y , n t ) are reshaped into 2D arrays (n x × n y , n t ). Because of this transformation, each column contains all pixel information from the same time.

SVD
To understand the SVD application in this study, we present the definitions illustrated in Fig. 1(c) [24]. Let matrix A of size (n x × n y , n t ) be a dataset with time step n t and n x × n y pixels of data. According to the definition of SVD, matrix A can be decomposed into three matrices as follows; * denotes the conjugate transpose: where matrix U of size (n x × n y , n t ) denotes the spatial eigenvector (U 1 , U 2 , . . . U k , . . . , U 1000 ), representing the specific spatial distribution of energy states of the signal corresponding to each frequency resolution. Matrix V of size (n t , n t ) denotes the temporal eigenvector (V 1 , V 2 , . . . V k , . . . , V 1000 ), representing the signal temporal variation corresponding to each spatial eigenvector in matrix U. Both U and V are weighted by the corresponding singular values λ k , which represent the energetic information in each eigenvector and are elements of square diagonal matrix S of size (n t , n t ). Additionally, singular values λ k in diagonal matrix S are sorted in descending order. In this study, high-energy states containing significant noise will be concentrated in the first few singular values λ k . Although subsequent singular values are loaded with less noise, the signal energy contributing to blood flow also decreases, thus making the noise and signal strength comparable and making it difficult to extract blood information. In summary, the SVD features are the signal energetic information provided by singular values λ k in Fig. 1(d), temporal behavior of signal fluctuation given by matrices V and singular vectors V k in Fig. 1(e), and spatial distribution encoded in matrices U and singular vectors U k in Fig. 1(f).
This study performed SVD to extract signals in the frequency ranges related to myogenic activity, neurogenic activity, and endothelial activity from other unwanted signals, because these types of activity contribute higher energies to blood flow in the vessels compared to respiratory and cardiac activities [17].
Similar to the Fourier method, for a measurement of the sampling frequency of Fs = 55.56 Hz and time window of 18 s (n t = 1000 frames), we can assume that each frequency window ∆f = 0.5 × Fs / n t = 0.027 Hz. Such frequency resolution provides sufficient precision to extract blood signals that contribute to myogenic activity. However, note that, in SVD, the temporal eigenvectors corresponding to the frequency band are ordered in the singular values (λ k ) of the signal energetic contributions, rather than the frequency order.
The choice of a threshold to truncate the desired signals for blood flow motion visualization is based on the theory above, and the results presented in Fig. 1(e) and 1(f). According to Fig. 1(e), temporal eigenvectors V 2 -V 6 filter out the DC portion of the signal contributed by tissue or clutter; hence, we can select V 2 for the bottom border of truncation. Additionally, according to Fig. 1(f), U 2 yields a clearer vein image (i.e., three vascular branches indicated by a red square) than U 1 , because U 2 reconstructs a vein image by extracting the signal in the endothelial activity frequency range (0.0095-0.021 Hz). This is consistent with the temporal eigenvector variation that is shown in Fig. 1(e), where V 2 exhibits a slow oscillation frequency over 1000 frames (18 s).
For the upper cutoff frequency region (approximately 0.145 Hz) of myogenic activity, we chose the fifth temporal eigenvector (V 5 ) because its upper cutoff frequency is 0.135 Hz. Therefore, we set n t = 2 and 5 as our lower and upper temporal eigenvector thresholds, respectively, to truncate the singular value matrix as our blood vessel visualization results.
Finally, we used the Eq. (2) to reconstruct a time-resolved movies of vascular myogenic activity. All computations were performed using MATLAB R2020a.

PPG signal decomposition analysis results
We first present the results of the reconstructed blood flow images (1000 frames, 55.56 fps) from each individual decomposed frequency band above the DC component. In Fig. 2(a), the reconstructed consecutive images exhibit no temporal variation. This is because the first spatial eigenvector (U 1 ) and first temporal eigenvector (V 1 ) only extract and reconstruct the lowest frequency component of the PPG signal. We can infer that this component of PPG contributes to the signal DC offset with the highest signal energy information among the singular values (λ 1 ). In Figs. 2(b) and 2(c), stable temporal variation can be observed with the finger vein disappearing and reappearing again. One can infer that this is caused by the endothelial, neurogenic, and myogenic activities of the finger. Based on this result, the second and third spatial (U 2 , U 3 ) and temporal (V 2 , V 3 ) eigenvectors mainly contribute to the localization and reconstruction of images for vein deformation during blood flow. Figure 2(d) shows that the images reconstructed from the fourth eigenvectors (U 4 , V 4 ) suffers from severe noise. This result is consistent with the result presented in Fig. 1(e), which shows that V 4 exhibits extremely fast oscillation. We suggest that this signal does not correspond to a biometric signal, and is generated by the imaging system (e.g., dark current, pixel crosstalk, and parasite capacitance noise). Figure 2(e) shows that two areas (indicated by arrows) fluctuate alternatingly. One can observe the fifth temporal eigenvector (V 5 ) corresponding to the upper region of the physiological bands of myogenic activity (0.1-0.13 Hz) as shown in frequency analysis result in Fig. 2(g), which contributes more energy to blood flow than vein deformation. In other words, we suggest that the fifth eigenvectors (U 5 , V 5 ) contribute to the blood flow portion of myogenic activity. Figure 2(f) depicts clear and stable respiratory-like temporal variation of 1 Hz in frequency analysis result shown in Fig. 2(g) during the time window of 18 s (1000 frames). In other words, we have proved that SVD is capable of extracting PPG components from respiratory physiological bands [25,26].

Imaging blood vein myogenic activity
We have demonstrated the superiority of the SVD filter over conventional Fourier bandpass filtering, as shown in Fig. 3. The consecutive frames in the original images without the filtering process, with Fourier bandpass filtering, and with SVD filtering are presented at the left top, left middle, and left bottom of Fig. 3, respectively. In Fig. 3(a), the finger vein biometric images are dominated by static signals, which prevent dynamic blood fluidic activity from being revealed. In Fig. 3(b), the conventional Fourier bandpass filtered images are strongly affected by noise, demonstrating Fourier filtering's weakness at extracting such a narrow frequency band signal (<2 Hz) to visualize blood flow information. In Fig. 3(c), SVD filtering isolates the vascular activity signal to highlight some vein branches that cannot be observed in Fig. 3(a), and detects myogenic activity obtained from vein endothelial activity and myogenic-activity-like behavior in the reconstructed images. These examples demonstrate that the SVD filtering method is powerful for isolating lowfrequency vascular activity signals from the DC components of biometric signals with higher performance than the standard Fourier method. Based on the definition of SVD for bio-signal extraction, clutter contributes higher spatial coherence related to singular values of higher energy, which are removed during the truncation process applied to the singular value matrix. In this study, we only focused on extracting vascular activity signals from the second-to fifth-lowest singular values to obtain vascular visualization with the highest contrast by removing clutter contributions and minimizing the effects of noise in the unwanted frequency band.

Imaging skin microcirculation
In this study, we also analyzed the capabilities of SVD to decompose vascular activity further, by performing reconstruction with an appropriate singular value ratio to improve the quality of visualizing blood fluidic behavior. First, we can conclude that the fourth eigenvector contributes to noise, rather than to the biometric signal. Therefore, we can filter it out by setting λ 4 = 0 during the reconstruction of the blood fluidic images. Similarly, the sixth eigenvector exhibits a signal oscillation frequency of approximately 1 Hz, which is far from the myogenic activity signal frequency; hence, we also filter out λ 6 to avoid a decrease in the quality of blood flow imaging reconstruction. In other words, according to our theory, perfusion images reconstructed from the second, third, and fifth eigenvectors should exhibit the clearest vascular activity, because the signal components of the fourth and sixth eigenvectors suffer from noise and unrelated signals that may be induced by factors such as the imaging system.
As shown in Fig. 1(d), after SVD, the singular value energy ratios of the second to fifth eigenvectors are in descending order of (17, 3, 3, 1) originally. By filtering out the fourth eigenvector, and then, modifying the singular value ratios for the remaining eigenvectors, the vascular activity in the frequency band can be enhanced and focused. For example, as shown in Fig. 3(c) (d) (e), by adjusting the ratios from (17, 3, 0, 1) to (17,3,0,6), the visualization of vessel pumping-and swelling-like behavior is weakened. In contrast, blood fluidic activity in the vertical direction is strengthened as the blood flows back and forth from the palm to the fingertip. This is because, as shown in Fig. 2(c), V 5 provides a stable signal that tracks blood flow, but contributes less energy (λ 5 ) than the endothelial activity signal and myogenic activity signal (λ 2 , λ 3 ) (in Fig. 1(d)). Therefore, by artificially adjusting singular value energy ratios, myogenic and blood fluidic activities can be made comparable in the observations. Another interesting finding in the time-resolved movies in the Supplement 1 corresponding to Fig. 3(d) is that, during the time window of 6-10 s, the vein is observed to shake horizontally during the transportation of blood (see Visualization 1). This demonstrates the powerful capabilities of SVD for extracting and reconstructing very small signals.
Additionally, by introducing the concept of phase shifting for blood volume oscillation, we can define the highest value as 90°and lowest value as −90°, because the extracted signal frequency band is so narrow that we can assume that it is close to a sinusoidal wave [27]. Similar results for blood movement can be observed easily by following the shifting of red regions, as shown in Fig. 3(f) (see Visualization 2).

Imaging palm skin microcirculation
In this section, we discuss the application of the SVD technique adopted in this study to enable the detection and visualization of skin microcirculation in a conformable imager [23]. By directly contacting the palm, as shown in Fig. 4(a), our conformable imager can directly capture hand palm vein images from static biometric signals without using an external optical system (e.g., lens or alignment stages) (see Visualization 3). By applying the SVD method, microcirculation under the palm skin can also be observed, as shown in Fig. 4(c). Two branches (indicated by green and blue arrows) of the vessels are visualized by the second and fourth spatial eigenvectors (U 2 , U 4 ), whereas clutter signals and noise are contained in the first, third (indicated by red arrows), and fifth spatial eigenvectors (U 1 , U 3 , U 5 ). For the perfusion images reconstruction in Fig. 4(d), extracting and combining only U 2 and U 4 yields the best visualization quality for a vessel image (see Visualization 4) that rejects clutter signals and noise (indicated by red arrows). This mechanism can be described as follows. Direct contact with an object overcomes the biometric signal loss issue during signal transport in biomedical sensing/imaging systems, and the high-quality visualization of static biometric signals indicates the possibility that dynamic biometric signals are included. By adjusting the sampling time window, we can adjust the sampling frequency window size in a manner similar to the Fourier method.
However, the Fourier method is not efficient for imaging processing, particularly in highresolution situations, because of its massive processing times (n x × n y iterations of the filtering process) and band filtering in an extremely low-frequency window size (approximately 0.1 Hz) incurs significant system requirements. By using the SVD method implemented in this study, we can halve the sampling frequency window size from 0.027 Hz (n t = 1000) to 0.014 Hz and 0.007 Hz (n t = 2000, 4000) by doubling the sampling time window, which enables us to double the precision of extracting biometric signals. By making such adjustments, one can easily visualize palm dynamic biometric signals at high resolution without incurring additional computational burdens on the system, unlike the standard Fourier method. The application of SVD in lens-less imaging provides a new approach for pattern recognition by determining the threshold in low eigenvectors owing to the signal energetic contributions [28,29], which is different from the traditional Fourier method (i.e., frequency), deep learning/machine learning (i.e., back-projection networks) [30,31]. The reason for its effectiveness when used in PPG is the strong temporal coherence property of a PPG signal in the low frequency domain (<2 Hz). Each signal energetic value corresponds to information contained in each portion of a biometric signal including blood flow induced relative vessel volume changes in different frequency windows, tissue clutter, and noise. While tissue clutter signal and noise do show higher singular energetic value than certain blood flow signals in some cases, it can be easily recognized by spatial vectors and filtered out by simply setting the singular value to zero during the image reconstruction process. Such processes can hardly be achieved with the Fourier method if an unwanted signal overwhelms neighboring target signals in the frequency domain. In other words, the ease of determining a singular value for image reconstruction provides higher efficiency than the Fourier method. Another important property of the SVD application in this study is the visibility of image reconstruction data processing steps (i.e., white box), which is different from that in the case of the deep learning/machine learning approach (i.e., black box). In this way, we can easily recognize different visualized biometric components from the image reconstruction of each singular value without using reference data (i.e., ground truth material), which simplifies the data preparation process and could be of high practical interests for developing remote PPG monitoring/imaging system [32].

Conclusion
In a PPG signal, the signal component that contains blood flow information is extremely close to the DC component in the frequency domain. Therefore, the capability of filtering a PPG signal in a 0.02-Hz window efficiently is crucial for reconstructing skin blood flow activity images. We demonstrated that the PPG signal components from 0.02 Hz to 0.1 Hz contribute necessary signal information for visualizing below-skin blood flow. Such narrow frequency band filtering cannot be achieved effectively using the standard Fourier method. Additionally, by analyzing decomposed spatial and temporal eigenvector matrices, we identified noise components that are caused by the system or other unknown sources. By filtering out unnecessary PPG signal and noise components and adjusting the singular value energy ratios, the reconstruction of high-quality blood perfusion images can be achieved for further applications such as blood flow monitoring.
Most importantly, our method exhibits much higher efficiency than the standard Fourier method. For the Fourier method, each pixel with n t consecutive time steps is used for the fast Fourier transform (FFT), bandpass filter, and inverse FFT processes. This required a tremendous amount of time (8 h) in our study. In contrast, the SVD method only required approximately 3 s to extract the required signals, thus taking a time span of less than 1 min for reconstructing perfusion images. This result demonstrated that our simple and easily implementable algorithm is a promising method for achieving below-skin blood flow perfusion imaging.