zoukankan      html  css  js  c++  java
  • opencv的频域滤波

    下面是频域滤波示例程序:

    在本程序中,共有五个自定义函数,分别是:
    1. myMagnitude(),在该函数中封装了Opencv中的magnitude函数,实现对于复数图像的幅值计算。
    2. dftshift(),该函数实现对图像四个象限的对角互换,相当于MatLab中 fftshift(),将频谱的原点(0,0)移到图像中心。示例1中采用了该函数实现了频谱图中心化。
    3. srcCentralized()用于傅里叶变换前的预处理,以便得到傅里叶频谱的原点(0,0)位于图像的中心。
    该函数与dftshift()目的一致,实现方法不同,一个是变换前预处理,一个是变换后处理。
    示例2中采用了该函数实现了频谱图中心化。
    4. displayDftSpectrum(),该函数用于显示复数图像。
    5. makeFilter(),用于制作频域滤波器,该函数利用了ptr()
       指针遍历图像的方法,实现了圆等滤波函数。
    下面是两个示例,分别采用了dftshift()、srcCentralized()实现频谱图的中心化。示例2与冈萨雷斯的《数字图像处理(第3版)》的4.6.7小结中的流程一致。
    在主程序内容,参见注释。
    示例1:在该程序中采用了第2个函数dftshift(),实现了频谱图中心化。
    
    
    #include <iostream>
    #include<opencv2/opencv.hpp>
    using namespace cv;
    using namespace std;
    
    void myMagnitude(Mat & complexImg,Mat & mI)
    {
        Mat planes[2];
        split(complexImg,planes);
        magnitude(planes[0],planes[1],mI);
    }
    void dftshift(Mat& ds)
    {
        int cx=ds.cols/2;//图像的中心点x 坐标
        int cy=ds.rows/2;//图像的中心点y 坐标
        Mat q0=ds(Rect(0,0,cx,cy));//左上
        Mat q1=ds(Rect(cx,0,cx,cy));//右上
        Mat q2=ds(Rect(0,cy,cx,cy));//左下
        Mat q3=ds(Rect(cx,cy,cx,cy));//右下
        Mat tmp;
        q0.copyTo(tmp);
        q3.copyTo(q0);
        tmp.copyTo(q3);
        q1.copyTo(tmp);
        q2.copyTo(q1);
        tmp.copyTo(q2);
    }
    void displayDftSpectrum(Mat&dftDst,String winName,bool inverseSpectrum)
    {
        Mat magI;
        myMagnitude(dftDst,magI);
        if(!inverseSpectrum)//如果是正向傅里叶变换谱
        {
            magI+=Scalar::all(1);
            log(magI,magI);
        }
        normalize(magI,magI,0,1,NORM_MINMAX);
        imshow(winName,magI);
    }
    void makeFilter(Mat&filter,int R,bool isLowPassFilter)
    {
        int cx=filter.cols/2,cy=filter.rows/2;
        int R2=R*R;
    
        for(int i=0;i<filter.rows;i++)
        {
            Vec2d* pf=filter.ptr<Vec2d>(i);
            for(int j=0;j<filter.cols;j++)
            {
                int r2=(j-cx)*(j-cx)+(i-cy)*(i-cy);
                if(r2<R2)
                {
                    pf[j][0]=1;
                    pf[j][1]=pf[j][0];
                }
            }
        }
        if(!isLowPassFilter)
        {
         filter =Scalar::all(1)-filter;
        }
        Mat displayFilter;
        extractChannel(filter,displayFilter,0);
        imshow("filter image",displayFilter);
    }
    
    int main()
    {
        ///读入灰度图像
        Mat src=imread("D:\Qt\MyImage\baboon.jpg",0);
        int ny=src.rows,nx=src.cols;
        cv::copyMakeBorder(src,src,0,ny,0,nx,BORDER_CONSTANT);
        imshow("original image",src);
        /***************傅里叶变换*************/
        src.convertTo(src,CV_64FC1);//转成浮点数据类型
        Mat dftDst(src.size(),CV_64FC2);//预设dft的输出结果矩阵
        dft(src,dftDst,DFT_COMPLEX_OUTPUT,0);//离散傅立叶变换
        dftshift(dftDst);//将傅里叶变换的结果,四象限对角互换
        displayDftSpectrum(dftDst,"dftSpectrum display",false);//显示傅里叶频谱图
    
        /***************频域滤波*************/
        Mat filter(src.size(),CV_64FC2,Scalar::all(0));
        makeFilter(filter,80,false);
        Mat fftTemp=dftDst.mul(filter);//复数频谱图与滤波器相乘,实现频域滤波
        Mat srcFiltered;
        dft(fftTemp,srcFiltered,DFT_INVERSE);//逆傅里叶变换
        //显示经过滤波后的图像
        displayDftSpectrum(srcFiltered,"filtered image",true);//显示逆傅里叶频谱图
        waitKey();
        return 0;
    }

     示例2:在该示例中采用了第3个函数srcCentralized(),实现频谱图中心化。程序流程即

    在本教材的6.2.3频域滤波步骤小结中,或者在冈萨雷斯的《数字图像处理》4.6.7节中,频域滤波流程总结如下:

      1. 给定一幅大小为M×N的输入图像f(x,y),从式(6.1-25)和式(6.1-26)得到填充参数P和Q。典型地,我们选择P=2M和Q=2N;
      2. 对f(x,y)添加必要数量的0,形成大小为P×Q填充后的图像
      3. 用(-1)(x+y)乘以fp(x,y),进行频谱中心化的预处理;
      4. 计算中心化预处理过的fp(x,y)的傅里叶变换,得到Fp(u,v);
      5. 生成一个实的、对称的滤波函数H(u,v),其大小为P×Q,频谱零点位于(P/2,Q/2)处。用阵列相乘形成乘积G(u,v)=F(u,v)H(u,v);
      6. 经过频域滤波后的图像:
              gp(x,y)={real[fft-1[G(u,v)]]}(-1)(x+y)
        其中,为忽略由于计算不准确导致的寄生复分量,选择了实部,下标P指出我们处理的是填充后的阵列;
      7. 通过从 gp(x,y)的左上象限提取M×N区域,得到最终处理结果g(x,y)。
    /*求取复数矩阵的幅值*/
    void myMagnitude(Mat & complexImg,Mat & mI)
    {
        Mat planes[2];
        split(complexImg,planes);
        magnitude(planes[0],planes[1],mI);
    }
    /*傅里叶变换后的频谱图后处理,将傅里叶普的原点(0,0)平移到图像的中心*/
    void dftshift(Mat& ds)
    {
        int cx=ds.cols/2;//图像的中心点x 坐标
        int cy=ds.rows/2;//图像的中心点y 坐标
        Mat q0=ds(Rect(0,0,cx,cy));//左上
        Mat q1=ds(Rect(cx,0,cx,cy));//右上
        Mat q2=ds(Rect(0,cy,cx,cy));//左下
        Mat q3=ds(Rect(cx,cy,cx,cy));//右下
        Mat tmp;
        q0.copyTo(tmp);
        q3.copyTo(q0);
        tmp.copyTo(q3);
        q1.copyTo(tmp);
        q2.copyTo(q1);
        tmp.copyTo(q2);
    }
    /*傅里叶变换前的预处理,以便频谱图的原点(0,0)移动到图像的中心*/
    void srcCentralized(Mat& src)
    {
        for(int i=0;i<src.rows;i++)
        {
            for(int j=0;j<src.cols;j++)
            {
                float* mv=CV_MAT_ELEM2(src,float,i,j);
                if((i+j)%2!=0)
                    mv[0]=-mv[0];//如果i+j为奇数,该像素值取负值
            }
        }
    }
    /*在窗口中显示复数图像,如果是正向傅里叶矩阵,需要取log才能显示更多频谱信息
     *如果是逆傅里叶变换,通过normalize归一化后,显示频谱图*/
    void displayDftSpectrum(Mat&dftDst,String winName,bool inverseSpectrum)
    {
        Mat magI;
        myMagnitude(dftDst,magI);
        if(!inverseSpectrum)//如果是正向傅里叶变换谱
        {
            magI+=Scalar::all(1);
            log(magI,magI);
        }
        normalize(magI,magI,0,1,NORM_MINMAX);
        imshow(winName,magI);
    }
    void makeFilter(Mat&filter,int R,bool isLowPassFilter)
    {
        int cx=filter.cols/2,cy=filter.rows/2;
        int R2=R*R;
    
        for(int i=0;i<filter.rows;i++)
        {
            Vec2d* pf=filter.ptr<Vec2d>(i);
            for(int j=0;j<filter.cols;j++)
            {
                int r2=(j-cx)*(j-cx)+(i-cy)*(i-cy);
                if(r2<R2)
                {
                    pf[j][0]=1;
                    pf[j][1]=pf[j][0];
                }
            }
        }
        if(!isLowPassFilter)
        {
            filter =Scalar::all(1)-filter;
        }
        Mat displayFilter;
        extractChannel(filter,displayFilter,0);
        imshow("filter image",displayFilter);
    }
    
    int main()
    {
        ///读入灰度图像
        Mat src=imread("D:\Qt\MyImage\interferometer2.jpg",0);
        int ny=src.rows,nx=src.cols;
        imshow("original image",src);
        src.convertTo(src,CV_32FC1);
        srcCentralized(src);
        cv::copyMakeBorder(src,src,0,ny,0,nx,BORDER_CONSTANT);
    
        /***************傅里叶变换*************/
        src.convertTo(src,CV_64FC1);//转成浮点数据类型
        Mat dftDst(src.size(),CV_64FC2);//预设dft的输出结果矩阵
        dft(src,dftDst,DFT_COMPLEX_OUTPUT,0);//离散傅立叶变换
        //    dftshift(dftDst);//将傅里叶变换的结果,四象限对角互换
        displayDftSpectrum(dftDst,"dftSpectrum display",false);//显示傅里叶频谱图
    
        /***************频域滤波*************/
        Mat filter(src.size(),CV_64FC2,Scalar::all(0));
        makeFilter(filter,30,true);
        Mat fftTemp=dftDst.mul(filter);//复数频谱图与滤波器相乘,实现频域滤波
        srcCentralized(fftTemp);
    
        Mat srcFiltered;
        dft(fftTemp,srcFiltered,DFT_INVERSE);//逆傅里叶变换
        //  显示经过滤波后的图像
        displayDftSpectrum(srcFiltered,"filtered image",true);//显示逆傅里叶频谱图
        waitKey();
        return 0;
    }


  • 相关阅读:
    redis状态与性能监控
    redis-stat 安装
    Redis-stat is not found
    查看Redis信息和状态
    查看、分析memcached使用状态
    Memcache内存分配策略
    memcached server LRU 深入分析
    Memcached常用命令及使用说明
    Web-超大文件上传-如何上传文件-大文件上传
    PHP-超大文件上传-如何上传文件-大文件上传
  • 原文地址:https://www.cnblogs.com/phoenixdsg/p/6949575.html
Copyright © 2011-2022 走看看