图像处理——同态滤波

1. 什么是同态滤波
   同态滤波（Homomorphic filter）是信号与图像处理中的一种常用技术，它采用了一种线性滤波在不同域中的非线性映射。

2. 首先介绍两个概念
   同态系统：是将非线性问题，转化为线性问题处理。即对非线性（乘性）混杂信号，通过某种数学运算（如对数变换），变成加性模型，而后采用线性滤波方法进行处理。
   同态滤波：是把频率滤波和空域灰度变换结合起来的一种图像处理方法，它根据图像的照度/反射率模型作为频域处理的基础，利用压缩亮度范围和增强对比度来改善图像的质量。

3. 同态滤波的原理
   一幅图像可看成由两部分组成，即
   fi代表随空间位置不同的光强（Illumination）分量，其特点是缓慢变化，集中在图像的低频部分。
   fr代表景物反射到人眼的反射（Reflectance）分量，其特点包含了景物各种信息，高频成分丰富。
   同态滤波过程，分为以下5个基本步骤：
   ① 原图做对数变换，得到如下两个加性分量，即
   ② 对数图像做傅里叶变换，得到其对应的频域表示为：
   ③ 设计一个频域滤波器H(u,v)，进行对数图像的频域滤波。
   ④ 傅里叶反变换，返回空域对数图像。
   ⑤ 取指数，得空域滤波结果。
   综上，同态滤波的基本步骤如下图所示。
   
   首先上面的分析一共有三点需要注意：
   ① 取对数，目的是基于我们假定的图像模型将低频和高频信号进行分离
   ② 然后频域处理，即减弱低频增强高频达到增强细节的作用
   ③ 反操作得到处理后图像
   可以看出，同态滤波的关键在于滤波器H的设计。对于一幅光照不均匀的图像，同态滤波可同时实现亮度调整和对比度提升，从而改善图像质量。为了压制低频的亮度分量，增强高频的反射分量，滤波器H应是一个高通滤波器，但又不能完全cut off 低频分量，仅作适当压制。
   因此，同态滤波器一般采用如下形式，即
   
   其中，rL< 1, rH >1，控制滤波器幅度的范围。Hhp通常为高通滤波器，如高斯（Gaussian）高通滤波器、巴特沃兹(Butterworth)高通滤波器、Laplacian滤波器等。
   如果Hhp采用Gaussian高通滤波器，则有：
   其中，D0为gaussian滤波器的截止频域，c为一个常数，控制滤波器的形态，即从低频到高频过渡段的陡度（斜率），其值越大，斜坡带越陡峭，见图2。
   其中我们需要关注的参数一共有4个，分别为：
   ① rh 高频权重
   ②rl 低频权重
   ③ c 高斯函数的陡峭程度(并不是很重要)
   ④d0 手动区分低频和高频的比例
   图2 同态滤波器幅频曲线

4. 同态滤波作用
   同态滤波利用去除乘性噪声(multiplicative noise)，可以同时增加对比度以及标准化亮度，借此达到图像增强的目的。
   一副图像可以表示为其照度(illumination)分量和反射(reflectance)分量的乘积，虽然在时域上这两者是不可分离的，但是经由傅立叶转换两者在频域中可以线性分离。由于照度可视为环境中的照明，相对变化很小，可以看作是图像的低频成分；而反射率相对变化较大，则可视为高频成分。通过分别处理照度和反射率对像元灰度值的影响，通常是借由高通滤波器(high-pass filter)，让图像的照明更加均匀，达到增强阴影区细节特征的目的。

5. 代码实现

import cv2
import numpy as np


def homomorphic_filter(src, d0=10, rl=0.5, rh=2.0, c=4, h=2.0, l=0.5):
    gray = src.copy()
    if len(src.shape) > 2:
        gray = cv2.cvtColor(src, cv2.COLOR_BGR2GRAY)  # 转换成灰度图像

    gray = np.log(1e-5 + gray)  # 取对数

    rows, cols = gray.shape
    gray_fft = np.fft.fft2(gray)  # FFT傅里叶变换
    gray_fftshift = np.fft.fftshift(gray_fft)  # FFT中心化

    M, N = np.meshgrid(np.arange(-cols // 2, cols // 2), np.arange(-rows // 2, rows // 2))
    D = np.sqrt(M ** 2 + N ** 2)  # 计算距离
    Z = (rh - rl) * (1 - np.exp(-c * (D ** 2 / d0 ** 2))) + rl  # H(u,v)传输函数

    dst_fftshift = Z * gray_fftshift
    dst_fftshift = (h - l) * dst_fftshift + l
    dst_ifftshift = np.fft.ifftshift(dst_fftshift)

    dst_ifft = np.fft.ifft2(dst_ifftshift)  # IFFT逆傅里叶变换
    dst = np.real(dst_ifft)  # IFFT取实部

    dst = np.exp(dst) - 1  # 还原
    dst = np.uint8(np.clip(dst, 0, 255))
    return dst


img = cv2.imread('hf1.jpg')
img_new = homomorphic_filter(img)

cv2.imshow("img", img)
cv2.imshow("img_new", img_new)
cv2.imwrite("img_new1.jpg", img_new)

key = cv2.waitKey(0)
cv2.destroyAllWindows()

参考链接：数字图像处理杂项-同态滤波