綜述：圖像濾波常用算法實現及原理解析

綜述：圖像濾波常用算法實現及原理解析

CV開發者都愛看的 計算機視覺工坊 昨天

點擊上方“計算機視覺工坊”，選擇“星標”

干貨第一時間送達

作者丨一支程序媛@知乎（已授權）來源丨https://zhuanlan.zhihu.com/p/279602383編輯丨極市平台

導讀

 

圖像濾波是一種非常重要的圖像處理技術，本文詳細介紹了四種常見的圖像濾波算法，並附上源碼，包括自適應中值濾波、高斯濾波、雙邊濾波和導向濾波。 

前言

本文介紹四種常見的圖像濾波算法，並附上源碼。圖像濾波是一種非常重要的圖像處理技術，現在大火的卷積神經網絡其實也是濾波的一種，都是用卷積核去提取圖像的特征模式。不過，傳統的濾波，使用的卷積核是固定的參數，是由經驗非常豐富的人去手動設計的，也稱為手工特征。而卷積神經網絡的卷積核參數初始時未知的，根據不同的任務由數據和神經網絡反向傳播算法去學習得到的參數，更能適應於不同的任務。

目錄

• 自適應中值濾波
• 高斯濾波
• 雙邊濾波
• 導向濾波

自適應中值濾波

中值濾波器

中值濾波器是一種常用的非線性濾波器，其基本原理是：選擇待處理像素的一個鄰域中各像素值的中值來代替待處理的像素。主要功能使某像素的灰度值與周圍領域內的像素比較接近，從而消除一些孤立的噪聲點，所以中值濾波器能夠很好的消除椒鹽噪聲。不僅如此，中值濾波器在消除噪聲的同時，還能有效的保護圖像的邊界信息，不會對圖像造成很大的模糊（相比於均值濾波器）。

中值濾波器的效果受濾波窗口尺寸的影響較大，在消除噪聲和保護圖像的細節存在着矛盾：濾波窗口較小，則能很好的保護圖像中的某些細節，但對噪聲的過濾效果就不是很好，因為實際中的噪聲不可能只占一個像素位置；反之，窗口尺寸較大有較好的噪聲過濾效果，但是會對圖像造成一定的模糊。另外，根據中值濾波器原理，如果在濾波窗口內的噪聲點的個數大於整個窗口內非噪聲像素的個數，則中值濾波就不能很好的過濾掉噪聲。

自適應中值濾波器

常規的中值濾波器，在噪聲的密度不是很大的情況下，效果不錯。但是當噪聲出現的概率較高時，常規的中值濾波的效果就不是很好了。有一個選擇就是增大濾波器的窗口大小，這雖然在一定程度上能解決上述的問題，但是會給圖像造成較大的模糊。

常規的中值濾波器的窗口尺寸是固定大小不變的，就不能同時兼顧去噪和保護圖像的細節。這時就要尋求一種改變，根據預先設定好的條件，在濾波的過程中，動態的改變濾波器的窗口尺寸大小，這就是自適應中值濾波器 Adaptive Median Filter。在濾波的過程中，自適應中值濾波器會根據預先設定好的條件，改變濾波窗口的尺寸大小，同時還會根據一定的條件判斷當前像素是不是噪聲，如果是則用鄰域中值替換掉當前像素；不是，則不作改變。

自適應中值濾波器有三個目的：

• 濾除椒鹽噪聲
• 平滑其他非脈沖噪聲
• 盡可能的保護圖像中細節信息，避免圖像邊緣的細化或者粗化。

自適應中值濾波算法描述

自適應濾波器不但能夠濾除概率較大的椒鹽噪聲，而且能夠更好的保護圖像的細節，這是常規的中值濾波器做不到的。自適應的中值濾波器也需要一個矩形的窗口 ，和常規中值濾波器不同的是這個窗口的大小會在濾波處理的過程中進行改變（增大）。需要注意的是，濾波器的輸出是一個像素值，該值用來替換點 處的像素值，點 是濾波窗口的中心位置。

