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

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');

Answer 1

我不是 100% 确定你想做什么。请参阅下面的代码并告诉我。运行 150 x 150 的案例我花了大约 30 秒。不确定是否足够快

[h,k] = meshgrid(0:150, 0:150);

a = 0.25;
b = 0.5;

phi = reshape( linspace( 0 , 2*pi , numel(h) ), size(h));

theta = linspace(0,2*pi,50)';

x = a*cos(theta);
y = b*sin(theta);

h = h(:);
k = k(:);
phi = phi(:);

ellipseX = arrayfun(@(H,K,F) x*cos(F)-y*sin(F) + H , h,k,phi, 'uni', 0);
ellipseY = arrayfun(@(H,K,F) x*sin(F)+y*cos(F) + K , h,k,phi, 'uni', 0);

ellipse = [ellipseX, ellipseY, repmat({'r'}, size(ellipseX))]';

fill(ellipse{:})

Answer 2

您可以做很多事情来加快计算速度。一个重要的方法是删除循环并用矢量化代码替换它们。 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 被多次计算，cos 和 sin 也是如此。您还分配 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

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

Why is this ellipse drawing program so very slow?

drawing

ellipse

octave