經典的同態濾波算法的優化及其應用參數配置。


  同態濾波,網絡上有很多文章提到過這個算法,我們摘取百度的一段文字簡要的說明了該算法的核心: 同態濾波是一種減少低頻增加高頻,從而減少光照變化並銳化邊緣或細節的圖像濾波方法。

  關於該算法,網絡上已經有很多資料了,也有很多給出了參考代碼,但是很痛心的是我看到的沒有一個是完全正確的,或多或少都存在瑕疵,有些雖然算法最后的效果是差不多正確的,實際上是和真正的算法是背道而馳的。

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

       對於一幅圖像f(x,y),可以表示為照射分量i(x,y)和反射分量r(x,y)的乘積。其中0<i(x,y)<∞,0<r(x,y)<1。i(x,y)描述景物的照明,變化緩慢,處於低頻成分。r(x,y)描述景物的細節,變化較快,處於高頻成分。因為該性質是乘性的,所以不能直接使用傅里葉變換對i(x,y)和r(x,y)進行控制,因此可以先對f(x,y)取對數,分離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],為了避免出現ln(0)的情況,故采用ln ( f(x,y) + 1 ) 來計算。

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

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

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

       最后,反對數(取指數),得到最后處理后的圖像。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)分別是處理后圖像的照射分量和入射分量。

       這個濾波器通常我們取如下形式:

              

       其中,

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

       如果Hhp采用Gaussian高通濾波器,則有:

              

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

                  

                               圖2 同態濾波器幅頻曲線

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

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

function I3 = homofilter(I)    %同態濾波函數 subplot(1,2,1),imshow(I);title('同態濾波之前原始圖像'); I=double(rgb2gray(I)); [M,N]=size(I); rL=0.5; rH=5;%可根據需要效果調整參數 c=3; d0=9; I1=log(I+1);%取對數 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;%高斯同態濾波 end end I2=ifft2(H.*FI);%傅里葉逆變換 I3=real(exp(I2)); subplot(1,2,2),imshow(I3,[]);title('同態濾波增強后');  

  我們看看其傳遞函數H那一行的代碼: 

  H(i,j)=(rH-rL).*(exp(c*(-D(i,j)./(d0^2))))+rL;%高斯同態濾波    

很明顯,這里有着嚴重的錯誤,1-exp操作,他直接寫成了exp。雖然他最后得到了一個似乎不錯的結果,但是這是違背科學的。

    再比如同態濾波及其實現 這篇文章其實也是不對的,其部分代碼如下:

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

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

  我們在看matlab—同態濾波的實現這篇文章的代碼:

clear all clc I =imread('moon128.bmp'); % tun.bmp figure(1),subplot(221),imshow(I); title('原圖') I=im2double(I);    %轉換數據類型為double型 [M,N]=size(I); P = 2*M; Q = 2*N; I2 = zeros(P,Q); for i = 1:M for j =1:N I2(i,j) = I(i,j);  %對圖像進行填充 end end figure(1),subplot(222),imshow(I2);title('填充后的圖像') I2=log(I2+1);    %取對數 FI=fft2(I2);    %傅里葉變換 rL=0.25; rH=2.2;     % 可根據需要效果調整參數 c=2.0;       %銳化參數 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; %高斯同態濾波 end end H=ifftshift(H);  %對H做反中心化 I3=ifft2(H.*FI);  %傅里葉逆變換 I4=real(I3); I5 =I4(1:M, 1:N);  %截取一部分 I6=exp(I5)-1;  %取指數 figure(1),subplot(223),imshow(I6,[]);title('同態濾波增強后'

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

  再找一篇文章的代碼:同態濾波用於光照不均勻校正,這里的貼的效果也不錯,代碼如下:

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

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

  通過以上代碼,我們發現有的代碼在做FFT變換前會將圖像擴大一倍,比如這個matlab的網站中關於同態濾波也提到了2倍尺寸:https://blogs.mathworks.com/steve/2013/06/25/homomorphic-filtering-part-1/,雖然說得很有道理,但是經過測試,我覺得真的有必要嗎,又占內存又減低了速度。

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

% 同態濾波器 % ImageIn   - 需要進行濾波的灰度圖像 % High      - 高頻增益,需要大於1 % Low       - 低頻增益,取值在0和1之間 % C         - 銳化系數 % Sigma     - 截止頻率,越大圖像越亮 % 輸出為進行濾波之后的灰度圖像 function [ImageOut] = HomoFilter(ImageIn, High, Low, C, Sigma) Img = double(ImageIn);      % 轉換圖像矩陣為雙精度型,不會改變數據本身 [Height, Width] = size(ImageIn);      % 返回的行數和列數 CenterX = floor(Width / 2);     % 中心點坐標 CenterY = floor(Height / 2); LogImg = log(Img + 1);      % 圖像對數數據 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;      % 同態濾波器函數 end end H = ifftshift(H);                   % 對H做反中心化 Log_FFT = H.* Log_FFT;              % 濾波,矩陣點乘 Log_FFT = ifft2(Log_FFT);           % 反傅立葉變換 Out = exp(Log_FFT)-1;               % 取指數 % 指數處理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); end end end 

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

  第一,我們沒有將圖像擴大2倍,實踐證明不需要。第二,我們沒有把圖像歸一化,如果使用歸一化的代碼,同樣的參數會有不同的效果。第三,必須對H做反中心化處理,如果不對H反中心化,就要對FFT的結果進行中心化,此時代碼如下所示:

function [ImageOut] = HomoFilter(ImageIn, High, Low, C, Sigma) Img = double(ImageIn);      % 轉換圖像矩陣為雙精度型,不會改變數據本身 [Height, Width] = size(ImageIn);      % 返回的行數和列數 CenterX = floor(Width / 2);     % 中心點坐標 CenterY = floor(Height / 2); LogImg = log(Img + 1);      % 圖像對數數據 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;      % 同態濾波器函數 end end Log_FFT = H.* Log_FFT;              % 濾波,矩陣點乘 Log_FFT = ifftshift(Log_FFT); Log_FFT = ifft2(Log_FFT);           % 反傅立葉變換 Out = exp(Log_FFT)-1;               % 取指數 % 指數處理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); end end end 

  算法使用場合及參數配置說明:

  (1)、光照不均勻圖像的均勻化。

                         

                     

   要實現此種效果建議的參數配置如下:High = 2, Low =0.2, C > 1,  50<Sigma < 200;

  (2) 偏暗圖像的增強

        

        

     要實現此種效果建議的參數配置如下:High = 2, Low =0.2, C = 0.1,  Sigma = max(Width, Height);

     以上是彩色的圖,彩色圖的處理方式有很多種,可以參考我以前所發的博文。

  鑒於算法的這個特性,這個算法應該也可以應用在16位圖像的HDR顯示上,有空做做試驗。

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

  

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

      


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM