老马迷图分享 http://blog.sciencenet.cn/u/zmpenguestc 网络空间做伏枥老马,志在千里育识途小驹!

博文

图像的同态滤波原理及实现

已有 32941 次阅读 2017-4-30 13:09 |个人分类:闻图思学|系统分类:教学心得| 图像处理, 同态滤波

态滤波(Homomorphic filter)是信号与图像处理中的一种常用技术,它采用了一种线性滤波在不同域中的非线性映射。这一技术是上世纪60年代由麻省理工学院(MIT)Thomas StockhamAlan V. Oppenheim Ronald W. Schafer 等几位学者提出。

如果您对原理部分不感兴趣,可直接跳至双虚线以下的关键代码实现部分

首先,介绍两个概念。

同态系统:是将非线性问题,转化为线性问题处理。即对非线性(乘性)混杂信号,通过某种数学运算(如对数变换),变成加性模型,而后采用线性滤波方法进行处理。

同态滤波:是把频率滤波和空域灰度变换结合起来的一种图像处理方法,它根据图像的照度/反射率模型作为频域处理的基础,利用压缩亮度范围和增强对比度来改善图像的质量。

下面介绍关于同态滤波的基本原理。

一幅图像可看成由两部分组成,即


其中,fi代表随空间位置不同的光强(Illumination)分量,其特点是缓慢变化,集中在图像的低频部分。fr代表景物反射到人眼的反射(Reflectance)分量。其特点包含了景物各种信息,高频成分丰富。

同态滤波过程,分为以下5个基本步骤:

原图做对数变换,得到如下两个加性分量,即


对数图像做傅里叶变换,得到其对应的频域表示为:


设计一个频域滤波器H(u,v),进行对数图像的频域滤波。

傅里叶反变换,返回空域对数图像。

取指数,得空域滤波结果。

综上,同态滤波的基本步骤如图1所示。


图1 同态滤波的基本步骤

可以看出,同态滤波的关键在于滤波器H的设计。对于一幅光照不均匀的图像,同态滤波可同时实现亮度调整和对比度提升,从而改善图像质量。为了压制低频的亮度分量,增强高频的反射分量,滤波器H应是一个高通滤波器,但又不能完全cut off 低频分量,仅作适当压制。

因此,同态滤波器一般采用如下形式,即

其中,gL< 1, gH >1,控制滤波器幅度的范围。Hhp常为高通滤波器,如高斯(Gaussian)高通滤波器、巴特沃兹(Butterworth)高通滤波器、Laplacian滤波器等。

如果Hhp采用Gaussian高通滤波器,则有:

其中,c为一个常数,控制滤波器的形态,即从低频到高频过渡段的陡度(斜率),其值越大,斜坡带越陡峭,见图2。


图2 同态滤波器幅频曲线

================================================================================

同态滤波的原理就这么简单,不过写代码实现起来,具体问题就来了。你不妨百度一下,网上一些论坛里,贴有大量可参考的代码,让您看了像雾像雨又像云,很可能会误导你。那么,只有吃透了原理,或自己动手测试一下,方能做出正确判断。

实现对一幅输入图像同态滤波的关键(MATLAB)代码如下:

输入图像取对数,一般取自然对数,即

                           f = log(1+I);                          % log transformation

这里图像I做对数变换前,需要转化为double型。针对灰度范围[0~255]的图像数据,log(I+1)是为了满足真数大于0,以防计算无意义。特别提醒,如果是归一化图像数据,则建议log(I+0.01)

fFFT,此步涉及后续的频域滤波。因此,FFT前数据f一般需零填充至本身的2倍尺寸。

                           Fp = fft2(f,P,Q);                   % FFT with zeros padding

这里,P = 2M,  Q = 2N,  M, N为数据f的阵列大小。

滤波器H设计,滤波器尺寸同样应为PQ

                          [v, u] = meshgrid(1:Q, 1:P);  % [y, x] = size(I),y为行数,x为列数。与习惯顺序[x, y]有点反,易混淆!

                          u = u - floor(P/2);                  % u方向 centralization

                          v = v - floor(Q/2);                  % v方向 centralization

                          D = u.^2 + v.^2;                      % compute the distances

                          H = 1-exp(-c*(D./D0^2));    % gaussian high-pass filter

                          H = (rH - rL)*H + rL;           % homomorphic filter

这里,D0gaussian滤波器的截止频域。因产生的为中心化滤波器,下一步滤波f的傅里叶变换也需要中心化处理;否则,H反中心化。

频域滤波。

                           H   = ifftshift(H);          % H反中心化,因Fp未中心化!

                           Gp = Fp.*H;                   % filtering in frequency domain

傅里叶反变换,返回空域的对数图像。

                           gp  = real(ifft2(Gp));    % IFFT取实部。网上有些代码直接取abs,计算误差会引起虚部不为零。No!  

                           g    = gp(1:M,1:N);        % 截取有效数据

取指数,得空域滤波结果。

                           g = exp(g)-1;                % 先前取对数时,为log(1+I)。 有代码先在频域做exp(Gp),再反变换。No!

显示g时,适当做数据的映射处理,以适合人眼观察。

以下为利用以上原理和关键代码实现的两个场景下的同态滤波实验。H参数的设置为:D0=80,  rL =0.25,  rH=2.2, c=2.0。其幅频曲线如图3所示。


图3 设置参数下的同态滤波器幅频图


图4 同态滤波结果,PET扫描图(上),隧洞口场景(下)

注:以上代码实现中,未考虑边界问题。关于边界处理问题,请参考前述博文:如何保持空域与频域滤波结果的一致性(续)

读者如需完整代码,可直接联系本人。

推荐参考:

[1]https://en.wikipedia.org/wiki/Homomorphic_filtering

[2]http://blogs.mathworks.com/steve/2013/06/25/homomorphic-filtering-part-1/

[3]http://blogs.mathworks.com/steve/2013/07/10/homomorphic-filtering-part-2/

[4]http://bbs.sjtu.edu.cn/bbstcon?board=DSP&reid=1316225259


扫一扫,可关注“老马迷图”微信公众号!

《教学后记》序列博文:

[1] 频域Laplacian图像锐化原理与实现

[2] 如何保持空域与频域滤波结果的一致性(续)

[3] 傅里叶变换的波形分辨率与频率分辨率

[4] 如何保持空域与频域滤波结果的一致性



https://wap.sciencenet.cn/blog-425437-1052070.html

上一篇:频域Laplacian图像锐化原理及实现
下一篇:最朴实的话语,却最温暖人心

10 徐令予 张云 信忠保 徐晓 李曙 魏焱明 haipengzhangdr styxmshy zouzoukankan yangjz2001

该博文允许注册用户评论 请点击登录 评论 (6 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2022-1-25 10:28

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部