1六、頻率域濾波 頻率域濾波基礎——傅里葉變換計算及應用基礎 基於OpenCV的同態濾波
一、頻率域與空間域之間的關係
在頻率域濾波基礎——傅里葉變換計算及應用基礎中，咱們已經知道了如何將圖像轉換到頻率域，以及如何將頻率域圖像經過傅里葉逆變換轉換回圖像，這樣一來，能夠利用空域圖像與頻譜之間的對應關係，嘗試將空域卷積濾波變換爲頻域濾波，然後再將頻域濾波處理後的圖像反變換回空間域，從而達到圖像加強的目的，這樣作的一個最主要的吸引力再域頻率域濾波具備直觀性的特色。html
根據著名的卷積定律，兩個二維連續函數再空間域的卷積可由其相應的兩個傅里葉變換乘積的反變換而獲得，反之，在頻率域中的卷積可由在空間域中乘積的傅里葉變換而獲得。即ios
$$
\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}
$$算法
二、頻率域濾波的基本步驟
本文重點就來說一講如何經過頻率域濾波來實現圖像的平滑和銳化。首先將頻率域濾波的步驟來作一個小結以下：ubuntu
- 給定一幅大小爲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.
$$less
其中，爲忽略因爲計算不許確致使的寄生覆成分，選擇實部，下標p指出咱們處理的是填充後的圖像。dom
- 從$g_{p}(x, y)$的左上象限提取MXN區域，獲得最終處理結果$g(x, y)$
由上述敘述可知，濾波的關鍵取決於濾波函數$H(u, v)$，經常稱爲濾波器，或濾波器傳遞函數，由於它在濾波中抑制或除去了頻譜中的某些份量，而保留其餘的一些頻率不受影響，從而達到濾波的目的。而本節將講解一些經常使用的濾波器。ide
三、使用頻率域濾波平滑圖像
在空間域咱們已經講過圖像平滑的目的及空間域平滑濾波，如今在頻率域濾波中咱們首先講解三種頻率平滑濾波器：理想濾波器、巴特沃斯濾波器、高斯濾波器。這三種濾波器涵蓋了從很是急劇（理想）到很是平滑的濾波範圍。函數
一、理想低通濾波器（ILPF)
首先最容易想到得衰減高頻成分得方法就是在一個稱爲截止頻率得位置截斷全部的高頻成分，將圖像頻譜中全部高於這一截止頻率的頻譜成分設爲0，低於截止頻率的成分保持不變。如圖所示（左圖是其透視圖，右圖是其圖像表示）：post
可以達到這種用效果的濾波器咱們稱之爲理想濾波器，用公式表示爲：優化
$$
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;
}
從左至右依次是：原圖，優化事後的圖、傅里葉變換圖、理想低通濾波器處理事後的高斯變換圖、理想低通濾波器處理事後的效果圖：
二、巴特沃思低通濾波器(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;
}
從左至右依次是：原圖，優化事後的圖、傅里葉變換圖、巴特沃斯濾波器平滑後的高斯變換圖、巴特沃斯濾波器平滑後的效果圖：
三、高斯低通濾波器(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;
}
從左至右依次是：原圖，優化事後的圖、傅里葉變換圖、高斯低通濾波後的高斯變換圖、高斯平滑事後的效果圖：
四、使用頻率域濾波銳化圖像
前面已經講過了經過衰減圖像的傅里葉變換的高頻成分能夠平滑圖像。由於邊緣和其餘灰度的急劇變化與高頻成分有關，因此圖像的銳化可在頻率域經過濾波來實現，高通濾波器會衰減傅里葉變換中的低頻成分而不擾亂高頻信息。故咱們能夠考慮用高通濾波器來進行圖像的銳化,一個高頻濾波器是從給定低通濾波器用下式獲得的:
$$
H_{\mathrm{HP}}(u, v)=1-H_{\mathrm{LP}}(u, v)
$$
與低通濾波器對應，這裏介紹理想高通濾波器，巴特沃斯高通濾波器，高斯高通濾波器三種高通濾波器。
一、理想高通濾波器 (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;
}
結果以下：
二、巴特沃思高通濾波器 （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;
}
結果以下：
三、高斯高通濾波器（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.一、帶阻濾波器、帶通濾波器
所謂的帶通濾波器便是將低頻濾波和高頻濾波相結合，最後咱們能夠利用帶阻濾波器濾掉週期性出現的噪聲，由前面可知，
理想帶通濾波器公式爲
$$
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.二、陷波濾波器
陷波器是一個相對於帶通濾波更有用的選擇濾波器，其拒絕事先定義好的一個關於頻率矩陣中心的一個領域的頻域。零相移濾波器必須是關於遠點對稱的，所以，一箇中心位於$\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鍵便可看到陷波結果：
七、同態濾波
圖像的造成是由光源的入射和反射到成像系統中來的到物體的外觀特徵的一個過程，假如咱們用形如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;
}
運行結果以下所示：
OpenCV 頻率域處理：離散傅里葉變換、頻率域濾波
CS425 Lab: Frequency Domain Processing
基於opencv的理想低通濾波器和巴特沃斯低通濾波器
opencv 頻率域濾波實例
OpenCV 頻率域實現鈍化模板、高提高濾波和高頻強調濾波 C++
【數字圖像處理】4.10:灰度圖像-頻域濾波 概論
【數字圖像處理】4.11:灰度圖像-頻域濾波 濾波器
【數字圖像處理】4.12:灰度圖像-頻域濾波 同態濾波
基於OpenCV的同態濾波
計算機視覺（六）：頻率域濾波器
OpenCV使用陷波濾波器減小摩爾波紋 C++
相關文章
- 1. 頻率域濾波
- 2. 頻域濾波(1)--基礎與低通濾波器
- 3. 頻率域濾波（2）
- 4. PIE SDK頻率域濾波
- 5. Matlab中常用的空間域濾波和頻率域濾波
- 6. 頻域濾波
- 7. 空間域濾波與頻率域濾波比較
- 8. 空間濾波&頻域濾波(2)
- 9. 空間濾波&頻域濾波(3)
- 10. [數字圖像處理]頻域濾波(1)--基礎與低通濾波器
- 更多相關文章...
- • R 基礎運算 - R 語言教程
- • Rust 基礎語法 - RUST 教程
- • 使用Rxjava計算圓周率
- • ☆基於Java Instrument的Agent實現
相關標籤/搜索
每日一句
-
每一个你不满意的现在，都有一个你没有努力的曾经。
歡迎關注本站公眾號,獲取更多信息
相關文章
- 1. 頻率域濾波
- 2. 頻域濾波(1)--基礎與低通濾波器
- 3. 頻率域濾波（2）
- 4. PIE SDK頻率域濾波
- 5. Matlab中常用的空間域濾波和頻率域濾波
- 6. 頻域濾波
- 7. 空間域濾波與頻率域濾波比較
- 8. 空間濾波&頻域濾波(2)
- 9. 空間濾波&頻域濾波(3)
- 10. [數字圖像處理]頻域濾波(1)--基礎與低通濾波器