写在前面
在实际应用中,效率是不得不考虑的问题。上一篇博客介绍了均值滤波原理,这一篇就写用积分图实现的快速均值滤波吧。
还是贴一下常规与快速的效率对比吧:
下图是常规均值滤波处理一张分辨率为485*528图像的时间(模板15*15):
下图是积分图快速均值滤波处理的时间(模板15*15):
可以说加速后,速度提升很多。而且最重要的是,用积分图的快速均值滤波受模板变化的影响不大!!!为深刻地体现这一特点,我把模板设置为151*151,可想而知,常规方法的计算量是多么的大!这个时候快速均值滤波的优势是多么的明显!请看下图:
下图是常规均值滤波处理一张分辨率为485*528图像的时间(模板151*151):
下图是积分图快速均值滤波处理的时间(模板151*151):
常规方法都达到了1738ms,快速方法还是3ms左右,惊艳有木有!!!
当然再加上sse优化、openmp+x86循环展开,速度到1ms都有可能吧。可惜我不会,以后再学sse优化吧。
积分图原理
积分图法也分为常规方法和改进方法。先介绍最原始的常规方法。
定义
积分图任意一个位置的值等于原图中该位置左上角所有值的和,即对于R*C图像矩阵的积分图可以定义计算:
举个栗子:
1 | 2 | 3 | 2 | 4 |
0 | 5 | 1 | 7 | 2 |
3 | 1 | 5 | 9 | 8 |
5 | 2 | 6 | 2 | 1 |
1 | 0 | 8 | 5 | 4 |
1 | 3 | 6 | 8 | 12 |
1 | 8 | 12 | 21 | 27 |
4 | 12 | 21 | 39 | 53 |
9 | 19 | 34 | 54 | 69 |
10 | 20 | 43 | 68 | 87 |
而且,积分图只需遍历一次图像即可有效地计算出来,因为积分图每一点(x,y)的值是:
所以,一旦积分图计算完毕,对任意矩形区域的和的计算就可以在常数时间内完成。如下图中,阴影矩形区域的和为:
举个栗子:
1 | 2 | 3 | 2 | 4 |
0 | 5 | 1 | 7 | 2 |
3 | 1 | 5 | 9 | 8 |
5 | 2 | 6 | 2 | 1 |
1 | 0 | 8 | 5 | 4 |
1 | 3 | 6 | 8 | 12 |
1 | 8 | 12 | 21 | 27 |
4 | 12 | 21 | 39 | 53 |
9 | 19 | 34 | 54 | 69 |
10 | 20 | 43 | 68 | 87 |
要求中间红色区域的和,只需对积分图进行如下计算:
也就是:5+1+7+1+5+9+2+6+2 = 54+1-8-9.
应用
积分图可用于快速计算Haar特征、二值化、均值滤波等等一切符合积分图法这一思想的算法。
由于均值滤波原理本质上是计算任意点的邻域内的平均值,而平均值由求和除以面积得到。所以无论如何改变模板大小,都可以利用积分图快速计算出每个点邻域内的和。
改进方法
对于积分图法的改进,有一位图像算法优化大佬给出了办法,图文并茂:13行代码实现最快速最高效的积分图像算法。思想就是不需要由三个位置的积分计算,只需由上方的Inegral(i-1,j)加上当前行的累加和即可,这样又减少了一次运算。这种改进办法比上面常规的方法快1.5倍左右!
代码
代码里包含了两种积分图的计算方法(常规和改进),暂时支持灰度图,RGB图就是每个通道分别计算就行。
#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgproc.hpp>
#include <ctime>
//
//积分图-常规方法
//由三个由三个位置的积分计算出来
//对于W*H图像:3*(W-1)*(H-1)次加减法
//
void Im_integral(cv::Mat& src,cv::Mat& dst){
int nr = src.rows;
int nc = src.cols;
dst = cv::Mat::zeros(nr + 1, nc + 1, CV_64F);
for (int i = 1; i < dst.rows; ++i){
for (int j = 1; j < dst.cols; ++j){
double top_left = dst.at<double>(i - 1, j - 1);
double top_right = dst.at<double>(i - 1, j);
double buttom_left = dst.at<double>(i, j - 1);
int buttom_right = src.at<uchar>(i-1, j-1);
dst.at<double>(i, j) = buttom_right + buttom_left + top_right - top_left;
}
}
}
//
//积分图-优化方法
//由上方negral(i-1,j)加上当前行的和即可
//对于W*H图像:2*(W-1)*(H-1)次加减法
//比常规方法快1.5倍左右
/
void Fast_integral(cv::Mat& src, cv::Mat& dst){
int nr = src.rows;
int nc = src.cols;
int sum_r = 0;
dst = cv::Mat::zeros(nr + 1, nc + 1, CV_64F);
for (int i = 1; i < dst.rows; ++i){
for (int j = 1, sum_r = 0; j < dst.cols; ++j){
//行累加,因为积分图相当于在原图上方加一行,左边加一列,所以积分图的(1,1)对应原图(0,0),(i,j)对应(i-1,j-1)
sum_r = src.at<uchar>(i-1 , j-1) + sum_r;
dst.at<double>(i, j) = dst.at<double>(i-1, j)+sum_r;
}
}
}
//积分图快速均值滤波
void Fast_MeanFilter(cv::Mat& src, cv::Mat& dst, cv::Size wsize){
//图像边界扩充
if (wsize.height % 2 == 0 || wsize.width % 2 == 0){
fprintf(stderr, "Please enter odd size!");
exit(-1);
}
int hh = (wsize.height - 1) / 2;
int hw = (wsize.width - 1) / 2;
cv::Mat Newsrc;
cv::copyMakeBorder(src, Newsrc, hh, hh, hw, hw, cv::BORDER_REFLECT_101);//以边缘为轴,对称
dst = cv::Mat::zeros(src.size(), src.type());
//计算积分图
cv::Mat inte;
Fast_integral(Newsrc, inte);
//均值滤波
double mean = 0;
for (int i = hh+1; i < src.rows + hh + 1;++i){ //积分图图像比原图(边界扩充后的)多一行和一列
for (int j = hw+1; j < src.cols + hw + 1; ++j){
double top_left = inte.at<double>(i - hh - 1, j - hw-1);
double top_right = inte.at<double>(i-hh-1,j+hw);
double buttom_left = inte.at<double>(i + hh, j - hw- 1);
double buttom_right = inte.at<double>(i+hh,j+hw);
mean = (buttom_right - top_right - buttom_left + top_left) / wsize.area();
//一定要进行判断和数据类型转换
if (mean < 0)
mean = 0;
else if (mean>255)
mean = 255;
dst.at<uchar>(i - hh - 1, j - hw - 1) = static_cast<uchar> (mean);
}
}
}
int main(){
cv::Mat src = cv::imread("I:\\Learning-and-Practice\\2019Change\\Image process algorithm\\Img\\Fig0334(a)(hubble-original).tif");
if (src.empty()){
return -1;
}
if (src.channels()>1)
cvtColor(src, src, CV_RGB2GRAY);
cv::Mat dst;
double t2 = (double)cv::getTickCount(); //测时间
Fast_MeanFilter(src, dst, cv::Size(151,151));//均值滤波
t2 = (double)cv::getTickCount() - t2;
double time2 = (t2 *1000.) / ((double)cv::getTickFrequency());
std::cout << "FASTmy_process=" << time2 << " ms. " << std::endl << std::endl;
cv::namedWindow("src");
cv::imshow("src", src);
cv::namedWindow("dst");
cv::imshow("dst", dst);
cv::waitKey(0);
}
效果
原图 快速积分图均值滤波(11*11) 经典均值滤波(11*11)
可见,效果上是一样的,但是效率大不一样,开篇已经介绍过了。