Python实现频域同态滤波增强医学影像灰度图像

由于转移函数H是二维的，所以需要对图像进行二维傅里叶变换。在进行傅里叶变换之前，需要将灰度图像进行归一化处理，使像素值范围在[0,1]之间。

代码实现如下：

import cv2
import numpy as np
import matplotlib.pyplot as plt

def homomorphic_filter(img, r_h, r_l, c, D_0, D_1):
    '''
    :param img: 原始灰度图像
    :param r_h: 高频增益
    :param r_l: 低频增益
    :param c: 控制滤波器的带宽
    :param D_0: 高斯低通滤波器截止频率
    :param D_1: 高斯高通滤波器截止频率
    :return: 增强后的灰度图像
    '''
    rows, cols = img.shape

    # 归一化处理
    img = img / 255

    # 中心化图像
    img = img * (-1) ** (np.arange(rows) + np.arange(cols)[:, np.newaxis])

    # 二维傅里叶变换
    img_fft = np.fft.fft2(img)

    # 高斯低通滤波器
    H_lp = np.zeros((rows, cols))
    for i in range(rows):
        for j in range(cols):
            H_lp[i, j] = np.exp(-((i - rows / 2) ** 2 + (j - cols / 2) ** 2) / (2 * D_0 ** 2))

    # 高斯高通滤波器
    H_hp = np.zeros((rows, cols))
    for i in range(rows):
        for j in range(cols):
            H_hp[i, j] = 1 - np.exp(-c * ((i - rows / 2) ** 2 + (j - cols / 2) ** 2) / (D_1 ** 2))

    # 同态滤波器
    H = r_h * H_hp + r_l * H_lp

    # 滤波
    img_fft_filtered = H * img_fft

    # 反变换
    img_filtered = np.real(np.fft.ifft2(img_fft_filtered))

    # 反中心化
    img_filtered = img_filtered * (-1) ** (np.arange(rows) + np.arange(cols)[:, np.newaxis])

    # 反归一化
    img_filtered = img_filtered * 255

    return img_filtered

# 读取灰度图像
img = cv2.imread('medical_image.jpg', cv2.IMREAD_GRAYSCALE)

# 同态滤波增强
img_filtered = homomorphic_filter(img, r_h=2, r_l=0.5, c=1, D_0=50, D_1=10)

# 显示增强后的图像
plt.imshow(img_filtered, cmap='gray')
plt.axis('off')
plt.show()