【opencv/core module】(七)Discrete Fourier Transform

说在前面

• opencv版本：4.0.1
• 操作系统：win10
• vs版本：2017
• 官方文档：Discrete Fourier Transform
• 其他说明：自学，记录，demo

Theory

$$F(k,l)=\sum _{i=0}^{N-1}\sum _{j=0}^{N-1}f(i,j)e^{-i2\pi (\frac{ki}{N}+\frac{lj}{N})}$$
$$e^{ix}=\cos x+i\sin x$$
离散傅里叶变换原理咱先缓缓 (可能后续补充) ，/(ㄒoㄒ)/~~
对一张图象进行傅里叶变换后，可以得到real/magnitude image&complex/phase image

• real/magnitude image
  包含了所有我们感兴趣的几何结构信息
• complex/phase image
  如果要进行逆傅里叶变换，该图像也需要

SourceCode

• Expand the image to an optimal size(将图像扩张至一个最佳的大小)

  • 离散傅里叶变换的性能取决于图像的大小。性能最佳时的图像大小往往是2、3、5的倍数。所以我们需要将图像的大小转成满足这种条件的大小。

  • 函数 getOptimalDFTSize(int vecsize) 可以获得满足下述关系的最小T
    $$T=2^p\cdot 3^q\cdot 5^r\ge vecsize$$
    其中$p、q、r$为非负整数。
    举个栗子：
    $$2^1\cdot 3^3\cdot 5^1=270>257$$
    

  • 函数 copyMakeBorder 给图像加边界(border)

    void cv::copyMakeBorder 
    (
    InputArray src, //输入图像
    OutputArray dst, //输出图像
    int top, //上边框宽度，top个像素
    int bottom, //下边框宽度
    int left, //左边框宽度
    int right, //右边框宽度
    int borderType, //边框类型
    const Scalar & value = Scalar() 
    //若边框类型为BORDER_CONSTANT，则代表边框的值
    )
    

    

    • code

    	Mat padded;                            //expand input image to optimal size
    	int m = getOptimalDFTSize( I.rows );
    	int n = getOptimalDFTSize( I.cols ); // on the border add zero values
    	copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols,
    									 BORDER_CONSTANT, Scalar::all(0));
    

    • 扩张后图像
      
    • 关于borderType

type	含义	图片
BORDER_CONSTANT	连续定值	见附
BORDER_REPLICATE	复制边缘像素	见附
BORDER_REFLECT	对称 (fedcba|abcdefgh|hgfedcb)	见附
BORDER_WRAP	重复原图	见附
BORDER_REFLECT_101	对称101 (fedcb|abcdefgh|gfedcb)
BORDER_TRANSPARENT	透明
BORDER_REFLECT101	对称101
BORDER_DEFAULT	默认(对称101)
BORDER_ISOLATED	(这个不懂)

• Make place for both the complex and the real values(扩张成两通道图像)

  • 使用函数merge()混合Mat，举个栗子
    矩阵A：
    $$\left[\begin{array}{ccc}1& 2& 3\\ 2& 3& 4\\ 4& 5& 6\end{array}\right]$$
    矩阵B：
    $$\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$$
    merge A B：
    $$\left[\begin{array}{cccccc}1& 0& 2& 0& 3& 0\\ 2& 0& 3& 0& 4& 0\\ 4& 0& 5& 0& 6& 0\end{array}\right]$$
  • merge函数

    void cv::merge 
    (
    const Mat * mv, //输入矩阵数组
    size_t count, //要处理矩阵数组中的前几个矩阵
    OutputArray dst //输出矩阵
    )
    

  • code

    //因为dft后的值比较大，因此我们使用浮点数
    Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
    Mat complexI;
    merge(planes, 2, complexI);// Add to the expanded another plane with zeros
    

• Make the Discrete Fourier Transform(进行离散傅里叶变换)

  • dft函数

    void cv::dft 
    (
    InputArray src, 
    OutputArray dst, 
    int flags = 0, 
    int nonzeroRows = 0 
    )
    

  • code

    dft(complexI, complexI); // this way the result may fit in the source matrix
    

• Transform the real and complex values to magnitude(求模)

  • theory
    复数求模
    $$M=\sqrt[2]{{Re(DFT(I))}^2+{Im(DFT(I))}^2}$$
  1. 使用dft函数后得到的图像(双通道)为实数图像+虚数图像，分别占用一个通道
  2. 使用split函数将实数图和虚数图分离开(即merge的逆操作)
  3. 使用magnitude函数求模即可
  • code

    void cv::split 
    (
    const Mat & src, //多通道的源图像
    Mat * mvbegin //与源图像通道数匹配的Mat数组
    )
    
    void cv::magnitude 
    (
    InputArray x, //向量在x方向的浮点数组
    InputArray y, //向量在y方向的浮点数组，与x数组大小相同
    OutputArray magnitude //输出数组
    )
    
    split(complexI, planes);   // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
    magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude
    Mat magI = planes[0];
    

