导向滤波

最新推荐文章于 2024-08-10 20:07:01 发布

todayq

最新推荐文章于 2024-08-10 20:07:01 发布

阅读量3.8k

点赞数 1

分类专栏：数字图像处理

数字图像处理专栏收录该内容

72 篇文章 0 订阅

订阅专栏

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

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

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

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

用来确定a,b的值

其中

论文所给算法如下：

matlab代码如下：

[plain]view plaincopy 
    
 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  

[plain]view plaincopy 
    
 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  

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

[cpp]view plaincopy 
    
 #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);  
     cvAddS(var_I,cvScalar(eps),var_I);  
     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);  
     cvAdd(a,mean_b,q);  
            
     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;  
 }  