A method for accurately extracting pulse waves based on human superficial artery video

Abstract

The present invention discloses a method for accurately extracting pulse waves based on human superficial artery video, which is used to solve the motion and light interference problems commonly existing in non-contact pulse wave measurement methods, and multiple filtering is performed on multiple noises to achieve accurate extraction of pulse waves. The present invention includes the following steps: S1: obtaining human superficial artery video; S2: preprocessing the video with adaptive gamma transform; S3: performing motion amplification on the video with an amplification algorithm combining Euler amplification and FRR filtering; S4: adaptive threshold inter-frame difference method and selecting multiple feature areas as regions of interest; S5: extracting the original pulse wave signal; S6: II-type Chebyshev, wavelet transform of adaptive threshold method, and adaptive notch filter filtering out noise in the signal; S7: modulation domain spectrum subtraction method enhancing the pulse wave signal. The present invention can effectively suppress the influence of motion artifacts and environmental brightness differences, and accurately select the region of interest to filter out noise, further improving the accuracy of pulse wave extraction.

Description

Technical Field

[0001] The present invention relates to a method for accurately extracting pulse waves based on human superficial artery videos, which can be used for non-contact measurement of pulse waves in multiple superficial arteries of the human body. Vital signs such as heart rate, blood pressure, and blood oxygen saturation can be analyzed through the detection and recording of pulse waves, and belongs to the field of medical equipment and computer vision technology. Background Art

[0002] The process of the heart completing a rhythmic contraction and relaxation is a cardiac cycle, and the process of completing a cardiac cycle is equivalent to the heart completing a blood pumping function. In the process of realizing the blood pumping function, the blood circulation system of the whole body plays a vital role. When the heart contracts, the human body's blood is ejected from the ventricle into the arterial vascular system, through the capillary network, to the venous system; when the heart relaxes, the blood reaches the arterial system through the venous system and then flows back to the atrium to prepare for the next blood ejection during the heart contraction. Therefore, the arterial system is the first link in the realization of the heart's blood pumping function, which can more clearly reflect the rhythmic movement of the cardiac cycle and the amount of blood pumped. The superficial arteries of the human body mainly include the superficial temporal artery, carotid artery, radial artery, brachial artery, femoral artery, dorsalis pedis artery, etc. These arteries have more obvious vascular pulsation signals caused by the contraction and relaxation of the heart.

[0003] Different tissues in the human body can absorb visible light to different degrees. The amount of blood in the blood vessels changes periodically due to the alternating contraction and relaxation of the heart. Therefore, the degree of absorption of incident light by blood in the arteries also changes periodically. However, various tissues such as muscles and bones are in a stable form in the human body, so the absorption state of incident light remains constant. Therefore, the fluctuation of the absorbance of incident light by the human body is mainly caused by the rhythmic contraction and relaxation of the heart.

[0004] Changes in blood volume in blood vessels, i.e. changes in blood volume, will cause changes in brightness in the video images of superficial human arteries. When the heart contracts, the blood volume in the blood vessels increases, and the brightness of the region of interest in the captured video image increases. When the heart relaxes, the blood volume in the blood vessels decreases, and the brightness of the region of interest in the captured video image decreases.

[0005] The existing non-contact pulse wave extraction schemes are inevitably affected by light and motion, which will reduce the accuracy of the extracted pulse wave. Among them, the light influence refers to the influence of weather conditions, different indoor and outdoor light intensity, different shooting angles, etc., which will affect the brightness difference of the reflected light caused by the image received during the video shooting process. Regarding the influence of movement, on the one hand, during the video shooting process, the subject will unconsciously shake slightly; on the other hand, when the human body's blood flows to various parts of the blood vessels, it will generate pressure waves of different sizes in different directions on the blood vessel wall. During the blood flow process, the pressure waves will affect the skin and cause periodic small movements on the skin surface.

[0006] In recent years, with the continuous rise of the medical industry and the continuous growth of the computer industry, the attention of computer vision in the medical field has gradually increased. In 1938, Hertzman proposed the photoplethysmography (PPG) technology, and this technology can analyze the heart rate of the subject from the acquired waveform. In 2000, ATWu et al. improved the traditional PPG technology and proposed the non-contact imaging photoplethysmography (IPPG) technology using computer knowledge, which has attracted in-depth research by a large number of scholars. Lei Hengbo et al.'s method for collecting photoplethysmographic pulse wave signals (patent number: 201610373750.2) obtains photoplethysmographic pulse wave signals by processing the skin video twice. Hao-YuWu et al. proposed the Euler amplification method to amplify the tiny movements of the face and radial artery caused by blood flow, providing a way to accurately extract the pulse wave signal converted from the video frame.

[0007] Li Jiangshan et al. visualized the pulse wave obtained after processing the radial artery video based on the Euler amplification video algorithm, selected the area with the highest signal-to-noise ratio as the region of interest in the signal after Euler amplification, used the independent component analysis algorithm to separate the effective pulse wave signal and noise, and finally used the blind source separation algorithm to extract the pulse wave. This scheme uses the independent component analysis method. Since the pulse signal of blood volume is not necessarily the first basis vector of the independent component component, there is randomness. Su Peiquan et al. also used the Euler algorithm to amplify the radial artery video in a non-contact heart rate measurement method and device based on Euler image amplification (patent number: 201710303637.1). The difference is that it uses the interception of the region of interest in the YCrCb color space, obtains the periodic pulse wave through the signal of skin brightness change, and calculates the heart rate. However, this method of selecting the region of interest has poor resistance to light interference.

