OpenCV4

Goals

  • 了解傅里叶变换
  • 了解如何使用 OpenCV 完成傅里叶变换
  • 了解 copyMakeBordermergedftgetOptimalDFTSizelognormalize

Code

C++:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
int (int argc, char *argv[])
{
std::string imagePath = "./car.jpg";
cv::Mat src = cv::imread(imagePath, cv::IMREAD_GRAYSCALE);
if (src.empty())
{
std::cout << "Can't load image." << std::endl;
return -1;
}


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

// 构建复平面
cv::Mat planes[] = { cv::Mat_<float>(padded), cv::Mat::zeros(padded.size(), CV_32F) };
cv::Mat complexImg;
cv::merge(planes, 2, complexImg);

// 离散傅里叶变换
cv::dft(complexImg, complexImg); // 支持 in-place 计算

// 计算幅值,并转换为log(1 + magnitude), magnitude = sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)
cv::split(complexImg, planes); // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))
cv::magnitude(planes[0], planes[1], planes[0]);
cv::Mat magnitudeImg = planes[0];
magnitudeImg += cv::Scalar::all(1);
cv::log(magnitudeImg, magnitudeImg);

magnitudeImg = magnitudeImg(cv::Rect(0, 0, magnitudeImg.cols & -2, magnitudeImg.rows & -2));

// 重新排列图像的象限,以便原点在图像中心
int cx = magnitudeImg.cols / 2;
int cy = magnitudeImg.rows / 2;
cv::Mat q0(magnitudeImg, cv::Rect(0, 0, cx, cy));
cv::Mat q1(magnitudeImg, cv::Rect(cx, 0, cx, cy));
cv::Mat q2(magnitudeImg, cv::Rect(0, cy, cx, cy));
cv::Mat q3(magnitudeImg, cv::Rect(cx, cy, cx, cy));

cv::Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);

cv::normalize(magnitudeImg, magnitudeImg, 0, 1, cv::NORM_MINMAX);

cv::imshow("src", src);
cv::imshow("magnitude", magnitudeImg);
cv::waitKey();
return 0;
}

Python:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
def ():
img_path = "./car.jpg"
src = cv2.imread(img_path, cv2.IMREAD_GRAYSCALE)
if src is None:
print("Can't load image!")
return -1
rows, cols = src.shape
m = cv2.getOptimalDFTSize(rows)
n = cv2.getOptimalDFTSize(cols)
# 0填充至DFT最佳大小,2/3/5的倍数
padded = cv2.copyMakeBorder(src, 0, m - rows, 0, n - cols,
cv2.BORDER_CONSTANT, value=[0,0,0])
planes = [np.float32(padded), np.zeros(padded.shape, np.float32)]
complex_img = cv2.merge(planes) # 合成双通道以容纳复值图像

complex_img = cv2.dft(complex_img)
planes = cv2.split(complex_img) # planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))
magnitude_img = cv2.magnitude(planes[0], planes[1]) # magnitude = sqrt(Re^2 + Im^2)

# 用log(1+mag)压缩,便于显示
magnitude_img += 1
magnitude_img = cv2.log(magnitude_img)

mag_rows, mag_cols = magnitude_img.shape
# crop
magnitude_img = magnitude_img[0:(mag_rows & -2), 0:(mag_cols & -2)]
cx = int(mag_rows/2)
cy = int(mag_cols/2)

# 重新排列使得原点在图像中心
q0 = magnitude_img[0:cx, 0:cy]
q1 = magnitude_img[cx:cx+cx, 0:cy]
q2 = magnitude_img[0:cx, cy:cy+cy]
q3 = magnitude_img[cx:cx+cx, cy:cy+cy]

tmp = q0.copy()
magnitude_img[0:cx, 0:cy] = q3
magnitude_img[cx:cx+cx, cy:cy+cy] = tmp
tmp = q1.copy()
magnitude_img[cx:cx+cx, 0:cy] = q2
magnitude_img[0:cx, cy:cy+cy] = tmp

cv2.normalize(magnitude_img, magnitude_img, 0, 1, cv2.NORM_MINMAX)

cv2.imshow("Src", src)
cv2.imshow("Spectrum magnitude", magnitude_img)
cv2.waitKey()


if __name__ == "__main__":
main()

Explanation

上述代码首先读取图像的灰度图,傅里叶变换的性能取决于图像大小,当其大小为2,3,5的倍数时执行效率更高,所以根据图像的大小使用 getOptimalDFTSize 来计算最佳大小,然后用0填充图像至相应大小,因为傅里叶变换结果是复数,所以结果包含两个图像,并且频域范围比空间域大,所以一般使用 float 类型,构建一个包含2通道的 float 图像来存储傅里叶变换结果,接着使用 dft 来完成离散傅里叶变换,使用 split 来将结果分离出实部和虚部两个图像,使用 magnitude 来计算幅值,M = sqrt(Re(I)^2 + Im(I)^2),再使用 log(1+M) 来将图像压缩,便于显示,最后重新排列图像使得原点处于图像中心,因为 OpenCV 规定的图像中心在左上角,并使用 normalize 对图像进行归一化。

更多傅里叶相关的知识,可百度或谷歌。

Reference