數(shù)字圖像處理中經(jīng)常會用到均值濾波芜茵，快速的均值方式通常會使用行求和，得到中間圖像然后進行列求和倡蝙。最后除以模版的面積九串。
基于積分圖的均值濾波可以在O(1)時間內(nèi)求出任意半徑的均值濾波結(jié)果。前提是求出積分圖寺鸥。
先貼結(jié)果猪钮，如果下圖積分圖，中間圖片是模版胆建，則原圖模版內(nèi)的像素和為
I1, I2, I3, I4為中間模板的四個頂點烤低，對應(yīng)的是積分圖中的像素值，中心點到矩形模版四條邊的距離為r, 則均值濾波后中間像素值為:
plot by ItchyHacker.jpeg
下面看看公式是如何得出的:
I1 的像素值是原圖A矩形內(nèi)像素值和笆载，I2 的像素值是原圖A + C矩形內(nèi)像素值和, I3 的像素值是原圖A+B矩形內(nèi)像素值和, I4 的像素值是原圖A+B+C+D矩形內(nèi)像素值和,
因此只要有積分圖之后扑馁，任意半徑的的均值濾波結(jié)果只需要取四個數(shù)加減然后除以均值濾波模版的面積便可。
OpenCV內(nèi)置積分算法凉驻，下面貼一下我自己實現(xiàn)的積分算法和內(nèi)置積分算法對比腻要，還有基于積分圖的均值濾波和兩次求和遍歷取平均的結(jié)果。程序邊緣區(qū)域沒有處理涝登，只是簡單的測試程序雄家。
//
// integralCompare.cpp
// FaceBeauty
//
// Created by yuhua.cheng on 2018/11/6.
// Copyright ? 2018 yuhua.cheng. All rights reserved.
//
#include <iostream>
#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
void myIntegral(const cv::Mat& img, cv::Mat& integral){
// 我自己實現(xiàn)的積分圖算法
int height = img.rows;
int width = img.cols;
cv::Mat tempImg;
img.convertTo(tempImg, CV_32FC3, 1/255.0);
integral = tempImg.clone();
// 第一行積分圖值計算
float* imgPtr = tempImg.ptr<float>(0);
float* interPrevPtr = integral.ptr<float>(0);
float* interPtr = integral.ptr<float>(0);
interPtr += 3;
for(int j = 1; j < width; j++){
interPtr[0] = interPrevPtr[0] + imgPtr[0];
interPtr[1] = interPrevPtr[1] + imgPtr[1];
interPtr[2] = interPrevPtr[2] + imgPtr[2];
imgPtr += 3;
interPrevPtr += 3;
interPtr += 3;
}
// 第一列積分圖值計算
for(int i = 1; i < height; i++){
interPrevPtr = integral.ptr<float>(i-1);
interPtr = integral.ptr<float>(i);
imgPtr = tempImg.ptr<float>(i);
interPtr[0] = interPrevPtr[0] + imgPtr[0];
interPtr[1] = interPrevPtr[1] + imgPtr[1];
interPtr[2] = interPrevPtr[2] + imgPtr[2];
}
// 利用第一行和第一列信息計算積分圖
for(int i = 1; i < height; i++){
interPrevPtr = integral.ptr<float>(i-1)+3;
interPtr = integral.ptr<float>(i)+3;
imgPtr = tempImg.ptr<float>(i)+3; // 這里可以改進
float rowSum[3] = {0};
for(int j = 1; j < width; j++){
rowSum[0] += (imgPtr-3)[0];
rowSum[1] += (imgPtr-3)[1];
rowSum[2] += (imgPtr-3)[2];
interPtr[0] = interPrevPtr[0] + imgPtr[0] + rowSum[0];
interPtr[1] = interPrevPtr[1] + imgPtr[1] + rowSum[1];
interPtr[2] = interPrevPtr[2] + imgPtr[2] + rowSum[2];
interPrevPtr += 3;
interPtr += 3;
imgPtr += 3;
}
}
}
void integralBoxBlur(const cv::Mat& img, cv::Mat& dst, int radius){
// 使用積分圖實現(xiàn)的boxBlur算法
int r = radius;
cv::Mat tempImg;
img.convertTo(tempImg, CV_32FC3, 1/255.0);
cv::Mat integral;
myIntegral(img, integral);
int weight = (2*r+1)*(2*r+1);
for(int row = r; row < img.rows - r; row++){
float* intePtr1 = integral.ptr<float>(row-r);
float* intePtr2 = integral.ptr<float>(row+r);
float* imgPtr = tempImg.ptr<float>(row);
intePtr1 += 3*r;
intePtr2 += 3*r;
imgPtr += 3*r;
for(int col = r; col < img.cols - r; col++){
imgPtr[0] = ((intePtr2+3*r)[0] + (intePtr1-3*r)[0] - (intePtr1+3*r)[0] - (intePtr2-3*r)[0])/weight;
imgPtr[1] = ((intePtr2+3*r)[1] + (intePtr1-3*r)[1] - (intePtr1+3*r)[1] - (intePtr2-3*r)[1])/weight;
imgPtr[2] = ((intePtr2+3*r)[2] + (intePtr1-3*r)[2] - (intePtr1+3*r)[2] - (intePtr2-3*r)[2])/weight;
intePtr1 += 3;
intePtr2 += 3;
imgPtr += 3;
}
}
tempImg.convertTo(dst, CV_8UC3, 255);
}
void sumBoxBlur(const cv::Mat& img, cv::Mat& dst, int radius){
// 普通求和平均的計算方式
int r = radius;
int weight = (2*r+1)*(2*r+1);
int height = img.rows;
int width = img.cols;
cv::Mat tempImg;
img.convertTo(tempImg, CV_32FC3, 1/255.0);
cv::Mat rowSum;
rowSum = tempImg.clone();
cv::Mat colSum;
colSum = tempImg.clone();
dst = tempImg.clone();
// 行求和
for(int i = 0; i < height; i++){
float *imgPtr = tempImg.ptr<float>(i);
float *colPtr = colSum.ptr<float>(i);
for(int j = 0; j < width; j++){
if(j < r){
imgPtr += 3;
colPtr += 3;
}else{
for(int bias = -r; bias <= r; bias++){
colPtr[0] += (imgPtr+3*bias)[0];
colPtr[1] += (imgPtr+3*bias)[1];
colPtr[2] += (imgPtr+3*bias)[2];
}
imgPtr +=3;
colPtr += 3;
}
}
}
cv::imshow("colSum",colSum);
// 列求和
for(int j = 0; j < width; j++){
for(int i = 0; i < height; i++){
if(i < r){
continue;
}else{
float *rowPtr = rowSum.ptr<float>(i)+3*j;
for(int bias = -r; bias <= r; bias++){
rowPtr[0] += (colSum.ptr<float>(i+bias)+3*j)[0];
rowPtr[1] += (colSum.ptr<float>(i+bias)+3*j)[1];
rowPtr[2] += (colSum.ptr<float>(i+bias)+3*j)[2];
}
}
}
}
dst = rowSum / weight;
dst.convertTo(dst, CV_8UC3, 255);
}
int main(){
cv::Mat img = cv::imread("/Users/yuhua.cheng/Documents/dataset/Face/10.jpg");
int height = img.rows;
int width = img.cols;
cv::Mat result1, result2, result3;
cv::Mat tempImg;
img.convertTo(tempImg, CV_32FC3, 1/255.0);
cv::Mat integral;
double t0, t1;
t0 = cv::getTickCount();
cv::integral(tempImg, integral, -1);
t1 = cv::getTickCount();
std::cout << "Opencv integral time: " << (t1 - t0) / cv::getTickFrequency() << std::endl;
t0 = cv::getTickCount();
myIntegral(img, integral);
t1 = cv::getTickCount();
std::cout << "My integral time: " << (t1 - t0) / cv::getTickFrequency() << std::endl;
t0 = cv::getTickCount();
integralBoxBlur(img, result1, 11);
t1 = cv::getTickCount();
std::cout << "Integral box Blur Time: " << (t1 - t0) / cv::getTickFrequency() << std::endl;
cv::imshow("my BoxBlur", result1);
t0 = cv::getTickCount();
cv::boxFilter(img, result2, -1, cv::Size(23,23));
t1 = cv::getTickCount();
std::cout << "Opencv box Filter Time: " << (t1 - t0) / cv::getTickFrequency() << std::endl;
cv::imshow("opencv BoxFilter", result2);
t0 = cv::getTickCount();
sumBoxBlur(img, result3, 11);
t1 = cv::getTickCount();
std::cout << "sumBoxFilter Time: " << (t1 - t0) / cv::getTickFrequency() << std::endl;
cv::imshow("sum BoxFilter", result3);
cv::waitKey(0);
}
結(jié)果:
Opencv integral time: 0.0155424
My integral time: 0.0254387
Integral box Blur Time: 0.0337875
Opencv box Filter Time: 0.00485557
sumBoxFilter Time: 0.319712
Reference
- https://blog.csdn.net/lijianlarry/article/details/78678297
- https://blog.csdn.net/u010839382/article/details/48241929
- Non-local Means Algorithm Introduction and Implementation: http://www.ipol.im/pub/art/2011/bcm_nlm/
- 積分圖的并行算法:https://blog.csdn.net/10km/article/details/51610735
- 快速均值濾波算法 : https://www.opengl.org/discussion_boards/showthread.php/167435-fastest-wide-width-box-filter
