欧拉运动放大算法实现心率检测

最新推荐文章于 2022-03-25 20:22:12 发布
最新推荐文章于 2022-03-25 20:22:12 发布
阅读量4k 收藏 73
点赞数 12
分类专栏： 数字图像处理 文章标签： c++ opencv 图像识别
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_43294951/article/details/106161760
版权

由于疫情的影响，一系列的非接触设备的到应用，如无人消毒机器人、红外测温。除了体温，心率也是衡量人体状况的重要指标。心率可以通过颈动脉、桡动脉还有皮肤颜色来进行提取。2012年，Wu等人基于欧拉视频放大(EVM)技术设计了一个非接触式心率检测系统在相机和目标人物本身没有明显晃动的情况下，且光照条件良好时，该系统测得心率的结果，能够获得临床级别的准确度。之后又有学者提出了基于相位的微小运动放大算法。近年来传统的运动放大技术还是EVM和PVM，自动识别和视觉增强方面的发展和探索空间还比较大。本文运用欧拉运动放大技术的原理，以增强胸骨上窝运动来获取心率。主要是通过非接触式拍摄一段视频，运用EVM放大运动，再提取心跳信号并分析得出心率。

先看一下效果吧
心率图
在这里插入图片描述放大胸骨上窝
在这里插入图片描述

本程序完全开源，网盘链接：https://pan.baidu.com/s/1XU0IhnPIKWYzQQN6bPG8pg
提取码：wj4k
如果对大家有用，希望可以给我的码云点个star~(github又被墙了)项目地址

那么什么是欧拉运动放大呢?简单的来说是在图像金字塔上进行操作，求出相邻两帧的图像变化，然后进行IIR滤波，还原通道，再进行放大和衰减。
由于处理域为YIQ，所以需要先将图像由RGB变换到YIQ，函数如下

Matx33f rgb2yiq_mat(0.299f, 0.587f, 0.114f,
	0.596f, -0.274f, -0.322f,
	0.211f, -0.523f, 0.312f);

Matx33f yiq2rgb_mat(1.0f, 0.956f, 0.621f,
	1.0f, -0.272f, -0.647f,
	1.0f, -1.106f, 1.703f);

//RGB2YIQ

Mat rgb2yiq(const Mat& img) {

	Mat img_out(img.size(), img.type());
	img_out = img.clone();
	//img_out.convertTo(img_out, CV_32F);

	for (int j = 0; j < img.rows; j++) {
		for (int i = 0; i < img.cols * 3; i += 3) {

			Vec3f pixel(img_out.at<float>(j, i + 2), img_out.at<float>(j, i + 1), img_out.at<float>(j, i));
			pixel = rgb2yiq_mat * pixel;


			for (int k = 0; k < 3; k++) {
				// if(pixel[k] > 255.0f) pixel[k] = 255.0f;
				// else if(pixel[k] < 0.0f) pixel[k] = 0.0f;

				img_out.at<float>(j, i + 2 - k) = pixel[k];
			}
		}
	}

	// img_out.convertTo(img_out, CV_8UC3);

	return img_out;
}

//YIQ2RGB
Mat yiq2rgb(const Mat& img) {

	Mat img_out(img.size(), img.type());
	img_out = img.clone();
	//img_out.convertTo(img_out, CV_32F);

	for (int j = 0; j < img.rows; j++) {
		for (int i = 0; i < img.cols * 3; i += 3) {

			Vec3f pixel(img_out.at<float>(j, i + 2), img_out.at<float>(j, i + 1), img_out.at<float>(j, i));
			pixel = yiq2rgb_mat * pixel;


			for (int k = 0; k < 3; k++) {
				if (pixel[k] > 255.0) pixel[k] = 255.0;
				else if (pixel[k] < 0.0) pixel[k] = 0.0;

				img_out.at<float>(j, i + 2 - k) = pixel[k];
			}
		}
	}

	//img_out.convertTo(img_out, CV_8UC3);

	return img_out;
}

然后构建拉普拉斯金字塔

bool laplacianpyramid(Mat img, int levels, vector<Mat_<Vec3f>>& pyramid)
{
    if (levels < 1)
    {
        cout << "level设置错误" << endl;
        return false;
    }
    else
    {

        Mat img1 = img.clone();
        for (int i = 0; i < levels - 1; i++)
        {
            Mat up, down;
            pyrDown(img1, down);
            pyrUp(down, up, img1.size());
            pyramid.push_back(img1 - up);
            img1 = down;
        }
        pyramid.push_back(img1);
    }
}

然后进行IIR滤波 注意代码中有一系列的文件操作，需要将这些文件代码进行注释(我是由于vector的copy导致数据出错，然后这样找bug的)

void IIR( vector<Mat_<Vec3f>>&src, vector<Mat_<Vec3f>>& dst, vector<Mat_<Vec3f>>& low, vector<Mat_<Vec3f>>& high)
{
    double fl = 0.05;
    double fh = 0.4;
    vector<Mat_<Vec3f>> low1, high1, dst1;
    for (int i = 0; i < src.size(); i++)
    {
        low1.push_back((1 - fl) * low[i] + fl * src[i]);
        high1.push_back((1 - fh) * high[i] + fh * src[i]);
        dst1.push_back(high[i]-low[i]);
    }
    low.assign(low1.begin(), low1.end());
    high.assign(high1.begin(), high1.end());
    dst.assign(dst1.begin(), dst1.end());
    ofstream file;
    file.open("dst.txt");
        for (int i = 0; i < dst[3].rows; i++)
        {
            for (int j = 0; j < dst[3].cols; j++)
            {
                file << float(dst[3].at<Vec3f>(i, j)[2]) << " ";
            }
            file << endl;
        }
    file.close();
    file.open("src.txt");
        for (int i = 0; i < dst[3].rows; i++)
        {
            for (int j = 0; j < dst[3].cols; j++)
            {
                file << src[3].at<Vec3f>(i, j)[2] << " ";
            }
            file << endl;
        }
    file.close();
    file.open("low.txt");
    int num = 3;
        for (int i = 0; i < dst[num].rows; i++)
        {
            for (int j = 0; j < dst[num].cols; j++)
            {
                file << low[num].at<Vec3f>(i, j)[2] << " ";
            }
            file << endl;
        }
    file.close();
    file.open("high.txt");
        for (int i = 0; i < dst[num].rows; i++)
        {
            for (int j = 0; j < dst[num].cols; j++)
            {
                file << high[num].at<Vec3f>(i, j)[2] << " ";
            }
            file << endl;
        }
    file.close();
}

