OpenCV 陷波滤波器消除周期性噪音 C++

观察下面图像：

有明显的近似水平线的噪音模式，我们希望他在频率域内的成分沿垂直轴集中；

下图是该图像的谱；

将垂直轴的分量去掉，构建的陷波器如下：

原图与结果图对比：

滤掉的空间噪音模式：

 

代码实现：

#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> &notch_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;
    }
}