OpenCV使用陷波滤波器减少摩尔波纹 C++

2023-11-09 05:50

本文主要是介绍OpenCV使用陷波滤波器减少摩尔波纹 C++,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

陷波滤波器是更有用的选择性滤波器。陷波滤波器拒绝事先定义的关于矩形中心的一个邻域的频率。

零相移滤波器必须是关于原点对称的,因此一个中心位于(u0,v0)的陷波在位置(-u0,-v0)必须有一个对应的陷波。

陷波带阻滤波器可以用中心已被平移到陷波滤波器中心的高通滤波器的乘积来构造。一般形式为:

                                                      H_{NR}(u,v) = \prod_{k=1}^{Q}H_{k}(u,v)H_{-k}(u,v)

其中,H_{k}(u,v)H_{-k}(u,v)是高通滤波器,他们的中心分别位于(u_k, v_k)和(-u_k, -v_k);这些中心是根据频率矩形的中心(M/2,N/2)确定的。对于每个滤波器,距离计算由下式执行:

                                                   D_{k}(u,v) = \left [ (u-M/2-u_{k})^{2} + (v-N/2-v_{k})^{2} \right ]^{1/2}

                                                   D_{-k}(u,v) = \left [ (u-M/2+u_{k})^{2} + (v-N/2+v_{k})^{2} \right ]^{1/2}

例如下面是一个n阶布特沃斯陷波带阻滤波器,它包含三个陷波对;

                                   H_{NR}(u,v) = \prod_{k=1}^{3}\left [ \frac{1}{1+\left [ D_{0k}/D_{k}(u,v) \right ]^{2n}} \right ]\left [ \frac{1}{1+\left [ D_{0k}/D_{-k}(u,v) \right ]^{2n}} \right ]

D_{k} 和 D_{-k} 由 上式计算得出; 常数D_{0k} 对每一个陷波对都是相同的。对于不同陷波对可以不同。其他陷波带阻滤波器可用相同方法构建。

前面说过,1-带阻滤波就是带通,因此陷波带通滤波器:

                              H_{NP}(u,v)=1-H_{NR}(u,v)

 

陷波滤波是选择性的修改DFT的局部区域。典型处理是交互式完成,它直接对DFT处理,

下面使用陷波滤波减少图像的摩尔波纹。

使用布特沃斯 n=4,D0=30处理

原图像             ----------------------    处理后图像

原图像频谱图    ---------------------------------------         处理后图像频谱图 

陷波滤波器

通过鼠标选择矩形区域,找出最大值点,当做(u,v);

代码实现:

#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::Point> &notch_pot, float D0 )
{
    cv::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<float>(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;
                int v_k = notch_pot.at(k).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;
//                std::cout <<  "notch_pot: " << notch_pot.at(k) << std::endl;
            }
            p[2*j]   = h_nr;
            p[2*j+1] = h_nr;

        }
    }

    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滤波器d0=" + std::to_string(D0);
    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::Point>  notch_point;

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;
    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;
        }

    }
    if( !notch_point.empty() ){
        for( unsigned int i = 0; i < notch_point.size(); i++ ){
            cv::circle( tempImage, notch_point.at(i), 3, cv::Scalar(1), 2);
            std::cout << notch_point.at(i) << std::endl;
        }
    }
    cv::imshow(name_win, tempImage);

    Mat ker = notch_kernel( complexIm,notch_point, 30 );
    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;

        if( g_rectangle.height > 0 && g_rectangle.width > 0 ){
            Mat imageROI = image(g_rectangle).clone();
            double  min, max;
            cv::Point   min_loc, 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( "ROI", imageROI );
            cv::imshow( "src", image );
        }

    }
        break;
    }
}

cv::Mat notchFilter_BTW(int rows,int cols,std::vector<cv::Point> np,
                float* D,int n=1,int cvtype=CV_32FC2)
{
    cv::Mat filt(rows,cols,cvtype,cv::Scalar::all(0));
    int cx=cols/2,cy=rows/2;
    int numNotch=np.size();
    float* D2=D;
    for(int i=0;i<numNotch;i++)
    {
        D2[i]=D[i]*D[i];
    }
    int uk[numNotch],vk[numNotch];//在画面上的实际陷波坐标点
    int u_k[numNotch],v_k[numNotch];//陷波共轭点
    float Dk[numNotch],D_k[numNotch];//陷波半径r
    float Hk[numNotch],H_k[numNotch];

    for(int i=0;i<rows;i++)
    {
        for(int j=0;j<cols;j++)
        {
            int u=cx-j,v=cy-i;//中心坐标
            for(int s=0;s<numNotch;s++)
            {
                uk[s]=u+np[s].x,vk[s]=v+np[s].y;
                Dk[s]=uk[s]*uk[s]+vk[s]*vk[s];//距离中心半径的平方
                Hk[s]=1-1/(1+pow(Dk[s]/D2[s],n));

                u_k[s]=u-np[s].x,v_k[s]=v-np[s].y;
                D_k[s]=u_k[s]*u_k[s]+v_k[s]*v_k[s];
                H_k[s]=1-1/(1+pow(D_k[s]/D2[s],n));
            }
            //.data返回的是uchar*型指针,所以要强制转换成浮点数型
            float* p=(float*)(filt.data+i*filt.step[0]+j*filt.step[1]);

            for(int c=0;c<filt.channels();c++)
            {
                p[c]=Hk[0]*H_k[0];
                for(int k=1;k<numNotch;k++)
                {
                    p[c]*=Hk[k]*H_k[k];
                }
            }

        }
    }
    return filt;
}


 

