Imaging system and method for contactless detection of pulse wave spatial and temporal distribution and characteristics

Abstract

An imaging system and method for contactless detection of the spatial and temporal distribution and characteristics of pulse waves combines an illumination light source capable of generating an incident light beam with a peak wavelength and wavelength bandwidth, an imaging device capable of measuring light scattered from human tissue, a polarization element capable of controlling the polarization states of the incident light beam and the light scattered from human tissue, a camera capable of capturing video at a high frame rate, and an analysis algorithm capable of processing iPPG video data composed of multiple images. This system obtains a human pulse wave signal with a high signal-to-noise ratio, as well as the spatial distribution of the pulse wave signal and the average pulse wave temporal waveform data at different pixel locations. The system first sets the video camera's image data acquisition parameters, then activates the camera to capture and store multiple iPPG video data frames. The system then calculates a human pulse wave signal with a high signal-to-noise ratio, as well as the spatial distribution of the pulse wave signal and the average pulse wave temporal waveform data at different spatial locations. This system can quantify the intensity of arterial wall sclerosis and thus quantitatively estimate the risk of cardiovascular and metabolic diseases.

Description

Technical Field

[0001] The present invention relates to a contactless imaging system and method, and more particularly to an imaging system and method for contactless detection of pulse wave spatial and temporal distribution and characteristics. Background Art

[0002] The beating of the human heart generates pressure changes in the blood, which propagate as waves within the blood vessels and cause changes in vascular volume, forming a pulse, referred to here as a pulse wave. Pulse waves drive blood flow within the circulatory system and are essential for maintaining life; therefore, pulse wave detection can provide key vital sign data. Furthermore, quantitative detection of the spatial and temporal distribution of pulse waves can provide big data on the status of the heart and circulatory system, useful for testing, monitoring, and caring for patients with metabolic and cardiovascular diseases. Traditional arterial imaging methods, including magnetic resonance imaging, X-rays, and computed tomography (CT), can image blood vessels in human tissues. These methods offer high spatial resolution and penetration depth, but these methods suffer from low temporal resolution, require the injection of contrast agents, are radioactive, and are expensive, making them unsuitable for large-scale population testing or patient treatment monitoring. Ultrasound arterial imaging is a safe imaging method, but its spatial and temporal resolution is limited. In contrast, label-free optical imaging offers the advantages of being radioactive and having high spatial and temporal resolution, enabling the detection of blood flow parameters within superficial tissues. For example, speckle and Doppler imaging, both methods using coherent light illumination, can measure the velocity of particles in the blood, such as red blood cells, and are already in clinical use. However, none of these methods can provide vital sign data related to pulse waves and heartbeats in arteries.

[0003] Photoplethysmography (PPG) is a photoelectric measurement technique that measures scattered light signals emitted from human tissue. After signal processing, it can be used to detect pulse waves propagating in blood within arteries. The PPG signal can be used to measure heart rate and other blood parameters, such as oxygen saturation. Pulse oximeters developed based on PPG technology have been widely used in clinical and home settings. These instruments utilize a dual-wavelength light-emitting diode and photodiode sensor configuration to accurately measure blood oxygen saturation and heart rate. The PPG principle has attracted the interest of many researchers due to its unique ability to capture pulse wave and blood parameters and its simplicity. In recent years, its application in wearable devices has received considerable attention. Although PPG technology has many advantages, its requirement for contact with human tissue can lead to the spread of pathogens among the subjects being tested, making it unsuitable for clinical applications such as wound monitoring. The spatial localization of PPG signals further limits its potential for detecting pulse wave propagation characteristics associated with the degree of arterial wall stiffness. Therefore, people hope to effectively utilize the convenience of PPG technology in measuring pulse waves and blood parameters and develop new methods that can measure the spatial and temporal distribution of pulse waves in a contactless manner.

[0004] In recent years, imaging photoplethysmography (iPPG) has emerged as an extension of PPG technology. By recording high-speed video signals, it measures the distribution of pulse waves within a region of interest and the time-domain waveform of pulse waves at different locations without contacting human tissue. To date, most iPPG research reports have selected one or more regions of interest within the field of view and then calculated the average pixel value within that region to improve the signal-to-noise ratio of very weak pulse wave signals. However, this average calculation of all pixels results in a complete or significant loss of spatial resolution of the pulse wave signal. To improve spatial resolution, researchers have reported using temporal correlation between iPPG video image data and a pulse wave reference signal to remove non-pulse wave noise from iPPG image data. However, the pulse wave reference signal requires measurement using other methods, such as PPG technology or electrocardiogram (ECG) signals, which requires sensors in contact with skin tissue and increases the difficulty of signal measurement. iPPG methods reported in the literature also include extracting pulse wave signals by calculating the average pixel value of all pixels within a region of interest (ROI) across different color channels within iPPG color video data as input, performing independent component analysis on this value. Another method involves analyzing the temporal correlation of image pixels within the field of view to determine the spatial distribution of the pulse wave. These methods select a region of interest (ROI) with a strong pulse wave signal and calculate the average pixel value of all pixels within that region. This value is then used as a pulse wave reference signal. The temporal correlation between each pixel within the ROI and the reference signal is then calculated, and the pulse wave signal component of this value, which varies over time, is extracted. Although these methods can extract pulse wave signals from iPPG video data, they generally still require calculating the average pixel value of all pixels within the ROI, or the spatial distribution signal-to-noise ratio of the extracted pulse wave signal remains low. Consequently, significant improvement in the spatial and temporal resolution of the pulse wave signal is difficult to achieve, hindering widespread application. Summary of the Invention

[0005] The technical problem to be solved by this invention is to overcome the shortcomings of existing technologies by providing an imaging system and method for contactless detection of the spatial and temporal distribution and characteristics of pulse waves. This system can extract human pulse wave signal components with high signal-to-noise ratio and high spatial resolution from imaging photoplethysmography video data, as well as the spatial distribution of this signal and pulse wave temporal waveform data at different spatial locations.

[0006] The technical solution adopted by the present invention is: an imaging system for contactlessly detecting the spatial and temporal distribution and characteristics of pulse waves, comprising a computer for storing, displaying, and processing image data, an imaging device connected to the computer via a data transmission line for acquiring M images of a human tissue region, and a light source assembly for providing incident light. The imaging device comprises a video camera, an imaging lens, and an imaging polarizer arranged sequentially in the same direction for receiving reflected light from the human tissue region. The imaging polarizer restricts the imaging lens to only receive reflected light with a linear polarization state and a polarization direction perpendicular to the polarization direction of the incident light. The M images acquired by the imaging device constitute video images from which iPPG video data of the human tissue region is extracted. The light source assembly comprises 1 to 10 light source groups arranged at equal angles and intervals around the imaging device for illuminating the human tissue region. Each light source group comprises a light source and a light source polarizer arranged coaxially with the light source. The light illuminating the human tissue region is linearly polarized incident light by adjusting the angle of rotation of the light source polarizer about the coaxial axis.

[0007] An image processing algorithm for calculating human pulse frequency from iPPG video data based on an imaging system that detects the spatial and temporal distribution and characteristics of pulse waves without contact, comprising the following steps:

[0008] 1) Determine whether each image in the iPPG video data collected from the M images contains saturated pixels. If each image does not contain saturated pixels, proceed to step 2); if more than one image contains saturated pixels, select the image with the most saturated pixels and calculate the number of saturated pixels p of this image. s The total number of pixels p of the image t and determine whether the ratio is greater than 2%. If so, the image processing process ends, otherwise, proceed to step 2);

[0009] 2) Based on the pixel value histogram of each image in the iPPG video data, calculate the boundary threshold of the region of interest with human blood vessel and tissue characteristics in the corresponding image. Based on the region of interest of each image, find the common region of interest of all images and intercept the corresponding common region of interest in each image;

[0010] 3) According to the common region of interest of each image, the average pixel value of the common region of interest of each image in the iPPG video data at the acquisition time t is calculated, which is recorded as I av (t;λ p ), where 0≤t≤T, T is the video acquisition time, λ p is the peak wavelength of the incident light, and calculates the average pixel value I av (t;λ p ) time average μav (λ p ) and standard deviation σ av (λ p ), calculate the average pixel standard score ΔI of the common area of ​​interest in each image according to the following formula av (t;λ p ),Right now

[0011]