• Switch to a logarithmic scale(缩放)

  • 由于傅里叶变换后的值过大(本例中超过255)，直接使用imshow的话显示出来是白色，所以我们取个自然对数，便于显示。
    $$M_1=\log (1+M)$$
  • 下图为取对数前，[0,0]像素的值
    
  • code

    magI += Scalar::all(1);                    // switch to logarithmic scale
    log(magI, magI);
    

• Crop and rearrange(剪裁&调整位置)

  • 由于在第一步我们将图像大小扩大了，因此我们需要去除这部分的值(不是很懂这部分的原理，有大佬吗？验证感觉并未实现这功能，见如下结果，宽度未变)
    
  • rearrange

    // crop the spectrum, if it has an odd number of rows or columns
    magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));
    //11 & -2 = 10
    //这里应该和下面的 magI.cols/2 有关，保证可以得到图像的中心
    //取比原值小的最接近的偶数
    
    // rearrange the quadrants of Fourier image  
    //so that the origin is at the image center
    int cx = magI.cols/2;
    int cy = magI.rows/2;
    
    //取ROI，之前有提到过，未发生数据的拷贝
    Mat q0(magI, Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant
    Mat q1(magI, Rect(cx, 0, cx, cy));  // Top-Right
    Mat q2(magI, Rect(0, cy, cx, cy));  // Bottom-Left
    Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right
    
    //使用copyTo函数将数据进行拷贝
    Mat tmp;                           // swap quadrants (Top-Left with Bottom-Right)
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);
    
    q1.copyTo(tmp);                    // swap quadrant (Top-Right with Bottom-Left)
    q2.copyTo(q1);
    tmp.copyTo(q2);
    

调换前：

调换后：

• Normalize(归一化)

  • 浮点类型的Mat可视部分貌似只有$[0,1]$，因此我们需要使用$[0,1]$来代替原来的值。
  • 原理也比较简单，找出最大值MAX、最小值MIN$$\frac{Value-MIN}{MAX-MIN}$$

    normalize(magI, magI, 0, 1, NORM_MINMAX); 
    // Transform the matrix with float values into a
    // viewable image form (float between values 0 and 1).
    

• Code

#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/imgcodecs.hpp"
#include "opencv2/highgui.hpp"

#include <iostream>

using namespace cv;
using namespace std;

int main()
{
	const char* filename = "gamma.jpg";

	Mat I = imread(filename, IMREAD_GRAYSCALE);
	if (I.empty()) {
		cout << "Error opening image" << endl;
		return -1;
	}

	Mat Mfloat = Mat_ <float>(I);
	imshow("float", Mfloat);

	Mat padded;                            //expand input image to optimal size
	int m = getOptimalDFTSize(I.rows);
	int n = getOptimalDFTSize(I.cols); // on the border add zero values
	/*cout << "input image rows:" << I.rows << endl;
	cout << "input image cols:" << I.cols << endl;

	cout << "optimal rows:" << m << endl;
	cout << "optimal cols:" << n << endl;*/

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

	imshow("padded", padded);


	Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
	Mat complexI;
	merge(planes, 2, complexI);     // Add to the expanded another plane with zeros

	dft(complexI, complexI);        // this way the result may fit in the source matrix

								// compute the magnitude and switch to 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))
	magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude
	Mat magI = planes[0];

	magI += Scalar::all(1);                    // switch to logarithmic scale
	log(magI, magI);

	cout << "after log Mat[0,0]:" << magI.at<Vec2f>(0, 0) << endl;

	// crop the spectrum, if it has an odd number of rows or columns
	magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));
	/*cout << "magI rows:" << magI.rows << endl;
	cout << "magI cols:" << magI.cols << endl;
	cout << "magI rows & -2:" << (magI.rows & -2) << endl;
	cout << "magI cols & -2:" << (magI.cols & -2) << endl;*/

	// rearrange the quadrants of Fourier image  
	//so that the origin is at the image center
	int cx = magI.cols / 2;
	int cy = magI.rows / 2;

	Mat q0(magI, Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant
	Mat q1(magI, Rect(cx, 0, cx, cy));  // Top-Right
	Mat q2(magI, Rect(0, cy, cx, cy));  // Bottom-Left
	Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right

	Mat tmp;                           // swap quadrants (Top-Left with Bottom-Right)
	q0.copyTo(tmp);
	q3.copyTo(q0);
	tmp.copyTo(q3);

	q1.copyTo(tmp);                    // swap quadrant (Top-Right with Bottom-Left)
	q2.copyTo(q1);
	tmp.copyTo(q2);

	normalize(magI, magI, 0, 1, NORM_MINMAX); 
	// Transform the matrix with float values int
	// viewable image form (float between values 0 and 1).
	cout << "after normal Mat[0,0]:" << magI.at<Vec2f>(0, 0) << endl;

	imshow("Input Image", I);    // Show the result
	imshow("spectrum magnitude", magI);
	cout << magI.depth() << endl;
	waitKey();

	return 0;
}

Result

• 原图1
  
• 原图2
  
• 输出结果1
  
• 输出结果2
  可以看到这个和原图2一样，倾斜了一点
  

附

• CONSTANT
  
• REPLICATE
  
• REFLECT
  
• WRAP
  

---

END-(CSDN)2019.6.30