ガウスフィルタリング/ガウス平滑化/ガウスファジイの実現とその高速アルゴリズム(Gaussian Filter,Gaussian Smooth,Gaussian Blur,Fast implementation)
ネット上では画像に対してガウスボケを行う文章が少なくなく、その原理は簡単で、ここではあまり紹介しません.ここでは、ガウスブラーを実現するためのいくつかのアルゴリズム(画像サイズがM*Nであると仮定し、filterの半径がrであることに注意し、多くの文章で使用されている用語はfilter sizeであり、半径サイズまたは直径サイズを指し、一般的にfilterの半径rは3倍または4倍のsigmaをとる):
1.最も原始的な実現方法は、ガウス関数を用いてガウステンプレートを生成し、その後、このテンプレートを用いて画像をボリューム化することである.この方法のアルゴリズム複雑度はO(M*N*r^2)であり,すなわち,各画素に対して複雑度はO(r^2)であり,フィルタサイズとquadratic関係にある.フィルタsizeを大きくすると、アルゴリズムは著しく遅くなります.
次はMatlab実装のコードです(コードの便宜上、以下のコードに現れるfilter sizeは半径ではなくフィルタ直径を指します!):
2.ガウスフィルタのkernelは分離可能であり、すなわち、2 Dのガウスkernelを2つの1 Dのkernelに分解して、まずx方向に画像を1 Dのガウスkernelの畳み込みを行い、次にy方向に画像を1 Dのガウスkernelの畳み込みを行い、最後の結果は1つの2 Dのガウスkernelを用いて画像の畳み込み効果と同じである.これにより、フィルタのアルゴリズム複雑度は、画素毎にO(r)に低下する.
3.FFT-based convolutionを使用します.画像のconvolutionはFFT変換後の周波数領域の積に等価であるため、ガウスkernelを画像と同程度の大きさに拡張し、画像とガウスkernelをそれぞれFFT変換し、そのまま両者のFFT結果画像をper pixelで乗算し、最終結果を逆FFT変換することができる.このアルゴリズムの全体的な複雑度はO(M*N*log(M*N))であり,フィルタの複雑度は画素毎にO(log(M*N))である.フィルタsizeが比較的小さい場合、このアルゴリズムは分離可能な1 D Gaussボリュームを使用するよりも、フィルタsizeが比較的大きい(r>log(M*N))場合、FFT−basedボリューム速度が比較的速いことがわかる.Matlabを用いてFFTベースのガウスファジイを実現し,ここでは非常に詳細な文章を説明する.
Matlabコードは以下の通りです(コードの便宜上、以下のコードに現れるfilter sizeは半径ではなくフィルタ直径を指します!):
4.再帰Gaussフィルタリングアルゴリズム.以上のアルゴリズムよりもさらに,approximate Gaussフィルタリングのアルゴリズムが提案されており,アルゴリズム速度をフィルタsizeに依存させず,各画素に対してフィルタのアルゴリズム複雑度はO(1)レベルである.この方面のアルゴリズム原理の紹介は比較的複雑で,信号処理に関する知識を理解しなければ理解できない(例えばz変換など,以下の参考文献を参照).幸いなことに、アルゴリズム自体はかなり簡単で、大体の構想は順方向にスキャンして画素値を蓄積し(蓄積時に特定の係数を使ってapproximateガウスkernelを使うことができます)、それから逆方向にスキャンして画素値を蓄積して、done!しかもアルゴリズムもseparableなので、x方向とy方向が別々にできるので、速度が非常に速くなります!関連するコードはここにあり、NVIDIAにもCUDAコードはここにあり、最新のGPU高速アルゴリズムはここを参照してください.原理はこことここを参考にすることができます.以下はいくつかのこの方面の参考文献である.
[1] Rachid Deriche - "Recursively implementing the Gaussian and its derivatives", 1993. [2] Lucas J. van Vliet, Ian T. Young and Piet W. Verbeek - "Recursive Gaussian derivative Filters", 1998 [3] Dave Hale, "Recursive Gaussian Filters", CWP-546 [4] Luis Alvarez, Rachid Deriche, Francisco Santana - "Recursivity and PDE's in Image Processing", 2000
1.最も原始的な実現方法は、ガウス関数を用いてガウステンプレートを生成し、その後、このテンプレートを用いて画像をボリューム化することである.この方法のアルゴリズム複雑度はO(M*N*r^2)であり,すなわち,各画素に対して複雑度はO(r^2)であり,フィルタサイズとquadratic関係にある.フィルタsizeを大きくすると、アルゴリズムは著しく遅くなります.
次はMatlab実装のコードです(コードの便宜上、以下のコードに現れるfilter sizeは半径ではなくフィルタ直径を指します!):
function [T] = GenGaussFilterKernel(size, sigma)
% here the size is the diameter of the filter
if mod(size, 2) == 0
size = size + 1;
end
T = zeros(size, size);
distMat = zeros(size, size);
centerIdx = ceil(size / 2);
for i = 1:size
for j = 1:size
distMat(i,j) = sqrt(abs(i - centerIdx) ^ 2 + abs(j - centerIdx) ^ 2);
end
end
T = exp(-distMat.^2 / (2 * sigma^2));
function [img_ext] = PaddingImgByMirrorEdge(img, r)
[h,w] = size(img);
% make a larger image, mirror the edges
img_ext = zeros(h + 2*r, w + 2*r);
img_ext(r+1:h+r, r+1:w+r) = img;
img_ext(1+r:h+r, 1:r) = flipdim(img(1:h, 2:r+1), 2); %add left border
img_ext(1+r:h+r, w+r+1:w+r+r) = flipdim(img(1:h, w-r:w-1), 2); %add right border
img_ext(1:r, 1:w+r+r) = flipdim(img_ext(r+2:r+r+1, 1:w+r+r), 1); %add up border
img_ext(h+r+1:h+r+r, 1:w+r+r) = flipdim(img_ext(h+r-r:h+r-1, 1:w+r+r), 1); %add bottom border
function [I] = ImageConvolution(img, filterKernel)
[h,w] = size(img);
filterSize = size(filterKernel);
r = floor(filterSize(1) / 2);
I = zeros(h, w);
% make a larger image, mirror the border to extend it
imgExt = PaddingImgByMirrorEdge(img, r);
% convolution
sumFilter = sum(sum(filterKernel));
tempRegion = zeros(filterSize);
for i = 1 : h
for j = 1 : w
tempRegion = imgExt(i:i+2*r, j:j+2*r);
I(i,j) = sum(sum(double(tempRegion) .* filterKernel)) / sumFilter;
end
end
I = uint8(I);
2.ガウスフィルタのkernelは分離可能であり、すなわち、2 Dのガウスkernelを2つの1 Dのkernelに分解して、まずx方向に画像を1 Dのガウスkernelの畳み込みを行い、次にy方向に画像を1 Dのガウスkernelの畳み込みを行い、最後の結果は1つの2 Dのガウスkernelを用いて画像の畳み込み効果と同じである.これにより、フィルタのアルゴリズム複雑度は、画素毎にO(r)に低下する.
3.FFT-based convolutionを使用します.画像のconvolutionはFFT変換後の周波数領域の積に等価であるため、ガウスkernelを画像と同程度の大きさに拡張し、画像とガウスkernelをそれぞれFFT変換し、そのまま両者のFFT結果画像をper pixelで乗算し、最終結果を逆FFT変換することができる.このアルゴリズムの全体的な複雑度はO(M*N*log(M*N))であり,フィルタの複雑度は画素毎にO(log(M*N))である.フィルタsizeが比較的小さい場合、このアルゴリズムは分離可能な1 D Gaussボリュームを使用するよりも、フィルタsizeが比較的大きい(r>log(M*N))場合、FFT−basedボリューム速度が比較的速いことがわかる.Matlabを用いてFFTベースのガウスファジイを実現し,ここでは非常に詳細な文章を説明する.
Matlabコードは以下の通りです(コードの便宜上、以下のコードに現れるfilter sizeは半径ではなくフィルタ直径を指します!):
function [I] = GaussianSmooth(img, filterSize, sigma, is_use_fft)
% here the size is the diameter of the filter
if nargin < 4
is_use_fft = 1;
end
T = GenGaussFilterKernel(filterSize, sigma); %generate gaussian kernel
[h, w] = size(img);
[d1, d2] = size(T);
r = floor(d1 / 2);
r1 = floor(h / 2);
r2 = floor(w / 2);
if is_use_fft == 0
% use convolution to implement
I = ImageConvolution(img, T);
else
% use fft to implement Gaussian smooth
extImg = PaddingImgByMirrorEdge(img, r);
extT = zeros(size(extImg)); %extend filter kernel to the same size
extT(1:end-h+1, 1:end-w+1) = T; %just place the original kernel on the topleft
fftImg = fft2(extImg);
fftT = fft2(extT);
resFFT = fftImg .* fftT; %multiply the FFT result
resI = ifft2(resFFT); %inverse FFT
% calculate coefficients of Gaussian smooth
extCoef = ones(size(extImg));
fftCoef = fft2(extCoef);
resCoef = fftCoef .* fftT;
resC = ifft2(resCoef);
% output image
I = uint8(resI ./ resC);
I = I(2*r+1:2*r+h, 2*r+1:2*r+w);
end
4.再帰Gaussフィルタリングアルゴリズム.以上のアルゴリズムよりもさらに,approximate Gaussフィルタリングのアルゴリズムが提案されており,アルゴリズム速度をフィルタsizeに依存させず,各画素に対してフィルタのアルゴリズム複雑度はO(1)レベルである.この方面のアルゴリズム原理の紹介は比較的複雑で,信号処理に関する知識を理解しなければ理解できない(例えばz変換など,以下の参考文献を参照).幸いなことに、アルゴリズム自体はかなり簡単で、大体の構想は順方向にスキャンして画素値を蓄積し(蓄積時に特定の係数を使ってapproximateガウスkernelを使うことができます)、それから逆方向にスキャンして画素値を蓄積して、done!しかもアルゴリズムもseparableなので、x方向とy方向が別々にできるので、速度が非常に速くなります!関連するコードはここにあり、NVIDIAにもCUDAコードはここにあり、最新のGPU高速アルゴリズムはここを参照してください.原理はこことここを参考にすることができます.以下はいくつかのこの方面の参考文献である.
[1] Rachid Deriche - "Recursively implementing the Gaussian and its derivatives", 1993. [2] Lucas J. van Vliet, Ian T. Young and Piet W. Verbeek - "Recursive Gaussian derivative Filters", 1998 [3] Dave Hale, "Recursive Gaussian Filters", CWP-546 [4] Luis Alvarez, Rachid Deriche, Francisco Santana - "Recursivity and PDE's in Image Processing", 2000