进行完这些操作后，就需要恢复图像并进行放大处理

//恢复源图像并进行放大
void returnsrc(vector<Mat_<Vec3f>>& src, Mat& dst,double alpha,double lambda_c)
{
    double lambda = sqrtf(src[0].cols * src[0].cols + src[0].rows * src[0].rows) / 3;
    int ext = 2;
    double delta;
    delta = lambda_c / 8 / (1 + alpha);
    vector<Mat_<Vec3f>>temp;
    for (int i = 0; i < src.size(); i++)
    {
        double curralpha = lambda / delta / 8 - 1;
        curralpha = curralpha * ext;
        if ((i == 0) || (i == src.size() - 1))
        {
            src[i] = 0;
        }
        else if (curralpha > alpha)
        {
            src[i] = src[i]*alpha;
        }
        else
        {
            src[i] = src[i] *curralpha;
        }
        lambda = lambda / 2;
    }
    Mat up= src[src.size()-1];
    for (int i = 0; i < src.size()-1; i++)
    {
        pyrUp(up, up);
        up = up + src[src.size() - 2-i];
    }
    Mat planes[3];
    split(up, planes);
    planes[0] = planes[0] * 0.1;
    planes[1] = planes[1] * 0.1;
    merge(planes, 3, up);
    //Mat result = yiq2rgb(up);
    //result.convertTo(result, CV_8UC3, 255.0, 1.0 / 255.0);
    dst = up.clone();
}

这样通过两帧就可以得到一帧放大后的图像
但是怎么得到心率呢？
一种简单的方法是获取ROI，并求响应；另外一种是独立成分分析ICA，这种我现在还没有搞明白，先不用。下面介绍ROI的响应并求出周期

我对输入图像进行了裁剪，大小为512*512，30fps
获取ROI 代码中的cost就是代价

for (int i = 0; i < 16; i++)
        {
            for (int j = 0; j < 16; j++)
            {
                Mat imgrio = img(Rect(i*32,j*32,32,32));
                cost[flag - 1][j + i*16] = count(imgrio);
            }
        }
float count(Mat& src)
{
    Mat a;
    //a = src.clone();
    cvtColor(src, a, CV_RGB2GRAY);
    a.convertTo(a, CV_32FC1);
    float sum = 0;
    for (int i = 0; i < src.rows; i++)
    {
        for (int j = 0; j < src.cols; j++)
        {
            sum += a.at<float>(i, j);
        }
    }
    return sum/32/32;
}

然后你需要找到响应值最大的，可以采用冒泡排序，下面的为升序，find函数是找到每个ROI随时间的响应最大

void BubbleSort(float* p, int length, int* ind_diff)
{
    for (int m = 0; m < length; m++)
    {
        ind_diff[m] = m;
    }

    for (int i = 0; i < length; i++)
    {
        for (int j = 0; j < length - i - 1; j++)
        {
            if (p[j] > p[j + 1])
            {
                float temp = p[j];
                p[j] = p[j + 1];
                p[j + 1] = temp;

                int ind_temp = ind_diff[j];
                ind_diff[j] = ind_diff[j + 1];
                ind_diff[j + 1] = ind_temp;
            }
        }
    }
}

int find(float X[210][256])
{
    float average[256];
    float s[256];
    int ind[256] = { 0 };
    for (int j = 0; j < 256; j++)
    {
        average[j] = 0;
        s[j] = 0;
        for (int i = 0; i < 210; i++)
        {
            average[j] += X[i][j];
        }
        average[j] /= 210;
        for (int i = 0; i < 210; i++)
        {
            s[j] += (X[i][j]-average[j])*(X[i][j] - average[j]);
        }
    }
    BubbleSort(s, 256, ind);
    for (int i = 0; i < 256; i++)
    {
        cout << "value: " << s[i] << "Index: " << ind[i] << endl;
    }
    //return distance(s, max_element(s, s + sizeof(s) / sizeof(s[0])));
    return ind[255];
}

这样我们的任务基本上结束了 全部代码如下，其中需要自建pdi文件

pdi.h

#ifndef PDI_H
#define PDI_H

#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
#include <algorithm>

using namespace cv;

Mat filter(const Mat& img, const Mat& mask);
Mat filterRGB(const Mat& img, const Mat& mask);
Mat filterGray(const Mat& img, const Mat& mask);

Mat rgb2yiq2rgb(const Mat& img);
Mat rgb2yiq(const Mat& img);
Mat yiq2rgb(const Mat& img);
Mat split(const Mat& img, int channel, bool mono = false);
Mat merge(const Mat& img_R, const Mat& img_G, const Mat& img_B);
uchar mean_y_uchar(const Mat& img);
float mean_y_float(const Mat& img);

