第5章 Python 数字图像处理(DIP) - 图像复原与重建14 - 逆滤波

逆滤波

逆滤波

图像的退化函数已知或者由前面的方法获取退化函数，则可以直接逆滤波

$$\hat{F}(u,v)=\frac{G(u,v)}{H(u,v)}\text{\hspace{1em}(5.78)}$$

$$\hat{F}(u,v)=F(u,v)+\frac{N(u,v)}{H(u,v)}\text{\hspace{1em}(5.79)}$$

即使知道退化函数，也不能准确地复原未退化的图像。当退化函数$H(u,v)$很小或是零值时时，比例很容易支配了$F(u,v)$。

解决的方法是，将滤波器频率限制到接近原点的值。

$$\hat{F}(u,v)=\frac{G(u,v)}{H(u,v)}\text{\hspace{1em}(5.78)}$$

def freq_cut_off(huv, radius=10):
    """
    set the frequency cut off of the degraded fuction 
    param: huv: input degraded fuction
    param: radius: the cut of radius of frequency, default is 10
    """
    M, N = huv.shape[:2]
    huv_00 = huv[M//2, N//2]
    D0 = radius
    U = np.arange(M)
    V = np.arange(N)
    u, v = np.meshgrid(U, V)
    D = np.sqrt((u - M//2)**2 + (v - N//2)**2)
    huv_1 = huv.copy()
    huv_1[D >= D0] = huv_00
    return huv_1

# 逆滤波，这个是自己想的方法，但效果不太一致
def get_degenerate_image(img, img_deg):
    """
    不填充图像做傅里叶变换后与退化函数做乘积，再反傅里叶变换
    """
    # FFT--------------------------------------------
    fft = np.fft.fft2(img)

    # FFT * H(u, v)----------------------------------
    fft_huv = fft * img_deg

    # IFFT-------------------------------------------
    ifft = np.fft.ifft2(fft_huv)
    
    return ifft

img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0525(a)(aerial_view_no_turb).tif', 0)
M, N = img_ori.shape[:2]
# 退化的图像
k = [0.0025, 0.001, 0.00025]
fp_cen = centralized_2d(img_ori)
huv = modeling_degrade(fp_cen, k=k[0])
ifft = get_degenerate_image(fp_cen, huv)
gxy = centralized_2d(ifft.real)
gxy = np.clip(gxy, 0, gxy.max())
gxy = np.uint8(normalize(gxy) * 255)

# 对退化的的图像中心化
gxy_cen = centralized_2d(gxy)
# 逆滤波
guv = np.fft.fft2(gxy_cen)
fuv = guv / huv
ifft_deg = np.fft.ifft2(fuv)
fxy = centralized_2d(ifft_deg.real)
fxy = np.clip(fxy, 0, fxy.max())
fxy = np.uint8(normalize(fxy) * 255)

# 半径40之外 H 截止
blpf = butterworth_low_pass_filter(huv, huv.shape, radius=40, n=10)
huv_40 = freq_cut_off(huv, radius=40)
fuv_40 = guv / huv_40
ifft_deg = np.fft.ifft2(fuv_40)
fxy_40 = centralized_2d(ifft_deg.real)
fxy_40 = np.clip(fxy_40, 0, fxy_40.max())
fxy_40 = np.uint8(normalize(fxy_40) * 255)

# 半径70之外 H 截止
huv_70 = freq_cut_off(huv, radius=70)
fuv_70 = guv / huv_70
ifft_deg = np.fft.ifft2(fuv_70)
fxy_70 = centralized_2d(ifft_deg.real)
fxy_70 = np.clip(fxy_70, 0, fxy_70.max())
fxy_70 = np.uint8(normalize(fxy_70) * 255)

# 半径85之外 H 截止
huv_85 = freq_cut_off(huv, radius=85)
fuv_85 = guv / huv_85
ifft_deg = np.fft.ifft2(fuv_85)
fxy_85 = centralized_2d(ifft_deg.real)
fxy_85 = np.clip(fxy_85, 0, fxy_85.max())
fxy_85 = np.uint8(normalize(fxy_85) * 255)

fig = plt.figure(figsize=(12, 12))
ax_1 = fig.add_subplot(2, 2, 1, xticks=[], yticks=[])
ax_1.imshow(fxy, 'gray'), ax_1.set_title("H(u, v)")
ax_2 = fig.add_subplot(2, 2, 2, xticks=[], yticks=[])
ax_2.imshow(fxy_40, 'gray'), ax_2.set_title("Radius = 40")
ax_3 = fig.add_subplot(2, 2, 3, xticks=[], yticks=[])
ax_3.imshow(fxy_70, 'gray'), ax_3.set_title("Radius = 70")
ax_4 = fig.add_subplot(2, 2, 4, xticks=[], yticks=[])
ax_4.imshow(fxy_85, 'gray'), ax_4.set_title("Radius = 85")
plt.tight_layout()
plt.show() 

# 逆滤波，这跟书上说的方法比较一致
def get_degenerate_image(img, img_deg):
    """
    不填充图像做傅里叶变换后与退化函数做乘积，再反傅里叶变换
    """
    # FFT--------------------------------------------
    fft = np.fft.fft2(img)

    # FFT * H(u, v)----------------------------------
    fft_huv = fft * img_deg

    # IFFT-------------------------------------------
    ifft = np.fft.ifft2(fft_huv)
    
    return ifft

img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0525(a)(aerial_view_no_turb).tif', 0)
M, N = img_ori.shape[:2]
# 退化的图像
k = [0.0025, 0.001, 0.00025]
fp_cen = centralized_2d(img_ori)
huv = modeling_degrade(fp_cen, k=k[0])
ifft = get_degenerate_image(fp_cen, huv)
gxy = centralized_2d(ifft.real)
gxy = np.clip(gxy, 0, gxy.max())
gxy = np.uint8(normalize(gxy) * 255)

# 对退化的的图像中心化
gxy_cen = centralized_2d(gxy)
# 逆滤波
guv = np.fft.fft2(gxy_cen)
fuv = guv / huv
ifft_deg = np.fft.ifft2(fuv)
fxy = centralized_2d(ifft_deg.real)
fxy = np.clip(fxy, 0, fxy.max())
fxy = np.uint8(normalize(fxy) * 255)

# 半径40之外 H 截止
blpf = butterworth_low_pass_filter(huv, huv.shape, radius=40, n=10)
# huv_40 = freq_cut_off(huv, radius=40)
fuv_40 = guv / huv * blpf
ifft_deg = np.fft.ifft2(fuv_40)
fxy_40 = centralized_2d(ifft_deg.real)
fxy_40 = np.clip(fxy_40, 0, fxy_40.max())
fxy_40 = np.uint8(normalize(fxy_40) * 255)

# 半径70之外 H 截止
blpf = butterworth_low_pass_filter(huv, huv.shape, radius=70, n=11)
# huv_70 = freq_cut_off(huv, radius=70)
fuv_70 = guv / huv * blpf
ifft_deg = np.fft.ifft2(fuv_70)
fxy_70 = centralized_2d(ifft_deg.real)
fxy_70 = np.clip(fxy_70, 0, fxy_70.max())
fxy_70 = np.uint8(normalize(fxy_70) * 255)

# 半径85之外 H 截止
blpf = butterworth_low_pass_filter(huv, huv.shape, radius=85, n=10)
# huv_85 = freq_cut_off(huv, radius=85)
fuv_85 = guv / huv * blpf
ifft_deg = np.fft.ifft2(fuv_85)
fxy_85 = centralized_2d(ifft_deg.real)
fxy_85 = np.clip(fxy_85, 0, fxy_85.max())
fxy_85 = np.uint8(normalize(fxy_85) * 255)

fig = plt.figure(figsize=(12, 12))
ax_1 = fig.add_subplot(2, 2, 1, xticks=[], yticks=[])
ax_1.imshow(fxy, 'gray'), ax_1.set_title("H(u, v)")
ax_2 = fig.add_subplot(2, 2, 2, xticks=[], yticks=[])
ax_2.imshow(fxy_40, 'gray'), ax_2.set_title("Radius = 40")
ax_3 = fig.add_subplot(2, 2, 3, xticks=[], yticks=[])
ax_3.imshow(fxy_70, 'gray'), ax_3.set_title("Radius = 70")
ax_4 = fig.add_subplot(2, 2, 4, xticks=[], yticks=[])
ax_4.imshow(fxy_85, 'gray'), ax_4.set_title("Radius = 85")
plt.tight_layout()
plt.show()