How to find vector indexes without for loops?

I'm trying to speed up my code, which essentially computes distances of 9 million points to 72k points and looks up values based on that distance.

I have one large mesh array (3,000 x 3,000 feet in X/Y coord.) which I can step through element by element. For each element, I need to compute the distance to 72001 points (xpos_vector and ypos_vector). I then round that distance and use a lookup table to get pre-computed values.

I tried looping 72001 times within the i/j loops, but that is taking an enormous amount of time. I was hoping to "vectorize" the distances, but the issue is that the find function won't allow me to find which index each element of the vector corresponds to, only the few indexes which happen to match perfectly.

I have on 72001 element vector and I need to return a 72001 element vector of indexes based on a lookup vector (of different size.) Is there any way to do this without a for loop? I.E. - it would be perfect if index=find(vector.==lookup) would run and return an array of size(vector).

[Image Analyst]

Do you not get an Out of Memory error? That's 2.6 terabytes of memory (9e6*72e3*4) and that's if you use single precision instead of double precision.

[Ed]

If I try to pre-allocate a large matrix ( matrix=zeros(3000,3000,72001);) then yes, I do get a memory error...

[Matt J]

Ed Commented:

ok well I've come a bit further, but now I'm stuck trying to figure out how to feed an array of indexes (72001 elements) to a smaller array (4500 elements) and return a larger array (72001 elements) of values corresponding to the indexes.

Without using an additional for-loop, I need to pass in an arbitrary number of indexes to a lookup array and return an array of the arbitrary size.

E.G. (smaller scale)

xpos_vector=rand(100,1);
xpos_vector=xpos_vector.*100;
ypos_vector=rand(100,1);
ypos_vector=ypos_vector.*100;
round_distance=zeros(100,1);
for i=1:1:100
for j=1:1:100
round_distance = round(sqrt((xpos_vector-j).^2 + (ypos_vector-i).^2));
end
end
MCNP_vector = rand(50,1); %this is pre-computed lookup values
mcnpvals_vector=zeros(100,1);
mcnpvals_vector=MCNP_vector(round_distance); %obviously this does not work

[Matt J]

It's not clear what your latest example is trying to accomplish because, within your double for-loop, round_distance never accumulates anything over i and j. In other words, this

    for i=1:1:100
     for j=1:1:100
      round_distance = round(sqrt((xpos_vector-j).^2 + (ypos_vector-i).^2));
     end
    end

is equivalent to this one line:

round_distance = round(sqrt((xpos_vector-100).^2 + (ypos_vector-100).^2));

It's also a bit confusing that you're looping over 100 different i,j pairs if round_distance is only 100x1.

I'm stuck trying to figure out how to feed an array of indexes (72001 elements) to a smaller array (4500 elements) and return a larger array (72001 elements) of values corresponding to the indexes.

The size of the indices versus the lookup data doesn't seem like it should really be the confounding issue. It is generally true in MATLAB that the output of an indexing operation will always be the same size of the indexing data and independent of the size of the lookup data. A trivially example with 1x1 lookup data would be

    >> lookup=10; indices=[1 1 1  1 1 ]; lookup(indices)

    ans =

        10    10    10    10    10

Answer 1

[Matt J]

You could use

[~,index]=histc(vector,binEdges)

where the binEdges are derived from your lookup data.

[Ed]

That's a clever idea, however as the Image Analyst pointed out above index would be a 3,000 x 3,000 x 72001 matrix I think.

I'm going to try a different method of pre-computing the distance lookup such that the distance in feet is the index.