Mat negative(const Mat& img, bool rgb = true);
Mat negative_rgb(const Mat& img);
Mat negative_y(const Mat& img);
Mat addBrightness(const Mat& img, int c);
Mat mulBrightness(const Mat& img, float c, bool rgb = true);
Mat mulBrightness_rgb(const Mat& img, float c);
Mat mulBrightness_y(const Mat& img, float c);
Mat thresholding(const Mat& img, uchar threshold, int type, bool mean = true);
Mat thresholdingGray(const Mat& img, uchar threshold, int type, bool mean = true);
Mat thresholdingYIQ(const Mat& img, float threshold, int type, bool mean = true);


Mat sobelFilterVer(const Mat& img);
Mat sobelFilterHor(const Mat& img);
Mat sobelFilter(const Mat& img);
Mat laplacianFilterVer(const Mat& img);
Mat laplacianFilterHor(const Mat& img);
Mat laplacianFilter(const Mat& img);

Mat noiseSaltPepper(const Mat& img, float density);
Mat noiseSaltPepperMono(const Mat& img, float density);
Mat noiseSaltPepperMono(const Mat& img, float density);

#endif

pdi.cpp

#include "pdi.hpp"

Matx33f rgb2yiq_mat(0.299f, 0.587f, 0.114f,
	0.596f, -0.274f, -0.322f,
	0.211f, -0.523f, 0.312f);

Matx33f yiq2rgb_mat(1.0f, 0.956f, 0.621f,
	1.0f, -0.272f, -0.647f,
	1.0f, -1.106f, 1.703f);

// 1.1
Mat rgb2yiq2rgb(const Mat& img) {

	return yiq2rgb(rgb2yiq(img));
}

Mat rgb2yiq(const Mat& img) {

	Mat img_out(img.size(), img.type());
	img_out = img.clone();
	//img_out.convertTo(img_out, CV_32F);

	for (int j = 0; j < img.rows; j++) {
		for (int i = 0; i < img.cols * 3; i += 3) {

			Vec3f pixel(img_out.at<float>(j, i + 2), img_out.at<float>(j, i + 1), img_out.at<float>(j, i));
			pixel = rgb2yiq_mat * pixel;


			for (int k = 0; k < 3; k++) {
				// if(pixel[k] > 255.0f) pixel[k] = 255.0f;
				// else if(pixel[k] < 0.0f) pixel[k] = 0.0f;

				img_out.at<float>(j, i + 2 - k) = pixel[k];
			}
		}
	}

	// img_out.convertTo(img_out, CV_8UC3);

	return img_out;
}

Mat yiq2rgb(const Mat& img) {

	Mat img_out(img.size(), img.type());
	img_out = img.clone();
	//img_out.convertTo(img_out, CV_32F);

	for (int j = 0; j < img.rows; j++) {
		for (int i = 0; i < img.cols * 3; i += 3) {

			Vec3f pixel(img_out.at<float>(j, i + 2), img_out.at<float>(j, i + 1), img_out.at<float>(j, i));
			pixel = yiq2rgb_mat * pixel;


			for (int k = 0; k < 3; k++) {
				if (pixel[k] > 255.0) pixel[k] = 255.0;
				else if (pixel[k] < 0.0) pixel[k] = 0.0;

				img_out.at<float>(j, i + 2 - k) = pixel[k];
			}
		}
	}

	//img_out.convertTo(img_out, CV_8UC3);

	return img_out;
}

// 2.1
Mat split(const Mat& img, int channel, bool mono) {

	int type;

	if (mono)
		type = CV_8UC1;
	else
		type = CV_8UC3;

	Mat img_out(img.size(), type);

	for (int j = 0; j < img.rows; j++) {
		for (int i = 0; i < img.cols; i++) {

			img_out.at<uchar>(j, i + 2 - channel + !mono * 2 * i + mono * (channel - 2)) = img.at<uchar>(j, 3 * i + 2 - channel);
		}
	}

	return img_out;
}

Mat merge(const Mat& img_R, const Mat& img_G, const Mat& img_B) {

	Mat img_out(img_R.size(), CV_8UC3);

	for (int j = 0; j < img_R.rows; j++) {
		for (int i = 0; i < img_R.cols; i++) {
			img_out.at<uchar>(j, 3 * i + 2) = img_R.at<uchar>(j, i);
			img_out.at<uchar>(j, 3 * i + 1) = img_G.at<uchar>(j, i);
			img_out.at<uchar>(j, 3 * i) = img_B.at<uchar>(j, i);
		}
	}

	return img_out;
}

// 1.3
Mat negative(const Mat& img, bool rgb) {

	if (rgb)
		return negative_rgb(img);
	else
		return negative_y(img);
}

// 1.3a
Mat negative_rgb(const Mat& img) {

	bool color = true;
	if (img.type() == CV_32FC1 || img.type() == CV_8UC1)
		color = false;

	Mat img_out(img.size(), img.type());

	for (int j = 0; j < img.rows; j++) {
		for (int i = 0; i < img.cols + color * img.cols * 2; i++) {

			img_out.at<uchar>(j, i) = 255 - img.at<uchar>(j, i);
		}
	}

	return img_out;
}

// 1.3b
Mat negative_y(const Mat& img) {

	Mat img_out = rgb2yiq(img);

	for (int j = 0; j < img.rows; j++) {
		for (int i = 2; i < 3 * img.cols; i += 3) {

			img_out.at<float>(j, i) = 255 - img_out.at<float>(j, i);
		}
	}

	return yiq2rgb(img_out);
}

