Vectorization of this loop

1 次查看（过去 30 天）
Dilunath Hareendranath
The following loop calculates the distance and angle values of every location from a point and stores in arrays named Radius and theta. This loop is called nearly 3600 times in the code. This loop is effecting the performance of the code. Please suggest some ways to vectorise this loop.
xwidth and ywidth varies from 500 to 750. So, memory needed is also very high. Please suggest ways to decrease the memory needed.
x1=0;
y1=1;
inj_x=round(xwidth/2.0);
inj_y=round(ywidth/2.0);
Radius=zeros(ywidth,xwidth);
theta=zeros(ywidth,xwidth);
for r=1:ywidth
for c=1:xwidth
x2=r-inj_y;
y2=c-inj_x;
Radius(r,c)=(x2^2+y2^2)^.5;
theta(r,c)=mod(atan2(x1*y2-x2*y1,x1*x2+y1*y2),2*pi);
end
end
Thanks in advance

请先登录，再回答此问题。

采纳的回答

Andrei Bobrov
Andrei Bobrov 2012-4-18
in your case
r = (1:ywidth).' - round(ywidth/2);
c = (1:xwidth) - round(xwidth/2);
Radius = bsxfun(@hypot,r,c);
theta = mod(bsxfun(@atan2,-r,c),2*pi);
  1 个评论
Dilunath Hareendranath
Dilunath Hareendranath 2012-4-18
Thanks andrei.. This code is taking only 0.09 seconds to process.. whereas previous code was taking 2.8 seconds.

请先登录，再进行评论。

更多回答（2 个）

Honglei Chen
Honglei Chen 2012-4-18
x1=0;
y1=1;
inj_x=round(xwidth/2.0);
inj_y=round(ywidth/2.0);
[x2,y2] = ndgrid((1:ywidth)'-inj_y,(1:xwidth)'-inj_x);
Radius=(x2.^2+y2.^2).^.5;
theta=mod(atan2(x1.*y2-x2.*y1,x1.*x2+y1.*y2),2*pi);
  1 个评论
Dilunath Hareendranath
Dilunath Hareendranath 2012-4-18
Thanks Honglei Chen for answering. This code is also taking less time.

请先登录，再进行评论。

Jan
Jan 2012-4-18
For a fair speed comparison cleanup the loops:
  • move all repeated calculation outside
  • SSQRT() is faster than ^0.5
twoPi = 2 * pi;
for r = 1:ywidth
x2 = r - inj_y;
x2_2 = x2 * x2;
x1x2 = x1 * x2;
y1x2 = y1 * x2;
for c = 1:xwidth
y2 = c-inj_x;
Radius(r,c) = sqrt(x2_2 + y2^2);
theta(r,c) = mod(atan2(x1*y2 - y1x2, x1x2 + y1*y2), twoPi);
end
end
Perhaps a partial vectorization is fastest:
twoPi = 2*pi;
for c = 1:xwidth
y2 = c-inj_x;
x2 = transpose(1-inj_y:ywidth - inj_y);
Radius(:,c) = sqrt(x2.^2 + y2^2);
theta(:,c) = mod(atan2(x1*y2-x2*y1, x1*x2+y1*y2), twoPi);
end
And fully vectorized:
x2 = transpose(1 - inj_y:ywidth - inj_y);
y2 = 1 - inj_x:xwidth - inj_x;
Radius = sqrt(bsxfun(@plus, x2 .^ 2 + y2 .^ 2);
k1 = bsxfun(@minus, x1 * y2, y1 * x2);
k2 = bsxfun(@plus, x1 * x2, y1 * y2);
theta = mod(bsxfun(@atan, k1, k2), 2*pi);
And if x1 and y1 are really fixed to 0 and 1 this can be simplified.

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by