数字图像处理-图像增强总结

图像增强技术：包括空域和频域两部分

空间域图像增强

• 直接对图像的像素本身进行操作
• 空域图像增强分为点处理和邻域处理。

点处理

• 在像素 (x, y) 处 g(x, y) 的值仅取决于 f(x, y) 的值
• 增强操作即为灰度级映射：s = E®
• 输出图像的像素值 g(x, y) 仅与输入图像中位于 (x, y) 的像素 f(x, y) 有关

邻域处理

• 像素 (x, y) 的邻域定义为中心位于像素 (x, y) 处的局部区域 Sxy
• 增强操作是以像素 (x, y) 为中心的邻域内所有像素值的函数：g(x, y) = E{f(s, t) | (s, t) ∈ Sxy}
• 输出图像的像素值 g(x, y) 与输入图像中位于 (x, y) 邻域内的像素都有关。

灰度级变换（点处理）

理论基础：灰度级变换将各像素灰度，都按同一函数进行变换。

技术要点：根据输入图像f(x,y)灰度值和输出图像g(x,y)灰度值之间的转换函数g(x，y)=T[f(x，y)]进行，包括线性变换和非线性变换，这种变换与像素的坐标无关，只和灰度级有关。

灰度变换可使图像动态范围增大，对比度得到扩展，使图像清晰、特征明显，是图像增强的重要手段之一。它主要利用图像的点运算来修正像素灰度，由输入像素点的灰度值确定相应输出像素点的灰度值，可以看作是“从像素到像素”的变换操作，不改变图像内的空间关系。

图像反转

理论基础：即反转图像灰度级，图像的每个像素点都进行翻转

技术要点：图像反转的公式如下，s为目标图像的像素点的像素值，r为原图像像素点的像素值

$$ s=L-1-r $$

实例：openCV实现

import cv2
import numpy as np

img = cv2.imread(r'C:\Users\xxx\Desktop\breast.tif',0)

reverse_img = 255 - img

cv2.imshow('srcimg',img)
cv2.imshow('reverse_img',reverse_img)
cv2.waitKey(0)

对数变换

理论基础：即对图像上的每一个像素点进行对数函数操作。可用于增强图像的暗部细节。

技术要点：对数变换的公式如下，s为目标图像的像素点的像素值，r为原图像像素点的像素值

$$ s=clog(1+r) $$

当c=1时，r与s关系图如下：

实例：

import cv2
import math
import numpy as np

def logTransform(c,img):

    #3通道RGB
    '''h,w,d = img.shape[0],img.shape[1],img.shape[2]
    new_img = np.zeros((h,w,d))
    for i in range(h):
        for j in range(w):
            for k in range(d):
                new_img[i,j,k] = c*(math.log(1.0+img[i,j,k]))'''

    #灰度图专属
    h,w = img.shape[0], img.shape[1]
    new_img = np.zeros((h, w))
    for i in range(h):
        for j in range(w):
            new_img[i, j] = c * (math.log(1.0 + img[i, j]))


    new_img = cv2.normalize(new_img,new_img,0,255,cv2.NORM_MINMAX)

    return new_img
    
#替换为你的图片路径
img = cv2.imread(r'C:\Users\xxx\Desktop\Fourier spectrum.tif',0)

log_img = logTransform(1.0,img)
cv2.imshow('log_img',log_img)
cv2.imwrite(r'C:\Users\xxx\Desktop\Fourier spectrum2.jpg',log_img)
cv2.waitKey(0)

下列图像分别是傅里叶频谱和应用上式中的对数变换(c=1)的结果：

由于对数曲线在像素值较低的区域斜率大，在像素值较高的区域斜率较小，所以图像经过对数变换后，较暗区域的对比度将有所提升。可用于增强图像的暗部细节。

幂律（伽马）变换

理论基础：即对图像上的每一个像素点进行伽马变换。伽马变换可以很好地拉伸图像的对比度，扩展灰度级。