// 1.4
Mat addBrightness(const Mat& img, int c) {

	bool color = true;
	if (img.type() == CV_32FC1 || img.type() == CV_8UC1)
		color = false;

	Mat img_out(img.size(), img.type());

	for (int j = 0; j < img.rows; j++) {
		for (int i = 0; i < img.cols + color * img.cols * 2; i++) {

			int pixel = img.at<uchar>(j, i) + c;

			if (pixel > 255) pixel = 255;
			else if (pixel < 0) pixel = 0;

			img_out.at<uchar>(j, i) = pixel;
		}
	}

	return img_out;
}

//1.5
Mat mulBrightness(const Mat& img, float c, bool rgb) {
	if (rgb)
		return mulBrightness_rgb(img, c);
	else
		return mulBrightness_y(img, c);
}

// 1.5a
Mat mulBrightness_rgb(const Mat& img, float c) {

	if (c < 0.0f) return img;

	bool color = true;
	if (img.type() == CV_32FC1 || img.type() == CV_8UC1)
		color = false;

	Mat img_out(img.size(), img.type());

	for (int j = 0; j < img.rows; j++) {
		for (int i = 0; i < img.cols + color * img.cols * 2; i++) {

			float pixel = img.at<uchar>(j, i) * c;

			if (pixel > 255.0f) pixel = 255.0f;

			img_out.at<uchar>(j, i) = pixel;
		}
	}

	return img_out;
}

// 1.5b
Mat mulBrightness_y(const Mat& img, float c) {

	if (c < 0.0f) return img;

	Mat img_out = rgb2yiq(img);

	for (int j = 0; j < img.rows; j++) {
		for (int i = 2; i < img.cols * 3; i += 3) {

			float pixel = img_out.at<float>(j, i) * c;

			if (pixel > 255.0f) pixel = 255.0f;

			img_out.at<float>(j, i) = pixel;
		}
	}

	return yiq2rgb(img_out);
}

// 1.6
Mat thresholding(const Mat& img, uchar threshold, int type, bool mean) {

	Mat img_out;
	bool color = true;

	if (img.type() == CV_32FC1 || img.type() == CV_8UC1) {
		color = false;
		img_out = img;
		return thresholdingGray(img_out, threshold, type, mean);
	}
	else {
		img_out = rgb2yiq(img);
		return thresholdingYIQ(img_out, float(threshold), type, mean);
	}
}

Mat thresholdingGray(const Mat& img, uchar threshold, int type, bool mean) {

	Mat img_out = img;
	if (mean) threshold = mean_y_uchar(img_out);

	for (int j = 0; j < img.rows; j++) {
		for (int i = 0; i < img.cols; i++) {
			if (type == 0)
				img_out.at<uchar>(j, i) = img_out.at<uchar>(j, i) > threshold ? 255 : 0;
			else if (type == 1)
				img_out.at<uchar>(j, i) = img_out.at<uchar>(j, i) > threshold ? 255 : img_out.at<uchar>(j, i);
			else
				img_out.at<uchar>(j, i) = img_out.at<uchar>(j, i) > threshold ? threshold : img_out.at<uchar>(j, i);
		}
	}

	return img_out;
}

Mat thresholdingYIQ(const Mat& img, float threshold, int type, bool mean) {

	Mat img_out = img;
	if (mean) threshold = mean_y_float(img_out);

	for (int j = 0; j < img.rows; j++) {
		for (int i = 2; i < img.cols * 3; i += 3) {
			if (type == 0)
				img_out.at<float>(j, i) = img_out.at<float>(j, i) > threshold ? 255 : 0;
			else if (type == 1)
				img_out.at<float>(j, i) = img_out.at<float>(j, i) > threshold ? 255 : img_out.at<float>(j, i);
			else
				img_out.at<float>(j, i) = img_out.at<float>(j, i) > threshold ? threshold : img_out.at<float>(j, i);
		}
	}

	return yiq2rgb(img_out);
}

uchar mean_y_uchar(const Mat& img) {

	long mean = 0;

	for (int j = 0; j < img.rows; j++)
		for (int i = 0; i < img.cols; i++)
			mean += img.at<uchar>(j, i);

	mean = mean / (img.rows * img.cols);
	return mean;
}

float mean_y_float(const Mat& img) {

	float mean = 0;

	for (int j = 0; j < img.rows; j++)
		for (int i = 2; i < img.cols * 3; i += 3)
			mean += img.at<float>(j, i);

	mean = mean / (img.rows * img.cols);
	return mean;
}

Mat filter(const Mat& img, const Mat& mask) {

	if (img.type() == CV_32FC3 || img.type() == CV_8UC3)
		return filterRGB(img, mask);

	if (img.type() == CV_32FC1 || img.type() == CV_8UC1)
		return filterGray(img, mask);
}

Mat filterGray(const Mat& img, const Mat& mask) {

	Mat img_out(img.size(), img.type());
	int raio = (mask.rows - 1) / 2;

	Mat img_pad;
	copyMakeBorder(img, img_pad, raio, raio, raio, raio, BORDER_REFLECT101);

	for (int j = raio; j < img.rows + raio; j++) {

		for (int i = raio; i < img.cols + raio; i++) {

			Rect region(i - raio, j - raio, mask.rows, mask.cols);
			Mat img_region = img_pad(region);

			img_region.convertTo(img_region, CV_32F);

			img_region = img_region.mul(mask);

			float pixel = sum(img_region)[0];

			if (pixel > 255) { pixel = 255; }
			if (pixel < 0) { pixel = 0; }

			img_out.at<uchar>(j - raio, i - raio) = pixel;
		}
	}

	return img_out;
}

Mat filterRGB(const Mat& img, const Mat& mask) {

	Mat channel[3];

	for (int i = 0; i < 3; i++) {
		channel[i] = filterGray(split(img, i, true), mask);
	}

	return merge(channel[0], channel[1], channel[2]);
}

