图像处理之离散傅里叶变换2

接着上一节继续讲述傅里叶变换

4 返回DFT最优尺寸大小：getOptimalDFTSize函数

getOptimalDFTSize函数返回给定向量尺寸的傅里叶最优尺寸大小。为了提高离散傅里叶变换的运行速度，需要扩充图像，而具体扩充多少，就由这个函数来计算得到。

C++

int getOptimalDFTSize(int vecsize)

此函数的唯一一个参数为int类型的vecsize，向量尺寸，即图片的rows、cols。

5 扩充图像边界：copyMakeBorder函数

copyMakeBorder函数的作用是扩充图像边界。

C++:

void copyMakeBorder (InputArray src, OutputArray dst,int top, int bottom, int left, 
int right, int borderType, const Scalar& value=Scalar())

• 第一个参数，InputArray类型的src，输入图像，即源图像，填Mat类的对象即可。
• 第二个参数，OutputArray 类型的dst，函数调用后的运算结果存在这里，即这个参数用于存放函数调用后的输出结果，需和源图片有一样的尺寸和类型，且size应该为Size（src.cols＋left＋right，src.rows＋top＋bottom）。
• 接下来的4个参数分别为int类型的top、bottom、left、right，分别表示在 源图像的四个方向上扩充多少像素，例如top＝2，bottom＝2，left＝2，right＝2 就意味着在源图像的上下左右各扩充两个像素宽度的边界。
• 第七个参数，borderType 类型的，边界类型，常见取值为BORDER＿CONSTANT，可参考 borderInterpolate（）得到更多的细节。
• 第八个参数，const Scalar＆类型的value，有默认值Scalar（），可以理解为默 认值为0。当 borderType取值为BORDER＿CONSTANT时，这个参数表示 边界值。

6 计算二维矢量的幅值：magnitude函数

magnitude()函数用于计算二维矢量的幅值。

C++:

void magnitude (InputArray x, InputArray y, OutputArray magnitude)

• 第一个参数，InputArray类型的x，表示矢量的浮点型X坐标值，也就是实部。
• 第二个参数，InputArray类型的y，表示矢量的浮点型Y坐标值，也就是虚部。
• 第三次参数，OutputArray 类型的 magnitude，输出的幅值，它和第一个参 数x有着同样的尺寸和类型。

下式可以表示 magnitude（）函数的原理：

7 计算自然对数：log函数

函数log()函数的功能是计算每个数组元素绝对值的自然对数。

C++:

void log(InputArray src, OutArray dst)

第一个参数为输入图像，第二个参数为得到的对数值。其原理如下。

8 矩阵归一化：normalize函数

normalize（）的作用是进行矩阵归一化。

C++:

void normalize(InputArray src, OutputArray dst, double alpha=1, double beta=0, 
    int norm_type=NORM_L2, int dtype = -1,InputArray mask=noArray())

• 第一个参数，InputArray类型的src。输入图像，即源图像，填Mat类的对象即可。
• 第二个参数，OutputArray 类型的dst。函数调用后的运算结果存在这里，和源图片有一样的尺寸和类型。
• 第三个参数，double类型的alpha。归一化后的最大值，有默认值1。
• 第四个参数，double类型的 beta。归一化后的最大值，有默认值0。
• 第五个参数，int类型的norm＿type。归一化类型，有NORM＿INF、NORM＿L1、 NORM＿L2和NORM＿MINMAX等参数可选，有默认值NORM＿L2。
• 第六个参数，int类型的dtype，有默认值—1。当此参数取负值时，输出矩阵和src有同样的类型，否则，它和src有同样的通道数，且此时图像深度为CV_MAT_DEPTH(dtype)。
• 第七个参数，InputArray类型的mask，可选的操作掩膜，有默认值noArray（）。

9 示例程序：离散傅里叶变换

这节中我们将学习一个以dft（函数为核心，对图像求傅里叶变换的，有详细注释的示例程序。

在此示例中，将展示如何计算以及显示傅里叶变换后的幅度图像。由于数字图像的离散性，像素值的取值范围也是有限的。比如在一张灰度图像中，像素灰度值一般在0到255之间。因此，我们这里讨论的也仅仅是离散傅里叶变换（DFT）。如果需要得到图像中的几何结构信息，那么就要用到离散傅里叶变换了。下面的步骤将以输入图像为单通道的灰度图像I为例，进行分步说明。

1．【第一步】载入原始图像