技术要点：对数变换的公式如下，其中c、γ 为常数，s为目标图像的像素点的像素值，r为原图像像素点的像素值

$$ s=cr^y $$

对于不同的 γ 值，s 与 r的关系曲线如下图所示：

• 当图像的整体灰度偏暗时，选择γ<1,可以使图像增亮；
• 当图像的整体灰度偏亮时，选择γ>1,可以使图像变暗，

实例：

import math
import numpy as np
import cv2

def gammaTranform(c,gamma,image):
    h,w,d = image.shape[0],image.shape[1],image.shape[2]
    new_img = np.zeros((h,w,d),dtype=np.float32)
    for i in range(h):
        for j in range(w):
            new_img[i,j,0] = c*math.pow(image[i, j, 0], gamma)
            new_img[i,j,1] = c*math.pow(image[i, j, 1], gamma)
            new_img[i,j,2] = c*math.pow(image[i, j, 2], gamma)
    cv2.normalize(new_img,new_img,0,255,cv2.NORM_MINMAX)
    new_img = cv2.convertScaleAbs(new_img)

    return new_img

img = cv2.imread(r'C:\Users\xxx\Desktop\gray.jpg',1)

new_img = gammaTranform(1,2.5,img)

cv2.imshow('x',new_img)
cv2.imwrite(r'C:\Users\xxx\Desktop\gray_2.5.jpg',new_img)
cv2.waitKey(0)

γ<1 增强亮度

以及c=1时的γ=0.6、0.4、0.3时应用上公示的结果

分段线性变换函数

理论基础：为突出感兴趣的目标或灰度区间，相对抑制不感兴趣的灰度区域，可采用分段线性变换。

技术要点：根据分段操作不同可以分为：对比度拉伸、灰度级分层、比特平面封层等等。下面以对比度拉伸演示。

实例：

import cv2
import imutils
import numpy as np

# image = cv2.imread('E:/city.PNG')

image = cv2.imread('E:/city.PNG')
gray_img = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

# 在灰度图进行分段线性对比度拉伸
# 此种方式变换函数把灰度级由原来的线性拉伸到整个范围[0, 255]
r_min, r_max = 255, 0
for i in range(gray_img.shape[0]):
    for j in range(gray_img.shape[1]):
        if gray_img[i, j] > r_max:
            r_max = gray_img[i, j]
        if gray_img[i, j] < r_min:
            r_min = gray_img[i, j]
r1, s1 = r_min, 0
r2, s2 = r_max, 255

precewise_img = np.zeros((gray_img.shape[0], gray_img.shape[1]), dtype=np.uint8)
k1 = s1 / r1
k3 = (255 - s2) / (255 - r2)
k2 = (s2 - s1) / (r2 - r1)
for i in range(gray_img.shape[0]):
    for j in range(gray_img.shape[1]):
        if r1 <= gray_img[i, j] <= r2:
            precewise_img[i, j] = k2 * (gray_img[i, j] - r1)
        elif gray_img[i, j] < r1:
            precewise_img[i, j] = k1 * gray_img[i, j]
        elif gray_img[i, j] > r2:
            precewise_img[i, j] = k3 * (gray_img[i, j] - r2)

cv2.imshow('origin image', imutils.resize(image, 480))
cv2.imshow('gray image', imutils.resize(gray_img, 480))
cv2.imshow('precewise image', imutils.resize(precewise_img, 480))
if cv2.waitKey(0) == 27:
    cv2.destroyAllWindows()

灰度图像：

拉伸后：

算数运算

理论基础：直接对图像每个点的像素的灰度值进行算运算。

技术要点：图像相加减 = 灰度值相加减

相加运算：

s(x, y) = f(x, y) + g(x, y)

主要应用

• 对同一场景的多幅序列图像求取平均值，降低加性随机噪声；
• 将一幅图像叠加到另一幅图像上，实现二次曝光效果。

