导向滤波

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

导向图像I，滤波输入图像p以及输出图像q。像素点 i 处的滤波结果是被表达成一个加权平均：


假设导向滤波器在导向图像I和滤波输出q之间是一个局部线性模型：

最小化下面的窗口Wk的代价函数：


用来确定a,b的值

其中

论文所给算法如下：

matlab代码如下：

[plain]  view plain copy

1. function q = guidedfilter(I, p, r, eps)  
2. %   GUIDEDFILTER   O(1) time implementation of guided filter.  
3. %  
4. %   - guidance image: I (should be a gray-scale/single channel image)  
5. %   - filtering input image: p (should be a gray-scale/single channel image)  
6. %   - local window radius: r  
7. %   - regularization parameter: eps  
8.   
9. [hei, wid] = size(I);  
10. N = boxfilter(ones(hei, wid), r); % the size of each local patch; N=(2r+1)^2 except for boundary pixels.  
11.   
12. mean_I = boxfilter(I, r) ./ N;  
13. mean_p = boxfilter(p, r) ./ N;  
14. mean_Ip = boxfilter(I.*p, r) ./ N;  
15. cov_Ip = mean_Ip - mean_I .* mean_p; % this is the covariance of (I, p) in each local patch.  
16.   
17. mean_II = boxfilter(I.*I, r) ./ N;  
18. var_I = mean_II - mean_I .* mean_I;  
19.   
20. a = cov_Ip ./ (var_I + eps); % Eqn. (5) in the paper;  
21. b = mean_p - a .* mean_I; % Eqn. (6) in the paper;  
22.   
23. mean_a = boxfilter(a, r) ./ N;  
24. mean_b = boxfilter(b, r) ./ N;  
25.   
26. q = mean_a .* I + mean_b; % Eqn. (8) in the paper;  
27. end  

[plain]  view plain copy

1. function imDst = boxfilter(imSrc, r)  
2.   
3. %   BOXFILTER   O(1) time box filtering using cumulative sum  
4. %  
5. %   - Definition imDst(x, y)=sum(sum(imSrc(x-r:x+r,y-r:y+r)));  
6. %   - Running time independent of r;   
7. %   - Equivalent to the function: colfilt(imSrc, [2*r+1, 2*r+1], 'sliding', @sum);  
8. %   - But much faster.  
9.   
10. [hei, wid] = size(imSrc);  
11. imDst = zeros(size(imSrc));  
12.   
13. %cumulative sum over Y axis  
14. imCum = cumsum(imSrc, 1);  
15. %difference over Y axis  
16. imDst(1:r+1, :) = imCum(1+r:2*r+1, :);  
17. imDst(r+2:hei-r, :) = imCum(2*r+2:hei, :) - imCum(1:hei-2*r-1, :);  
18. imDst(hei-r+1:hei, :) = repmat(imCum(hei, :), [r, 1]) - imCum(hei-2*r:hei-r-1, :);  
19.   
20. %cumulative sum over X axis  
21. imCum = cumsum(imDst, 2);  
22. %difference over Y axis  
23. imDst(:, 1:r+1) = imCum(:, 1+r:2*r+1);  
24. imDst(:, r+2:wid-r) = imCum(:, 2*r+2:wid) - imCum(:, 1:wid-2*r-1);  
25. imDst(:, wid-r+1:wid) = repmat(imCum(:, wid), [1, r]) - imCum(:, wid-2*r:wid-r-1);  
26. end  

以下为单通道图像导向滤波opencv实现：

[cpp]  view plain copy

