1、頻率域與空間域之間的關系
在頻率域濾波基礎——傅里葉變換計算及應用基礎中，我們已經知道了如何將圖像轉換到頻率域，以及如何將頻率域圖像通過傅里葉逆變換轉換回圖像，這樣一來，可以利用空域圖像與頻譜之間的對應關系，嘗試將空域卷積濾波變換為頻域濾波，而后再將頻域濾波處理后的圖像反變換回空間域，從而達到圖像增強的目的，這樣做的一個最主要的吸引力再域頻率域濾波具有直觀性的特點。
根據著名的卷積定律，兩個二維連續函數再空間域的卷積可由其相應的兩個傅里葉變換乘積的反變換而得到，反之，在頻率域中的卷積可由在空間域中乘積的傅里葉變換而得到。即
$$
\begin{array}{l}{f_{1}(t) \leftrightarrow F_{1}(\omega), f_{2}(t) \leftrightarrow F_{2}(\omega)} \\ {f_{1}(t) * f_{2}(t) \leftrightarrow F_{1}(\omega) \cdot F_{2}(\omega)}\end{array}
$$
2、頻率域濾波的基本步驟
本文重點就來講一講如何通過頻率域濾波來實現圖像的平滑和銳化。首先將頻率域濾波的步驟來做一個小結如下：
- 給定一幅大小為MXN的輸入圖像f(x,y),選擇填充參數P=2M.Q=2N.
- 對f(x,y)添加必要數量的0，形成大小為PXQ的填充圖像$f_{p}(x, y)$。
- 用$(-1)^{x+y}$乘以$f_{p}(x, y)$，移到其變換中心。
- 計算上圖的DFT，得到F(u,v).
- 生成一個實的，對稱的濾波函數$H(u, v)$，其大小為PXQ，中心在$(P / 2, Q / 2)$處。用陣列相乘得到乘積$G(u, v)=H(u, v) F(u, v)$；即$G(i, k)=H(i, k) F(i, k)$
- 得到處理后的圖像：
$$
\mathrm{g}_{p}(x, y)=\left\{\text {real}\left[\mathfrak{J}^{-1}[G(u, v)]\right\}(-1)^{x+y}\right.
$$
其中，為忽略由於計算不准確導致的寄生復成分，選擇實部，下標p指出我們處理的是填充后的圖像。
- 從$g_{p}(x, y)$的左上象限提取MXN區域，得到最終處理結果$g(x, y)$
由上述敘述可知，濾波的關鍵取決於濾波函數$H(u, v)$，常常稱為濾波器，或濾波器傳遞函數，因為它在濾波中抑制或除去了頻譜中的某些分量，而保留其他的一些頻率不受影響，從而達到濾波的目的。而本節將講解一些常用的濾波器。
3、使用頻率域濾波平滑圖像
在空間域我們已經講過圖像平滑的目的及空間域平滑濾波，現在在頻率域濾波中我們首先講解三種頻率平滑濾波器：理想濾波器、巴特沃斯濾波器、高斯濾波器。這三種濾波器涵蓋了從非常急劇（理想）到非常平滑的濾波范圍。
1、理想低通濾波器（ILPF)
首先最容易想到得衰減高頻成分得方法就是在一個稱為截止頻率得位置截斷所有的高頻成分，將圖像頻譜中所有高於這一截止頻率的頻譜成分設為0，低於截止頻率的成分保持不變。如圖所示（左圖是其透視圖，右圖是其圖像表示）：
能夠達到這種用效果的濾波器我們稱之為理想濾波器，用公式表示為：
$$
H(u, v)=\left\{\begin{array}{ll}{1} & {\text { if } D(u, v) \leq D_{0}} \\ {0} & {\text { if } D(u, v)>D_{0}}\end{array}\right.
$$
其中，D0表示通帶的半徑。D(u,v)的計算方式也就是兩點間的距離，很簡單就能得到。
$$
D(u, v)=\sqrt{\left(u-\frac{P}{2}\right)^{2}+\left(v-\frac{Q}{2}\right)^{2}}
$$
算法在OpenCV中的實現如下所示：
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>
using namespace cv;
using namespace std;
int main()
{
//Mat input = imread(p[1], CV_LOAD_IMAGE_GRAYSCALE);
//這是在ubuntu上運行的，p[1]為控制台給出的參數，即圖片路徑
//如果不知道怎么傳入參數，可以直接改為
Mat input=imread("2222.jpg",CV_LOAD_IMAGE_GRAYSCALE);
//image.jpg必須放在當前目錄下，image.jpg即為輸入圖片
imshow("input", input);
int w = getOptimalDFTSize(input.cols);
int h = getOptimalDFTSize(input.rows);
Mat padded;
copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));
padded.convertTo(padded, CV_32FC1);
imshow("padded", padded);
for (int i = 0; i < padded.rows; i++)//中心化操作，其余操作和上一篇博客的介紹一樣
{
float *ptr = padded.ptr<float>(i);
for (int j = 0; j < padded.cols; j++) ptr[j] *= pow(-1, i + j);
}
Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
Mat complexImg;
merge(plane, 2, complexImg);
dft(complexImg, complexImg);
//****************************************************
Mat mag = complexImg;
mag = mag(Rect(0, 0, mag.cols & -2, mag.rows & -2));//這里為什么&上-2具體查看opencv文檔
//其實是為了把行和列變成偶數 -2的二進制是11111111.......10 最后一位是0
//獲取中心點坐標
int cx = mag.cols / 2;
int cy = mag.rows / 2;
float D0 = 20;
for (int y = 0; y < mag.rows; y++)
{
double* data = mag.ptr<double>(y);
for (int x = 0; x < mag.cols; x++)
{
double d = sqrt(pow((y - cy), 2) + pow((x - cx), 2));
if (d <= D0) {
//小於截止頻率不做操作
}
else {
data[x] = 0;
}
}
}
//******************************************************************
split(mag, plane);
magnitude(plane[0], plane[1], plane[0]);
plane[0] += Scalar::all(1);
log(plane[0], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("mag", plane[0]);
//******************************************************************
//*************************idft*************************************
idft(mag, mag);
split(mag, plane);
magnitude(plane[0], plane[1], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("idft-mag", plane[0]);
waitKey();
return 0;
}
從左至右依次是：原圖，優化過后的圖、傅里葉變換圖、理想低通濾波器處理過后的高斯變換圖、理想低通濾波器處理過后的效果圖：
2、巴特沃思低通濾波器(BLPF）
截止頻率位於距離原點$D_{0}$處的n階巴特沃斯濾波器的傳遞函數定義為：
$$
H(u, v)=\frac{1}{1+\left[D(u, v) / D_{0}\right]^{2 n}}
$$
式中
$$
D(u, v)=\sqrt{\left(u-\frac{P}{2}\right)^{2}+\left(v-\frac{Q}{2}\right)^{2}}
$$
下圖即為巴特沃斯濾波器的透視圖及其圖像顯示：
與理想低通濾波不同，巴特沃斯低通濾波器並沒有再通過頻率和濾除頻率之間給出明顯截止的急劇不確定性。對於具有平滑傳遞函數的濾波器，可在這樣一點上定義截止頻率，即使H(u,v)下降為其最大值的某個百分比的點。
同理，我們可以在OpenCV下實現相關算法：
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>
using namespace cv;
using namespace std;
int main()
{
//Mat input = imread(p[1], CV_LOAD_IMAGE_GRAYSCALE);
//這是在ubuntu上運行的，p[1]為控制台給出的參數，即圖片路徑
//如果不知道怎么傳入參數，可以直接改為
Mat input=imread("2222.jpg",CV_LOAD_IMAGE_GRAYSCALE);
//image.jpg必須放在當前目錄下，image.jpg即為輸入圖片
imshow("input", input);
//Butterworth_Low_Paass_Filter(input);
int w = getOptimalDFTSize(input.cols);
int h = getOptimalDFTSize(input.rows);
Mat padded;
copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));
padded.convertTo(padded, CV_32FC1);
imshow("padded", padded);
for (int i = 0; i < padded.rows; i++)//中心化操作，其余操作和上一篇博客的介紹一樣
{
float *ptr = padded.ptr<float>(i);
for (int j = 0; j < padded.cols; j++) ptr[j] *= pow(-1, i + j);
}
Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
Mat complexImg;
merge(plane, 2, complexImg);
//****************************************************
dft(complexImg, complexImg);
split(complexImg, plane);
magnitude(plane[0], plane[1], plane[0]);
plane[0] += Scalar::all(1);
log(plane[0], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("complexImg", plane[0]);
//****************************************************
Mat mag = complexImg;
mag = mag(Rect(0, 0, mag.cols & -2, mag.rows & -2));//這里為什么&上-2具體查看opencv文檔
//其實是為了把行和列變成偶數 -2的二進制是11111111.......10 最后一位是0
//獲取中心點坐標
int cx = mag.cols / 2;
int cy = mag.rows / 2;
int n = 10;//表示巴特沃斯濾波器的次數
//H = 1 / (1+(D/D0)^2n)
double D0 = 10;
for (int y = 0; y < mag.rows; y++) {
double* data = mag.ptr<double>(y);
for (int x = 0; x < mag.cols; x++) {
//cout << data[x] << endl;
double d = sqrt(pow((y - cy), 2) + pow((x - cx), 2));
//cout << d << endl;
double h = 1.0 / (1 + pow(d / D0, 2 * n));
if (h <= 0.5) {
data[x] = 0;
}
else {
//data[x] = data[x]*0.5;
//cout << h << endl;
}
//cout << data[x] << endl;
}
}
//******************************************************************
split(mag, plane);
magnitude(plane[0], plane[1], plane[0]);
plane[0] += Scalar::all(1);
log(plane[0], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("mag", plane[0]);
//******************************************************************
//*************************idft*************************************
idft(mag, mag);
split(mag, plane);
magnitude(plane[0], plane[1], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("idft-mag", plane[0]);
waitKey();
return 0;
}
從左至右依次是：原圖，優化過后的圖、傅里葉變換圖、巴特沃斯濾波器平滑后的高斯變換圖、巴特沃斯濾波器平滑后的效果圖：
3、高斯低通濾波器(GLPF)
高斯低通濾波的頻率域二維形式由下式給出：
$$
H(u, v)=e^{-D^{2}(u, v) / 2 D_{0}^{2}}
$$
式中
$$
D(u, v)=\sqrt{\left(u-\frac{P}{2}\right)^{2}+\left(v-\frac{Q}{2}\right)^{2}}
$$
是離頻率矩陣中心的距離。$D_{0}$是截止頻率，當$D(u, v)=D_{0}$時，GLPF下降到其最大值的0.607處。
高斯濾波器的透視圖及其圖像顯示如下圖所示：
在OpenCV下實現高斯低通濾波算法如下：
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>
using namespace cv;
using namespace std;
int main()
{
//Mat input = imread(p[1], CV_LOAD_IMAGE_GRAYSCALE);
//這是在ubuntu上運行的，p[1]為控制台給出的參數，即圖片路徑
//如果不知道怎么傳入參數，可以直接改為
Mat input=imread("2222.jpg",CV_LOAD_IMAGE_GRAYSCALE);
//image.jpg必須放在當前目錄下，image.jpg即為輸入圖片
imshow("input", input);
int w = getOptimalDFTSize(input.cols);
int h = getOptimalDFTSize(input.rows);
Mat padded;
copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));
padded.convertTo(padded, CV_32FC1);
imshow("padded", padded);
for (int i = 0; i < padded.rows; i++)//中心化操作，其余操作和上一篇博客的介紹一樣
{
float *ptr = padded.ptr<float>(i);
for (int j = 0; j < padded.cols; j++) ptr[j] *= pow(-1, i + j);
}
Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
Mat complexImg;
merge(plane, 2, complexImg);
dft(complexImg, complexImg);
//************************gaussian****************************
Mat gaussianBlur(padded.size(), CV_32FC2);
float D0 = 2 * 10 * 10;
for (int i = 0; i < padded.rows; i++)
{
float*p = gaussianBlur.ptr<float>(i);
for (int j = 0; j < padded.cols; j++)
{
float d = pow(i - padded.rows / 2, 2) + pow(j - padded.cols / 2, 2);
p[2 * j] = expf(-d / D0);
p[2 * j + 1] = expf(-d / D0);
}
}
multiply(complexImg, gaussianBlur, gaussianBlur);//矩陣元素對應相乘法，注意，和矩陣相乘區分
//*****************************************************************
split(complexImg, plane);
magnitude(plane[0], plane[1], plane[0]);
plane[0] += Scalar::all(1);
log(plane[0], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("dft", plane[0]);
//******************************************************************
split(gaussianBlur, plane);
magnitude(plane[0], plane[1], plane[0]);
plane[0] += Scalar::all(1);
log(plane[0], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("gaussianBlur", plane[0]);
//******************************************************************
//*************************idft*************************************
idft(gaussianBlur, gaussianBlur);
split(gaussianBlur, plane);
magnitude(plane[0], plane[1], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("idft-gaussianBlur", plane[0]);
waitKey();
return 0;
}
從左至右依次是：原圖，優化過后的圖、傅里葉變換圖、高斯低通濾波后的高斯變換圖、高斯平滑過后的效果圖：
4、使用頻率域濾波銳化圖像
前面已經講過了通過衰減圖像的傅里葉變換的高頻成分可以平滑圖像。因為邊緣和其他灰度的急劇變化與高頻成分有關，所以圖像的銳化可在頻率域通過濾波來實現，高通濾波器會衰減傅里葉變換中的低頻成分而不擾亂高頻信息。故我們可以考慮用高通濾波器來進行圖像的銳化,一個高頻濾波器是從給定低通濾波器用下式得到的:
$$
H_{\mathrm{HP}}(u, v)=1-H_{\mathrm{LP}}(u, v)
$$
與低通濾波器對應，這里介紹理想高通濾波器，巴特沃斯高通濾波器，高斯高通濾波器三種高通濾波器。
1、理想高通濾波器 (IHPF)
與低通對應，理想高通濾波器定義為：
$$
H(u, v)=\left\{\begin{array}{ll}{0} & {\text { if } D(u, v) \leq D_{0}} \\ {1} & {\text { if } D(u, v)>D_{0}}\end{array}\right.
$$
透視圖和圖像顯示如下：
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>
using namespace cv;
using namespace std;
int main()
{
//Mat input = imread(p[1], CV_LOAD_IMAGE_GRAYSCALE);
//這是在ubuntu上運行的，p[1]為控制台給出的參數，即圖片路徑
//如果不知道怎么傳入參數，可以直接改為
Mat input = imread("2222.jpg", CV_LOAD_IMAGE_GRAYSCALE);
//image.jpg必須放在當前目錄下，image.jpg即為輸入圖片
imshow("input", input);
int w = getOptimalDFTSize(input.cols);
int h = getOptimalDFTSize(input.rows);
Mat padded;
copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));
padded.convertTo(padded, CV_32FC1);
imshow("padded", padded);
for (int i = 0; i < padded.rows; i++)//中心化操作，其余操作和上一篇博客的介紹一樣
{
float *ptr = padded.ptr<float>(i);
for (int j = 0; j < padded.cols; j++) ptr[j] *= pow(-1, i + j);
}
Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
Mat complexImg;
merge(plane, 2, complexImg);
dft(complexImg, complexImg);
//****************************************************
Mat mag = complexImg;
mag = mag(Rect(0, 0, mag.cols & -2, mag.rows & -2));//這里為什么&上-2具體查看opencv文檔
//其實是為了把行和列變成偶數 -2的二進制是11111111.......10 最后一位是0
//獲取中心點坐標
int cx = mag.cols / 2;
int cy = mag.rows / 2;
float D0 = 20;
for (int y = 0; y < mag.rows; y++)
{
double* data = mag.ptr<double>(y);
for (int x = 0; x < mag.cols; x++)
{
double d = sqrt(pow((y - cy), 2) + pow((x - cx), 2));
if (d > D0) {
//小於截止頻率不做操作
}
else {
data[x] = 0;
}
}
}
//******************************************************************
split(mag, plane);
magnitude(plane[0], plane[1], plane[0]);
plane[0] += Scalar::all(1);
log(plane[0], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("mag", plane[0]);
//******************************************************************
//*************************idft*************************************
idft(mag, mag);
split(mag, plane);
magnitude(plane[0], plane[1], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("idft-mag", plane[0]);
waitKey();
return 0;
}
結果如下：
2、巴特沃思高通濾波器 （BHPF)
同理n階巴特沃斯高通濾波器的定義為
$$
H(u, v)=\frac{1}{1+\left[D_{0} / D(u, v)\right]^{2 n}}
$$
透視圖和圖像顯示如下：
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>
using namespace cv;
using namespace std;
int main()
{
//Mat input = imread(p[1], CV_LOAD_IMAGE_GRAYSCALE);
//這是在ubuntu上運行的，p[1]為控制台給出的參數，即圖片路徑
//如果不知道怎么傳入參數，可以直接改為
Mat input = imread("2222.jpg", CV_LOAD_IMAGE_GRAYSCALE);
//image.jpg必須放在當前目錄下，image.jpg即為輸入圖片
imshow("input", input);
//Butterworth_Low_Paass_Filter(input);
int w = getOptimalDFTSize(input.cols);
int h = getOptimalDFTSize(input.rows);
Mat padded;
copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));
padded.convertTo(padded, CV_32FC1);
imshow("padded", padded);
for (int i = 0; i < padded.rows; i++)//中心化操作，其余操作和上一篇博客的介紹一樣
{
float *ptr = padded.ptr<float>(i);
for (int j = 0; j < padded.cols; j++) ptr[j] *= pow(-1, i + j);
}
Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
Mat complexImg;
merge(plane, 2, complexImg);
//****************************************************
dft(complexImg, complexImg);
split(complexImg, plane);
magnitude(plane[0], plane[1], plane[0]);
plane[0] += Scalar::all(1);
log(plane[0], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("complexImg", plane[0]);
//****************************************************
Mat mag = complexImg;
mag = mag(Rect(0, 0, mag.cols & -2, mag.rows & -2));//這里為什么&上-2具體查看opencv文檔
//其實是為了把行和列變成偶數 -2的二進制是11111111.......10 最后一位是0
//獲取中心點坐標
int cx = mag.cols / 2;
int cy = mag.rows / 2;
int n = 1;//表示巴特沃斯濾波器的次數
//H = 1 / (1+(D/D0)^2n)
double D0 = 30;
for (int y = 0; y < mag.rows; y++) {
double* data = mag.ptr<double>(y);
for (int x = 0; x < mag.cols; x++) {
//cout << data[x] << endl;
double d = sqrt(pow((y - cy), 2) + pow((x - cx), 2));
//cout << d << endl;
double h = 1.0 / (1 + pow(d / D0, 2 * n));
if (h > 0.5) {//低通變高通
data[x] = 0;
}
else {
//data[x] = data[x]*0.5;
//cout << h << endl;
}
//cout << data[x] << endl;
}
}
//******************************************************************
split(mag, plane);
magnitude(plane[0], plane[1], plane[0]);
plane[0] += Scalar::all(1);
log(plane[0], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("mag", plane[0]);
//******************************************************************
//*************************idft*************************************
idft(mag, mag);
split(mag, plane);
magnitude(plane[0], plane[1], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("idft-mag", plane[0]);
waitKey();
return 0;
}
結果如下：
3、高斯高通濾波器（GHPF)
高斯高通濾波器定義為：
$$
H(u, v)=1-\mathrm{e}^{-D^{2}(u, v) / 2 D_{0}^{2}}
$$
透視圖和圖像顯示如下：
OpenCV實現為：
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>
using namespace cv;
using namespace std;
int main()
{
//Mat input = imread(p[1], CV_LOAD_IMAGE_GRAYSCALE);
//這是在ubuntu上運行的，p[1]為控制台給出的參數，即圖片路徑
//如果不知道怎么傳入參數，可以直接改為
Mat input=imread("2222.jpg",CV_LOAD_IMAGE_GRAYSCALE);
//image.jpg必須放在當前目錄下，image.jpg即為輸入圖片
imshow("input", input);
int w = getOptimalDFTSize(input.cols);
int h = getOptimalDFTSize(input.rows);
Mat padded;
copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));
padded.convertTo(padded, CV_32FC1);
imshow("padded", padded);
for (int i = 0; i < padded.rows; i++)//中心化操作，其余操作和上一篇博客的介紹一樣
{
float *ptr = padded.ptr<float>(i);
for (int j = 0; j < padded.cols; j++) ptr[j] *= pow(-1, i + j);
}
Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
Mat complexImg;
merge(plane, 2, complexImg);
dft(complexImg, complexImg);
//************************gaussian****************************
Mat gaussianSharpen(padded.size(), CV_32FC2);
float D0 = 2 * 10 * 10;
for (int i = 0; i < padded.rows; i++)
{
float*q = gaussianSharpen.ptr<float>(i);
for (int j = 0; j < padded.cols; j++)
{
float d = pow(i - padded.rows / 2, 2) + pow(j - padded.cols / 2, 2);
q[2 * j] = 1 - expf(-d / D0);
q[2 * j + 1] = 1 - expf(-d / D0);
}
}
multiply(complexImg, gaussianSharpen, gaussianSharpen);
//*****************************************************************
split(complexImg, plane);
magnitude(plane[0], plane[1], plane[0]);
plane[0] += Scalar::all(1);
log(plane[0], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("dft", plane[0]);
//******************************************************************
split(gaussianSharpen, plane);
magnitude(plane[0], plane[1], plane[0]);
plane[0] += Scalar::all(1);
log(plane[0], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("gaussianSharpen", plane[0]);
//******************************************************************
//*************************idft*************************************
idft(gaussianSharpen, gaussianSharpen);
split(gaussianSharpen, plane);
magnitude(plane[0], plane[1], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("idft-gaussianSharpen", plane[0]);
waitKey();
return 0;
}
從左至右依次是：原圖，優化過后的圖、傅里葉變換圖、高斯銳化后的高斯變換圖、高斯銳化過后的效果圖：
6、選擇率濾波器
在很多應用中，圖像處理的目的是處理指定頻段或頻率矩陣的小區域，此時我們需要使用選擇濾波器對其進行處理，選擇濾波器分為帶通濾波器和帶阻濾波器和陷波器。
6.1、帶阻濾波器、帶通濾波器
所謂的帶通濾波器即是將低頻濾波和高頻濾波相結合，最后我們可以利用帶阻濾波器濾掉周期性出現的噪聲，由前面可知，
理想帶通濾波器公式為
$$
H(\mathrm{u}, v)=\left\{\begin{array}{l}{0 i f D_{0}-\frac{W}{2} \leq D \leq D_{0}+\frac{W}{2}} \\ {1}\end{array}\right.
$$
巴特沃斯帶通濾波器公式為
$$
H(u, v)=\frac{1}{1+\left[\frac{D W}{D^{2}-D_{0}^{2}}\right]^{2 n}}
$$
高斯帶通濾波器公式為：
$$
H(u, v)=1-\mathrm{e}^{-\left[\frac{D^{2}-D_{0}^{2}}{D W}\right]^{2}}
$$
與由低通濾波器得到高通濾波器相同的方法，由帶阻濾波器我們可以得到帶通濾波器
$$
H_{\mathrm{BP}}(u, v)=1-H_{\mathrm{BR}}(u, v)
$$
示例
以下就是理想帶通濾波器的OpenCV實現及效果圖
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>
using namespace cv;
using namespace std;
int main()
{
//Mat input = imread(p[1], CV_LOAD_IMAGE_GRAYSCALE);
//這是在ubuntu上運行的，p[1]為控制台給出的參數，即圖片路徑
//如果不知道怎么傳入參數，可以直接改為
Mat input = imread("2222.jpg", CV_LOAD_IMAGE_GRAYSCALE);
//image.jpg必須放在當前目錄下，image.jpg即為輸入圖片
imshow("input", input);
int w = getOptimalDFTSize(input.cols);
int h = getOptimalDFTSize(input.rows);
Mat padded;
copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, BORDER_CONSTANT, Scalar::all(0));
padded.convertTo(padded, CV_32FC1);
imshow("padded", padded);
for (int i = 0; i < padded.rows; i++)//中心化操作，其余操作和上一篇博客的介紹一樣
{
float *ptr = padded.ptr<float>(i);
for (int j = 0; j < padded.cols; j++) ptr[j] *= pow(-1, i + j);
}
Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
Mat complexImg;
merge(plane, 2, complexImg);
dft(complexImg, complexImg);
//****************************************************
Mat mag = complexImg;
mag = mag(Rect(0, 0, mag.cols & -2, mag.rows & -2));//這里為什么&上-2具體查看opencv文檔
//其實是為了把行和列變成偶數 -2的二進制是11111111.......10 最后一位是0
//獲取中心點坐標
int cx = mag.cols / 2;
int cy = mag.rows / 2;
float D0 = 20;
for (int y = 0; y < mag.rows; y++)
{
double* data = mag.ptr<double>(y);
for (int x = 0; x < mag.cols; x++)
{
double d = sqrt(pow((y - cy), 2) + pow((x - cx), 2));
if ( D0<= d&& d <= D0+30 ) {
}
else {
data[x] = 0;
}
}
}
//******************************************************************
split(mag, plane);
magnitude(plane[0], plane[1], plane[0]);
plane[0] += Scalar::all(1);
log(plane[0], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("mag", plane[0]);
//******************************************************************
//*************************idft*************************************
idft(mag, mag);
split(mag, plane);
magnitude(plane[0], plane[1], plane[0]);
normalize(plane[0], plane[0], 1, 0, CV_MINMAX);
imshow("idft-mag", plane[0]);
waitKey();
return 0;
}
6.2、陷波濾波器
陷波器是一個相對於帶通濾波更有用的選擇濾波器，其拒絕事先定義好的一個關於頻率矩陣中心的一個領域的頻域。零相移濾波器必須是關於遠點對稱的，因此，一個中心位於$\left(u_{0}, v_{0}\right)$的陷波在位置$\left(-u_{0},-v_{0}\right)$必須有一個對應的陷波。陷波帶阻濾波器可以用中心已被平移到陷波濾波器中心的高通濾波器的乘積來構造。一般形式為：
$$
H_{N R}(u, v)=\prod_{k=1}^{Q} H_{k}(u, v) H_{-k}(u, v)
$$
其中，$H_{k}(u, v) H_{-k}(u, v)$是高通濾波器，他們的中心分別位於$\left(u_{k}, v_{k}\right)$和$\left(-u_{k},-v_{k}\right)$處些中心是根據頻率矩形的中心(M/2,N/2)確定的。對於每個濾波器，距離計算由下式執行：
$$
D_{k}(u, v)=\left[\left(u-M / 2-u_{k}\right)^{2}+\left(v-N / 2-v_{k}\right)^{2}\right]^{1 / 2}
$$
$$
D_{-k}(u, v)=\left[\left(u-M / 2+u_{k}\right)^{2}+\left(v-N / 2+v_{k}\right)^{2}\right]^{1 / 2}
$$
關於摩爾紋，簡單來說是差拍原理的一種表現。從數學上講，兩個頻率接近的等幅正弦波疊加，合成信號的幅度將按照兩個頻率之差變化。如果感光元件CCD（CMOS）像素的空間頻率與影像中條紋的空間頻率接近，就會產生摩爾紋。我們手機對着電腦拍照的時候看到 的波紋即是摩爾紋，現在手機都自帶摩爾紋消除了，所以這里無法用手機得到帶摩爾紋的圖片，這里給出數字圖像處理的摩爾紋圖片鏈接，帶摩爾紋圖片下載的鏈接為：http://www.imageprocessingplace.com/DIP-3E/dip3e_book_images_downloads.htm第四章的Fig0464(a)(car_75DPI_Moire).tif
下面通過OpenCV來完成摩爾紋消除實驗
示例：
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>
using namespace cv;
using namespace std;
typedef cv::Mat Mat;
Mat image_add_border(Mat &src)
{
int w = 2 * src.cols;
int h = 2 * src.rows;
std::cout << "src: " << src.cols << "*" << src.rows << std::endl;
cv::Mat padded;
copyMakeBorder(src, padded, 0, h - src.rows, 0, w - src.cols,
cv::BORDER_CONSTANT, cv::Scalar::all(0));
padded.convertTo(padded, CV_32FC1);
std::cout << "opt: " << padded.cols << "*" << padded.rows << std::endl;
return padded;
}
//transform to center 中心化
void center_transform(cv::Mat &src)
{
for (int i = 0; i < src.rows; i++) {
float *p = src.ptr<float>(i);
for (int j = 0; j < src.cols; j++) {
p[j] = p[j] * pow(-1, i + j);
}
}
}
//對角線交換內容
void zero_to_center(cv::Mat &freq_plane)
{
// freq_plane = freq_plane(Rect(0, 0, freq_plane.cols & -2, freq_plane.rows & -2));
//這里為什么&上-2具體查看opencv文檔
//其實是為了把行和列變成偶數 -2的二進制是11111111.......10 最后一位是0
int cx = freq_plane.cols / 2; int cy = freq_plane.rows / 2;//以下的操作是移動圖像 (零頻移到中心)
cv::Mat part1_r(freq_plane, cv::Rect(0, 0, cx, cy));//元素坐標表示為(cx,cy)
cv::Mat part2_r(freq_plane, cv::Rect(cx, 0, cx, cy));
cv::Mat part3_r(freq_plane, cv::Rect(0, cy, cx, cy));
cv::Mat part4_r(freq_plane, cv::Rect(cx, cy, cx, cy));
cv::Mat tmp;
part1_r.copyTo(tmp);//左上與右下交換位置(實部)
part4_r.copyTo(part1_r);
tmp.copyTo(part4_r);
part2_r.copyTo(tmp);//右上與左下交換位置(實部)
part3_r.copyTo(part2_r);
tmp.copyTo(part3_r);
}
void show_spectrum(cv::Mat &complexI)
{
cv::Mat temp[] = { cv::Mat::zeros(complexI.size(),CV_32FC1),
cv::Mat::zeros(complexI.size(),CV_32FC1) };
//顯示頻譜圖
cv::split(complexI, temp);
cv::Mat aa;
cv::magnitude(temp[0], temp[1], aa);
// zero_to_center(aa);
cv::divide(aa, aa.cols*aa.rows, aa);
double min, max;
cv::Point min_loc, max_loc;
cv::minMaxLoc(aa, &min, &max,&min_loc, &max_loc);
std::cout << "min: " << min << " max: " << max << std::endl;
std::cout << "min_loc: " << min_loc << " max_loc: " << max_loc << std::endl;
cv::circle(aa, min_loc, 20, cv::Scalar::all(1), 3);
cv::circle(aa, max_loc, 20, cv::Scalar::all(1), 3);
cv::imshow("src_img_spectrum", aa);
}
//頻率域濾波
cv::Mat frequency_filter(cv::Mat &padded, cv::Mat &blur)
{
cv::Mat plane[] = { padded, cv::Mat::zeros(padded.size(), CV_32FC1) };
cv::Mat complexIm;
cv::merge(plane, 2, complexIm);
cv::dft(complexIm, complexIm);//fourior transform
show_spectrum(complexIm);
cv::multiply(complexIm, blur, complexIm);
cv::idft(complexIm, complexIm, CV_DXT_INVERSE);//idft
cv::Mat dst_plane[2];
cv::split(complexIm,dst_plane);
center_transform(dst_plane[0]);
// center_transform(dst_plane[1]);
cv::magnitude(dst_plane[0], dst_plane[1], dst_plane[0]);//求幅值(模)
// center_transform(dst_plane[0]); //center transform
return dst_plane[0];
}
//陷波濾波器
cv::Mat notch_kernel(cv::Mat &scr, std::vector<cv::Point> ¬ch_pot, float D0)
{
cv::Mat notch_pass(scr.size(), CV_32FC2);
int row_num = scr.rows;
int col_num = scr.cols;
int n = 4;
for (int i = 0; i < row_num; i++) {
float *p = notch_pass.ptr<float>(i);
for (int j = 0; j < col_num; j++) {
float h_nr = 1.0;
for (unsigned int k = 0; k < notch_pot.size(); k++) {
int u_k = notch_pot.at(k).y;
int v_k = notch_pot.at(k).x;
float dk = sqrt(pow((i - u_k), 2) +
pow((j - v_k), 2));
float d_k = sqrt(pow((i + u_k), 2) +
pow((j + v_k), 2));
float dst_dk = 1.0 / (1.0 + pow(D0 / dk, 2 * n));
float dst_d_k = 1.0 / (1.0 + pow(D0 / d_k, 2 * n));
h_nr = dst_dk * dst_d_k * h_nr;
// std::cout << "notch_pot: " << notch_pot.at(k) << std::endl;
}
p[2 * j] = h_nr;
p[2 * j + 1] = h_nr;
}
}
cv::Mat temp[] = { cv::Mat::zeros(scr.size(), CV_32FC1),
cv::Mat::zeros(scr.size(), CV_32FC1) };
cv::split(notch_pass, temp);
std::string name = "notch濾波器d0=" + std::to_string(D0);
cv::Mat show;
cv::normalize(temp[0], show, 1, 0, CV_MINMAX);
cv::imshow(name, show);
return notch_pass;
}
std::string name_win("Notch_filter");
cv::Rect g_rectangle;
bool g_bDrawingBox = false;//是否進行繪制
cv::RNG g_rng(12345);
std::vector<cv::Point>notch_point;
void on_MouseHandle(int event, int x, int y, int flags, void*param);
void DrawRectangle(cv::Mat& img, cv::Rect box);
int main()
{
Mat srcImage = cv::imread("Fig0464(a)(car_75DPI_Moire).tif", cv::IMREAD_GRAYSCALE);
if (srcImage.empty())
return -1;
imshow("src_img", srcImage);
Mat padded = image_add_border(srcImage);
center_transform(padded);
Mat plane[] = { padded, cv::Mat::zeros(padded.size(), CV_32FC1) };
Mat complexIm;
merge(plane, 2, complexIm);
cv::dft(complexIm, complexIm);//fourior transform
Mat temp[] = { cv::Mat::zeros(complexIm.size(),CV_32FC1),
cv::Mat::zeros(complexIm.size(),CV_32FC1) };
//顯示頻譜圖
split(complexIm, temp);
Mat aa;
magnitude(temp[0], temp[1], aa);
divide(aa, aa.cols*aa.rows, aa);
cv::namedWindow(name_win);
cv::imshow(name_win, aa);
cv::setMouseCallback(name_win, on_MouseHandle, (void*)&aa);
Mat tempImage = aa.clone();
int key_value = -1;
while (1) {
key_value = cv::waitKey(10);
if (key_value == 27) {//按下ESC鍵查看濾波后的效果
break;
}
}
if (!notch_point.empty()) {
for (unsigned int i = 0; i < notch_point.size(); i++) {
cv::circle(tempImage, notch_point.at(i), 3, cv::Scalar(1), 2);
std::cout << notch_point.at(i) << std::endl;
}
}
cv::imshow(name_win, tempImage);
Mat ker = notch_kernel(complexIm, notch_point, 30);
cv::multiply(complexIm, ker, complexIm);
split(complexIm, temp);
magnitude(temp[0], temp[1], aa);
divide(aa, aa.cols*aa.rows, aa);
imshow("aa", aa);
cv::idft(complexIm, complexIm, CV_DXT_INVERSE);//idft
cv::Mat dst_plane[2];
cv::split(complexIm, dst_plane);
center_transform(dst_plane[0]);
//center_transform(dst_plane[1]);
// cv::magnitude(dst_plane[0],dst_plane[1],dst_plane[0]); //求幅值(模)
cv::normalize(dst_plane[0], dst_plane[0], 1, 0, CV_MINMAX);
imshow("dest", dst_plane[0]);
cv::waitKey(0);
return 1;
}
void on_MouseHandle(int event, int x, int y, int falgs, void* param)
{
Mat& image = *(cv::Mat*)param;
switch (event) {//鼠標移動消息
case cv::EVENT_MOUSEMOVE: {
if (g_bDrawingBox) {//標識符為真，則記錄下長和寬到Rect型變量中
g_rectangle.width = x - g_rectangle.x;
g_rectangle.height = y - g_rectangle.y;
}
}
break;
case cv::EVENT_LBUTTONDOWN: {
g_bDrawingBox = true;
g_rectangle = cv::Rect(x, y, 0, 0);//記錄起點
std::cout << "start point( " << g_rectangle.x << "," << g_rectangle.y << ")" << std::endl;
}
break;
case cv::EVENT_LBUTTONUP: {
g_bDrawingBox = false;
bool w_less_0 = false, h_less_0 = false;
if (g_rectangle.width < 0) {//對寬高小於0的處理
g_rectangle.x += g_rectangle.width;
g_rectangle.width *= -1;
w_less_0 = true;
}
if (g_rectangle.height < 0) {
g_rectangle.y += g_rectangle.height;
g_rectangle.height *= -1;
h_less_0 = true;
}
std::cout << "end Rect[ " << g_rectangle.x << "," << g_rectangle.y << "," << g_rectangle.width << ","
<< g_rectangle.height << "]" << std::endl;
if (g_rectangle.height > 0 && g_rectangle.width > 0) {
Mat imageROI = image(g_rectangle).clone();
double min, max;
cv::Point min_loc, max_loc;
cv::minMaxLoc(imageROI, &min, &max, &min_loc, &max_loc);
cv::circle(imageROI, max_loc, 10, 1);
max_loc.x += g_rectangle.x;
max_loc.y += g_rectangle.y;
notch_point.push_back(max_loc);
cv::circle(image, max_loc, 10, 1);
// cv::imshow( "ROI", imageROI );
cv::imshow("src", image);
}
}
break;
}
}
運行程序后，用鼠標在頻譜圖上畫圈，將四個陷波點標記后，在命令框內按下ESC鍵即可看到陷波結果：
7、同態濾波
圖像的形成是由光源的入射和反射到成像系統中來的到物體的外觀特征的一個過程，假如我們用形如f(x,y)的二維函數來表示圖像，則f(x,y)可以由兩個分量來表示：（1）入射到被觀察場景的光源照射總量（2）場景中物體所反射光的總量。這兩個分量分別被稱為入射分量和反射分量，可以分別表示為i(x,y)和r(x,y)。兩個函數的乘積為f(x,y)
$$
f(x, y)=i(x, y) r(x, y)
$$
其中$0<i(x, y)<\infty ,0<r(x, y)<1$，由於兩個函數的乘積的傅里葉變換是不可分離的，我們不能輕易地使用上式分別對照明和反射的頻率成分進行操作，即
$$
\mathcal{F}(F(x, y)) \neq \mathcal{F}(i(x, y)) \mathcal{F}(r(x, y))
$$
但是，假設我們定義
$$
\begin{aligned} z(x, y) &=\ln F(x, y) \\ &=\ln i(x, y)+\ln r(x, y) \end{aligned}
$$
則有
$$
\begin{aligned} \mathcal{F}(z(x, y)) &=\mathcal{F}(\ln F(x, y)) \\ &=\mathcal{F}(\ln i(x, y))+\mathcal{F}(\ln r(x, y)) \end{aligned}
$$
或
$$
Z(\omega, \nu)=I(\omega, \nu)+R(\omega, \nu)
$$
其中Z,I,R分別是z,$\ln i(x, y)$和$\ln r(x, y)$的傅里葉變換。函數Z表示低頻照明圖像和高頻反射圖像兩個圖像之和的傅立葉變換,其傳遞函數如下所示：
如果我們現在應用具有抑制低頻分量和增強高頻分量的傳遞函數的濾波器，那么我們可以抑制照明分量並增強反射分量。從而有
$$
\begin{aligned} S(\omega, \nu) &=H(\omega, \nu) Z(\omega, \nu) \\ &=H(\omega, \nu) I(\omega, \nu)+H(\omega, \nu) R(\omega, \nu) \end{aligned}
$$
其中S是傅里葉變換的結果。在空間領域
$$
\begin{aligned} s(x, y) &=\mathcal{F}^{-1}(S(\omega, \nu)) \\ &=\mathcal{F}^{-1}(H(\omega, \nu) I(\omega, \nu))+\mathcal{F}^{-1}(H(\omega, \nu) R(\omega, \nu)) \end{aligned}
$$
由定義：
$$
i^{\prime}(x, y)=\mathcal{F}^{-1}(H(\omega, \nu) I(\omega, \nu))
$$
$$
r^{\prime}(x, y)=\mathcal{F}^{-1}(H(\omega, \nu) R(\omega, \nu))
$$
我們有
$$
s(x, y)=i^{\prime}(x, y)+r^{\prime}(x, y)
$$
最后，因為z(x,y)是通過去輸入圖像的對數形成的，我們可以通過取濾波后的而結果的指數這一反處理來形成輸出圖像：
$$
\begin{aligned} \hat{F}(x, y) &=\exp [s(x, y)] \\ &=\exp \left[i^{\prime}(x, y)\right] \exp \left[r^{\prime}(x, y)\right] \\ &=i_{0}(x, y) r_{0}(x, y) \end{aligned}
$$
以上便是同態濾波的整個流程，可以總結如下：
圖像的照射分量通常由慢的空間變化來表征，而反射分量往往引起突變，特別是在不同物體的連接部分。這些特性導致圖像取對數后的傅里葉變換的低頻成分與照射分量相聯系，而高頻成分與反射分量相聯系。雖然這些聯系只是粗略的近似，但它們用在圖像濾波中是有益的。
使用同態濾波器可更好的控制照射成分和反射成分。這種控制需要指定一個濾波器H(u，v)，它可用不同的可控方法影響傅里葉變換的低頻和高頻成分。例如，采用形式稍微變化一下的高斯高通濾波器可得到函數：
$$
H(u, v)=\left(\gamma_{H}-\gamma_{L}\right)\left[1-e^{-c\left[D^{2}(u, v) / D_{0}^{2}\right]}\right]+\gamma_{L}
$$
其中$D(u, v)=\sqrt{\left(u-\frac{P}{2}\right)^{2}+\left(v-\frac{Q}{2}\right)^{2}}$，常數c控制函數邊坡的銳利度。
利用OpenCV的相關實現如下所示：
示例：
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>
using namespace cv;
using namespace std;
Mat Butterworth_Homomorphic_Filter(Size sz, int D, int n, float high_h_v_TB, float low_h_v_TB, Mat& realIm)
{
Mat single(sz.height, sz.width, CV_32F);
cout << sz.width << " " << sz.height << endl;
Point centre = Point(sz.height / 2, sz.width / 2);
double radius;
float upper = (high_h_v_TB * 0.1);
float lower = (low_h_v_TB * 0.1);
long dpow = D * D;
float W = (upper - lower);
for (int i = 0; i < sz.height; i++)
{
for (int j = 0; j < sz.width; j++)
{
radius = pow((float)(i - centre.x), 2) + pow((float)(j - centre.y), 2);
float r = exp(-n * radius / dpow);
if (radius < 0)
single.at<float>(i, j) = upper;
else
single.at<float>(i, j) = W * (1 - r) + lower;
}
}
single.copyTo(realIm);
Mat butterworth_complex;
//make two channels to match complex
Mat butterworth_channels[] = { Mat_<float>(single), Mat::zeros(sz, CV_32F) };
merge(butterworth_channels, 2, butterworth_complex);
return butterworth_complex;
}
//DFT 返回功率譜Power
Mat Fourier_Transform(Mat frame_bw, Mat& image_complex, Mat &image_phase, Mat &image_mag)
{
Mat frame_log;
frame_bw.convertTo(frame_log, CV_32F);
frame_log = frame_log / 255;
/*Take the natural log of the input (compute log(1 + Mag)*/
frame_log += 1;
log(frame_log, frame_log); // log(1 + Mag)
/*2. Expand the image to an optimal size
The performance of the DFT depends of the image size. It tends to be the fastest for image sizes that are multiple of 2, 3 or 5.
We can use the copyMakeBorder() function to expand the borders of an image.*/
Mat padded;
int M = getOptimalDFTSize(frame_log.rows);
int N = getOptimalDFTSize(frame_log.cols);
copyMakeBorder(frame_log, padded, 0, M - frame_log.rows, 0, N - frame_log.cols, BORDER_CONSTANT, Scalar::all(0));
/*Make place for both the complex and real values
The result of the DFT is a complex. Then the result is 2 images (Imaginary + Real), and the frequency domains range is much larger than the spatial one. Therefore we need to store in float !
That's why we will convert our input image "padded" to float and expand it to another channel to hold the complex values.
Planes is an arrow of 2 matrix (planes[0] = Real part, planes[1] = Imaginary part)*/
Mat image_planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
/*Creates one multichannel array out of several single-channel ones.*/
merge(image_planes, 2, image_complex);
/*Make the DFT
The result of thee DFT is a complex image : "image_complex"*/
dft(image_complex, image_complex);
/***************************/
//Create spectrum magnitude//
/***************************/
/*Transform the real and complex values to magnitude
NB: We separe Real part to Imaginary part*/
split(image_complex, image_planes);
//Starting with this part we have the real part of the image in planes[0] and the imaginary in planes[1]
phase(image_planes[0], image_planes[1], image_phase);
magnitude(image_planes[0], image_planes[1], image_mag);
//Power
pow(image_planes[0], 2, image_planes[0]);
pow(image_planes[1], 2, image_planes[1]);
Mat Power = image_planes[0] + image_planes[1];
return Power;
}
void Inv_Fourier_Transform(Mat input, Mat& inverseTransform)
{
/*Make the IDFT*/
Mat result;
idft(input, result, DFT_SCALE);
/*Take the exponential*/
exp(result, result);
vector<Mat> planes;
split(result, planes);
magnitude(planes[0], planes[1], planes[0]);
planes[0] = planes[0] - 1.0;
normalize(planes[0], planes[0], 0, 255, CV_MINMAX);
planes[0].convertTo(inverseTransform, CV_8U);
}
//SHIFT
void Shifting_DFT(Mat &fImage)
{
//For visualization purposes we may also rearrange the quadrants of the result, so that the origin (0,0), corresponds to the image center.
Mat tmp, q0, q1, q2, q3;
/*First crop the image, if it has an odd number of rows or columns.
Operator & bit to bit by -2 (two's complement : -2 = 111111111....10) to eliminate the first bit 2^0 (In case of odd number on row or col, we take the even number in below)*/
fImage = fImage(Rect(0, 0, fImage.cols & -2, fImage.rows & -2));
int cx = fImage.cols / 2;
int cy = fImage.rows / 2;
/*Rearrange the quadrants of Fourier image so that the origin is at the image center*/
q0 = fImage(Rect(0, 0, cx, cy));
q1 = fImage(Rect(cx, 0, cx, cy));
q2 = fImage(Rect(0, cy, cx, cy));
q3 = fImage(Rect(cx, cy, cx, cy));
/*We reverse each quadrant of the frame with its other quadrant diagonally opposite*/
/*We reverse q0 and q3*/
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
/*We reverse q1 and q2*/
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
}
void homomorphicFiltering(Mat src, Mat& dst)
{
Mat img;
Mat imgHls;
vector<Mat> vHls;
if (src.channels() == 3)
{
cvtColor(src, imgHls, CV_BGR2HSV);
split(imgHls, vHls);
vHls[2].copyTo(img);
}
else
src.copyTo(img);
//DFT
//cout<<"DFT "<<endl;
Mat img_complex, img_mag, img_phase;
Mat fpower = Fourier_Transform(img, img_complex, img_phase, img_mag);
Shifting_DFT(img_complex);
Shifting_DFT(fpower);
//int D0 = getRadius(fpower,0.15);
int D0 = 10;
int n = 2;
int w = img_complex.cols;
int h = img_complex.rows;
//BHPF
// Mat filter,filter_complex;
// filter = BHPF(h,w,D0,n);
// Mat planes[] = {Mat_<float>(filter), Mat::zeros(filter.size(), CV_32F)};
// merge(planes,2,filter_complex);
int rH = 150;
int rL = 50;
Mat filter, filter_complex;
filter_complex = Butterworth_Homomorphic_Filter(Size(w, h), D0, n, rH, rL, filter);
//do mulSpectrums on the full dft
mulSpectrums(img_complex, filter_complex, filter_complex, 0);
Shifting_DFT(filter_complex);
//IDFT
Mat result;
Inv_Fourier_Transform(filter_complex, result);
if (src.channels() == 3)
{
vHls.at(2) = result(Rect(0, 0, src.cols, src.rows));
merge(vHls, imgHls);
cvtColor(imgHls, dst, CV_HSV2BGR);
}
else
result.copyTo(dst);
}
int main()
{
Mat img = imread("test.jpg");
imshow("rogin img", img);
Mat dst;
homomorphicFiltering(img, dst);
imshow("homoMorphic filter", (img+dst)/2);
waitKey();
return 0;
}
運行結果如下所示：