相减运算：

d(x, y) = f(x, y) - g(x, y)

主要应用

• 显示两幅图像的差异，检测同一场景两幅图像之间的变化，如：视频中镜头边界的检测，帧间变化检测
• 图像分割：如分割运动的车辆，减法去掉静止部分，剩余的是运动像素和噪声（背景减除法），数字减影血管造影

不做过多实例演示了吧

直方图处理

直方图均衡化

理论基础：对在图像中像素个数多的灰度值（即对画面起主要作用的灰度值）进行展宽，而对像素个数少的灰度值（即对画面不起主要作用的灰度值）进行归并，从而增大对比度，使图像清晰，达到增强的目的。

技术要点（计算过程）：

• 得到原始图片的灰度直方图
• 得到各个灰度级对应的概率密度函数
• 通过概率密度函数得到累积分布函数
• 累计分布函数乘以255，得到每一个灰度级对应的新的灰度
• 将旧灰度映射得到新的灰度，即更新整张图片的灰度

实例：

import cv2
import matplotlib.pyplot as plt

# 1.以灰度图的形式读入
img = cv2.imread(path_02, 0)

# 2.均衡化处理
dst = cv2.equalizeHist(img)

# 3.计算原图和均衡化后的直方图
cul_img = cv2.calcHist([img], [0], None, [256], [0, 256])
cul_dis = cv2.calcHist([dst], [0], None, [256], [0, 256])

# 3.结果展示
fig, axex = plt.subplots(nrows=2, ncols=2, figsize=[10, 8], dpi=100)
axex[0][0].imshow(img, cmap=plt.cm.gray)
axex[0][0].set_title("原图")
axex[0][1].imshow(dst, cmap=plt.cm.gray)
axex[0][1].set_title("均衡化的结果")
axex[1][0].plot(cul_img)
axex[1][0].grid()
axex[1][1].plot(cul_dis)
axex[1][1].grid()

直方图规定化

从上面可以看出，直方图的均衡化自动的确定了变换函数，可以很方便的得到变换后的图像，但是在有些应用中这种自动的增强并不是最好的方法。有时候，需要图像具有某一特定的直方图形状（也就是灰度分布），而不是均匀分布的直方图，这时候可以使用直方图规定化。

理论基础：直方图规定化，也叫做直方图匹配，用于将图像变换为某一特定的灰度分布，也就是其目的的灰度直方图是已知的。规定化操作能够有目的的增强某个灰度区间，相比于，均衡化操作，规定化多了一个输入，但是其变换后的结果也更灵活。

技术要点：

• 将原始图像的灰度直方图进行均衡化，得到一个变换函数s=T(r)s=T(r)，其中s是均衡化后的像素，r是原始像素
• 对规定的直方图进行均衡化，得到一个变换函数v=G(z)v=G(z)，其中v是均衡化后的像素，z是规定化的像素
• 上面都是对同一图像的均衡化，其结果应该是相等的,s=v,且z=G−1(v)=G−1(T(r))s=v,且z=G−1(v)=G−1(T(r))

实例：

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

"""
直方图匹配： 让一张图参考另一张图， 让他们的色调保持一致
    步骤：
        计算原图累计直方图
        计算参考图的累计直方图
        计算两个累计直方图的差异
        生成原图和参考图之间的颜色映射
"""
# 计算单通道图像 累计概率
def getAccumulateRatios(img):

    height = img.shape[0]
    width = img.shape[1]

    # 1.计算直方图
    hist = cv.calcHist([img],[0],None,[256],[0,255])
    # print(hist)
    # plt.plot(hist)
    # plt.show()

    # 2.计算每个灰度值出现的概率
    ratios = hist/(height*width);

    # 3. 计算累计概率
    sumRatios = np.zeros(256,np.float)
    sum1=0;
    for i,r in enumerate(ratios):
        sum1 = sum1 + r;
        sumRatios[i] = sum1;

    # 4. 绘制累计概率直方图
    # x = np.linspace(0,255,256);
    # plt.bar(x,sumRatios)
    # plt.show()

    return sumRatios;


