opencv学 之图像傅里叶变换dft

//*************

在学习信号与系统或通信原理等课程里面可能对傅里叶变换有了一定的了解。我们知道傅里叶变换是把一个信号从时域变换到其对应的频域进行分析。如果有小伙伴还对傅里叶变换处于很迷糊的状态，请戳这里，非常通俗易懂。而在图像处理中也有傅里叶分析的概念，我这里给出在其官方指导文件opencv_tutorials中给出的解释。 
傅里叶变换可以将一幅图片分解为正弦和余弦两个分量，换而言之，他可以将一幅图像从其空间域（spatial domain）转换为频域（frequency domain）。这种变换的思想是任何函数可以很精确的接近无穷个sin()函数和cos()函数的和。傅里叶变换提供了这种方法来达到这种效果。对于二位图像其傅里叶变换公式如下： 
 
式中f(i, j)是图像空间域的值而F是频域的值。傅里叶转换的结果是复数，这也显示出了傅里叶变换是一副实数图像（real image）和虚数图像（complex image）叠加或者是幅度图像（magitude image）和相位图像（phase image）叠加的结果。在实际的图像处理算法中仅有幅度图像（magnitude image）图像能够用到，因为幅度图像包含了我们所需要的所有图像几何结构的信息。但是，如果想通过修改幅度图像或者相位图像来间接修改原空间图像，需要保留幅度图像和相位图像来进行傅里叶逆变换，从而得到修改后图像。

1.dft()

首先看一下opencv提供的傅里叶变换函数dft()，其定义如下：

C++: void dft(InputArray src, OutputArray dst, int flags=0, int nonzeroRows=0);
1
参数解释： 
. InputArray src: 输入图像，可以是实数或虚数 
. OutputArray dst: 输出图像，其大小和类型取决于第三个参数flags 
. int flags = 0: 转换的标识符，有默认值0.其可取的值如下所示： 
。DFT_INVERSE: 用一维或二维逆变换取代默认的正向变换 
。DFT_SCALE: 缩放比例标识符，根据数据元素个数平均求出其缩放结果，如有N个元素，则输出结果以1/N缩放输出，常与DFT_INVERSE搭配使用。 
。DFT_ROWS: 对输入矩阵的每行进行正向或反向的傅里叶变换；此标识符可在处理多种适量的的时候用于减小资源的开销，这些处理常常是三维或高维变换等复杂操作。 
。DFT_COMPLEX_OUTPUT: 对一维或二维的实数数组进行正向变换，这样的结果虽然是复数阵列，但拥有复数的共轭对称性（CCS），可以以一个和原数组尺寸大小相同的实数数组进行填充，这是最快的选择也是函数默认的方法。你可能想要得到一个全尺寸的复数数组（像简单光谱分析等等），通过设置标志位可以使函数生成一个全尺寸的复数输出数组。 
。DFT_REAL_OUTPUT: 对一维二维复数数组进行逆向变换，这样的结果通常是一个尺寸相同的复数矩阵，但是如果输入矩阵有复数的共轭对称性（比如是一个带有DFT_COMPLEX_OUTPUT标识符的正变换结果），便会输出实数矩阵。 
. int nonzeroRows = 0: 当这个参数不为0，函数会假设只有输入数组（没有设置DFT_INVERSE）的第一行或第一个输出数组（设置了DFT_INVERSE）包含非零值。这样的话函数就可以对其他的行进行更高效的处理节省一些时间，这项技术尤其是在采用DFT计算矩阵卷积时非常有效。

2. getOptimalDFTSize()

返回给定向量尺寸经过DFT变换后结果的最优尺寸大小。其函数定义如下：

C++: int getOptimalDFTSize(int vecsize);
1
参数解释： 
int vecsize: 输入向量尺寸大小(vector size) 
DFT变换在一个向量尺寸上不是一个单调函数，当计算两个数组卷积或对一个数组进行光学分析，它常常会用0扩充一些数组来得到稍微大点的数组以达到比原来数组计算更快的目的。一个尺寸是2阶指数（2,4,8,16,32…）的数组计算速度最快，一个数组尺寸是2、3、5的倍数（例如：300 = 5*5*3*2*2）同样有很高的处理效率。 
getOptimalDFTSize()函数返回大于或等于vecsize的最小数值N，这样尺寸为N的向量进行DFT变换能得到更高的处理效率。在当前N通过p,q,r等一些整数得出N = 2^p*3^q*5^r. 
这个函数不能直接用于DCT（离散余弦变换）最优尺寸的估计，可以通过getOptimalDFTSize((vecsize+1)/2)*2得到。

