C++版本的OpenCV实现二维图像的卷积定理(通过傅里叶变换实现二维图像的卷积过程，附代码！！)

前言

工作中用到许多卷积过程，需要转成C++代码的实现，使用OpenCV库自带的二维卷积过程所耗费的时间比较久，为了提升代码的运行效率可以考虑使用卷积定理实现二维图像的卷积过程。

一、卷积定理简单介绍

卷积定理是傅立叶变换满足的一个重要性质。卷积定理指出，函数卷积的傅立叶变换是函数傅立叶变换的乘积。具体分为时域卷积定理和频域卷积定理，时域卷积定理即时域内的卷积对应频域内的乘积；频域卷积定理即频域内的卷积对应时域内的乘积，两者具有对偶关系。

二、不同卷积过程对应的傅里叶变换过程

1、“Same”卷积

假设原始图像的大小为mxm，卷积核的大小为nxn。此时进行same卷积，则卷积后生成的图像大小为mxm，此时进行卷积是需要对原始图像进行padding操作，需要对原始图像周围进行n-1个补零操作，此时被卷积图像大小padding为(m+n-1)x(m+n-1)，然后进行卷积操作，进行傅里叶变化是需要先将被卷积图像以及卷积核的大小都padding为(m+n-1)x(m+n-1)，因为要进行傅里叶变换后的乘积操作，因此被卷积图像以及卷积核的大小需要相等。经过卷积定理的操作之后，将生成的(m+n-1)x(m+n-1)大小的图像按照padding的逆操作进行裁剪，得到mxm的图像即为same卷积得到的卷积图像。
Python代码验证

import numpy as np
from scipy import signal

# 原始图像 f(x)
gray = np.uint16(np.random.randint(100, size=(7, 7)))
# 卷积核 g(x)
kenel = np.ones((3, 3))/9
# ----- Conv = f(x)*g(x) ----- #
# f(x)*g(x)
Conv = signal.convolve2d(gray, kenel, mode='same') # 使用full卷积类型，得到(M+N-1)X(M+N-1)大小

#--------- ifft{ F(f(x))·F(g(x)) } ---------#
# 傅里叶变换前进行 padding 填充，图像和卷积核都补零到 (M+N-1)x(M+N-1)大小。
img_pad = np.pad(gray, ((1, 1), (1, 1)), 'constant')
kenel_pad = np.pad(kenel, ((3, 3),(3, 3)), 'constant')

# F(f(x))
img_fft = np.fft.fftshift(np.fft.fft2(img_pad))
# F(g(x))
kenel_fft = np.fft.fftshift(np.fft.fft2(kenel_pad))
# ifft( F(f(x))·F(g(x)) )
FFT = np.fft.ifftshift(np.fft.ifft2(np.fft.fftshift(img_fft*kenel_fft)))

#--------- 打印结果 ---------#
print(" f(x) ↓")
print(gray)
print(" g(x) ↓")
print(kenel)

print("\n\n f(x)*g(x) ↓")
print(np.uint8(Conv))
print("\n ifft[F·G] ↓")
print(np.uint8(np.abs(FFT)))
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

结果展示
上述代码的结果如下图所示，可见卷积定理得到的结果经过裁剪后的红框中的内容与same卷积得到的结果一致。

2、“Full”卷积

假设原始图像的大小为mxm，卷积核的大小为nxn。此时进行full卷积，则卷积后生成的图像大小为(m+n-1)x(m+n-1)，此时进行卷积是需要对原始图像进行padding操作，需要对原始图像周围进行2n-2个补零操作，此时被卷积图像大小padding为(m+2n-2)x(m+2n-12)，然后进行卷积操作，进行傅里叶变化是需要先将被卷积图像以及卷积核的大小都padding为(m+n-1)x(m+n-1)，因为要进行傅里叶变换后的乘积操作，因此被卷积图像以及卷积核的大小需要相等。经过卷积定理的操作之后，得到(m+n-1)x(m+n-1)的图像即为full卷积得到的卷积图像。
Python代码验证