"""
 传入进来的数据仍然上单通道数据
    1. 计算累计直方图的差异
        找到最小的差异
    2. 获取颜色映射    
"""
def map_color(src_channel,refer_channel):
    # src 单通道累计直方图
    src_sumRatios = getAccumulateRatios(src_channel);
    # refer 单通道累计直方图
    refer_sumRatios = getAccumulateRatios(refer_channel);

    colorMaps = np.zeros(256,np.uint8);

    # 遍历原图每一个灰度的累计概率
    for i,srcRatio in enumerate(src_sumRatios):
        # 得到 i: 灰度值  ratio ：累计概率  0  0.19
        min = 100;
        referColor = 0;
        for j,referRatio in enumerate(refer_sumRatios):
            # j: 表示参考图的灰度值， referRatio累计概率
            diff = np.abs(srcRatio - referRatio)
            if diff < min:
                # 更新最小值
                min = diff;
                referColor = j;

        # 0 ---> referColor
        colorMaps[i] = referColor;

    # print(colorMaps)
    return colorMaps;


""" 完成单通道直方图匹配  """
def oneChannelMatch(src_channel,refer_channel):
    # 单通道颜色值映射
    colorMaps = map_color(src_channel, refer_channel);

    # 绘制原图直方图
    # cv.imshow("src",src_channel)
    # plt.hist(src_channel.ravel(),bins=256,color="blue")
    # plt.show()
    # 绘制参考图直方图
    # cv.imshow("refer", refer_channel)
    # plt.hist(refer_channel.ravel(), bins=256,color="green")
    # plt.show()

    one_channel = src_channel
    height = one_channel.shape[0]
    width = one_channel.shape[1]
    for row in range(height):
        for col in range(width):
            # 获取原来的灰度值
            gray = one_channel[row, col];
            # 去颜色映射表中查找新的颜色值
            referColor = colorMaps[gray];
            # 替换为新的颜色值
            one_channel[row, col] = referColor;

    # 绘制生成之后的直方图
    # cv.imshow("dst", one_channel)
    # plt.hist(refer_channel.ravel(), bins=256, color="red")
    # plt.show()

    return one_channel

if __name__ == '__main__':
    src = cv.imread("../img/1.jpg",cv.IMREAD_COLOR);
    cv.imshow("src",src)
    refer = cv.imread("../img/2.jpg",cv.IMREAD_COLOR);
    cv.imshow("refer",refer)

    # 先计算src的累计直方图 , 计算单通道
    src_channels = cv.split(src);

    # 先计算refer的累计直方图 , 计算单通道
    refer_channels = cv.split(refer);
    
    #dst_channel0 = oneChannelMatch(src_channels[0],refer_channels[0])
    #dst_channel1 = oneChannelMatch(src_channels[1],refer_channels[1])
    #dst_channel2 = oneChannelMatch(src_channels[2],refer_channels[2])
    results = []
    fro i in range(3):
        res = oneChannelMatch(src_channels[i],refer_channels[i])
        reslut.append(res)
    dst = cv.merge(results)
    cv.imshow("dst",dst);


    cv.waitKey(0)
    cv.destroyAllWindows()

原图：

规定话过后：

空间域滤波（领域处理）

理论基础：空间滤波，就是直接在灰度值上，做一些滤波操作。滤波一词，其实来源于频域，将某个频率成分滤除的意思。大部分线性的空间滤波器（比如均值滤波器），是在空间上进行一些灰度值上的操作，达到平滑或锐化图像的作用。包括平滑滤波器和锐化滤波器。

