本篇博文來自博主Imageshop，打賞或想要查閱更多內(nèi)容可以移步至Imageshop。

轉(zhuǎn)載自：https://www.cnblogs.com/Imageshop/p/9766056.html ??侵刪

　　同態(tài)濾波，網(wǎng)絡(luò)上有很多文章提到過這個算法，我們摘取百度的一段文字簡要的說明了該算法的核心：?同態(tài)濾波是一種減少低頻增加高頻，從而減少光照變化并銳化邊緣或細(xì)節(jié)的圖像濾波方法。

　　關(guān)于該算法，網(wǎng)絡(luò)上已經(jīng)有很多資料了，也有很多給出了參考代碼，但是很痛心的是我看到的沒有一個是完全正確的，或多或少都存在瑕疵，有些雖然算法最后的效果是差不多正確的，實際上是和真正的算法是背道而馳的。

　　我們在這里只有簡單的語句來描述下該算法的過程。

? ? ? ?對于一幅圖像f(x,y)，可以表示為照射分量i(x,y)和反射分量r(x,y)的乘積。其中0<i(x,y)<∞，0<r(x,y)<1。i(x,y)描述景物的照明，變化緩慢，處于低頻成分。r(x,y)描述景物的細(xì)節(jié)，變化較快，處于高頻成分。因為該性質(zhì)是乘性的，所以不能直接使用傅里葉變換對i(x,y)和r(x,y)進(jìn)行控制，因此可以先對f(x,y)取對數(shù)，分離i(x,y)和r(x,y)。令z(x,y) = ln f(x,y) = ln i(x,y) + ln r(x,y)。由于f(x,y)的取值范圍為[0, L-1]，為了避免出現(xiàn)ln(0)的情況，故采用ln ( f(x,y) + 1 ) 來計算。

? ? ? 然后取傅里葉變換，得到 Z(u,v) = Fi(u,v) + Fr(u,v)。

? ? ? ?然后使用一個濾波器，對Z(u,v)進(jìn)行濾波，有 S(u,v) = H(u,v) Z(u,v) = H(u,v)Fi(u,v) + H(u,v)Fr(u,v)。

? ? ? ?濾波后，進(jìn)行反傅里葉變換，有 s(x, y) = IDFT( S(u,v) )。

? ? ? ?最后，反對數(shù)（取指數(shù)），得到最后處理后的圖像。g(x,y) = exp^(s(x,y)) = i0(x,y)+r0(x,y)。由于我們之前使用ln ( f(x,y)+1)，因此此處使用exp^(s(x,y)) - 1。? i0(x,y)和r0(x,y)分別是處理后圖像的照射分量和入射分量。

? ? ? ?這個濾波器通常我們?nèi)∪缦滦问?#xff1a;

? ? ? ? ? ? ??

? ? ? ?其中，

? ? ? ? ? ?γL< 1，γH?>1，控制濾波器幅度的范圍。Hhp通常為高通濾波器，如高斯（Gaussian）高通濾波器、巴特沃茲(Butterworth)高通濾波器、Laplacian濾波器等。

? ? ? ?如果Hhp采用Gaussian高通濾波器，則有：

? ? ? ? ? ? ??

? ? ? ?其中，c為一個常數(shù)，控制濾波器的形態(tài)，即從低頻到高頻過渡段的陡度（斜率），其值越大，斜坡帶越陡峭，見下圖。

? ? ? ? ? ? ? ? ??

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?圖2 同態(tài)濾波器幅頻曲線

　　如果英文可以的，直接看http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/OWENS/LECT5/node4.html這篇文章。

　　其實過程很簡單，但是很多博文都給出了錯誤的代碼，比如在圖像增強(qiáng)處理之：同態(tài)濾波與Retinex算法(一)同態(tài)濾波一文中，其給出的主要代碼如下：

