导向滤波（GuidedFilter）

何凯明去雾算法中的导向滤波实现，原文地址导向滤波

#include "cv.h"
#include "highgui.h"

IplImage * cumsum(IplImage *src,int rc);
IplImage * boxFilter(IplImage *src,int r);
IplImage * myGuidedFilter(IplImage * I,IplImage *img_pp,int r, double eps);

#include "myGuidedFilter.h"

IplImage * cumsum(IplImage *src,int rc)
{
	IplImage *Imdst = cvCloneImage(src);
	if (rc==1)
	{
		for(int y=1;y<src->height;y++)
		{
			for(int x=0;x<src->width;x++)
			{
				cvSetReal2D(Imdst,y,x,cvGetReal2D(Imdst,y-1,x)+cvGetReal2D(Imdst,y,x));
			}
		}
	}
	else if (rc==2)
	{
		for(int y=0;y<src->height;y++)
		{
			for(int x=1;x<src->width;x++)
			{
				cvSetReal2D(Imdst,y,x,cvGetReal2D(Imdst,y,x-1)+cvGetReal2D(Imdst,y,x));
			}
		}
	}
	return Imdst;
}
IplImage * boxFilter(IplImage *src,int r)
{
	IplImage *Imdst = cvCloneImage(src);
	int depth=src->depth;
	IplImage *subImage;

	//imCum = cumsum(imSrc, 1);
	IplImage *imCum = cumsum(Imdst,1);

	//imDst(1:r+1, :) = imCum(1+r:2*r+1, :);
	for (int y = 0;y<r;y++)
	{
		for(int x = 0;x<Imdst->width;x++)
		{
			cvSetReal2D(Imdst,y,x,cvGetReal2D(imCum,y+r,x));
		}
	}

	//imDst(r+2:hei-r, :) = imCum(2*r+2:hei, :) - imCum(1:hei-2*r-1, :);
	for (int y = r+1;y<Imdst->height-r-1;y++)
	{
		for(int x = 0;x<Imdst->width;x++)
		{
			cvSetReal2D(Imdst,y,x,(cvGetReal2D(imCum,y+r,x)-cvGetReal2D(imCum,y-r-1,x)));
		}
	}

	//imDst(hei-r+1:hei, :) = repmat(imCum(hei, :), [r, 1]) - imCum(hei-2*r:hei-r-1, :);
	subImage = cvCreateImage(cvSize(Imdst->width,r),depth,1);
	for(int y=0;y<r;y++)
	{
		for(int x=0;x<Imdst->width;x++)
		{
			cvSetReal2D(subImage,y,x,cvGetReal2D(imCum,Imdst->height-1,x));
		}
	}
	for (int y = Imdst->height-r;y<Imdst->height;y++)
	{
		for(int x = 0;x<Imdst->width;x++)
		{
			cvSetReal2D(Imdst,y,x,cvGetReal2D(subImage,y-Imdst->height+r,x)-cvGetReal2D(imCum,y-r-1,x));
		}
	}
	cvReleaseImage(&subImage);
 
	imCum = cumsum(Imdst, 2);
	//imDst(:, 1:r+1) = imCum(:, 1+r:2*r+1);
	for (int y = 0;y<Imdst->height;y++)
	{
		for(int x = 0;x<r;x++)
		{
			cvSetReal2D(Imdst,y,x,cvGetReal2D(imCum,y,x+r));
		}
	}

	//imDst(:, r+2:wid-r) = imCum(:, 2*r+2:wid) - imCum(:, 1:wid-2*r-1);
	for (int y = 0;y<Imdst->height;y++)
	{
		for(int x = r+1;x<Imdst->width-r-1;x++)
		{
			cvSetReal2D(Imdst,y,x,(cvGetReal2D(imCum,y,x+r)-cvGetReal2D(imCum,y,x-r-1)));
		}
	}
	cvReleaseImage(&subImage);
 
	//imDst(:, wid-r+1:wid) = repmat(imCum(:, wid), [1, r]) - imCum(:, wid-2*r:wid-r-1);
	subImage = cvCreateImage(cvSize(r,Imdst->height),depth,1);
	for(int y=0;y<Imdst->height;y++)
	{
		for(int x=0;x<r;x++)
		{
			cvSetReal2D(subImage,y,x,cvGetReal2D(imCum,y,Imdst->width-1));
		}
	}
	for (int y = 0;y<Imdst->height;y++)
	{
		for(int x = Imdst->width-r;x<Imdst->width;x++)
		{
			cvSetReal2D(Imdst,y,x,cvGetReal2D(subImage,y,x-Imdst->width+r)-cvGetReal2D(imCum,y,x-r-1));
		}
	}
	cvReleaseImage(&subImage);

	return Imdst;
}

