#include #include #include #include #include //VideoCapture #include //cvtColor /* Compute Mean Square Error - slow version */ double computeMSE0(cv::Mat m1, cv::Mat m2) { assert(m1.size() == m2.size() && m1.type() == m2.type() && m1.type() == CV_8UC1); unsigned long sum = 0; for(int i = 0; i < m1.rows; ++i) { for(int j = 0; j < m1.cols; ++j) { const int diff = m1.at(i, j) - m2.at(i, j); sum += diff * diff; } } const double MSE = sum / (double)(m1.rows*m1.cols); return MSE; } /* Compute Mean Square Error - fast version */ double computeMSE1(const cv::Mat &m1, const cv::Mat &m2) { assert(m1.size() == m2.size() && m1.type() == m2.type() && m1.type() == CV_8UC1); unsigned long sum = 0; int rows = m1.rows; int cols = m1.cols; if (m1.isContinuous() && m2.isContinuous()) { cols *= rows; rows = 1; } for (int y = 0; y < rows; ++y) { const unsigned char *p1 = m1.ptr(y); const unsigned char *p2 = m2.ptr(y); for(int x = 0; x < cols; ++x) { const int diff = p1[x]-p2[x]; sum += diff * diff; } } const double MSE = sum / (double)(m1.rows*m1.cols); return MSE; } /* Compute Mean Square Error - with matrix operations */ double computeMSE2(cv::Mat m1, cv::Mat m2) { assert(m1.size() == m2.size() && m1.type() == m2.type() && m1.type() == CV_8UC1); cv::Mat m; cv::absdiff(m1, m2, m); //stay on 8 bits m.convertTo(m, CV_32F); //convert to float (32 bits) to hold square values m = m.mul(m); const double sum = (cv::sum(m))[0]; const double MSE = sum / m.total(); return MSE; } /* Compute Mean Square Error - with matrix operations (earlier conversion to float) */ double computeMSE3(cv::Mat m1, cv::Mat m2) { assert(m1.size() == m2.size() && m1.type() == m2.type() && m1.type() == CV_8UC1); m1.convertTo(m1, CV_32FC1); m2.convertTo(m2, CV_32FC1); cv::Mat e; cv::subtract(m1, m2, e); e = e.mul(e); const double sum = (cv::sum(e))[0]; const double MSE = sum / e.total(); return MSE; } /* Compute Mean Square Error - with matrix operations */ double computeMSE5(cv::Mat m1, cv::Mat m2) { const double diff_frames = cv::pow(cv::norm(m1, m2, cv::NORM_L2), 2.0); const double MSE = diff_frames/(m1.rows*m1.cols); return MSE; } /* double computeMSE_WRONG(cv::Mat m1, cv::Mat m2) { assert(m1.size() == m2.size() && m1.type() == m2.type() && m1.type() == CV_8UC1); cv::Mat m; cv::pow(m1-m2, 2, m); //WRONG ! const double MSE = cv::sum(m)[0] / m.total(); return MSE; } */ /* Compute Peak Signal-to-Noise Ratio */ double computePSNR(double MSE) { const double PSNR = 10 * log10(255*255 / MSE); return PSNR; } /* Compute Entropy - slow version */ double computeEntropy0(cv::Mat m) { assert(m.type() == CV_8UC1); const int size = 256; long histogram[size] = {0}; for (int i = 0; i < m.rows; ++i) { for(int j = 0; j < m.cols; ++j) { histogram[m.at(i, j)] += 1; } } double entropy = 0.0; for (int i=0; i 0) { const double pi = (histogram[i]) / (double)(m.rows*m.cols); entropy -= pi * log2(pi); } } return entropy; } /* Compute Entropy - fast version */ double computeEntropy1(cv::Mat m) { assert(m.type() == CV_8UC1); const int size = 256; long histogram[size] = {0}; int rows = m.rows; int cols = m.cols; if (m.isContinuous()) { cols *= rows; rows = 1; } for (int i = 0; i < rows; ++i) { const unsigned char *p = m.ptr(i); for(int j = 0; j < cols; ++j) { histogram[p[j]] += 1; } } const double norm = 1. /(double)(rows*cols); double entropy = 0.0; for (int i=0; i 0) { const double pi = (histogram[i]) * norm; entropy -= pi * log2(pi); } } return entropy; } /* Compute Entropy - with matrix operations */ double computeEntropy2(cv::Mat m) { const int histSize = 256; const float range[] = { 0, 256 } ; cv::Mat hist; const float* histRange = { range }; const bool uniform = true; const bool accumulate = false; cv::calcHist(&m, 1, 0, cv::Mat(), hist, 1, &histSize, &histRange, uniform, accumulate); //hist is a matrix of size 256x1 (i.e., a column vector) and of type CV_32F hist /= m.rows*m.cols; cv::Mat logP; cv::log(hist, logP); cv::divide(logP, log(2), logP); const double entropy = - (cv::sum(hist.mul(logP)))[0]; return entropy; } /* Compute error image - slow version */ void computeErrorImage0(const cv::Mat &im, const cv::Mat &imC, cv::Mat &imErr) { assert(im.cols == imC.cols && im.rows == imC.rows && im.type() == CV_8UC1 && imC.type() == CV_8UC1); imErr = cv::Mat(imC.rows, imC.cols, CV_8UC1); for(int y = 0; y < imErr.rows; ++y) { for(int x = 0; x < imErr.cols; ++x) { imErr.at(y, x) = std::min(255, std::max(0, im.at(y, x)-imC.at(y, x)+128)); } } } /* Compute Error Image - faster version */ void computeErrorImage1(const cv::Mat &im, const cv::Mat &imC, cv::Mat &imErr) { assert(im.cols == imC.cols && im.rows == imC.rows && im.type() == CV_8UC1 && imC.type() == CV_8UC1); imErr.create(im.rows, im.cols, CV_8UC1); //nothing done if matrix already has the desired size int rows = im.rows; int cols = im.cols; if (im.isContinuous() && imC.isContinuous() && imErr.isContinuous()) { cols *= rows; rows = 1; } for (int y = 0; y < rows; ++y) { const unsigned char *p1 = im.ptr(y); const unsigned char *p2 = imC.ptr(y); unsigned char *pe = imErr.ptr(y); for(int x = 0; x < cols; ++x) { pe[x] = cv::saturate_cast(p1[x]-p2[x] + 128); } } } void computeErrorImage2(const cv::Mat &im, const cv::Mat &imC, cv::Mat &imErr) { assert(im.cols == imC.cols && im.rows == imC.rows && im.type() == CV_8UC1 && imC.type() == CV_8UC1); imErr = cv::min(255, cv::max(0, im-imC+128)); } void computeErrorImage3(const cv::Mat &im, const cv::Mat &imC, cv::Mat &imErr) { assert(im.cols == imC.cols && im.rows == imC.rows && im.type() == CV_8UC1 && imC.type() == CV_8UC1); imErr = im-imC+128; } int main(int argc, char **argv) { if(argc != 3) { std::cerr<<"Usage: "< previousFrames; unsigned long frameNumber = 0; for ( ; ; ) { cv::Mat frameBGR; cap >> frameBGR; if (frameBGR.empty()) { break; } //convert from BGR to YCrCb colorspace cv::Mat frameYCbCr; cv::cvtColor(frameBGR, frameYCbCr, CV_BGR2YCrCb); //get Y plane std::vector planes(3); cv::split(frameYCbCr, planes); const cv::Mat Y = planes[0]; if (previousFrames.size() >= deltaT) { cv::Mat prevY = previousFrames.front(); previousFrames.pop(); const double MSE = computeMSE1(Y, prevY); const double PSNR = computePSNR(MSE); const double ENT = computeEntropy1(Y); std::cout<