3.magnitude()

计算二维矢量的幅值，其定义如下：

C++: void magnitude(InputArray x, InputArray y, OutputArray magnitude)；
1
参数解释： 
. InputArray x: 浮点型数组的x坐标矢量，也就是实部 
. InputArray y: 浮点型数组的y坐标矢量，必须和x尺寸相同 
. OutputArray magnitude: 与x类型和尺寸相同的输出数组 
其计算公式如下： 

4. copyMakeBorder() 
扩充图像边界，其函数定义如下：

C++: void copyMakeBorder(InputArray src, OutputArray dst, int top, int bottom, int left, int right, int borderType, const Scalar& value=Scalar() )；
1
参数解释： 
. InputArray src: 输入图像 
. OutputArray dst: 输出图像，与src图像有相同的类型，其尺寸应为Size(src.cols+left+right, src.rows+top+bottom) 
. int类型的top、bottom、left、right: 在图像的四个方向上扩充像素的值 
. int borderType: 边界类型，由borderInterpolate()来定义，常见的取值为BORDER_CONSTANT 
. const Scalar& value = Scalar(): 如果边界类型为BORDER_CONSTANT则表示为边界值

5. normalize() 
归一化就是把要处理的数据经过某种算法的处理限制在所需要的范围内。首先归一化是为了后面数据处理的方便，其次归一化能够保证程序运行时收敛加快。归一化的具体作用是归纳同意样本的统计分布性，归一化在0-1之间是统计的概率分布，归一化在某个区间上是统计的坐标分布，在机器学习算法的数据预处理阶段，归一化也是非常重要的步骤。其定义如下：

C++: void normalize(InputArray src, OutputArray dst, double alpha=1, double beta=0, int norm_type=NORM_L2, int dtype=-1, InputArray mask=noArray() )
1
参数解释： 
. InputArray src: 输入图像 
. OutputArray dst: 输出图像，尺寸大小和src相同 
. double alpha = 1: range normalization模式的最小值 
. double beta = 0: range normalization模式的最大值，不用于norm normalization(范数归一化)模式 
. int norm_type = NORM_L2: 归一化的类型，主要有 
。NORM_INF: 归一化数组的C-范数（绝对值的最大值） 
。NORM_L1: 归一化数组的L1-范数（绝对值的和） 
。NORM_L2: 归一化数组的L2-范数（欧几里得） 
。NORM_MINMAX: 数组的数值被平移或缩放到一个指定的范围，线性归一化，一般较常用。 
. int dtype = -1: 当该参数为负数时，输出数组的类型与输入数组的类型相同，否则输出数组与输入数组只是通道数相同，而depth = CV_MAT_DEPTH(dtype) 
. InputArray mask = noArray(): 操作掩膜版，用于指示函数是否仅仅对指定的元素进行操作。

示例程序：

#include
#include
#include
#include

using namespace std;
using namespace cv;

int main()
{
    Mat I = imread("lena.jpg", IMREAD_GRAYSCALE);       //读入图像灰度图

    //判断图像是否加载成功
    if (I.empty())
    {
        cout << "图像加载失败!" << endl;
        return -1;
    }
    else
        cout << "图像加载成功!" << endl << endl;

    Mat padded;                 //以0填充输入图像矩阵
    int m = getOptimalDFTSize(I.rows);
    int n = getOptimalDFTSize(I.cols);

    //填充输入图像I，输入矩阵为padded，上方和左方不做填充处理
    copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));

    Mat planes[] = { Mat_(padded), Mat::zeros(padded.size(),CV_32F) };
    Mat complexI;
    merge(planes, 2, complexI);     //将planes融合合并成一个多通道数组complexI

    dft(complexI, complexI);        //进行傅里叶变换

    //计算幅值，转换到对数尺度(logarithmic scale)
    //=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
    split(complexI, planes);        //planes[0] = Re(DFT(I),planes[1] = Im(DFT(I))
                                    //即planes[0]为实部,planes[1]为虚部
    magnitude(planes[0], planes[1], planes[0]);     //planes[0] = magnitude
    Mat magI = planes[0];

    magI += Scalar::all(1);
    log(magI, magI);                //转换到对数尺度(logarithmic scale)

    //如果有奇数行或列，则对频谱进行裁剪
    magI = magI(Rect(0, 0, magI.cols&-2, magI.rows&-2));

    //重新排列傅里叶图像中的象限，使得原点位于图像中心
    int cx = magI.cols / 2;
    int cy = magI.rows / 2;

    Mat q0(magI, Rect(0, 0, cx, cy));       //左上角图像划定ROI区域
    Mat q1(magI, Rect(cx, 0, cx, cy));      //右上角图像
    Mat q2(magI, Rect(0, cy, cx, cy));      //左下角图像
    Mat q3(magI, 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-1之间的浮点数将矩阵变换为可视的图像格式
    normalize(magI, magI, 0, 1, CV_MINMAX);

    imshow("输入图像", I);
    imshow("频谱图", magI);
    waitKey(0);

    return 0;
}
程序分析： 
1.图像填充：

