From eea31bd5624756783a37d95c6db78499853ba326 Mon Sep 17 00:00:00 2001 From: stephan schulz Date: Fri, 13 Nov 2020 09:13:04 -0500 Subject: [PATCH] added peak detection using this c++ implementation of the using z-scores peak detection algorithm https://stackoverflow.com/questions/22583391/peak-signal-detection-in-realtime-timeseries-data/46998001#46998001 --- RPPG.cpp | 100 +++++++++++++++++++++++++++++++++++++++++++++++++++---- RPPG.hpp | 82 ++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 174 insertions(+), 8 deletions(-) mode change 100644 => 100755 RPPG.cpp mode change 100644 => 100755 RPPG.hpp diff --git a/RPPG.cpp b/RPPG.cpp old mode 100644 new mode 100755 index 4390aaa..f4f30ea --- a/RPPG.cpp +++ b/RPPG.cpp @@ -90,7 +90,7 @@ void RPPG::processFrame(Mat &frameRGB, Mat &frameGray, int time) { // Set time this->time = time; - + cout<<"time "<= fps * minSignalSize) { - +// cout<<"s.rows >= fps * minSignalSize"; // Filtering switch (rPPGAlg) { case g: @@ -287,6 +291,12 @@ void RPPG::trackFace(Mat &frameGray) { vector cornersFound_0; Mat err; +// cout<<"frameGray() "<roi = Rect(Point(box.tl().x + 0.3 * box.width, box.tl().y + 0.1 * box.height), Point(box.tl().x + 0.7 * box.width, box.tl().y + 0.25 * box.height)); + +// cout<<"updateROI() roi "< inputY; + inputY.push_back(s_f.at(0, 0)); + // Draw signal double vmin, vmax; Point pmin, pmax; @@ -624,8 +643,31 @@ void RPPG::draw(cv::Mat &frameRGB) { p2 = Point(drawAreaTlX + i * widthMult, drawAreaTlY + (vmax - s_f.at(i, 0))*heightMult); line(frameRGB, p1, p2, RED, 2); p1 = p2; + + inputY.push_back(s_f.at(i, 0)); } + //draw zScore peak data + int lag = 5; //30; + ld threshold = 3.5; //5.0; + ld influence = 0.5; //0.0; + unordered_map> output = z_score_thresholding(inputY, lag, threshold, influence); + cout << output["signals"] << endl; + + cv::Point p3(drawAreaTlX, drawAreaTlY + (vmax - output["signals"][0])*heightMult); + cv::Point p4; + for (int i = 1; i < output["signals"].size(); i++) { + p4 = cv::Point(drawAreaTlX + i * widthMult, drawAreaTlY + (vmax - output["signals"][i])*heightMult); + line(frameRGB, p3, p4, WHITE, 2); + p3 = p4; + } + + if(output["signals"][output["signals"].size()-1] == 1){ + circle( frameRGB, cv::Point(50,50),30,Scalar( 0, 0, 255 ),FILLED,LINE_8 ); + } else{ + circle( frameRGB, cv::Point(50,50),30,Scalar( 0, 0, 255 ),1,LINE_8 ); + } + // Draw powerSpectrum const int total = s_f.rows; Mat bandMask = Mat::zeros(s_f.size(), CV_8U); @@ -664,3 +706,47 @@ void RPPG::draw(cv::Mat &frameRGB) { line(frameRGB, Point(corners[i].x,corners[i].y-5), Point(corners[i].x,corners[i].y+5), GREEN, 1); } } + +/** + * This is the implementation of the Smoothed Z-Score Algorithm. + * This is direction translation of https://stackoverflow.com/a/22640362/1461896. + * + * @param input - input signal + * @param lag - the lag of the moving window + * @param threshold - the z-score at which the algorithm signals + * @param influence - the influence (between 0 and 1) of new signals on the mean and standard deviation + * @return a hashmap containing the filtered signal and corresponding mean and standard deviation. + */ +//https://stackoverflow.com/questions/22583391/peak-signal-detection-in-realtime-timeseries-data/46998001#46998001 +unordered_map> RPPG::z_score_thresholding(vector input, int lag, ld threshold, ld influence) { + unordered_map> output; + + uint n = (uint) input.size(); + vector signals(input.size()); + vector filtered_input(input.begin(), input.end()); + vector filtered_mean(input.size()); + vector filtered_stddev(input.size()); + + VectorStats lag_subvector_stats(input.begin(), input.begin() + lag); + filtered_mean[lag - 1] = lag_subvector_stats.mean(); + filtered_stddev[lag - 1] = lag_subvector_stats.standard_deviation(); + + for (int i = lag; i < n; i++) { + if (abs(input[i] - filtered_mean[i - 1]) > threshold * filtered_stddev[i - 1]) { + signals[i] = (input[i] > filtered_mean[i - 1]) ? 1.0 : -1.0; + filtered_input[i] = influence * input[i] + (1 - influence) * filtered_input[i - 1]; + } else { + signals[i] = 0.0; + filtered_input[i] = input[i]; + } + VectorStats lag_subvector_stats(filtered_input.begin() + (i - lag), filtered_input.begin() + i); + filtered_mean[i] = lag_subvector_stats.mean(); + filtered_stddev[i] = lag_subvector_stats.standard_deviation(); + } + + output["signals"] = signals; + output["filtered_mean"] = filtered_mean; + output["filtered_stddev"] = filtered_stddev; + + return output; +} diff --git a/RPPG.hpp b/RPPG.hpp old mode 100644 new mode 100755 index f8011a1..07c5f49 --- a/RPPG.hpp +++ b/RPPG.hpp @@ -16,6 +16,84 @@ #include +//---added zScore code +#include +#include +#include +#include +#include +#include + +typedef long double ld; +typedef unsigned int uint; +typedef std::vector::iterator vec_iter_ld; + +/** + * Overriding the ostream operator for pretty printing vectors. + */ +template +std::ostream &operator<<(std::ostream &os, std::vector vec) { + os << "["; + if (vec.size() != 0) { + std::copy(vec.begin(), vec.end() - 1, std::ostream_iterator(os, " ")); + os << vec.back(); + } + os << "]"; + return os; +} + +/** + * This class calculates mean and standard deviation of a subvector. + * This is basically stats computation of a subvector of a window size qual to "lag". + */ + +class VectorStats { +public: + /** + * Constructor for VectorStats class. + * + * @param start - This is the iterator position of the start of the window, + * @param end - This is the iterator position of the end of the window, + */ + VectorStats(vec_iter_ld start, vec_iter_ld end) { + this->start = start; + this->end = end; + this->compute(); + } + + /** + * This method calculates the mean and standard deviation using STL function. + * This is the Two-Pass implementation of the Mean & Variance calculation. + */ + void compute() { + ld sum = std::accumulate(start, end, 0.0); + uint slice_size = std::distance(start, end); + ld mean = sum / slice_size; + std::vector diff(slice_size); + std::transform(start, end, diff.begin(), [mean](ld x) { return x - mean; }); + ld sq_sum = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0); + ld std_dev = std::sqrt(sq_sum / slice_size); + + this->m1 = mean; + this->m2 = std_dev; + } + + ld mean() { + return m1; + } + + ld standard_deviation() { + return m2; + } + +private: + vec_iter_ld start; + vec_iter_ld end; + ld m1; + ld m2; +}; + +//----original RPPG code using namespace cv; using namespace dnn; using namespace std; @@ -45,6 +123,8 @@ class RPPG { typedef vector Contour2f; + unordered_map> z_score_thresholding(vector input, int lag, ld threshold, ld influence); + private: void detectFace(Mat &frameRGB, Mat &frameGray); @@ -86,7 +166,7 @@ class RPPG { int64_t lastSamplingTime; int64_t lastScanTime; int low; - int64_t now; +// int64_t now; bool faceValid; bool rescanFlag;