TWI653027B - Algorithmic method for extracting human pulse rate from compressed video data of a human face - Google Patents

Abstract

The present invention includes a first overlapping cutting step, a first processing step, a first overlapping addition step, a second overlapping cutting step, a second processing step, and a second overlapping addition step. By the foregoing steps, a method of separating single channel signals is performed to perform physiological signal extraction. Then, the green channel with the least influence of the compression algorithm is selected for processing, and then processed by singular spectrum analysis, double relationship filtering, frequency mask screening, etc., and finally the processed final heart rate signal is obtained. Therefore, the heart rate extraction algorithm that has the face compression image is quite novel, the application range is wide, and the single channel signal separation method greatly reduces the video image transmission amount.

Description

Heart rate extraction algorithm for face compressed image

The invention relates to a heart rate extraction calculation method for a face compressed image, in particular to a heart rate extraction algorithm with a face compressed image which is quite novel, has a wide application range and a single channel signal separation method substantially reduces the amount of video image transmission. Heart rate extraction calculation method for compressed images.

In 2008, Wim Verkruysse and others pointed out that heart rate waveforms can be detected by natural light and consumer cameras, so that physiological information can be analyzed remotely.

In 2010, Ming-Zher Poh et al. used blind source separation technology to analyze the color signals in the video and extract the heart rate signals from the human body.

In 2012, Hao-Yu Wu and others proposed the Euler image enlargement algorithm to amplify the small changes in skin color in the video stream, so that changes in skin color can be observed.

In 2013, Guha Balakrishnan and others pointed out that the rhythmic beat of the heart can cause small shaking of the head, which can be detected by detecting the slight shaking of the head in the video.

In 2015, the City University of Hong Kong and the Eindhoven University of Technology in the Netherlands each developed a measurement technique that restricted the subject's limbs in a non-sports state, such as the subject being able to turn or shake the head. In the case of the case, a reliable heart rate value is detected.

In 2016, Tulyakov et al. proposed an adaptive matrix filling method to automatically detect which region contains heart rate signals while detecting heart rate signals, and use these regions for heart rate detection.

All of the above methods use uncompressed video data as the processing object, and the three-channel signal separation technology is used to separate the heart rate signal from the original signal. When processing compressed video, because the compression algorithm has a great influence on the heart rate signal, the above methods cannot effectively obtain accurate heart rate information.

In 2016, Hanfland and others compressed the original video and compared the compressed video with the original video. The results show that the heart rate signal still exists in the compressed video, but the overall quality has been greatly reduced.

In 2017, McDuff et al. used two compression methods (x264 and x265) to compress the original video to different bit rates. The results show that video compression significantly reduces the signal-to-noise ratio of the heart rate signal.

In addition, the current remote light volume change tracing method (Rppg) uses uncompressed image as the processing object, that is, extracts the heart rate information implied by the face on the uncompressed image data; one of the negative effects of this method It will cause storage difficulties. As we all know, the amount of data in uncompressed images is huge. For example, a resolution of 640*480, 30fps one minute image requires about 1.7GB of storage space. Such a huge storage requirement will inevitably result in waste of resources; another problem is that uncompressed image data cannot be transmitted at a long distance, which brings great limitations to the application of rPPG technology. At present, the network transmission capability cannot realize the instant transmission of uncompressed images; the existing rPPG technology cannot be applied to the occasions requiring real-time transmission of long-distance images.

In short, today's images are transmitted over long distances, and image compression is almost always performed, reducing the amount of data transmitted. The more common compression methods are the following four types: x264, x265, vp8, and vp9. However, the compressed image is almost ineffective in obtaining accurate heart rate information.

In view of this, it is necessary to develop a technique that can solve the above disadvantages.

The object of the present invention is to provide a heart rate extraction calculation method for a face compressed image, which has a novel heart rate extraction algorithm for a face compressed image, a wide application range and a single channel signal separation method. The width reduces the amount of video image transmission and so on. In particular, the problem to be solved by the present invention is that uncompressed image data cannot be transmitted over long distances.