1. #include "myGuidedFilter_Mat.h"  
2.   
3. CvMat * cumsum(CvMat *src,int rc)  
4. {  
5.     CvMat *Imdst = cvCreateMat(src->rows,src->cols,CV_64FC1);  
6.     Imdst=cvCloneMat(src);  
7.     if (rc==1)  
8.     {  
9.         for(int y=1;y<src->height;y++)  
10.         {  
11.             double *ptr0=(double *)(Imdst->data.ptr+(y-1)*Imdst->step);  
12.             double *ptr=(double *)(Imdst->data.ptr+y*Imdst->step);  
13.             for(int x=0;x<src->width;x++)  
14.             {  
15.                 ptr[x]=ptr0[x]+ptr[x];  
16.                 //cvSetReal2D(Imdst,y,x,cvGetReal2D(Imdst,y-1,x)+cvGetReal2D(Imdst,y,x));  
17.             }  
18.         }  
19.     }  
20.     else if (rc==2)  
21.     {  
22.         for(int y=0;y<src->height;y++)  
23.         {  
24.             double *ptr=(double *)(Imdst->data.ptr+y*Imdst->step);  
25.             for(int x=1;x<src->width;x++)  
26.             {  
27.                 ptr[x]=ptr[x-1]+ptr[x];  
28.                 //cvSetReal2D(Imdst,y,x,cvGetReal2D(Imdst,y,x-1)+cvGetReal2D(Imdst,y,x));  
29.             }  
30.         }  
31.     }  
32.     return Imdst;  
33. }  
34. CvMat * boxFilter(CvMat *src,int r)  
35. {  
36.     CvMat *Imdst = cvCreateMat(src->rows,src->cols,CV_64FC1);  
37.     Imdst=cvCloneMat(src);  
38.     CvMat *subImage;  
39.   
40.     //imCum = cumsum(imSrc, 1);  
41.     CvMat *imCum = cumsum(Imdst,1);  
42.   
43.     //imDst(1:r+1, :) = imCum(1+r:2*r+1, :);  
44.     for (int y = 0;y<r;y++)  
45.     {  
46.         //double *ptrDst=(double *)Imdst->data.ptr+y*Imdst->step;  
47.         //double *ptrCum=(double *)imCum->data.ptr+(y+r)*imCum->step;  
48.         for(int x = 0;x<Imdst->width;x++)  
49.         {  
50.             //ptrDst[x]=ptrCum[x];  
51.             cvSetReal2D(Imdst,y,x,cvGetReal2D(imCum,y+r,x));  
52.         }  
53.     }  
54.   
55.     //imDst(r+2:hei-r, :) = imCum(2*r+2:hei, :) - imCum(1:hei-2*r-1, :);  
56.     for (int y = r+1;y<Imdst->height-r-1;y++)  
57.     {  
58.         for(int x = 0;x<Imdst->width;x++)  
59.         {  
60.             cvSetReal2D(Imdst,y,x,(cvGetReal2D(imCum,y+r,x)-cvGetReal2D(imCum,y-r-1,x)));  
61.         }  
62.     }  
63.   
64.     //imDst(hei-r+1:hei, :) = repmat(imCum(hei, :), [r, 1]) - imCum(hei-2*r:hei-r-1, :);  
65.     subImage = cvCreateMat(r,Imdst->width,CV_64FC1);  
66.     CvMat *tem=cvCreateMat(1,Imdst->width,CV_64FC1);  
67.     cvGetRow(imCum,tem,imCum->height-1);  
68.     cvRepeat(tem,subImage);  
69.     /*for(int y=0;y<r;y++) 
70.     { 
71.         for(int x=0;x<Imdst->width;x++) 
72.         { 
73.             cvSetReal2D(subImage,y,x,cvGetReal2D(imCum,Imdst->height-1,x)); 
74.         } 
75.     }*/  
76.     for (int y = Imdst->height-r;y<Imdst->height;y++)  
77.     {  
78.         for(int x = 0;x<Imdst->width;x++)  
79.         {  
80.             cvSetReal2D(Imdst,y,x,cvGetReal2D(subImage,y-Imdst->height+r,x)-cvGetReal2D(imCum,y-r-1,x));  
81.         }  
82.     }  
83.     cvReleaseMat(&subImage);  
84.     cvReleaseMat(&tem);  
85.    
86.     imCum = cumsum(Imdst, 2);  
87.     //imDst(:, 1:r+1) = imCum(:, 1+r:2*r+1);  
88.     for (int y = 0;y<Imdst->height;y++)  
89.     {  
90.         for(int x = 0;x<r;x++)  
91.         {  
92.             cvSetReal2D(Imdst,y,x,cvGetReal2D(imCum,y,x+r));  
93.         }  
94.     }  
95.   
96.     //imDst(:, r+2:wid-r) = imCum(:, 2*r+2:wid) - imCum(:, 1:wid-2*r-1);  
97.     for (int y = 0;y<Imdst->height;y++)  
98.     {  
99.         for(int x = r+1;x<Imdst->width-r-1;x++)  
100.         {  
101.             cvSetReal2D(Imdst,y,x,(cvGetReal2D(imCum,y,x+r)-cvGetReal2D(imCum,y,x-r-1)));  
102.         }  
103.     }  
104.    
105.     //imDst(:, wid-r+1:wid) = repmat(imCum(:, wid), [1, r]) - imCum(:, wid-2*r:wid-r-1);  
106.     subImage = cvCreateMat(Imdst->height,r,CV_64FC1);  
107.     tem=cvCreateMat(Imdst->height,1,CV_64FC1);  
108.     cvGetCol(imCum,tem,imCum->width-1);  
109.     cvRepeat(tem,subImage);  
110.     /*for(int y=0;y<Imdst->height;y++) 
111.     { 
112.         for(int x=0;x<r;x++) 
113.         { 
114.             cvSetReal2D(subImage,y,x,cvGetReal2D(imCum,y,Imdst->width-1)); 
115.         } 
116.     }*/  
117.     for (int y = 0;y<Imdst->height;y++)  
118.     {  
119.         for(int x = Imdst->width-r;x<Imdst->width;x++)  
120.         {  
121.             cvSetReal2D(Imdst,y,x,cvGetReal2D(subImage,y,x-Imdst->width+r)-cvGetReal2D(imCum,y,x-r-1));  
122.         }  
123.     }  
124.     cvReleaseMat(&subImage);  
125.   
126.     return Imdst;  
127. }  
128.   
129. CvMat * myGuidedFilter_Mat(CvMat * I,CvMat *img_pp,int r, double eps)  
130. {  
131.     int height = img_pp->height;  
132.     int width = img_pp->width;  
133.     int type = CV_64FC1;  
134.   
135.     CvMat *ones = cvCreateMat(height,width,type);  
136.     cvSet(ones,cvRealScalar(1));  
137.     CvMat * N = boxFilter(ones,r);  
138.   
139.     //求I的均值  
140.     CvMat * mean_I = cvCreateMat(height,width,type);  
141.     cvDiv(boxFilter(I,r),N,mean_I);  
142.     //求P的均值  
143.     CvMat * mean_p = cvCreateMat(height,width,type);  
144.     cvDiv(boxFilter(img_pp,r),N,mean_p);  
145.     //求I*P的均值  
146.     CvMat * pr = cvCreateMat(height,width,type);  
147.     cvMul(I,img_pp,pr);  
148.     CvMat * mean_Ip = cvCreateMat(height,width,type);  
149.     cvDiv(boxFilter(pr,r),N,mean_Ip);   
150.   
151.     //求I与p协方差  
152.     cvMul(mean_I,mean_p,pr);  
153.     CvMat * cov_Ip = cvCreateMat(height,width,type);  
154.     cvSub(mean_Ip,pr,cov_Ip);  
155.   
156.     //求I的方差  
157.     CvMat * var_I  = cvCreateMat(height,width,type);  
158.     cvMul(I,I,pr);  
159.     cvDiv(boxFilter(pr,r),N,var_I);  
160.     cvMul(mean_I,mean_I,pr);  
161.     cvSub(var_I,pr,var_I);  
162.       
163.     //求a  
164.     CvMat * a = cvCreateMat(height,width,type);  
165.     cvAddS(var_I,cvScalar(eps),var_I);  
166.     cvDiv(cov_Ip,var_I,a);  
167.     //求b  
168.     CvMat * b = cvCreateMat(height,width,type);  
169.     cvMul(a,mean_I,pr);  
170.     cvSub(mean_p,pr,b);  
171.            
172.     //求a的均值  
173.     CvMat * mean_a = cvCreateMat(height,width,type);  
174.     cvDiv(boxFilter(a,r),N,mean_a);  
175.     //求b的均值  
176.     CvMat * mean_b = cvCreateMat(height,width,type);  
177.     cvDiv(boxFilter(b,r),N,mean_b);  
178.            
179.     //求Q  
180.     CvMat * q = cvCreateMat(height,width,type);  
181.     cvMul(mean_a,I,a);  
182.     cvAdd(a,mean_b,q);  
183.            
184.     cvReleaseMat(&ones);  
185.     cvReleaseMat(&mean_I);  
186.     cvReleaseMat(&mean_p);  
187.     cvReleaseMat(&pr);  
188.     cvReleaseMat(&mean_Ip);  
189.     cvReleaseMat(&cov_Ip);  
190.     cvReleaseMat(&var_I);  
191.     cvReleaseMat(&a);  
192.     cvReleaseMat(&b);  
193.     cvReleaseMat(&mean_a);  
194.     cvReleaseMat(&mean_b);  
195.     return q;  
196. }