OpenCV导向滤波实现（Guided Filter）代码，以及使用颜色先验算法去雾

本文主要介绍导向滤波，但是在网上看这算法还能去雾，不知道是具体是怎么利用导向滤波实现去雾的，希望过来人指点迷津，这快主要是重写了导向滤波应用与彩色图像的部分代码，希望与大家共同交流。

论文主要如下：

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代码：

<span style="font-size:24px;">　　function q = guidedfilter(I, p, r, 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

　　[hei, wid] = size(I);
　　N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.

　　% imwrite(uint8(N), 'N.jpg');
　　% figure,imshow(N,[]),title('N');
　　

　　mean_I = boxfilter(I, r) ./ N;
　　mean_p = boxfilter(p, r) ./ N;
　　mean_Ip = boxfilter(I.*p, r) ./ N;
　　cov_Ip = mean_Ip - mean_I .* mean_p; % this is the covariance of (I, p) in each local patch.

　　mean_II = boxfilter(I.*I, r) ./ N;
　　var_I = mean_II - mean_I .* mean_I;

　　a = cov_Ip ./ (var_I + eps); % Eqn. (5) in the paper;
　　b = mean_p - a .* mean_I; % Eqn. (6) in the paper;

　　mean_a = boxfilter(a, r) ./ N;
　　mean_b = boxfilter(b, r) ./ N;

　　q = mean_a .* I + mean_b; % Eqn. (8) in the paper;
　　end</span>

去雾代码2：

#include "stdafx.h"
#include <opencv2\opencv.hpp>
#include "cv.h"
#include <cxcore.h>
#include "highgui.h"
#include <windows.h>
#include <math.h>
using namespace cv;
using namespace std; 

//求三个通道中最小值的最小值时调用的函数
double min(double b,double g, double r)
{
	double result = b;
	if(result>g)
		result = g;
	if(result>r)
		result = r;
	return result;
};
double max(double a,double b)
{
	double MAX;
	if (a<b)
		MAX = b;
	else
		MAX = a;
	return MAX;
};
double min2(double a,double b)//比较两个数值中的最小值并返回
{
	double MIN;
	if (a<b)
		MIN = a;
	else
		MIN = b;
	return MIN;
};
//这个函数相当于darkchannel的功能，但在padarray时。使用的是将边缘像素复制的方法，不是matlab的将边缘处镜像复制,计算出darkchannel后有计算了最大值A
double doDarkChannel(IplImage* in,int patchsize)
{

	int height,width,step,channels;//图像的宽，高，等信息，height，width是输入图像的尺寸，也是输出图像的尺寸，step是输出图像jout的（j对应matlab代码中的darkchannel的输出图像J）
	int i,j,k;//用于循环的变量
	uchar *data2;//输出的结果图像的指针
	height = in->height;//获取输入图像的宽高等信息
	width = in->width;
	int patch = patchsize/2;//图像要延拓的边缘的宽度
	IplImage* mout=cvCreateImage(cvSize(in->width+patchsize,in->height+patchsize),in->depth,in->nChannels);   //存放图像被镜像延拓后图像的空图像
	cvCopyMakeBorder(in,mout,cvPoint(patch,patch),IPL_BORDER_REPLICATE);//这个函数相当于padarray，mout中存放padarrry后的图像
	IplImage* jout = cvCreateImage(cvSize(in->width,in->height),in->depth,1);//darkchannel 的输出结果，J

	step = jout->widthStep/sizeof(uchar);//step是单通道输出图像jout的widthstep
	data2 = (uchar *)jout->imageData;//指向输出图像的数据头

	for(i=0;i<height;i++)
		{
       for(j=0;j<width;j++)
          {
				cvSetImageROI(mout, cvRect(j, i, patchsize, patchsize));//操作输入图像的（i，j）点处patchsize大小的图像块
				IplImage*  patch_out=cvCreateImage(cvSize(patchsize,patchsize),in->depth,in->nChannels);//存储三通道图像块的临时内存区，循环体里用到的内存区域再循环体里申请，在循环体里释放
				cvCopy(mout,patch_out);//将patchsize大小的图像块存入临时图像块patch_out
				cvResetImageROI(mout); //释放mout

				//以下内容是利用cnMinMaxloc分别计算三个通道中的最小值
				double MinValue;
				double MaxValue;
				double B_Min,G_Min,R_Min;
				CvPoint MinLocation;
				CvPoint MaxLocation;
				cvSetImageCOI(patch_out,1);
				cvMinMaxLoc(patch_out,& MinValue,& MaxValue,& MinLocation,& MaxLocation);
				B_Min = MinValue;
				cvSetImageCOI(patch_out,2);
				cvMinMaxLoc(patch_out,& MinValue,& MaxValue,& MinLocation,& MaxLocation);
				G_Min = MinValue;
				cvSetImageCOI(patch_out,3);
				cvMinMaxLoc(patch_out,& MinValue,& MaxValue,& MinLocation,& MaxLocation);
				R_Min = MinValue;
				int dark_point = (int)min(B_Min,G_Min,R_Min);
				//三个通道的最小值都已经被分别提取出来了
				data2[i*step+j] = dark_point;//step 是jout的step，是单通道的
				cvReleaseImage(&patch_out);
			};
	};
	double MinValue;
	double MaxValue;
	double B_Min,G_Min,R_Min;
	CvPoint MinLocation;
	CvPoint MaxLocation;
	cvSetImageCOI(jout,1);
	cvMinMaxLoc(jout,& MinValue,& MaxValue,& MinLocation,& MaxLocation);
	cvReleaseImage(&mout);
	cout<<"计算暗通道函数运行成功"<<"\n";
	return MaxValue;

};
//该函数的作用相当于matlab代码中求取三个通道中最小值，然后以最小值组成一幅灰度图
IplImage* doMinChannel(IplImage* in)
{
	IplImage* b = cvCreateImage(cvSize(in->width,in->height),in->depth,1);
	IplImage* g = cvCreateImage(cvSize(in->width,in->height),in->depth,1);
	IplImage* r = cvCreateImage(cvSize(in->width,in->height),in->depth,1);//创建保存读入图像三个通道的的内存区域
	IplImage* w = cvCreateImage(cvSize(in->width,in->height),in->depth,1);//创建保存输出图像的内存区域（三个通道中最小值组成的一幅灰度图）
	cvSetImageCOI(in,1);
	cvCopy(in,b);
	cvSetImageCOI(in,2);
	cvCopy(in,g);
	cvSetImageCOI(in,3);
	cvCopy(in,r);//将三个通道的的值分别存入r，g，b三块内存区域中

	//cvSplit(src,dst1,dst2,dst3,NULL);

	int height = in->height;//获取输入图像的宽高等信息
	int width = in->width;
	int i,j,k;//用于循环的变量
	uchar *data_w;
	uchar *data_b;
	uchar *data_g;
	uchar *data_r;
	int step = b->widthStep/sizeof(uchar);
	data_w = (uchar *)w->imageData;//指向输出图像的数据头
	data_b = (uchar *)b->imageData;//指向输出图像的数据头
	data_g = (uchar *)g->imageData;//指向输出图像的数据头
	data_r = (uchar *)r->imageData;//指向输出图像的数据头
	for(i=0;i<height;i++)
		{
       for(j=0;j<width;j++)
          {
			  double b,g,r;
			  int MIN;//b,g,r三个通道的最小值
			  b = data_b[i*step+j];
			  g = data_g[i*step+j];
			  r = data_r[i*step+j];
			  MIN = (int)min(b,g,r);
			  data_w[i*step+j] = MIN;
		  };
		};
	cout<<"计算三个通道最小值并组成一幅新灰度图的函数运行成功"<<"\n";//表示该函数运行成功
	return w;
}	;
IplImage* doCalculateV(IplImage* w,IplImage* diff,IplImage* smooth)
{
    IplImage* b = cvCreateImage(cvSize(w->width,w->height),w->depth,1);
	IplImage* v = cvCreateImage(cvSize(w->width,w->height),w->depth,1);
	int height = w->height;//获取输入图像的宽高等信息
	int width = w->width;
	int i,j,k;//用于循环的变量
	uchar *data_w;
	uchar *data_diff;
	uchar *data_v;
	uchar *data_b;
	uchar *data_smooth;
	int step = w->widthStep/sizeof(uchar);
	data_w = (uchar *)w->imageData;//指向输出图像的数据头
	data_diff = (uchar *)diff->imageData;//指向输出图像的数据头
	data_v = (uchar *)v->imageData;//指向输出图像的数据头
	data_b = (uchar *)b->imageData;//指向输出图像的数据头
	data_smooth = (uchar *)smooth->imageData;
	for(i=0;i<height;i++)
		{
       for(j=0;j<width;j++)
          {
			  double W;
			  double DIFF;
			  double B;
			  double SMOOTH;
			  double p = 0.78;//p = 0.78
			  double MIN,MAX;
			  W = data_w[i*step+j];
			  DIFF = data_diff[i*step+j];
			  SMOOTH = data_smooth[i*step+j];
			  B = W-DIFF;
			  MIN = min2(B,SMOOTH);
			  MAX = max(MIN,0);
			  data_v[i*step+j] = p*MAX;
		  };
		};
	cout<<"计算v函数运行成功"<<"\n";//表示该函数运行成功
	return v;
};
//计算最终的去雾图像的函数
IplImage* doFinally(IplImage* in,IplImage* v,double A)
{
	IplImage* b = cvCreateImage(cvSize(in->width,in->height),in->depth,1);
	IplImage* g = cvCreateImage(cvSize(in->width,in->height),in->depth,1);
	IplImage* r = cvCreateImage(cvSize(in->width,in->height),in->depth,1);
	IplImage* result = cvCreateImage(cvSize(in->width,in->height),in->depth,3);//创建存储输出图像的内存区域
	int height = in->height;//获取输入图像的宽高等信息
	int width = in->width;
	int i,j;//用于循环的变量
	cvSetImageCOI(in,1);
	cvCopy(in,b);
	cvSetImageCOI(in,2);
	cvCopy(in,g);
	cvSetImageCOI(in,3);
	cvCopy(in,r);//将三个通道的的值分别存入r，g，b三块内存区域中
	//cvSplit(in,b,g,r,NULL);//将图像拆分为三个通道
	uchar *data_b;
	uchar *data_g;
	uchar *data_r;
	uchar *data_v;
	int step = b->widthStep/sizeof(uchar);
	//data_w = (uchar *)w->imageData;//指向输出图像的数据头
	data_b = (uchar *)b->imageData;//指向蓝色通道的数据头
	data_g = (uchar *)g->imageData;//指向绿色通道的数据头
	data_r = (uchar *)r->imageData;//指向红色通道的数据头
	data_v = (uchar *)v->imageData;
	//计算蓝色通道的去雾结果
	for(i=0;i<height;i++)
		{
       for(j=0;j<width;j++)
          {
			  double B,G,R,V,VAB,VAG,VAR;
			  V = data_v[i*step+j];
			  B = data_b[i*step+j];
			  VAB = fabs(B-V)/(fabs(1-V/A));   //会有一些值大于256，需要进行归一化
			  if(VAB>255)
				  VAB = 255;
			  else
				  VAB = VAB;
			  data_b[i*step+j] = VAB;
			  G = data_g[i*step+j];
			  VAG = fabs(G-V)/(fabs(1-V/A));
			  if(VAG>255)
				  VAG = 255;
			  else
				  VAG = VAG;
			  data_g[i*step+j] = VAG;
			  R = data_r[i*step+j];
			  VAR = fabs(R-V)/(fabs(1-V/A));
			  if(VAR>255)
				  VAR = 255;
			  else
				  VAR = VAR;
			  data_r[i*step+j] = VAR;
		  };
		};
	cvMerge(b,g,r,NULL,result);//这个函数可能也有问题~
	cout<<"最终去雾算法运行成功"<<"\n";//表示该函数运行成功
	return result;

}

int main(int argc, char** argv)
{
	cvNamedWindow("Source Image");
	cvNamedWindow("Result Image");
	IplImage* image = cvLoadImage("D:/4.bmp",1);  //input a image，0表示以灰度图形式读入图像，-1表示以图像自身真实通道数读入图像，1表示以三通道读入图像
													//此处可改成自己的图片路径
	cvShowImage("Source Image",image);//显示源图像
	int patchsize = 20;
	//IplImage* out = doDarkChannel(image,patchsize);//不能直接将返回的图像数据赋给一个未指定大小的指针，
	IplImage* out = cvCreateImage(cvSize(image->width,image->height),image->depth,3);//创建存储输出图像的内存区域
	IplImage* w = cvCreateImage(cvSize(image->width,image->height),image->depth,1);//创建存储输出图像的内存区域
	//cvCopy(doDarkChannel(image,patchsize),out);//将patchsize大小的图像块存入临时图像块patch_out
	IplImage* smooth = cvCreateImage(cvSize(image->width,image->height),image->depth,1);//创建存储输出图像的内存区域
	IplImage* diff = cvCreateImage(cvSize(image->width,image->height),image->depth,1);//存储w与I_smooth差值绝对值的的内存区域
	IplImage* v = cvCreateImage(cvSize(image->width,image->height),image->depth,1);//存储w与I_smooth差值绝对值的的内存区域
	int A_MAX = doDarkChannel(image,patchsize);//求取暗通道的最大值，A_MAX相当于Matlab中的A，传入doFinally
	cvCopy(doMinChannel(image),w);//计算三个通道的最小值并以最小值组成一副灰度图，进行下一步高斯平滑
	//w计算没问题

	cvSaveImage("D://result//w.bmp",w);

	cvSmooth(w,smooth,CV_GAUSSIAN,39,39,4.5,4.5);//39x39;不使用相关而使用卷积进行计算，将边界点复制得到拓展边界

	cvSaveImage("D://result//smooth.bmp",smooth);

	cvAbsDiff(smooth,w,diff);
	//diff有问题，应该是由于smooth导致的
	cvSaveImage("D://result//diff.bmp",diff);

	cvCopy(doCalculateV(w,diff,smooth),v);//计算v，v的含义从matlab代码中可找到，v传入doFinally进行最终的结果计算
	//v有问题；由于smooth有问题，w没问题，diff有问题，导致v有问题
	cvSaveImage("D://result//v.bmp",v);
	cvCopy(doFinally(image,v,A_MAX),out);//计算最终的去雾结果的函数的调用
	//cvSaveImage("D://v.bmp",v);//测试能否顺利产生图像v的代码
	cout<<"A_MAX="<<A_MAX<<"\n";

	cvSaveImage("D://result//finally.bmp",out);

	cvShowImage("Result Image",out);//imshow the result
	cvWaitKey(0);
	cvReleaseImage(&image);//release the storage space
	cvReleaseImage(&out);//release the storage space
	cvReleaseImage(&w);//release the storage space
	cvReleaseImage(&smooth);//release the storage space
	cvReleaseImage(&diff);//release the storage space
	cvReleaseImage(&v);//release the storage space
	cvDestroyWindow("Source Image");
	cvDestroyWindow("Result Image");
	//system("pause"); //避免一闪而过
	return 0;
}

参考文献：

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

时间： 2024-10-25 14:53:03

OpenCV导向滤波实现（Guided Filter）代码，以及使用颜色先验算法去雾的相关文章

CUDA加opencv复现导向滤波算法

CUDA是GPU通用计算的一种,其中现在大热的深度学习底层GPU计算差不多都选择的CUDA,在这我们先简单了解下其中的一些概念,为了好理解,我们先用DX11里的Compute shader来和CUDA比较下,这二者都可用于GPU通用计算. 先上一张微软MSDN上的图. Compute shader: 线程块: Dispatch(x,y,z), 索引SV_GroupID 线程组: [numthreads(SIZE_X, SIZE_Y, 1)], 索引SV_GroupThreadID. 组内索引:

用 OpenCV 实现 Guided Filter

最近发现 OpenCV 的 C++ 接口写得和 Matlab 实在是太相似了.我用 OpenCV 2.49 来实现 Guided Filter 的时候几乎是直接照抄作者公开的 Matlab 代码.废话不多直接贴我的代码: cv::Mat guidedFilter(cv::Mat I, cv::Mat p, int r, double eps) { /* % GUIDEDFILTER O(1) time implementation of guided filter. % % - guidance

超高速导向滤波实现过程纪要（欢迎挑战）

自从何凯明提出导向滤波后,因为其算法的简单性和有效性,该算法得到了广泛的应用,以至于新版的matlab都将其作为标准自带的函数之一了,利用他可以解决的所有的保边滤波器的能解决的问题,比如细节增强.HDR压缩.细节羽化.去雾.风格化,而且由于其保边特性,如果很多传统函数中使用高斯滤波或者均值滤波的地方用他代替,能很好解决一些强边缘的过渡不自然问题,比如retinex.Highlight/shadow等应用中,因此,快速的实现该算法具有很强的适用意义. 本文简要的记录了本人在优化导向滤波实现的过程中

学习OpenCV——粒子滤波（网上两篇文章总结）

粒子滤波的理论实在是太美妙了,用一组不同权重的随机状态来逼近复杂的概率密度函数.其再非线性.非高斯系统中具有优良的特性.opencv给出了一个实现,但是没有给出范例,学习过程中发现网络上也找不到.learning opencv一书中有介绍,但距离直接使用还是有些距离.在经过一番坎坷后,终于可以用了,希望对你有帮助. 本文中给出的例子跟我的另一篇博文是同一个应用例子,都是对二维坐标进行平滑.预测使用方法: 1.创建并初始化 const int stateNum=4;//状态数 const in

FFmpeg的HEVC解码器源代码简单分析：环路滤波（Loop Filter）

===================================================== HEVC源代码分析文章列表: [解码 -libavcodec HEVC 解码器] FFmpeg的HEVC解码器源代码简单分析:概述 FFmpeg的HEVC解码器源代码简单分析:解析器(Parser)部分 FFmpeg的HEVC解码器源代码简单分析:解码器主干部分 FFmpeg的HEVC解码器源代码简单分析:CTU解码(CTU Decode)部分-PU FFmpeg的HEVC解码器源代码简单

FFmpeg的H.264解码器源代码简单分析：环路滤波（Loop Filter）部分

本文分析FFmpeg的H.264解码器的环路滤波(Loop Filter)部分.FFmpeg的H.264解码器调用decode_slice()函数完成了解码工作.这些解码工作可以大体上分为3个步骤:熵解码,宏块解码以及环路滤波.本文分析这3个步骤中的第3个步骤. 函数调用关系图环路滤波(Loop Filter)部分的源代码在整个H.264解码器中的位置如下图所示. 单击查看更清晰的图片环路滤波(Loop Filter)部分的源代码的调用关系如下图所示. 单击查看更清晰的图片环路滤波主要用于

FilterEngine 类解析——OpenCV图像滤波核心引擎(zz)

<2>FilterEngine 类解析——OpenCV图像滤波核心引擎 FilterEngine类是OpenCV关于图像滤波的主力军类,OpenCV图像滤波功能的核心引擎.各种滤波函数比如blur, GaussianBlur,到头来其实是就是在函数末尾处定义了一个Ptr<FilterEngine>类型的f,然后f->apply( src, dst )了一下而已. 这个类可以把几乎是所有的滤波操作施加到图像上.它包含了所有必要的中间缓存器.有很多和滤波相关的create系函数的

springboot扫描自定义的servlet和filter代码详解_java - JAVA

文章来源:嗨学网敏而好学论坛www.piaodoo.com 欢迎大家相互学习这几天使用spring boot编写公司一个应用,在编写了一个filter,用于指定编码的filter,如下: /** * Created by xiaxuan on 16/11/1. */ @WebFilter(urlPatterns = "/*",filterName="CharacterEncodeFilter", initParams={ @WebInitParam(name=&

图像去雾利用最小值滤波(待改进)

图像去雾我的个神啊~ 调了这么久终于阶段性"胜利了" -- 利用最小值滤波正确实现图像去雾. 因为是C语言实现,很多细节都要注意,很砸时间--本来看起来很简单的程序,就是有bug,debug的时间很长了,这也是一种锻炼吧,难得的实战机会.磨C 这把刀这么久了,是该做做具体应用了!刀子是否锋利,要接受考验! 程序很大程度上没有利用OpenCV的一些数据操作的接口,除了图像显示和读取之外,其他的数据处理部分都是自己实现的.这几天的时间全砸进去了~ 暑假哇~时间过的慢点可好~ 由于项目正在