数字图像处理中经常会用到均值滤波，快速的均值方式通常会使用行求和，得到中间图像然后进行列求和。最后除以模版的面积。
基于积分图的均值滤波可以在O(1)时间内求出任意半径的均值滤波结果。前提是求出积分图。
先贴结果，如果下图积分图，中间图片是模版，则原图模版内的像素和为
I1, I2, I3, I4为中间模板的四个顶点，对应的是积分图中的像素值，中心点到矩形模版四条边的距离为r, 则均值滤波后中间像素值为:
plot by ItchyHacker.jpeg
下面看看公式是如何得出的:
I1 的像素值是原图A矩形内像素值和，I2 的像素值是原图A + C矩形内像素值和, I3 的像素值是原图A+B矩形内像素值和, I4 的像素值是原图A+B+C+D矩形内像素值和,
因此只要有积分图之后，任意半径的的均值滤波结果只需要取四个数加减然后除以均值滤波模版的面积便可。
OpenCV内置积分算法，下面贴一下我自己实现的积分算法和内置积分算法对比，还有基于积分图的均值滤波和两次求和遍历取平均的结果。程序边缘区域没有处理，只是简单的测试程序。
//
// 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){
// 我自己实现的积分图算法
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){
// 使用积分图实现的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);
}
结果:
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
