写在前面： 刚开始接触数字图像处理频率域滤波时，很是陌生，感觉这个技术使用范围很窄，不如空域直接处理来的实在，最近看书发现有些情况又不得不在频率域中进行操作，个人感觉图像的复原与重建就是最大的应用点。特此实现一些基本的频率域滤波操作为后学习打下基础…
前处理： 包括对图像边界填充，使之达到OpenCV傅里叶变换最佳尺寸，然后就是将
f
(
x
,
y
)
f(x,y)
f(x,y)乘以
−
1
(
x
+
y
)
-1^{(x+y)}
−1(x+y)，使傅里叶变换位于填充后图像大小的频率矩形的中心。
后处理： 包括从反变换回来的图的左上象限提取一个大小为MXN(原图尺寸)的区域，得到与输入图像大小相同的滤波后的结果 g ( x , y ) g(x,y) g(x,y)。
具体操作步骤详见《数字图像处理》第4章。
//***************************************tool***************************************//
//function: 使傅里叶频谱中心化
void fftshift(const cv::Mat& src, cv::Mat& dst)
{
dst = src.clone();
if (dst.type() != CV_32F) dst.convertTo(dst, CV_32F);
for (int r = 0; r < dst.rows; ++r)
{
for (int c = 0; c < dst.cols; ++c)
{
if ((r + c) % 2 != 0) //r + c 不为偶数，当前点变为负数
dst.at<float>(r, c) = -dst.at<float>(r, c);
}
}
}
//function: 创建最佳傅里叶变换尺寸
cv::Size Make_BetterSize(const cv::Mat src)
{
int w = cv::getOptimalDFTSize(src.cols);
int h = cv::getOptimalDFTSize(src.rows);
//width, height
return Size(w, h);
}
//function: 将图像中负数截断为0
//opencv中应该有现成的工具函数，需要查一下文档
void MinusToZero(const cv::Mat& src, cv::Mat& dst)
{
dst = src.clone();
for (int r = 0; r < dst.rows; ++r)
{
for (int c = 0; c < dst.cols; ++c)
{
if (dst.at<float>(r, c) < 0)
dst.at<float>(r, c) = 0.0;
}
}
}
//*****************************calculate filtering kernel*****************************//
//function: 创建理想低通滤波器
void calcILPF(const cv::Mat& src, int D0, Mat& dst)
{
Size filterSize = Make_BetterSize(src);
Mat H(filterSize, CV_32FC1, Scalar::all(0));
int cy = H.rows / 2;
int cx = H.cols / 2;
for (int y = 0; y < H.rows; ++y)
{
float* H_data = H.ptr<float>(y);
for (int x = 0; x < H.cols; ++x)
{
double d = sqrt(pow(y - cy, 2) + pow(x - cx, 2));
//当距离小于截止频率时，设置为1
if (d <= (double)D0)
H_data[x] = 1.f;
}
}
//改为双通道图像便于与F(u,v)相乘
Mat planes[] = { H.clone(), H.clone() };
Mat complexH;
merge(planes, 2, complexH);
complexH.copyTo(dst);
}
//function: 巴特沃斯低通滤波器
void calcButterWorth(const cv::Mat& src, int D0, Mat& dst, float n)
{
Size filterSize = Make_BetterSize(src);
cv::Mat H = cv::Mat::zeros(filterSize, CV_32FC1);
int cy = H.rows / 2;
int cx = H.cols / 2;
for (int y = 0; y < H.rows; ++y)
{
float* H_data = H.ptr<float>(y); //指向第y行首个元素的指针
for (int x = 0; x < H.cols; ++x)
{
double d = sqrt(pow(y - cy, 2) + pow(x - cx, 2));
//根据巴特沃斯低通滤波器计算H(u, v)
H_data[x] = 1.0 / (1 + pow(d / D0, 2 * n));
}
}
Mat planes[] = { H.clone(), H.clone() };
Mat complexH;
merge(planes, 2, complexH);
complexH.copyTo(dst);
}
//function: 创建高斯低通滤波器
void calcGLPF(const cv::Mat& src, int D0, Mat& dst)
{
Size filterSize = Make_BetterSize(src);
cv::Mat H = cv::Mat::zeros(filterSize, CV_32FC1);
int cy = H.rows / 2;
int cx = H.cols / 2;
for (int y = 0; y < H.rows; ++y)
{
float* H_data = H.ptr<float>(y); //指向第y行首个元素的指针
for (int x = 0; x < H.cols; ++x)
{
double d = sqrt(pow(y - cy, 2) + pow(x - cx, 2));
//根据gauss低通滤波器计算H(u, v)
H_data[x] = expf((-d * d) / (2 * D0 * D0));
}
}
Mat planes[] = { H.clone(), H.clone() };
Mat complexH;
merge(planes, 2, complexH);
complexH.copyTo(dst);
}
//function: 创建高斯高通滤波器
//这里未做简化，只需要用1减去高斯低通核即可创建高通核
void calcGHPF(const cv::Mat& src, int D0, Mat& dst)
{
Size filterSize = Make_BetterSize(src);
cv::Mat H = cv::Mat::ones(filterSize, CV_32FC1);
int cy = H.rows / 2;
int cx = H.cols / 2;
for (int y = 0; y < H.rows; ++y)
{
float* H_data = H.ptr<float>(y); //指向第y行首个元素的指针
for (int x = 0; x < H.cols; ++x)
{
double d = sqrt(pow(y - cy, 2) + pow(x - cx, 2));
//根据gauss低通滤波器计算H(u, v)
H_data[x] = 1.f - expf((-d * d) / (2 * D0 * D0));
}
}
Mat planes[] = { H.clone(), H.clone() };
Mat complexH;
merge(planes, 2, complexH);
complexH.copyTo(dst);
}
//********************************执行G(u,v) = F(u,v) * H(u,v)*******************************//
//function: 将输入图像在频率域进行处理，然后输出
//InputImage: 输入单通道图像
//H: 频率域滤波核
void frequency_filter(const cv::Mat& InputImage, const cv::Mat& H, cv::Mat& OutputImage)
{
cv::Mat padded;
int m = getOptimalDFTSize(InputImage.rows);
int n = getOptimalDFTSize(InputImage.cols);
//on the border add zero values
cv::copyMakeBorder(InputImage, padded, 0, m - InputImage.rows, 0, n - InputImage.cols, BORDER_CONSTANT, Scalar::all(0));
cv::Mat planes[] = { Mat_<float>(padded), cv::Mat::zeros(padded.size(), CV_32F) };
//使傅里叶变换中心化
fftshift(planes[0], planes[0]);
//创建一个复数图像F(u, v)
cv::Mat complexF;
cv::merge(planes, 2, complexF);
//执行傅里叶变换
cv::dft(complexF, complexF);
cv::Mat complexFH;
cv::Mat ifft;
//采用对应像素相乘G(u, v) = H(u, v) * F(u, v)
cv::multiply(complexF, H, complexFH);
cv::idft(complexFH, ifft, DFT_REAL_OUTPUT);
//最后不要忘了再完成一次移动
fftshift(ifft, ifft);
//为了显示效果将负数截断为0
MinusToZero(ifft, ifft);
normalize(ifft, ifft, 0, 1, NORM_MINMAX);
ifft.convertTo(ifft, CV_8UC1, 255);
//裁剪回原尺寸输出
ifft(cv::Rect(0, 0, InputImage.cols, InputImage.rows)).copyTo(OutputImage);
}
//*********************************main()*******************************//
int main()
{
string path = "blurring_test.tif";
Mat SrcImage = imread(path, IMREAD_GRAYSCALE);
if (!SrcImage.data) {
std::cout << "Could not open or find the image" << std::endl;
return -1;
}
//高通滤波结果和高通核
Mat H_outimg, H_kernel_H;
//低通滤波结果和低通核
Mat L_outimg, L_kernel_H;
calcGLPF(SrcImage, 60, L_kernel_H);
frequency_filter(SrcImage, L_kernel_H, L_outimg);
calcGHPF(SrcImage, 80, H_kernel_H);
frequency_filter(SrcImage, H_kernel_H, H_outimg);
imwrite("low_blurring_test.png", L_outimg);
imwrite("high_blurring_test.png", H_outimg);
cv::waitKey(0);
return 0;
}
原图：
高斯低通滤波：
高斯高通滤波：
