Prev: Reduce elements in a point cloud
Next: Warnings engOpen - Handle Array ..., MCOS initial...
From: Pete on 21 Dec 2009 16:42
Hello all
I calculate a number of random walks, generating x and y coordinates. I wish to map these to a spatial grid (X,Y) of a different resolution. What I want to do is count the number of 'stops' (x,y pairs) that 'land' in each cell of the spatial grid (X,Y). The code below might make it a bit easier to see what I mean. There is a severe bottleneck in the last part of the code where i calculate euclidean distances. If I increase the size and number of the walks (which I really need to do, by a great deal), then this part of the code slows everything down to a crawl. Can anyone think of a way to speed that part up, or replace it?

--------------------------------------------------------------------------

n = 100; % number of steps in each walk
N = 100; % number of separate walks
minstep = 10; % minimum length of each step
mu = 2; % exponent of a power law distribution
x = zeros(N,n); % x coordinates of the walks
y = zeros(N,n); % y coordinates of the walks

angle = 2*pi*rand(N,n); % random trajectory angles
length = minstep*rand(N,n).^(1/(1-mu)); % random trajectory lengths
x = length.*cos(angle); % steps
y = length.*sin(angle); % steps
X = cumsum(x,2); % x coordinates of walk
Y = cumsum(y,2); % y coordinates of walk

% Plot walk coordinates
shg
figure(1),h = plot(X,Y,'b.');
axis square
box on
hold off;
zzz=1;

% Cartesian grid
ls = -5000; % [m]. extent to left of origin
rs = 5000; % [m]. extent to right of origin
ts = 5000; % [m]. extent above origin
bs = -5000; % [m]. extent below origin
xnp = 100; % number of x coordinates
ynp = 100; % number of y coordinates
xk = linspace(ls,rs,xnp); % vector of x coordinates
yk = linspace(ts,bs,ynp); % vector of y coordinates
[Xk,Yk] = meshgrid(xk,yk); % X and Y coordinates

% Storage array for counts
grid = zeros(size(Xk)); % [#]. Will contain a count of the
% number of walk coordinates that 'land'
% in each cell defined by the cartesian grid

% Remove random points that lie outside cartesian grid boundaries
[i] = find(X>rs);X(i) = [];Y(i) = [];
[j] = find(X<ls);X(j) = [];Y(j) = [];
[k] = find(Y>ts);X(k) = [];Y(k) = [];
[l] = find(Y<bs);X(l) = [];Y(l) = [];
nrp = numel(X);

% Assign each random point to its nearest cartesian coordinate
for i = 1:nrp
euc = sqrt((X(i)-Xk).^2+(Y(i)-Yk).^2);
grid(euc==min(min(euc))) = grid(euc==min(min(euc)))+1;
end
---------------------------------------------------------------
From: Pete on 21 Dec 2009 17:35
I have tried the following also, but the arrays get too large when I increase n and N to the values I want (e.g., n = 1000, N = 1000):

% Assign each random point to its nearest cartesian coordinate
Xrep = repmat(reshape(X,[1 1 nrp]),[ynp xnp 1]);
Yrep = repmat(reshape(Y,[1 1 nrp]),[ynp xnp 1]);
Xkrep = repmat(Xk,[1 1 nrp]);
Ykrep = repmat(Yk,[1 1 nrp]);
euc = sqrt((Xrep-Xkrep).^2+(Yrep-Ykrep).^2);

for i = 1:nrp
euk = euc(:,:,i);
grid(euk==min(min(euk))) = ...
grid(euk==min(min(euk)))+1;
end
From: ImageAnalyst on 21 Dec 2009 23:59
I really don't understand that last for loop. You have X and Y which
is a list of the X and Y coordinates. These are floating point
numbers so you just have to figure out what index they correspond to
and increment it. Something like this
for k =1:nrp
indexX = X(k);
indexY = Y(k);
grid(indexY, indexX) = grid(indexY, indexX) + 1
end

You'd need to modify that if your x and y were in the range (-5000 to
+5000) instead of (1 to nrp) but that's a simple mapping, for example
just add 5000 and divide by 100 or something like that to get the
index in the range of 1-100. This would be many, many fewer
computations than what you are doing.

From: Pete on 22 Dec 2009 06:07
ImageAnalyst <imageanalyst(a)mailinator.com> wrote in message <13e0efdd-f5d1-4dbb-9a14-0bffd471b4a3(a)m25g2000yqc.googlegroups.com>...
> I really don't understand that last for loop. You have X and Y which
> is a list of the X and Y coordinates. These are floating point
> numbers so you just have to figure out what index they correspond to
> and increment it. Something like this
> for k =1:nrp
> indexX = X(k);
> indexY = Y(k);
> grid(indexY, indexX) = grid(indexY, indexX) + 1
> end
>
> You'd need to modify that if your x and y were in the range (-5000 to
> +5000) instead of (1 to nrp) but that's a simple mapping, for example
> just add 5000 and divide by 100 or something like that to get the
> index in the range of 1-100. This would be many, many fewer
> computations than what you are doing.

Hey ImageAnalyst
The reality is that X and Y are not created using meshgrid - I just put that in to try and make my example a bit clearer. X and Y are non-uniformly spaced coordinates - covering a wide range of resolutions and a very large extent.........so it does not seem possible to create a finer grid that encompasses all X,Y pairs (to assist in indexing). I was hoping someone would know a nice trick or command that could help...........
From: ImageAnalyst on 22 Dec 2009 10:41
Pete:
Yes I realized that. And I thought my last paragraph was the nice
"trick" that would speed it up for you.

Why calculate the min of a huge bunch of numbers, and then compare it
to that same bunch of numbers just to find an index that you already
have? And then you do this unnecessary stuff not once, but twice.
And it's not only twice, it's twice for every one of your ten thousand
numbers - what a huge waste of time and effort for something that can
be done so simply and efficiently.
-ImageAnalyst
 |  Next  |  Last
Pages: 1 2
Prev: Reduce elements in a point cloud
Next: Warnings engOpen - Handle Array ..., MCOS initial...