[0012] 4) Use Fourier transform algorithm to calculate the average pixel standard score value ΔI av (t;λ p ) of the frequency domain signal Δi av (f;λ p ), determine the pulse frequency f of a given human body h The pulse frequency measurement value f of the strongest component within the range c , the pulse frequency measurement value f c Refers to the given human pulse frequency f h The frequency with the largest Fourier component value within the range, the human pulse frequency f h The range is 0.6 Hz to 4.0 Hz.

[0013] An image processing algorithm based on calculating human pulse frequency, and a pixel independent component analysis method for calculating the spatial distribution image and time domain distribution data of human pulse wave signals, comprising the following steps:

[0014] 1) Based on the standard deviation normalization method, the single-pixel standard score value I of the common area of ​​interest in each image in the iPPG video data is calculated according to the following formula: n (x,y;t;λ p ),Right now

[0015]

[0016] Where x and y are the coordinates of a single pixel position; I(x, y; t; λ p ) is the single pixel value with coordinates (x, y); μ I (x,y;λ p ) and σ I (x,y;λ p ) are the temporal average and standard deviation of the single pixel value with coordinates (x, y); λ p is the peak wavelength of the incident light;

[0017] 2) Select a single pixel with coordinates (x, y), divide the surrounding pixels centered on the pixel into K-1 surrounding pixel groups, and calculate the average pixel standard score value I of each surrounding pixel group nz(x, y; t; λ0), using the average pixel standard score value I of the surrounding pixel group nz The values ​​of (x, y; t; λ0) at different times are used as components, expressed as the average pixel standard score value time vector I nz (x,y;λ p ), where z = 2, ..., K, where the value of K is the average pixel standard score value I affecting each pixel group nz The number of signal sources (x, y; t; λ0) that change with time is the same;

[0018] 3) Calculate the signal source S of each single pixel with coordinates (x, y) k (x,y;t;λ p )’s signal source time vector S k (x,y;λ p ), k=1,…,K, wherein the signal source time vector S k (x,y;λ p ) must satisfy the condition of maximum temporal statistical independence between the time vectors of each signal source and is defined as the output time vector of the independent component analysis algorithm of the pixel group, called the independent component time vector;

[0019] 4) Using Fourier transform algorithm, calculate each independent component time vector S k (x,y;λ p ) or the corresponding signal source S k (x,y;t;λ p ) spectrum distribution, determine the signal source S k (x,y;t;λ p ) The maximum component frequency value f0 within the human pulse frequency range, the K independent component time vectors have the closest pulse frequency measurement value f c The independent component time vector of the f0 value is the pulse wave independent component time vector, denoted as S c (x,y;λ p ), the pulse wave independent component time domain signal is S c (x,y;t;λ p ), is the pulse wave independent component time domain signal of the pixel value of a single pixel with coordinates (x, y);

[0020] 5) Using filtering algorithm, calculate the pulse wave independent component time domain signal S c (x,y;t;λ p ) The filtered signal is recorded as the pulse wave time domain signal S cf (x,y;t;λ p ) or pulse wave time vector S cf (x,y;λ p );

[0021] The filtering algorithm is based on Fourier transform in the frequency domain, or adopts a digital filtering method such as a finite impulse response filtering algorithm;

[0022] 6) Based on Fourier transform, calculate the pulse wave time domain signal S of each single pixel with coordinates (x, y) cf (x,y;t;λ p ) is the frequency distribution signal, recorded as the pulse wave frequency domain signal s cf (x,y;f;λ p );

[0023] 7) According to the pulse wave frequency domain signal s cf (x,y;f;λ p ) and the following formula to calculate the pulse wave index r of the single pixel value with coordinates (x, y) one by one cf (x,y;λ p ),Right now

[0024]

[0025] Where Δf = F / 2d is the frequency domain step size, F is the frame rate of the video camera, d is the smallest integer that satisfies the condition M≤2d≤2M, M is the number of images, N≥2, W≥N are the number of terms in the sum of the numerator and denominator in formula (4), respectively;

[0026] 8) Repeat steps 4) to 7) until the pulse wave time domain signal S of all pixels in the common region of interest in each image is completed. cf (x,y;t;λ p ) and the pulse wave index r of all pixels cf (x,y;λ p ) calculation;

[0027] 9) The pulse wave index r of all pixels in the common region of interest in each image obtained in step 8) cf (x,y;λ p ), calculate the maximum and minimum values ​​of the pulse wave index in the common area of ​​interest, and then calculate the pulse wave index r based on the maximum and minimum values cf (x,y;λ p ) normalized value, according to the pulse wave index r cf (x,y;λ p )'s (x, y) coordinates and normalized values, generate and output each image in the common region of interest and with the incident light peak wavelength λ p Pulse wave signal spatial distribution image data;

[0028] 10) According to the pulse wave signal spatial distribution image data in the common region of interest of each image obtained in step 9), select the image with coordinates (x, y) in the common region of interest of each image and r cf (x,y;λ p ) is greater than the set value, forming a high signal-to-noise ratio pixel set. The single pixel with coordinates (x, y) and the surrounding pixels centered on the single pixel are divided into K pixel groups. The number of pixels in each pixel group is greater than 1. The average pixel value I of each pixel group at the acquisition time t is calculated. ka (x,y;t;λ p ), k=1,…,K;

[0029] 11) The average pixel value I of the K pixel groups of high signal-to-noise ratio pixels located at (x, y) obtained in step 10) ka (x,y;t;λ p ), calculate the average pixel standard score value I of each pixel group according to the following formula nka (x,y;t;λ p ):

[0030]

[0031] Among them, μ Ika (x,y;λ p ) and σ Ika (x,y;λ p ) are the average pixel values ​​I ka (x,y;t;λ p )'s time average and standard deviation values; the average pixel standard score value I nka (x,y;t;λ p ) corresponds to the time vector I nka (x,y;λ p ) is used as the input time vector in formula (3), and the average signal source, that is, the average independent component time vector, is calculated according to steps 12) to 15), which is recorded as S ka (x,y;λ p ), in the average independent component time vector S ka (x,y;t;λ p ) to determine the average pulse wave independent component signal S ca (x,y;t;λ p ), using filtering algorithm, calculate the average pulse wave time domain signal S of the pixel at (x, y) caf (x,y;t;λ p ), calculate and compare the time domain waveform phases of pixels at different positions in the high signal-to-noise ratio pixel set, and exclude high signal-to-noise ratio pixels whose phase changes exceed 90 degrees compared with adjacent high signal-to-noise ratio pixels;

[0032] 12) The average pulse wave time domain signal S of each pixel in the high signal-to-noise ratio pixel set obtained according to steps 10) and 11) caf (x,y;t;λ p ), calculating the propagation time and propagation time delay of two pixels located at different coordinate positions with the same time domain waveform peak phase, calculating the spatial distance between the two coordinate positions, and calculating the pulse wave propagation velocity at the corresponding waveform peak between the two positions based on the ratio of the spatial distance to the time delay and the arterial blood flow and pressure wave propagation direction;

[0033] 13) For each combination of two pixels in the high signal-to-noise ratio pixel set obtained in steps 10) and 11), calculate the pulse wave propagation velocity of all peaks of each combination within the video acquisition time T according to step 12), thereby obtaining the spatial distribution, temporal distribution, mean values ​​of the spatial distribution and temporal distribution, standard deviations of the spatial distribution and temporal distribution, and related statistical characteristic parameters of the pulse wave propagation velocity within the common region of interest.

[0034] The imaging system and method for contactlessly detecting the spatial and temporal distribution and characteristics of pulse waves disclosed herein are an imaging detection method for contactlessly capturing iPPG video data. Combining an illumination light source capable of generating incident light with a peak wavelength and wavelength bandwidth, a polarization element and imaging device capable of controlling the polarization state of the incident light and scattered light from human tissue, and a video camera with a frame rate of F, the system captures iPPG video data comprising multiple images, extracts a pulse wave spatial distribution map and temporal signal data, obtains and outputs a pulse wave signal spatial distribution map, selects multiple locations with a high signal-to-noise ratio, calculates the pixel group time vector and average pulse wave temporal signal waveform data for each selected location, thereby significantly improving the spatial and temporal resolution of the pulse wave signal. Furthermore, the system accurately calculates and compares the average pulse wave temporal signal waveform phase at different pixel locations in the iPPG video data, and calculates the pulse wave propagation velocity and related statistical characteristic parameters between pixels with the same phase but different locations. These parameters can be used to quantify the intensity of arterial wall sclerosis, thereby quantitatively estimating the risk of cardiovascular and metabolic diseases. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 Schematic diagram of the imaging system for contactless detection of pulse wave spatial and temporal distribution and characteristics of the present invention;

