OpenCV 频率域处理:离散傅里叶变换、频率域滤波

数字图像本质就是一种信号,那信号自然就有频率。在数字图像中,频率就是指灰度变化的速度。并且数字图像信号是离散的,那么要分析频率域时,就要用到离散傅里叶变换及其逆变换之类的。本文公式主要来自《冈萨雷斯》。


某知乎大神做的图

前言

如果读者跟我一样,是EE专业出身的学生(电子、电气、自动化、通信),学过复变函数、信号与系统之类的课程,那么理论部分可以直接去看《冈萨雷斯》的第四章频率域滤波。因为有相关基础,就算不能完全看懂也能半懂。

如果读者是其它专业出身,比如计算机类和机械类专业,那对这方面的理解可能就会不太容易(信号之类的课程应该是选修?),因为对傅里叶变换的理解可能会有所欠缺。当年要用到傅里叶变换的时候(不是用在图像上),还没有学到相关课程,那个时候可是让我抓耳挠腮。不过很幸运阅读了知乎大神Heinrich的文章,很快就对傅里叶变换有了初步了解,在这里推荐一下。

傅里叶变换

在连续的情况下,傅里叶变换和逆变换很简洁,公式如下


傅里叶变换

反变换

由欧拉公式,傅里叶变换可以写成这样。其中,正弦项的频率由μ决定


欧拉公式 + 傅里叶变换

在离散的情况下,公式如下,也很好理解


正变换

逆变换

上面的都是一维的情况。显然,图像是二维的,于是要推广到二维的傅里叶变换。类似,连续情况下傅里叶变换公式如下


二维傅里叶变换

二维傅里叶逆变换

离散的情况如下。其中图像为M*N的图像。


DFT

IDFT

二维DFT一般是复函数,用极坐标来表示如下


极坐标形式

取幅度就得到了被称为傅里叶谱或者频谱的玩意

幅度

计算一下反正切就得到了相角。


相角

因为是二维信号,所以频谱和相角都是二维的。其中相角包含了频率的位置信息。总的来说,图像经过DFT之后得到了频谱和相角。我们在分析的时候,一般只看频谱。但如果要进行逆变换,就同时需要频谱和相角的信息,才能正确还原图像。

频率域滤波

频率域滤波的基本公式如下。H为频率域上的滤波函数。


频率域滤波公式

常见的频率域滤波函数有理想高低通滤波器、布特沃斯滤波器、高斯滤波器。其中前面两者的滤波效果会发生振铃现象。这里的代码实现只搞下高斯滤波。

高斯低通滤波
中间高四周低,高频被滤除
高斯高通滤波

低频被滤除

编程实现

opencv有用于离散傅里叶变换和逆变换的函数,但还是要稍微处理一些细节,并且频率域的滤波需要自己去编写。

/********************************************************************
 * Created by 杨帮杰 on 12/8/2018
 * Right to use this code in any way you want without
 * warranty, support or any guarantee of it working
 * E-mail: yangbangjie1998@qq.com
 * Association: SCAU 华南农业大学
 ********************************************************************/

#include <iostream>
#include <vector>
#include <opencv2/core.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/features2d.hpp>
#include <opencv2/xfeatures2d.hpp>
#include <opencv2/calib3d.hpp>

#define IMG_PATH "//home//jacob//图片//lenna.jpg"

using namespace std;
using namespace cv;