The technical means for solving the above problem is to provide a heart rate extraction calculation method for a face compressed image, which comprises the following steps: 1. The first overlapping cutting step; two. The first processing step: [a] pre-processing step; [b] band-pass filtering processing step; [c] first singular spectrum analysis step; [d] double-relationship screening step; [e] reconstruction step; The first overlap addition step; four. The second overlapping cutting step; The second processing step: [f] frequency mask construction step; [g] second singular spectrum analysis step; [h] frequency mask screening step; The second overlap addition step.

The invention will be described in detail in the following examples in conjunction with the drawings:

10‧‧‧Video installation

20‧‧‧Network connection device

S1‧‧‧ first overlapping cutting step

S2‧‧‧ first processing steps

S21‧‧‧Pretreatment steps

S22‧‧‧ Bandpass Filtering Procedure

S23‧‧‧First singular spectrum analysis step

S24‧‧‧Two-fold relationship screening step

S25‧‧‧Reconstruction steps

S3‧‧‧ first overlap addition step

S4‧‧‧Second overlapping cutting step

S5‧‧‧Second processing steps

S51‧‧‧ Frequency mask construction steps

S52‧‧‧Second singular spectrum analysis step

S53‧‧‧frequency mask screening step

S6‧‧‧Second overlap addition step

M‧‧‧Face compressed video

M1‧‧‧Small video

T1, T2‧‧‧ length of time

T3‧‧ ‧ long time

K‧‧‧ original image

P‧‧‧Face area

G0(t)‧‧‧ original signal

G1(t)‧‧‧ first signal

G1m‧‧‧ first subsequence

G1f‧‧‧ first main frequency

G2(t)‧‧‧second signal

G22‧‧‧second signal after overlapping and adding

G3(t)‧‧‧third signal

G3j‧‧‧ center frequency

G3ju‧‧‧pass upper limit frequency

G3jd‧‧‧pass lower limit frequency

G3m‧‧‧Second subsequence

G3f‧‧‧second main frequency

S34‧‧‧frequency mask screening step

G4(t)‧‧‧4th signal

G44‧‧‧4th signal after overlapping and adding

L1, LA‧‧‧ first curve

L2, LB‧‧‧ second curve

L3, LC‧‧‧ third curve

L4, LD‧‧‧ fourth curve

Figure 1 is a flow chart of the calculation method of the present invention

Figure 2 is a schematic view of the present invention

Figure 3 is a schematic view of the first overlapping cutting process of the present invention.

Figure 4 is a schematic diagram of each of the short films of the present invention containing N祯 original images.

5A and 5B are respectively schematic diagrams of tracking of a face region according to the present invention.

Figure 6 is a schematic diagram of the original signal of the present invention.

7 and 8 are respectively schematic diagrams of the pre-processing and the post-processing of the band-pass filtering process of the original signal of the present invention.

Figure 9 is a schematic diagram of the first signal of the present invention.

Figure 10 is a schematic diagram of the first singular spectrum analysis step of the present invention.

Figure 11 is a schematic view of the second signal of the present invention

Figure 12 is a schematic diagram of the first overlap addition process of the present invention.

Figure 13 is a schematic view of the second overlapping cutting process of the present invention

Figure 14 is a schematic diagram of the construction process of the frequency mask of the present invention

Figure 15 is a schematic diagram of the second singular spectrum analysis process of the present invention.

Figure 16 is a schematic diagram of the second overlap addition process of the present invention.

Figure 17A is a schematic diagram of a static film of the present invention after performing a band pass filtering process

Figure 17B is a schematic diagram of the first singular spectrum analysis step of Figure 17A and the double relationship screening processing step

Figure 17C is a schematic diagram of the second singular spectrum analysis and frequency mask screening step of Figure 17B

Figure 17D is a comparison chart of 17A, 17B and 17C

Figure 18A is a schematic diagram of a dynamic film of the present invention after performing a band pass filtering process

Figure 18B is a schematic diagram of the first singular spectrum analysis step of Figure 18A and the double relationship screening processing step

Figure 18C is a schematic diagram of the second singular spectrum analysis and frequency mask screening step of Figure 18B

