图像处理之频率域应用案例(C++)
前言
本文主要介绍几个比较典型的图像从空间域转换到频率域的案例。
一、原理
频率特征是图像的灰度变化特征。一个重要的经验结论:低频代表图像整体轮廓,高频代表了图像噪声,中频代表图像边缘、纹理等细节。因此,我们可以设计不同的滤波器对频率域图像进行滤波,从而实现我们想获得的频率特征。
注意:
1.在做傅里叶变换(FFT)时,通常只使用灰度图像数据,所以需要将彩色图像数据转换为灰度图像数据。同时灰度图像数据作为傅里叶变换复数单元的实数部分,虚数部分直接将值设置为0即可。
2.变换到频域后的数据,低频中心(直流分量)分散在uv坐标系图像的四个角落上,为了方便后续频域处理,通常需要将低频中心位置偏移到uv坐标系图像的中心,需要做一个位置变换处理,然后进行频域滤波,滤波完成再将低频中心重新偏移到四个角落上做IFFT变换。
二、代码实现
1.案例介绍
在复杂背景却又有一些纹理背景中,想要检测出缺陷(位于中心处),在空间域可能比较难以直接检测,可以转换到频率域,设计相应的滤波器,最后再转换到空间域进行后处理。
2.代码实现
#include <iostream>
#include <opencv.hpp>
using namespace std;
class Filter
{
public:
Filter(){ }
virtual cv::Mat getFilter()=0; //获得不同的滤波器矩阵
virtual ~Filter() { }
};
class idealLowPassFilter :public Filter //理想低通滤波器
{
int h, w;
double D0;
public:
idealLowPassFilter(int H, int W, double d0) : h(H), w(W), D0(d0)
{
}
//h,w--为图像的长宽、D0为带通半径
cv::Mat getFilter()
{
//1. 计算频谱中心
int cx = w / 2;
int cy = h / 2;
//2. 构建滤波器
cv::Mat filter(cv::Size(w, h), CV_32FC1);
for(int i=0;i<h;i++)
for (int j = 0; j < w; j++)
{
double D = pow(i - cy, 2) + pow(j - cx, 2);
if (D > D0 * D0)
filter.at<float>(i, j) = 0;
else
filter.at<float>(i, j) = 1;
}
//3. 中心化
cv::Mat temp = filter.clone();
temp(cv::Rect(0, 0, cx, cy)).copyTo(filter(cv::Rect(cx, cy, cx, cy)));//左上到右下
temp(cv::Rect(0, cy, cx, cy)).copyTo(filter(cv::Rect(cx, 0, cx, cy)));//左下到右上
temp(cv::Rect(cx, 0, cx, cy)).copyTo(filter(cv::Rect(0, cy, cx, cy)));//右上到左下
temp(cv::Rect(cx, cy, cx, cy)).copyTo(filter(cv::Rect(0, 0, cx, cy)));//右下到左上
return filter;
}
};
class idealHighPassFilter :public Filter //理想高通滤波器
{
int h, w;
double D0;
public:
idealHighPassFilter(int H, int W, double d0) : h(H), w(W), D0(d0)
{
}
//h,w--为图像的长宽、D0为带通半径
cv::Mat getFilter()
{
//1. 计算频谱中心
int cx = w / 2;
int cy = h / 2;
//2. 构建滤波器
cv::Mat filter(cv::Size(w, h), CV_32FC1);
for (int i = 0; i < h; i++)
for (int j = 0; j < w; j++)
{
double D = pow(i - cy, 2) + pow(j - cx, 2);
if (D > D0 * D0)
filter.at<float>(i, j) = 1;
else
filter.at<float>(i, j) = 0;
}
//3. 中心化
cv::Mat temp = filter.clone();
temp(cv::Rect(0, 0, cx, cy)).copyTo(filter(cv::Rect(cx, cy, cx, cy)));//左上到右下
temp(cv::Rect(0, cy, cx, cy)).copyTo(filter(cv::Rect(cx, 0, cx, cy)));//左下到右上
temp(cv::Rect(cx, 0, cx, cy)).copyTo(filter(cv::Rect(0, cy, cx, cy)));//右上到左下
temp(cv::Rect(cx, cy, cx, cy)).copyTo(filter(cv::Rect(0, 0, cx, cy)));//右下到左上
return filter;
}
};
class GaussLowPassFilter :public Filter //高斯低通滤波器
{
int h, w;
double D0;
public:
GaussLowPassFilter(int H, int W, double d0) :h(H), w(W), D0(d0)
{
}
//h,w--为图像的长宽、D0为带通半径
cv::Mat getFilter()
{
//1. 计算频谱中心
int cx = w / 2;
int cy = h / 2;
//2. 构建滤波器
cv::Mat filter(cv::Size(w, h), CV_32FC1);
for (int i = 0; i < h; i++)
for (int j = 0; j < w; j++)
{
double D = pow(i - cy, 2) + pow(j - cx, 2);
filter.at<float>(i, j) = exp(-1 * D / (2 * D0 * D0));
}
//3. 中心化
cv::Mat temp = filter.clone();
temp(cv::Rect(0, 0, cx, cy)).copyTo(filter(cv::Rect(cx, cy, cx, cy)));//左上到右下
temp(cv::Rect(0, cy, cx, cy)).copyTo(filter(cv::Rect(cx, 0, cx, cy)));//左下到右上
temp(cv::Rect(cx, 0, cx, cy)).copyTo(filter(cv::Rect(0, cy, cx, cy)));//右上到左下
temp(cv::Rect(cx, cy, cx, cy)).copyTo(filter(cv::Rect(0, 0, cx, cy)));//右下到左上
return filter;
}
};
class GaussHighPassFilter :public Filter //高斯高通滤波器
{
int h, w;
double D0;
public:
GaussHighPassFilter(int H, int W, double d0) :h(H), w(W), D0(d0)
{
}
//h,w--为图像的长宽、D0为带通半径
cv::Mat getFilter()
{
//1. 计算频谱中心
int cx = w / 2;
int cy = h / 2;
//2. 构建滤波器
cv::Mat filter(cv::Size(w, h), CV_32FC1);
for (int i = 0; i < h; i++)
for (int j = 0; j < w; j++)
{
double D = pow(i - cy, 2) + pow(j - cx, 2);
filter.at<float>(i, j) = 1-exp(-1*D/(2*D0*D0));
}
//3. 中心化
cv::Mat temp = filter.clone();
temp(cv::Rect(0, 0, cx, cy)).copyTo(filter(cv::Rect(cx, cy, cx, cy)));//左上到右下
temp(cv::Rect(0, cy, cx, cy)).copyTo(filter(cv::Rect(cx, 0, cx, cy)));//左下到右上
temp(cv::Rect(cx, 0, cx, cy)).copyTo(filter(cv::Rect(0, cy, cx, cy)));//右上到左下
temp(cv::Rect(cx, cy, cx, cy)).copyTo(filter(cv::Rect(0, 0, cx, cy)));//右下到左上
return filter;
}
};
class ButterworthLowPassFilter :public Filter //巴特沃斯低通滤波器
{
int h, w;
double D0; //截止频率
int n; //阶数
public:
ButterworthLowPassFilter(int H, int W, double d0, int N):h(H),w(W),D0(d0), n(N)
{
}
//h,w--为图像的长宽、D0为带通半径
cv::Mat getFilter()
{
//1. 计算频谱中心
int cx = w / 2;
int cy = h / 2;
//2. 构建滤波器
cv::Mat filter(cv::Size(w, h), CV_32FC1);
for (int i = 0; i < h; i++)
for (int j = 0; j < w; j++)
{
double D = pow(i - cy, 2) + pow(j - cx, 2);
filter.at<float>(i, j) = 1/(pow(D / (D0 * D0),n)+1);
}
//3. 中心化
cv::Mat temp = filter.clone();
temp(cv::Rect(0, 0, cx, cy)).copyTo(filter(cv::Rect(cx, cy, cx, cy)));//左上到右下
temp(cv::Rect(0, cy, cx, cy)).copyTo(filter(cv::Rect(cx, 0, cx, cy)));//左下到右上
temp(cv::Rect(cx, 0, cx, cy)).copyTo(filter(cv::Rect(0, cy, cx, cy)));//右上到左下
temp(cv::Rect(cx, cy, cx, cy)).copyTo(filter(cv::Rect(0, 0, cx, cy)));//右下到左上
return filter;
}
};
class ButterworthHighPassFilter :public Filter //巴特沃斯高通滤波器
{
int h, w;
double D0; //截止频率
int n; //阶数
public:
ButterworthHighPassFilter(int H, int W, double d0,int N) :h(H), w(W), D0(d0),n(N)
{
}
//h,w--为图像的长宽、D0为带通半径
cv::Mat getFilter()
{
//1. 计算频谱中心
int cx = w / 2;
int cy = h / 2;
//2. 构建滤波器
cv::Mat filter(cv::Size(w, h), CV_32FC1);
for (int i = 0; i < h; i++)
for (int j = 0; j < w; j++)
{
double D = pow(i - cy, 2) + pow(j - cx, 2);
filter.at<float>(i, j) = 1 / (pow((D0 * D0)/D, n) + 1);
}
//3. 中心化
cv::Mat temp = filter.clone();
temp(cv::Rect(0, 0, cx, cy)).copyTo(filter(cv::Rect(cx, cy, cx, cy)));//左上到右下
temp(cv::Rect(0, cy, cx, cy)).copyTo(filter(cv::Rect(cx, 0, cx, cy)));//左下到右上
temp(cv::Rect(cx, 0, cx, cy)).copyTo(filter(cv::Rect(0, cy, cx, cy)));//右上到左下
temp(cv::Rect(cx, cy, cx, cy)).copyTo(filter(cv::Rect(0, 0, cx, cy)));//右下到左上
return filter;
}
};
//自定义滤波器
class myPassFilter :public Filter //理想高通滤波器
{
int h, w;
double D0;
public:
myPassFilter(int H, int W, double d0) : h(H), w(W), D0(d0)
{
}
//h,w--为图像的长宽、D0为带通半径
cv::Mat getFilter()
{
//1. 计算频谱中心
int cx = w / 2;
int cy = h / 2;
//2. 构建滤波器
cv::Mat filter(cv::Size(w, h), CV_32FC1);
for (int i = 0; i < h; i++)
for (int j = 0; j < w; j++)
{
double D = pow(i - cy, 2) + pow(j - cx, 2);
if (D > D0 * D0)
filter.at<float>(i, j) = 1;
else
filter.at<float>(i, j) = 0.5;
}
//3. 中心化
cv::Mat temp = filter.clone();
temp(cv::Rect(0, 0, cx, cy)).copyTo(filter(cv::Rect(cx, cy, cx, cy)));//左上到右下
temp(cv::Rect(0, cy, cx, cy)).copyTo(filter(cv::Rect(cx, 0, cx, cy)));//左下到右上
temp(cv::Rect(cx, 0, cx, cy)).copyTo(filter(cv::Rect(0, cy, cx, cy)));//右上到左下
temp(cv::Rect(cx, cy, cx, cy)).copyTo(filter(cv::Rect(0, 0, cx, cy)));//右下到左上
return filter;
}
};
/*
* @param cv::Mat src 输入图片,单通道灰度图像
* @param cv::Mat dst 输出图片
* @brief 空间域转换到频率域,进行频率域的滤波,然后再返回空间域分析变化
*/
void Spatial2Frequency(const cv::Mat& src, cv::Mat& dst)
{
//1. 获得dft的最优尺寸
int h = cv::getOptimalDFTSize(src.rows);
int w = cv::getOptimalDFTSize(src.cols);
//2. 对图像进行填充
cv::Mat padImg;
cv::copyMakeBorder(src, padImg, 0, h - src.rows, 0, w - src.cols, cv::BORDER_CONSTANT, cv::Scalar(0));
//3. 创建矩阵存储实数和虚数
cv::Mat planes[] = { cv::Mat_<float>(padImg),cv::Mat::zeros(padImg.size(),CV_32FC1) };
cv::Mat dftImg;
cv::merge(planes, 2, dftImg);
//4. 执行傅里叶变换,转换到频率域
cv::dft(dftImg, dftImg);
//5. 构建不同的滤波器
double D0 = 10; //带宽
//要获得不同的滤波器,只需要修改此处的类即可。
//Filter* filter = new idealLowPassFilter(h,w,D0); //理想低通滤波器
//Filter* filter = new idealHighPassFilter(h,w,D0); //理想高通滤波器
//Filter* filter = new GaussLowPassFilter(h, w, D0); //高斯低通滤波器
//Filter* filter = new GaussHighPassFilter(h, w, D0); //高斯高通滤波器
//Filter* filter = new ButterworthLowPassFilter(h, w, D0,1); //巴特沃斯低通滤波器
//Filter* filter = new ButterworthHighPassFilter(h, w, D0, 1); //巴特沃斯高通滤波器
Filter* filter = new myPassFilter(h, w, D0); //自定义滤波器
cv::Mat PassFilter=filter->getFilter().clone();
//6. 将滤波器放到一个两通道的Mat矩阵中
cv::Mat filterPlanes[] = { cv::Mat_<float>(PassFilter),cv::Mat::zeros(PassFilter.size(),CV_32FC1) };
cv::merge(filterPlanes, 2,PassFilter);
//7. 相应的频谱进行滤波
cv::mulSpectrums(dftImg, PassFilter, dftImg,0);
//8. 反变换回空间域
cv::Mat spatial;
std::vector<cv::Mat> channels;
cv::idft(dftImg, dftImg,cv::DFT_SCALE);
cv::split(dftImg, channels);
cv::magnitude(channels[0], channels[1], spatial);
//9. 归一化
cv::normalize(spatial, spatial, 0, 255, cv::NORM_MINMAX, CV_8UC1);
spatial(cv::Rect(0, 0, src.cols, src.rows)).copyTo(dst);
//cv::imshow("dst", dst);
//释放指针
delete filter;
}
/*
* @param cv::Mat src 输入图片,单通道灰度图像
* @param cv::Mat dst 输出图片,三通道的结果图像
* @brief 后处理函数
*/
void postProcess(cv::Mat& src, cv::Mat& dst)
{
//通过观察图片,发现缺陷一半白,一半黑,采用阈值处理,然后膨胀再采用与操作,提取相同区域
cv::Mat light, dark;
cv::threshold(src, dark, 10, 255, cv::THRESH_BINARY_INV);
cv::threshold(src, light, 230, 255, cv::THRESH_BINARY);
cv::Mat kernel = cv::getStructuringElement(cv::MORPH_RECT, cv::Size(7, 7));
cv::dilate(dark, dark, kernel);
cv::dilate(light, light, kernel);
cv::bitwise_and(light, dark, src);
//cv::imshow("dst_1", src);
//绘制轮廓
std::vector<std::vector<cv::Point>> contours;
std::vector<cv::Vec4i> hierarchy;
cv::findContours(src, contours, hierarchy, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_SIMPLE);
for (int i = 0; i < contours.size(); i++)
{
cv::Rect rect = cv::boundingRect(contours[i]);
cv::rectangle(dst, cv::Rect(rect.x, rect.y, 10, 10), cv::Scalar(0, 0, 255));
}
}
int main()
{
//读取图片
string filepath = "F://work_study//algorithm_demo//hard_blob.bmp";
cv::Mat src = cv::imread(filepath, cv::IMREAD_GRAYSCALE);
if (src.empty())
{
std::cout << "imread error" << std::endl;
return -1;
}
cv::Mat dst(src.size(), src.type());
//空间域转换到频率域,设计滤波器进行处理,提取想要的特征
Spatial2Frequency(src, dst);
//可视化结果,将单通道灰度图转换为三通道彩色图像
cv::cvtColor(src, src, cv::COLOR_GRAY2BGR);
//后处理
postProcess(dst, src);
//显示图片
imwrite("dst.bmp", dst);
cv::waitKey(0);
return 0;
}
3.结果展示
总结
本文介绍了图像在频率域中的典型应用,并且对滤波器类进行了重构,方便后续读者根据自己的需要进行继承设计自己的滤波器,代码中已经实现了理想高/低通滤波器、高斯高/低通滤波器、巴特沃斯高/低通滤波器。
因为笔者水平有限,有错误欢迎指出,代码本人均在本地运行实验正确,大家放心使用。