python實現直方圖均衡化，理想高通濾波與高斯低通濾波

寫在前面

HIT大三上學期視聽覺信號處理課程中視覺部分的實驗二，經過和學長們實驗的對比發現每一級實驗要求都不一樣，因此這里標明了是2019年秋季學期的視覺實驗二。

由於時間緊張，代碼沒有進行任何優化，實驗算法僅供參考。

實驗要求

1. 實現圖像直方圖均衡化，要求顯示均衡化前、后直方圖以及均衡化后圖像。

2. 對單通道圖像進行DFT變換，要求顯示幅度圖和相位圖，並設計理想高通濾波器和高斯低通濾波器對圖像進行頻域濾波，並顯示濾波之后的圖像。 注：除DFT和IDFT外，不允許調庫

實驗代碼

代碼首先貼在這里，僅供參考

源代碼

首先是實現圖像直方圖均衡化的代碼

from ReadBmp import ReadBmp
import matplotlib.pyplot as plt
import  numpy as np

filename1 = "1.bmp"
bmp = ReadBmp(filename1)
bmp.gray()

# 統計各像素點數
h = np.array([0 for i in range(256)])
h1 = []
for pixel in bmp.data:
    h[pixel[0]] = h[pixel[0]] + 1
    h1.append(pixel[0])
# 畫出原先的直方圖
plt.subplot(1,2,1)
plt.hist(h1, bins = 256)

# 歸一化
hs = h / len(bmp.data)
# 計算累計分布
hp = np.array([0.0 for i in range(256)])
for i in range(256):
    hp[i] = np.round(np.sum(hs[0:i+1]) * 255)
T = hp.astype('uint8')

# 創建新圖像，並統計新圖像的各個像素點的個數
hn = np.array([0 for i in range(256)])
h2 = []
for pixel in bmp.data:
    s = T[pixel[0]]
    pixel[0] = s
    pixel[1] = s
    pixel[2] = s
    hn[pixel[0]] = hn[pixel[0]] + 1
    h2.append(s)
bmp.creataBmp("2.bmp")
# 畫出新圖像的直方圖
plt.subplot(1,2,2)
plt.hist(h2, bins = 256)
plt.show()

其中讀取bmp圖像的程序是我自己寫的，這里不再贅述（直接調用了實驗一寫好的bmp文件讀取程序，具體代碼見github）

然后是理想高通濾波和高斯低通濾波的程序

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

def IdealHighPassFiltering(f_shift):
    # 設置濾波半徑
    D0 = 8
    # 初始化
    m = f_shift.shape[0]
    n = f_shift.shape[1]
    h1 = np.zeros((m, n))
    x0 = np.floor(m/2)
    y0 = np.floor(n/2)
    for i in range(m):
        for j in range(n):
            D = np.sqrt((i - x0)**2 + (j - y0)**2)
            if D >= D0:
                h1[i][j] = 1
    result = np.multiply(f_shift, h1)
    return result

def GaussLowPassFiltering(f_shift):
    # 設置濾波半徑
    D0 = 8
    # 初始化
    m = f_shift.shape[0]
    n = f_shift.shape[1]
    h1 = np.zeros((m, n))
    x0 = np.floor(m/2)
    y0 = np.floor(n/2)
    for i in range(m):
        for j in range(n):
            D = np.sqrt((i - x0)**2 + (j - y0)**2)
            h1[i][j] = np.exp((-1)*D**2/2/(D0**2))
    result = np.multiply(f_shift, h1)
    return result

img =cv2.imread('1.bmp',0)
f=np.fft.fft2(img)
f_shift=np.fft.fftshift(f)
# 幅度圖
s= np.log(abs(f_shift))
# 相位圖
p= abs(np.angle(f_shift))
plt.subplot(2,2,1)
plt.imshow(s, 'gray')
plt.subplot(2,2,2)
plt.imshow(p, 'gray')

# 理想高通濾波
IHPF = IdealHighPassFiltering(f_shift)
new_f1 = np.fft.ifftshift(IHPF)
new_image1 = np.uint8(np.abs(np.fft.ifft2(new_f1)))
plt.subplot(2,2,3)
plt.imshow(new_image1, 'gray')
# 高斯低通濾波
GLPF = GaussLowPassFiltering(f_shift)
new_f2 = np.fft.ifftshift(GLPF)
new_image2 = np.uint8(np.abs(np.fft.ifft2(new_f2)))
plt.subplot(2,2,4)
plt.imshow(new_image2, 'gray')
plt.show()

實驗結果

實驗結果的圖像如下圖所示

這是直方圖均衡化的運行結果，其中左圖是做直方圖均衡化前對像素統計的直方圖，右圖是做直方圖均衡化后對像素統計的直方圖

這是畫出圖像的幅度圖和相位圖，對圖像做理想高通濾波和高斯低通濾波的結果，其中1是幅度圖，2是相位圖，3是理想高通濾波的結果，4是高斯低通濾波的結果

要點講解

關於實驗算法的大致流程，已在程序中用注釋標明，這里只記錄需要注意的地方

1. 在直方圖均衡化中，要注意直方圖函數的用法，該函數可以對輸入的矩陣中的數值進行自動的分類，而不用手動分好之后再用直方圖函數去畫，之前沒有仔細看該函數的用法，導致找了好長時間沒找到問題。

2. 關於直方圖均衡化的原理，需要一定的數學公式的推導，然而其結論十分簡單，只有一個公式，對理解公式推導無需求的直接拿公式來用就可以了。（即程序中的歸一化加計算累積分布部分）

3. 關於理想高通濾波和高斯低通濾波，可以看到其代碼大體框架是一致的，只有中間if語句中的判斷條件不一樣，理解原理之后直接替換if語句即可。

4. 最終畫出理想高通濾波的圖像時，網上給出的一些運行結果可能和我們運行的不一致，類似於下圖中圖一的圖像，原因在於這一句話中：new_image1 = np.uint8(np.abs(np.fft.ifft2(new_f1)))，注意到中間調用了np.abs()函數，由於經過DFT變換之后所得的數值是復數，而不是實數，若將np.abs()替換為np.real()，即只取其實部，則得到的是圖一，而使用np.abs()，得到的則是圖二。個人感覺圖二才是相對正確一些的結果，圖一相當於將虛部舍棄掉了，這是不太好的。

以下是圖一

以下是圖二

知識總結

關於直方圖均衡化的數學原理，我參考的是這篇博客。直方圖均衡化的公式推導較為復雜，然而其結果卻十分簡單，因此在python中一個循環就實現了。
直方圖均衡化的數學原理

關於如何畫直方圖，參考的是這兩篇博客
python之matplotlib.pyplot直方圖

關於相位圖和振幅圖的畫法，參考的是下面這兩篇博客
python opencv dft
進行傅里葉變換，畫出幅值譜圖、相位譜圖以及實現位移后的幅值譜圖和相位譜圖
案例解釋圖像傅里葉變換的幅度譜和相位譜的以及反變換

關於理想高通和高斯低通濾波器的實現，參考的是以下這幾篇博客
理想高通濾波器
理想低通濾波器、Butterworth濾波器和高斯濾波器(matlab)

其中在理想高通和高斯低通濾波器中，需要用到矩陣乘法，關於幾種不同的矩陣乘法參考的是這篇博客
numpy 三種矩陣乘法

實驗總結

很簡單的一個實驗，總時間加起來大致4個小時左右，做這么快的主要原因是，本來同學來哈爾濱找我玩，正玩得開心，結果突然告訴我這個實驗要今天晚上交（之前以為這個實驗還有一周的ddl），所以潛能爆發，從三點寫到六點差不多就寫完了