[0008] Many scholars have based their research on contactless pulse wave detection technology on the video results after Euler amplification, but Euler amplification is sensitive to noise. While amplifying useful signals, it also amplifies noise, and endpoint noise may be generated during use. Gao Guanqun et al. proposed a phase-based pulse wave extraction method, replacing the Euler amplification algorithm with complex controllable pyramid decomposition, and using a time-domain bandpass filter to select the phase difference signal of the horizontal phase spectrum in the region of interest, and then filtering out the noise to extract the pulse wave. This complex controllable pyramid decomposition algorithm does not amplify noise, but can only amplify the edge position, and the calculation is complex. Meng Wen's method for real-time display of pulse waveforms based on face video (patent number: CN201911249939.0) combines the Euler amplification and adversarial network methods to express the pulse wave state of the subject as accurately as possible through model training, but in the process of model establishment and adversarial network construction, a large amount of data is required for training to ensure the accuracy of the adversarial network. Yu Haotian et al. proposed the EVM-FRR amplification algorithm, which is an extension and improvement of the Euler amplification algorithm. The one-dimensional time signal extracted by it no longer has phase distortion, that is, while amplifying the useful signal of the video frame, it almost does not amplify the noise. However, the impact of endpoint noise is not solved, and the impact of light interference is not considered. Summary of the invention

[0009] The purpose of the present invention is to provide a method for accurately extracting pulse waves based on human superficial artery video. The method can reduce the influence of light interference and motion interference to a certain extent, and greatly improve the signal-to-noise ratio to obtain accurate pulse wave signals.

[0010] The technical solution of the present invention is as follows:

[0011] The method for accurately extracting pulse waves based on human superficial artery video includes the following steps:

[0012] S1: Get the video of the superficial artery of the human body;

[0013] S2: Preprocess the video using adaptive gamma transform;

[0014] S3: Video motion amplification of human superficial arteries using the EVM-FRR amplification algorithm;

[0015] S4: Adopting the adaptive threshold frame difference method and selecting multiple feature regions to determine the region of interest (ROI);

[0016] S5: extracting the original pulse wave signal;

[0017] S6: Use type II Chebyshev filter, wavelet transform with adaptive threshold method, and adaptive notch filter to filter out the noise that may exist in the signal;

[0018] S7: Enhance the pulse wave signal using modulation domain spectral subtraction.

[0019] Step S1 obtains the position of superficial arteries of the human body, mainly including the superficial temporal artery, carotid artery, radial artery, brachial artery, femoral artery, dorsalis pedis artery, etc. The obvious pulsation positions are selected as ideal video acquisition positions.

[0020] Step S2 Adaptive gamma transformation is used for video preprocessing system to adjust the video output caused by the difference in light brightness by changing the γ parameter, which can eliminate the impact of image quality caused by uneven lighting. The γ parameter adjusts the following formula:

[0021]

[0022] Among them, X(x,y) is the input image; O(x,y) is the output image; γ is the gamma correction parameter; the size of the γ parameter is determined by the illumination component and the brightness mean of the image. The formula for determining the γ parameter is:

[0023]

[0024] Where m is the mean brightness of the illumination image; L(x,y) is the illumination component value at point (x,y).

[0025] In the video preprocessing system of adaptive gamma transform, the γ correction parameter corresponding to a certain point can be obtained from the brightness mean value of the image and the illumination component value at that point. When the correction parameter γ value is greater than 1, the local brightness of the input image is attenuated; when the correction parameter γ value is less than 1, the local brightness of the input image is enhanced; when the correction parameter γ value is equal to 1, the local brightness of the input image remains unchanged. By searching and judging each pixel, local brightness compensation can be completed to solve the impact of uneven illumination on image quality.

[0026] Step S3: The method for amplifying the motion of human superficial artery video based on the EVM-FRR amplification algorithm comprises the following steps:

[0027] S31: performing spatial filtering on the video frame, i.e. constructing a Laplacian pyramid to decompose the video frame into frequency space bands of different scales;

[0028] S32: performing time-domain bandpass filtering on the image layers with the same spatial resolution in the multiple spatial bands obtained after spatial filtering by using the following formula, wherein the formula is:

[0029]

[0030] S33: I(x, t) is the brightness of a frame image in the video when the spatial point position is x and the time point is t. f(x) is the brightness of the corresponding spatial point position x when time t = 0, that is, f(x) = I(x, 0). δ(t) is the motion displacement of the spatial point after time t.

[0031] S34: The tiny pulse signal to be amplified is determined by a time domain bandpass filter, that is, the brightness change signal B(x, t) corresponding to the tiny pulse signal when the spatial point position of the video frame image is x and the time point is t, wherein the brightness change signal B(x, t) is determined by the following formula:

[0032]

[0033] S35: The brightness change corresponding to the pulse beat signal of the video frame image at the spatial point x and the time point t is amplified by α times, and then superimposed on the image brightness signal I(x,t) at the spatial point x and the time point t in the original video frame, as shown in the following formula:

[0034]

[0035] Where α is the amplification factor, It is the image brightness signal that is amplified once and superimposed at the spatial point x and the time point t.

[0036] S36: For the image brightness signal of the first amplification superposition Perform time reversal at the spatial point x to obtain a reversed time signal at the spatial point x and time -t;

[0037] S37: taking the inverted time signal as input, amplifying the image brightness change signal corresponding to the pulse beating signal at the spatial position x and the time -t by α times through the amplification factor relationship;

[0038] S38: The tiny signal at the spatial point position x and the time point -t after the second amplification is time-reversed again to obtain a video frame signal at the spatial position x and the time point t That is the video frame signal after EVM-FRR amplification.