在描述自適應中值濾波器時需要用到如下的符號：

• 窗口中的最小灰度值
• 窗口中的最大灰度值
• 窗口中的灰度值的中值
• 表示坐標 處的灰度值
• 允許的最大窗口尺寸

自適應中值濾波器有兩個處理過程，分別記為：和。

A :

 

 

如果A1 > 0 且 A2 < 0，跳轉到 B；

否則，增大窗口的尺寸 如果增大后窗口的尺寸 ，則重復A過程。否則，輸出 𝑍𝑚𝑒𝑑

B:

 

 

如果B1 > 0 且 B2 <0則輸出 ,否則輸出 

自適應中值濾波原理說明

過程A的目的是確定當前窗口內得到中值 𝑍𝑚𝑒𝑑 是否是噪聲。如果 𝑍𝑚𝑖𝑛<𝑍𝑚𝑒𝑑<𝑍𝑚𝑎𝑥，則𝑍𝑚𝑒𝑑不是噪聲，這時轉到過程B測試當前窗口的中心位置的像素 𝑍𝑥𝑦 是否是一個噪聲點。如果 𝑍𝑚𝑖𝑛<𝑍𝑥𝑦<𝑍𝑚𝑎𝑥，則 𝑍𝑥𝑦 不是一個噪聲，此時輸出𝑍𝑥𝑦 ;如果不滿足上述條件，則可判定 𝑍𝑥𝑦 是噪聲，這是輸出中值 𝑍𝑚𝑒𝑑 （在A中已經判斷出 𝑍𝑚𝑒𝑑 不是噪聲）。

如果在過程A中，得到則 𝑍𝑚𝑒𝑑 不符合條件 𝑍𝑚𝑖𝑛<𝑍𝑚𝑒𝑑<𝑍𝑚𝑎𝑥 ，則可判斷得到的中值 𝑍𝑚𝑒𝑑 是一個噪聲。在這種情況下，需要增大濾波器的窗口尺寸，在一個更大的范圍內尋找一個非噪聲點的中值，直到找到一個非噪聲的中值，跳轉到B；或者，窗口的尺寸達到了最大值，這時返回找到的中值，退出。

從上面分析可知，噪聲出現的概率較低，自適應中值濾波器可以較快的得出結果，不需要去增加窗口的尺寸；反之，噪聲的出現的概率較高，則需要增大濾波器的窗口尺寸，這也符合種中值濾波器的特點：噪聲點比較多時，需要更大的濾波器窗口尺寸。

算法實現

有了算法的詳細描述，借助於OpenCV對圖像的讀寫，自適應中值濾波器實現起來也不是很困難。首先定義濾波器最小的窗口尺寸以及最大的窗口尺寸。要進行濾波處理，首先要擴展圖像的邊界，以便對圖像的邊界像素進行處理。copyMakeBorder根據選擇的BorderTypes使用不同的值擴充圖像的邊界像素，具體可參考OpenCV的文檔信息。下面就是遍歷圖像的像素，對每個像素進行濾波處理。需要注意一點，不論濾波器多么的復雜，其每次的濾波過程，都是值返回一個值，來替換掉當前窗口的中心的像素值。函數adpativeProcess就是對當前像素的濾波過程，其代碼如下：

