導向濾波小結：從導向濾波（guided filter）到快速導向濾波（fast guide filter）的原理，應用及opencv實現程式碼
阿新 • • 發佈：2019-02-07
1. 導向濾波簡介
導向濾波是何凱明在學生時代提出的一個保邊濾波（edge-preserving smoothing）演算法。何凱明在cv圈應該算是名人了，學生時代關於影象去霧的研究就以第一作者的身份獲得Best Paper Award（CVPR 2009），而且今年剛剛又斬獲Marr Prize（ICCV 2017）。更多關於何凱明的最新研究動態可以點選以下連結何凱明。
導向濾波顧名思義，就是有選擇（導向）性的濾波，其與我們經常提及的高斯濾波、雙邊濾波相比，它具有導向性，說具體點就是，它通過輸入一副影象（矩陣）作為導向圖，這樣濾波器就知道什麼地方是邊緣，這樣就可以更好的保護邊緣，最終達到在濾波的同時，保持邊緣細節。所以有個說法是導向濾波是各向異性的濾波器，而高斯濾波、雙邊濾波這些是各向同性濾波器，我覺得也是很貼切。
導向濾波作為一種保邊濾波，可以運用在很多場合，比如美顏，去霧，三維重建等。
如果你僅僅只是需要運用這個演算法，現在opencv 3.0和MATLAB 14都已經添加了guided filter的API，可以直接呼叫。
opencv中的API如下void cv::ximgproc::guidedFilter()，具體的可以參考opencv的幫助文件關於導向濾波的介紹guidedFilter。
但是需要提醒的是，opencv中guidedFilter()函式包含在ximgproc模組下，但是從官方下載的標準的opencv.exe程式中並不包含該模組，需要分別下載opencv的source檔案和contrib模組的source檔案，然後自己編譯，具體可以參考 opencv3.1.0+contrib模組編譯總結。
2. 導向濾波的原理
查看了很多相關的資料，覺得白馬負金羈的導向濾波(Guided Filter)的解析與實現一文將其原理解釋的非常通俗易懂了，這裡就不再贅述。僅給出最後的推導結果，其中fmean為一個視窗半徑為r的均值濾波器（對應的視窗大小為2*r+1），corr為相關，var為方差，cov為協方差。
cv::Mat GuidedFilter(cv::Mat I, cv::Mat p, int r, double eps)
{
/*
% GUIDEDFILTER O(N) 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
*/
cv::Mat _I;
I.convertTo(_I, CV_64FC1,1.0/255);
I = _I;
cv::Mat _p;
p.convertTo(_p, CV_64FC1,1.0/255);
p = _p;
//[hei, wid] = size(I);
int hei = I.rows;
int wid = I.cols;
r=2*r+1;//因為opencv自帶的boxFilter（）中的Size,比如9x9,我們說半徑為4
//mean_I = boxfilter(I, r) ./ N;
cv::Mat mean_I;
cv::boxFilter(I, mean_I, CV_64FC1, cv::Size(r, r));
//mean_p = boxfilter(p, r) ./ N;
cv::Mat mean_p;
cv::boxFilter(p, mean_p, CV_64FC1, cv::Size(r, r));
//mean_Ip = boxfilter(I.*p, r) ./ N;
cv::Mat mean_Ip;
cv::boxFilter(I.mul(p), mean_Ip, CV_64FC1, cv::Size(r, r));
//cov_Ip = mean_Ip - mean_I .* mean_p; % this is the covariance of (I, p) in each local patch.
cv::Mat cov_Ip = mean_Ip - mean_I.mul(mean_p);
//mean_II = boxfilter(I.*I, r) ./ N;
cv::Mat mean_II;
cv::boxFilter(I.mul(I), mean_II, CV_64FC1, cv::Size(r, r));
//var_I = mean_II - mean_I .* mean_I;
cv::Mat var_I = mean_II - mean_I.mul(mean_I);
//a = cov_Ip ./ (var_I + eps); % Eqn. (5) in the paper;
cv::Mat a = cov_Ip / (var_I + eps);
//b = mean_p - a .* mean_I; % Eqn. (6) in the paper;
cv::Mat b = mean_p - a.mul(mean_I);
//mean_a = boxfilter(a, r) ./ N;
cv::Mat mean_a;
cv::boxFilter(a, mean_a, CV_64FC1, cv::Size(r, r));
//mean_b = boxfilter(b, r) ./ N;
cv::Mat mean_b;
cv::boxFilter(b, mean_b, CV_64FC1, cv::Size(r, r));
//q = mean_a .* I + mean_b; % Eqn. (8) in the paper;
cv::Mat q = mean_a.mul(I) + mean_b;
return q;
}
cv::Mat fastGuidedFilter(cv::Mat I_org, cv::Mat p_org, int r, double eps, int s)
{
/*
% GUIDEDFILTER O(N) 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
*/
cv::Mat I,_I;
I_org.convertTo(_I, CV_64FC1, 1.0 / 255);
resize(_I,I,Size(),1.0/s,1.0/s,1);
cv::Mat p,_p;
p_org.convertTo(_p, CV_64FC1, 1.0 / 255);
//p = _p;
resize(_p, p, Size(),1.0/s,1.0/s,1);
//[hei, wid] = size(I);
int hei = I.rows;
int wid = I.cols;
r = (2 * r + 1)/s+1;//因為opencv自帶的boxFilter（）中的Size,比如9x9,我們說半徑為4
//mean_I = boxfilter(I, r) ./ N;
cv::Mat mean_I;
cv::boxFilter(I, mean_I, CV_64FC1, cv::Size(r, r));
//mean_p = boxfilter(p, r) ./ N;
cv::Mat mean_p;
cv::boxFilter(p, mean_p, CV_64FC1, cv::Size(r, r));
//mean_Ip = boxfilter(I.*p, r) ./ N;
cv::Mat mean_Ip;
cv::boxFilter(I.mul(p), mean_Ip, CV_64FC1, cv::Size(r, r));
//cov_Ip = mean_Ip - mean_I .* mean_p; % this is the covariance of (I, p) in each local patch.
cv::Mat cov_Ip = mean_Ip - mean_I.mul(mean_p);
//mean_II = boxfilter(I.*I, r) ./ N;
cv::Mat mean_II;
cv::boxFilter(I.mul(I), mean_II, CV_64FC1, cv::Size(r, r));
//var_I = mean_II - mean_I .* mean_I;
cv::Mat var_I = mean_II - mean_I.mul(mean_I);
//a = cov_Ip ./ (var_I + eps); % Eqn. (5) in the paper;
cv::Mat a = cov_Ip / (var_I + eps);
//b = mean_p - a .* mean_I; % Eqn. (6) in the paper;
cv::Mat b = mean_p - a.mul(mean_I);
//mean_a = boxfilter(a, r) ./ N;
cv::Mat mean_a;
cv::boxFilter(a, mean_a, CV_64FC1, cv::Size(r, r));
Mat rmean_a;
resize(mean_a, rmean_a, Size(I_org.cols, I_org.rows),1);
//mean_b = boxfilter(b, r) ./ N;
cv::Mat mean_b;
cv::boxFilter(b, mean_b, CV_64FC1, cv::Size(r, r));
Mat rmean_b;
resize(mean_b, rmean_b, Size(I_org.cols, I_org.rows),1);
//q = mean_a .* I + mean_b; % Eqn. (8) in the paper;
cv::Mat q = rmean_a.mul(_I) + rmean_b;
return q;
}取s=8，計算結果和之前的guidedFilter()的計算結果如下所示：
#include<opencv.hpp>
#include<ximgproc.hpp>
int main(void)
{
cv::Mat src = cv::imread("d:/Opencv Picture/bilateral filter.png", 1);
cv::imshow("src", src);
cv::Mat dst(src.size(), src.type());
float eps = 0.02 * 255 * 255;//eps的取值很關鍵（乘於255的平方）
cv::ximgproc::guidedFilter(src,src,dst,16,eps,-1);
cv::imshow("dst", dst);
cvWaitKey();
return 0;
}