C++ MatlabOctave对conv2函数使用哪种方法

conv2实现直接方法(即4个循环).这不是O(n⁴)，而是O(nk)，其中n是图像中的像素数，k是内核中的非零像素数.

FFT实现是O(nlogn)，独立于k，但常量非常大，因此您需要一个较大的k才能使其比直接实现更高效.

您可以简单地通过比较计算fft2(img)和conv2(img,ones(3,3))所需的时间来测试这一点:

img = ones(1024,1024);
k = ones(3,3);
timeit(@()fft2(img))
timeit(@()conv2(img,k,'same'))

ans =
    0.0059

ans =
    0.0015

通过FFT计算卷积需要调用fft2次的3倍，外加复值乘法，因此比本例中的直接实现慢12倍以上.

您还可以看到conv2次随着内核大小的增加而增加:

>> k = ones(7,7); timeit(@()conv2(img,k,'same'))
ans =
    0.0066

>> k = ones(19,19); timeit(@()conv2(img,k,'same'))
ans =
    0.0152

如果它使用FFT卷积，情况就不会是这样.

请注意，MATLAB的卷积实现非常高效，很难将其与您自己的C++代码相匹配.尽管如此，即使使用简单的实现，您也不希望使用FFT，除非您的内核非常大.如果您确实想自己实现这一点，请参阅this other answer of mine获取有关如何有效地实现卷积的一些提示.

如果使用FFT实现卷积，则需要接受周期性边界条件，或者充分填充图像以避免边界效应.不过，如果您不想像conv2那样在图像外假定为零(这根本没有用)，那么在Direct方法中也需要这种填充.

替代方案

对于某些内核，有比FFT或conv2方法更有效的实现:

• 如果内核是可分离的，您可以沿图像行应用1D卷积，然后沿图像列应用1D卷积，以在极大地减少计算工作量的情况下获得最终结果(如果内核是n乘m，则您将进行n+m次而不是n*m次乘加.

• 如果内核具有统一的权重(就像我在上面使用的示例中一样)，那么您可以使用与内核中的像素数无关的成本来应用它.在移动内核时，会减go 左侧退出内核的像素，然后添加右侧进入内核的像素.

• 如果内核是高斯或Gabor内核，您还可以使用众所周知的IIR实现(其成本与内核的大小无关).IIR=无限脉冲响应.

高斯型滤波器

DIPlib有three implementations of the Gaussian filter:FIR(直接方法，每个轴上都有一系列一维滤波器)、IIR和FT(使用FFT).以下是它们对2048x2048图像、不同西格玛(使用直接调用DIPlib功能的DIPImage工具箱)的计时的比较:

img = randn(2048, 2048, 'single');
sigma = [1, 2, 3, 5, 7, 10, 15, 20];
t = zeros(numel(sigma), 3);
for ii = 1:numel(sigma)
    t(ii, 1) = timeit(@()gaussf(img, sigma(ii), 'fir'));
    t(ii, 2) = timeit(@()gaussf(img, sigma(ii), 'iir'));
    t(ii, 3) = timeit(@()gaussf(img, sigma(ii), 'ft'));
end

clf
set(gca, 'FontSize', 12)
plot(sigma, t * 1000, '.-', 'LineWidth', 1)
legend({'FIR', 'IIR', 'FT'})
xticks(sigma)
xlabel('Sigma of Gaussian (kernel is 2\lceil 3\sigma\rceil + 1)')
ylabel('Time (ms)')

• FIR:
  • 一维内核大小为2*ceil(3*sigma)+1.
  • 这个实现利用了内核是对称的这一事实，所以它不做k次乘法和加法，而是做k次加法和k/2次乘法.但除此之外，代码非常简单(例如，它是符合标准的C++，没有显式使用SIMD指令，尽管编译器应该在适当的时候使用这些指令).
  • 正如您所预期的那样，这在时间上不是线性的，因为在沿列应用过滤器时，缓存的使用不是最优的.因此，对于较大的西格玛，存储器访问的成本相对较大.
• IIR:
  • 这个实现很老了，从1998年开始，纯C语言.
  • IIR过滤器是沿每个轴依次应用的一维过滤器，与FIR过滤器类似.
• 《金融时报》:
  • 对于FFT，这使用了PocketFFT(我可以使用FFTW构建DIPlib，但在我的机器M1Mac上，这没有太大区别，因为FFTW没有特定于M1芯片的代码).
  • 它只使用了两个FFT:一个用于图像，一个用于结果.内核本身不进行变换，而是直接在频域中生成.由于高斯是可分离的，您可以生成1D内核并将其乘以FFT的列和行.
  • 对于较小的sigma，此方法速度较慢，因为它在频域中生成核，而较大的sigma会生成较小的核(与之相乘的零数更多).

Sobel滤波器

Sobel过滤器是3x3，但只有6个非零权重.6乘法和加法的应用非常便宜，计算一次FFT从来都不便宜，更不用说三次了.请注意，遍历图像是这里最昂贵的部分(它需要从内存获取数据，这需要比6次乘法和加法更多的时间).如果您想要计算这两个导数，那么将它们组合到图像上的相同循环中是值得的，以最好地利用缓存.