opencv實現導向濾波（GuidedFilter）
阿新 • • 發佈：2019-01-04
何凱明去霧演算法中的導向濾波實現，原文地址導向濾波。
導向影象I，濾波輸入影象p以及輸出影象q。畫素點 i 處的濾波結果是被表達成一個加權平均：
假設導向濾波器在導向影象I和濾波輸出q之間是一個區域性線性模型：
最小化下面的視窗Wk的代價函式：
用來確定a,b的值
其中
論文所給演算法如下：
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
以下為單通道影象導向濾波opencv實現：
#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; }