int main()
{
    Mat inputImg = imread(IMG_PATH, IMREAD_GRAYSCALE);

    if(inputImg.empty())
    {
        cout << "图片没读到,傻逼!" << endl;
        return -1;
    }

    //得到DFT的最佳尺寸(2的指数),以加速计算
    Mat paddedImg;
    int m = getOptimalDFTSize(inputImg.rows);
    int n = getOptimalDFTSize(inputImg.cols);

    cout << "图片原始尺寸为:" << inputImg.cols << "*" << inputImg.rows <<endl;
    cout << "DFT最佳尺寸为:" << n << "*" << m <<endl;

    //填充图像的下端和右端
    copyMakeBorder(inputImg, paddedImg, 0, m - inputImg.rows,
                   0, n - inputImg.cols, BORDER_CONSTANT, Scalar::all(0));

    //将填充的图像组成一个复数的二维数组(两个通道的Mat),用于DFT
    Mat matArray[] = {Mat_<float>(paddedImg), Mat::zeros(paddedImg.size(), CV_32F)};
    Mat complexInput, complexOutput;
    merge(matArray, 2, complexInput);

    dft(complexInput, complexOutput);

    //计算幅度谱(傅里叶谱)
    split(complexOutput, matArray);
    Mat magImg;
    magnitude(matArray[0], matArray[1], magImg);

    //转换到对数坐标
    magImg += Scalar::all(1);
    log(magImg, magImg);

    //将频谱图像裁剪成偶数并将低频部分放到图像中心
    int width = (magImg.cols / 2)*2;
    int height = (magImg.cols / 2)*2;
    magImg = magImg(Rect(0, 0, width, height));

    int colToCut = magImg.cols/2;
    int rowToCut = magImg.rows/2;

    //获取图像,分别为左上右上左下右下
    //注意这种方式得到的是magImg的ROI的引用
    //对下面四幅图像进行修改就是直接对magImg进行了修改
    Mat topLeftImg(magImg, Rect(0, 0, colToCut, rowToCut));
    Mat topRightImg(magImg, Rect(colToCut, 0, colToCut, rowToCut));
    Mat bottomLeftImg(magImg, Rect(0, rowToCut, colToCut, rowToCut));
    Mat bottomRightImg(magImg, Rect(colToCut, rowToCut, colToCut, rowToCut));

    //第二象限和第四象限进行交换
    Mat tmpImg;
    topLeftImg.copyTo(tmpImg);
    bottomRightImg.copyTo(topLeftImg);
    tmpImg.copyTo(bottomRightImg);

    //第一象限和第三象限进行交换
    bottomLeftImg.copyTo(tmpImg);
    topRightImg.copyTo(bottomLeftImg);
    tmpImg.copyTo(topRightImg);

    //归一化图像
    normalize(magImg, magImg, 0, 1, CV_MINMAX);

    //傅里叶反变换
    Mat complexIDFT, IDFTImg;
    idft(complexOutput, complexIDFT);
    split(complexIDFT, matArray);
    magnitude(matArray[0], matArray[1], IDFTImg);
    normalize(IDFTImg, IDFTImg, 0, 1, CV_MINMAX);

    imshow("输入图像", inputImg);
    imshow("频谱图像", magImg);
    imshow("反变换图像", IDFTImg);

    /***********************频率域滤波部分*************************/
    //高斯低通滤波函数(中间高两边低)
    Mat gaussianBlur(paddedImg.size(),CV_32FC2);
    //高斯高通滤波函数(中间低两边高)
    Mat gaussianSharpen(paddedImg.size(),CV_32FC2);
    double D0=2*10*10*10;
    for(int i=0;i<paddedImg.rows;i++)
    {
        float*p=gaussianBlur.ptr<float>(i);
        float*q=gaussianSharpen.ptr<float>(i);
        for(int j=0;j<paddedImg.cols;j++)
        {
            double d=pow(i-paddedImg.rows/2,2)+pow(j-paddedImg.cols/2,2);
            p[2*j]=expf(-d/D0);
            p[2*j+1]=expf(-d/D0);

            q[2*j]=1-expf(-d/D0);
            q[2*j+1]=1-expf(-d/D0);
        }
    }

    //将实部和虚部按照频谱图的方式换位
    //低频在图像中心,用于滤波
    split(complexOutput, matArray);

    //matArray[0]表示DFT的实部
    Mat q1(matArray[0], Rect(0, 0, colToCut, rowToCut));
    Mat q2(matArray[0], Rect(colToCut, 0, colToCut, rowToCut));
    Mat q3(matArray[0], Rect(0, rowToCut, colToCut, rowToCut));
    Mat q4(matArray[0], Rect(colToCut, rowToCut, colToCut, rowToCut));

    //第二象限和第四象限进行交换
    q1.copyTo(tmpImg);
    q4.copyTo(q1);
    tmpImg.copyTo(q4);

    //第一象限和第三象限进行交换
    q2.copyTo(tmpImg);
    q3.copyTo(q2);
    tmpImg.copyTo(q3);

    //matArray[1]表示DFT的虚部
    Mat qq1(matArray[1], Rect(0, 0, colToCut, rowToCut));
    Mat qq2(matArray[1], Rect(colToCut, 0, colToCut, rowToCut));
    Mat qq3(matArray[1], Rect(0, rowToCut, colToCut, rowToCut));
    Mat qq4(matArray[1], Rect(colToCut, rowToCut, colToCut, rowToCut));

    //第二象限和第四象限进行交换
    qq1.copyTo(tmpImg);
    qq4.copyTo(qq1);
    tmpImg.copyTo(qq4);

    //第一象限和第三象限进行交换
    qq2.copyTo(tmpImg);
    qq3.copyTo(qq2);
    tmpImg.copyTo(qq3);

    merge(matArray, 2, complexOutput);

    //滤波器和DFT结果相乘(矩阵内积)
    multiply(complexOutput,gaussianBlur,gaussianBlur);
    multiply(complexOutput,gaussianSharpen,gaussianSharpen);

    //计算频谱
    split(gaussianBlur,matArray);
    magnitude(matArray[0],matArray[1],magImg);
    magImg+=Scalar::all(1);
    log(magImg,magImg);
    normalize(magImg,magImg,1,0,CV_MINMAX);
    imshow("高斯低通滤波频谱",magImg);

    split(gaussianSharpen,matArray);
    magnitude(matArray[0],matArray[1],magImg);
    magImg+=Scalar::all(1);
    log(magImg,magImg);
    normalize(magImg,magImg,1,0,CV_MINMAX);
    imshow("高斯高通滤波频谱",magImg);

    //IDFT得到滤波结果
    Mat gaussianBlurImg;
    idft(gaussianBlur, complexIDFT);
    split(complexIDFT, matArray);
    magnitude(matArray[0], matArray[1], gaussianBlurImg);
    normalize(gaussianBlurImg, gaussianBlurImg, 0, 1, CV_MINMAX);

    Mat gaussianSharpenImg;
    idft(gaussianSharpen, complexIDFT);
    split(complexIDFT, matArray);
    magnitude(matArray[0], matArray[1], gaussianSharpenImg);
    normalize(gaussianSharpenImg, gaussianSharpenImg, 0, 1, CV_MINMAX);

    imshow("高斯低通滤波", gaussianBlurImg);
    imshow("高斯高通滤波", gaussianSharpenImg);


    waitKey(0);

    return 0;
}

输入-DFT-IDFT

低通滤波

高通滤波

Reference:
opencv学习(十五)之图像傅里叶变换dft
opencv 频率域滤波实例
opencv官方例程 :dtf.cpp
《数字图像处理》 ——冈萨雷斯

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 158,736评论 4 362
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 67,167评论 1 291
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 108,442评论 0 243
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 43,902评论 0 204
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 52,302评论 3 287
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,573评论 1 216
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,847评论 2 312
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,562评论 0 197
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,260评论 1 241
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,531评论 2 245
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 32,021评论 1 258
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,367评论 2 253
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 33,016评论 3 235
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,068评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,827评论 0 194
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,610评论 2 274
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,514评论 2 269

推荐阅读更多精彩内容