[0039] The EVM-FRR amplification algorithm combines the Euler amplification algorithm with the FRR filter. It uses two time reversals and two amplification schemes based on the IIR filter, which improves the noise sensitivity problem in the Euler algorithm. In the process of video signal amplification, no distorted waveform caused by phase distortion is generated.

[0040] Step S4 adopts the inter-frame difference method with adaptive threshold and mainly includes the following steps:

[0041] S41: Perform differential operation on the pixel grayscale of adjacent frames of the collected arterial video or two images separated by several frames. Since the target is constantly moving, the feature points in adjacent video frames have obvious changes relative to the static background. Therefore, it is easier to obtain the absolute value of the image pixel difference between the video frames, that is, the video frame image after differential.

[0042] S42: Select the video frame image threshold using an adaptive threshold algorithm. Assume that there are L gray levels in the video frame image and the binarization threshold is T. The video frame image is divided into two regions, K0 and K1. The gray values ​​of the pixels in the K0 region are in the range of [0, T], and the gray values ​​of the pixels in the K1 region are in the range of [T+1, L-1].

[0043] Therefore, the probability of gray level j appearing is:

[0044] Among them, n j Represents the number of pixels with gray level j, and N is the total number of pixels in a single video frame image.

[0045] The probability of the number of pixels in the K0 area appearing in the entire video frame image is:

[0046]

[0047] The average grayscale of the K0 region is:

[0048] The probability of the number of pixels in the K1 area appearing in the entire video frame image is:

[0049]

[0050] The average grayscale of the K1 area is:

[0051] The inter-class variance of the moving area and the static area in the video frame image is:

[0052] S 2 (T) = P0P1(E0-E1) 2 (11)

[0053] The threshold T when the inter-class variance is the largest is selected as the optimal binarization threshold.

[0054] S43: Binarization of video frame images. Compare the grayscale value of each pixel of the video frame image with the determined optimal binarization threshold. The pixel value of the point greater than the optimal threshold becomes 255, which is determined as a moving part; the pixel value of the point less than the optimal threshold becomes 0, which is determined as a stationary part.

[0055] S44: Perform connectivity analysis on the processed binary image, and then connect each frame of the image to form a new video, that is, obtain a moving part video.

[0056] Step S4 of selecting multiple feature regions includes the following steps:

[0057] S45: Using the SIFT algorithm to extract feature points of the moving part video frame image obtained by the inter-frame difference method, and determine the coordinates of the feature points.

[0058] S46: Among the determined feature points, at least three feature points are randomly selected to form a feature point group, the feature point group forms a polygon, and the centroid coordinates of the feature point group are calculated. The centroid calculation formula is as follows:

[0059]

[0060]

[0061] Among them, x and y are the horizontal and vertical coordinates of the center of mass, respectively. i ,y i They are the horizontal and vertical coordinates of each feature point (i=1,2,3)

[0062] S47: A rectangular area is intercepted with the centroid as the center as a single region of interest of the pulse wave signal in the video frame image, and six regions of interest are determined as regions of interest of the final pulse wave video frame image in the same way.

[0063] Step S5 extracts the original pulse wave signal, which mainly includes:

[0064] S51: Convert each frame of the video signal slice image of the region of interest from the RGB color space to the standard chromaticity model CIE-XYZ. The conversion formula is as follows:

[0065]

[0066] S52: Convert the slice of the region of interest of the CIE-XYZ standard colorimetric system to the Adobe RGB color space. The conversion formula is as follows:

[0067]

[0068] S53: Separate the RGB three channels of the region of interest slice in the Adobe RGB color space.

[0069] S54: performing mean processing on the video frame signal pixels of the region of interest in each channel, and combining them in time sequence to obtain three groups of one-dimensional time series.

[0070] S55: Convert each channel one-dimensional time series x i (t)(i=1,2,3) is normalized, and its expression is as follows:

[0071] Among them, μ i For x i The pixel mean of (t), σ i For x i (t) is the pixel standard deviation.

[0072] After the normalized signal is subjected to fast Fourier transform to convert the time domain signal into a frequency domain signal, a bandpass filter is used to perform a primary filtering on the frequency domain signal in the range of 0.8 to 3.3 Hz to obtain the original pulse wave.

[0073] Step S6 uses a type II Chebyshev filter, a wavelet transform with an adaptive threshold method, and an adaptive notch filter to filter out possible noise in the signal. The type II Chebyshev filter uses a monotonous passband to ensure that the data in the region of interest within the passband is not affected by the filter, while the noise in the stopband is quickly filtered out to reduce phase loss. The wavelet transform with an adaptive threshold method changes the reference variable And n, the adjustment function switches between soft threshold and hard threshold to filter out Gaussian white noise. For the adaptive notch filter, the frequency distribution is obtained by fast Fourier transform, the frequency points exceeding the amplitude threshold are matched with the notch filter center frequency, and the linear input signal of the input video frame image is separated from the noise signal part using the weight coefficient to achieve the purpose of eliminating certain specific frequency noises in the video frame signal.

[0074] The type II Chebyshev filter is filtered by a transfer function, where the transfer function formula is as follows:

[0075] Where W0 is the cut-off frequency, is a Chebyshev polynomial of order n.

[0076] The formula for the attenuation ∈ at different decibels χ in the transfer function is as follows:

[0077]

[0078] Regarding adaptive threshold wavelet transform, when When is 0 and n is 0, the wavelet transform is the wavelet transform of the soft threshold function. When set to a sufficiently large value, the wavelet transform is transformed into a wavelet transform of a hard threshold function. The wavelet transform switching function of the soft and hard threshold is as follows:

[0079]

[0080] in, is the wavelet decomposition coefficient after processing; m j,k is the wavelet decomposition coefficient before processing; and n are both positive numbers. The value of λ is determined by the following formula:

