快一个月时间没有写博客了,为啥嘞?一方面专业课的考试可不能挂,另一方面的话最近泡了个女朋友,而且是十分粘人那种 = =!。废话不多说,进入正题。
数字图像本质就是一种信号,那信号自然就有频率。在数字图像中,频率就是指灰度变化的速度。并且数字图像信号是离散的,那么要分析频率域时,就要用到离散傅里叶变换及其逆变换之类的。本文公式主要来自《冈萨雷斯》。
某知乎大神做的图
前言
如果读者跟我一样,是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
低通滤波
高通滤波
作者:Jacob杨帮帮
链接:https://www.jianshu.com/p/286050594d9e