[0036] Figure 2 is a flow chart of an image processing method for calculating human pulse frequency based on iPPG video data according to the present invention;

[0037] Figure 3 It is the distribution diagram of the common region of interest and the average pixel value of the region and the normalized value of the average pixel value in the time domain and frequency domain;

[0038] Figure 4 is a flowchart of the pixel group independent component analysis method of the present invention);

[0039] Figure 5 Taking K=3 as an example, the image processing calculation flow based on the pixel group independent component analysis algorithm is shown;

[0040] Figure 6 Taking K=3 as an example, a schematic diagram of a group of three pixels with the center position (x, y) is selected;

[0041] Figure 7 are the 10 different incident light peak wavelengths λ obtained after imaging the left palm in Example 1 p Pulse wave spatial distribution diagram;

[0042] Figure 8 are the 10 different incident light peak wavelengths λ obtained after imaging the left palm of Example 2 p Pulse wave spatial distribution diagram;

[0043] Figure 9 It is the high signal-to-noise ratio pixel-averaged pulse wave time domain signal at four different locations when the incident light peak wavelength is 530 nm and 850 nm;

[0044] Figure 10 It is the average pulse wave time domain signal of high signal-to-noise ratio pixels located at 12 different positions when the peak wavelength of the incident light is 850 nanometers;

[0045] Figure 11 yes Figure 10 The delay time of the average pulse wave time domain signal of high signal-to-noise ratio pixels at three different positions is shown. DETAILED DESCRIPTION

[0046] The imaging system and method for contactless detection of pulse wave spatial and temporal distribution and characteristics are described in detail below with reference to the embodiments and drawings.

[0047] like Figure 1As shown, the imaging system for contactlessly detecting the spatial and temporal distribution and characteristics of pulse waves of the present invention includes a computer 8 for storing, displaying, and processing image data, an imaging device connected to the computer 8 via a data transmission line 7 for acquiring M images of a human tissue region 3, and a light source assembly for providing incident light. The imaging device includes a video camera 6, an imaging lens 5, and an imaging polarizer 4 arranged sequentially in the same direction for receiving reflected light from the human tissue region 3. By appropriately setting the optical axis direction of the imaging polarizer 4, the imaging lens 5 can be restricted to only receive diffusely reflected light having a linear polarization state and a polarization direction perpendicular to the polarization direction of the incident light. The M images acquired by the imaging device constitute video images from which iPPG video data of the human tissue region 3 is extracted. The light source assembly includes 1 to 10 light source groups arranged at equal angles and intervals around the imaging device for illuminating the human tissue region 3. Each light source group includes a light source 1 and a light source polarizer 2 arranged coaxially with the light source 1. The light illuminating the human tissue region 3 is incident light having a linear polarization state by adjusting the angle of rotation of the light source polarizer 2 about the coaxial axis.

[0048] The frame rate F of the video camera 6 is greater than the human pulse frequency f h , 10HZ ~ 1000HZ; pixels are 250,000 ~ 10 million; video acquisition time T value is greater than 1 / f h , i.e., greater than 2 seconds. A video camera model FS-1600D-10GE produced by JAI can be used, with a frame rate F of up to 250 Hz, to collect iPPG video data having M images.

[0049] The diameters of the imaging polarizer 4 and the light source polarizer 2 range from 5 mm to 200 mm. For example, the light source polarizer 2 can be a Thorlabs-branded LPNIRE-B near-infrared polarizer, which forms a linearly polarized incident light beam that is incident on the human tissue 3 to be examined. The imaging polarizer 4 can also be a Thorlabs-branded LPNIRE-B near-infrared polarizer.

[0050] The light source 1 is one of an incandescent lamp, a light emitting diode, a fluorescent lamp, and a laser, with a wavelength range of 400nm to 1100nm and a bandwidth of 0.01nm to 1000nm. The angle between the central axis of the incident light generated by each light source group and the vertical direction is 5° to 85°. For example, a plurality of high-power, peak wavelength λ p It is a ring light source composed of 850 nanometer light-emitting diodes.

[0051] The imaging lens 5 can be selected from existing camera lenses with different focal lengths according to the required field of view, such as the camera lens model #58-000 with a focal length of 8 mm sold by Edmund, or a camera lens made of other optical elements.

[0052] The human tissue area 3 is, for example, palm skin tissue, face or ear.

[0053] like Figure 2 As shown, the image processing algorithm of the present invention for calculating human pulse frequency from iPPG video data based on an imaging system for contactless detection of pulse wave spatial and temporal distribution and characteristics includes the following steps:

[0054] 1) Determine whether each image in the iPPG video data collected from the M images contains saturated pixels. If each image does not contain saturated pixels, proceed to step 2); if more than one image contains saturated pixels, select the image with the most saturated pixels and calculate the number of saturated pixels p of this image. s The total number of pixels p of the image t and determine whether the ratio is greater than 2%. If so, the image processing process ends, otherwise, proceed to step 2);

[0055] 2) Based on the pixel value histogram of each image in the iPPG video data, calculate the boundary threshold of the region of interest with human blood vessel and tissue characteristics in the corresponding image. Based on the region of interest of each image, find the common region of interest of all images and intercept the corresponding common region of interest in each image;

[0056] 3) According to the common region of interest of each image, the average pixel value of the common region of interest of each image in the iPPG video data at the acquisition time t is calculated, which is recorded as I av (t;λ p ), where 0≤t≤T, T is the video acquisition time, λ p is the peak wavelength of the incident light, and calculates the average pixel value I av (t;λ p ) time average μ av (λ p ) and standard deviation σ av (λ p ), calculate the average pixel standard score ΔI of the common area of ​​interest in each image according to the following formula av (t;λ p ),Right now

[0057]

[0058] 4) Use Fourier transform algorithm to calculate the average pixel standard score value ΔI av(t;λ p ) of the frequency domain signal Δi av (f;λ p ), determine the pulse frequency f of a given human body h The pulse frequency measurement value f of the strongest component within the range c , the pulse frequency measurement value f c Refers to the given human pulse frequency f h The frequency with the largest Fourier component value within the range, the human pulse frequency f h The range is 0.6 Hz to 4.0 Hz.

[0059] The present invention can use a variety of methods to calculate the common region of interest in the iPPG video data. One implementation method can be based on a histogram data analysis algorithm of pixel values ​​(see N. Otsu, "AThreshold Selection Method from Gray-Level Histograms," IEEE Trans. Syst. Man Cybern. 9, 62-66 (1979)) to obtain a pixel value threshold that can determine the boundary of the common region of interest, such as Figure 3 As shown in (A), as another implementation method, the pixel value threshold that can determine the boundary of the common region of interest can also be obtained based on the fuzzy clustering algorithm (see D.Graves and W.Pedrycz, "Kernel-based fuzzy clustering and fuzzyclustering: A comparative experimental study," Fuzzy Sets and Systems, 161(4), 522-543(2010)).

[0060] Figure 3 is the distribution diagram of the common region of interest and the average pixel value and the normalized value of the average pixel value in the time domain and frequency domain;

[0061] Figure 3 (A) Reads all images in the iPPG video data and then calculates the common region of interest (ROI) within all M images using a boundary threshold calculation method. The ROI is the human tissue area within the field of view, such as palm skin tissue, marked as the white area in the black and white image.

[0062] Figure 3 (B) is the average pixel value of the image in the common area of ​​interest I av (t;λ p ) and the image acquisition time t;

[0063] Figure 3(C) is the average pixel standard score value ΔI of the common area of ​​interest av (t;λ p ) and the image acquisition time t;

[0064] Figure 3 (D) is ΔI av (t;λ p ) of the spectrum distribution signal Δi av (f;λ p An example of the functional relationship between ) and signal frequency f, f c It is the human pulse frequency value determined within a given human pulse frequency range.

[0065] Figure 3 (B) to Figure 3 (D) Comparison of the average pixel value of each image I av (t;λ p ), the average pixel standard score value ΔI of each image calculated according to formula (1) av (t;λ p ) and the frequency domain distribution signal Δi calculated according to Fourier transform av (f;λ p ) data, from which we can see that the standard score calculation defined by formula (1) can significantly enhance the frequency domain distribution signal Δi av (f;λ p ) The frequency peak value within the human heart rate range (e.g., between 0.6 Hz and 4.0 Hz) corresponds to the human heart rate measurement value f c ,exist Figure 3 (D) shows 1.40 Hz or 84 beats per minute. In the subsequent independent component analysis calculation of the pixel group, the measured value f c It can be used to determine the pulse wave independent component vector from the output K independent component time vectors.