[0081]

[0082] Among them, MAD is the median value of the absolute value of the first layer wavelet, and N is the length of the video frame signal.

[0083] The adaptive notch filter makes a corresponding frequency distribution by finding the part exceeding the amplitude threshold in the frequency domain range of 0.8 to 3.3 Hz of the wavelet transformed signal, and determines the number k of notch points exceeding the amplitude threshold.

[0084] The adaptive notch filter uses the least mean square algorithm (LMS) as the adaptive algorithm. The output of the filtered signal

[0085] is: e(n)=d(n)-y(n) (21)

[0086] d(n) is the input signal, and y(n) is the output of the adaptive filter. The formula corresponding to y(n) is:

[0087] y(n)=ω1(n)x1(n)+ω2(n)x2(n)+…+ω 2k-1 (n)x 2k-1 (n)+ω 2k (n)x 2k (n) (22)

[0088] Among them, ω i (n) is the weight coefficient of the reference signal, x i (n) is the input of the reference signal;

[0089] The adaptive filtering algorithm uses the minimum mean square error algorithm, so that the weight coefficient is:

[0090] ω i (n+1)=ω i (n)+2μe(n)x i (n) (23)

[0091] Among them, μ is the length of the adjustment step, which controls the speed of adaptive convergence.

[0092] The reference signal input is an orthogonal signal of the center frequency. The reference signal formula is:

[0093]

[0094] Among them, ω0 is the center frequency, A is the filter amplitude, which can be calculated based on μA 2 To adjust the bandwidth when the signal amplitude reaches -3dB, so that the shoulder of the notch of the notch can be arbitrarily narrow. Different filter frequencies ω0 are configured according to the frequency distribution exceeding the threshold, and the bandwidth is adjusted by configuring the filter amplitude A.

[0095] Step S7 uses modulation domain spectrum subtraction to enhance the effective pulse wave signal, and its main steps include:

[0096] First, obtain the amplitude spectrum |X(p,r)| and phase spectrum ∠X(p,r) (where r is the discrete frequency) of the pulse wave signal after filtering; perform fast Fourier transform along time at each frequency point of the frequency domain amplitude spectrum of the pulse wave signal to obtain the modulation amplitude spectrum |X(δ,r,l)| and modulation phase spectrum ∠X(δ,r,l) of the pulse wave signal (where δ is the modulation frame and l is the modulation frequency). The modulation domain amplitude spectrum obtained by modulation spectrum subtraction is:

[0097]

[0098] Among them, β is used to determine the type of spectrum subtraction. When β=1, it is determined to be the amplitude spectrum subtraction type; when β=2, it is determined to be the energy spectrum subtraction type. η is the over-subtraction factor in the modulation domain, ξ is the spectrum suppression factor, that is, the lower limit of the enhanced modulation spectrum, which is between 0 and 1; |E(δ,r,l)| is the modulation amplitude spectrum estimation of the noise. By comparing the signal-to-noise ratio of the signal with the preset effective signal threshold, the part of the signal-to-noise ratio that exceeds the threshold is the modulation amplitude spectrum of the signal, otherwise it is the modulation spectrum of the noise.

[0099] in,

[0100] Wherein, N is the number of frames corresponding to the no-interest region (noise region) of the leading video frame.

[0101] Modulation phase spectrum after phase spectrum compensation enhancement:

[0102] ∠G(δ,r,l)=arg[Z′(δ,r,l)] (27)

[0103] Among them, Z′(δ,r,l) is the modulation spectrum after modulation domain spectrum subtraction and phase spectrum compensation enhancement, and its formula is as follows:

[0104] Z′(δ,r,l)=G g (δ,r,l)+Γ(δ,r,l) (28)

[0105] Among them, G g (δ,r,l) is the modulation domain amplitude estimate after spectral subtraction, and Γ(δ,r,l) is the modulation domain phase compensation.

[0106] The modulation domain amplitude estimate after spectral subtraction is: G g (δ,r,l)=|G(δ,r,l)|e j∠X(δ,r,l) (29)

[0107] Phase compensation degree: Γ(δ,r,l)=λφ(l)|E(δ,r,l)| (30)

[0108] Among them, λ is an empirical constant, φ(l) is an antisymmetric function, and its formula is as follows:

[0109]

[0110] Therefore, the final modulation spectrum is: Z(δ,r,l)=|G g (δ,r,l)|e j∠G(δ,r,l) (32)

[0111] The final modulation spectrum is subjected to an inverse fast Fourier transform to obtain the frequency domain amplitude spectrum after signal enhancement. This frequency domain amplitude spectrum is combined with the frequency domain phase spectrum and then subjected to another inverse fast Fourier transform to obtain the pulse wave signal waveform after spectrum subtraction and phase spectrum compensation enhancement.

[0112] It can be seen from the above technical solutions that the beneficial effects of the present invention are:

[0113] 1. The method proposed in the present invention can collect any superficial artery video of the human body through a camera and analyze a more accurate pulse wave.

[0114] 2. The present invention adopts adaptive gamma transformation in the preprocessing stage, and adjusts the brightness of the original video frame image by adaptively adjusting the size of the γ value to perform image illumination compensation, thereby eliminating the random illumination noise in the high-frequency part.

[0115] 3. The present invention uses the EVM-FRR amplification algorithm to overcome the problem of the traditional Euler amplification algorithm being sensitive to noise. For small motion amplification, the waveform of the original signal will not be distorted due to phase distortion in the extracted one-dimensional time signal sequence.