其中,M/2-u_k,N/2- v_k就是计算出来的坐标值,不需要用到M/2, N/2。  

 

这篇关于OpenCV使用陷波滤波器减少摩尔波纹 C++的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/374300

相关文章

中文分词jieba库的使用与实景应用(一)

知识星球:https://articles.zsxq.com/id_fxvgc803qmr2.html 目录 一.定义: 精确模式(默认模式): 全模式: 搜索引擎模式: paddle 模式(基于深度学习的分词模式): 二 自定义词典 三.文本解析   调整词出现的频率 四. 关键词提取 A. 基于TF-IDF算法的关键词提取 B. 基于TextRank算法的关键词提取

使用SecondaryNameNode恢复NameNode的数据

1)需求: NameNode进程挂了并且存储的数据也丢失了,如何恢复NameNode 此种方式恢复的数据可能存在小部分数据的丢失。 2)故障模拟 (1)kill -9 NameNode进程 [lytfly@hadoop102 current]$ kill -9 19886 (2)删除NameNode存储的数据(/opt/module/hadoop-3.1.4/data/tmp/dfs/na

Hadoop数据压缩使用介绍

一、压缩原则 (1)运算密集型的Job,少用压缩 (2)IO密集型的Job,多用压缩 二、压缩算法比较 三、压缩位置选择 四、压缩参数配置 1)为了支持多种压缩/解压缩算法,Hadoop引入了编码/解码器 2)要在Hadoop中启用压缩,可以配置如下参数

Makefile简明使用教程

文章目录 规则makefile文件的基本语法:加在命令前的特殊符号:.PHONY伪目标: Makefilev1 直观写法v2 加上中间过程v3 伪目标v4 变量 make 选项-f-n-C Make 是一种流行的构建工具,常用于将源代码转换成可执行文件或者其他形式的输出文件(如库文件、文档等)。Make 可以自动化地执行编译、链接等一系列操作。 规则 makefile文件

【C++ Primer Plus习题】13.4

大家好,这里是国中之林! ❥前些天发现了一个巨牛的人工智能学习网站,通俗易懂,风趣幽默,忍不住分享一下给大家。点击跳转到网站。有兴趣的可以点点进去看看← 问题: 解答: main.cpp #include <iostream>#include "port.h"int main() {Port p1;Port p2("Abc", "Bcc", 30);std::cout <<

使用opencv优化图片(画面变清晰)

文章目录 需求影响照片清晰度的因素 实现降噪测试代码 锐化空间锐化Unsharp Masking频率域锐化对比测试 对比度增强常用算法对比测试 需求 对图像进行优化,使其看起来更清晰,同时保持尺寸不变,通常涉及到图像处理技术如锐化、降噪、对比度增强等 影响照片清晰度的因素 影响照片清晰度的因素有很多,主要可以从以下几个方面来分析 1. 拍摄设备 相机传感器:相机传

C++包装器

包装器 在 C++ 中,“包装器”通常指的是一种设计模式或编程技巧,用于封装其他代码或对象,使其更易于使用、管理或扩展。包装器的概念在编程中非常普遍,可以用于函数、类、库等多个方面。下面是几个常见的 “包装器” 类型: 1. 函数包装器 函数包装器用于封装一个或多个函数,使其接口更统一或更便于调用。例如,std::function 是一个通用的函数包装器,它可以存储任意可调用对象(函数、函数

C++11第三弹:lambda表达式 | 新的类功能 | 模板的可变参数

🌈个人主页: 南桥几晴秋 🌈C++专栏: 南桥谈C++ 🌈C语言专栏: C语言学习系列 🌈Linux学习专栏: 南桥谈Linux 🌈数据结构学习专栏: 数据结构杂谈 🌈数据库学习专栏: 南桥谈MySQL 🌈Qt学习专栏: 南桥谈Qt 🌈菜鸡代码练习: 练习随想记录 🌈git学习: 南桥谈Git 🌈🌈🌈🌈🌈🌈🌈🌈🌈🌈🌈🌈🌈�

【C++】_list常用方法解析及模拟实现

相信自己的力量,只要对自己始终保持信心,尽自己最大努力去完成任何事,就算事情最终结果是失败了,努力了也不留遗憾。💓💓💓 目录   ✨说在前面 🍋知识点一:什么是list? •🌰1.list的定义 •🌰2.list的基本特性 •🌰3.常用接口介绍 🍋知识点二:list常用接口 •🌰1.默认成员函数 🔥构造函数(⭐) 🔥析构函数 •🌰2.list对象

06 C++Lambda表达式

lambda表达式的定义 没有显式模版形参的lambda表达式 [捕获] 前属性 (形参列表) 说明符 异常 后属性 尾随类型 约束 {函数体} 有显式模版形参的lambda表达式 [捕获] <模版形参> 模版约束 前属性 (形参列表) 说明符 异常 后属性 尾随类型 约束 {函数体} 含义 捕获:包含零个或者多个捕获符的逗号分隔列表 模板形参:用于泛型lambda提供个模板形参的名