使用陷波滤波减少莫尔 (波纹) 模式（C++）

        处理图片中的莫尔波纹，需要用到陷波滤波器，本文使用巴特沃斯滤波器设计滤波器，并且可以通过opencv库调用鼠标callback操作，捕捉频域中的噪音点。在频域对图片进行处理，达到很好的滤波效果。 

#include 
#include 
#include 
//#include 
//#include 
using namespace cv;
using namespace std;
void DFT(Mat& input_image,Mat& output_image,Mat& transform_image)   //傅立叶变化函数
{
    for (int i = 0; i < input_image.rows; i++)
    {
        float* p = input_image.ptr(i);
        for (int j = 0; j < input_image.cols; j++)
        {
            p[j]=p[j] * pow(-1 , i+j);
        }
    }
    //创建一个双通道矩阵planes，用来储存复数的实部与虚部
    Mat planes[] = { Mat_(input_image), Mat::zeros(input_image.size(), CV_32FC1)  }; //创建一个双通道矩阵planes，用来储存复数的实部与虚部， CV_32F表示OPENCV里面的32位浮点数
    merge(planes, 2, transform_image);//从多个单通道数组中创建一个多通道数组:transform_image。函数Merge将几个数组合并为一个多通道阵列，即输出数组的每个元素将是输入数组元素的级联
    dft(transform_image, transform_image);//傅立叶变换
    split(transform_image, planes); // 将双通道分为两个单通道，一个表示实部，一个表示虚部
    magnitude(planes[0], planes[1], output_image); //计算复数的幅值，保存在output_image（频谱图）
    //前面得到的频谱图数级过大，不好显示，因此转换
    output_image += Scalar(1);   // 取对数前将所有的像素都加1，防止log0
    log(output_image, output_image);   // 取对数
    normalize(output_image, output_image, 0, 1, NORM_MINMAX); //归一化
    output_image = output_image(Rect(0, 0, output_image.cols & -2, output_image.rows & -2));      //如果有奇数行或列，则对频谱进行裁剪
}
Mat butterworth_filter_kernel(Mat& scr, std::vector& notch_pot, float D0)
{
    //vector表示C++里的一种容器，Point是opencv里一种数据结构表示一个二维点的坐标
    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(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;//陷点的y坐标
                int v_k = notch_pot.at(k).x;//陷点的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;
            }
            p[2 * j] = h_nr;
            p[2 * j + 1] = h_nr;
            //双通道数组，图像频率矩阵的实部和虚部都要乘
        }
    }
    Mat temp[] = { cv::Mat::zeros(scr.size(), CV_32FC1),cv::Mat::zeros(scr.size(), CV_32FC1) };
    split(notch_pass, temp);//双通道分开
    Mat show;   //用来显示幅值
    normalize(temp[0], show, 1, 0, NORM_MINMAX);
    imshow("butterworth_filter_kernel",show);
    return notch_pass;
}
Rect g_rectangle(60, 60, 100, 100);
bool g_bDrawingBox = false; //是否进行绘制
std::vector  notch_point;
void on_MouseHandle(int event, int x, int y, int flags, void* param)
{
    Mat& image = *(cv::Mat*)param;
    switch (event)
    {
        //鼠标右键按下事件
        case cv::EVENT_MOUSEMOVE:
        {   //鼠标键按住时根据起始点画生成的矩形并绘制图片
            if (g_bDrawingBox)
            {
                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);//定义矩形的起点
        }
        break;
        case cv::EVENT_LBUTTONUP:
        {    //鼠标左键松开事件
            g_bDrawingBox = false;
        //    bool w_less_0 = false, h_less_0 = false;
            if (g_rectangle.width<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;
            }
            if (g_rectangle.height > 0 && g_rectangle.width > 0)
            {
                //图像的大小大于一
                Mat imageROI = image(g_rectangle).clone();
                double min;//矩阵最小值
                double max;//矩阵最大值
                Point   min_loc;//矩阵最小值坐标
                Point   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("滤波器陷点所在位置", image);
            }
        }
        break;
    }
}
void esc()// 按下esc键退出
{
    int key_value = -1;
    while (1)
    {
        key_value = cv::waitKey(10);
        if (key_value == 27)
        {
            break;
        }
    }
}
void IDFT(Mat& output_image,Mat& transform_image)  //逆傅立叶
{
    //最终傅立叶逆变换的双通道矩阵
    cv::idft(transform_image, transform_image, DFT_INVERSE);
    Mat dst_plane[2];
    cv::split(transform_image, dst_plane);
    output_image=dst_plane[0];//第一通道的实部
    for (int i = 0; i < output_image.rows; i++)
    {
        float* p = output_image.ptr(i);
        for (int j = 0; j < output_image.cols; j++)
        {
            p[j]=p[j] * pow(-1 , i+j);
        }
    }
    cv::normalize(output_image, output_image, 1, 0, NORM_MINMAX);
}
void amplitude_display(Mat& input_image,Mat& transform_image)  //显示双通道频谱的幅度图像
{
    Mat amplitude;
    Mat planes[] = { Mat_(input_image), Mat::zeros(input_image.size(), CV_32FC1)  };
    split(transform_image, planes); // 将双通道分为两个单通道，一个表示实部，一个表示虚部
    magnitude(planes[0], planes[1], amplitude);
    amplitude += Scalar(1);   // 取对数前将所有的像素都加1，防止log0
    log(amplitude, amplitude);   // 取对数
    normalize(amplitude, amplitude, 0, 1, NORM_MINMAX); //归一化
    imshow("amplitude",amplitude);
    waitKey(0);
}
int main(int argc, char** argv) {
    Mat image;
    image = imread("/Users/huangziyu/Desktop/zhn.jpg", 0);
    Mat transform_image;     //傅立叶变换结果
    Mat image_output;     //傅立叶变换幅度结果
    imshow("image",image);
    waitKey(0);
    int m = getOptimalDFTSize(2*image.rows);
    int n = getOptimalDFTSize(2*image.cols);
    copyMakeBorder(image, image, 0, m - image.rows, 0, n - image.cols, BORDER_CONSTANT, Scalar::all(0));        //扩展图像矩阵，为2，3，5的倍数时运算速度快
    image.convertTo(image, CV_32FC1);
    Mat A;
    normalize(image,A, 1 , 0, NORM_MINMAX);
    imshow("image",A);
    waitKey(0);
    DFT(image, image_output,transform_image);    //傅立叶变换
    imshow("image_output", image_output);
    cv::setMouseCallback("image_output", on_MouseHandle, (void*)& image_output);  //鼠标返回函数
    esc();
    Mat filter = butterworth_filter_kernel(transform_image, notch_point, 30);
    cv::multiply(transform_image, filter, transform_image);  //图像频谱和滤波器相乘
    amplitude_display(image,transform_image);
    Mat final_image;  //最终图片
    IDFT(final_image,transform_image);
    imshow("final_image",final_image(Rect(0,0,235,340)));
    waitKey(0);
    return 0;
}

        处理结果如下：

频率图：

频率中噪音点：

 巴特沃斯滤波器过滤后的频谱：

最终结果： 

zhn.jpg图片处理： 完美！