Mat padded;
int m = getOptimalDFTSize(I.rows);
int n = getOptimalDFTSize(I.cols);
copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));

根据前面的理论介绍可以知道当图像尺寸为2、3、5的倍数时可以得到最快的处理速度，所以通过getOptimalDFTSize()函数获取最佳DFT变换尺寸，之后再结合copyMakeBorder()函数对图像进行扩充。

2.为实部和虚部分配存储空间

Mat planes[] = { Mat_(padded), Mat::zeros(padded.size(),CV_32F) };
Mat complexI;
merge(planes, 2, complexI);

傅里叶变换的结果是复数，这就意味着经过傅里叶变换每个图像值都会变成两个值，此外其频域(frequency domains)范围比空间域(spatial counterpart)范围大很多。我们通常以浮点型数据格式对结果进行存储。因此我们将输入图像转换为这种类型，通过另外的通道扩充图像。

3.傅里叶变换

dft(complexI, complexI);
1
傅里叶变换函数，对图像进行傅里叶变换。

4.将实数和复数的值转换为幅度值

split(complexI, planes);
magnitude(planes[0], planes[1], planes[0]); 
Mat magI = planes[0];
1
2
3
复数包含实部和虚部两个部分，傅里叶变换的结果是一个复数，傅里叶变换的幅度计算公式是： 
 
5.转换为对数尺度(Switch to a logarithmic scale)

magI += Scalar::all(1);
log(magI, magI);
1
2
之所以要进行对数转换是因为傅里叶变换后的结果对于在显示器显示来讲范围比较大，这样的话对于一些小的变化或者是高的变换值不能进行观察。因此高的变化值将会转变成白点，而较小的变化值则会变成黑点。为了能够获得可视化的效果，可以利用灰度值将我们的线性尺度(linear scale)转变为对数尺度(logarithmic scale)，其计算公式如下： 

6.剪切和象限变换

magI = magI(Rect(0, 0, magI.cols&-2, magI.rows&-2));
int cx = magI.cols / 2;
int cy = magI.rows / 2;

Mat q0(magI, Rect(0, 0, cx, cy));
Mat q1(magI, Rect(cx, 0, cx, cy));  
Mat q2(magI, Rect(0, cy, cx, cy));
Mat q3(magI, 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);

在进行傅里叶变换时，为了取得更快的计算效果，对图像进行了扩充，现在就需要对新增加的行列进行裁剪了。为了可视化的需要，我们同样需要对显示的结果图像像素进行调整，如果不进行调整，最后显示的结果是这样的： 
 
可以看到四周的角上时高频分量，现在我们通过重新调整象限将高频分量调整到图像正中间。

7.归一化

normalize(magI, magI, 0, 1, CV_MINMAX);
1
对结果进行归一化处理同样是处于可视化的目的。现在我们得到了幅度值，但是这仍然超出了0-1的显示范围。这就需要利用normalize()函数对数据进行归一化处理。

程序运行结果如图： 

————————————————
版权声明：本文为CSDN博主「梧桐栖鸦」的原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/keith_bb/article/details/53389819

//***********

傅里叶变换将图像分解成其正弦和余弦分量，它将图像由空域转换为时域。任何函数都可以近似的表示为无数正弦和余弦函数的和，傅里叶变换就是实现这一步的，数学上一个二维图像的傅里叶变换为：

正在上传…重新上传取消转存失败重新上传取消

公式中，f是图像在空域的值，F是频域的值。转换的结果是复数，但是不可能通过一个真实图像和一个复杂的图像或通过大小和相位图像去显示这样的一个图像。然而，在整个图像处理算法只对大小图像是感兴趣的，因为这包含了所有我们需要的图像几何结构的信息。

可通过以下几步显示一副傅里叶变换后的图像

1、将图像扩展到它的最佳尺寸，DFT(直接傅里叶变换)的性能依赖于图片的尺寸，当图像是2，3，5的倍数时往往是最快的。因此，为了达到最优性能通常采用垫边界值的方法，得到一个最佳的尺寸。