import numpy as np
from scipy import signal

# 原始图像 f(x)
gray = np.uint16(np.random.randint(100, size=(9, 9)))
# 卷积核 g(x)
kenel = np.ones((5, 5))/9
# ----- Conv = f(x)*g(x) ----- #
# f(x)*g(x)
Conv = signal.convolve2d(gray, kenel, mode='full') # 使用full卷积类型，得到(M+N-1)X(M+N-1)大小


#--------- ifft{ F(f(x))·F(g(x)) } ---------#
# 傅里叶变换前进行 padding 填充，图像和卷积核都补零到 (M+N-1)x(M+N-1)大小。
img_pad = np.pad(gray, ((2, 2), (2, 2)), 'constant')
kenel_pad = np.pad(kenel, ((4, 4), (4, 4)), 'constant')

# F(f(x))
img_fft = np.fft.fftshift(np.fft.fft2(img_pad))
# F(g(x))
kenel_fft = np.fft.fftshift(np.fft.fft2(kenel_pad))
# ifft( F(f(x))·F(g(x)) )
FFT = np.fft.ifftshift(np.fft.ifft2(np.fft.fftshift(img_fft*kenel_fft)))

#--------- 打印结果 ---------#
print(" f(x) ↓")
print(gray)
print(" g(x) ↓")
print(kenel)

print("\n\n f(x)*g(x) ↓")
print(np.uint8(Conv))
print("\n ifft[F·G] ↓")
print(np.uint8(np.abs(FFT)))
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

结果展示
上述代码的结果如下图所示，可见卷积定理得到的结果与full卷积得到的结果一致。

3、“Valid”卷积

假设原始图像的大小为mxm，卷积核的大小为nxn。此时进行valid卷积，则卷积后生成的图像大小为(m-n+1)x(m-n+1)，此时不需要对被卷积图像进行padding操作，直接进行卷积操作，进行傅里叶变化是需要先将卷积核的大小padding为mxm，因为要进行傅里叶变换后的乘积操作，因此被卷积图像以及卷积核的大小需要相等。经过卷积定理的操作之后，将生成的mxm的图像进行裁剪，得到(m-n+1)x(m-n+1)的图像即为valid卷积得到的卷积图像。
Python代码验证

import numpy as np
from scipy import signal

# 原始图像 f(x)
gray = np.uint16(np.random.randint(100, size=(9, 9)))
# 卷积核 g(x)
kenel = np.ones((5, 5))/9
# ----- Conv = f(x)*g(x) ----- #
# f(x)*g(x)
Conv = signal.convolve2d(gray, kenel, mode='valid') # 使用full卷积类型，得到(M+N-1)X(M+N-1)大小

#--------- ifft{ F(f(x))·F(g(x)) } ---------#
# 傅里叶变换前进行 padding 填充，图像和卷积核都补零到 (M+N-1)x(M+N-1)大小。
img_pad = np.pad(gray, ((0, 0), (0, 0)), 'constant')
kenel_pad = np.pad(kenel, ((2, 2), (2, 2)), 'constant')

# F(f(x))
img_fft = np.fft.fftshift(np.fft.fft2(img_pad))
# F(g(x))
kenel_fft = np.fft.fftshift(np.fft.fft2(kenel_pad))
# ifft( F(f(x))·F(g(x)) )
FFT = np.fft.ifftshift(np.fft.ifft2(np.fft.fftshift(img_fft*kenel_fft)))

#--------- 打印结果 ---------#
print(" f(x) ↓")
print(gray)
print(" g(x) ↓")
print(kenel)

print("\n\n f(x)*g(x) ↓")
print(np.uint8(Conv))
print("\n ifft[F·G] ↓")
print(np.uint8(np.abs(FFT)))
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

结果展示
上述代码的结果如下图所示，可见卷积定理得到的结果经过裁剪后的红框中的内容与valid卷积得到的结果一致。

三、基于OpenCV库实现的二维图像卷积定理