[0116] 4. The present invention adopts the inter-frame difference method and multiple feature regions to select the region of interest. Among them, the inter-frame difference method is used to improve the operation speed while applying its strong adaptability to dynamic environments, which will assist the illumination compensation in the adaptive gamma transform to obtain a clearer video of the moving area. In addition, by using the multi-feature region selection method, there are overlapping parts in the multiple intercepted regions of interest, which can increase the coverage of the feature region, and to a certain extent, the pulse wave signals to be extracted can be complementary in space, improve the accuracy of feature region tracking, and thus reduce motion artifacts. The two methods are combined to jointly determine the region of interest, which can assist the pre-processing stage in reducing the influence of illumination and eliminate motion artifacts to a certain extent.

[0117] 5. In terms of noise reduction, the present invention adopts a method combining type II Chebyshev filter, wavelet transform of adaptive threshold method and adaptive notch filter to perform multi-dimensional noise reduction on the extraction of pulse wave. Among them, type II Chebyshev filter is mainly used to filter out high-frequency noise, and ensure that the passband signal is not affected by filtering, while the stopband signal is filtered out and the influence of the phase on the original pulse wave waveform is reduced; the wavelet transform of adaptive threshold method filters out Gaussian white noise by adjusting the soft and hard thresholds; the adaptive notch filter is for residual noise interference of some specific frequencies.

[0118] 6. The present invention uses modulation domain spectrum subtraction to enhance weak effective pulse wave signals and improve the extraction of pulse wave signal details in video frames. BRIEF DESCRIPTION OF THE DRAWINGS

[0119] Figure 1 This is a flow chart of the method for accurately extracting pulse waves based on human superficial artery video according to the present invention.

[0120] Figure 2 This is a schematic diagram of a method for accurately extracting pulse waves based on carotid artery video according to an embodiment of the present invention.

[0121] Figure 3 The figure is a flow chart of the inter-frame difference method described in the present invention.

[0122] Figure 4 A schematic diagram of a system for obtaining an enhanced modulation domain by using the modulation domain spectrum subtraction method described in the present invention.

[0123] Figure 5 This is a comparison chart of the color representation of the Adobe RGB color space mentioned in the present invention and the color representation of the standard RGB color space.

[0124] Figure 6 This is a comparison diagram of the carotid artery RGB slice and the Adobe RGB slice according to the embodiment of the present invention. DETAILED DESCRIPTION

[0125] The present invention will be further described below in conjunction with the accompanying drawings.

[0126] Reference Figure 1 The present invention provides a method for accurately extracting pulse waves based on human superficial artery video, taking the carotid artery as an example, comprising the following steps:

[0127] S1: Obtaining a video of the superficial artery of the human body: Use a high-definition camera to capture a 30-second video of the subject's carotid artery, and use natural light sources in the environment to ensure that the shooting area is illuminated as much as possible to ensure the universality and accuracy of the experiment. In order to verify the accuracy of the pulse wave waveform obtained by the method of the present invention, a finger-clip oximeter is used to obtain a pulse wave for comparison while collecting the carotid artery video.

[0128] S2: Preprocess the video using adaptive gamma transform: The collected carotid artery video is used as the input signal, and the guide video is also input. The guided filtering algorithm is used to perform a filtering process that can maintain the edge of the video frame image through the frame-by-frame image of the input signal and the downsampled image of the frame-by-frame image of the guide video. Then, the brightness correction process of the video frame is performed by adaptive gamma transform frame by frame.

[0129] The guided filtering algorithm expression is:

[0130]

[0131] Where D is the frame-by-frame image of the guide video. q is the output image after the frame-by-frame image of the video passes through the filter window, that is, the illumination component of the image. k 、b k is the coefficient of the linear transformation of the filter function. k is the filter window.

[0132] Determine the loss function in the filter window, thereby obtaining the optimal solution of the linear transformation coefficient, so as to achieve the purpose of minimizing the difference between the input image and the guide image. The loss function is:

[0133]

[0134] Among them, X is the original input of the carotid artery video frame by frame, and ε is the regularization parameter. k ,b k ) takes the minimum value, the linear transformation coefficient obtains the optimal solution:

[0135]

[0136]

[0137] Among them, |ω| represents ω kThe number of pixels in the k represents the mean of the guide image D in the filter window; represents the variance of the guide image D in the filter window; is the mean value of the original input image X of the carotid artery video frame by frame in the filtering window.

[0138] The expression of the output frame-by-frame image after guided filtering is obtained as follows:

[0139]

[0140] The illumination value group on the pixel points of the carotid artery video frame image is used as a parameter for adaptive gamma preprocessing, and its expression is:

[0141]

[0142] Among them, O(x, y) is the output image of the carotid artery video frame by frame after gamma transformation, γ is the gamma correction parameter, and the expression of γ is as follows:

[0143]

[0144] Among them, m is the mean brightness of the illumination image; L is the illumination component value at the pixel point.

[0145] S3: EVM-FRR amplification algorithm processes the video frame image after gamma transformation: the traditional Euler amplification method is combined with the FRR filtering method to reduce phase loss while amplifying the signal.

[0146] The EVM-FRR amplification process includes the following steps:

[0147] (1) Spatial decomposition: Convolution and sampling are performed on the rows and columns of each frame of the carotid artery video, and even columns and rows are removed (even columns are removed by convolution and sampling on rows, and even rows are removed by convolution and sampling on columns) to construct a Laplacian pyramid to decompose it into frequency space bands of different scales;

[0148] (2) Time domain filtering: Use a time domain bandpass filter to filter the collected carotid artery video frame image. The time domain bandpass filter formula is as follows:

[0149]