uchar adaptiveProcess(const Mat &im, int row,int col,int kernelSize,int maxSize){ vector<uchar> pixels; for (int a = -kernelSize / 2; a <= kernelSize / 2; a++) for (int b = -kernelSize / 2; b <= kernelSize / 2; b++) { pixels.push_back(im.at<uchar>(row + a, col + b)); } sort(pixels.begin(), pixels.end()); auto min = pixels[0]; auto max = pixels[kernelSize * kernelSize - 1]; auto med = pixels[kernelSize * kernelSize / 2]; auto zxy = im.at<uchar>(row, col); if (med > min && med < max) { // to B if (zxy > min && zxy < max) return zxy; else return med; } else { kernelSize += 2; if (kernelSize <= maxSize) return adpativeProcess(im, row, col, kernelSize, maxSize); // 增大窗口尺寸，繼續A過程。 else return med; }}

有了上面這個函數，剩下的只需要對全部像素做一個遍歷即可，更為完整的代碼，請見我的Github地址：

https://github.com/zhangqizky/common-image-filteringgithub.com

高斯濾波

高斯濾波也是一種非常常見的濾波方法，其核的形式為：

其中是圖像中的點的坐標,  為標准差，高斯模板就是利用這個函數來計算的，x和y都是代表，以核中心點為坐標原點的坐標值。這里想說一下 的作用，當 比較小的時候，生成的高斯模板中心的系數比較大，而周圍的系數比較小，這樣對圖像的平滑效果不明顯。而當 比較大時，生成的模板的各個系數相差就不是很大，比較類似於均值模板，對圖像的平滑效果比較明顯。

高斯濾波沒有特別多可說的，最主要的作用是濾除高斯噪聲，即符合正態分布的噪聲。

實現的方式有兩種，第一種是按照公式暴力實現，代碼如下:

//O(m * n * ksize^2)void GaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma){ CV_Assert(src.channels() || src.channels() == 3); //只處理3通道或單通道的圖片 double **GaussianTemplate = new double *[ksize]; for(int i = 0; i < ksize; i++){ GaussianTemplate[i] = new double [ksize]; } generateGaussianTemplate(GaussianTemplate, ksize, sigma); //padding int border = ksize / 2; copyMakeBorder(src, dst, border, border, border, border, BORDER_CONSTANT); int channels = dst.channels(); int rows = dst.rows - border; int cols = dst.cols - border; for(int i = border; i < rows; i++){ for(int j = border; j< cols; j++){ double sum[3] = {0}; for(int a = -border; a <= border; a++){ for(int b = -border; b <= border; b++){ if(channels == 1){ sum[0] += GaussianTemplate[border+a][border+b] * dst.at<uchar>(i+a, j+b); }else if(channels == 3){ Vec3b rgb = dst.at<Vec3b>(i+a, j+b); auto k = GaussianTemplate[border+a][border+b]; sum[0] += k * rgb[0]; sum[1] += k * rgb[1]; sum[2] += k * rgb[2]; } } } for(int k = 0; k < channels; k++){ if(sum[k] < 0) sum[k] = 0; else if(sum[k] > 255) sum[k] = 255; } if(channels == 1){ dst.at<uchar >(i, j) = static_cast<uchar >(sum[0]); }else if(channels == 3){ Vec3b rgb = {static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2])}; dst.at<Vec3b>(i, j) = rgb; } } } for(int i = 0; i < ksize; i++) delete[] GaussianTemplate[i]; delete[] GaussianTemplate;}

當核比較大時，高斯濾波會比較費時，此時可以使用分離X和Y通道的形式來實現可分離的高斯濾波。為什么可以這么做？因為高斯函數中沒有這樣的耦合項，即x和y是相對獨立的，此時就可以將兩個維度分離開來。

//分離實現高斯濾波//O(m*n*k)void separateGaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma){ assert(src.channels()==1 || src.channels() == 3); //只處理單通道或者三通道圖像 //生成一維的 double *matrix = new double[ksize]; double sum = 0; int origin = ksize / 2; for(int i = 0; i < ksize; i++){ double g = exp(-(i-origin) * (i-origin) / (2 * sigma * sigma)); sum += g; matrix[i] = g; } for(int i = 0; i < ksize; i++) matrix[i] /= sum; int border = ksize / 2; copyMakeBorder(src, dst, border, border, border, border, BORDER_CONSTANT); int channels = dst.channels(); int rows = dst.rows - border; int cols = dst.cols - border; //水平方向 for(int i = border; i < rows; i++){ for(int j = border; j < cols; j++){ double sum[3] = {0}; for(int k = -border; k<=border; k++){ if(channels == 1){ sum[0] += matrix[border + k] * dst.at<uchar>(i, j+k); }else if(channels == 3){ Vec3b rgb = dst.at<Vec3b>(i, j+k); sum[0] += matrix[border+k] * rgb[0]; sum[1] += matrix[border+k] * rgb[1]; sum[2] += matrix[border+k] * rgb[2]; } } for(int k = 0; k < channels; k++){ if(sum[k] < 0) sum[k] = 0; else if(sum[k] > 255) sum[k] = 255; } if(channels == 1) dst.at<Vec3b>(i, j) = static_cast<uchar>(sum[0]); else if(channels == 3){ Vec3b rgb = {static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2])}; dst.at<Vec3b>(i, j) = rgb; } } } //豎直方向 for(int i = border; i < rows; i++){ for(int j = border; j < cols; j++){ double sum[3] = {0}; for(int k = -border; k<=border; k++){ if(channels == 1){ sum[0] += matrix[border + k] * dst.at<uchar>(i+k, j); }else if(channels == 3){ Vec3b rgb = dst.at<Vec3b>(i+k, j); sum[0] += matrix[border+k] * rgb[0]; sum[1] += matrix[border+k] * rgb[1]; sum[2] += matrix[border+k] * rgb[2]; } } for(int k = 0; k < channels; k++){ if(sum[k] < 0) sum[k] = 0; else if(sum[k] > 255) sum[k] = 255; } if(channels == 1) dst.at<Vec3b>(i, j) = static_cast<uchar>(sum[0]); else if(channels == 3){ Vec3b rgb = {static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2])}; dst.at<Vec3b>(i, j) = rgb; } } } delete [] matrix;}

同樣，完整的代碼請見：

https://github.com/zhangqizky/common-image-filtering/tree/maingithub.com

雙邊濾波

雙邊濾波是一種非線性濾波方法，是結合了圖像的鄰近度和像素值相似度的一種折中，在濾除噪聲的同時可以保留原圖的邊緣信息。整個雙邊濾波是由兩個函數構成：一個函數是由空間距離決定的濾波器系數，另外一個詩由像素差值決定的濾波器系數。整個雙邊濾波的公式如下：

其中權重系數 取決於定義域核：

和值域核

的乘積。其中定義域核影響的是空間位置，如果把圖像看成一個二維函數，那么定義域就是圖像的坐標，值域就是該坐標處對應的像素值。定義域核就是普通的高斯核，全局使用一個就可以。但值域核是需要對每個像素點滑動進行計算的。

那么如何理解雙邊濾波呢

高斯濾波的濾波核的意義是，濾波后的像素值等於窗口內的像素值的加權平均值，權值系數是符合高斯分布，距離該點越近，權值越大。但是沒有考慮像素值與當前點的差距。現在加上值域核，意義就在，濾波后當前點的像素值還會受到領域內像素值與自身的像素值差異的影響，不僅僅是距離來決定。這樣，在平緩的區域里，由於像素值差異非常小，則值域的權重趨向於1，所以雙邊濾波就近似為高斯濾波。而在邊緣區域中，由於像素值的差異比較大，則值域核趨向於0，權重下降，即當前像素受到領域內像素影響比較小，從而保留了邊緣信息。

雙邊濾波的代碼

opencv中提供了bilateralFilter()函數來實現雙邊濾波操作，其原型如下：

void cv::bilateralFilter(InputArray src,
OutputArray 	dst,
int 	d,
double 	sigmaColor,
double 	sigmaSpace,
int 	borderType = BORDER_DEFAULT 
)		

• InputArray src: 輸入圖像，可以是Mat類型，圖像必須是8位整型或浮點型單通道、三通道的圖像。
• OutputArray dst: 輸出圖像，和原圖像有相同的尺寸和類型。
• int d: 表示在過濾過程中每個像素鄰域的直徑范圍。如果這個值是非正數，則函數會從第五個參數sigmaSpace計算該值。
• double sigmaColor: 顏色空間過濾器的值，這個參數的值月大，表明該像素鄰域內有越寬廣的顏色會被混合到一起，產生較大的半相等顏色區域。（這個參數可以理解為值域核的 和 )
• double sigmaSpace: 坐標空間中濾波器的sigma值，如果該值較大，則意味着越遠的像素將相互影響，從而使更大的區域中足夠相似的顏色獲取相同的顏色。當d>0時，d指定了鄰域大小且與sigmaSpace無關，否則d正比於sigmaSpace. （這個參數可以理解為空間域核的 和 ）
• int borderType=BORDER_DEFAULT: 用於推斷圖像外部像素的某種邊界模式，有默認值BORDER_DEFAULT.

具體代碼如下：

#include <iostream>#include <opencv2/opencv.hpp>
using namespace std;using namespace cv;
//定義全局變量const int g_ndMaxValue = 100;const int g_nsigmaColorMaxValue = 200;const int g_nsigmaSpaceMaxValue = 200;int g_ndValue;int g_nsigmaColorValue;int g_nsigmaSpaceValue;
Mat g_srcImage;Mat g_dstImage;
//定義回調函數void on_bilateralFilterTrackbar(int, void*);
int main(){ g_srcImage = imread("lena.jpg");
 //判斷圖像是否加載成功 if(g_srcImage.empty()) { cout << "圖像加載失敗!" << endl; return -1; } else cout << "圖像加載成功!" << endl << endl;
 namedWindow("原圖像", WINDOW_AUTOSIZE); imshow("原圖像", g_srcImage);
 //定義輸出圖像窗口屬性和軌跡條屬性 namedWindow("雙邊濾波圖像", WINDOW_AUTOSIZE); g_ndValue = 10; g_nsigmaColorValue = 10; g_nsigmaSpaceValue = 10;
 char dName[20]; sprintf(dName, "鄰域直徑 %d", g_ndMaxValue);
 char sigmaColorName[20]; sprintf(sigmaColorName, "sigmaColor %d", g_nsigmaColorMaxValue);
 char sigmaSpaceName[20]; sprintf(sigmaSpaceName, "sigmaSpace %d", g_nsigmaSpaceMaxValue);
 //創建軌跡條 createTrackbar(dName, "雙邊濾波圖像", &g_ndValue, g_ndMaxValue, on_bilateralFilterTrackbar); on_bilateralFilterTrackbar(g_ndValue, 0);
 createTrackbar(sigmaColorName, "雙邊濾波圖像", &g_nsigmaColorValue, g_nsigmaColorMaxValue, on_bilateralFilterTrackbar); on_bilateralFilterTrackbar(g_nsigmaColorValue, 0);
 createTrackbar(sigmaSpaceName, "雙邊濾波圖像", &g_nsigmaSpaceValue, g_nsigmaSpaceMaxValue, on_bilateralFilterTrackbar); on_bilateralFilterTrackbar(g_nsigmaSpaceValue, 0);
 waitKey(0);
 return 0;}
void on_bilateralFilterTrackbar(int, void*){ bilateralFilter(g_srcImage, g_dstImage, g_ndValue, g_nsigmaColorValue, g_nsigmaSpaceValue); imshow("雙邊濾波圖像", g_dstImage);}

導向濾波(Guide Filter)

需要有高斯濾波和雙邊濾波的相關知識背景才能更好的理解導向濾波。在導向濾波中，首先利用了局部線性模型。這個模型認為某函數上一點與其近鄰部分的點成線性關系，一個復雜的函數就可以用很多局部的線性函數來表示，當需要求該函數上某一點的值時，只需要計算所有包含該點的線性函數的值並取平均值即可。這種模型，在表示非解析函數上，非常有用。

同理，我們可以認為圖像是一個二維函數，並且假設該函數的輸出與輸入在一個二維窗口內滿足線性關系，如下：

其中，是輸出像素的值，是輸入圖像的值，和是像素索引，和是當窗口中心位於k時該線性函數的系數。其實，輸入圖像不一定是待濾波的圖像本身，也可以是其他圖像即引導圖像，這也是為何稱為引導濾波的原因。對上式兩邊取梯度，可以得到:

即當輸入圖像有梯度時，輸出也有類似的梯度，現在可以解釋為什么引導濾波有邊緣保持特性了。下一步是求出線性函數的系數，也就是線性回歸，即希望擬合函數的輸出值與真實值之間的差距最小，也就是讓下式最小:

這里只能是待濾波圖像，並不像那樣可以是其他圖像。同時，之前的系數用於防止求得的過大，也是調節濾波器濾波效果的重要參數(相當於L2正則化的權重懲罰)。接下來利用最小二乘法的原理令  和 得到2個二元一次方程，求解得到：

，

其中 是 在窗口的平均值， 是 在窗口 的方差， 是窗口 中的像素個數，  是待濾波圖像在窗口 中的均值。在計算每個窗口的線性系數時，我們可以發現一個像素會被多個窗口包含，也就是說，每個像素都由多個線性函數所描述。因此，如之前所說，要具體求某一點的輸出值時，只需將所有包含該點的線性函數值平均即可，如下:

這里， 是所有包含像素 的窗口， 是其中心位置。

當把引導濾波用作邊緣保持濾波器時，往往有  ，如果 ，顯然是為最小值的解，從上式可以看出，這時的濾波器沒有任何作用，將輸入原封不動的輸出。如果 ，在像素強度變化小的區域（或單色區域），有近似於（或等於0，而近似於（或等於） ，即做了一個加權均值濾波；而在變化大的區域，近似於1，近似於0，對圖像的濾波效果很弱，有助於保持邊緣。而 的作用就是界定什么是變化大，什么是變化小。在窗口大小不變的情況下，隨着的增大，濾波效果越明顯。

在濾波效果上，引導濾波和雙邊濾波差不多，然后在一些細節上，引導濾波較好(在PS的磨皮美白中，經過實踐，效果更好)。引導濾波最大的優勢在於，可以寫出時間復雜度與窗口大小無關的算法，因此在使用大窗口處理圖片時，其效率更高。

同樣，OpenCV中也有導向濾波的接口。具體代碼如下：

void cv::ximgproc::guidedFilter	(	InputArray 	guide,
InputArray 	src,
OutputArray 	dst,
int 	radius,
double 	eps,
int 	dDepth = -1 
)	

guide	引導圖像，3通道，如果大於3通道則只有前三個會被用到
src	待濾波圖像
dst	輸出圖像
radius	導向濾波的窗口
eps	正則化參數
dDepth	可選，圖像的深度參數

這邊有個基於scipy實現的python代碼，可以參考一下：

import numpy as npimport scipy as spimport scipy.ndimage

def box(img, r): """ O(1) box filter img - >= 2d image r - radius of box filter """ (rows, cols) = img.shape[:2] imDst = np.zeros_like(img)

 tile = [1] * img.ndim tile[0] = r imCum = np.cumsum(img, 0) imDst[0:r+1, :, ...] = imCum[r:2*r+1, :, ...] imDst[r+1:rows-r, :, ...] = imCum[2*r+1:rows, :, ...] - imCum[0:rows-2*r-1, :, ...] imDst[rows-r:rows, :, ...] = np.tile(imCum[rows-1:rows, :, ...], tile) - imCum[rows-2*r-1:rows-r-1, :, ...]
 tile = [1] * img.ndim tile[1] = r imCum = np.cumsum(imDst, 1) imDst[:, 0:r+1, ...] = imCum[:, r:2*r+1, ...] imDst[:, r+1:cols-r, ...] = imCum[:, 2*r+1 : cols, ...] - imCum[:, 0 : cols-2*r-1, ...] imDst[:, cols-r: cols, ...] = np.tile(imCum[:, cols-1:cols, ...], tile) - imCum[:, cols-2*r-1 : cols-r-1, ...]
 return imDst
def _gf_color(I, p, r, eps, s=None): """ Color guided filter I - guide image (rgb) p - filtering input (single channel) r - window radius eps - regularization (roughly, variance of non-edge noise) s - subsampling factor for fast guided filter """ fullI = I fullP = p if s is not None: I = sp.ndimage.zoom(fullI, [1/s, 1/s, 1], order=1) p = sp.ndimage.zoom(fullP, [1/s, 1/s], order=1) r = round(r / s)
 h, w = p.shape[:2] N = box(np.ones((h, w)), r)
 mI_r = box(I[:,:,0], r) / N mI_g = box(I[:,:,1], r) / N mI_b = box(I[:,:,2], r) / N
 mP = box(p, r) / N
 # mean of I * p mIp_r = box(I[:,:,0]*p, r) / N mIp_g = box(I[:,:,1]*p, r) / N mIp_b = box(I[:,:,2]*p, r) / N
 # per-patch covariance of (I, p) covIp_r = mIp_r - mI_r * mP covIp_g = mIp_g - mI_g * mP covIp_b = mIp_b - mI_b * mP
 # symmetric covariance matrix of I in each patch: # rr rg rb # rg gg gb # rb gb bb var_I_rr = box(I[:,:,0] * I[:,:,0], r) / N - mI_r * mI_r; var_I_rg = box(I[:,:,0] * I[:,:,1], r) / N - mI_r * mI_g; var_I_rb = box(I[:,:,0] * I[:,:,2], r) / N - mI_r * mI_b;
 var_I_gg = box(I[:,:,1] * I[:,:,1], r) / N - mI_g * mI_g; var_I_gb = box(I[:,:,1] * I[:,:,2], r) / N - mI_g * mI_b;
 var_I_bb = box(I[:,:,2] * I[:,:,2], r) / N - mI_b * mI_b;
 a = np.zeros((h, w, 3)) for i in range(h): for j in range(w): sig = np.array([ [var_I_rr[i,j], var_I_rg[i,j], var_I_rb[i,j]], [var_I_rg[i,j], var_I_gg[i,j], var_I_gb[i,j]], [var_I_rb[i,j], var_I_gb[i,j], var_I_bb[i,j]] ]) covIp = np.array([covIp_r[i,j], covIp_g[i,j], covIp_b[i,j]]) a[i,j,:] = np.linalg.solve(sig + eps * np.eye(3), covIp)
 b = mP - a[:,:,0] * mI_r - a[:,:,1] * mI_g - a[:,:,2] * mI_b
 meanA = box(a, r) / N[...,np.newaxis] meanB = box(b, r) / N
 if s is not None: meanA = sp.ndimage.zoom(meanA, [s, s, 1], order=1) meanB = sp.ndimage.zoom(meanB, [s, s], order=1)
 q = np.sum(meanA * fullI, axis=2) + meanB
 return q

def _gf_gray(I, p, r, eps, s=None): """ grayscale (fast) guided filter I - guide image (1 channel) p - filter input (1 channel) r - window raidus eps - regularization (roughly, allowable variance of non-edge noise) s - subsampling factor for fast guided filter """ if s is not None: Isub = sp.ndimage.zoom(I, 1/s, order=1) Psub = sp.ndimage.zoom(p, 1/s, order=1) r = round(r / s) else: Isub = I Psub = p

 (rows, cols) = Isub.shape
 N = box(np.ones([rows, cols]), r)
 meanI = box(Isub, r) / N meanP = box(Psub, r) / N corrI = box(Isub * Isub, r) / N corrIp = box(Isub * Psub, r) / N varI = corrI - meanI * meanI covIp = corrIp - meanI * meanP

 a = covIp / (varI + eps) b = meanP - a * meanI
 meanA = box(a, r) / N meanB = box(b, r) / N
 if s is not None: meanA = sp.ndimage.zoom(meanA, s, order=1) meanB = sp.ndimage.zoom(meanB, s, order=1)
 q = meanA * I + meanB return q

def _gf_colorgray(I, p, r, eps, s=None): """ automatically choose color or gray guided filter based on I's shape """ if I.ndim == 2 or I.shape[2] == 1: return _gf_gray(I, p, r, eps, s) elif I.ndim == 3 and I.shape[2] == 3: return _gf_color(I, p, r, eps, s) else: print("Invalid guide dimensions:", I.shape)

def guided_filter(I, p, r, eps, s=None): """ run a guided filter per-channel on filtering input p I - guide image (1 or 3 channel) p - filter input (n channel) r - window raidus eps - regularization (roughly, allowable variance of non-edge noise) s - subsampling factor for fast guided filter """ if p.ndim == 2: p3 = p[:,:,np.newaxis]
 out = np.zeros_like(p3) for ch in range(p3.shape[2]): out[:,:,ch] = _gf_colorgray(I, p3[:,:,ch], r, eps, s) return np.squeeze(out) if p.ndim == 2 else out

def test_gf(): import imageio cat = imageio.imread('cat.bmp').astype(np.float32) / 255 tulips = imageio.imread('tulips.bmp').astype(np.float32) / 255
 r = 8 eps = 0.05
 cat_smoothed = guided_filter(cat, cat, r, eps) cat_detail = cat / cat_smoothed print(cat_detail.shape) cat_smoothed_s4 = guided_filter(cat, cat, r, eps, s=4)
 imageio.imwrite('cat_smoothed.png', cat_smoothed) imageio.imwrite('cat_smoothed_s4.png', cat_smoothed_s4) imageio.imwrite('cat_smoothed_detailed.png',cat_detail)
 tulips_smoothed4s = np.zeros_like(tulips) tulips_detailed = np.zeros_like(tulips) for i in range(3): tulips_smoothed4s[:,:,i] = guided_filter(tulips, tulips[:,:,i], r, eps, s=4)
 tulips_detailed = tulips / tulips_smoothed4s imageio.imwrite('tulips_detailed.png',tulips_detailed) imageio.imwrite('tulips_smoothed4s.png', tulips_smoothed4s)
 tulips_smoothed = np.zeros_like(tulips) for i in range(3): tulips_smoothed[:,:,i] = guided_filter(tulips, tulips[:,:,i], r, eps) imageio.imwrite('tulips_smoothed.png', tulips_smoothed)
if __name__ == '__main__': test_gf()

一副圖像，經過mask是圖像本身的導向濾波之后，得到一張細節圖和一張濾波圖。下面從左到右分別是原圖，細節圖和濾波圖。其實這是現在很多low-level領域的預處理步驟。拿到細節圖之后可以用卷積神經網絡做下面的處理。

 

 

 

這里還推薦一個很好的輪子，C++實現的。

https://github.com/atilimcetin/guided-filter

結語

以上就是常見的四種濾波算法的介紹。

 

推薦閱讀

• 一文圖解卡爾曼濾波(Kalman Filter)

• 消除Aliasing！加州大學&英偉達提出深度學習下采樣新思路：自適應低通濾波器層

• Transformer在CV領域有可能替代CNN嗎？還有哪些應用前景？

添加極市小助手微信（ID : cvmart2），備注：姓名-學校/公司-研究方向-城市（如：小極-北大-目標檢測-深圳），即可申請加入極市目標檢測/圖像分割/工業檢測/人臉/醫學影像/3D/SLAM/自動駕駛/超分辨率/姿態估計/ReID/GAN/圖像增強/OCR/視頻理解等技術交流群：每月大咖直播分享、真實項目需求對接、求職內推、算法競賽、干貨資訊匯總、與 10000+來自港科大、北大、清華、中科院、CMU、騰訊、百度等名校名企視覺開發者互動交流~