Figure 18D is a comparison diagram of the extracted final heart rate signal of the present invention after the 18A, 18B and 18C pictures.

Referring to FIGS. 1 and 2, the present invention is a heart rate extraction calculation method for a face compressed image, which includes the following steps:

One. The first overlapping cutting step S1: referring to FIG. 3, obtaining a face compressed movie M after video compression processing, the length of time is T1, and the face compressed movie M is superimposed and cut into a plurality of small segments of film M1. The length of each small piece of film M1 is T2; the method of overlapping cutting is defined as the beginning of the face compression film M, taking the first piece of the film M1, and then every other step for a long time T3 Then, another small piece of film M1 is retrieved, and the foregoing action is repeated until the face compression movie M ends, so that a plurality of the small pieces of the movie M1 can be obtained. For example, the face compression movie M is a movie of 600 seconds (that is, the length of time T1), and each of the small films M1 is 3 seconds (that is, the length of time T2), and the step T3 is 1.5 seconds for a long time. It is overlapped and cut into 399 small pieces of film M1.

two. In the first processing step S2, the following steps are performed one by one for each of the small segments of the movie M1:

[a] Pre-processing step S21: Each of the small-length films M1 includes N祯 original images K (see FIG. 4); face region tracking is performed for each of the original images K (see FIGS. 5A and 5B, In the prior art, a face region P having X by Y pixels is averaged, and the green values of the complex pixels of the X by Y are averaged to obtain a scalar quantity, which is defined as a single祯 average green value; and so on, the average green value of N祯 can be obtained, and an original signal G0(t) is further obtained from the average green value of N祯 (refer to Fig. 6, the longitudinal coordinate axis GV is the average green value, and GV is the abbreviation of English Green Value, and G0(t) is a signal curve formed by the average green value of N single ticks over a period of time), where t=1 to N.

[b] band pass filtering processing step S22: performing band pass filtering on the original signal G0(t), and selecting a signal between frequencies of 0.8 Hz to 2.0 Hz (see FIGS. 7 and 8) to obtain a first A signal G1(t) (see Figure 9), where t = 1 to N.

[c] first singular spectrum analysis (SSA) step S23: decomposing the first signal G1(t) into a plurality of first sub-sequences G1m (see FIG. 10), and for all the One subsequence The column G1m performs fast Fourier transform to obtain a spectrum of the plurality of first sub-sequences G1m, and the frequency having the largest amplitude in the spectrum of each of the first sub-sequences G1m is taken as its first main frequency G1f. The foregoing singular spectrum analysis is a known technology. In practice, the application software that is generally distributed in the market (for example, MATLAB) can be implemented, so the detailed process of the singular spectrum analysis will not be described here.

[d] double relationship screening step S24: comparing the first sub-sequence G1m two or two, if the two first main frequencies G1f are doubled after the comparison, the two first main frequencies G1f are retained, otherwise The two first main frequencies G1f are discarded. If none of the two or two alignments satisfy the double relationship, then all of the first subsequences G1m are retained; and at least two of the first reserved subsequences are obtained.

[e] reconstruction step S25: reconstructing all of the first reserved subsequences of the previous step to obtain a second signal G2(t) (refer to FIG. 11), where t=1 to N, the second signal The horizontal axis of G2(t) is time, the vertical axis is the intensity of green value, and the second signal G2(t) has the length of time T2.

three. First overlap addition step S3: Referring to FIG. 12, the plurality of second signals G2(t) obtained after the first processing step S2 are processed by the overlap-add technique of the existing raised cosine window method to obtain one The overlapped second signal G22 has the length of time T1; the overlap and add technique is defined as being adjacent between the two second signals G2(t), and the overlap is added by the long time T3. That is, the above processing is performed for each small segment signal (for example, the second signal G2(t)), the processing result is multiplied by the raised cosine window (Hanning window), and the processing results are added to obtain the processed long time. The signal (for example, the second signal G22 after the overlap is added). In practice, the application software that can be circulated through the general market (for example, MATLAB) can be implemented accordingly, so the detailed process will not be described here.

four. The second overlap cutting step S4: refer to FIG. 13 to obtain the overlapped second signal G22 having the time length T1, and the second signal G22 is overlapped and cut into a plurality of small segments after the overlap is added. Each of the small segments of the signal has the length of time T2; the manner of the overlapping cut is defined as the beginning of the second signal G22 after the overlap is added, and a small segment of the signal is extracted, and then every other step is T3 for a long time. Then take another small segment of the signal and repeat the above actions until the overlap is added. After the second signal G22 ends, a plurality of the small signals can be obtained, and each of the small signals is defined as a third signal G3(t), where t=1 to N.

Fives. The second processing step S5: performing the following steps one by one for each of the third signals G3(t):

[f] Frequency mask construction step S51: the frequency at which the amplitude of the third signal G3(t) is the largest is taken as the center frequency G3j (refer to FIG. 14), and the upper limit is set with the center frequency G3j as the center. The frequency G3ju and one pass the lower limit frequency G3jd to obtain a frequency mask range.

[g] second singular spectrum analysis (SSA) step S52: decomposing the third signal G3(t) into a plurality of second sub-sequences G3m (see Fig. 15), and performing fast on all second sub-sequences G3m The Fourier transform converts the spectrum of the plurality of second sub-sequences G3m, and the frequency having the largest amplitude in the spectrum of each of the second sub-sequences G3m is taken as its second main frequency G3f.

[h] frequency mask screening step S53: only retaining the second sub-sequence G3m of the second main frequency G3f within the frequency mask range, then screening at least one second reserved sub-sequence, which becomes one The fourth signal G4(t), where t=1 to N, the horizontal axis of the fourth signal G4(t) is time, and the vertical axis is the green value intensity. If the screening does not obtain any of the second reserved subsequences, the fourth signal G4(t) is directly equal to the third signal G3(t).

six. Second overlap addition step S6: Referring to FIG. 16, the plurality of fourth signals G4(t) obtained after the second processing step S5 are processed by the overlap-add technique of the existing raised cosine window method. The fourth signal G44 is obtained after the overlap and addition, and has the length of time T1, which is the final heart rate signal.

In practice, in the original image capturing step S1, two sets of video devices 10 and one network connecting device 20 are provided. The two sets of video devices 10 are connected to the network connection device 20 to achieve a video contact.

The compression method of the face compressed image is selected from one of x264, x265, vp8, and vp9.

Regarding the aforementioned singular spectrum analysis (SSA), which is a well-known technique, it is explained as follows: Input: a vector of length N = [ x ₁ , x ₂ , ..., x _N ] ^T .

The first step: Hankel matrix embedding (English is Hankel matrix embedding), which is to convert y into matrix X as follows.

Where K = N - L +1, L and K are the rows and columns of the matrix, which are set by the user.

The second step: the aforementioned matrix X performs Singular value decomposition (SVD):

Where σ _i is the singular value of X, u _i and v _i are the corresponding singular vectors, and r is the rank of X. make For a submatrix of rank 1, then: X = X ₁ + X ₂ +... + X _r .

The third step: reconstruction. The method is Diagonal averaging, which is more complicated to refer to. The function is to convert the matrix addition form of the above formula into a vector addition form: x = y ₁ + y ₂ +... + y _r .

Each y _{i is} called a reconstructed component (or RC for short).

Step 4: Select the RC that meets the needs in y _i as needed to get the final time series.

For example: In our algorithm, x is the time series obtained by a time window of the green channel. For example, if the time window (ie, the length of time T2) is 3s and the frame rate is 30fps, the length of x is N=90. K is generally set to half of N, which is 45. Calculated K = 46. After SVD, r is generally the smallest of K and N, ie r=45. In the fourth step, there are r=45 y _i , we can only consider the first 20 to speed up the calculation.

The focus of the invention is to fully consider the influence of the compressed image on the physiological signal of the human body, and adopt the method of single channel signal separation to perform physiological signal extraction. Then select the channel (G channel) that is least affected by the compression algorithm for processing. That is, only the green value of the X-by-th pixel of the face region is averaged. The advantage of this is to fully avoid the impact of the compression algorithm on the heart rate signal, so that the heart rate signal can be accurately and stably extracted from the compressed video. In addition, with the singular spectrum analysis (SSA), the frequency structure of the heart rate signal can be used to effectively identify the heart rate signal in the mixed signal containing the noise. Then through the frequency mask, the noise is further filtered to obtain a more accurate heart rate signal. The invention extracts the heart rate signal from the mixed signal by three methods, can greatly ensure the accuracy of the heart rate signal, and can effectively avoid the influence of the compression algorithm on the heart rate signal calculation.

For the experimental results of this case, please refer to Figures 17A to 18D.

In a static movie (the compression mode is vp8, the bit rate is 100 kb/3), when only the band pass filtering process (referred to as filter) of the third step is performed, the waveform of the result in the frequency domain is as shown in FIG. 17A; If the first singular spectrum analysis and the double relationship filtering process (referred to as filter+SSA) are continued, the waveform of the result in the frequency domain is as shown in Fig. 17B; if it continues to the second singular spectrum analysis and The frequency mask screening step (referred to as filter + SSA + refine; that is, the case), the result of the waveform in the frequency domain is shown in Figure 17C. The comparison between the three and the actual heart rate signal in the time domain is as shown in Fig. 17D (where the first curve L1, the second curve L2, the third curve L3, and the fourth curve L4 respectively represent the actual heart rate, band pass filtering, band As shown by the pass filter + SSA, band pass filter + SSA + frequency mask, it can be seen that the result of this case is closest to the actual actual heart rate signal.

In addition, in a dynamic movie (the compression mode is x264, the bit rate is 584 kb/3), when only the band pass filtering process (referred to as filter) in the third step is performed, the waveform in the frequency domain is as shown in FIG. 18A. If the first singular spectrum analysis and the double relationship filtering process (referred to as filter+SSA) are continued, the waveform of the result in the frequency domain is as shown in Fig. 18B; if it continues to the second singular spectrum Analysis and frequency mask screen The selection step (referred to as filter + SSA + refine; that is, the case), the result of the waveform in the frequency domain is shown in Figure 18C. The comparison between the three and the actual heart rate signal in the time domain is as shown in Fig. 18D (where the first curve LA, the second curve LB, the third curve LC, and the fourth curve LD respectively represent the actual heart rate, band pass filtering, band Through filtering + SSA, band pass filtering + SSA + frequency mask), it can be seen that the result of this case is closest to the actual actual heart rate signal.

The advantages and functions of the present invention are as follows:

[1] The heart rate extraction algorithm for face compressed images is quite novel. The invention adopts unique singular spectrum analysis, double relationship screening, frequency mask screening and other processing steps, and can extract the human heart rate signal in the compressed video data, which is an unprecedented technical means. Therefore, the heart rate extraction algorithm for the face compressed image of this case is quite novel.

[2] A wide range of applications. The invention can be applied to all occasions where image compression is required. For example, in telemedicine, the patient's video data is compressed and transmitted to the hospital for further analysis. In mobile applications, video captured by users is transmitted over the wireless network to the cloud for heart rate measurement and analysis. By the invention, the technology can be applied to the fields of telemedicine, home care or fitness training, especially in the field of home care in the field of telemedicine. Therefore, the application range is wide.

[3] The single channel signal separation method greatly reduces the amount of video image transmission. The invention adopts the method of single channel signal separation to perform physiological signal extraction, and then selects the channel (G channel) which is least affected by the compression algorithm for processing. That is, only the green value of the plurality of pixels of X by Y in the face area is averaged, and the video image transmission amount can be greatly reduced. Therefore, the single channel signal separation method greatly reduces the amount of video image transmission.

The present invention has been described in detail with reference to the preferred embodiments of the present invention, without departing from the spirit and scope of the invention.

Claims (3)

A heart rate extraction calculation method for a face compressed image includes the following steps: 1. The first overlapping cutting step: obtaining a face compressed movie subjected to video compression processing, and superimposing the face compressed film into a plurality of short films; the overlapping cutting manner is defined as compressing the movie by the face At the beginning, the first piece of the film is captured, and then every other step is taken for another long time, and then the other piece of the film is captured, and the action is repeated until the face compression movie ends, and then the plurality of pieces of the film can be obtained; two. In the first processing step, the following steps are performed one by one for each of the short films: [a] Preprocessing step: each of the short films contains N祯 original images; face region tracking is performed for each original image. Obtaining a face region having X by Y pixels, averaging the green values of the plurality of pixels of the X by Y, to obtain a scalar quantity, which is defined as a single 祯 average green value; and so on , an N祯 average green value is obtained, and an original signal G0(t) is obtained from the N祯 average green value, wherein t=1 to N; [b] band pass filtering processing step: the original signal G0 (for the foregoing) t) performing band pass filtering processing, selecting a signal between frequencies of 0.8 Hz to 2.0 Hz, and obtaining a first signal G1(t), where t=1 to N; [c] first singular spectrum analysis step: Decoding the first signal G1(t) into a plurality of first sub-sequences, and performing fast Fourier transform on all the first sub-sequences to obtain a plurality of spectra of the first sub-sequence, and then each of the first sub-sequences The frequency with the largest amplitude in the spectrum as its first dominant frequency; [d] Relationship screening step: comparing the first sub-sequences to two or two, if the two first main frequencies are doubled after the comparison, the two first main frequencies are reserved, otherwise the two first main frequencies are discarded; If none of the two or two comparisons satisfy the double relationship, then all the first subsequences are retained; and at least two first reserved subsequences are obtained; [e] reconstruction step: all of the previous steps The first reserved subsequence is reconstructed to obtain a second signal G2(t), where t=1 to N, the horizontal axis of the second signal G2(t) is time, the vertical axis is green value intensity, and the second The signal G2(t) has the length of time; The first overlapping addition step: the plurality of second signals G2(t) obtained after the first processing step are processed by the overlap addition technique of the existing raised cosine window method to obtain an overlapped second signal , having the length of time; the overlap-adding technique is defined as being adjacent between the two second signals G2(t), and the steps are overlapped and added for a long time; The second overlapping cutting step is: obtaining the second overlapping signal after the overlap, and having the length of time, after the overlapping is added, the second signal is overlapped and cut into a plurality of small segments, and each of the small segments has the length of time The method of overlapping cutting is defined as the beginning of the second signal after the overlap is added, taking a small segment of the signal, and then every other step for a long time, and then capturing another small segment of the signal, repeating the foregoing After the overlapping, the second signal ends, and then the plurality of small signals can be obtained, and each of the small signals is defined as a third signal G3(t), where t=1 to N; The second processing step: performing the following steps one by one for each of the third signals G3(t): [f] Frequency mask construction step: the frequency with the largest amplitude of the third signal G3(t) is taken as its center frequency And centering on the center frequency, setting a pass frequency and a pass frequency to obtain a frequency mask range; [g] second singular spectrum analysis step: decomposing the third signal G3(t) into a plurality of second sub-sequences, and performing fast Fourier transform on all second sub-sequences to obtain a spectrum of the plurality of second sub-sequences, and then using the maximum amplitude of the spectrum in each second sub-sequence as a second main frequency; [h] a frequency mask screening step: only the second sub-sequence whose second main frequency is within the frequency mask is retained, and at least one second reserved sub-sequence can be selected. It becomes a fourth signal G4(t), where t=1 to N, the horizontal axis of the fourth signal G4(t) is time, and the vertical axis is the green value intensity; if the screening does not obtain any of the second reserved objects Sequence, the fourth signal G4(t) is directly equal to the third signal G3(t); . The second overlapping addition step: the plurality of the fourth signals G4(t) obtained after the second processing step are processed by the overlap addition technique of the existing raised cosine window method to obtain an overlap addition The four signal, which has the length of time, is the final heart rate signal. The method for calculating a heart rate extraction algorithm for a face compressed image according to claim 1, wherein in the step of acquiring the original image, two sets of video devices and a network connecting device are disposed, and the two sets of video devices pass through the Network connection device to reach a video contact. The heart rate extraction calculation method for the face compressed image according to the first aspect of the invention, wherein the compression method of the face compressed image is selected from one of x264, x265, vp8, and vp9.