頻域濾波器
阿新 • • 發佈：2020-07-26
傅立葉變換可以得到影象的頻率資訊，頻率資訊代表了影象畫素變化的快慢，通過對影象頻率資訊的處理，可以實現影象濾波的效果。
頻域濾波的步驟如下
（1）給定一幅M×N大小的影象$f(x,y)$，選擇填充引數$P=2M$，$Q=2N$。
（2）對$f(x,y)$填充必要的0，形成大小為M×N的填充影象$f_p(x,y)$。
（3）用$(-1)^{x+y}$乘以$f_p(x,y)$移到其變換的中心。
（4）計算來自步驟3影象的DFT，得到$F(u,v)。
（5）生成一個實的對稱的濾波函式$H(u,v)$，其大小為$P*Q$，中心在$(\frac{M}{2},\frac{N}{2})$。用陣列相乘得到$G(u,v)=F(u,v)*H(u,v)$。
（6）得到處理後的影象
（7）通過從$f_p(x,y)$左上象限提取c處理結果$g(x,y)$。
進行影象填充和頻率中心平移的自定義函式如下
# 將影象填充到size大小，預設為填充到原影象的兩倍 def my_fill(image, size=(-1,-1)): raw = 2 * image.shape[0] column = 2 * image.shape[1] if size[0] > 0: raw, column = size # 這裡的資料結構可能更改 filled_image = np.zeros((raw, column), dtype=np.int64)for i in range(image.shape[0]): for j in range(image.shape[1]): filled_image[i][j] = image[i][j] return filled_image # 影象從offset擷取size大小，預設從左上角擷取一半 def my_ifill(image, size=(-1, -1), offset=(0, 0)): raw = int(image.shape[0] / 2) column = int(image.shape[1] / 2) if size[0] > 0: raw, column= size ifilled_image = np.zeros((raw, column), dtype=np.int64) for i in range(raw): for j in range(column): ifilled_image[i][j] = round(image[i+offset[0]][j+offset[1]]) return ifilled_image
# 對畫素乘以(-1)^(x+y) def my_get_fp(image): image_fp = image.copy() raw, column = image.shape for i in range(raw): for j in range(column): image_fp[i][j] = image[i][j] * pow(-1, i + j) return image_fp
按照上述處理步驟自定義的頻域濾波函式如下
# 自定義模板頻域濾波函式 def myfunc_seldifine(image, model, output_offset=(0, 0)): """ 用於得到頻域濾波的結果，濾波模板需要自己設計 :param image: 輸入影象，不用填充 :param model: 頻域濾波模板，需要將自己設計好的模板輸入，尺寸大於影象尺寸。在函式中影象會根據模板大小進行填充 :param output_offset: 在進行反變換後的大尺寸影象擷取時的偏移量 :return: 濾波後的輸出影象，尺寸與原影象相同 """ # 填充影象 image_fill = my_fill(image, size=model.shape) # 將頻率中心移到中點 image_fp = my_get_fp(image_fill) # 求解fft image_fft = np.fft.fft2(image_fp) fft_cal = model * image_fft # 進行傅立葉反變換 image_ifft = np.fft.ifft2(fft_cal) # 將影象移到中心點 image_ifft_real = image_ifft.real image_ifp = my_get_fp(image_ifft_real) result = my_ifill(image_ifp, size=image.shape, offset=output_offset) return result # 典型頻域濾波函式 def myfunc_filter(image, d0=60, filetype = filter_type.butterworth_low_pass, n=1): """ 用於得到頻域濾波的結果，可以選擇典型的濾波型別。 :param image: 輸入源影象 :param d0: 濾波範圍 :param filetype: 濾波型別，有以下幾種 filter_type.butterworth_low_pass：巴特沃斯低通 filter_type.butterworth_high_pass：巴特沃斯高通 filter_type.gaussion_low_pass：高斯低通 filter_type.gaussion_high_pass：高斯高通 filter_type.laplace_pass：Laplace濾波 :param n: 濾波器的階數 :return: 濾波後的輸出影象，和原影象尺寸相同。 """ # 填充影象 image_fill = my_fill(image) # 將頻率中心移到中點 image_fp = my_get_fp(image_fill) # 求解fft image_fft = np.fft.fft2(image_fp) # 產生巴特沃斯濾波模板 model = my_get_butterworth_model(image_fft.shape, d0, n) if filetype == filter_type.butterworth_high_pass: model = my_get_butterworth_highpass_model(image_fft.shape, d0, n) elif filetype == filter_type.gaussion_low_pass: model = my_get_gaussion_model(image_fft.shape, d0) elif filetype == filter_type.gaussion_high_pass: model = my_get_gaussion_highpass_model(image_fft.shape, d0) elif filetype == filter_type.laplace_pass: model = my_get_laplace_model(image_fft.shape) # 進行頻域濾波 fft_cal = model * image_fft # 進行傅立葉反變換 image_ifft = np.fft.ifft2(fft_cal) # 將影象移到中心點 image_ifft_real = image_ifft.real image_ifp = my_get_fp(image_ifft_real) result = my_ifill(image_ifp) if filetype == filter_type.laplace_pass: result = result / np.max(result) * 255 # 這裡有一個問題，就是計算出的影象有沒有負的畫素 result = cv.convertScaleAbs(result) return result
進行頻域濾波的一個關鍵步驟是確當恰當的截止頻率，可以通過規定影象總功率$P_T$的圓來建立一組標準的截止頻率軌跡。
低通濾波器
理想低通濾波器
在以原點為中心的一個圓內，無衰減地通過所有頻率，而圓外的頻率全部濾除掉的濾波器為理想低通濾波器。理想低通濾波器公式如下。
理想低通濾波器雖然在離散域容易實現，但是濾波後的影象會有明顯的振鈴現象，所以實際中不常使用。理想高通濾波器類似。
布特沃斯低通濾波器
n階布特沃斯濾波器定義為
由於布特沃斯濾波器在截止頻率處為平滑過渡，所以基本上沒有振鈴現象。
根據定義計算布特沃斯濾模板函式為
# 得到巴特沃斯濾波模板 def my_get_butterworth_model(size, d0, n=1): model = np.zeros(size, dtype=np.float64) wide, height = size for i in range(wide): for j in range(height): length = math.sqrt((i - wide / 2) ** 2 + (j - height / 2) ** 2) model[i][j] = 1 / (1 + pow(length / d0, 2 * n)) return model
高斯低通濾波器
二維高斯低通濾波器定義為
其中$D_0$是截止頻率。高斯低通濾波器在截止頻率出同樣是平滑過渡，所以也沒有振鈴現象。
高斯低通濾波器模板函式為
# 得到高斯濾波模板 def my_get_gaussion_model(size, d0): model = np.zeros(size, dtype=np.float64) wide, height = size for i in range(wide): for j in range(height): length2 = (i - wide / 2) ** 2 + (j - height / 2) ** 2 model[i][j] = math.exp(-1 * length2 / (2 * d0 ** 2)) return model
布特沃斯高通濾波器
截止頻率為$D_0$的n階布特沃斯高通濾波器定義為
與布特沃斯低通濾波器相比$D_0$和$D(u,v)$的位置相反。模板函式為
# 得到巴特沃斯高通濾波模板 def my_get_butterworth_highpass_model(size, d0, n=1): model = np.zeros(size, dtype=np.float64) wide, height = size for i in range(wide): for j in range(height): length = math.sqrt((i - wide / 2) ** 2 + (j - height / 2) ** 2) + 0.00000001 model[i][j] = 1 / (1 + pow(d0 / length, 2 * n)) return model
高斯高通濾波器
高斯高通濾波器定義為
模板函式
# 得到高斯高通濾波模板 def my_get_gaussion_highpass_model(size, d0): model = np.zeros(size, dtype=np.float64) wide, height = size for i in range(wide): for j in range(height): length2 = (i - wide / 2) ** 2 + (j - height / 2) ** 2 model[i][j] = 1 - math.exp(-1 * length2 / (2 * d0 ** 2)) return model
頻域Laplace濾波器
通過變換可以將空域的Laplace運算元變換到頻域，得到頻域Laplace運算元如下
在對將低頻平移到中心的影象進行處理時，用以下公式計算
計算Laplace濾波模板的函式為
# 得到拉普拉斯濾波模板 def my_get_laplace_model(size): model = np.zeros(size, dtype=np.float64) wide, height = size for i in range(wide): for j in range(height): length2 = (i - wide / 2) ** 2 + (j - height / 2) ** 2 model[i][j] = -4*pow(math.pi, 2) * length2 return model
注意到影象經過Laplace濾波之後，得到的是影象的邊緣資訊，而且資訊畫素值有正有負。如果想要顯示邊緣資訊，需要將得到的畫素值取絕對值。如果想要通過邊緣資訊對影象進行銳化，則不需要取絕對值，直接用原影象減去邊緣資訊的 一定倍數即可。
進行Laplace銳化的函式為
# 拉普拉斯銳化 def myfunc_laplace_sharpen(image): laplace_result = myfunc_seldifine(image, my_get_laplace_model((2 * image.shape[0], 2 * image.shape[1]))) laplace_result = laplace_result / np.max(laplace_result) * 255 result = np.zeros(image.shape, dtype=np.int64) for i in range(image.shape[0]): for j in range(image.shape[1]): result[i][j] = int(image[i][j]) - int(laplace_result[i][j]) sharpen_image = cv.convertScaleAbs(result) return sharpen_image