图像处理算法100问python实现-part3

32-35 傅里叶变换

32.傅里叶变换（Fourier Transform)

傅里叶变换通常用来分离模拟信号或音频信号等连续一维信号。但是，数字图像表示为[0,255]内的离散值，且用H*W的二维矩阵来表示，所以使用二维的离散傅里叶变换处理图像。
二维离散傅里叶变换公式如下：I为输入图像，G为频谱图。

二维离散傅里叶逆变换从频率分量按照下式复原图像：其中exp(j)为复数，编程时需用绝对值形态。

注意：用for语句非常耗时，善用numpy来减少计算量

33.傅里叶变换-低通滤波

保留傅里叶变换后的频谱图中的低频成分
图像中的低频分量，指的是图像强度（亮度/灰度）变换平缓的地方，也就是大片色块的地方。 人眼对图像中的高频信号更为敏感
通过离散傅立叶变换得到的频率在左上、右上、左下、右下等地方频率较低，在中心位置频率较高。需要将图像划分为四个区域进行对调，才可转换为低频在中间区，高频在边缘区的标准形式，如下图：构建低通滤波器时低频中心在图像中心

34.傅里叶变换-高通滤波

图像中的高频分量，指的是图像强度（亮度/灰度）变化剧烈的地方，也就是我们常说的边缘（轮廓）
图像中的低频分量，指的是图像强度（亮度/灰度）变换平缓的地方，也就是大片色块的地方。 人眼对图像中的高频信号更为敏感。

35.傅里叶变换-带通滤波

可以保留介于低频成分和高频成分之间的分量的带通滤波器
其他置为0.

总结:跑起来太耗时，不如用opencv自带的傅里叶变换函数

代码

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

def bgr2gray(img):
	gray = 0.2126 * img[..., 2] + 0.7152 * img[..., 1] + 0.0722 * img[..., 0]
	return gray
	
#傅里叶变换
def dft(img):
    H, W, C = img.shape
    K, L = H, W

    G = np.zeros((K, L, C), dtype=np.complex)

    # 利用numpy生成各像素的坐标矩阵
    # 生成行方向映射矩阵，并沿列方向复制
    new_w = np.tile(np.arrange(W), (H, 1))
    # 生成列方向映射矩阵，并沿行方向复制
    new_h = np.tile(np.arrange(H).reshape(H, -1), (1, W))

    for l in range(L):
        for k in range(K):
            for c in range(C):
                G[l, k, c] = np.sum(img[..., c] * exp(-2j * np.pi * (new_w * k / K + new_h * l / L))) / np.sqrt(K * L)
    return G

#使用numpy循环做傅里叶变换时间复杂度高，可使用opencv中的dft函数进行傅里叶变换
dft = cv2.dft(np.float32(img),flags = cv2.DFT_COMPLEX_OUTPUT)#傅里叶变换
dft_shift = np.fft.fftshift(dft)#低频平移到中心

# 傅里叶逆变换
def idft(G):
    # prepare out image
    H, W, _ = G.shape
    out = np.zeros((H, W, channel), dtype=np.float32)

    # prepare processed index corresponding to original image positions
    x = np.tile(np.arange(W), (H, 1))
    y = np.arange(H).repeat(W).reshape(H, -1)

    # idft
    for c in range(channel):
        for l in range(H):
            for k in range(W):
            #exp(j)是个复数，实际编程的时候请务必使用下式中的绝对值形态
                out[l, k, c] = np.abs(np.sum(G[..., c] * np.exp(2j * np.pi * (x * k / W + y * l / H)))) / np.sqrt(W * H)

    out = np.clip(out, 0, 255)
    out = out.astype(np.uint8)

    return out