[0066] like Figure 4 As shown, the pixel group independent component analysis method for calculating the spatial distribution image and time domain distribution data of the human pulse wave signal of the present invention includes the following steps:

[0067] 1) Based on the standard deviation normalization method, the single-pixel standard score value I of the common area of ​​interest in each image in the iPPG video data is calculated according to the following formula: n (x,y;t;λ p ),Right now

[0068]

[0069] Where x and y are the coordinates of a single pixel position; I(x, y; t; λ p) is the single pixel value with coordinates (x, y); μ I (x,y;λ p ) and σ I (x,y;λ p ) are the temporal average and standard deviation of the single pixel value with coordinates (x, y); λ p is the peak wavelength of the incident light;

[0070] 2) Select a single pixel with coordinates (x, y), divide the surrounding pixels centered on the pixel into K-1 surrounding pixel groups, and calculate the average pixel standard score value I of each surrounding pixel group nz (x, y; t; λ0), using the average pixel standard score value I of the surrounding pixel group nz The values ​​of (x, y; t; λ0) at different times are used as components, expressed as the average pixel standard score value time vector I nz (x,y;λ p ), where z = 2, ..., K, where the value of K is the average pixel standard score value I affecting each pixel group nz The number of signal sources (x, y; t; λ0) that change with time is the same;

[0071] The average pixel standard score value I of each pixel group nz (x,y;t;λ p ) is defined according to the distance from a single pixel with coordinates (x, y), wherein the peripheral pixel group with z=2 is composed of all pixels in a square formed with a set length as the side length and centered on the single pixel with coordinates (x, y) excluding the single pixel with coordinates (x, y); the peripheral pixel group with z=3 is composed of all pixels in a square ring with a set width formed along the outer periphery of the peripheral pixel group with z=2 and centered on the single pixel with coordinates (x, y); .....and the peripheral pixel group with z=K is composed of all pixels in a square ring with a set width formed along the outer periphery of the peripheral pixel group with z=K-1 and centered on the single pixel with coordinates (x, y).

[0072] 3) Calculate the signal source S of each single pixel with coordinates (x, y) k (x,y;t;λ p )’s signal source time vector S k (x,y;λ p ), k=1,…,K, wherein the signal source time vector S k (x,y;λ p ) must satisfy the condition of maximum temporal statistical independence between the time vectors of each signal source and is defined as the output time vector of the independent component analysis algorithm of the pixel group, called the independent component time vector;

[0073] The calculation coordinates of each signal source S in a single pixel (x, y) k (x,y;t;λ p )’s signal source time vector S k (x,y;λ p ), firstly, the standard score value of a single pixel with coordinates (x, y) I n (x,y;t;λ p ) is represented as a single pixel standard score value time vector I n1 (x,y;λ p ), and the average pixel standard score value time vector I of the K-1 surrounding pixel groups nz (x,y;λ p ) are combined into K time vectors I nk (x,y;λ p ), where k = 1, ..., K, set I nk (x,y;λ p ) and the time vector S of the signal source that affects the time variation k (x,y;λ p ) is a linear relationship, which can be expressed by the following matrix formula:

[0074]

[0075] Where [A] is an unknown matrix, also called a mixing matrix; the K time vectors I are transformed into nk (x,y;λ p ) is used as the input time vector of the pixel group independent component analysis algorithm, and the matrix [A] and the signal source time vector S that satisfy formula (3) are iteratively solved. k (x,y;λ p ).

[0076] 4) Using Fourier transform algorithm, calculate each independent component time vector S k (x,y;λ p ) or the corresponding signal source S k (x,y;t;λ p ) spectrum distribution, determine the signal source S k (x,y;t;λ p ) The maximum component frequency value f0 within the human pulse frequency range, the independent component time vector with the f0 value closest to the pulse frequency measurement value fc among the K independent component time vectors is the pulse wave independent component time vector, denoted by S c (x,y;λ p ), the pulse wave independent component time domain signal is S c(x, y; t; λp) is the pulse wave independent component time domain signal of the pixel value of a single pixel with coordinates (x, y);

[0077] 5) Using filtering algorithm, calculate the pulse wave independent component time domain signal S c (x,y;t;λ p ) The filtered signal is recorded as the pulse wave time domain signal S cf (x,y;t;λ p ) or pulse wave time vector S cf (x,y;λ p );

[0078] Wherein, the filtering algorithm is based on Fourier transform in the frequency domain, or adopts a digital filtering method such as a finite impulse response filtering algorithm;

[0079] 6) Based on Fourier transform, calculate the pulse wave time domain signal S of each single pixel with coordinates (x, y) cf (x,y;t;λ p ) is the frequency distribution signal, recorded as the pulse wave frequency domain signal s cf (x,y;f;λ p );

[0080] 7) According to the pulse wave frequency domain signal s cf (x,y;f;λ p ) and the following formula to calculate the pulse wave index r of the single pixel value with coordinates (x, y) one by one cf (x,y;λ p ),Right now

[0081]

[0082] Where Δf = F / 2d is the frequency domain step size, F is the frame rate of the video camera, d is the smallest integer that satisfies the condition M≤2d≤2M, M is the number of images, N≥2, W≥N are the number of terms in the sum of the numerator and denominator in formula (4), respectively;

[0083] 8) Repeat steps 4) to 7) until the pulse wave time domain signal S of all pixels in the common region of interest in each image is completed. cf (x,y;t;λ p ) and the pulse wave index r of all pixels cf (x,y;λ p ) calculation;

[0084] 9) The pulse wave index r of all pixels in the common region of interest in each image obtained in step 8) cf (x,y;λ p), calculate the maximum and minimum values ​​of the pulse wave index in the common area of ​​interest, and then calculate the pulse wave index r based on the maximum and minimum values cf (x,y;λ p ) normalized value, according to the pulse wave index r cf (x,y;λ p )'s (x, y) coordinates and normalized values, generate and output each image in the common region of interest and with the incident light peak wavelength λ p Pulse wave signal spatial distribution image data;

[0085] 10) According to the pulse wave signal spatial distribution image data in the common region of interest of each image obtained in step 9), select the image with coordinates (x, y) in the common region of interest of each image and r cf (x,y;λ p ) is greater than the set value, forming a high signal-to-noise ratio pixel set. The single pixel with coordinates (x, y) and the surrounding pixels centered on the single pixel are divided into K pixel groups. The number of pixels in each pixel group is greater than 1. The average pixel value I of each pixel group at the acquisition time t is calculated. ka (x,y;t;λ p ), k=1,…,K;

[0086] The k value of the pixel group is defined according to the distance between the pixels in the common region of interest and the single pixel with coordinates (x, y), wherein the pixel group of k=1 is composed of all pixels in a square formed with a single pixel with coordinates (x, y) as the center and a set length as the side length; the pixel group of k=2 is composed of all pixels in a square ring with a set width formed along the outer periphery of the pixel group of k=1 with a single pixel with coordinates (x, y) as the center; the pixel group of k=3 is composed of all pixels in a square ring with a set width formed along the outer periphery of the pixel group of k=2 with a single pixel with coordinates (x, y) as the center; ... the pixel group of k=K is composed of all pixels in a square ring with a set width formed along the outer periphery of the pixel group of k=K-1 with a single pixel with coordinates (x, y) as the center.

[0087] 11) The average pixel value I of the K pixel groups of high signal-to-noise ratio pixels located at (x, y) obtained in step 10) ka (x,y;t;λ p ), calculate the average pixel standard score value I of each pixel group according to the following formula nka (x,y;t;λ p ):

[0088]

[0089] Among them, μIka (x,y;λ p ) and σ Ika (x,y;λ p ) are the average pixel values ​​I ka (x,y;t;λ p )'s time average and standard deviation values; the average pixel standard score value I nka (x,y;t;λ p ) corresponds to the time vector I nka (x,y;λ p ) is used as the input time vector in formula (3), and the average signal source, that is, the average independent component time vector, is calculated according to steps 12) to 15), which is recorded as S ka (x,y;λ p ), in the average independent component time vector S ka (x,y;t;λ p ) to determine the average pulse wave independent component signal S ca (x,y;t;λ p ), using filtering algorithm, calculate the average pulse wave time domain signal S of the pixel at (x, y) caf (x,y;t;λ p ), calculate and compare the time domain waveform phases of pixels at different positions in the high signal-to-noise ratio pixel set, and exclude high signal-to-noise ratio pixels whose phase changes exceed 90 degrees compared with adjacent high signal-to-noise ratio pixels;