技术要点：空间滤波工作机理如上图。模板或者说核，掩膜逐个点移动，如果是线性滤波的话，那么（x,y）这个点经过滤波后的像素值由原来的a,变为模板系数和该系数对应的像素值的乘积和。

平滑滤波器

理论基础：平滑滤波器用于模糊处理和降低噪声. 常应用于预处理任务中, 例如在大目标提取之前去除图像中的一些琐碎细节, 以及桥接直线或曲线的缝隙, 通过线性滤波和非线性滤波模糊处理, 可以降低噪声。

技术要点：平滑线性滤波器的输出是包含滤波器模板邻域内的像素的简单平均值, 这些滤波器也被成为均值滤波器。

根据卷积核的不同也可以分为不同的滤波器：盒式滤波器、高斯核滤波器、统计排序滤波器

盒式滤波器的滤波核

高斯滤波器的核

统计排序滤波器：

基于滤波器邻域中的像素值的顺序，排序结果决定了滤波器的输出。

实例：统计排序滤波器-中值滤波

# 9.12: 统计排序滤波器 (Statistical sorting filter)
    img = cv2.imread("../images/Fig0508a.tif", 0)  # flags=0 读取为灰度图像
    img_h = img.shape[0]
    img_w = img.shape[1]

    m, n = 3, 3
    kernalMean = np.ones((m, n), np.float32)  # 生成盒式核

    # 边缘填充
    hPad = int((m-1) / 2)
    wPad = int((n-1) / 2)
    imgPad = np.pad(img.copy(), ((hPad, m-hPad-1), (wPad, n-wPad-1)), mode="edge")

    imgMedianFilter = np.zeros(img.shape)  # 中值滤波器
    for i in range(img_h):
        for j in range(img_w):
            # # 1. 中值滤波器 (median filter)
            pad = imgPad[i:i+m, j:j+n]
            imgMedianFilter[i, j] = np.median(pad)

    plt.figure(figsize=(9, 7))
    plt.subplot(221), plt.axis('off'), plt.title("median filter")
    plt.imshow(imgMedianFilter, cmap='gray', vmin=0, vmax=255)
    plt.tight_layout()
    plt.show()

锐化滤波器

理论基础： 锐化滤波器在去噪的同时凸显物体的边缘信息，保持图像边缘信息不变或者是增强边缘信息。图像锐化的实质上增强原图像的高频分量，图像锐化滤波器为高通滤波器，边缘和轮廓一般位于灰度突变的地方，因此可以使用灰度差分提取图像边缘和轮廓。

技术要点：以锐化算子区分有二阶微分的图像锐化：拉普拉斯锐化、一阶微分的图像锐化：梯度锐化两种

二阶微分的图像锐化：拉普拉斯锐化

理论基础： 用算子产生的图像暗色背景叠加浅灰色边及突变点明显。由于拉普拉斯是一种微分算子，拉普拉斯图像强调原图中的灰度突变区域，衰减灰度变化慢区域，恒定区域变为0。

技术要点：将原始图像和拉普拉斯图像叠加在一起的简单方法可以保护拉普拉斯锐化处理的效果，同时又能复原背景信息。所以适用拉普拉斯变换对图像锐化增强的基本方法可表示为：

实例：

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

def Laplace(img, kernel):

    des_8U = cv2.filter2D(img, -1, kernel=kernel, borderType=cv2.BORDER_DEFAULT)
    des_16S = cv2.filter2D(img, ddepth=cv2.CV_16SC1, kernel=kernel, borderType=cv2.BORDER_DEFAULT)

    g = img - des_16S
    g[g<0] = 0
    g[g>255] = 255

    plt.figure(figsize=(10,14))

    # origin, des_8U, des_16S, filtered
    plt.subplot(221)
    plt.imshow(img, cmap='gray')
    plt.title('origin')

    plt.subplot(222)
    plt.imshow(des_8U, cmap='gray')
    plt.title('des-8U')

    plt.subplot(223)
    plt.imshow(des_16S, cmap='gray')
    plt.title('des-16S')

    plt.subplot(224)
    plt.imshow(g, cmap='gray')
    plt.title('g')
    plt.show()

