I have a program to draw a grid of ellipses with a uniform phase distribution. However, it is very slow.

I'd like my code to be faster, so that I can use, for example, N = 150, M = 150. How can I speed up this code?

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

This is what the code produces:

Picture

Updated:

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

2 Answers

1
Cris Luengo On Best Solutions

There are many things you can do to speed up the computation. One important one is to remove loops and replace them with vectorized code. Octave works much faster when doing many computations as once, as opposed to one at a time.

For example, instead of

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

one can write

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

However, here we still do too much work because we compute index at locations that we know will be outside the extend of the ellipse. Ideally we'd find a smaller region around each point (x,y) that we know the ellipse fits into.

Another important thing to do is preallocate the array. f grows within the loop iterations. Instead, one should create f with its final size before the loop starts.

There are also many redundant computations being made. For example w.*pi is computed multiple times, and the cos and sin of it too. You also assign cos(0.001)+i.*sin(0.001) to output pixels over and over again, this could be a constant computed once.

The following code runs in MATLAB in a tiny fraction of a second (though Octave will be a lot slower). I've also separated N and M properly (so the output is not always square) and made the step sizes a variable for improved understanding of the code. I'm assigning 1 to the ellipse elements, you can replace them by your constant by multiplying 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

output of algorithm

0
user1543042 On

I'm not 100% sure what you're trying to do. See the code below and let me know. It took me about 30 s to run the 150 x 150 case. Not sure if that's fast enough

[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{:})