Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

added peak detection #30

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

added peak detection #30

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
100 changes: 93 additions & 7 deletions RPPG.cpp
100644 → 100755
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@ void RPPG::processFrame(Mat &frameRGB, Mat &frameGray, int time) {

// Set time
this->time = time;

cout<<"time "<<time<<endl;
if (!faceValid) {

cout << "Not valid, finding a new face" << endl;
Expand All @@ -109,9 +109,11 @@ void RPPG::processFrame(Mat &frameRGB, Mat &frameGray, int time) {
} else {

cout << "Tracking face" << endl;

// cout<<"corners.size() "<<corners.size()<<endl;
trackFace(frameGray);
}

// cout<<"frameGray() "<<frameGray.rows<<" , "<<frameGray.cols<<endl;

if (faceValid) {

Expand All @@ -124,11 +126,12 @@ void RPPG::processFrame(Mat &frameRGB, Mat &frameGray, int time) {
push(t);
push(re);
}

cout<<"s.rows "<<s.rows<<" t.rows "<<t.rows<<" re.rows "<<re.rows<<endl;
assert(s.rows == t.rows && s.rows == re.rows);

// New values
Scalar means = mean(frameRGB, mask);
// cout<<"means "<<means<<endl;
// Add new values to raw signal buffer
double values[] = {means(0), means(1), means(2)};
s.push_back(Mat(1, 3, CV_64F, values));
Expand All @@ -143,10 +146,11 @@ void RPPG::processFrame(Mat &frameRGB, Mat &frameGray, int time) {
// Update band spectrum limits
low = (int)(s.rows * LOW_BPM / SEC_PER_MIN / fps);
high = (int)(s.rows * HIGH_BPM / SEC_PER_MIN / fps) + 1;
// cout<<"band spectrum limits "<<low<<" , "<<high;

// If valid signal is large enough: estimate
if (s.rows >= fps * minSignalSize) {

// cout<<"s.rows >= fps * minSignalSize";
// Filtering
switch (rPPGAlg) {
case g:
Expand Down Expand Up @@ -287,6 +291,12 @@ void RPPG::trackFace(Mat &frameGray) {
vector<uchar> cornersFound_0;
Mat err;

// cout<<"frameGray() "<<frameGray.rows<<" , "<<frameGray.cols<<endl;
// cout<<"frameGray.type() "<<frameGray.type();
//
// cout<<"lastFrameGray() "<<lastFrameGray.rows<<" , "<<lastFrameGray.cols<<" , "<<lastFrameGray.type()<<endl;


// Track face features with Kanade-Lucas-Tomasi (KLT) algorithm
calcOpticalFlowPyrLK(lastFrameGray, frameGray, corners, corners_1, cornersFound_1, err);

Expand Down Expand Up @@ -345,6 +355,8 @@ void RPPG::trackFace(Mat &frameGray) {
void RPPG::updateROI() {
this->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 "<<roi<<endl;
}

void RPPG::updateMask(Mat &frameGray) {
Expand All @@ -366,11 +378,15 @@ void RPPG::invalidateFace() {
}

void RPPG::extractSignal_g() {

cout<<"extractSignal_g()"<<endl;

// Denoise
Mat s_den = Mat(s.rows, 1, CV_64F);
// cout<<"D s.rows "<<s.rows<<" s_den "<<s_den<<endl;

denoise(s.col(1), re, s_den);

// cout<<"E s.rows "<<s.rows<<" s_den "<<s_den<<endl;

// Normalise
normalization(s_den, s_den);

Expand Down Expand Up @@ -601,7 +617,7 @@ void RPPG::draw(cv::Mat &frameRGB) {
rectangle(frameRGB, roi, GREEN);

// Draw bounding box
rectangle(frameRGB, box, RED);
rectangle(frameRGB, box, WHITE);

// Draw signal
if (!s_f.empty() && !powerSpectrum.empty()) {
Expand All @@ -610,6 +626,9 @@ void RPPG::draw(cv::Mat &frameRGB) {
double displayHeight = box.height/2.0;
double displayWidth = box.width*0.8;

vector<ld> inputY;
inputY.push_back(s_f.at<double>(0, 0));

// Draw signal
double vmin, vmax;
Point pmin, pmax;
Expand All @@ -624,8 +643,31 @@ void RPPG::draw(cv::Mat &frameRGB) {
p2 = Point(drawAreaTlX + i * widthMult, drawAreaTlY + (vmax - s_f.at<double>(i, 0))*heightMult);
line(frameRGB, p1, p2, RED, 2);
p1 = p2;

inputY.push_back(s_f.at<double>(i, 0));
}

//draw zScore peak data
int lag = 5; //30;
ld threshold = 3.5; //5.0;
ld influence = 0.5; //0.0;
unordered_map<string, vector<ld>> 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);
Expand Down Expand Up @@ -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<string, vector<ld>> RPPG::z_score_thresholding(vector<ld> input, int lag, ld threshold, ld influence) {
unordered_map<string, vector<ld>> output;

uint n = (uint) input.size();
vector<ld> signals(input.size());
vector<ld> filtered_input(input.begin(), input.end());
vector<ld> filtered_mean(input.size());
vector<ld> 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;
}
82 changes: 81 additions & 1 deletion RPPG.hpp
100644 → 100755
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,84 @@

#include <stdio.h>

//---added zScore code
#include <vector>
#include <algorithm>
#include <unordered_map>
#include <cmath>
#include <iterator>
#include <numeric>

typedef long double ld;
typedef unsigned int uint;
typedef std::vector<ld>::iterator vec_iter_ld;

/**
* Overriding the ostream operator for pretty printing vectors.
*/
template<typename T>
std::ostream &operator<<(std::ostream &os, std::vector<T> vec) {
os << "[";
if (vec.size() != 0) {
std::copy(vec.begin(), vec.end() - 1, std::ostream_iterator<T>(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<ld> 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;
Expand Down Expand Up @@ -45,6 +123,8 @@ class RPPG {

typedef vector<Point2f> Contour2f;

unordered_map<string, vector<ld>> z_score_thresholding(vector<ld> input, int lag, ld threshold, ld influence);

private:

void detectFace(Mat &frameRGB, Mat &frameGray);
Expand Down Expand Up @@ -86,7 +166,7 @@ class RPPG {
int64_t lastSamplingTime;
int64_t lastScanTime;
int low;
int64_t now;
// int64_t now;
bool faceValid;
bool rescanFlag;

Expand Down