img0 = 'moon.tif'

f = cv2.imread(img0, cv2.IMREAD_GRAYSCALE)

kernel1 = np.asarray([[0, 1, 0],
                      [1, -4, 1],
                      [0, 1, 0]])

Laplace(f, kernel1)

kernel2 = np.asarray([[1, 1, 1],
                      [1, -8, 1],
                      [1, 1, 1]])

Laplace(f, kernel2)

一阶微分的图像锐化：梯度锐化

理论基础：梯度向量的分量是微分，是线性算子，但梯度的幅值不是线性算子，是做了平方和方根。梯度的偏微分不是旋转不变的，梯度向量的幅值是旋转不变的。有时计算用绝对值来近似平方和方根。

技术要点：Sobel算子：

作用：

• 去除变化慢的背景
• 在灰度平坦区域中增强细小的突变
• 突出灰度图像中看不见的斑点

实例：

import cv2
import numpy as np
import copy

OriginalImg = cv2.imread('cns.jpg')
R ,G, B = cv2.split(OriginalImg)
#RGB图像拆分成R，G，B三通道数据


G_pad = np.pad(G,((1, 1), (1, 1)),'edge')
#G通道到数据填充
newImg = np.full((700,700),np.nan)
#创建空np数组，用于存储sobel算子处理后的数据
for i in range(1,700):
    for j in range(1,700):
        tmp = copy.copy(G_pad[i-1:i+2,j-1:j+2])
        x = np.multiply(tmp, np.array([[-1,2,-1],[0,0,0],[1,2,1]]))
        y = np.multiply(tmp, np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]]))
        dF = abs(np.sum(x)) + abs(np.sum(y))
        newImg[i-1,j-1] = dF

newImg = np.uint8(newImg)
#图像边缘填充

cv2.namedWindow('G')
cv2.resizeWindow('G',(280,280))
cv2.imshow('G',G)

cv2.namedWindow('newImg')
cv2.resizeWindow('newImg',(280,280))
cv2.imshow('newImg',newImg)
cv2.imwrite('sobelImg.jpg',newImg)
cv2.waitKey(0)

梯度锐化后：

原图：

频率域图像增强

从空间域变换到频率域，傅里叶变换可以做到转换过程不丢失任何信息。

一维傅里叶变换：

单变量连续函数f(x)的傅里叶变换F(u)定义为等式：

反变换可以定义为：

上述这两个式子，就反映了通过一个函数可以做变换求到其傅里叶变换，或已知傅里叶变换可以完全地求出原始函数。单变量离散函数f(x)（其中x=0,1,2,…,m-1）的傅里叶变换由以下等式给出：

同样，若给出F(u)，能用反DFT来获得：

F(u)的值的范围覆盖的域（u的值）称为频率域，因为u决定了变换的频率成分。F(u)的M项中的每一个被称为变换的频率分量。

二维傅里叶变换：

二维连续函数f(x, y)的傅里叶变换F(u, v)定义为：

由此可以得它的反变换：

对于图像尺寸为M*N的函数f(x, y)的二维离散傅里叶变换为：