#低通滤波器
#设从低频的中心到高频的距离为 ，我们保留0.5r的低频分量。
def lpf(G,ratio=0.5): 
    H,W,C=G.shape

    #四个区对调
    _G=np.zeros_like(G)
    _G[:H//2,:W//2]=G[H//2:,W/2:]
    _G[:H // 2, W // 2:] = G[H // 2:, 0:W / 2]
    _G[H // 2:, :W // 2] = G[:H // 2, w / 2:]
    _G[H // 2:, W // 2:] = G[:H // 2, :W // 2]
    # 生成各像素的坐标矩阵
    w_new = np.tile(np.arange(W), (H, 1))
    h_new = np.tile(np.arange(H).reshape(H, -1), (1, W))
    # 生成低通滤波器
    mask = np.ones((H, W), dtype=np.float32)
    r1 = w_new - W // 2
    r2 = h_new - H // 2
    r = np.sqrt(r1 ** 2 + r2 ** 2)
    mask[r >= (ratio * W//2)] = 0
    # 扩展为三维
    mask = np.repeat(mask, C).reshape(H, W, C)
    result *= mask
    mid = result.copy()
    # 1，4区互换
    result[:H // 2, :W // 2] = mid[H // 2:, W // 2:]
    result[H // 2:, W // 2:] = mid[:H // 2, :W // 2]
    # 2，3区互换
    result[H // 2:, :W // 2] = mid[:H // 2, W // 2:]
    result[:H // 2, W // 2:] = mid[H // 2:, :W // 2]
    return result

##高通
def HPF(image, range=0.2):
    H, W, C = image.shape
    # 傅里叶变换后矩阵含有虚部，因此result矩阵类型也应含有虚部
    result = np.zeros((H, W, C), dtype=np.complex)
    # 1，4区互换
    result[:H // 2, :W // 2] = image[H // 2:, W // 2:]
    result[H // 2:, W // 2:] = image[:H // 2, :W // 2]
    # 2，3区互换
    result[H // 2:, :W // 2] = image[:H // 2, W // 2:]
    result[:H // 2, W // 2:] = image[H // 2:, :W // 2]
    # 生成各像素的坐标矩阵
    w_new = np.tile(np.arange(W), (H, 1))
    h_new = np.tile(np.arange(H).reshape(H, -1), (1, W))
    # 生成高通滤波器
    mask = np.ones((H, W), dtype=np.float32)
    r1 = w_new - W//2
    r2 = h_new - H//2
    r = np.sqrt(r1 ** 2 + r2 ** 2)
    mask[r <= (range * W // 2)] = 0
    #扩展为三维
    mask = np.repeat(mask, C).reshape(H, W, C)
    result *= mask
    mid = result.copy()
    # 1，4区互换
    result[:H // 2, :W // 2] = mid[H // 2:, W // 2:]
    result[H // 2:, W // 2:] = mid[:H // 2, :W // 2]
    # 2，3区互换
    result[H // 2:, :W // 2] = mid[:H // 2, W // 2:]
    result[:H // 2, W // 2:] = mid[H // 2:, :W // 2]
    return result

##带通滤波
def MPF(image, range_bottle=0.1, range_top=0.5): #保留0.1r-0.5r间的频率
    H, W, C = image.shape
    result = np.zeros((H, W, C), dtype=np.complex)
    result[:H // 2, :W // 2] = image[H // 2:, W // 2:]
    result[H // 2:, W // 2:] = image[:H // 2, :W // 2]
    result[H // 2:, :W // 2] = image[:H // 2, W // 2:]
    result[:H // 2, W // 2:] = image[H // 2:, :W // 2]
    w_new = np.tile(np.arange(W), (H, 1))
    h_new = np.tile(np.arange(H).reshape(H, -1), (1, W))
    mask = np.ones((H, W), dtype=np.float32)
    r1 = w_new - W//2
    r2 = h_new - H//2
    r = np.sqrt(r1 ** 2 + r2 ** 2)
    mask[r >= (range_top * W // 2)] = 0
    mask[r <= (range_bottle * W // 2)] = 0
    mask = np.repeat(mask, C).reshape(H, W, C)
    result *= mask
    mid = result.copy()
    result[:H // 2, :W // 2] = mid[H // 2:, W // 2:]
    result[H // 2:, W // 2:] = mid[:H // 2, :W // 2]
    result[H // 2:, :W // 2] = mid[:H // 2, W // 2:]
    result[:H // 2, W // 2:] = mid[H // 2:, :W // 2]
    return r
 
# Read image
img = cv2.imread("imori.jpg").astype(np.float32)
H, W, C = img.shape

# Gray scale
gray = bgr2gray(img)

# DFT
G = dft(img)

# HPF
G = hpf(G)

# IDFT
out = idft(G)

# Save result
cv2.imshow("result", out)
cv2.waitKey(0)
cv2.imwrite("out.jpg", out)

#使用numpy循环做傅里叶变换时间复杂度高，可使用opencv中的dft函数进行傅里叶变换
dft = cv2.dft(np.float32(img),flags = cv2.DFT_COMPLEX_OUTPUT)#傅里叶变换
dft_shift = np.fft.fftshift(dft)#低频平移到中心

36-40 JPEG压缩算法