上一个教程 : 改变图像的对比度和亮度
下一个教程 : 使用 XML 和 YAML 文件输入和输出文件
原作者 | Bernát Gábor |
---|---|
兼容性 | OpenCV >= 3.0 |
目标
我们将寻求以下问题的答案:
- 什么是傅立叶变换,为什么要使用它?
- 如何在 OpenCV 中进行傅立叶变换?
- 函数的使用,如:copyMakeBorder() 、merge() 、dft() 、getOptimalDFTSize() 、log() 和 normalize()。
源代码
- C++
您可以从此处下载,或在 OpenCV 源代码库的samples/cpp/tutorial_code/core/discrete_fourier_transform/discrete_fourier_transform.cpp
中找到它。
下面是 dft() 的使用示例:
#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/imgcodecs.hpp"
#include "opencv2/highgui.hpp"
#include <iostream>
using namespace cv.NET
using namespace std;
static void help(char ** argv)
{
cout << endl
<< "This program demonstrated the use of the discrete Fourier transform (DFT). " << endl
<< "The dft of an image is taken and it's power spectrum is displayed." << endl << endl
<< "Usage:" << endl
<< argv[0] << " [image_name -- 默认为 lena.jpg]" << endl << endl;
}
int main(int argc, char ** argv)
{
help(argv);
const char* filename = argc >=2 ? argv[1] : "lena.jpg";
Mat I = imread( samples::findFile( filename ), IMREAD_GRAYSCALE);
if( I.empty()){
cout << "Error opening image" << endl;
return EXIT_FAILURE;
}
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));
Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexI;
merge(planes, 2, complexI); // 在扩展后的平面上再添加一个零平面
dft(complexI, complexI); // 这样得到的结果才可能适合源矩阵
// 计算幅值并转换为对数标度
// => 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); // 切换到对数标度
log(magI, magI);
// 如果光谱的行数或列数为奇数,则裁剪光谱
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)); // Bottom-Right
Mat tmp; // 交换象限(左上角和右下角)
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp); // 交换象限(右上角和左下角)
q2.copyTo(q1);
tmp.copyTo(q2);
normalize(magI,magI,0,1,NORM_MINMAX);//将带有浮点数值的矩阵转换为可查看的图像形式(浮点数值介于 0 和 1 之间)。
// 可查看的图像形式(数值在 0 和 1 之间的浮点)。
imshow("Input Image" , I ); // 显示结果
imshow("spectrum magnitude", magI);
waitKey();
return EXIT_SUCCESS;
}
- Java
您可以从此处下载,或在 OpenCV 源代码库的samples/java/tutorial_code/core/discrete_fourier_transform/DiscreteFourierTransform.java
中找到它。
下面是 dft() 的使用示例:
import org.opencv.core.*;
import org.opencv.highgui.HighGui;
import org.opencv.imgcodecs.Imgcodecs;
import java.util.List;
import java.util.*;
class DiscreteFourierTransformRun{
private void help() {
System.out.println("" +
"This program demonstrated the use of the discrete Fourier transform (DFT). \n" +
"The dft of an image is taken and it's power spectrum is displayed.\n" +
"Usage:\n" +
"./DiscreteFourierTransform [image_name -- 默认为 ../data/lena.jpg]");
}
public void run(String[] args){
help();
String filename = ((args.length > 0) ? args[0] : "../data/lena.jpg");
Mat I = Imgcodecs.imread(filename, Imgcodecs.IMREAD_GRAYSCALE);
if( I.empty() ) {
System.out.println("Error opening image");
System.exit(-1);
}
Mat padded = new Mat(); // 将输入图像扩展到最佳尺寸
int m = Core.getOptimalDFTSize( I.rows() );
int n = Core.getOptimalDFTSize( I.cols() ); // 在边框上添加零值
Core.copyMakeBorder(I, padded, 0, m - I.rows(), 0, n - I.cols(), Core.BORDER_CONSTANT, Scalar.all(0));
List<Mat> planes = new ArrayList<Mat>();
PADEDED.CONVERT TO(PADEDED, CvType.CV_32F);
planes.add(padded);
planes.add(Mat.zeros(padded.size(), CvType.CV_32F));
Mat complexI = new Mat();
Core.merge(planes, complexI); // 在扩展后的平面中添加另一个零平面
Core.dft(complexI, complexI); // 这样得到的结果才可能适合源矩阵
// 计算幅值并转换为对数标度
// => log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
Core.split(complexI, planes); // planes.get(0) = Re(DFT(I)
// planes.get(1) = Im(DFT(I))
Core.magnitude(planes.get(0), planes.get(1), planes.get(0));// planes.get(0) = magnitude
Mat magI = planes.get(0);
Mat matOfOnes = Mat.ones(magI.size(), magI.type());
Core.add(matOfOnes, magI, magI); // 切换到对数标度
Core.log(magI, magI);
// 如果光谱的行数或列数为奇数,则裁剪光谱
magI = magI.submat(new Rect(0, 0, magI.cols() & -2, magI.rows() & -2));
// 重新排列傅立叶图像的象限,使原点位于图像中心
int cx = magI.cols()/2;
int cy = magI.rows()/2;
Mat q0 = new Mat(magI, new Rect(0, 0, cx, cy)); // Top-Left - 每个象限创建一个 ROI
Mat q1 = new Mat(magI, new Rect(cx, 0, cx, cy)); // 右上角
Mat q2 = new Mat(magI, new Rect(0, cy, cx, cy)); // 左下角
Mat q3 = new Mat(magI, new Rect(cx, cy, cx, cy)); // Bottom-Right
Mat tmp = new Mat(); // 交换象限(左上角和右下角)
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp); // 交换象限(右上角和左下角)
q2.copyTo(q1);
tmp.copyTo(q2);
magI.convertTo(magI, CvType.CV_8UC1);
Core.normalize(magI,magI,0,255,Core.NORM_MINMAX,CvType.CV_8UC1);//将带有浮点数值的矩阵转换为可视图像形式(CV_8UC1
// 转换为可查看的图像形式(浮点值介于
// 值 0 和 255 之间的浮点)。
HighGui.imshow("Input Image" , I ); // 显示结果
HighGui.imshow("Spectrum Magnitude", magI );
HighGui.waitKey();
System.exit(0);
}
}
public class DiscreteFourierTransform {
public static void main(String[] args) {
// 加载本地库。
System.loadLibrary(Core.NATIVE_LIBRARY_NAME);
new DiscreteFourierTransformRun().run(args);
}
}
- Python
您可以从此处下载,或在 OpenCV 源代码库的samples/python/tutorial_code/core/discrete_fourier_transform/discrete_fourier_transform.py
中找到它。
下面是 dft() 的用法示例:
from __future__ import print_function
import sys
import cv2 as cv
import numpy as np
def print_help():
print('''
This program demonstrated the use of the discrete Fourier transform (DFT).
The dft of an image is taken and it's power spectrum is displayed.
Usage:
discrete_fourier_transform.py [image_name -- default lena.jpg]''')
def main(argv):
print_help()
filename = argv[0] if len(argv) > 0 else 'lena.jpg'.
I = cv.imread(cv.samples.findFile(filename), cv.IMREAD_GRAYSCALE)
if I 为空:
print('Error opening image')
return -1
rows,cols = I.shape
m = cv.getOptimalDFTSize( rows )
n = cv.getOptimalDFTSize( cols )
padded = cv.copyMakeBorder(I, 0, m - rows, 0, n - cols, cv.BORDER_CONSTANT, value=[0, 0, 0])
planes = [np.float32(padded), np.zeros(padded.shape, np.float32)]
complexI = cv.merge(planes) # 在扩展后的平面上再添加一个零平面
cv.dft(complexI, complexI) # 这样,结果就可能适合源矩阵了
cv.split(complexI, planes) # planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
cv.magnitude(planes[0], planes[1], planes[0])# planes[0] = magnitude
magI = planes[0]
matOfOnes = np.ones(magI.shape, dtype=magI.dtype)
cv.add(matOfOnes, magI, magI) # 切换到对数标度
cv.log(magI, magI)
magI_rows, magI_cols = magI.shape
# 如果光谱的行数或列数为奇数,则裁剪光谱
magI = magI[0:(magI_rows & -2),0:(magI_cols & -2)] cx = int(magI_rows & -2)
cx = int(magI_rows/2)
cy = int(magI_cols/2)
q0 = magI[0:cx, 0:cy] # 左上角 - 每个象限创建一个 ROI
q1 = magI[cx:cx+cx, 0:cy] # 右上角
q2 = magI[0:cx, cy:cy+cy] # 左下角
q3 = magI[cx:cx+cx, cy:cy+cy] # 右下方
tmp = np.copy(q0) # 交换象限(左上角和右下角)
magI[0:cx, 0:cy] = q3
magI[cx:cx + cx, cy:cy + cy] = tmp
tmp = np.copy(q1) # 交换象限(从右上角到左下角)
magI[cx:cx + cx, 0:cy] = q2
magI[0:cx, cy:cy + cy] = tmp
cv.normalize(magI, magI, 0, 1, cv.NORM_MINMAX) # 将包含浮点数值的矩阵转换成一个
cv.imshow("Input Image" , I ) # 显示结果
cv.imshow("spectrum magnitude", magI)
cv.waitKey()
if __name__ == "__main__":
main(sys.argv[1:])
解释
傅立叶变换将图像分解为正弦和余弦成分。换句话说,它将图像从空间域转换到频率域。其原理是,任何函数都可以用无限个正弦和余弦函数之和来近似。傅立叶变换就是实现这一目的的一种方法。从数学上讲,二维图像的傅里叶变换是:
F
(
k
,
l
)
=
∑
i
=
0
N
−
1
∑
j
=
0
N
−
1
f
(
i
,
j
)
e
−
i
2
π
(
k
i
N
+
l
j
N
)
F\left( k,l \right) =\sum_{i=0}^{N-1}{\sum_{j=0}^{N-1}{f\left( i,j \right) e^{-i2\pi \left( \frac{ki}{N}+\frac{lj}{N} \right)}}}
F(k,l)=i=0∑N−1j=0∑N−1f(i,j)e−i2π(Nki+Nlj)
e
i
x
=
cos
x
+
i
sin
x
e^{ix}=\cos x+i\sin x
eix=cosx+isinx
其中,f 是空间域的图像值,F 是频率域的图像值。变换的结果是复数。可以通过实数图像和复数图像,或通过幅值图像和相位图像来显示。不过,在整个图像处理算法中,只有幅值图像才是有意义的,因为它包含了我们所需的有关图像几何结构的所有信息。不过,如果您打算以这些形式对图像进行一些修改,然后需要对其进行重新变换,则需要保留这两种形式的图像。
在本示例中,我将演示如何计算和显示傅立叶变换的幅值图像。数字图像是离散的。这意味着它们可以从给定的域值中取值。例如,在基本灰度图像中,数值通常介于 0 和 255 之间。因此,傅立叶变换也需要是离散类型的,从而产生离散傅立叶变换 (DFT)。当您需要从几何角度确定图像结构时,就会用到这种方法。以下是操作步骤(在输入图像 I 为灰度图像的情况下):
将图像扩展到最佳尺寸
DFT 的性能取决于图像的大小。当图像大小是 2、3 和 5 的倍数时,DFT 的速度最快。因此,要获得最佳性能,通常最好在图像中填充边框值,以获得具有这些特征的图像大小。getOptimalDFTSize() 会返回这种最佳尺寸,我们可以使用 copyMakeBorder() 函数来扩展图像的边框(添加的像素初始化为零):
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));
同时放置复数值和实数值
傅立叶变换的结果是复值。这意味着每个图像值的结果都是两个图像值(每个分量一个)。此外,频域范围远大于空间范围。因此,我们通常至少以浮点格式存储这些值。因此,我们要将输入图像转换为这种类型,并用另一个通道来扩展它,以保存复数值:
Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexI;
merge(planes, 2, complexI); // 在扩展后的另一个平面上添加零点
进行离散傅立叶变换
可以进行就地计算(输入与输出相同):
dft(complexI, complexI); // 这样计算结果就可以放入源矩阵中了
将实数和复数值转换为幅度
复数有实部(Re)和虚部(Im)。DFT 的结果就是复数。DFT 的幅值为
M
=
Re
(
D
F
T
(
I
)
)
2
+
Im
(
D
F
T
(
I
)
)
2
M=\sqrt{\text{Re}\left( DFT\left( I \right) \right) ^2+\text{Im}\left( DFT\left( I \right) \right) ^2}
M=Re(DFT(I))2+Im(DFT(I))2
翻译为 OpenCV 代码:
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];
切换到对数标度
事实证明,傅立叶系数的动态范围太大,无法在屏幕上显示。我们有一些变化较小的值和一些变化较大的值,无法像这样观察。因此,数值大的会显示为白点,数值小的会显示为黑点。为了将灰度值用于可视化,我们可以将线性刻度转换为对数刻度:
M 1 = l o g ( 1 + M ) M_1=log\left( 1+M \right) M1=log(1+M)
转换为 OpenCV 代码:
magI += Scalar::all(1); // 切换到对数标度
log(magI, magI);
裁剪并重新排列
还记得在第一步,我们展开了图像吗?现在是时候丢掉新引入的值了。出于可视化目的,我们还可以重新排列结果的象限,使原点(零,零)与图像中心相对应。
// 如果光谱的行数或列数为奇数,则裁剪光谱
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)); // Bottom-Right
Mat tmp; // 交换象限(左上角和右下角)
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp); // 交换象限(右上角和左下角)
q2.copyTo(q1);
tmp.copyTo(q2);
规范化
为了实现可视化,我们再次这样做。现在我们有了幅度,但仍超出了图像显示范围(0 到 1)。我们使用 cv::normalize() 函数将数值归一化到这个范围。
normalize(magI, magI, 0, 1, NORM_MINMAX); // 将带有浮点数值的矩阵转换为可视图像形式(浮点数值介于 0 和 1 之间)。
// 可查看的图像形式(数值在 0 和 1 之间的浮点)。
结果
一种应用思路是确定图像中存在的几何方向。例如,让我们找出文字是否水平?通过观察一些文字,你会发现文字的线条有点像水平线,而字母则有点像垂直线。在傅立叶变换中也可以看到文本片段的这两个主要组成部分。让我们用水平图像和旋转图像来表示文字。
如果是水平文本,则
如果是旋转文字
您可以看到,频域中影响最大的分量(幅度图像上最亮的点)随着图像上物体的几何旋转而变化。由此,我们可以计算出偏移量,并进行图像旋转,以纠正最终的错位。