Matlab/Octave 使用哪种方法来实现 conv2 函数？

Question

Matlab/Octave 使用哪种方法conv2函数？

是吗：

• 快速傅立叶变换
• 3 个或更多 for 循环，例如经典迭代
• 其他？

我正在寻找最快的方法conv2。我将为它编写 C 代码。

Answer 1

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

FFT 实现的时间复杂度为 O( n log n )，与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 计算卷积需要 3 倍的fft2调用，加上复数值乘法，因此在这种情况下比直接实现慢 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，除非您的内核非常大。如果您确实想自己实现这一点，请参阅我的另一个答案，了解有关如何有效实现卷积的一些技巧。

如果使用FFT来实现卷积，则需要接受周期性边界条件，或者充分填充图像以避免边界效应。尽管如果您不想在图像之外假设零，那么您也需要在直接方法中进行这种填充conv2（这根本没有用）。

备择方案

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

• 如果内核是可分离的，您可以沿着图像行应用一维卷积，然后沿着图像列应用一维卷积，以获得最终结果，从而大大减少计算工作量（如果内核是n-by- m，您可以这样做n+m而不是n*m乘法和加法。

• 如果内核具有统一的权重（如我上面使用的示例），那么您可以以独立于内核中像素数量的成本来应用它。当您移动内核时，您会减去左侧退出内核的像素，并添加进入右侧内核的像素。

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

高斯滤波器

DIPlib 具有高斯滤波器的三种实现：FIR（直接方法，作为沿每个轴的一系列一维滤波器）、IIR 和 FT（使用 FFT）。以下是针对不同 sigma 的 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)')

• 冷杉：
  • 一维内核大小为2*ceil(3*sigma)+1。
  • 该实现利用了内核是对称的这一事实，因此k它不进行乘法和加法，而是进行k加法和k/2乘法。但除此之外，代码相当简单（例如，它是符合标准的 C++，没有显式使用 SIMD 指令，尽管编译器应该在适当的时候使用这些指令）。
  • 正如您所期望的，这在时间上不是线性的，因为当沿列应用过滤器时，缓存的使用效果较差。因此，对于较大的西格玛，内存访问的成本相对较高。
• IIR：
  • 这个实现相当古老，从 1998 年开始，纯 C 实现。
  • IIR 滤波器是沿每个轴依次应用的一维滤波器，与 FIR 滤波器类似。
• 英国金融时报：
  • 对于 FFT，它使用 PocketFFT（我可以使用 FFTW 构建 DIPlib，但在我的机器（M1 Mac）上，这没有太大区别，因为 FFTW 没有特定于 M1 芯片的代码）。
  • 它仅使用两种 FFT：一种用于图像前向，一种用于结果后向。内核本身没有经过变换，而是直接在频域中生成。高斯是可分离的，您可以生成一维内核并将其乘以 FFT 的列和行。
  • 对于较小的 sigma，此方法速度较慢，因为它在频域中生成内核，其中较大的 sigma 生成较小的内核（需要乘以更多的零）。

索贝尔过滤器

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