[0150] Among them, I(x,t) is the brightness of a frame image in the video when the spatial point position is x and the time point is t. f(x) is the brightness of the corresponding spatial point position x when time t = 0, that is, f(x) = I(x,0). δ(t) is the motion displacement generated by the spatial point after time t.

[0151] The brightness change signal B(x, t) corresponding to the small signal in the carotid artery video frame image is:

[0152]

[0153] (3) Motion amplification: The brightness change signal B(x, t) corresponding to the tiny signal in the carotid artery video frame image is amplified by α times and then superimposed on the corresponding image brightness signal I(x, t), as shown in the following formula:

[0154]

[0155] Where α is the amplification factor, It is the image brightness signal at the spatial point x and the time point t after the video frame image signal is amplified and superimposed once.

[0156] (4) Time reversal: Time-reversing at the spatial point x obtains a reversed time signal at time -t;

[0157] (5) Secondary motion amplification: With the inverted time signal as input, the image brightness change signal corresponding to the carotid artery pulsation signal at the spatial position x and time -t is amplified again by α times through the amplification factor relationship, and superimposed on the time-reversed signal sequence of step (4);

[0158] (6) Secondary time reversal and output: The carotid artery pulsation signal after secondary amplification is time-reversed again to obtain an output video frame signal with a spatial position of x and a time point of t. That is the carotid artery video frame signal after EVM-FRR amplification.

[0159] S4: Inter-frame difference method and multi-feature region selection, the carotid artery video amplified by EVM-FRR is used as the input of the signal selection region of interest module. Determine the grayscale level L in each frame of the carotid artery video image, give an initial threshold T based on experience in advance, and partition each video frame image with the initial threshold T as the boundary. The part with grayscale value less than T is region K0, and the part with grayscale value greater than T is region K1.

[0160] The probability P and average grayscale E of the pixel points in the area K0 and area K1 in the entire video frame image are calculated respectively, and the expressions are as follows:

[0161] Area K0:

[0162]

[0163] Area K1:

[0164]

[0165] The threshold T corresponding to the maximum inter-class variance is selected as the optimal binarization threshold, where the inter-class variance formula is

[0166] For: S 2 (T) = P0P1(E0-E1) 2 (47)

[0167] The pixel grayscale values ​​of the carotid artery video frame are compared with the optimal binarization threshold one by one. The pixel values ​​of points greater than the optimal threshold are changed to 255, which are determined as moving parts; the pixel values ​​of points less than the optimal threshold are changed to 0, which are determined as static parts. In this way, the carotid artery pulsation effective area and other non-pulsation areas are divided.

[0168] The SIFT algorithm is used in the carotid artery pulsation area to mark the feature points of the carotid artery, and among the marked feature points, three feature points are grouped into one group, and 6 feature point groups are selected.

[0169] Determine the centroid coordinates of each feature point group, and cut a rectangular area with the centroid as the center as a single region of interest of the pulse wave signal in the video frame image. Complete the six regions of interest in the same way. The centroid calculation formula is as follows:

[0170]

[0171]

[0172] Among them, x and y are the horizontal and vertical coordinates of the center of mass, respectively. i ,y i are the horizontal and vertical coordinates of each feature point (i=1, 2, 3).

[0173] S5: Extract the original pulse wave signal. Slice the region of interest of each video frame image and convert the slice image from RGB color space to Adobe RGB color space. The conversion formula is as follows:

[0174]

[0175] For the slices of the region of interest in the Adobe RGB color space, the RGB three channels are separated, and each video frame image obtains three separated single-channel grayscale images. Each video frame image in each channel is averaged, and each video frame single-channel image can obtain a pixel average value. The obtained pixel average values ​​are combined in the time sequence of the video frames to obtain a one-dimensional time series of a single channel. Similarly, the one-dimensional time series of the three channels are obtained, represented by x1(t), x2(t) and x3(t).

[0176] For a one-dimensional time series xi (t)(i=1,2,3) is normalized, and its expression is as follows:

[0177]

[0178] Among them, μ i For x i The pixel mean of (t), σ i For x i (t) is the pixel standard deviation.

[0179] After the normalized signal is transformed into a frequency domain signal by fast Fourier transform, the frequency domain signal is initially filtered in the range of 0.8 to 3.3 Hz using a bandpass filter, and the relatively pure waveform channel waveform is selected as the original pulse wave.

[0180] S6: Signal denoising is accomplished by combining type II Chebyshev filter, wavelet transform with adaptive threshold method, and adaptive notch filter. The extracted original pulse wave signal is used as the input of the signal denoising module, and the attenuation of the filter function ∈ is determined according to different decibels χ. The formula is as follows:

[0181]

[0182] The formula of the original pulse wave signal through the type II Chebyshev filter is:

[0183]

[0184] Where W0 is the cut-off frequency, is a Chebyshev polynomial of order n.

[0185] The signal waveform filtered by the II-type Chebyshev filter is used as input for wavelet decomposition to determine the wavelet coefficients of the signal decomposition. The wavelet coefficients are then subjected to adaptive threshold processing through a threshold function, where the formula for the adaptive threshold function of the wavelet transform is as follows:

[0186]

[0187] in, is the wavelet decomposition coefficient after processing; m j,k is the wavelet decomposition coefficient before processing; and n are both positive numbers, when When is 0 and n is 0, the wavelet transform is the wavelet transform of the soft threshold function. When set to a sufficiently large value, the wavelet transform is transformed into a wavelet transform of a hard threshold function. The formula for determining the λ value is as follows:

[0188]

[0189] Among them, MAD is the median value of the absolute value of the first layer wavelet, and N is the length of the video frame signal.

[0190] The wavelet coefficients determined by the adaptive threshold wavelet transform function are used to perform wavelet reconstruction on the pulse wave signal to obtain a reconstructed pulse wave, thereby removing the Gaussian white noise in the original pulse wave signal.