[0090] 12) The average pulse wave time domain signal S of each pixel in the high signal-to-noise ratio pixel set obtained according to steps 10) and 11) caf (x,y;t;λ p ), calculating the propagation time and propagation time delay of two pixels located at different coordinate positions with the same time domain waveform peak phase, calculating the spatial distance between the two coordinate positions, and calculating the pulse wave propagation velocity at the corresponding waveform peak between the two positions based on the ratio of the spatial distance to the time delay and the arterial blood flow and pressure wave propagation direction;

[0091] 13) For each pair of pixels in the high signal-to-noise ratio pixel sets obtained in steps 10) and 11), calculate the pulse wave propagation velocity of all peaks for each pair within the video acquisition time T according to step 12), thereby obtaining the spatial distribution, temporal distribution, mean value of the spatial and temporal distribution, standard deviation of the spatial and temporal distribution, and related statistical characteristic parameters of the pulse wave propagation velocity within the common region of interest. The relevant statistical characteristic parameters include mean value, standard deviation, skewness coefficient, and kurtosis.

[0092] The pixel group independent component analysis method described in the present invention is derived from the independent component analysis algorithm well known in the industry, that is, the independent component analysis algorithm is applied to the data analysis of the iPPG video signal. The functional relationship between each pixel value and the imaging time variable t in the extracted common interest region is analyzed. According to the independent component analysis, the time-varying component corresponding to the pulse wave frequency in the pixel value signal is extracted and recorded as the pulse wave independent component S c (x,y;t;λ p ). Formula (3) is the matrix expression of the independent component analysis method of the pixel group, where the left side of the equation is composed of K time vectors I nk (x,y;λ p ), the right side of the equation is the mixing matrix [A] and K time vectors S k (x,y;λ p ), where K is an integer greater than 1 and k is any integer between 1 and K. If the spatial distribution image data of the pulse wave signal with high spatial resolution is calculated, the standard fraction value of the single pixel at (x, y) is taken as the input time vector I n1 (x,y;λ p ), while I n2 (x,y;λ p ) to I nK (x,y;λ p ) is taken as the average pixel standard score time vector of the K-1 surrounding pixel groups centered at the pixel at position (x, y) obtained according to formula (2). If the average pulse wave time domain signal data with a high signal-to-noise ratio is calculated, then the K input time vectors I nk (x,y;λ p ) are the time vectors of the average pixel standard score value. S k (x,y;λ p ) represents the time vector of the kth signal source, such as the independent component S1(x,y;λ p ) can represent the time change of pixel value caused by heartbeat or pulse as the signal source, S2(x, y; λ p ) can represent the temporal variation of pixel values ​​caused by breathing as a signal source. The mixing matrix [A] is related to the interaction between human tissue and incident light in the irradiated area and represents the temporal variation of pixel values ​​caused by the K signal sources through this interaction process.

[0093] The input signal time vector I of the pixel group independent component analysis method of the present invention is nk (x,y;λ p The number of ) is K. The setting of K value needs to consider the number of possible signal sources. If it is unknown, K value cannot be too small. However, K value should not be too large. If it is greater than the number of possible signal sources, it will lead to the iterative solution calculation of [A] and S in the independent component analysis of pixel groups.k (x,y;λ p ) takes too long, which increases the computational cost. One implementation of the pixel group independent component analysis method of the present invention may set K=3, but K may also be set to other positive integer values ​​such as 2 or 4. Figure 5 The image processing calculation flow based on the independent component analysis method of pixel groups with a K value of 3 is shown. The time vector I of each pixel in the common area of ​​interest is nk (x,y;λ p ) starts, obtains the input time vector required for independent component analysis of the pixel group, and calculates the pulse wave independent component time vector S point by point c (x,y;λ p ) and, after Fourier transformation and filtering, obtain the pulse wave spatial distribution map. Then, a set of high signal-to-noise ratio pixels is selected, and the average pulse wave time-domain signal waveform data is calculated after excluding high signal-to-noise ratio pixels with a phase shift of more than 90 degrees compared to adjacent high signal-to-noise ratio pixels.

[0094] The pixel group independent component analysis algorithm of the present invention can obtain the spatial distribution image data of the pulse wave signal with high spatial resolution; an implementation method of calculating the time vector of the input pixel group can be seen Figure 6 , Figure 6 Taking K=3 as an example, a single pixel standard score value at position (x, y) and two average pixel standard scores consisting of surrounding pixels centered at the pixel at (x, y) are selected. Figure 6 The right image is an enlarged view of the area within the rectangular box in the left image, showing the three pixel value time vectors I of one pixel at the wrist. nk (x,y;λ p ), where the time vector k = 1 is taken from the standard score value of a single pixel at (x, y), given by Figure 6 The small black square in the right figure indicates that the time vector k = 2 is taken from the average pixel standard score of the 24 nearest neighboring pixels centered at the pixel at (x, y). Figure 6 The light-colored square ring in the right figure indicates that the time vector k = 3 is taken from the average pixel standard score of the 24 next-nearest neighboring pixels centered at the pixel at (x, y) and outside the 24 nearest neighboring pixels, and is given by Figure 6 This is represented by the large black square ring frame in the figure on the right.

[0095] The pixel group independent component analysis algorithm of the present invention can also obtain an average pulse wave time domain signal with a high signal-to-noise ratio by increasing the number of pixels included in the average pixel value calculation. If K = 3, calculate the average pixel standard score value time vector I of the input pixel group nka (x,y;λ p) is implemented as follows: the time vector k=1 is taken from the average pixel standard score value of 81 pixels within a square with a side length of 9 pixels centered at the pixel at (x, y), the time vector k=2 is taken from the average pixel standard score value of the 40 nearest neighboring pixels centered at the pixel at (x, y) outside the above square with a side length of 9 pixels, and the time vector k=3 is taken from the average pixel standard score value of the 48 next nearest neighboring pixels centered at the pixel at (x, y) outside the above 121 surrounding pixels.

[0096] The first of the two implementation methods mentioned above is used to perform pixel group independent component analysis calculation, and the pulse wave independent component S c (x,y;t;λ p ) has a higher spatial resolution, but a lower signal-to-noise ratio; while using the second method, the obtained pulse wave independent component S ca (x,y;t;λ p ) has a high signal-to-noise ratio, but a relatively low spatial resolution. The method for calculating the time vector of the input pixel group described in the present invention can also be implemented in other ways according to the characteristics of the iPPG video signal.

[0097] Figure 7 are the 10 different incident light peak wavelengths λ obtained after imaging the left palm of subject 1 p Pulse wave spatial distribution diagram;

[0098] Ten light sources made of different light emitting diode arrays are used to provide 10 beams of λ p The incident light with different values ​​and average wavelength bandwidth Δλ value of about 36 nanometers is used to collect iPPG video data, and the pulse wave spatial distribution map is extracted by image processing software based on pixel group independent component analysis algorithm, where the first input time vector I n1 (x,y;λ p ) is the standard score value of a single pixel at (x, y), the second and third input time vectors I of the pixel group independent component analysis n2 (x,y;λ p ) and I n3 (x,y;λ p ) are the average pixel standard score value time vectors of the 24 nearest neighbor pixels and the 24 next nearest neighbor pixels centered at the (x, y) pixel position, and the numbers marked in the figure are the peak wavelengths λ p The spatial distribution of all 10 pulse waves is calculated using the pulse wave index r cf (x,y;λ p ) grayscale image after normalization of maximum and minimum values, pulse wave index r after standard score cf(x,y;λ p The corresponding relationship between the ) value and the grayscale is shown in the grayscale bar on the right side of the figure.

[0099] Figure 8 are the 10 different incident light peak wavelengths λ obtained after imaging the left palm of subject 2 p The spatial distribution diagram of the pulse wave, its data collection and processing are similar to Figure 7 same;

