经典的同态滤波算法的优化及其应用参数配置。

　　同态滤波，网络上有不少文章提到过这个算法，咱们摘取百度的一段文字简要的说明了该算法的核心： 同态滤波是一种减小低频增长高频，从而减小光照变化并锐化边缘或细节的图像滤波方法。

　　关于该算法，网络上已经有不少资料了，也有不少给出了参考代码，可是很痛心的是我看到的没有一个是彻底正确的，或多或少都存在瑕疵，有些虽然算法最后的效果是差很少正确的，其实是和真正的算法是背道而驰的。

　　咱们在这里只有简单的语句来描述下该算法的过程。

       对于一幅图像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，γ_H >1，控制滤波器幅度的范围。H_hp一般为高通滤波器，如高斯（Gaussian）高通滤波器、巴特沃兹(Butterworth)高通滤波器、Laplacian滤波器等。

       若是H_hp采用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