drawing - 为什么这个椭圆绘图程序这么慢？

Question

我有一个程序来绘制具有均匀相位分布的椭圆网格。但是，它非常缓慢。

我希望我的代码更快，以便我可以使用，例如N = 150，M = 150. 如何加快此代码的速度？

N = 10;
M = 10;
y = 1;
x = 1;
a = 1;
b = 2;
for k = 1:N
  for m = 1:N
    w = rand(1,1);
    for l = 1:N
      for s = 1:N
        if(((l-x)*cos(w*pi)+(s-y)*sin(w*pi)).^2/a^2 + (-(l-x)*sin(w*pi) + (s-y)*cos(w*pi)).^2/b.^2 <= 1)
          f(l,s) = 1*(cos(0.001)+i*sin(0.001));
        end
      end
    end
    y = y+4;
  end
  y = 1;
  x = x+5;
end
image(arg(f),'CDataMapping','scaled');

这是代码产生的：

更新：

N = 10;
M = 10;
y = 1;
x = 1;
a = 1;
b = 2;
for x = 1:5:N
  for y = 1:4:N
    w = rand(1);
    for l = 1:N
      for s = 1:N
        if(((l-x).*cos(w.*pi)+(s-y).*sin(w.*pi)).^2/a.^2 + (-(l-x).*sin(w.*pi) + (s-y).*cos(w.*pi)).^2/b.^2 <= 1)
          f(l,s) = cos(0.001)+i.*sin(0.001);
        end
      end
    end
  end
  y = 1;
end
image(arg(f),'CDataMapping','scaled');

标签： drawingoctaveellipse

Answer 1

您可以做很多事情来加快计算速度。一个重要的方法是删除循环并用矢量化代码替换它们。当一次执行多个计算时，Octave 的工作速度要快得多，而不是一次一个。

例如，而不是

for l = 1:N
  for s = 1:N
    if(((l-x).*cos(w.*pi)+(s-y).*sin(w.*pi)).^2/a.^2 + (-(l-x).*sin(w.*pi) + (s-y).*cos(w.*pi)).^2/b.^2 <= 1)
      f(l,s) = cos(0.001)+i.*sin(0.001);
    end
  end
end

一个可以写

l = 1:N;
s = (1:N).';
index = ((l-x).*cos(w.*pi)+(s-y).*sin(w.*pi)).^2/a.^2 + (-(l-x).*sin(w.*pi) + (s-y).*cos(w.*pi)).^2/b.^2 <= 1;
f(index) = cos(0.001)+i.*sin(0.001);

然而，在这里我们仍然做太多的工作，因为我们index在我们知道将在椭圆延伸之外的位置进行计算。(x,y)理想情况下，我们会在我们知道椭圆适合的每个点周围找到一个较小的区域。

另一件重要的事情是预分配数组。f在循环迭代中增长。相反，应该f在循环开始之前创建其最终大小。

还进行了许多冗余计算。例如w.*pi被多次计算，cosand sinof it 也是。您还cos(0.001)+i.*sin(0.001)可以一遍又一遍地分配给输出像素，这可能是一次计算的常数。

以下代码在 MATLAB 中运行只需几分之一秒（尽管 Octave 会慢很多）。我也正确地分开N了M（所以输出并不总是方形的）并将步长设置为变量以提高对代码的理解。我分配1给椭圆元素，您可以通过乘以用常量替换它们f = f * (cos(0.001)+i*sin(0.001))。

N = 150;
M = 200;
a = 5;
b = 10;
x_step = 25;
y_step = 25;
f = zeros(N,M);
for x = x_step/2:x_step:M
  for y = y_step/2:y_step:N
    phi = rand(1)*pi;
    cosphi = cos(phi);
    sinphi = sin(phi);
    l = (1:M)-x;
    s = (1:N).'-y;
    index = (l*cosphi+s*sinphi).^2/a.^2 + (-l*sinphi + s*cosphi).^2/b.^2 <= 1;
    f(index) = 1;
  end
end