[0100] Figure 7 and Figure 8 The iPPG video data collected from the left palms of two different subjects using incident light of 10 different peak wavelengths were compared, and the pulse wave spatial distribution maps obtained after processing using the pixel group independent component analysis method described in this invention were compared. Comparing the spatial distribution of pulse wave signals with incident light of different peak wavelengths, it can be seen that the pulse wave signal intensity distribution is basically consistent with the arterial anatomical structure in the palm skin tissue. Comparing the results of different peak wavelengths; it can be seen that the pulse wave signal intensity distribution is at the peak wavelength of the incident light λ p The enhancement effect is more pronounced at wavelengths of 530 and 850 nanometers. The pulse wave signal is strongest at the radial and ulnar arteries at the wrist, the superficial palmar arch in the middle of the palm, and the arterioles and capillaries at the fingertips. By analyzing the spatial distribution of the pulse wave signal, we can select locations with strong pulse wave signals for further analysis of their time domain waveforms. Figure 9 Shows the incident light peak wavelength λ p Figure 1 shows the spatial and temporal distribution of the average pulse wave signal at 530 nm and 850 nm for Subject 1. The time domain signal waveforms were compared using different implementations of the pixel group time vector calculation method at four different pixel locations on the radial and ulnar arteries of the wrist. Figure 9 (B) The average pixel standard score value of a larger area is used to calculate the input pixel group time vector of independent component analysis, while Figure 9 The corresponding calculation of (C) uses the average pixel standard score value of the smaller area. Figure 9 (B) and Figure 9 (C), it can be seen that using the average pixel standard score value of a larger area to calculate the input pixel group time vector can effectively improve the signal-to-noise ratio of the pulse wave signal, and the incident light peak wavelength λ p When the wavelength is 850 nm, the signal-to-noise ratio of the pulse wave signal is higher. Figure 9 (B) The average pulse wave time domain signal waveform on the right side also shows that the pulse wave signals at different positions have a fixed phase relationship with the electrocardiogram signal, which proves that the independent component signal S obtained by the pixel group independent component analysis algorithm of the present invention is c (x,y;t;λ p) is a pulse wave signal, that is, its signal source is a pulse wave driven by heartbeat.

[0101] Figure 9 It is the high signal-to-noise ratio pixel-averaged pulse wave time domain signal at four different positions when the incident light peak wavelength is 530 nanometers and 850 nanometers;

[0102] Figure 9 (A) shows two different incident light peak wavelengths λ obtained after imaging the left palm of subject 1. p The pulse wave spatial distribution diagram, its data collection and processing steps are the same as Figure 7 The number marked in the lower left corner of the figure is λ p Value, other numbers are pixel position numbers, pulse wave index r after standard fraction cf (x,y;λ p The corresponding relationship between the ) value and the grayscale can be seen in the grayscale bar between the two figures.

[0103] Figure 9 (B) The left side of the figure shows the peak wavelength λ of the incident light p The average pulse wave time domain signal waveform of high signal-to-noise ratio pixels at four different locations on the radial artery and ulnar artery of subject 1's left wrist at a wavelength of 530 nanometers. Figure 9 (B) The right side of the figure shows the peak wavelength λ of the incident light p The average pulse wave time domain waveform of high signal-to-noise ratio pixels at four different locations on the left palm of subject 1 with a value of 850 nm is marked with different lines. Each pixel position number corresponds to Figure 9 The corresponding position numbers of the left and right images in (A) are the same as those of the pixel group independent component analysis algorithm of the two image series, that is, the first input time vector I n1 (x,y;λ p ) is the average pixel standard fractional time vector of the 81 central and nearest neighboring pixels centered at the pixel (x, y), and the second input time vector I n2 (x,y;λ p ) is the average pixel standard fractional time vector of the 40 nearest neighboring pixels centered at the pixel at (x, y), and the third input time vector I n3 (x,y;λ p ) is the average pixel standard fractional time vector of the 48 next-nearest neighboring pixels centered at the pixel at (x, y), Figure 9 (B) The right side of the figure also includes the electrocardiogram time domain signal measured simultaneously by inventor 1, marked with a black solid line, referred to as ECG.

[0104] Figure 9 (C) The left side of the figure shows the peak wavelength λ of the incident light pThe average pulse wave time domain waveform of high signal-to-noise ratio pixels at four different locations on the radial artery and ulnar artery of subject 1's left wrist at a value of 532 nanometers. Figure 9 (C) The right side of the figure shows the peak wavelength λ of the incident light p The average pulse wave time domain signal waveform of high signal-to-noise ratio pixels at four different locations on the radial artery and ulnar artery of subject 1's left wrist with a value of 850 nanometers. Each pixel position number corresponds to Figure 9 The corresponding position numbers of the left and right images in (A) are the same as those in the image processing method for obtaining the average pulse wave time domain signal waveform at different pixel positions. Figure 9 (B) The image processing and annotation methods are the same;

[0105] Figure 10 It is the average pulse wave time domain signal of high signal-to-noise ratio pixels located at 12 different positions on the left palm of subject 1 under the condition of incident light peak wavelength of 850 nanometers;

[0106] Figure 10 (A) is the calculated pulse wave spatial distribution diagram, and its data acquisition and processing steps are the same as Figure 6 The number marked in the lower left corner of the figure is λ p value, the other numbers are pixel position numbers, pulse wave index r cf (x,y;λ p ) value and grayscale correspondence Figure 8 same,;

[0107] Figure 10 (B) is the average pulse wave time domain waveform of high signal-to-noise ratio pixels at four different locations of the radial artery at the wrist. Each pixel position number corresponds to Figure 10 The corresponding wrist position numbers in (A) are the same as those in (B). The image processing methods for obtaining the pulse wave time domain signal waveforms at different pixel positions are the same as those in (C). Figure 9 (B) The image processing method is the same as the one in the annotation;

[0108] Figure 10 (C) is the average pulse wave time domain signal waveform of 8 high signal-to-noise ratio pixels at different locations in the palmar arch and fingertip arterioles. Each pixel position number corresponds to Figure 10 (A) The corresponding position numbers of the middle palm and fingertips are the same as the image processing method to obtain the average pulse wave time domain signal waveform at different pixel positions. Figure 9 (B) The image processing method is the same as that of the annotated image.

[0109] observe Figure 9 and Figure 10 It can also be found that Figure 9 (B) and Figure 9(C) 4 high signal-to-noise ratio pixels at different locations are selected, 2 of which are close to each other, either at the ulnar artery of the wrist (locations 1 and 2) or at the radial artery (locations 3 and 4). The average pulse wave time domain signal waveform phases of the two close locations are almost opposite, while the waveform phases of locations 1 and 4, and locations 2 and 3 at different arteries are almost the same. The same phase relationship can also be found in Figure 10 (B) is seen, that is, Figure 10 (A) The phases between positions 1 and 2, and 3 and 4 at the radial artery of the wrist are almost opposite, while the phases between positions 1 and 3 and 2 and 4 are almost the same. Figure 10 The waveform phases at other palm positions shown in (C), namely positions 5 through 12, do not exhibit nearly opposite phases. These data can be understood using the following model: the average pulse wave time-domain signal waveform extracted from iPPG video data is not only related to the interaction between the heartbeat and light and human tissue, but also considers the volume changes in blood vessels during pulse wave propagation, namely, vascular expansion and contraction. The associated vascular wall movement causes movement of adjacent skin tissue. The radial and ulnar arteries at the wrist have larger diameters and are close to the skin surface. Therefore, vascular wall movement results in significant skin surface movement. Because the skin surface is the interface between air and skin tissue, the optical refractive index difference between the two sides is the largest. According to the Fresnel equation in classical electromagnetic wave theory, an optically rough skin surface produces strong diffuse reflections. Consequently, the significant skin surface movement of the radial and ulnar arteries at the wrist causes corresponding changes in pixel values, resulting in nearly opposite phases of pulse waves at different adjacent positions. The relationship between this phase change and vascular volume change is extremely complex, and existing hemodynamic and tissue optical models are unable to quantify it. Therefore, according to the above model, when comparing the pulse wave phase relationship of high signal-to-noise ratio pixels at different locations near the artery and calculating the pulse wave propagation velocity, high signal-to-noise ratio pixels with phase differences greater than 90 degrees or even almost opposite phase waveforms should be excluded, otherwise the calculation error of time delay will occur. If the diameter of the blood vessels in the artery is small, such as the capillaries or arterioles in the middle of the palm and fingertips, or the blood vessels are deep under the skin, the skin surface movement will be very small and therefore will not affect the pulse wave phase, such as Figure 10 (C) shows that there is no need to consider the above issues.