function I3 = homofilter(I) %同態(tài)濾波函數(shù) subplot(1,2,1),imshow(I);title('同態(tài)濾波之前原始圖像'); I=double(rgb2gray(I)); [M,N]=size(I); rL=0.5; rH=5;%可根據(jù)需要效果調(diào)整參數(shù) c=3; d0=9; I1=log(I+1);%取對數(shù) FI=fft2(I1);%傅里葉變換 n1=floor(M/2); n2=floor(N/2); for i=1:M for j=1:N D(i,j)=((i-n1).^2+(j-n2).^2); H(i,j)=(rH-rL).*(exp(c*(-D(i,j)./(d0^2))))+rL;%高斯同態(tài)濾波 end end I2=ifft2(H.*FI);%傅里葉逆變換 I3=real(exp(I2)); subplot(1,2,2),imshow(I3,[]);title('同態(tài)濾波增強(qiáng)后');

　　我們看看其傳遞函數(shù)H那一行的代碼：?

H(i,j)=(rH-rL).*(exp(c*(-D(i,j)./(d0^2))))+rL;%高斯同態(tài)濾波 很明顯，這里有著嚴(yán)重的錯誤，1-exp操作，他直接寫成了exp。雖然他最后得到了一個似乎不錯的結(jié)果，但是這是違背科學(xué)的。

　 ? 再比如同態(tài)濾波及其實現(xiàn)?這篇文章其實也是不對的，其部分代碼如下：

% 生成同態(tài)濾波函數(shù)，中心在(m+1,n+1) Homo = zeros(P, Q); a = D0^2; % 計算一些不變的中間參數(shù) r = rh-rl; for u = 1 : Pfor v = 1 : Qtemp = (u-(m+1.0))^2 + (v-(n+1.0))^2;Homo(u, v) = r * (1-exp((-c)*(temp/a))) + rl;end end%進(jìn)行濾波 G = F1 .* Homo;% 反傅里葉變換 gp = ifft2(G);

　 由于之前計算的fft結(jié)果沒有中心化，所以進(jìn)行濾波前需要對Homo進(jìn)行反中心化，所以這里的代碼也不對（這里是我錯了，他代碼前面有* (-1)^(i+j)，有這個的話后面不需要ifftshift了）。

?　我們在看matlab—同態(tài)濾波的實現(xiàn)這篇文章的代碼：

