【3.15】matlab 高斯濾波，快速濾波

Gauss濾波快速實現方法（轉）

　　二維高斯函數具備旋轉對稱性，處理後不會對哪個方向上的邊緣進行了過多的濾波，所以相對其餘濾波器，具備沒法比擬的優越性。可是傳統Gauss濾波隨着圖像尺寸的增長，運算複雜度呈平方上漲，所以須要對其優化改進。下面，分別介紹傳統型，分解型和遞歸迭代型三種實現方法。

　　

1 傳統型

　　Gauss濾波首先須要構建一個Gauss濾波核，公式爲：

Matlab實現代碼：

dSigma =0.8;
fK1=1.0/(2*dSigma*dSigma);
fK2=fK1/pi;
iSize = 5;
step = floor(iSize/2 + 0.5);
for i = 1 : iSize
    x=i-step;
    fTemp=fK2*exp(-x*x*fK1);
    for j = 1 : iSize
        y=j-step;
        model(x+step,y+step)=fTemp*exp(-y*y*fK1);
    end
end
dSum = sum(sum(model));
model = model/dSum;                     %Gauss核數值歸一化

　　

　　

圖1 Gauss濾波核（5*5大小）

　　接下來就是將輸入圖像和濾波核進行卷積操做。其實質就是對原始圖像進行加權求和，把這個「和」賦給中心像素。對於一個2048*2048的圖像，須要進行104734756次乘法運算，和104734756次加法運算，運算複雜度是很高的。

2 分解型

　　

　　咱們能夠把一個二維Gauss核分解爲兩個一維高斯核，而後先對行作一次一維卷積，再對這個卷積結果作一次一維列卷積，獲得的結果徹底如出一轍，而開銷會小不少。

一維高斯核函數:

　　Matlab代碼實現：

dSigma =0.8;
fK1=1.0/(2*dSigma*dSigma);
fK2=fK1/pi;
iSize = 5;
step = floor(iSize/2 + 0.5);
for i = 1 : iSize
    x=i-step;
    fTemp=fK2*exp(-x*x*fK1);
    model(1,x+step) = fTemp;
end
dM = sum(model);
model = model / dM;

　　

圖2 一維高斯濾波核（1*5大小）

　　一維卷積原理和二維卷積同樣，只不過咱們只須要將同一行或同一列上的數據，按位置一一加權求和，再把「和」賦給中心元素。

　　對於一個2048*2048的圖像，須要進行41918464次乘法運算，和41918464次加法運算。相比傳統運算量，只是前者的1/2.4985。若是遇到頻繁計算Gauss濾波的算法，後者明顯比前者速度快不少。

3 遞歸迭代型

　　第二種方法較第一種方法，雖然有了較大改善，可是任然複雜度較高。這裏再介紹一種更快速的逼近Gauss濾波方法。

     具體步驟分爲兩步：首先對圖像作一次前向濾波，其次，對圖像再作一次後向濾波。

　　Forward：

　　Backward:

qFactor = 5;
b0Coeff = 1.57825 + (2.44413 * qFactor) + (1.4281 * qFactor * qFactor) + (0.422205 * qFactor * qFactor * qFactor);
b1Coeff = (2.44413 * qFactor) + (2.85619 * qFactor * qFactor) + (1.26661 * qFactor * qFactor * qFactor);
b2Coeff = (-1.4281 * qFactor * qFactor) + (-1.26661 * qFactor * qFactor * qFactor);
b3Coeff = 0.422205 * qFactor * qFactor * qFactor;
normalizationCoeff = 1 - ((b1Coeff + b2Coeff + b3Coeff) / b0Coeff);
vDenCoeff = [b0Coeff, -b1Coeff, -b2Coeff, -b3Coeff] / b0Coeff;
vXSignal = zeros(61, 1);
vXSignal(31) = 10;
vYSignal = filter(normalizationCoeff, vDenCoeff, vXSignal);
vYSignal = filter(normalizationCoeff, vDenCoeff, vYSignal(end:-1:1));
figure();
plot(vYSignal);

　　

參考資料：Recursive implementation of the Gaussian filter。 Ian T. Young,1995

　　　　　http://dsp.stackexchange.com/questions/22075/recursive-implementation-of-the-gaussian-filter