子矩阵（和子向量）上的操作可以向量化吗？

Question

我目前正在研究 octave 的边缘检测器。来自 Java 和 Python 等其他编程语言，我习惯于在 for 循环中进行迭代，而不是对整个矩阵执行操作。现在在八度音程中，这会导致严重的性能下降，而且我很难弄清楚如何对我的代码进行矢量化。我有以下两段代码：

1)

function zc = ZeroCrossings(img, T=0.9257)
  zc = zeros(size(img));

  # Iterate over central positions of all 3x3 submatrices
  for y = 2:rows(img) - 1
    for x = 2:columns(img) - 1
      ndiff = 0;

      # Check all necessary pairs of elements of the submatrix (W/E, N/S, NW/SE, NE/SW)
      for d = [1, 0; 0, 1; 1, 1; 1, -1]'
        p1 = img(y-d(2), x-d(1));
        p2 = img(y+d(2), x+d(1));
        if sign(p1) != sign(p2) && abs(p1 - p2) >= T
          ndiff++;
        end
      end

      # If at least two pairs fit the requirements, these coordinates are a zero crossing
      if ndiff >= 2
        zc(y, x) = 1;
      end
    end
  end
end

2)

function g = LinkGaps(img, k=5)
  g = zeros(size(img));

  for i = 1:rows(img)
    g(i, :) = link(img(i, :), k);
  end
end

function row = link(row, k)
  # Find first 1
  i = 1;
  while i <= length(row) && row(i) == 0
    i++;
  end

  # Iterate over gaps
  while true
    # Determine gap start
    while i <= length(row) && row(i) == 1
      i++;
    end
    start = i;

    # Determine gap stop
    while i <= length(row) && row(i) == 0
      i++;
    end

    # If stop wasn't reached, exit loop
    if i > length(row)
      break
    end

    # If gap is short enough, fill it with 1s
    if i - start <= k
      row(start:i-1) = 1;
    end
  end
end

这两个函数都迭代子矩阵（或第二种情况下的行和子行），尤其是第一个似乎大大减慢了我的程序。

1. 此函数采用像素矩阵 (img) 并返回二进制 (0/1) 矩阵，其中找到零交叉点（对应的 3x3 邻域满足特定要求的像素）为 1。

   外部的 2 个 for 循环似乎应该可以以某种方式矢量化。我可以将主体放入它自己的函数中（将必要的子矩阵作为参数），但我不知道如何在所有子矩阵上调用此函数，将它们对应的（中心）位置设置为返回值。

   如果内部 for 循环也可以向量化，则加分。

2. 此函数从前一个输出中获取二进制矩阵，并填充其行中的空白（即将它们设置为 1）。间隙定义为一系列长度

   现在我确定至少外循环（LinkGaps 中的那个）是可矢量化的。但是，link 中的 while 循环再次对子向量进行操作，而不是单个元素，所以我不确定如何对其进行向量化。

Answer 1

不是一个完整的解决方案，但这里有一个想法，你可以如何在没有任何循环的情况下完成第一个：

% W/E
I1 = I(2:end-1,1:end-2);
I2 = I(2:end-1,3:end  );
C = (I1 .* I2 < 0) .* (abs(I1 - I2)>=T);

% N/S
I1 = I(1:end-2,2:end-1);
I2 = I(3:end,  2:end-1);
C = C + (I1 .* I2 < 0) .* (abs(I1 - I2)>=T);

% proceed similarly with NW/SE and NE/SW
% ...

% zero-crossings where count is at least 2
ZC = C>=2;

想法：形成两个适当移动的子图像，检查符号差异（乘积负数）并阈值差异。两个测试都返回一个逻辑 (0/1) 矩阵，逐元素乘积执行逻辑，结果是一个 0/​​1 矩阵，其中两个测试都成功了。可以添加这些矩阵来跟踪计数 (ndiff)。