我们在这一步以灰度模式读取原始图像，进行是否读取成功的检测，并显示出读取到的图像。代码如下。

2．【第二步】将图像扩大到合适的尺寸

离散傅里叶变换的运行速度与图片的尺寸有很大关系。当图像的尺寸是2、3、5 的整数倍时，计算速度最快。因此，为了达到快速计算的目的，经常通过添凑新的边缘像素的方法获取最佳图像尺寸。函数 getOptimalDFTSize（）用于返回最佳尺寸，而函数copyMakeBorder（）用于填充边缘像素，这一步代码如下。

3．【第三步】为傅里叶变换的结果（实部和虚部）分配存储空间

傅里叶变换的结果是复数，这就是说对于每个原图像值，结果会有两个图像值。此外，频域值范围远远超过空间值范围，因此至少要将频域储存在 float 格式中。所以我们将输入图像转换成浮点类型，并多加一个额外通道来储存复数部分。

4．【第四步】进行离散傅里叶变换

这里的离散傅里叶变换为图像就地计算模式（in—place，输入输出为同一图像）。

5．【第五步】将复数转换为幅值

复数包含实数部分（Re）和虚数部分（imaginary—Im）。离散傅里叶变换的结果是复数，对应的幅度可以表示为：

那么，转化为OpenCV代码，就是下面这样的。

6．【第六步】进行对数尺度（logarithmic scale）缩放．

傅里叶变换的幅度值范围大到不适合在屏幕上显示。高值在屏幕上显示为白点，而低值为黑点，高低值的变化无法有效分辨。为了在屏幕上凸显出高低变化的连续性，我们可以用对数尺度来替换线性尺度，公式如下。

而写成OpenCV代码，就是下面这样的。

7．【第七步】剪切和重分布幅度图象限

因为在第二步中延扩了图像，那现在是时候将新添加的像素剔除了。为了方便显示，也可以重新分布幅度图象限位置（注：将第五步得到的幅度图从中间划开，得到4张1／4子图像，将每张子图像看成幅度图的一个象限，重新分布，即将4个角点重叠到图片中心）。这样的话原点（0，0）就位移到图像中心了。OpenCV代码如下。

8．【第八步】归一化

这一步仍然是为了显示。现在有了重分布后的幅度图，但是幅度值仍然超过可显示范围［0，1］。我们使用 normalize（）函数将幅度归一化到可显示范围。

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

处理完成，最后进行显示即可。

将上述代码组合到一起，变得到了程序的完整的代码 。

//---------------------------------【头文件、命名空间包含部分】-----------------------------
//    描述：包含程序所使用的头文件和命名空间
//-------------------------------------------------------------------------------------------------
#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
using namespace cv;


//-----------------------------------【ShowHelpText( )函数】----------------------------------
//     描述：输出一些帮助信息
//----------------------------------------------------------------------------------------------
void ShowHelpText()
{
  //输出欢迎信息和OpenCV版本
  printf("\n\n\t\t\t非常感谢购买《OpenCV3编程入门》一书！\n");
  printf("\n\n\t\t\t此为本书OpenCV2版的第28个配套示例程序\n");
  printf("\n\n\t\t\t   当前使用的OpenCV版本为：" CV_VERSION );
  printf("\n\n  ----------------------------------------------------------------------------\n");
}



//--------------------------------------【main( )函数】-----------------------------------------
//          描述：控制台应用程序的入口函数，我们的程序从这里开始执行
//-------------------------------------------------------------------------------------------------
int main( )
{

  //【1】以灰度模式读取原始图像并显示
  Mat srcImage = imread("1.jpg", 0);
  if(!srcImage.data ) { printf("读取图片错误，请确定目录下是否有imread函数指定图片存在~！ \n"); return false; } 
  imshow("原始图像" , srcImage);   

  ShowHelpText();

  //【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】将复数转换为幅值，即=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
  split(complexI, planes); // 将多通道数组complexI分离成几个单通道数组，planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
  magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude  
  Mat magnitudeImage = planes[0];

  //【6】进行对数尺度(logarithmic scale)缩放
  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));  // ROI区域的右上
  Mat q2(magnitudeImage, Rect(0, cy, cx, cy));  // ROI区域的左下
  Mat q3(magnitudeImage, Rect(cx, cy, cx, cy)); // ROI区域的右下
  //交换象限（左上与右下进行交换）
  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); 

  //【9】显示效果图
  imshow("频谱幅值", magnitudeImage);    
  waitKey();

  return 0;
}

原图

效果图