[0191] The frequency distribution exceeding the amplitude threshold in the frequency domain segment of 0.8 to 3.3 Hz is found in the reconstructed pulse wave signal, and the number of notch points is determined by the frequency distribution. Each notch point corresponds to two mutually orthogonal input parameter signals.

[0192]

[0193] Where ω0 is the center frequency and A is the filter amplitude. Different filter frequencies ω0 are configured by frequency distribution exceeding the threshold, and the bandwidth is adjusted by configuring the filter amplitude A.

[0194] The output of the filtered signal is: e(n) = d(n) - y(n) (57)

[0195] Where d(n) is the input signal and y(n) is the output of the adaptive filter. The formula for y(n) is:

[0196] y(n)=ω1(n)x1(n)+ω2(n)x2(n)+…+ω 2k-1 (n)x 2k-1 (n)+ω 2k (n)x 2k (n) (58)

[0197] The weight coefficient ω of the reference signal i (n) is determined by the following formula:

[0198] ω i (n+1)=ω i (n)+2μe(n)x i (n) (59)

[0199] Among them, μ is the length of the adjustment step, and the weight coefficient of the reference signal is adjusted through e(n) feedback.

[0200] S7: Perform fast Fourier transform on the filtered signal to obtain the modulation domain amplitude spectrum |X(δ,r,l)| and modulation domain phase spectrum ∠X(δ,r,l) of the signal, where r is the discrete frequency, δ is the modulation frame, and l is the modulation frequency.

[0201] Perform modulation domain spectrum subtraction on the modulation amplitude spectrum of the pulse wave signal to obtain the improved modulation domain amplitude value:

[0202]

[0203] Among them, β is used to determine the type of spectrum subtraction. When β = 1, it is determined to be the amplitude spectrum subtraction type; when β = 2, it is determined to be the energy spectrum subtraction type. η is the over-subtraction factor in the modulation domain, ξ is the spectrum suppression factor, that is, the lower limit of the enhanced modulation spectrum, between 0 and 1; |E(δ,r,l)| is the modulation amplitude spectrum estimate of the noise.

[0204]

[0205] Wherein, N is the number of frames corresponding to the no-interest region (noise region) of the leading video frame.

[0206] The modulation domain amplitude spectrum improved by modulation domain spectrum subtraction and the unimproved modulation domain phase spectrum ∠X(δ,r,l) together constitute the modulation domain amplitude spectrum estimate G after modulation domain spectrum subtraction. g (δ,r,l):

[0207] G g (δ,r,l)=|G(δ,r,l)|e j∠X(δ,r,l) (62)

[0208] The phase compensation degree Γ(δ, r, l) is determined by estimating the modulation amplitude spectrum of the noise:

[0209] Γ(δ,r,l)=λφ(l)|E(δ,r,l)| (63)

[0210] Among them, λ is an empirical constant, φ(l) is an antisymmetric function, and its formula is as follows:

[0211]

[0212] Modulation domain amplitude estimation G after spectral subtraction g (δ, r, l) and the phase compensation degree Γ(δ, r, l) jointly synthesize the modulation spectrum Z′(δ, r, l) after modulation domain spectrum subtraction enhancement:

[0213] Z′(δ,r,l)=G g (δ,r,l)+Γ(δ,r,l) (65)

[0214] Therefore, the modulation phase spectrum after modulation spectrum subtraction enhancement is:

[0215] ∠G(δ,r,l)=arg[Z′(δ,r,l)] (66)

[0216] The final modulation spectrum is obtained: Z(δ,r,l)=|G g (δ,r,l)|e j∠G(δ,r,l) (67)

[0217] The final modulation spectrum is subjected to an inverse fast Fourier transform to obtain the frequency domain amplitude spectrum after signal enhancement. After combining the frequency domain amplitude spectrum with the frequency domain phase spectrum, an inverse fast Fourier transform is performed again to obtain the pulse wave signal waveform after spectrum subtraction and phase spectrum compensation enhancement.