2、为傅立叶变换结果的实部和虚部分配存储空间。傅里叶变换的结果是一个复数，这意味着每幅图的结果都有一个实部和虚部，此外，频域范围远远大于它对应的空间范围。因此，我们这些通常至少以一个浮点数格式存储这些数值。因此，我们会将我们的输入图像转换为这种类型并且扩展它与另一通道存放复数值

3、进行傅里叶变换。

4、将复数转换为幅值，DFT的幅值由以下公式得出：正在上传…重新上传取消转存失败重新上传取消

5、切换到对数刻度。对图像进行对数尺度的缩放，结果证明，傅立叶系数矩阵的动态范围太大，无法显示在屏幕上，我们无法通过这样去观察一些小的和高的变化值。因此那些高的数值将转化成白点而小的数值会变成黑点，使用灰度值进行可视化，我们可以将线性刻度转换为对数刻度，以便于观察。

转存失败重新上传取消转存失败重新上传取消

6、剪切和重分布幅度图象，第一步我们扩展了图像，这里我们去掉扩展的那部分值，基于可视化的目的，我们还可以重新排列结果的象限，使原点(0，0)对应于与图像中心

7、归一化。目前得到的幅值图像仍然太大，超出了显示的范围，归一化这范围内的值，可以进一步达到可视化的目的

实现程序

?

正在上传…重新上传取消转存失败重新上传取消

傅里叶变换后的图片

正在上传…重新上传取消转存失败重新上传取消

我是做Tracking 的，对于速度要求非常高。发现傅里叶变换能够使用。

于是学习之。

核心： 最根本的一点就是将时域内的信号转移到频域里面。这样时域里的卷积能够转换为频域内的乘积！

      在分析图像信号的频率特性时，对于一幅图像，直流分量表示预想的平均灰度。低频分量代表了大面积背景区域和缓慢变化部分，高频部分代表了它的边缘，细节，跳跃部分以及颗粒噪声.  因此，我们能够做对应的锐化和模糊的处理：提出当中的高频分量做傅里叶逆变换得到的就是锐化的结果。

提出当中的低频分量做傅里叶逆变换得到的就是模糊的结果。

           最不能理解的应该是：截取频域图中的不论什么一个区域相应的都是原来的整张图的区域。而不是相应的局部。

 由于频域内的各个点都反映的是整张图的一个状态。

我们能够用时间和频率来理解：当你走完一段单位路程的时候。如果你花了100秒，那么你的频率就是0.01HZ。

这个0.01HZ显然体现的是一个总体的结果。而不是局部。

我们再由公式来看：

能够非常明显的知道频域内的每个点的值都是由整个图像求出来的。当然以上得出的结果，我们一般仅仅关注幅值频谱图。

也就是说真正起作用的就是前面的那个cos x而已. 于是我们能够知道。在整个范围内（0<k <N, 0<l <N）,低频分量集中于四个角。

且其它地方的值仅仅可能比这个小。

在原点的傅里叶变换即等于图像的平均灰度级。由于 在原点处经常为零，F(0，0)有时称做 频率谱的直流成分。 

使用：

当图像的尺寸是2，3，5的整数倍时，计算速度最快。因此opencv里面有一个函数：

它能够使得图片的尺寸能够满足这个要求。

可是这样就须要对原来的图像进行大小的处理，因此使用函数：CopyMakeBorder复制图像而且制作边界。

（处理边界卷积）

将原始的图像I 扩充为理想的大小放在padded里面。

接下来我们须要给计算出来的结果分配空间：

然后便能够进行傅里叶变换了：

得到的结果有两部分。实数部分和虚数部分，你能够分别对这两部分进行操作：

当然还能够进行：归一化：

另外重要的一个应用是： convolveDFT。

 当中的 *代表的是 卷积。我认为这也是我们进行离散傅里叶变换的目的。

使得计算的速度大大的添加。

先来说一下卷积在图像中的意义：

如果图像f(x),模板是g(x),然后将模版g(x)在模版中移动,每到一个位置,就把f(x)与g(x)的定义域相交的元素进行乘积而且求和,得出新的图像一点,就是被卷积后的图像. 模版又称为卷积核.卷积核做一个矩阵的形状.（当然边缘点可能须要特殊的处理，同一时候这个操作和滤波也非常像，或许就是一回事）。

当中：

C 为什么是这种勒？想想一个特殊的样例就知道了：当A，B尺寸相等的时候，这个时候的高斯滤波得到的也就是中心点的那一个值（卷积核滤波的区别在于须要绕中心180度旋转）。