同态滤波

Python图像的频域图像增强-同态滤波器黑白图片处理黑边处理2有注释代码
数字图像空间域 频域
同态滤波是一种在频域中同时能够压缩图像的亮度范围和增强图像对比度的方法。下面推导同态滤波的流程:
基于图像成像模型:一幅图像f(x, y)可以表示成它的照度分量i(x, y)与反射分量r(x, y)的乘积。
成像模型(照度和反射) :
f(x,y)=i(x,y)r(x,y)
(1)两边取对数:
lnf(x,y)= lni(x, y) +Inr(x,y)
(2)两边取傅立叶变换:
F(u,v)= I(u,v) + R(u,v)
(3)用一频域函数H(u, v)处理F(u, v):
H(u,v)F(u,v)= H(u,v)I(u,v)+ H(u,v)R(u,v)
算法解释
对于一幅由物理过程产生的图像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)。取完对数后先进行归一化再进行傅里叶变换。在这个过程中，由于f(x,y)的取值范围为[0, L-1]，为了避免出现ln(0)的情况，故采用ln ( f(x,y) + 1 ) 来表示。
然后取二维傅里叶变换，得到 Z(u,v) = Fi(u,v) + Fr(u,v)，并使用np.fft.fftshift()把转换后的频域图的低频部分放到图像中间，起到便于观察的作用。
再使用一个滤波器，对Z(u,v)进行滤波，有 S(u,v) = H(u,v) Z(u,v) = H(u,v)Fi(u,v) + H(u,v)Fr(u,v)。
滤波后，进行依次进行np.fft.ifftshift()和二维反傅里叶变换，有 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。由于之前使用了归一化这里在取完指数后再乘以255进行反归一化。 i0(x,y)和r0(x,y)分别是处理后图像的照射分量和入射分量。
黑白图片的同态滤波代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
import cv2
import numpy as np


def homomorphic_filter(src, d0=2, r1=1.0, rh=2.0, c=4, h=2.0, l=0.5):

    # 图像灰度化处理，BGR转为GRAY
    print(src.shape)
    gray = src.copy()

    if len(src.shape) > 2:  # 维度>2
        gray = cv2.cvtColor(src, cv2.COLOR_BGR2GRAY)

    # 图像格式处理
    # 转为np.float64格式
    gray = np.float64(gray)
    print(gray.dtype)
    print(gray.max(), gray.min())

    # 对数域
    gray = np.log(gray + 1.0)
    gray = gray / np.log(256)  # 归一化
    print(gray.max(), gray.min())

    # 傅里叶变换
    #np.fft.fft2() 二维傅里叶变换
    gray_fft = np.fft.fft2(gray)
    #np.fft.fftshift() 把转换后的频域图的低频部分放到图像中间，仅起到便于观察的作用
    gray_fftshift = np.fft.fftshift(gray_fft)

    # arange函数用于创建等差数组
    rows, cols = gray.shape
    print(rows,cols)
    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 - r1) * (1 - np.exp(-c * (D ** 2 / d0 ** 2))) + r1  # filter
    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)

    # 选取元素的模
    dst = np.abs(dst_ifft)
    print(dst.min(), dst.max())

    # 对数反变换
    dst = np.exp(dst) - 1
    print(dst.min(), dst.max())
    dst = (dst - dst.min()) / (dst.max() - dst.min())  # 归一化
    print(dst.min(), dst.max())
    dst *= 255  #反归一化？
    print(dst.min(), dst.max())

    # dst中，比0小的都会变成0，比255大的都变成255
    # uint8是专门用于存储各种图像的（包括RGB，灰度图像等），范围是从0–255
    """
    numpy中clip函数用法详解
    numpy.clip(a, a_min, a_max, out=None)
    参数说明
    a : 输入的数组
    a_min: 限定的最小值 也可以是数组 如果为数组时 shape必须和a一样
    a_max:限定的最大值 也可以是数组 shape和a一样
    out：剪裁后的数组存入的数组
    """
    dst = np.uint8(np.clip(dst, 0,255))
    return dst


if __name__ == "__main__":
    # 把图片变成了(400, 400)丢失了RGB通道信息
    img = cv2.resize(cv2.imread('./img/ys.jpg', 0), (400, 400))
    img_new = homomorphic_filter(img)
    """
    按水平方向（列顺序）堆叠数组构成一个新的数组
    堆叠的数组需要具有相同的维度
    """
    contrast = np.hstack((img, img_new))
    cv2.imshow('contrast', contrast)
    cv2.waitKey()
    cv2.destroyAllWindows()

彩色图片的同态滤波代码：python版本暂未发现

1
#python版本暂未发现