Vectorizing code to calculate radial distance for all points

5 vues (au cours des 30 derniers jours)
Cong Bang Huynh
Cong Bang Huynh le 12 Jan 2013
Hi all!
I have currently written this code to calculate the radial distance for all points with their x, y and z coordinates.
for x = xmin:step:xmax
for y = ymin:step:ymax
for z = zmin:step:zmax
r((x-xmin)/step+1,(y-ymin)/step+1,(z-zmin)/step+1) = sqrt(x^2 + y^2 + z^2);
end
end
end
I am wondering if there is a way to vectorize this code so that I don't have to run the nested loops and hence speed up my calculation time, because I have to run the code for quite a number of times to analyze different data sets.
I am new to matlab so this question may seem quite trivial but please do help me :) Thank you in advance.
  1 commentaire
Roger Stafford
Roger Stafford le 12 Jan 2013
Aside from the question of vectorizing code, it is unwise to set up for-loops in the way you show here. The problem is that expressions such as (x-xmin)/step+1 used as indices may not turn out to be exact positive integers due to tiny round-off errors and thereby would produce error messages. For example, with xmin = 0, step = 0.1, and xmax = 1, my machine produces a non-integer for (x-xmin)/step+1 at the seventh step and would complain if I used it as an index. With for-loops you should instead do something like the following so that you are guaranteed to have integer values for your indices:
x = xmin:step:xmax;
y = ymin:step:ymax;
z = zmin:step:zmax;
for kx = 1:length(x)
for ky = 1:length(y)
for kz = 1:length(z)
r(kx,ky,kz) = sqrt(x(kx)^2 + y(ky)^2 + z(kz)^2);
end
end
end

Connectez-vous pour commenter.

Réponse acceptée

Matt J
Matt J le 12 Jan 2013
x = (xmin:step:xmax).';
y = ymin:step:ymax;
z = reshape(zmin:step:zmax,1,1,[]);
r=sqrt(bsxfun(@plus, bsxfun(@plus,x.^2,y.^2),z.^2 ));

Plus de réponses (1)

Roger Stafford
Roger Stafford le 12 Jan 2013
You can also do:
[X,Y,Z] = ndgrid(xmin:step:xmax,ymin:step:ymax,zmin:step:zmax);
r = sqrt(X.^2+Y.^2+Z.^2);

Catégories

En savoir plus sur Spline Postprocessing dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by