OpenCV中的图像处理——离散傅里叶变换实例(11)

我们通过这个示例来学习 一个以dft()为核心，对图像求傅里叶变换的过程。

程序示例

如何计算和显示傅里叶变换后的幅度图像。由于数字图像的离散性，我门也仅仅实现的是离散傅里叶变换，通过傅里叶变换得到图像中的几何结构信息。下面以输入图像的单通道灰度图为例。
新函数部分参考上讲的内容： https://editor.csdn.net/md/?articleId=124288761
1.【第一步】载入图像
以灰度图模式读取原始图像，进行是否读取成功的判断，并显示图像。

Mat srcImage = imread("1.jpg", 0);
if(!srcImage.data) {
	printf("获取图像失败！\n");
	return false;
}
imshow("原始图像", srcImage);

2 .【第二步】将图像扩大到合适的尺寸。
离散傅里叶变换的运行速度与图片的尺寸有很大的关系。当图像尺寸是2，3，5的整数倍时，计算速度最快。因此，为了达到快速计算的目的，通常会添凑新的边缘像素的幅方法来获取最佳的图像尺寸。函数getOptimalDFTSize()用于返回最佳尺寸，而函数copyMakeBorder()用于填充边缘像素。

//[2]将输入图像延扩到最佳的尺寸，边界用0补充
int m = getOptimalDFTSize(I.rows);
int n = getOptimalDFTSize(I.cols);
//将添加的像素初始化为0
Mat padded;
copyMakeBorder(I, padded, 0, m-I.rows, 0, n-I.cols, BORDER_CONSTANT, Scalar::all(0));

3 .【第三步】为傅里叶变换的结果(实部和虚部)分配存储空间。
傅里叶变换的结果是复数，对每个图像值会有实和虚部两个图像值。频域值的范围大与时域值范围，因此要将频域值存储在float格式中，所以我们此时就需要将输入图像转为浮点类型,并加一个额为的通道来存储复数部分,

//为傅里叶变换的结果分配(实部和虚部)存储空间
//将palanes数组组合合并成一个多通道的数组complexI
Mat palnes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexI;
merge(palnes, 2, complexI);

4 .【第四步】进行离散傅里叶变换
图像就地计算模式(in-place),输入输出为同一幅图像。

//进行就地离散傅里叶变换
dft(complexI, complecI);

5 .【第五步】将复数转变为幅值
离散傅里叶变换的结果是复数，对应的幅值：$M=\sqrt{Re(DFT(I){)}^2+Im(DFT(I){)}^2}$

//将复数转变为幅值，即log(1+ M)
split(complexI, palnes);   //将多通道数组complex分离成单通道数组
//palnes[0] = Re(DFT(I)), palnes[1] = Im(DFT(I))
magnitude(palnes[0], palnes[1], palnes[0]); // palnes[0] = magnitude
Mat magnitudeImage = planes[0];

6 .【第六步】进行对数尺度缩放
离散傅里叶的幅度范围不适合在屏幕上显示。高值呈现为白点，低值为黑点，高低值的变化无法有效的分辨，为了在屏幕上显示高低变化的连续性，可用用对数尺度替换线性尺度。
$M_1=log(1+M)$

//进行对数尺度的缩放
magnitudeImage += Scalar::all(1);
log(magnitudeImage, magnitudeImage); //求自然对数

7.【第七步】剪切和重分布幅度图象限
在第二步中扩展了图像，将新添加的像素剔除，为了方便显示，可以重新分布幅度图象限的位置(将第五步中的图从中间划开得到四个子图，每张子图就是幅度图的象限，重新分布，即将四个角点重叠到图像中心)。这样原点就位移到了图像的中心。

//剪切和重新分布幅度图象限
//若有奇数行或奇数列，进行频谱裁减
 magnitudeImage = magnitudeImage(Rect(0,0,magnitudeImage.cols & -2, magnitudeImage.rows & -1));
 //重新排列傅里叶图像中的象限，使得原点位于图像中心
 int cv = magnitudeImage.cols/2;
 int cy = magnitudeImage.rows/2;
 Mat q0(magnitudeImage,Rect(0,0,cx,cy));  //ROI左上角
 Mat q1(magnitudeImage, Rect(cx, 0, cx,cy));  // 右下角
 Mat q2(magnitudeImage, Rect(0,cy, cx, cy));
 Mat q3(magnitudeImage, 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);

8 . 【第八步】归一化
重新分布后的幅度图，其幅度值仍然超过了可示范围[0,1]，我们使用normalize()函数将幅度归一化到可显示的范围

//归一化，用0-1之间的浮点值将矩阵变换到可视的图像格式
//normalize(magnitudeImage, magnitudeImage, 0 , 1, CV_MINMAX);//opencv2
normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX);//opencv3

9.【第九步】显示效果图

//显示频谱图
imshow("频谱图",magnitudeImage);

分步式的综合示例如下：

#include "opencv2/opencv.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include<iostream>
using namespace std;
using namespace cv;

int main()
{
    //[1]以灰度图读取原始图像并显示
    Mat srcImage = imread("../1.jpg", 0);
    if(!srcImage.data) {
        printf("获取图像失败\n");
        return false;
    }
    imshow("原始图", srcImage);

    //[2]将输入图像扩展到最佳尺寸，边界填0
    int m= getOptimalDFTSize(srcImage.rows);
    int n = getOptimalDFTSize(srcImage.cols);
    //将添加的像素初始化为0
    Mat padded;
    copyMakeBorder(srcImage, padded, 0, m - srcImage.rows, 0, n - srcImage.cols, BORDER_CONSTANT, Scalar::all(0));

    //[3]为傅里叶变换的结果分配存储空间(实部和虚部)
    //将planes数组合并成一个多通道的数组complexI
    Mat planes[] = {Mat_<float>(padded), Mat ::zeros(padded.size(), CV_32F)};
    Mat complexI;
    merge(planes, 2, complexI);

    //[4]进行就地离散傅里叶变换
    dft(complexI, complexI);

    //[5]将复数转换为幅度值
    split(complexI, planes);  //将多通道数组分离成几个单通道数组
    // planes[0] = Re(DFT(I));
    // planes[1] = Im(DFT(I));

    magnitude(planes[0], planes[1], planes[0]);

    Mat magnitudeImage = planes[0];

    //[6]进行对数尺度缩放
    magnitudeImage += Scalar::all(1);
    log(magnitudeImage, magnitudeImage);

    //[7]剪切和重分幅度图象限
    //若有奇数行或列，进行频谱裁剪
    magnitudeImage = magnitudeImage(Rect(0, 0, magnitudeImage.cols & -2, magnitudeImage.rows & -2));
    //重新排列傅里叶变换中的象限，使得原点位于图像中心
    int cx = magnitudeImage.cols/2;
    int cy = magnitudeImage.rows/2;
    Mat q0(magnitudeImage,Rect(0,0,cx,cy));  //ROI左上角
    Mat q1(magnitudeImage, Rect(cx, 0, cx,cy));  // 右下角
    Mat q2(magnitudeImage, Rect(0,cy, cx, cy));
    Mat q3(magnitudeImage, 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);

    //[8]归一化，用0-1之间的浮点值将矩阵变换到可视的图像格式
    //normalize(magnitudeImage, magnitudeImage, 0 , 1, CV_MINMAX);//opencv2
    normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX);//opencv3

    imshow("频谱图", magnitudeImage);
    // waitKey(0);
    while(char (waitKey(1)) != 'q'){}
    return 0;

}

运行结果：