观察下面图像：
有明显的近似水平线的噪音模式，我们希望他在频率域内的成分沿垂直轴集中；
下图是该图像的谱；
将垂直轴的分量去掉，构建的陷波器如下：
原图与结果图对比：
滤掉的空间噪音模式：
代码实现：
#include "opencv2/opencv.hpp"
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::Rect> ¬ch_Rect )
{
// cv::Mat notch_pass = cv::Mat::zeros(scr.size(),CV_32FC2);
cv::Mat notch_pass = cv::Mat::ones(scr.size(),CV_32FC2);
int row_num = scr.rows;
int col_num = scr.cols;
for( unsigned int k = 0; k < notch_Rect.size(); k++ ){
for(int i=0; i<row_num; i++ ){
float *p = notch_pass.ptr<float>(i);
for(int j=0; j<col_num; j++ ){
if( i > notch_Rect.at(k).y
&& i < notch_Rect.at(k).y + notch_Rect.at(k).height
&& j > notch_Rect.at(k).x
&& j < notch_Rect.at(k).x + notch_Rect.at(k).width ){
p[2*j] = 0; //=1 与上面的zeros 共同获取空间噪音模式
p[2*j+1] = 0; //=1
}
}
}
}
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滤波器";
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::Rect> notch_rect;
void on_MouseHandle(int event, int x, int y, int flags, void*param);
void DrawRectangle(cv::Mat& img, cv::Rect box);
int main(int argc, char * argv[])
{
if( argc != 2){
std::cerr << "Usage: " << argv[0] << " <img_name> " << std::endl;
return -1;
}
Mat srcImage = cv::imread(argv[1], cv::IMREAD_GRAYSCALE);
if( srcImage.empty() )
return -1;
cv::resize( srcImage, srcImage, cv::Size(), 0.5, 0.5);
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 key,break
break;
}
}
cv::imshow(name_win, tempImage);
Mat ker = notch_kernel( complexIm, notch_rect );
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;
notch_rect.push_back(g_rectangle);
}
break;
}
}
如果觉得《OpenCV 陷波滤波器消除周期性噪音 C++》对你有帮助，请点赞、收藏，并留下你的观点哦！