Mat sobelFilterHor(const Mat& img) {

	Mat	mask = (Mat_<float>(3, 3) <<
		-1, 0, 1,
		-2, 0, 2,
		-1, 0, 1);

	return filter(img, mask);
}

Mat sobelFilterVer(const Mat& img) {

	Mat	mask = (Mat_<float>(3, 3) <<
		1, 2, 1,
		0, 0, 0,
		-1, -2, -1);

	return filter(img, mask);
}

Mat sobelFilter(const Mat& img) {

	Mat img_out = sobelFilterVer(img) + sobelFilterHor(img);

	bool color = true;
	if (img_out.type() == CV_32FC1 || img_out.type() == CV_8UC1)
		color = false;

	for (int j = 0; j < img_out.rows; j++) {
		for (int i = 0; i < img_out.cols + color * img.cols * 2; i++) {

			uchar pixel = img_out.at<uchar>(j, i);

			if (pixel > 255) pixel = 255;
			else if (pixel < 0) pixel = 0;

			img_out.at<uchar>(j, i) = pixel;
		}
	}

	return img_out;
}

Mat laplacianFilterHor(const Mat& img) {

	Mat	mask = (Mat_<float>(3, 3) <<
		0, 1, 0,
		1, -4, 1,
		0, 1, 0);

	return filter(img, mask);
}

Mat laplacianFilterVer(const Mat& img) {

	Mat	mask = (Mat_<float>(3, 3) <<
		1, 1, 1,
		1, -8, 1,
		1, 1, 1);

	return filter(img, mask);
}

Mat laplacianFilter(const Mat& img) {

	// Mat img_out = laplacianFilterVer(img) + laplacianFilterHor(img);

	// bool color = true;
	// if(img_out.type() == CV_32FC1 || img_out.type() == CV_8UC1)
	// 	color = false;

	// for(int j = 0; j < img_out.rows; j++){
	// 	for(int i = 0; i < img_out.cols + color*img.cols*2; i++){

	// 		uchar pixel = img_out.at<uchar>(j, i);

	// 		if(pixel > 255) pixel = 255;
	// 		else if(pixel < 0) pixel = 0;

	// 		img_out.at<uchar>(j, i) = pixel;
	// 	}
	// }

	// return img_out;

	Mat	mask = (Mat_<float>(3, 3) <<
		0, 1, 0,
		1, -4, 1,
		0, 1, 0);

	return filter(img, mask);
}

Mat noiseSaltPepperMono(const Mat& img, float density) {

	// Mat black(img.size(), IPL_DEPTH_1U);
	// Mat white(img.size(), IPL_DEPTH_1U);
	uchar pixel;
	Mat out_img(img.size(), img.type());

	/*for (int j = 0; j < img.rows; j++) {
		for (int i = 0; i < img.cols; i++) {

			pixel = uchar(rand() % 256);
			out_img.at<uchar>(j, i) = pixel < uint(127 * density) ? 0 : img.at<uchar>(j, i);
			out_img.at<uchar>(j, i) = pixel > uint(127 * (2.0f - density)) ? 255 : out_img.at<uchar>(j, i);
		}
	}*/

	//randu(saltpepper_noise, 0, 255);

	// Mat black = saltpepper_noise < uint(127*density);
	// Mat white = saltpepper_noise > uint(127*(2.0-density));

	// Mat saltpepper_img = img.clone();

	// saltpepper_img.setTo(255, white);
	// saltpepper_img.setTo(0, black);

	return out_img;
}

Mat noiseSaltPepperRGB(const Mat& img, float density) {

	Mat channel[3];

	for (int i = 0; i < 3; i++) {
		channel[i] = noiseSaltPepperMono(split(img, i, true), density);
	}

	return merge(channel[0], channel[1], channel[2]);

}

Mat noiseSaltPepper(const Mat& img, float density) {

	if (img.type() == CV_32FC3 || img.type() == CV_8UC3)
		return noiseSaltPepperRGB(img, density);

	if (img.type() == CV_32FC1 || img.type() == CV_8UC1)
		return noiseSaltPepperMono(img, density);
}

main.cpp

#include <opencv2/opencv.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <vector>
#include <fstream>
#include<algorithm>
#include <iterator>
#include <valarray>
#include "pdi.hpp"
#include <graphics.h>      // 引用图形库头文件
#include <conio.h>


using namespace std;
using namespace cv;

