Blood pulsation imaging (BPI) is a non-invasive optical method based on photoplethysmography (PPG). It is used for the visualization of changes in the spatial distribution of blood in the microvascular bed. BPI specifically allows measurements of the relative phase of blood pulsations and using it we detected a novel type of PPG fast waveforms, which were observable in limited areas with asynchronous regional blood supply. In all subjects studied, these fast waveforms coexisted with traditional slow waveforms of PPG. We are therefore presenting a novel lock-in image processing technique of blood pulsation imaging, which can be used for detailed temporal characterization of peripheral microcirculation.
© 2014 Optical Society of America
The non-invasive optical method, photoplethysmography (PPG), is a traditional and widely used approach allowing the detection of blood volume changes in the peripheral tissues of the body [1,2]. Application of the PPG technology, based on pulse oximetry, allows noninvasive measuring of oxygen saturation in blood. Implementation of PPG is very simple since it requires only two optoelectronic components: a light source for illumination and a photoreceiver for measuring the light intensity either transmitted or reflected from part of a living body [2–4]. The output signal from the photoreceiver varies in time proportionally to the changes in peripheral blood volume at the site of measurement showing modulation corresponding both to heartbeats and respiratory rate [2,4–8].
Conventional PPG devices such as the pulse oximeter provide measurements only in a single point. To overcome this limitation several research groups reported the development of an imaging system for remote oximetry based on the reflection-mode PPG and operating under illumination at two or more wavelengths [8–10]. In these systems a two-dimensional (2D) photosensitive matrix is used instead of a single photoreceiver. We proposed recently a PPG imaging of high spatial resolution, called Blood Pulsation Imaging (BPI), which exploits illuminations at the single wavelength . Data processing of digital images in BPI is based on the detection of time-varying modulation of each pixel in the recorded video frames synchronously with the heartbeats [11–14]. The most important feature of this approach is that, apart from traditional measurement of amplitudes, BPI allows to visualize the distribution of the relative phase of the blood pulsations. As we more recently showed , the relative phase of blood pulsations is a very important channel of biomedical information. It can be used for instance, as a biomarker for dysfunction of the autonomic vascular control which may occur in such neurological diseases as migraine . However, as the BPI method is new, the reliability of the blood pulsation phase (BPP) measurements needs to be further improved.
Therefore, the aim of this study was to elaborate a criterion for the quantitative estimation of data processing resulting in more adequate calculation of the 2D distribution of the blood-pulsation phase. We found that the shape of the PPG waveform is strongly variable over the palm and wrist areas in all studied subjects, which leads to an ambiguity in the measurement of the relative phase of blood pulsations. Moreover, in areas of the subject’s body with asynchronous regional blood supply we detected new fast type waveforms of PPG, which are different from the traditional slow type waveforms. To take into consideration this uneven distribution in the shape of the waveform, we propose a new algorithm for data processing, which includes a proper choice of the reference area and dual definition of borders of cardiac cycles. Two shapes of the reference function (sinusoidal and triangle) were used for the calculation of BPP maps. By comparing the calculated spatial distributions of BPP with the relative phase between PPG waveforms measured in different regions of interest (ROI), we found that the shape of the reference function could be used for a more reliable BPP mapping.
2. Methods and materials
2.1 Experimental setup
The experimental setup for video recording required for BPI is shown in Fig. 1. The video recordings of the subject (part of the palmar and wrist) were done by a digital monochrome 8bit CMOS camera (EO-1312, Edmund Optics, Barrington, NJ, USA) with attached lens (18 - 108 mm FL, Canon, Tokyo, Japan) in a reflection mode (with a minimal angle between the optical axis of the lens and illuminator). The video frames with focused images of the illuminated subject were downloaded to a personal computer via the universal serial bus (USB). The light source was powered from another USB port providing subject’s illumination with very stable light intensity. The frame rate of video recording was 50 frames per second (fps) and the frame size was 640 × 480 pixels.
To illuminate the subject we used conventional light-emitting diode (LED) (H2A3-530, Roithner Lasertechnik GmbH, Wien, Austria) operating at the wavelength of 525 nm. We used illumination by the green light because the maximal temporal modulation of the light intensity back reflected from the vascular bed was observed in this range of wavelengths . The optical power of the LED was 80 mW, the wavelength bandwidth was 60 nm. The illuminator provided uniform illumination of the area of 10 × 10 cm2 at a distance between the subject and video camera of about 50 cm. Note that during video recording the lens diaphragm and the camera exposure time were adjusted so that the recorded frames did not contain saturated pixels (the maximal pixel value is 255 for 8-bit camera) but were as high as possible. During the experiment, the ambient illumination was much lower than that provided by the illuminator.
The illuminating light was linearly polarized by means of a polarizer P1 (Fig. 1) attached to the illuminator. The polarizer P2 was attached to the camera lens so that its transmission axis was oriented orthogonal to the polarization vector of the illuminating light. By using crossed polarizations we reject the light reflected from the skin surface-air interface.
In order to estimate the signal-to-noise ratio (SNR) of the calculated parameters of blood pulsations we used an approach in which the data processing was simultaneously applied to the subject’s area and to an inanimate object. To this end we attached a piece of cardboard (Fig. 1) close to the palm so that they did not contact mechanically and the reflected light from the palm did not touch the cardboard.
2.2 Experimental protocol
Before measurement the subject placed his palm on the sponge support. Taking a comfort position the subject was asked to avoid hand movements during one minute of video recording. The measurements were carried out with seventeen healthy subjects. Age of the subjects was from 28 to 89 years. The group consisted of 9 women and 8 men. All subjects gave their informed consent of participation in the experiment and in the publication of the results in the written form.
This study was conducted in accordance with the ethical standards laid down in the 1964 Declaration of Helsinki. The study plan was approved by the Research Ethics Committee, Hospital District of Northern Savo.
2.3 Previous algorithm of data processing
The recorded video frames were processed by using custom BPI software implemented in the MATLAB® platform. Previous algorithm of BPI consisted of two steps: formation of the reference function and evaluation of the correlation matrices contained the 2D distribution of the blood pulsation amplitudes (BPA) and phase (BPP) for each cardiac cycle.
In our previous algorithm of BPI data processing  we created the reference function from the PPG waveform obtained by spatial averaging of the pixel values within a large (50% of the frame size) reference area in the center of the recorded frame. The resulted waveform was divided into even parts with the length equal to the mean period of cardiac cycles. Each part of the PPG waveform was approximated by a sinusoidal function of the same length. Natural variability of the heartbeats rate  was taken into account by matching the minimum of the sinusoidal reference with the respective minimum of the PPG waveform in each cardiac cycle [11–13]. The resulted combined sinusoidal waveform was further used as a reference function needed to calculate the correlation matrix.
In this report, the above described algorithm for formation of the reference function was used only for preliminary evaluation of the BPA and BPP maps to choose the optimal reference area. Thereafter, the PPG waveform was calculated by averaging the pixel values within the chosen area. In contrast with the previous algorithm here the reference function was not adjusted to be in the phase with the individual PPG pulse. It was only the time duration of a cardiac pulse evaluated from the temporal coordinates of the extremes of the PPG signal, which is needed for the formation of the reference function.
3. Algorithms for synchronous detection of blood pulsations
3.1 Preliminary BPI
In the first step of the new algorithm we calculate initial BPA and BPP maps using previous BPI algorithm with the sinusoidal reference function defined from the PPG waveform measured in the large reference area. An example of typical BPA and BPP maps is shown in Fig. 2. It is seen that both the amplitude and phase of blood pulsations is unevenly distributed over the subject’s area, which was previously reported in our papers devoted to BPI [11–14]. One can clearly see several areas with the big phase difference coded by different colors (blue-green versus red) in Fig. 2(b) (BPP map). It means that the blood is asynchronously supplied to these areas. Existence of such asynchronicity was reported in our recent paper . However, we did not previously check the shape of the PPG waveform in these areas.
Unexpectedly, we found that PPG waveforms are not the same in different areas of a subject’s body. Three typical temporal shapes of the waveform are shown in Fig. 3. They were obtained by averaging the pixel values within regions of 10 × 10 pixels situated in areas with elevated amplitude of blood pulsations. Positions of chosen ROIs are marked by the squares in Fig. 2. ROI-1 and ROI-2 were placed in the areas with asynchronous blood pulsations. Maximal and minimal values of the PPG signal in the reflection mode of the experiment are associated with the diastole and systole phases of the vascular system, respectively [9,10]. The waveform in Fig. 3(a) is characterized by the slower systole-to-diastole transition, i.e. by longer rising time of a PPG pulse. Here we call this type of transition as a slow type. This type of the PPG waveform was observed and reported in all numerous papers devoted to photoplethysmography [9,10,17–19]. Note that Fig. 3 shows the optical intensity which changes in the opposite direction of the blood pressure: the larger blood volume corresponds to the larger attenuation of light , and consequently, to the smaller intensity returned to the CMOS camera. Therefore, the rise in the blood pressure corresponds to the fall in the light intensity. Similar waveforms of slow diastole were observed in a larger area of all subjects studied in our experiment. However, an inversed waveform with the sharp systole-to-diastole transition (fast type) [Fig. 3(b)] was observed in areas with asynchronous blood pulsations, which are visualized by the hues of red color in the BPP map in Fig. 2(b). Notably, such fast type waveforms were never observed before. In addition, there are places (for example, ROI-3 in Fig. 2) in which the PPG waveform has an intermediate symmetrical shape [Fig. 3(c)].
Observed diversity of PPG waveforms shows that the random choice of the reference area (as it was in the previous BPI algorithm) could lead to incorrect calculation of BPA/BPP maps if the areas with different PPG waveforms are joined in a large size ROI from which we generate the reference function. In addition, some regions of microcirculation with atypical fast type PPG waveforms would be overlooked. Since the reference function is a key parameter of any algorithm of synchronous detection, one should carefully choose both the position and size of the ROI for reference function formation.
3.2 Choosing of the reference area
On the one hand, the reference area should be of small size to exclude an interference of adjacent zones with inversed waveforms because it can result in incorrect estimation of the cardiac cycle duration. On the other hand, time-varying modulation of pixel values is very small (typically, less than one per cent), which requires data averaging over large number of pixels to increase SNR of the pulsatile component. In our experiments, the reference-area size of 10 × 10 pixels was found to be the compromise. This size of the reference area corresponds to the physical size of 3 × 3 mm2 as measured at the subject’s hand. For such a ROI, the pulsatile component of the PPG waveform has high enough SNR as seen in Fig. 3, and this size fits the minimal dimension of areas in which the blood pulsates with the same phase as it can be estimated in Fig. 2(b).
Position of the reference area was chosen considering the above requirements for correct estimation of duration of each cardiac cycle: it should be placed inside the areas (i) with bigger amplitude of blood pulsations to achieve the high SNR, and (ii) with uniform phase distribution to exclude the interference of different waveforms. In the example of preliminary BPA/BPP maps shown in Fig. 2, the position of ROI-1 satisfies these conditions.
To estimate, which shape of the PPG waveform is more common in the palm and wrist, we carried out a statistical analysis within the cohort of 17 subjects participated in our experiment. For all studied subjects the slow PPG waveform [with the sharp diastole-to-systole transition, such as in Fig. 3(a)] was observed at the largest part of the palm and wrist. Therefore, it would be useful for further analysis to place the reference area in the zone with the dominating shape of the PPG waveform.
3.3 Reference functions formation
The PPG waveform extracted from the chosen reference area was further used for the formation of the reference function within each cardiac cycle. Considering the diversity of the PPG waveforms shown in Fig. 3, mapping of the relative phase of blood pulsations may seriously depend on the shape of the reference function. Moreover, as we show below, the diversity of PPG waveforms leads to ambiguity of the relative phase magnitude because the time difference between the moments of reaching the systole phase in two spatially separated points may be different from this for the diastole phase. To estimate how the reference function does affect BPP and BPA maps, we calculated these maps using two types of reference functions: symmetrical sinusoidal and triangle with variable asymmetry. Formation of any reference function was started from defining the beginning and duration of each cardiac cycle, which was done by finding positions on the time scale of either maxima or minima of the PPG waveform extracted from the reference area. When the cardiac cycle is defined by the minima of the PPG waveform, the relative phase of blood pulsations in different pixels is calculated as the difference in time between the moments at which the vascular system reaches the diastole phase in these pixels. Alternatively, it is mapping of the time difference in reaching the systole phase if the cardiac cycle is defined by the maxima of the PPG waveform. Consequently, for each type of the reference functions we define two mathematically different shapes, S and D, for calculations of spatial distributions of time differences for the systole and diastole phases, respectively.
3.3.1 Sinusoidal reference function
Schematic of formation of the sinusoidal reference function of D-type, , is shown in Fig. 4(a). The blue curve is the PPG waveform from the chosen reference region of 10 × 10 pixels within j-cardiac cycle. Circles show the beginning and the end of each cardiac cycle corresponded to D-type of the reference function, which temporal coordinates, and , respectively. Duration is varying from one cycle to another and its particular magnitude is depicted above the graph. Here j is a sequential number of the cardiac cycle. The real part of the reference function was defined so that its minima coincide with the beginning and end of each cardiac cycle. The imaginary part was just 90°-phase shifted from the real part. Therefore, the reference function of the k- cardiac cycle is a complex one:
The normalization coefficient was defined separately for each cardiac cycle from the condition that the convolution of the PPG waveform with the respective reference function should be equal to the peak-to-peak amplitude of:Fig. 4 with its scale on the right while the scale of the PPG waveform is shown on the left.
The sinusoidal reference function of S-type was formed similarly as in Eq. (1):Eq. (1)] which has the negative sign. The normalization coefficient was defined as
3.3.2 Triangular reference function
The scheme of formation of the triangular reference function of D-type, , is shown in Fig. 5(a). The blue curve is the same PPG waveform as in Fig. 4 with beginning and the end of each cardiac cycle corresponded to the D-type reference function. The base of the reference function is a triangular shape with three vertices which positions were defined inside j- cardiac cycle as
Thereafter, the reference function based on the shape of the defined by Eq. (5) for j-cardiac cycle can be written as
The normalization coefficient was defined separately for each j-cardiac cycle from the condition that the convolution of the PPG waveform () with the respective triangular shape should be equal to the peak-to-peak amplitude of:Eq. (2), and is one of the shapes , obtained by the matching of the temporal coordinate with the maximum of the PPG waveform in the j-cardiac cycle, which corresponded to the best correlation of shape with the PPG waveform. The summation is done over the cycle duration ().
The S-type triangular reference function was defined as:Eq. (4), and is one of the shapes obtained by matching the temporal coordinate with the minimum of the PPG waveform in the k-cardiac cycle. The summation is done over the cycle duration.
Triangle reference functions and are also alternating-sign functions with the mean equal to zero. In each graph of Fig. 5, the scale of the reference function is on the right while the scale of the PPG waveform is shown on the left.
3.4 Image preprocessing
Before calculating BPA/BPP maps we applied digital processing of recorded frames to compensate for accidental motion of subject’s arm. The custom software for this aim was implemented using the MATLAB® platform (of the Mathworks). The whole frame was split into small regions for detecting the mean gradient of the pixel values within the subject’s area in the frames sequence starting from the first frame of each cardiac cycle. Once the motion of any region has being detected, we stabilized exactly the region with moved gradient while the surrounding area remained stable and hence not modified.
After motion compensation the pixel values in the sequence of the frames during a cardiac cycle comprise two components: a pulsatile component attributed to interrogated blood volume, and a slowly varying baseline (DC) caused by respiration and sympathetic nervous system activity. Before applying the BPI algorithm we filtered out the DC component by calculating
The video frames modified according to Eq. (10) for each subject were further processed for calculation of blood pulsation parameters with four reference functions , , , and , resulting in four arrays of the correlation matrices.
3.5 Calculation of correlation matrices
3.5.1 BPI with sinusoidal reference function
Synchronous detection of the pixel value modulation using the sinusoidal reference function of D-type, , defined by Eq. (1) results in a correlation matrix evaluated for j-cardiac cycle asEq. (10)] spatio-temporal distribution of pixel values in the frames during the j- cardiac cycle in which the reference function is defined.
The correlation matrix calculated in the BPI algorithm with S-type of the reference function was similarly defined asFig. 4). matrix contains information about the moments of reaching the diastole phase, whereas matrix allows estimation of the moments for achievement the systole phase. Correlation matrices in Eqs. (11) and (12) are complex because of complexity of the and [see Eqs. (1) and (3)]. Since the sinusoidal reference function is defined in the quadrature, the modulus of the correlation matrices and describes the spatial distribution of blood pulsation amplitude (i.e. BPA map), while the argument of matrices value corresponds to the relative phase of blood pulsations (BPP maps). The values of the BPA maps represent the amplitudes of the pixels variations in and , measured in per cent of values in and , respectively. The values of the BPP maps represent the time difference in reaching either the diastole or systole phase in the pixel with coordinates and the reference area (ROI-1). Values of BPP maps are measured in degrees after normalization of the time difference with the duration of the cardiac cycle.
3.5.2 BPI with triangular reference function
Calculation of the blood pulsations parameters using triangular reference function of D-type, , was executed for the same recorded video frames as with the sinusoidal reference . First we find the temporal position at which each pixel in these frames reaches its maximal value during the j- cardiac cycle. To this end we calculate a series of the correlation matrices of the filtered-out spatio-temporal distribution given by Eq. (10) with different triangular reference functions of . The function is defined by Eq. (5) with the position of its maximum coinciding with one of the recorded frames within the j- cardiac cycle. Formation of the series of reference functions is schematically depicted in Fig. 6 by red dotted lines. The number N of these functions in the series is equal to the number of the recorded frames within the j-cycle minus one. The m-correlation matrix in the series is given by
By finding the index m at which one of the matrices has the maximum, we define the moment when the value of the pixel with spatial coordinates of reaches the diastole phase. An example of a pixel value modulation is shown in Fig. 6 by the blue line. In this particular case, the correlation matrix with the reference function is the maximal one. Therefore, defines the diastole moment in this case. Note that the moment of reaching the diastole may vary with the coordinates . The phase of blood pulsations was estimated as the ratio calculated in degrees in respect to the average phase in the reference area (ROI-1).
Second, we estimate the amplitude of blood pulsations in the pixel with the coordinates as the maximal correlation value in the series of . This allows us to calculate the BPA map for the j-cardiac cycles.
For calculations of BPP/BPA maps using S-type of triangular reference functions we defined the borders of each cardiac cycle by the maxima of the PPG waveform in the reference area, whereas the series of the references functions has variable position of the triangle minimum unlike for the [see difference between (a) and (b) of Fig. 5]. Therefore, like with the sinusoidal reference functions, a BPP map calculated with the D-type reference represents spatial distribution of the time difference for reaching the diastole phase while with the S-type we calculate the time difference for the systole phase.
4. Blood pulsation imaging
4.1 Fast type PPG waveforms
By applying the algorithm of video-data processing described in the previous section we calculated BPP and BPA maps for each cardiac cycle of every subject. Typical distribution of the blood pulsation phase and amplitude calculated with D-type sinusoidal reference function is shown in Fig. 7. Both BPP and BPA maps are overlaid in this figure with the initial camera image to provide an anatomic reference of the respective subject. As one can see in Fig. 7(a), there are adjacent limited areas in which the blood pulsates asynchronously, i.e. with big a relative phase. These asynchronous areas are marked in the BPP map by the red and blue colors. They are usually situated at the wrist of each subject alongside the radial and ulnar arteries. However, sometimes these areas are extended perpendicular to the artery [see Fig. 10(b)]. Moreover, in 11 of 17 subjects the asynchronous areas were also observed in the middle of the palm. We found that the atypical fast type of the PPG waveform [such as shown in Fig. 3(b)] was always observed in the areas marked in our BPP maps by the red color.
The quantitative evaluation of these areas was done by calculating the amount of the pixels inside the area where the BPP values were more than 90°. The partial area with the fast type PPG waveform was calculated as per cent of the whole area of the subject's tissue within the video frame. The partial area, from which the new fast type of the PPG waveform was observed, was variable from one subject to another. The results of statistical analysis of this partial area distribution in the subjects of different age are presented in Fig. 8. The averaged partial area within all the subjects was 3.2% ± 1.3% (n = 17, P<0.0001, one-sample t test).
The fast-type partial area varied from 0.7% to 5.5%, where the highest value was observed in the 28-years old subject, while the lowest one – in the 64-years old subject. However, no correlation of the fast-type area was found either with age or gender of the studied subjects (P > 0.05, n = 17).
4.2 BPI with four different reference functions
Proper choice of the reference area allowed us to increase the reliability of estimation both the beginning and duration of each cardiac cycle. These two parameters are sufficient for formation of four different reference functions: sinusoidal and triangular of both D- and S-type defined in Eqs. (1), (3), (6) and (8). Using these references we calculated BPA and BPP maps applying the algorithms of Sect. 3.5 for every cardiac cycle. An example of both BPA and BPP maps for a random cycle of one of the subjects is shown in Fig. 9. The reference area of 10 × 10 pixels from which we extracted the PPG waveform for formation of the reference functions was chosen in the ROI-1. The phase of blood pulsations of each pixel in the BPP maps was calculated relatively to the phase of pulsations in the reference area (ROI-1).
At first glance, all distributions of blood pulsations amplitude in the upper row of Fig. 9 seem to be almost identical. However, there are few differences indicated by red arrows in the BPA maps. For example, in the top left corner of the maps (a) and (c) calculated with the D-type reference functions an area of the elevated amplitude is weaker and shorter than it is in the maps (b) and (d) of the S-type. In addition, the bright spot which is clearly seen in the lower left corner of maps (a), (b) and (d), completely disappears in Fig. 9(c) calculated with reference function. In BPP maps shown in the lower row of Fig. 9 we observe not only variations of the relative phase magnitude (which is visualized by variation of pseudo colors), but also a change of the topology in maps calculated with D-type reference functions compare with S-type in the region marked by the blue arrow.
In order to understand what shape of the reference function describes more reliably the microvascular dynamics, the corresponded criterion is required. Such a criterion is proposed in the next section.
4.3 PPG signal as a criterion for BPP maps reliability
Since there is no alternative technique which could measure the relative phase of blood pulsations of a live body, we decided to use raw PPG signal measured in different areas as an adequacy criterion for data processing with different reference functions. To this end, we calculated the PPG signal in ROIs of 10 × 10 pixels placed into several areas with high amplitude of blood pulsations. Positions of such ROIs are shown in Fig. 10(a), 10(b) using BPP maps of two subjects. In both subjects the PPG signal within ROI-1 was chosen for estimation of each cardiac cycle beginning/duration for formation of the reference function. Time traces of the PPG signals are shown in adjacent graphs in Fig. 10(c), 10(d) within ROIs 1-3 to compare their relative evolution from one cardiac cycle to another. It is seen that the blood pulsates in ROI-3 almost in the counter phase with ROI-1. However, due to different shape of the PPG waveforms from these ROIs, quantitative estimation of the relative phase between them is difficult. Nevertheless, we can always estimate the time difference between the moments of systole or diastole from the pair of PPG waveforms within the reference area (ROI-1) and ROIs 2-3, which is convenient to recalculate to the phase difference by dividing the former by the cycle's duration.
As it follows from the time traces in Fig. 10, the time difference between the moments of reaching the diastole phase in ROI-1 and ROI-3 (measured between the positions of the maxima in blue and red curves) is not the same as for the systole phase (between the positions of minimum). Therefore, we obtained two types of the relative phase difference: for the systole and diastole. These two types should be separately compared with the phase difference calculated by BPI algorithm with D- and S-types reference functions, respectively.
4.4 Blood pulsation phase comparison
Both the BPI algorithm and raw PPG waveforms allows us to estimate the relative phase of blood pulsations for every cardiac cycle. However, these quantities vary in time because of the heart rate variability which is inherent in any subject . Therefore, for correct comparison of BPI results we first averaged the values of the BPP maps within the same ROIs as for extracted PPG waveforms, and then averaged them in time over 60 cardiac cycles. In all cases (for systole and diastole) the phase of blood pulsations was calculated in respect to the phase measured in the chosen reference area (ROI-1). For each subject we chosen 28 ROIs situated in areas with high enough amplitude of blood pulsations. The size of each ROI was 10 × 10 pixels. Examples of the ROI choice for two different subjects are shown in Fig. 11(a) and 11(b) where we show only 7 of 28 chosen ROIs for clarity of presentation.
We found that the BPI algorithm with the sinusoidal reference function gives more reliable values of the D-type relative phase in 93, 71, 71, 64, and 50 per cent of ROIs for five studied subjects. Representative example of D-type-phase comparison is shown in Fig. 11(c). The green bar shows the phase difference calculated from row PPG waveforms while the blue and red bars are the relative phase calculated with the sinusoidal and triangular reference functions, respectively. Whiskers show the standard deviation of the magnitude during averaging over 60 cardiac cycles. It is seen that the phase calculations obtained with the function is much closer to the estimations from PPG signals than that with the function.
In contrast, for the S-type relative phase, calculations with the triangular reference function give more reliable results in 75, 64, 61, and 57 per cent of chosen ROIs for four of five studied subjects. In one subject, phase calculations with the sinusoidal reference function were more reliable even for the S-type (64%). An example of S-type-phase comparison for the case of more reliable calculations with is shown in Fig. 11(e). Fragments of the PPG time-traces measured in ROI-1 (blue lines) and ROI-2 (red lines) for subject (a) and (b) are shown in Fig. 11(d) and 11(f), respectively. Diversity of PPG waveforms in different regions of the body and in different subjects does not allow us to find a universal reference function for reliable estimations of the relative phase distribution.
4.5 Pulsation amplitudes comparison
We also used the PPG waveforms measured in different regions of subject to estimate the BPI-algorithm reliability in calculation of the blood pulsation amplitude with sinusoidal and triangle reference functions of D- and S-types. Special attention was paid to comparison of BPA mapping in regions with asynchronous blood pulsations, i.e., with big relative phase of pulsations. An example of such comparison is shown in Fig. 12 for one of the studied subjects. Similarly with comparison of BPP calculations, the ROIs size was 10 × 10 pixels.
It is seen that both sinusoidal and triangular reference functions of both types provide reliable estimation of the blood pulsation amplitude in ROI-2 [Fig. 12(a)]. Note that this ROI was selected in the region in which the blood pulsates synchronously with the reference area. Such types of ROI are dominating in all studied subjects (see Sect. 3.2). In contrast, we observed that the D-type triangular function results in underestimated BPA in the ROI-3 characterized by asynchronous blood pulsations [Fig. 12(b)]. This underestimation of the pulsations amplitude by the BPI algorithm with function was found in all five subjects in areas with asynchronous blood supply. Nevertheless, sinusoidal reference functions provide reliable BPA calculations even in the regions with asynchronous pulsations.
In this study by using an advance mode of PPG we found new temporal readouts for the variable peripheral microcirculation. The amplitude of the pulsatile component of the PPG signal is known to be influenced by respiration, sympathetic nervous system activity as well as other factors [2,4,19–21]. The shape of the recorded signal in the traditional modes of PPG, however, remains approximately constant [22–24]. It is characterized by a slow transition time from systole to diastole [9,10,17–19]. Using an advance mode of PPG called BPI, we observed for the first time that there are limited regions in healthy adult subjects in which the PPG fast type waveform shape becomes inverted with the faster transition from systole to diastole. This fast type of the PPG waveform co-exists with the prevailing traditional slow type in all studied subjects. Such heterogeneity of the PPG waveforms leads to ambiguity in the measurements of spatial distribution of the blood pulsation phase: in a pair of spatially separated regions, the time intervals between systole phases can be different from those of the diastole phase.
A possible explanation for the co-existence of slow and fast types of waveforms is the dominating presence of the vascular tone provided by the sympathetic nerve activity  leading to the slow growth of PPGs (slow-type waveform) in most of the studied regions. In contrast, in a limited number of atypical regions, this vascular tone was attenuated resulting in fast-type PPGs. Even though detectable in fewer regions of the palm and wrist, such fast-type signals can potentially serve as novel biomarkers of the microcirculation, as one can predict that the presence of regions with fast-type signals will be altered during pathological states associated with attenuated skin sympathetic nerve activity induced by emotional stress , exposure to low or high temperatures , or exercise training . It also remains to be established if there are other body regions where this fast type waveforms are more abundant.
For a complete 2D description of the complex processes of microcirculation we are proposing a new algorithm of Blood Pulsation Imaging (BPI), which includes the proper choice of the reference area and a synchronous detection of the blood pulsations with two types of reference functions, defined for both the systole and diastole phases. We have compared the results of BPP mapping using the BPI algorithm with sinusoidal and triangular reference functions of D- and S- types, and found that there is not a single reference function for the reliable estimation of the pulsation phase in all regions. Statistically, D-type sinusoidal reference function provides a more reliable BPP map for the diastole phase, while the S-type triangular reference function provides more reliable calculations for the systole phase. Both these functions can also be used for the calculations of blood pulsation amplitude, BPA mapping.
It is worth noting that the atypical fast PPG waveforms were observed in regions to which the blood is supplied asynchronously compared with the dominating regions with typical, slow type of PPG signal. In addition to allowing the observation of the microcirculation heterogeneity, using the analysis of combined BPA and BPP maps allowed us to detect ‘hot spots’  and regions of asynchronous blood supply [12,13] but also to characterize the non-uniform shape of waveforms in the neighboring skin regions.
VT acknowledges the Department of Applied Physics, University of Eastern Finland, Kuopio for the partial financial support. Research of RG was supported by the Program of Competitive Growth of the Kazan Federal University. The authors are grateful to Genevieve Bart for the reading of the manuscript.
References and links
1. A. B. Hertzman and C. R. Spealman, “Observations on the finger volume pulse recorded photoelectrically,” Am. J. Physiol. 119, 334–335 (1937).
4. K. Nakajima, T. Tamura, and H. Miike, “Monitoring of heart and respiratory rates by photoplethysmography using a digital filtering technique,” Med. Eng. Phys. 18(5), 365–372 (1996). [CrossRef] [PubMed]
7. N. Selvaraj, A. K. Jaryal, J. Santhosh, K. K. Deepak, and S. S. Anand, “Assessment of heart rate variability derived from finger-tip photoplethysmography as compared to electrocardiography,” J. Med. Eng. Technol. 32(6), 479–484 (2008). [CrossRef] [PubMed]
8. F. P. Wieringa, F. Mastik, and A. F. W. van der Steen, “Contactless multiple wavelength photoplethysmographic imaging: a first step toward “SpO2 camera” technology,” Ann. Biomed. Eng. 33(8), 1034–1041 (2005). [CrossRef] [PubMed]
9. K. Humphreys, T. Ward, and C. Markham, “Noncontact simultaneous dual wavelength photoplethysmography: a further step toward noncontact pulse oximetry,” Rev. Sci. Instrum. 78(4), 044304 (2007). [CrossRef] [PubMed]
11. A. A. Kamshilin, S. V. Miridonov, V. Teplov, R. Saarenheimo, and E. Nippolainen, “Photoplethysmographic imaging of high spatial resolution,” Biomed. Opt. Express 2(4), 996–1006 (2011). [CrossRef] [PubMed]
12. A. A. Kamshilin, V. Teplov, E. Nippolainen, S. V. Miridonov, and R. Giniatullin, “Variability of microcirculation detected by blood pulsation imaging,” PLoS ONE 8(2), e57117 (2013). [CrossRef] [PubMed]
13. N. Zaproudina, V. Teplov, E. Nippolainen, J. A. Lipponen, A. A. Kamshilin, M. Närhi, P. A. Karjalainen, and R. Giniatullin, “Asynchronicity of facial blood perfusion in migraine,” PLoS ONE 8(12), e80189 (2013). [CrossRef] [PubMed]
14. V. Teplov, A. Shatillo, E. Nippolainen, O. Gröhn, R. Giniatullin, and A. A. Kamshilin, “Fast vascular component of cortical spreading depression revealed in rats by blood pulsation imaging,” J. Biomed. Opt. 19(4), 046011 (2014). [CrossRef] [PubMed]
16. G. G. Berntson, J. T. Bigger Jr, D. L. Eckberg, P. Grossman, P. G. Kaufmann, M. Malik, H. N. Nagaraja, S. W. Porges, J. P. Saul, P. H. Stone, and M. W. van der Molen, “Heart rate variability: Origins, methods, and interpretive caveats,” Psychophysiology 34(6), 623–648 (1997). [CrossRef] [PubMed]
17. B. Khanokh, Y. Slovik, D. Landau, and M. Nitzan, “Sympathetically induced spontaneous fluctuations of the photoplethysmographic signal,” Med. Biol. Eng. Comput. 42(1), 80–85 (2004). [CrossRef] [PubMed]
19. C. G. Scully, J. Lee, J. Meyer, A. M. Gorbach, D. Granquist-Fraser, Y. Mendelson, and K. H. Chon, “Physiological parameter monitoring from optical recordings with a mobile phone,” IEEE Trans. Biomed. Eng. 59(2), 303–306 (2012). [CrossRef] [PubMed]
20. A. Johansson and P. Å. Öberg, “Estimation of respiratory volumes from the photoplethysmographic signal. Part I: experimental results,” Med. Biol. Eng. Comput. 37(1), 42–47 (1999). [CrossRef] [PubMed]
21. M. Nitzan, A. Babchenko, B. Khanokh, and D. Landau, “The variability of the photoplethysmographic signal - a potential method for the evaluation of the autonomic nervous system,” Physiol. Meas. 19(1), 93–102 (1998). [CrossRef] [PubMed]
22. A. B. Hertzman, “The blood supply of various skin areas as estimated by the photoelectric plethysmograph,” Am. J. Physiol. 124, 328–340 (1938).
23. H. Lax, A. W. Feinberg, and B. M. Cohen, “Studies of the arterial pulse wave. I. The normal pulse wave and its modification in the presence of human arteriosclerosis,” J. Chronic Dis. 3(6), 618–631 (1956). [CrossRef] [PubMed]
24. S. C. Millasseau, J. M. Ritter, K. Takazawa, and P. J. Chowienczyk, “Contour analysis of the photoplethysmographic pulse measured at the finger,” J. Hypertens. 24(8), 1449–1456 (2006). [CrossRef] [PubMed]
25. V. G. Macefield, “Sympathetic microneurography,” in Autonomic Nervous System, R. M. Buijs and D. F. Swaab, eds., (Elsevier B.V., 2013), pp. 353–364.
26. R. Brown, C. James, L. A. Henderson, and V. G. Macefield, “Autonomic markers of emotional processing: skin sympathetic nerve activity in humans during exposure to emotionally charged images,” Front Physiol 3, 394 (2012). [CrossRef] [PubMed]