IplImage * myGuidedFilter(IplImage * I,IplImage *img_pp,int r, double eps)
{
	int height = img_pp->height;
	int width = img_pp->width;
	int depth = IPL_DEPTH_64F;

	IplImage *ones = cvCreateImage(cvSize(width,height),depth,1);
	cvSet(ones,cvRealScalar(1));
	IplImage * N = boxFilter(ones,r);

	//求I的均值
	IplImage * mean_I = cvCreateImage(cvSize(width,height),depth,1);
	cvDiv(boxFilter(I,r),N,mean_I);
	//求P的均值
	IplImage * mean_p = cvCreateImage(cvSize(width,height),depth,1);
	cvDiv(boxFilter(img_pp,r),N,mean_p);
	//求I*P的均值
	IplImage * pr = cvCreateImage(cvSize(width,height),depth,1);
	cvMul(I,img_pp,pr);
	IplImage * mean_Ip = cvCreateImage(cvSize(width,height),depth,1);
	cvDiv(boxFilter(pr,r),N,mean_Ip); 

	//求I与p协方差
	cvMul(mean_I,mean_p,pr);
	IplImage * cov_Ip = cvCreateImage(cvSize(width,height),depth,1);
	cvSub(mean_Ip,pr,cov_Ip);

	//求I的方差
	IplImage * var_I  = cvCreateImage(cvSize(width,height),depth,1);
	cvMul(I,I,pr);
	cvDiv(boxFilter(pr,r),N,var_I);
	cvMul(mean_I,mean_I,pr);
	cvSub(var_I,pr,var_I);
	
	//求a
	IplImage * a = cvCreateImage(cvSize(width,height),depth,1);
	cvAddS(var_I,cvScalar(eps),var_I);
	cvDiv(cov_Ip,var_I,a);
	//求b
	IplImage * b = cvCreateImage(cvSize(width,height),depth,1);
	cvMul(a,mean_I,pr);
	cvSub(mean_p,pr,b);
		 
	//求a的均值
	IplImage * mean_a = cvCreateImage(cvSize(width,height),depth,1);
	cvDiv(boxFilter(a,r),N,mean_a);
	//求b的均值
	IplImage * mean_b = cvCreateImage(cvSize(width,height),depth,1);
	cvDiv(boxFilter(b,r),N,mean_b);
		 
	//求Q
	IplImage * q = cvCreateImage(cvSize(width,height),depth,1);
	cvMul(mean_a,I,a);
	cvAdd(a,mean_b,q);
		 
	cvReleaseImage(&ones);
	cvReleaseImage(&mean_I);
	cvReleaseImage(&mean_p);
	cvReleaseImage(&pr);
	cvReleaseImage(&mean_Ip);
	cvReleaseImage(&cov_Ip);
	cvReleaseImage(&var_I);
	cvReleaseImage(&a);
	cvReleaseImage(&b);
	cvReleaseImage(&mean_a);
	cvReleaseImage(&mean_b);
	return q;
}

导向滤波（GuidedFilter）,布布扣,bubuko.com

导向滤波（GuidedFilter）

原文：http://blog.csdn.net/wds555/article/details/23176313