# opencv实现导向滤波（GuidedFilter）

matlab代码如下：

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.

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
function imDst = boxfilter(imSrc, r)

%   BOXFILTER   O(1) time box filtering using cumulative sum
%
%   - Definition imDst(x, y)=sum(sum(imSrc(x-r:x+r,y-r:y+r)));
%   - Running time independent of r;
%   - Equivalent to the function: colfilt(imSrc, [2*r+1, 2*r+1], 'sliding', @sum);
%   - But much faster.

[hei, wid] = size(imSrc);
imDst = zeros(size(imSrc));

%cumulative sum over Y axis
imCum = cumsum(imSrc, 1);
%difference over Y axis
imDst(1:r+1, :) = imCum(1+r:2*r+1, :);
imDst(r+2:hei-r, :) = imCum(2*r+2:hei, :) - imCum(1:hei-2*r-1, :);
imDst(hei-r+1:hei, :) = repmat(imCum(hei, :), [r, 1]) - imCum(hei-2*r:hei-r-1, :);

%cumulative sum over X axis
imCum = cumsum(imDst, 2);
%difference over Y axis
imDst(:, 1:r+1) = imCum(:, 1+r:2*r+1);
imDst(:, r+2:wid-r) = imCum(:, 2*r+2:wid) - imCum(:, 1:wid-2*r-1);
imDst(:, wid-r+1:wid) = repmat(imCum(:, wid), [1, r]) - imCum(:, wid-2*r:wid-r-1);
end


#include "myGuidedFilter_Mat.h"

CvMat * cumsum(CvMat *src,int rc)
{
CvMat *Imdst = cvCreateMat(src->rows,src->cols,CV_64FC1);
Imdst=cvCloneMat(src);
if (rc==1)
{
for(int y=1;y<src->height;y++)
{
double *ptr0=(double *)(Imdst->data.ptr+(y-1)*Imdst->step);
double *ptr=(double *)(Imdst->data.ptr+y*Imdst->step);
for(int x=0;x<src->width;x++)
{
ptr[x]=ptr0[x]+ptr[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++)
{
double *ptr=(double *)(Imdst->data.ptr+y*Imdst->step);
for(int x=1;x<src->width;x++)
{
ptr[x]=ptr[x-1]+ptr[x];
//cvSetReal2D(Imdst,y,x,cvGetReal2D(Imdst,y,x-1)+cvGetReal2D(Imdst,y,x));
}
}
}
return Imdst;
}
CvMat * boxFilter(CvMat *src,int r)
{
CvMat *Imdst = cvCreateMat(src->rows,src->cols,CV_64FC1);
Imdst=cvCloneMat(src);
CvMat *subImage;

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

//imDst(1:r+1, :) = imCum(1+r:2*r+1, :);
for (int y = 0;y<r;y++)
{
//double *ptrDst=(double *)Imdst->data.ptr+y*Imdst->step;
//double *ptrCum=(double *)imCum->data.ptr+(y+r)*imCum->step;
for(int x = 0;x<Imdst->width;x++)
{
//ptrDst[x]=ptrCum[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 = cvCreateMat(r,Imdst->width,CV_64FC1);
CvMat *tem=cvCreateMat(1,Imdst->width,CV_64FC1);
cvGetRow(imCum,tem,imCum->height-1);
cvRepeat(tem,subImage);
/*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));
}
}
cvReleaseMat(&subImage);
cvReleaseMat(&tem);

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)));
}
}

//imDst(:, wid-r+1:wid) = repmat(imCum(:, wid), [1, r]) - imCum(:, wid-2*r:wid-r-1);
subImage = cvCreateMat(Imdst->height,r,CV_64FC1);
tem=cvCreateMat(Imdst->height,1,CV_64FC1);
cvGetCol(imCum,tem,imCum->width-1);
cvRepeat(tem,subImage);
/*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));
}
}
cvReleaseMat(&subImage);

return Imdst;
}

CvMat * myGuidedFilter_Mat(CvMat * I,CvMat *img_pp,int r, double eps)
{
int height = img_pp->height;
int width = img_pp->width;
int type = CV_64FC1;

CvMat *ones = cvCreateMat(height,width,type);
cvSet(ones,cvRealScalar(1));
CvMat * N = boxFilter(ones,r);

//求I的均值
CvMat * mean_I = cvCreateMat(height,width,type);
cvDiv(boxFilter(I,r),N,mean_I);
//求P的均值
CvMat * mean_p = cvCreateMat(height,width,type);
cvDiv(boxFilter(img_pp,r),N,mean_p);
//求I*P的均值
CvMat * pr = cvCreateMat(height,width,type);
cvMul(I,img_pp,pr);
CvMat * mean_Ip = cvCreateMat(height,width,type);
cvDiv(boxFilter(pr,r),N,mean_Ip);

//求I与p协方差
cvMul(mean_I,mean_p,pr);
CvMat * cov_Ip = cvCreateMat(height,width,type);
cvSub(mean_Ip,pr,cov_Ip);

//求I的方差
CvMat * var_I  = cvCreateMat(height,width,type);
cvMul(I,I,pr);
cvDiv(boxFilter(pr,r),N,var_I);
cvMul(mean_I,mean_I,pr);
cvSub(var_I,pr,var_I);

//求a
CvMat * a = cvCreateMat(height,width,type);
cvDiv(cov_Ip,var_I,a);
//求b
CvMat * b = cvCreateMat(height,width,type);
cvMul(a,mean_I,pr);
cvSub(mean_p,pr,b);

//求a的均值
CvMat * mean_a = cvCreateMat(height,width,type);
cvDiv(boxFilter(a,r),N,mean_a);
//求b的均值
CvMat * mean_b = cvCreateMat(height,width,type);
cvDiv(boxFilter(b,r),N,mean_b);

//求Q
CvMat * q = cvCreateMat(height,width,type);
cvMul(mean_a,I,a);

cvReleaseMat(&ones);
cvReleaseMat(&mean_I);
cvReleaseMat(&mean_p);
cvReleaseMat(&pr);
cvReleaseMat(&mean_Ip);
cvReleaseMat(&cov_Ip);
cvReleaseMat(&var_I);
cvReleaseMat(&a);
cvReleaseMat(&b);
cvReleaseMat(&mean_a);
cvReleaseMat(&mean_b);
return q;
}

05-05 2万+
10-28 3万+
04-20 732
10-17 2558
07-10 6万+
09-18 1493
12-12 959
11-18 996
04-21 3700
06-06 1万+
10-18
04-17 1136