//空域滤波--构建拉普拉斯金字塔
bool laplacianpyramid(Mat img, int levels, vector<Mat_<Vec3f>>& pyramid)
{
    if (levels < 1)
    {
        cout << "level设置错误" << endl;
        return false;
    }
    else
    {

        Mat img1 = img.clone();
        for (int i = 0; i < levels - 1; i++)
        {
            Mat up, down;
            pyrDown(img1, down);
            pyrUp(down, up, img1.size());
            pyramid.push_back(img1 - up);
            img1 = down;
        }
        pyramid.push_back(img1);
    }
}
//时域滤波--二阶IIR
void IIR( vector<Mat_<Vec3f>>&src, vector<Mat_<Vec3f>>& dst, vector<Mat_<Vec3f>>& low, vector<Mat_<Vec3f>>& high)
{
    double fl = 0.05;
    double fh = 0.4;
    vector<Mat_<Vec3f>> low1, high1, dst1;
    for (int i = 0; i < src.size(); i++)
    {
        low1.push_back((1 - fl) * low[i] + fl * src[i]);
        high1.push_back((1 - fh) * high[i] + fh * src[i]);
        dst1.push_back(high[i]-low[i]);
    }
    low.assign(low1.begin(), low1.end());
    high.assign(high1.begin(), high1.end());
    dst.assign(dst1.begin(), dst1.end());
    ofstream file;
    file.open("dst.txt");
        for (int i = 0; i < dst[3].rows; i++)
        {
            for (int j = 0; j < dst[3].cols; j++)
            {
                file << float(dst[3].at<Vec3f>(i, j)[2]) << " ";
            }
            file << endl;
        }
    file.close();
    file.open("src.txt");
        for (int i = 0; i < dst[3].rows; i++)
        {
            for (int j = 0; j < dst[3].cols; j++)
            {
                file << src[3].at<Vec3f>(i, j)[2] << " ";
            }
            file << endl;
        }
    file.close();
    file.open("low.txt");
    int num = 3;
        for (int i = 0; i < dst[num].rows; i++)
        {
            for (int j = 0; j < dst[num].cols; j++)
            {
                file << low[num].at<Vec3f>(i, j)[2] << " ";
            }
            file << endl;
        }
    file.close();
    file.open("high.txt");
        for (int i = 0; i < dst[num].rows; i++)
        {
            for (int j = 0; j < dst[num].cols; j++)
            {
                file << high[num].at<Vec3f>(i, j)[2] << " ";
            }
            file << endl;
        }
    file.close();
}
//恢复源图像并进行放大
void returnsrc(vector<Mat_<Vec3f>>& src, Mat& dst,double alpha,double lambda_c)
{
    double lambda = sqrtf(src[0].cols * src[0].cols + src[0].rows * src[0].rows) / 3;
    int ext = 2;
    double delta;
    delta = lambda_c / 8 / (1 + alpha);
    vector<Mat_<Vec3f>>temp;
    for (int i = 0; i < src.size(); i++)
    {
        double curralpha = lambda / delta / 8 - 1;
        curralpha = curralpha * ext;
        if ((i == 0) || (i == src.size() - 1))
        {
            src[i] = 0;
        }
        else if (curralpha > alpha)
        {
            src[i] = src[i]*alpha;
        }
        else
        {
            src[i] = src[i] *curralpha;
        }
        lambda = lambda / 2;
    }
    Mat up= src[src.size()-1];
    for (int i = 0; i < src.size()-1; i++)
    {
        pyrUp(up, up);
        up = up + src[src.size() - 2-i];
    }
    Mat planes[3];
    split(up, planes);
    planes[0] = planes[0] * 0.1;
    planes[1] = planes[1] * 0.1;
    merge(planes, 3, up);
    //Mat result = yiq2rgb(up);
    //result.convertTo(result, CV_8UC3, 255.0, 1.0 / 255.0);
    dst = up.clone();
}

//void returnsrc(vector<Mat_<Vec3f>>& src, Mat& dst)
//{
//    vector<Mat_<Vec3f>>temp;
//    for (int i = 0; i < src.size(); i++)
//    {
//        Mat a = src[i];
//        Mat p = yiq2rgb(a);
//        Mat planes[3];
//        split(p, planes);
//        planes[1] = planes[1] * 100;
//        planes[2] = planes[2] * 100;
//        merge(planes, 3, p);
//        temp.push_back(p);
//    }
//    Mat up = temp[src.size() - 1];
//    for (int i = 0; i < src.size() - 1; i++)
//    {
//        pyrUp(up, up);
//        up = up + temp[src.size() - 2 - i];
//    }
//    up.convertTo(up, CV_8UC3, 255.0, 1.0 / 255.0);
//    dst = up.clone();
//}
float count(Mat& src)
{
    Mat a;
    //a = src.clone();
    cvtColor(src, a, CV_RGB2GRAY);
    a.convertTo(a, CV_32FC1);
    float sum = 0;
    for (int i = 0; i < src.rows; i++)
    {
        for (int j = 0; j < src.cols; j++)
        {
            sum += a.at<float>(i, j);
        }
    }
    return sum/32/32;
}
void BubbleSort(float* p, int length, int* ind_diff)
{
    for (int m = 0; m < length; m++)
    {
        ind_diff[m] = m;
    }

    for (int i = 0; i < length; i++)
    {
        for (int j = 0; j < length - i - 1; j++)
        {
            if (p[j] > p[j + 1])
            {
                float temp = p[j];
                p[j] = p[j + 1];
                p[j + 1] = temp;

                int ind_temp = ind_diff[j];
                ind_diff[j] = ind_diff[j + 1];
                ind_diff[j + 1] = ind_temp;
            }
        }
    }
}
int find(float X[210][256])
{
    float average[256];
    float s[256];
    int ind[256] = { 0 };
    for (int j = 0; j < 256; j++)
    {
        average[j] = 0;
        s[j] = 0;
        for (int i = 0; i < 210; i++)
        {
            average[j] += X[i][j];
        }
        average[j] /= 210;
        for (int i = 0; i < 210; i++)
        {
            s[j] += (X[i][j]-average[j])*(X[i][j] - average[j]);
        }
    }
    BubbleSort(s, 256, ind);
    for (int i = 0; i < 256; i++)
    {
        cout << "value: " << s[i] << "Index: " << ind[i] << endl;
    }
    //return distance(s, max_element(s, s + sizeof(s) / sizeof(s[0])));
    return ind[255];
}