[0218] The above embodiments are only used to illustrate the technical solutions of the present invention, but not to limit the same. Although the above embodiments describe the present invention in detail, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or replace some of the technical features therein with equivalents. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for accurately extracting pulse waves based on human superficial artery video, characterized in that: It includes the steps: Step S1: obtaining a video of a superficial artery of a human body; Step S2: pre-processing the video using adaptive gamma transform; Step S3: using the EVM-FRR amplification algorithm to amplify the motion of the human superficial artery video; Step S4: Adopting the adaptive threshold frame difference method and selecting multiple feature regions to determine the region of interest (ROI); Step S5: extracting the original pulse wave signal; Step S6: using a type II Chebyshev filter, a wavelet transform using an adaptive threshold method, and an adaptive notch filter to filter out possible noise in the signal; Step S7: enhancing the filtered pulse wave signal using modulation domain spectrum subtraction; The EVM-FRR amplification algorithm of step S3 is an algorithm combining the Euler video amplification algorithm and the FRR filter, which includes the following steps: A: Spatial decomposition: Construct a Laplacian pyramid to spatially sample the video frames and decompose the spatial domain into frequency-space bands of different scales; B: Time domain filtering: Use a time domain bandpass filter to perform time domain filtering on the frequency space band after Laplace decomposition to obtain a valid signal in the frequency range of 0.8 to 3.3 Hz; C: Signal amplification: Through the amplification factor relationship, the acquired effective signal is amplified by α times; D: Amplification and reconstruction: the amplified video frame is expanded to the size of the original video and superimposed with the effective signal at the spatial point x and the time point t in the original video channel; E: Time reversal: reverse the time of the superimposed effective signal and pass it through the filter in the reverse direction; F: Amplification enhancement: Repeat the spatial decomposition, time domain filtering, signal amplification, and amplification reconstruction in steps AD for the signal that passes through the filter in the reverse direction to obtain the effect of double amplification of the video frame; G: Signal output: For the amplified and enhanced video frame signal at the spatial point x, it is output after time reversal; Among them, the adaptive threshold frame difference method of step S4 includes the following steps: subtracting the grayscale value of the image pixel point of the i-1 frame from the i-th frame in the sampled video, and calculating the absolute value of the image pixel difference to obtain the video frame image after the difference; binarizing the pixel points of the video frame image; comparing the binarization result with the set threshold; determining the pixel area with a grayscale value greater than the threshold as a moving part, otherwise it is a static area; after processing all the video frames, each frame is connected to form a new video through connectivity analysis, that is, a moving part video is obtained; wherein i=2,3,4... Among them, the multi-feature area selection method in step S4 includes the following steps: edge detection is performed on the enlarged motion area to obtain feature points and marked coordinates of the feature points; three feature points are selected from the marked feature points to form a feature point group, and the area centroid of the feature point group is used as the center point of the feature area. By means of area positioning, the coordinates of the center position of the area surrounded by the feature point group are calculated, and a rectangular area is intercepted as a region of interest (ROI), thereby obtaining a single ROI area; multiple ROI areas are selected by the same method; multiple ROI areas are obtained by using multiple feature area selection methods, so that the coverage of the feature area is increased, and the extracted multiple blood volume pulse wave signals can be complementary in space, so as to remove motion interference, further remove system errors, and better restore the pulse wave signal.

2. The method for accurately extracting pulse waves based on human superficial artery video as claimed in claim 1, characterized in that: Step S2 is an adaptive gamma transform for the video preprocessing system. The main steps are to first perform fast guided filtering on the downsampled images of the input image and the guide image to obtain the illumination component at each pixel in the video frame, and use the pixel value of each pixel (x, y) in the video frame as input. The illumination value at the pixel point and the illumination mean of the image are used as two input parameters of the point, and the variables in the gamma function are adjusted to eliminate the influence of uneven illumination on the image quality.

3. The method for accurately extracting pulse waves based on human superficial artery video as claimed in claim 1, characterized in that: Step S5 extracts the original pulse wave signal and mainly includes: firstly, slicing the region of interest in the video frame image frame by frame, and converting the sliced ​​image from the RGB color space to the Adobe RGB color space; then separating the three channels of R, G, and B in the Adobe RGB color space; performing pixel mean processing on each separated channel to obtain a one-dimensional time series of the three channels; then performing normalization processing, Fourier transform processing, and bandpass filtering on the one-dimensional time series of each channel to select a channel with a higher signal-to-noise ratio to obtain the original pulse wave signal.

4. The method for accurately extracting pulse waves based on human superficial artery video as claimed in claim 1, characterized in that: Step S6 uses a type II Chebyshev filter to perform filtering using the first formula as a transfer function, and determines the attenuation of different decibels using the second formula; wherein the first formula is: Where W0 is the cut-off frequency, is the nth order Chebyshev polynomial; the second formula is: Among them, ∈ is the attenuation degree and χ is the decibel number.

5. The method for accurately extracting pulse waves based on human superficial artery video according to claim 1, characterized in that: Step S6 The wavelet transform of the adaptive threshold method is mainly achieved by changing the variable parameters in the third formula and n to switch between soft threshold function and hard threshold function; when When is 0 and n is 0, the wavelet transform is the wavelet transform of the soft threshold function. When it is set to a sufficiently large value, the wavelet transform is transformed into a wavelet transform of a hard threshold function; the third formula is: in, is the wavelet decomposition coefficient after processing; m j,k is the wavelet decomposition coefficient before processing; the λ value is determined by the fourth formula, and n are both positive numbers; wherein the fourth formula is: Among them, MAD is the median value of the absolute value of the first layer wavelet, and N is the length of the video frame signal.

6. The method for accurately extracting pulse waves based on human superficial artery video according to claim 1, characterized in that: Step S6 adaptive notch filtering limits the waveform after wavelet transformation to a frequency domain segment of 0.8 to 3.3 Hz, and records the corresponding frequency distribution of the part exceeding the amplitude threshold within this frequency range to obtain the number of notch points exceeding the amplitude threshold and the filtering center frequency; the orthogonal signal of the center frequency is used as a reference signal, and then the weight coefficient of the linear combination is continuously adjusted through the residual generated corresponding to each frequency point of the video signal, and the linear input signal of the input video frame image is separated from the noise signal part according to the weight coefficient.

7. The method for accurately extracting pulse waves based on human superficial artery video according to claim 1, characterized in that: Step S7 modulation domain spectrum subtraction suppresses residual noise while enhancing the effective signal of the pulse wave; it includes the following steps: firstly, obtaining the amplitude spectrum and phase spectrum from the filtered signal, and performing fast Fourier transform on the amplitude spectrum at each frequency point along the time to obtain the modulation amplitude spectrum and the modulation phase spectrum; The phase spectrum is compensated by using an antisymmetric function to change the angular relationship between the conjugate terms of the amplitude spectrum; The modulation phase spectrum and the modulation amplitude spectrum after phase compensation together constitute the modulation spectrum, and then the frequency domain amplitude spectrum is obtained through inverse Fourier transform. This frequency domain amplitude spectrum is combined with the frequency domain phase spectrum; the pulse wave signal in the video frame image can be enhanced by performing inverse Fourier transform again.

Images