[0110] Figure 11 Compared Figure 10 The delay time of the average pulse wave time domain signal of high signal-to-noise ratio pixels at three different positions is shown;

[0111] Figure 11 (A) shows the peak wavelength λ of the incident light pThe average pulse wave time domain signal waveform of high signal-to-noise ratio pixels at three different positions on the left palm of subject 1 was obtained after imaging the left palm of the subject at an optical wavelength of 850 nanometers. Each pixel position number corresponds to Figure 10 (A) The corresponding positions of the palm are numbered, and the pulse wave time delays of the three palm pixel positions are much smaller than the average heartbeat pulse period. Pixel position 1 is located at the radial artery of the wrist, position 7 is located at the superficial palmar arch in the middle of the palm, and position 11 is located at the arteriole of the fingertip. The image processing method for obtaining the average pulse wave time domain signal waveform at different pixel positions is the same, all of which are consistent with Figure 9 (B) The image processing method is the same as the one in the annotation;

[0112] Figure 11 (B) with Figure 11 (A) The pulse wave average time domain signal waveform data is the same, but the horizontal axis time t value range is reduced from 10 seconds to 0.4 seconds to accurately display Figure 11 (A) Waveforms near the second peak and valley of the pulse wave at three palm positions. The three vertical dashed lines represent the time positions of the second peak and valley of the averaged pulse wave time-domain waveform at the three different positions. The figure shows that the pulse wave peak and valley arrive at position 7 on the superficial palmar arch in the middle of the palm later than at position 1 on the radial artery at the wrist, and the pulse wave peak and valley arrive at position 11 on the arterioles at the fingertips later than at position 7.

[0113] Figures 7 to 10 The experimental data shown in the figure proves that the pulse wave spatial distribution map calculated by the pixel group independent component analysis method described in the present invention can first calculate the high signal-to-noise ratio pixel set and exclude the pixels in which the pulse wave appears to have an almost inverted waveform. Then, the delay and propagation speed caused by the pulse wave propagating in the artery are calculated between the other pixel positions. Figure 11 (A) The average pulse wave time domain signal waveform compared is selected from Figure 10 Among the three positions in (A), position 1 is located at the radial artery at the wrist, and its waveform is similar to that of position 7 located at the superficial palmar arch of the middle palm artery and position 11 at the fingertip artery, rather than being almost inversely proportional. Figure 11 (B) It can be seen that the pulse wave time domain signal waveforms at these three locations have a time delay caused by the finite propagation speed. Comparing the trough times in the figure, it can be seen that the pulse wave first reaches position 1 at the radial artery at the wrist, then delays to position 7 at the superficial palmar arch of the middle artery, and finally reaches position 11 at the arteriole of the fingertip. This is consistent with the direction of blood flow and pulse wave propagation, which starts from the heart and flows through the arteries to the limbs.

[0114] according to Figure 10The pulse wave spatial distribution diagram of subject 1 and the average pulse wave time domain signal waveform of the high signal-to-noise ratio pixels at the selected 12 positions are shown. The pulse wave propagation velocity between the two areas of the palm can be calculated. First, positions 2 and 3, where almost antiphase waveforms appear at the radial artery of the wrist, are excluded. Then, all different pairwise high signal-to-noise ratio pixel combinations between the two areas are selected, such as position 1 of the radial artery at the wrist and position 5 of the superficial palmar arch of the middle artery, position 1 and position 6, position 1 and position 7, and so on. The time delay of the same peak of a certain two-position combination is calculated. Then, the time delay of all peaks is calculated. Finally, the average and standard deviation of the time delay of all peaks of all pairwise position combinations between the two areas are calculated, and recorded as Δt a,m and Δt a,s .according to Figure 10 (A) The average straight-line distance between the two selected positions is calculated by calibrating the pixel size, which is recorded as d. The average pulse wave velocity between the two positions is v = d / Δt a,m Table 1 gives the Figure 10 The pulse wave signal of subject 1 shown here shows the average distance between different pixels in two different areas of the palm, the average pulse wave time delay, and the average wave velocity data. These data range from 1 to 10 meters per second, which is consistent with the average wave velocity range of the human pulse wave measured using other methods.

[0115] Table 1. Average straight-line distance, average pulse wave time delay, and average wave velocity between two different areas of the palm

[0116]

[0117] according to Figure 10 The calculated data is obtained from the pulse wave signal shown.

Claims

1. An image processing algorithm for calculating human pulse frequency from iPPG video data based on an imaging system for contactless detection of pulse wave spatial and temporal distribution and characteristics, the imaging system for contactless detection of pulse wave spatial and temporal distribution and characteristics comprising a computer (8) for storing, displaying and processing image data, an imaging device connected to the computer (8) via a data transmission line (7) for acquiring M images of a human tissue region (3), and a light source assembly for providing incident light, characterized in that: The imaging device comprises a video camera (6), an imaging lens (5) and an imaging polarizer (4) for receiving the reflected light of the human tissue area (3), which are sequentially arranged in the same direction. The imaging polarizer (4) limits the imaging lens (5) to only receive reflected light with a linear polarization state and a polarization direction perpendicular to the polarization direction of the incident light. The M images acquired by the imaging device are video images constituting iPPG video data of the human tissue area (3) extracted therefrom. The light source assembly comprises 1 to 10 light source groups arranged at equal angles and intervals around the imaging device for irradiating the human tissue area (3). Each light source group comprises a light source (1) and a light source polarizer (2) coaxially arranged with the light source (1). The light irradiating the human tissue area (3) is made into incident light with a linear polarization state by adjusting the angle of rotation of the light source polarizer (2) around the coaxial axis. The image processing algorithm is characterized in that the image processing algorithm comprises the following steps: 1) Determine whether each image in the iPPG video data collected from the M images contains saturated pixels. If each image does not contain saturated pixels, proceed to step 2); if more than one image contains saturated pixels, select the image with the most saturated pixels and calculate the number of saturated pixels p of this image. s The total number of pixels p of the image t and determine whether the ratio is greater than 2%. If so, the image processing process ends, otherwise, proceed to step 2); 2) Based on the pixel value histogram of each image in the iPPG video data, calculate the boundary threshold of the region of interest with human blood vessel and tissue characteristics in the corresponding image. Based on the region of interest of each image, find the common region of interest of all images and intercept the corresponding common region of interest in each image; 3) According to the common region of interest of each image, the average pixel value of the common region of interest of each image in the iPPG video data at the acquisition time t is calculated, which is recorded as I av (t;λ p ), where 0≤t≤T, T is the video acquisition time, λ p is the peak wavelength of the incident light, and calculates the average pixel value I av (t;λ p ) time average μ av (λ p ) and standard deviation σ av (λ p ), calculate the average pixel standard score ΔI of the common area of ​​interest in each image according to the following formula av (t;λ p ),Right now 4) Use Fourier transform algorithm to calculate the average pixel standard score value ΔI av (t;λ p ) of the frequency domain signal Δi av (f;λ p ), determine the pulse frequency f of a given human body h The pulse frequency measurement value f of the strongest component within the range c , the pulse frequency measurement value f c Refers to the given human pulse frequency f h The frequency with the largest Fourier component value within the range, the human pulse frequency f h The range is 0.6 Hz to 4.0 Hz.

