论文下载地址：http://research.microsoft.com/en-us/um/people/jiansun/papers/GuidedFilter_ECCV10.pdf
本文主要介绍导向滤波，但是在网上看这算法还能去雾，不知道是具体是怎么利用导向滤波实现去雾的，希望过来人指点迷津，这块主要是重写了导向滤波应用于彩色图像的部分代码，希望与大家共同交流。
论文主要如下：
Kaiming He, Jian Sun, Xiaoou Tang. Single Image Haze Removal Using Dark Channel Prior
大致内容是提出了一个叫做暗原色先验的东西来对有雾图像进行处理，十分巧妙，有兴趣者可以看看。这里使用OpenCV实现文中的去雾算法，然而论文提到的soft matting未在本程序中实现。
原理如下：
滤波效果：
单通道效果：
方法1效果：
方法2效果：
效果----为何要滤波：
guied filter滤波代码：使用了两种方法，代码来源后面参考文献中。我做了一些修改和比对工作。
// Guided Filter.cpp : 定义控制台应用程序的入口点。
//
#include "stdafx.h"
#include <iostream>
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#pragma comment(lib,"opencv_core2410d.lib")
#pragma comment(lib,"opencv_highgui2410d.lib")
#pragma comment(lib,"opencv_imgproc2410d.lib")
using namespace std;
using namespace cv;
Mat getimage(Mat &a)
{
int hei =a.rows;
int wid = a.cols;
Mat I(hei, wid, CV_64FC1);
//convert image depth to CV_64F
a.convertTo(I, CV_64FC1,1.0/255.0);
//normalize the pixel to 0~1
/*
for( int i = 0; i< hei; i++){
double *p = I.ptr<double>(i);
for( int j = 0; j< wid; j++){
p[j] = p[j]/255.0;
}
}
*/
return I;
}
Mat cumsum(Mat &imSrc, int rc)
{
if(!imSrc.data)
{
cout << "no data input!\n" << endl;
}
int hei = imSrc.rows;
int wid = imSrc.cols;
Mat imCum = imSrc.clone();
if( rc == 1)
{
for( int i =1;i < hei; i++)
{
for( int j = 0; j< wid; j++)
{
imCum.at<double>(i,j) += imCum.at<double>(i-1,j);
}
}
}
if( rc == 2)
{
for( int i =0;i < hei; i++)
{
for( int j = 1; j< wid; j++)
{
imCum.at<double>(i,j) += imCum.at<double>(i,j-1);
}
}
}
return imCum;
}
Mat boxfilter(Mat &imSrc, int r)
{
int hei = imSrc.rows;
int wid = imSrc.cols;
Mat imDst = Mat::zeros( hei, wid, CV_64FC1);
//imCum = cumsum(imSrc, 1);
Mat imCum = cumsum(imSrc,1);
//imDst(1:r+1, :) = imCum(1+r:2*r+1, :);
for( int i = 0; i<r+1; i++)
{
for( int j=0; j<wid; j++ )
{
imDst.at<double>(i,j) = imCum.at<double>(i+r,j);
}
}
//imDst(r+2:hei-r, :) = imCum(2*r+2:hei, :) - imCum(1:hei-2*r-1, :);
for( int i =r+1; i<hei-r;i++)
{
for( int j = 0; j<wid;j++)
{
imDst.at<double>(i,j) = imCum.at<double>(i+r,j)-imCum.at<double>(i-r-1,j);
}
}
//imDst(hei-r+1:hei, :) = repmat(imCum(hei, :), [r, 1]) - imCum(hei-2*r:hei-r-1, :);
for( int i = hei-r; i< hei; i++)
{
for( int j = 0; j< wid; j++)
{
imDst.at<double>(i,j) = imCum.at<double>(hei-1,j)-imCum.at<double>(i-r-1,j);
}
}
imCum = cumsum(imDst, 2);
//imDst(:, 1:r+1) = imCum(:, 1+r:2*r+1);
for( int i = 0; i<hei; i++)
{
for( int j=0; j<r+1; j++ )
{
imDst.at<double>(i,j) = imCum.at<double>(i,j+r);
}
}
//imDst(:, r+2:wid-r) = imCum(:, 2*r+2:wid) - imCum(:, 1:wid-2*r-1);
for( int i =0 ; i<hei;i++)
{
for( int j = r+1; j<wid-r ;j++ )
{
imDst.at<double>(i,j) = imCum.at<double>(i,j+r)-imCum.at<double>(i,j-r-1);
}
}
//imDst(:, wid-r+1:wid) = repmat(imCum(:, wid), [1, r]) - imCum(:, wid-2*r:wid-r-1);
for( int i = 0; i< hei; i++)
{
for( int j = wid-r; j<wid; j++)
{
imDst.at<double>(i,j) = imCum.at<double>(i,wid-1)-imCum.at<double>(i,j-r-1);
}
}
return imDst;
}
Mat guidedfilter( Mat &I, Mat &p, int r, double eps )
{
int hei = I.rows;
int wid = I.cols;
//N = boxfilter(ones(hei, wid), r);
Mat one = Mat::ones(hei, wid, CV_64FC1);
Mat N = boxfilter(one, r);
//mean_I = boxfilter(I, r) ./ N;
Mat mean_I(hei, wid, CV_64FC1);
divide(boxfilter(I, r), N, mean_I);
//mean_p = boxfilter(p, r) ./ N;
Mat mean_p(hei, wid, CV_64FC1);
divide(boxfilter(p, r), N, mean_p);
//mean_Ip = boxfilter(I.*p, r) ./ N;
Mat mul_Ip(hei, wid, CV_64FC1);
Mat mean_Ip(hei, wid, CV_64FC1);
multiply(I,p,mul_Ip);
divide(boxfilter(mul_Ip, r), N, mean_Ip);
//cov_Ip = mean_Ip - mean_I .* mean_p
//this is the covariance of (I, p) in each local patch.
Mat mul_mean_Ip(hei, wid, CV_64FC1);
Mat cov_Ip(hei, wid, CV_64FC1);
multiply(mean_I, mean_p, mul_mean_Ip);
subtract(mean_Ip, mul_mean_Ip, cov_Ip);
//mean_II = boxfilter(I.*I, r) ./ N;
Mat mul_II(hei, wid, CV_64FC1);
Mat mean_II(hei, wid, CV_64FC1);
multiply(I, I, mul_II);
divide(boxfilter(mul_II, r), N, mean_II);
//var_I = mean_II - mean_I .* mean_I;
Mat mul_mean_II(hei, wid, CV_64FC1);
Mat var_I(hei, wid, CV_64FC1);
multiply(mean_I, mean_I, mul_mean_II);
subtract(mean_II, mul_mean_II, var_I);
//a = cov_Ip ./ (var_I + eps);
Mat a(hei, wid, CV_64FC1);
for( int i = 0; i< hei; i++){
double *p = var_I.ptr<double>(i);
for( int j = 0; j< wid; j++){
p[j] = p[j] +eps;
}
}
divide(cov_Ip, var_I, a);
//b = mean_p - a .* mean_I;
Mat a_mean_I(hei ,wid, CV_64FC1);
Mat b(hei ,wid, CV_64FC1);
multiply(a, mean_I, a_mean_I);
subtract(mean_p, a_mean_I, b);
//mean_a = boxfilter(a, r) ./ N;
Mat mean_a(hei, wid, CV_64FC1);
divide(boxfilter(a, r), N, mean_a);
//mean_b = boxfilter(b, r) ./ N;
Mat mean_b(hei, wid, CV_64FC1);
divide(boxfilter(b, r), N, mean_b);
//q = mean_a .* I + mean_b;
Mat mean_a_I(hei, wid, CV_64FC1);
Mat q(hei, wid, CV_64FC1);
multiply(mean_a, I, mean_a_I);
add(mean_a_I, mean_b, q);
return q;
}
/*****************
http://research.microsoft.com/en-us/um/people/kahe/eccv10/
推酷上的一篇文章：
http://www.tuicool.com/articles/Mv2iiu
************************/
cv::Mat guidedFilter2(cv::Mat I, cv::Mat p, int r, double eps)
{
/*
% GUIDEDFILTER O(1) time implementation of guided filter.
%
% - guidance image: I (should be a gray-scale/single channel image)
% - filtering input image: p (should be a gray-scale/single channel image)
% - local window radius: r
% - regularization parameter: eps
*/
cv::Mat _I;
I.convertTo(_I, CV_64FC1);
I = _I;
cv::Mat _p;
p.convertTo(_p, CV_64FC1);
p = _p;
//[hei, wid] = size(I);
int hei = I.rows;
int wid = I.cols;
//N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.
cv::Mat N;
cv::boxFilter(cv::Mat::ones(hei, wid, I.type()), N, CV_64FC1, cv::Size(r, r));
//mean_I = boxfilter(I, r) ./ N;
cv::Mat mean_I;
cv::boxFilter(I, mean_I, CV_64FC1, cv::Size(r, r));
//mean_p = boxfilter(p, r) ./ N;
cv::Mat mean_p;
cv::boxFilter(p, mean_p, CV_64FC1, cv::Size(r, r));
//mean_Ip = boxfilter(I.*p, r) ./ N;
cv::Mat mean_Ip;
cv::boxFilter(I.mul(p), mean_Ip, CV_64FC1, cv::Size(r, r));
//cov_Ip = mean_Ip - mean_I .* mean_p; % this is the covariance of (I, p) in each local patch.
cv::Mat cov_Ip = mean_Ip - mean_I.mul(mean_p);
//mean_II = boxfilter(I.*I, r) ./ N;
cv::Mat mean_II;
cv::boxFilter(I.mul(I), mean_II, CV_64FC1, cv::Size(r, r));
//var_I = mean_II - mean_I .* mean_I;
cv::Mat var_I = mean_II - mean_I.mul(mean_I);
//a = cov_Ip ./ (var_I + eps); % Eqn. (5) in the paper;
cv::Mat a = cov_Ip/(var_I + eps);
//b = mean_p - a .* mean_I; % Eqn. (6) in the paper;
cv::Mat b = mean_p - a.mul(mean_I);
//mean_a = boxfilter(a, r) ./ N;
cv::Mat mean_a;
cv::boxFilter(a, mean_a, CV_64FC1, cv::Size(r, r));
mean_a = mean_a/N;
//mean_b = boxfilter(b, r) ./ N;
cv::Mat mean_b;
cv::boxFilter(b, mean_b, CV_64FC1, cv::Size(r, r));
mean_b = mean_b/N;
//q = mean_a .* I + mean_b; % Eqn. (8) in the paper;
cv::Mat q = mean_a.mul(I) + mean_b;
return q;
}
int _tmain(int argc, _TCHAR* argv[])
{
int r = 4;
double eps = 0.01;
string image_name ;
cout<<"input name:"<<endl;
cin>>image_name;
/*
CV_LOAD_IMAGE_ANYDEPTH - If set, return 16-bit/32-bit image when the input has the corresponding depth,
otherwise convert it to 8-bit.
CV_LOAD_IMAGE_COLOR - If set, always convert image to the color one
CV_LOAD_IMAGE_GRAYSCALE - If set, always convert image to the grayscale one
>0 Return a 3-channel color image.
Note:
In the current implementation the alpha channel, if any, is stripped from the output image.
Use negative value if you need the alpha channel.
=0 Return a grayscale image.
<0 Return the loaded image as is (with alpha channel).
*/
Mat image_src = imread(image_name,CV_LOAD_IMAGE_COLOR);
Mat image_gray(image_src.size(),CV_8UC1);
cvtColor(image_src,image_gray,CV_BGR2GRAY);
vector<Mat> bgr_src,bgr_dst;
split(image_src,bgr_src);//分解每个通道
Mat dst_color;
double time;
time = (double)getTickCount();
for(int i=0;i<3;i++)
{
Mat I = getimage(bgr_src[i]);
Mat p = I.clone();
Mat q = guidedfilter(I, p, r, eps);
//string number ;
//sprintf((char *)number.c_str(),"%d",i);
//imshow(number,q);
//imshow("方法1：", q);
bgr_dst.push_back(q);
//cv::merge(q,dst_color);
}
merge(bgr_dst,dst_color);
//imwrite("filtered.bmp", q*255);
time = 1000*((double)getTickCount() - time)/getTickFrequency();
cout <<endl<<"Time of guided filter for runs: " << time << " milliseconds."<< endl;
imshow("原图像的灰度图", image_gray);
imshow("方法1：", dst_color);
imwrite("result.jpg",dst_color*255);
double time2 = 0;
time2 = (double)getTickCount();
Mat I = getimage(image_gray);
Mat p = I.clone();
//int r = 8;
//double eps = 0.04;
//Mat q = guidedfilter(I, p, r, eps);
//imwrite("filtered.bmp", q*255);
//*/
/*imshow("原图像的灰度图", image_gray);
imshow("方法1：", q);*/
imshow("方法2：",guidedFilter2(I, p, r, eps));
time2 = 1000*((double)getTickCount() - time2)/getTickFrequency();
cout <<endl<<"Time of guided filter2 for runs: " << time2 << " milliseconds."<< endl;
waitKey(0);
return 0;
}
下面的代码还没有真正的调试，只是找到了，先放在这里，后面有空再看看研究一下。
去雾代码1：
#include<iostream.h>
#include<cv.h>
#include<highgui.h>
char tbarname1[] = "调节block";
//定义两个滑动条,用于调节参数
char tbarname2[] = "调节w";
//w是为了保留一部分的雾
int block=5;
int w1=80;
double w;
IplImage *src=NULL;
IplImage *dst=NULL;
//定义去雾函数如下
IplImage *quw(IplImage *src,int block,double w)
{
//图像分别有三个颜色通道
IplImage *dst1=NULL;
IplImage *dst2=NULL;
IplImage *dst3=NULL;
IplImage *imgroi1;
//dst1的ROI
IplImage *imgroi2;
//dst2的ROI
IplImage *imgroi3;
//dst3的ROI
IplImage *roidark;
//dark channel的ROI
IplImage *dark_channel=NULL;
//暗原色先验的指针
IplImage *toushelv=NULL;
//透射率
//去雾算法运算后的三个通道
IplImage *j1=NULL;
IplImage *j2=NULL;
IplImage *j3=NULL;
//去雾后的图像,三通道合并成
IplImage *dst=NULL;
//源图像ROI位置以及大小
CvRect ROI_rect;
//分离的三个通道
dst1=cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,1);
dst2=cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,1);
dst3=cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,1);
//为各个ROI分配内存
imgroi1=cvCreateImage(cvSize(block,block),IPL_DEPTH_8U,1);
imgroi2=cvCreateImage(cvSize(block,block),IPL_DEPTH_8U,1);
imgroi3=cvCreateImage(cvSize(block,block),IPL_DEPTH_8U,1);
roidark=cvCreateImage(cvSize(block,block),IPL_DEPTH_8U,1);
//为j1 j2 j3分配大小
j1=cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,1);
j2=cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,1);
j3=cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,1);
//为暗原色先验指针分配大小
dark_channel=cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,1);
//为透射率指针分配大小
toushelv=cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,1);
//dst分配大小
dst=cvCreateImage(cvSize(src->width,src->height),IPL_DEPTH_8U,3);
//将原彩色图像分离成三通道
cvSplit(src,dst1,dst2,dst3,NULL);
//求暗原色
ROI_rect.width=block;
ROI_rect.height=block;
ROI_rect.x=0;
ROI_rect.y=0;
int i;
int j;
double min1=0;
double max1=0;
double min2=0;
double max2=0;
double min3=0;
double max3=0;
double min=0;
CvScalar value;
for(i=0;i<src->width/block;i++)
{ for(j=0;j<src->height/block;j++)
{
//分别计算三个通道内ROI的最小值
cvSetImageROI(dst1,ROI_rect);
cvCopy(dst1,imgroi1,NULL);
cvMinMaxLoc(imgroi1,&min1,&max1,NULL,NULL);
cvSetImageROI(dst2,ROI_rect);
cvCopy(dst2,imgroi2,NULL);
cvMinMaxLoc(imgroi2,&min2,&max2,NULL,NULL);
cvSetImageROI(dst3,ROI_rect);
cvCopy(dst3,imgroi3,NULL);
cvMinMaxLoc(imgroi3,&min3,&max3,NULL,NULL);
//求三个通道内最小值的最小值
if(min1<min2)
min=min1;
else
min=min2;
if(min>min3)
min=min3;//min为这个ROI中暗原色
value=cvScalar(min,min,min,min);//min放在value中
//min赋予dark_channel中相应的ROI
cvSetImageROI(dark_channel,ROI_rect);
cvSet(roidark,value,NULL);
cvCopy(roidark,dark_channel,NULL);
//释放各个ROI
cvResetImageROI(dst1);
cvResetImageROI(dst2);
cvResetImageROI(dst3);
cvResetImageROI(dark_channel);
//转入下一个ROI
ROI_rect.x=block*i;
ROI_rect.y=block*j;
}
}
//保存暗原色先验的图像
cvSaveImage("f:/dark_channel_prior.jpg",dark_channel);
//利用得到的暗原色先验dark_channel_prior.jpg求大气光强
double min_dark;
double max_dark;
CvPoint min_loc;
CvPoint max_loc;//max_loc是暗原色先验最亮一小块的原坐标
cvMinMaxLoc(dark_channel,&min_dark,&max_dark,&min_loc,&max_loc,NULL);
cout<<max_loc.x<<" "<<max_loc.y<<endl;
ROI_rect.x=max_loc.x;
ROI_rect.y=max_loc.y;
double A_dst1;//定义大气光成分的估计值
double dst1_min;
double A_dst2;
double dst2_min;
double A_dst3;
double dst3_min;
cvSetImageROI(dst1,ROI_rect);
//按照论文方法求大气光强估计值
cvCopy(dst1,imgroi1,NULL);
cvMinMaxLoc(imgroi1,&dst1_min,&A_dst1,NULL,NULL);
cvSetImageROI(dst2,ROI_rect);
cvCopy(dst2,imgroi2,NULL);
cvMinMaxLoc(imgroi2,&dst2_min,&A_dst2,NULL,NULL);
cvSetImageROI(dst3,ROI_rect);
cvCopy(dst3,imgroi3,NULL);
cvMinMaxLoc(imgroi3,&dst3_min,&A_dst3,NULL,NULL);
cout<<A_dst1<<" "<<A_dst2<<" "<<A_dst3<<endl;//这三值为大气光强度估计值
//求透射率
int k;
int l;
CvScalar m;
CvScalar n;//暗原色先验各元素值
for(k=0;k<src->height;k++)
{
for(l=0;l<src->width;l++)
{
m=cvGet2D(dark_channel,k,l);
n=cvScalar(255-w*m.val[0]);
//w目的是保留一部分的雾,使图像看起来真实些
cvSet2D(toushelv,k,l,n);
}
}
cvSaveImage("f:/toushelv.jpg",toushelv);
//求无雾图像
int p,q;
double tx;
double jj1,jj2,jj3;
CvScalar ix,jx;
for(p=0;p<src->height;p++)
{
for(q=0;q<src->width;q++)
{
tx=cvGetReal2D(toushelv,p,q);
tx=tx/255;
if(tx<0.1)
tx=0.1;
ix=cvGet2D(src,p,q);
jj1=(ix.val[0]-A_dst1)/tx+A_dst1;//根据雾产生模型运算,还原出无雾图像
jj2=(ix.val[1]-A_dst2)/tx+A_dst2;
jj3=(ix.val[2]-A_dst3)/tx+A_dst3;
jx=cvScalar(jj1,jj2,jj3,0.0);
cvSet2D(dst,p,q,jx);
}
}
cvSaveImage("f:/removed_haze.jpg",dst);
//释放指针
cvReleaseImage(&dst1);
cvReleaseImage(&dst2);
cvReleaseImage(&dst3);
cvReleaseImage(&imgroi1);
cvReleaseImage(&imgroi2);
cvReleaseImage(&imgroi3);
cvReleaseImage(&roidark);
cvReleaseImage(&dark_channel);
cvReleaseImage(&toushelv);
cvReleaseImage(&j1);
cvReleaseImage(&j2);
cvReleaseImage(&j3);
return dst;
}
void on_trackbar1(int h)
{
dst=quw(src,block,w);
cvShowImage("目的图像",dst);
// cvWaitKey(0);
}
void on_trackbar2(int h)
{
w=(double)w1/100;
dst=quw(src,block,w);
cvShowImage("目的图像",dst);
// cvWaitKey(0);
}
//主函数如下
void main()
{
//打开图像
src=cvLoadImage("8.jpg",-1);
//创造窗口
cvNamedWindow("有雾图像",CV_WINDOW_AUTOSIZE);
cvShowImage("有雾图像",src);
cvNamedWindow("目的图像",CV_WINDOW_AUTOSIZE);
cvCreateTrackbar(tbarname1, "目的图像", &block, 15, on_trackbar1);
cvCreateTrackbar(tbarname2, "目的图像", &w1, 100, on_trackbar2);
cvWaitKey(0);
cvReleaseImage(&src);
cvReleaseImage(&dst);
}
去雾matlab代码：
去雾代码2：
参考文献：
http://www.tuicool.com/articles/Mv2iiu
http://blog.csdn.net/holybang/article/details/28093305
http://www.tuicool.com/articles/MJZr2e
http://blog.sina.com.cn/s/blog_4d8730df0100m8lz.html