//本代码实现的是使用193x193大小的卷积核对193x193大小的被卷积图像进行卷积操作
Conv2Kernel(kernelImage, inputImage, outputImage, n, n);
cv::Mat matPaded1;
cv::Mat matPaded2;
cv::Mat kernelPaded;
cv::Mat fftKernel;
cv::Mat fftMat1;
cv::Mat fftMat2;
cv::Mat result1;
cv::Mat result2;

cv::copyMakeBorder(matKernel, kernelPaded, 192, 215, 192, 215, cv::BORDER_CONSTANT, cv::Scalar(0.f));
cv::copyMakeBorder(matConv1, matPaded1, 192, 215, 192, 215, cv::BORDER_CONSTANT, cv::Scalar(0.f));

cv::Mat planes1[] = { cv::Mat_<float>(kernelPaded),cv::Mat::zeros(kernelPaded.size(),CV_32F)
};
cv::Mat planes2[] = { cv::Mat_<float>(matPaded1),cv::Mat::zeros(matPaded1.size(),CV_32F)
};
cv::merge(planes1, 2, fftKernel);
cv::merge(planes2, 2, fftMat1);

cv::dft(fftKernel, fftKernel, cv::DFT_COMPLEX_OUTPUT);
cv::dft(fftMat1, fftMat1, cv::DFT_COMPLEX_OUTPUT);

fftshift(fftKernel);
fftshift(fftMat1);
cv::Mat fftMultiplication1;
cv::mulSpectrums(fftKernel, fftMat1, fftMultiplication1, cv::DFT_ROWS);
cv::idft(fftMultiplication1, result1, cv::DFT_INVERSE + cv::DFT_SCALE + cv::DFT_COMPLEX_OUTPUT);

ifftshift(result1);
cv::split(result1, planes1);
cv::magnitude(planes1[0], planes1[1], planes1[0]);

cv::Mat matConv3 = planes1[0](cv::Rect(192, 192, 193, 193)).clone();
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

代码详解
1、对被卷积图像和卷积核执行padding操作
由于使用的是same卷积，按照上文same卷积的过程对其进行padding操作，图像上下左右四个方向都补上192个0。下述代码右边和下面之所以补零数不是192，是因为在进行傅里叶变换时，特定长度的矩阵速度更快，因此补了更多的0。但是我在实验过程中并没有发现速度有变快
2、对padding后的图像矩阵进行傅里叶变换
对padding后的卷积核以及被卷积图像执行傅里叶变换，然后进行fftshift操作。
3、对傅里叶变换后的图像矩阵进行点乘操作
4、对点乘后的举证进行傅里叶逆变换
对点乘后的矩阵执行傅里叶逆变换，然后进行ifftshift操作。
5、截取合适位置的图像，具体位置根据padding过程确定，得到的图像矩阵即为卷积后的图像矩阵。

四、基于FFTW库实现的二维图像卷积定理

                fftw_complex* pImgOut = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * 600 * 600);
                fftw_plan pRef;
                pRef = fftw_plan_dft_2d(600, 600, pImgIn, pImgOut, FFTW_FORWARD, FFTW_ESTIMATE);

                fftw_execute(pRef);
                fftw_destroy_plan(pRef);
                fftw_free(pImgIn);
                fftw_free(pImgOut);
1
2
3
4
5
6
7
8

五、总结与讨论

经过不懈努力终于使用C++实现了卷积定理，想要进一步提升卷积过程的速度，优化代码性能。但是天不遂人愿，在我的任务中使用卷积定义实现卷积过程和使用OpenCV库中的卷积操作接口速度相差无几，几乎没有看到在计算速度上有提升，忙了半天白忙活了，难受~。分析原因可能是由于本次任务的卷积过程比较特殊，卷积核、被卷积图像的大小一致，导致进行padding操作时，padding后的图像几乎是原图像的9倍大小，造成傅里叶变换以及傅里叶逆变换的过程速度变慢，影响了整体流程的速度。