同态滤波, 网络上有很多文章提到过这个算法, 我们摘取百度的一段文字简要的说明了该算法的核心: 同态滤波是一种减少低频增加高频, 从而减少光照变化并锐化边缘或细节的图像滤波方法.
关于该算法, 网络上已经有很多资料了, 也有很多给出了参考代码, 但是很痛心的是我看到的没有一个是完全正确的, 或多或少都存在瑕疵, 有些虽然算法最后的效果是差不多正确的, 实际上是和真正的算法是背道而驰的.
我们在这里只有简单的语句来描述下该算法的过程.
对于一幅图像 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, 控制滤波器幅度的范围. 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 进行反中心化, 所以这里的代码也不对.
我们在看 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
来源: https://www.cnblogs.com/Imageshop/p/9766056.html