int main()
{
//#define find_plot
#define creat
#ifdef creat
    VideoCapture frame("04032.mp4");
    VideoWriter writer("result_ht_1.avi", CV_FOURCC('M', 'J', 'P', 'G'), 30, Size(512, 512));
    Mat img,out;
    vector<Mat_<Vec3f>>result;
    vector<Mat_<Vec3f>>low;
    vector<Mat_<Vec3f>>high;
    int number = 0;
#define level 7    
    while (1)
    {
        vector<Mat_<Vec3f>>pyramid;
        frame >> img;
        string a = to_string(number+1)+".jpg";
        imwrite(a, img);
        if (img.empty())
        {
            break;
        }
        else
        {
            img.convertTo(img, CV_32FC3, 1.0 / 255.0f);
            Mat imgyiq;
            imgyiq = rgb2yiq(img);
            Mat imgba = imgyiq.clone();
            laplacianpyramid(imgyiq, level, pyramid);
            if (number == 0)
            {
                low.assign(pyramid.begin(), pyramid.end());
                high.assign(pyramid.begin(), pyramid.end());
                result.assign(pyramid.begin(), pyramid.end());
            }
            else
            {
                IIR(pyramid, result, low, high);
                returnsrc(result, out,15,24);
                //imshow("", out);
                //imshow("s", imgba);
                out = (out + imgba);
                Mat finnal = yiq2rgb(out);
                finnal.convertTo(finnal, CV_8UC3, 255.0);
                writer << finnal;
                imshow("", finnal);
                waitKey(1);
            }
            number++;
        }
    }
#endif
#ifdef find_plot
    VideoCapture frame("result_ht_1.avi");
    Mat img;
    int flag = 0;
    float cost[210][256];
    while (true)
    {
        frame >> img;
        //img.convertTo(img, CV_32FC3, 1.0 / 255.0f);
        //Mat imgyiq = rgb2yiq(img);
        //Mat planes[3];
        //split(imgyiq, planes);
        flag++;
        if (img.empty())
        {
            break;
        }
        if (flag > 210)
        {
            break;
        }
        for (int i = 0; i < 16; i++)
        {
            for (int j = 0; j < 16; j++)
            {
                Mat imgrio = img(Rect(i*32,j*32,32,32));
                cost[flag - 1][j + i*16] = count(imgrio);
            }
        }
        
    }

    int loc = find(cost);
    cout << loc << endl;
    ofstream file;
    file.open("zhenfu.txt");
    for (int i = 0; i < 210; i++)
    {
        file << float(cost[i][loc]) << ",";
    }
    file.close();
    /*initgraph(420,200);
    setbkcolor(WHITE);
    cleardevice();
    for (int i = 0; i < 210; i++)
    {
        putpixel(2 * i, 200-int(float(cost[i][loc])), BLACK);
    }
    saveimage(_T("test.bmp"));
    _getch();
    closegraph();*/
#endif
}
鹏鹏诶
关注
专栏目录
基于面部视频的实时心率监测系统 day eight
周子博的博客
03-16 1217
今天实现对面部视频的欧拉放大 将感兴趣区域进行欧放大 ，欧拉放大的几个步骤，参考CSDN上guifeng93 对欧拉放大的实现步骤 1,RGB转换为 YIQ 。YIQ 是 NTSC 电视机系统所采用的颜色空间，Y 是提供黑白电视及彩色电视的亮度信号，I 代表 In-phase，色彩从橙色到青色，Q 代表 Quadrature-phase，色彩从紫色到黄绿色。 https://blog.csdn....
微小运动检测（一）欧拉运动放大算法
lianghong1564868的博客
09-19 6074
一、关于微小运动      自然界中包含很多形式的运动，很多是肉眼不可见的，比如人脸面部血液的流动，然而这种微小运动却包含了很多信息，比如我们可以通过分析面部血液流动的频率，进而得到心率。MIT提出的欧拉算法可以放大这种不可见的运动，使之清晰的呈现在我们面前。 二、欧拉运动放大算法的基本原理     &nbs...
eulerian_video_magnification:实现欧拉视频放大并用于心率检测等
04-27
eulerian_video_magnification 实现欧拉视频放大并用于心率检测等 依赖库 opencv2.4.9及以上（后续版本会将opencv中的傅立叶变换替换为FFTW） 功能 根据中的欧拉视频放大原理实现了运动放大和颜色放大两种算法。目前只支持读入视频来实现放大。 即将加入的功能 对人脸进行检测和跟踪。 将程序中用到的第三方GraphUtils画图库重构，主要是因为这个库用的是比较老版本的opencv，数据结构和函数接口都是C版本的，很容易出现内存方面的问题，用起来也不方便。 实现心率的检测。 MIT之后的论文有实现基于相位的放大，比这个版本更好，因为会把放大的信号加入一个相位平移，从而不会放大噪声，所以之后会考虑实现新版本的程序。 效果 原视频 运动放大 颜色放大
基于人脸视频的实时心率监测系统 DAY 15
周子博的博客
05-01 3806
欧拉放大已经实现 在网上查找了主流的心率特征提取方式 1.对脉搏跳动视频帧的像素点亮度值在时域上进行亮度方差统计，同时在Yc吒b颜色空间中分割出皮肤区域；其次，根据亮度方差统计和皮肤分割结果，结合图像形态学处理方法提取视频中桡动脉跳动区域；最后，对所提取桡动脉部位时域上亮度信号采用傅里叶变换进行时频分析，实现心率非接触式测量。 将图片帧转换为HSV ...
PyEVM:EVM（欧拉视频放大率）的Python实现
05-11
EVM（欧拉视频放大率）的Python实现 这是欧拉视频放大倍数《的欧拉视频放大倍率》的python实现。 我们的目标是揭示用肉眼难以或无法看到的视频中的时间变化，并以指示性的方式进行显示。 我们的方法（称为欧拉视频放大率）采用标准视频序列作为输入，并应用空间分解，然后对帧进行时间滤波。 然后将得到的信号放大以显示隐藏的信息。使用我们的方法，我们可以可视化血液在脸部填充时的流动，并且可以放大和显示微小的运动。 我们的技术可以实时运行，以显示在用户选择的时间频率处发生的现象。 安装OpenCV3 由于OpenCV3.X不支持Python3，因此您需要手动安装opencv3。 首先，下载适用于python3的opencv3： 适用于Python3的OpenCV3： ：//www.lfd.uci.edu/~gohlke/pythonlibs/#opencv 然后使用pip安装open
视频运动放大一：欧拉视频放大，《Eulerian video magnification for revealing subtle changes in the world》
Alocus的博客
03-24 8923
《Eulerian video magnification for revealing subtle changes in the world》（EVM） 目录 ​ 《Eulerian video magnification for revealing subtle changes in the world》（EVM） 背景 方法 实验 总结 背景 视频的时空处理相结合可以放大细微的变化，从而揭示周围世界的一些重要方面。如人类皮肤的颜色随血液循环而略有变化，细微表情和微小...
欧拉影像放大算法（Eulerian Video Magnification）的原理和实现
Weng Shihao's blog
03-25 4151
原文链接 因为原文网站不能翻译LaTeXLaTeXLaTeX公式了，后续在本站的一些转载也没有公式，故在此重新整理。 0.引言 人类的视觉感知存在有限的感知域。对于超出感知域的变化，我们无法感知。然而，这类信号却可能蕴藏着惊人的秘密。 比如，血液循环使得人体的皮肤发生细微的周期性变化，这个裸眼无法感知的变化却和人的心率非常吻合。2011 年，MIT 的一个亚裔学生 Mingzhe Poh 就利用这个微弱的信号设计了一个“魔镜” —— 不仅能照出你的模样，还能测出你的心率。 Mingzhe Poh 的这面神奇
心率算法
08-25 2720
class HeartUtils { private var lastMaxPos = 0 private var running = true private var index =0 private var maxPos = 0 private var maxValue = 0f fun init(index:Int,value:Float):Boolean{ if(maxValue<value){ m.
C++实现的欧拉放大代码
10-12
欧拉放大算法是一种在数值分析领域中用于模拟动力系统行为的方法。它是由数学家莱昂哈德·欧拉在18世纪提出的，是初等数值积分的一种基础技术。在这个特定的上下文中，"C++实现的欧拉放大代码" 提供了一个用C++编程...
心率算法_C语言_心率_
09-29
心率算法。通过心动描记法在逐个心律基础上估算心率-一种混合方法。
欧拉视频放大EVM（matlab）.zip
06-06
这是利用欧拉算法进行视频放大的源程序，关于算法的原理可以参考我的博客
心率算法（参照网上资料提供的算法编程）
05-20
参照网上资料提供的算法编程，参数是心率的值的buf，返回值是心率值。
实时心率算法
01-25
数据源为下位机采集到的心电电压值，取5s数据作为适应性数据，然后采用每秒数据迭代法计算实时心率值，1s刷新一次结果
QtEVM, 在OpenCV和Qt上，基于 C 的EVM ( 欧拉视频放大) 实现.zip
09-18
QtEVM, 在OpenCV和Qt上，基于 C 的EVM ( 欧拉视频放大) 实现 QtEVM基于OpenCV和Qt的EVM ( Eulerian视频放大器)的C 实现。支持：运动放大空间滤波器：拉普拉斯金字塔时间滤波器：IIR带通滤波器颜色放大空间滤波器：高斯金字塔时间滤波器：理想带通滤波
EVM算法的Python实现
热门推荐
flyingzhao
09-27 1万+
EVM(Eulerian Video Magnification,欧拉视频放大或欧拉影像放大)是一种将视频中的微小变化进行放大的算法，正如论文标题所说《Eulerian Video Magnification for Revealing Subtle Changes in the World》，该算法可以将视频中的微小变化转变为肉眼可以
视频的运动放大：欧拉方法与加速度联系
050100150200250RawPhase-basedOurs5290视频加速放大0Yichao Zhang，Silvia L. Pintea和Jan C. vanGemert视觉实验室，代尔夫特理工大学，代尔夫特，荷兰0100 120 ...（b）基于相位的运动放大[26]。（c）我们提
化秋毫为波澜：运动放大算法(深度学习版)
SIGAI_CSDN的博客
10-17 3528
什么是运动放大(Motion Magnification)？ 将视频中对应位置的运动进行放大，简单理解的话，就是找到时间段内的运动矢量，进行放大，然后权值叠加回去。 为什么需要运动放大？ 因为很多自然界或者生物的 subtle behaviour 不易被肉眼察觉（如飞机翼的震动，受风影响摇晃的建筑，生物皮肤变化等等），这些微变化只有通过运动放大，才能更好地被机器或者人类来做后续的视频...
欧拉图算法
Avenger_Tao的专栏
04-20 2159
首先考虑是欧拉图。若从一个顶点出发，利用DFS，将新搜索到的顶点保存到栈顶，并且将走过的边给删除，那么一定能访问到这个初始顶点。因为除起始顶点外其他的顶点若被访问到那么一定能去访问别的顶点。但是起始顶点不同，它是被直接访问的，即没有通过别的顶点去访问它，因此我们一定能找到这样的回路，使得起始顶点和终点是同一顶点，设为v，并且v的所有边都在这条回路上。我们用栈来保存这些顶点（按访问的顺序），如果此时
欧拉放大算法 python
最新发布
11-10
在Python中，我们可以使用OpenCV库来实现欧拉放大算法。具体步骤如下： 首先，导入需要的库： ```python import cv2 import numpy as np from scipy import interpolate ``` 然后，读取输入图像或视频序列： ```...
评论 14
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值