给出F(u, v)，可通过反DFT得到f(x, y）如下：

注：u和v是频率变量，x和y是空间或图像变量。

频率与滤波

傅里叶变换的结果，与图像中的强度变化模式具有一定的联系。 例如，变化最慢的频率成分(u=v=0)对应一幅图像的平均灰度级。

反映的是对图像每个像素点的灰度级求和然后除以MN，得到平均灰度级。低频对应着图像的总体灰度级的显示，高频对应图像中变化快的分量（图像的细节）。

步骤：

• 用(-1)^(x+y)乘以输入图像来进行中心变换
• 由(1)计算图像的DFT，得到F(u，v)
• 用滤波器函数H(u, v)乘以F(u, v)
• 计算(3)中结果的反DFT
• 得到(4)中结果的实部
• 用(-1)^(x+y)乘以(5)中的结果

低通滤波器

理论基础：“截断”傅里叶变换中的所有高频成分，这些成分处在距变换原点的距离比指定距离D0要远的多的位置。

技术要点：其变换函数为：

其中，D0是指定的非负数值，D(u,v)是(u,v)点距频率矩形原点的距离。由于变换被中心化了，如果要研究的图像尺寸为M*N，从点(u,v)到傅里叶变换中心的距离为：

滤波器h(x,y)有两个主要特性：在原点处的一个主要成分，及中心成分周围呈周期性的成分。中心成分主要决定模糊，周期性的成分主要决定了理想滤波器振铃现象的特性。 中心成分的半径和距原点每单位距离上周期的数量都与理想滤波器的截止频率成反比。

实例：

>> Fimg=fft2(double(f));
>> Fimg=fftshift(Fimg);
>> [M,N]=size(f);
>> dist1=5; 
>> z1=zeros(M,N);  
>> for i=1:M 
    for j=i:N 
       if(sqrt(((i-M/2)^2+(j-N/2)^2))<dist1) 
           z1(i,j)=1; 
        end 
    end 
end  
>> g1=Fimg.*z1;  
>> g1=ifftshift(g1); 
>> img1=real(ifft2(g1)); 
>> imshow(img1);

高通滤波器

理论基础：让高频率通过而滤掉或衰减低频,其作用是使图像得到锐化处理,突出图像的边界。

技术要点：其变换函数为：

其中D0为理想高通滤波器的截止频率，注意与理想高通滤波器的表达式区别在于将0和1交换而条件不变。

实例：matlab

close all;
clear all;
clc;
I = imread('coins.png');

subplot(121),imshow(I);
title('原始图像');

% 函数fft2()用于计算二维傅立叶变换
% 函数fftshift()是对函数fft2()作傅里叶变换后得到的频谱进行平移,将变换后的图像频谱中心从矩阵的原点移到矩阵的中心
% 作二维傅里叶变换前一定要用函数im2double()把原始图像的数据类型由uint8转化为double类型
% 否则会因为unit8数据类型只能表示0-255的整数而出现数据截断,进而出现错误结果
s=fftshift(fft2(im2double(I)));
[a,b]=size(s);
a0=round(a/2);
b0=round(b/2);
d0=50; % 将理想高通滤波器的截止频率D0设置为50
for i=1:a %双重for循环计算频率点(i,j)与频域中心的距离D(i,j)=sqrt((i-round(a/2)^2+(j-round(b/2)^2))
    for j=1:b 
        distance=sqrt((i-a0)^2+(j-b0)^2);
        if distance<=d0  % 根据理想高通滤波器产生公式,当D(i,j)<=D0,置为0
            h=0;
        else
            h=1;        % 根据理想高通滤波器产生公式,当D(i,j)>D0,置为1
        end
        s(i,j)=h*s(i,j);% 频域图像乘以滤波器的系数
    end
end
% real函数取元素的实部
s=real(ifft2(ifftshift(s)));% 最后进行二维傅里叶反变换转换为时域图像
subplot(122),imshow(s,[]);
title('理想高通滤波所得图像');

同态滤波器

理论基础：一种线性滤波在不同域中的非线性映射。

技术要点：一幅图像可看成由两部分组成，即

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

同态滤波过程，分为以下5个基本步骤：
① 原图做对数变换，得到如下两个加性分量，即

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

③ 设计一个频域滤波器H(u,v)，进行对数图像的频域滤波。
④ 傅里叶反变换，返回空域对数图像。
⑤ 取指数，得空域滤波结果。

实例：

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()