clear all clc I =imread('moon128.bmp'); % tun.bmp figure(1),subplot(221),imshow(I); title('原圖') I=im2double(I); %轉(zhuǎn)換數(shù)據(jù)類型為double型 [M,N]=size(I); P = 2*M; Q = 2*N; I2 = zeros(P,Q); for i = 1:Mfor j =1:NI2(i,j) = I(i,j); %對圖像進(jìn)行填充end end figure(1),subplot(222),imshow(I2);title('填充后的圖像') I2=log(I2+1); %取對數(shù) FI=fft2(I2); %傅里葉變換 rL=0.25; rH=2.2; % 可根據(jù)需要效果調(diào)整參數(shù) c=2.0; %銳化參數(shù) D0=20; n1=floor(P/2); n2=floor(Q/2); for u=1:P for v=1:Q D(u,v)=sqrt(((u-n1).^2+(v-n2).^2)); %頻率域中點（u，v）與頻率矩形中心的距離 H(u,v)=(rH-rL).*(1-exp(-c.*(D(u,v)^2./D0^2)))+rL; %高斯同態(tài)濾波 end end H=ifftshift(H); %對H做反中心化 I3=ifft2(H.*FI); %傅里葉逆變換 I4=real(I3); I5 =I4(1:M, 1:N); %截取一部分 I6=exp(I5)-1; %取指數(shù) figure(1),subplot(223),imshow(I6,[]);title('同態(tài)濾波增強(qiáng)后'）

　　個人覺得這個比較接近正確的結(jié)果了，但是我覺得還是有點問題，im2double函數(shù)會將圖像數(shù)據(jù)歸一化到【0,1】之間，和大部分的實現(xiàn)方式都不同。而且這個時候如果避免log后面參數(shù)為0，也不易加上1了，數(shù)量級太大了，加上0.0001之類的就可以了。

　　再找一篇文章的代碼：同態(tài)濾波用于光照不均勻校正，這里的貼的效果也不錯，代碼如下：

clear all;[filename,pathname]=uigetfile('*.*','選擇圖像文件');%通過瀏覽文件夾來讀取圖片 if isequal(filename,0) %判斷是否選擇msgbox('沒有選擇任何圖片'); elseimage=imread(strcat(pathname,filename));%獲取圖像路徑path，然后讀取圖片file end figure();subplot(2,2,1); imshow(abs(image)); title('原始圖像'); img=im2double(image); %轉(zhuǎn)換圖像矩陣為雙精度型 lnimg=log(img+0.000001); %取對數(shù) Fimg=fft2(lnimg); %傅里葉變換 P=fftshift(Fimg); %將頻域原點移到圖像中心； [M,N]=size(P); %返回的行數(shù)和列數(shù)在P作為單獨的輸出變量 subplot(2,2,2);imshow(uint8(abs(P)),[]);title('濾波前的頻譜圖像'); %顯示無符號8位數(shù)，即256級的灰度圖像 x0=floor(M/2); y0=floor(N/2);%表示將向量M和N每個元素與2作除法后取整 %同態(tài)濾波參數(shù)設(shè)置 D0=100;%截止頻率 c=1.50;%銳化系數(shù) Hh=0.5;Hl=2; %Hh<1,Hl>1，Hh為高頻增益，Hl為低頻增益， %通過改變這兩個參數(shù)，得到不同的濾波效果 for u=1:M for v=1:N D(u,v)=sqrt((u-x0)^2+(v-y0)^2);%點（u,v）到頻率平面原點的距離H(u,v)=(Hh-Hl)*(1-exp(-c*(D(u,v)^2/D0^2)))+Hl;%同態(tài)濾波器函數(shù)end end hImg=Fimg.*H(u,v);%濾波，矩陣點乘 Q=fftshift(hImg);%傅里葉逆變換 subplot(2,2,3),imshow(uint8(abs(Q))),title('濾波后的頻譜圖像') gImg=ifft2(hImg);%反傅立葉變換 Y=exp(gImg); %取指數(shù) J=im2uint8(Y);%轉(zhuǎn)換圖像矩陣為無符號8位數(shù)，即256級的灰度圖像 subplot(2,2,4),imshow(J),title(' 濾波后的增強(qiáng)圖像')

　　這里的代碼很明顯的Hh還比Hl小，我是在搞不懂作者這里是怎么考慮的，同樣的道理少了反中心化那一步。

　　通過以上代碼，我們發(fā)現(xiàn)有的代碼在做FFT變換前會將圖像擴(kuò)大一倍，比如這個matlab的網(wǎng)站中關(guān)于同態(tài)濾波也提到了2倍尺寸：https://blogs.mathworks.com/steve/2013/06/25/homomorphic-filtering-part-1/，雖然說得很有道理，但是經(jīng)過測試，我覺得真的有必要嗎，又占內(nèi)存又減低了速度。

　　綜合各篇文章的代碼，通過實驗我整理后的matlab代碼如下：

% 同態(tài)濾波器 % ImageIn - 需要進(jìn)行濾波的灰度圖像 % High - 高頻增益,需要大于1 % Low - 低頻增益,取值在0和1之間 % C - 銳化系數(shù) % Sigma - 截止頻率，越大圖像越亮 % 輸出為進(jìn)行濾波之后的灰度圖像 function [ImageOut] = HomoFilter(ImageIn, High, Low, C, Sigma)Img = double(ImageIn); % 轉(zhuǎn)換圖像矩陣為雙精度型,不會改變數(shù)據(jù)本身[Height, Width] = size(ImageIn); % 返回的行數(shù)和列數(shù)CenterX = floor(Width / 2); % 中心點坐標(biāo)CenterY = floor(Height / 2);LogImg = log(Img + 1); % 圖像對數(shù)數(shù)據(jù)Log_FFT = fft2(LogImg); % 傅里葉變換for Y = 1 : Height for X = 1 : Width Dist= (X - CenterX) * (X - CenterX) + (Y - CenterY) * (Y - CenterY); % 點（X,Y）到頻率平面原點的距離H(Y, X)=(High - Low) * (1 - exp(-C * (Dist / (2 * Sigma * Sigma)))) + Low; % 同態(tài)濾波器函數(shù)end endH = ifftshift(H); % 對H做反中心化 Log_FFT = H.* Log_FFT; % 濾波，矩陣點乘 Log_FFT = ifft2(Log_FFT); % 反傅立葉變換 Out = exp(Log_FFT)-1; % 取指數(shù) % 指數(shù)處理ge = exp(g)-1;% 歸一化到[0, L-1]Max = max(Out(:));Min = min(Out(:));Range = Max - Min;for Y = 1 : Height for X = 1 : Width ImageOut(Y, X) = uint8(255 * (Out(Y, X) - Min) / Range); endend end

　　以上代碼只針對灰度圖像。

　　第一，我們沒有將圖像擴(kuò)大2倍，實踐證明不需要。第二，我們沒有把圖像歸一化，如果使用歸一化的代碼，同樣的參數(shù)會有不同的效果。第三，必須對H做反中心化處理，如果不對H反中心化，就要對FFT的結(jié)果進(jìn)行中心化，此時代碼如下所示：

function [ImageOut] = HomoFilter(ImageIn, High, Low, C, Sigma)Img = double(ImageIn); % 轉(zhuǎn)換圖像矩陣為雙精度型,不會改變數(shù)據(jù)本身[Height, Width] = size(ImageIn); % 返回的行數(shù)和列數(shù)CenterX = floor(Width / 2); % 中心點坐標(biāo)CenterY = floor(Height / 2);LogImg = log(Img + 1); % 圖像對數(shù)數(shù)據(jù)Log_FFT = fft2(LogImg); % 傅里葉變換Log_FFT = fftshift(Log_FFT);for Y = 1 : Height for X = 1 : Width Dist= (X - CenterX) * (X - CenterX) + (Y - CenterY) * (Y - CenterY); % 點（X,Y）到頻率平面原點的距離H(Y, X)=(High - Low) * (1 - exp(-C * (Dist / (2 * Sigma * Sigma)))) + Low; % 同態(tài)濾波器函數(shù)end end Log_FFT = H.* Log_FFT; % 濾波，矩陣點乘Log_FFT = ifftshift(Log_FFT);Log_FFT = ifft2(Log_FFT); % 反傅立葉變換 Out = exp(Log_FFT)-1; % 取指數(shù) % 指數(shù)處理ge = exp(g)-1;% 歸一化到[0, L-1]Max = max(Out(:));Min = min(Out(:));Range = Max - Min;for Y = 1 : Height for X = 1 : Width ImageOut(Y, X) = uint8(255 * (Out(Y, X) - Min) / Range); endend end

　　算法使用場合及參數(shù)配置說明：

　　（1）、光照不均勻圖像的均勻化。

? ? ? ? ? ? ?? ? ? ?? ? ?

? ? ? ? ? ? ?? ?? ? ?

?　　要實現(xiàn)此種效果建議的參數(shù)配置如下：High = 2, Low =0.2, C > 1, ?50<Sigma < 200;

　　(2) 偏暗圖像的增強(qiáng)

? ? ????

? ? ?? ?

? ? ?要實現(xiàn)此種效果建議的參數(shù)配置如下：High = 2, Low =0.2, C = 0.1, ?Sigma = max(Width, Height)；

? ? ?以上是彩色的圖，彩色圖的處理方式有很多種，可以參考我以前所發(fā)的博文。

　 鑒于算法的這個特性，這個算法應(yīng)該也可以應(yīng)用在16位圖像的HDR顯示上，有空做做試驗。

? ? ? matlab的速度還是很慢的，我已經(jīng)用C++結(jié)合SSE優(yōu)化對上述過程進(jìn)行了編碼優(yōu)化，其中的FFT使用了opencv的代碼，目前opencv最新版本的FFT已經(jīng)支持任意長度的序列了，但是為了速度，一般還是調(diào)用GetOptimalDftSize獲取最佳的序列長度，然后用SSE優(yōu)化一下其他的輔助處理，速度上對800*600的灰度圖30ms左右。詳見附件的SSE_Optimization_Demo的enhance菜單。

　　

? ? ? ?Demo下載地址：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar

總結(jié)

以上是生活随笔為你收集整理的经典的同态滤波算法的优化及其应用参数配置的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。