2. A method for calculating the image processing algorithm for calculating the human pulse frequency according to claim 1, and a method for calculating the independent component analysis of pixel groups of the spatial distribution image and time domain distribution data of the human pulse wave signal, characterized in that: The steps include: 1) Based on the standard deviation normalization method, the single-pixel standard score value I of the common area of ​​interest in each image in the iPPG video data is calculated according to the following formula: n (x,y;t;λ p ),Right now Where x and y are the coordinates of a single pixel position; I(x, y; t; λ p ) is the single pixel value with coordinates (x, y); μ I (x,y;λ p ) and σ I (x,y;λ p ) are the temporal average and standard deviation of the single pixel value with coordinates (x, y); λ p is the peak wavelength of the incident light; 2) Select a single pixel with coordinates (x, y), divide the surrounding pixels centered on the pixel into K-1 surrounding pixel groups, and calculate the average pixel standard score value I of each surrounding pixel group nz (x, y; t; λ0), using the average pixel standard score value I of the surrounding pixel group nz The values ​​of (x, y; t; λ0) at different times are used as components, expressed as the average pixel standard score value time vector I nz (x,y;λ p ), where z = 2, ..., K, where the value of K is the average pixel standard score value I affecting each pixel group nz The number of signal sources (x, y; t; λ0) that change with time is the same; 3) Calculate the signal source S of each single pixel with coordinates (x, y) k (x,y;t;λ p )’s signal source time vector S k (x,y;λ p ), k=1,…,K, wherein the signal source time vector S k (x,y;λ p ) must satisfy the condition of maximum temporal statistical independence between the time vectors of each signal source and is defined as the output time vector of the independent component analysis algorithm of the pixel group, called the independent component time vector; First, the single pixel standard score value I with coordinates (x, y) n (x,y;t;λ p ) is represented as a single pixel standard score value time vector I n1 (x,y;λ p ), and the average pixel standard score value time vector I of the K-1 surrounding pixel groups nz (x,y;λ p ) are combined into K time vectors I nk (x,y;λ p ), where k = 1, ..., K, set I nk (x,y;λ p ) and the time vector S of the signal source that affects the time variation k (x,y;λ p ) is a linear relationship, which can be expressed by the following matrix formula: Where [A] is an unknown matrix, also called a mixing matrix; the average pixel standard score value time vector I of the surrounding pixel group is converted by the pixel group independent component analysis algorithm. nz (x,y;λ p ) and single pixel standard score value I n (x,y;t;λ p ) is used as the input time vector of the pixel group independent component analysis algorithm, and the matrix [A] and the signal source time vector S that satisfy formula (3) are iteratively solved. k (x,y;λ p ); 4) Using Fourier transform algorithm, calculate each independent component time vector S k (x,y;λ p ) or the corresponding signal source S k (x,y;t;λ p ) spectrum distribution, determine the signal source S k (x,y;t;λ p ) The maximum component frequency value f0 within the human pulse frequency range, the K independent component time vectors have the closest pulse frequency measurement value f c The independent component time vector of the f0 value is the pulse wave independent component time vector, denoted as S c (x,y;λ p ), the pulse wave independent component time domain signal is S c (x,y;t;λ p ), is the pulse wave independent component time domain signal of the pixel value of a single pixel with coordinates (x, y); 5) Using filtering algorithm, calculate the pulse wave independent component time domain signal S c (x,y;t;λ p ) The filtered signal is recorded as the pulse wave time domain signal S cf (x,y;t;λ p ) or pulse wave time vector S cf (x,y;λ p ); The filtering algorithm is based on Fourier transform in the frequency domain, or adopts a digital filtering method such as a finite impulse response filtering algorithm; 6) Based on Fourier transform, calculate the pulse wave time domain signal S of each single pixel with coordinates (x, y) cf (x,y;t;λ p ) is the frequency distribution signal, recorded as the pulse wave frequency domain signal s cf (x,y;f;λ p ); 7) According to the pulse wave frequency domain signal s cf (x,y;f;λ p ) and the following formula to calculate the pulse wave index r of the single pixel value with coordinates (x, y) one by one cf (x,y;λ p ),Right now Where Δf = F / 2d is the frequency domain step size, F is the frame rate of the video camera, d is the smallest integer that satisfies the condition M≤2d≤2M, M is the number of images, N≥2, W≥N are the number of terms in the sum of the numerator and denominator in formula (4), respectively; 8) Repeat steps 4) to 7) until the pulse wave time domain signal S of all pixels in the common region of interest in each image is completed. cf (x,y;t;λ p ) and the pulse wave index r of all pixels cf (x,y;λ p ) calculation; 9) The pulse wave index r of all pixels in the common region of interest in each image obtained in step 8) cf (x,y;λ p ), calculate the maximum and minimum values ​​of the pulse wave index in the common area of ​​interest, and then calculate the pulse wave index r based on the maximum and minimum values cf (x,y;λ p ) normalized value, according to the pulse wave index r cf (x,y;λ p )'s (x, y) coordinates and normalized values, generate and output each image in the common region of interest and with the incident light peak wavelength λ p Pulse wave signal spatial distribution image data; 10) According to the pulse wave signal spatial distribution image data in the common region of interest of each image obtained in step 9), select the image with coordinates (x, y) in the common region of interest of each image and r cf (x,y;λ p ) is greater than the set value, forming a high signal-to-noise ratio pixel set. The single pixel with coordinates (x, y) and the surrounding pixels centered on the single pixel are divided into K pixel groups. The number of pixels in each pixel group is greater than 1. The average pixel value I of each pixel group at the acquisition time t is calculated. ka (x,y;t;λ p ), k=1,…,K; 11) The average pixel value I of the K pixel groups of high signal-to-noise ratio pixels located at (x, y) obtained in step 10) ka (x,y;t;λ p ), calculate the average pixel standard score value I of each pixel group according to the following formula nka (x,y;t;λ p ): Among them, μ Ika (x,y;λ p ) and σ Ika (x,y;λ p ) are the average pixel values ​​I ka (x,y;t;λ p )'s time average and standard deviation values; the average pixel standard score value I nka (x,y;t;λ p ) corresponds to the time vector I nka (x,y;λ p ) is used as the input time vector in formula (3), and the average signal source, that is, the average independent component time vector, is calculated according to steps 3) to 6), which is recorded as S ka (x,y;λ p ), in the average independent component time vector S ka (x,y;t;λ p ) to determine the average pulse wave independent component signal S ca (x,y;t;λ p ), using filtering algorithm, calculate the average pulse wave time domain signal S of the pixel at (x, y) caf (x,y;t;λ p ), calculate and compare the time domain waveform phases of pixels at different positions in the high signal-to-noise ratio pixel set, and exclude high signal-to-noise ratio pixels whose phase changes exceed 90 degrees compared with adjacent high signal-to-noise ratio pixels; 12) The average pulse wave time domain signal S of each pixel in the high signal-to-noise ratio pixel set obtained according to steps 10) and 11) caf (x,y;t;λ p ), calculating the propagation time and propagation time delay of two pixels located at different coordinate positions with the same time domain waveform peak phase, calculating the spatial distance between the two coordinate positions, and calculating the pulse wave propagation velocity at the corresponding waveform peak between the two positions based on the ratio of the spatial distance to the time delay and the arterial blood flow and pressure wave propagation direction; 13) For each combination of two pixels in the high signal-to-noise ratio pixel set obtained in steps 10) and 11), calculate the pulse wave propagation velocity of all peaks of each combination within the video acquisition time T according to step 12), thereby obtaining the spatial distribution, temporal distribution, mean values ​​of the spatial distribution and temporal distribution, standard deviations of the spatial distribution and temporal distribution, and related statistical characteristic parameters of the pulse wave propagation velocity within the common region of interest.

3. The pixel group independent component analysis method according to claim 2, characterized in that: The average pixel standard score value I of each pixel group in step 2) nz (x,y;t;λ p ) is defined according to the distance from a single pixel with coordinates (x, y), wherein the peripheral pixel group with z=2 is composed of all pixels in a square formed with a set length as the side length and centered on the single pixel with coordinates (x, y) excluding the single pixel with coordinates (x, y); the peripheral pixel group with z=3 is composed of all pixels in a square ring with a set width formed along the outer periphery of the peripheral pixel group with z=2 and centered on the single pixel with coordinates (x, y); .....and the peripheral pixel group with z=K is composed of all pixels in a square ring with a set width formed along the outer periphery of the peripheral pixel group with z=K-1 and centered on the single pixel with coordinates (x, y).

4. The pixel group independent component analysis method according to claim 2, characterized in that: In step 10), the pixel group k value is defined according to the distance between the pixels in the common region of interest and the single pixel with coordinates (x, y), wherein the pixel group of k=1 is composed of all pixels in a square formed with a single pixel with coordinates (x, y) as the center and a set length as the side length; the pixel group of k=2 is composed of all pixels in a square ring with a set width formed along the outer periphery of the pixel group of k=1 with a single pixel with coordinates (x, y) as the center; the pixel group of k=3 is composed of all pixels in a square ring with a set width formed along the outer periphery of the pixel group of k=2 with a single pixel with coordinates (x, y) as the center; ... the pixel group of k=K is composed of all pixels in a square ring with a set width formed along the outer periphery of the pixel group of k=K-1 with a single pixel with coordinates (x, y) as the center.

5. The pixel group independent component analysis method according to claim 2, characterized in that: The relevant statistical characteristic parameters of the pulse wave propagation velocity in the common region of interest described in step 13) include: mean value, standard deviation value, skewness coefficient, and kurtosis.

Images