beliavsky@aol.com wrote:

> David.Paterson@csiro.au wrote:
>
>>I'm looking for Fortran source code if possible.
>>
>>Given set A with roughly 10,000 (x,y,z) triples and set B with roughly
>>1,000,000 (x,y,z) triples.
>>
>>For each point in set A find the nearest point from set B.
>>
>>What's the quickest algorithm?
>>
>>-------------------------------
>>
>>A very slow algorithm would be to check all 1,000,000 triples in B for
>>every point in A.
>
>
> Suppose v1 and v2 are both 3-d vectors of REALs. I wonder if Fortran
> compilers, given the expression
>
> if ((v1-v2)**2 > c)
>
> are "smart" enough to recognize that if (v1(1)-v2(1))**2 > c,
> computations involving v1(2:3) and v2(2:3) can be skipped.
>
> I also wonder if coding this logic explicitly would result in a slower
> program than just doing the comparison as above. I may check later when
> I have time.
>
> Btw a similar question was asked recently in comp.lang.fortran -- one
> can Google "kd tree" there.
>
> (comp.lang.fortran added to follow-ups)
>

I've attacked this problem in the past in the following way (2D problem
in my case).
First divide the region into cubes, given, say, by C(i,j,k), i=1,..,Nx,
j=1,..,Ny, k=1,..,Nz, where C(i,j,k) spans the region
Xmin + (i-1)*delta < x <= Xmin + i*delta
Ymin + (j-1)*delta < y <= Ymin + j*delta
Zmin + (k-1)*delta < z <= Zmin + k*delta
and delta, Nx, Ny, Nz are suitably chosen.

Then for each point in B determine the cube it falls in, i.e. the values
of (i,j,k). This is very fast. Create a list of B points for each
cube, i.e. for each triple (i,j,k).

Then when looking for the nearest point in set B to a given point in set
A that is in the cube (ia,ja,ka), you just need to examine the B points
in the cube (ia,ja,ka) and the adjacent cubes (27 altogether, in
general). If delta is badly chosen, or if your points are very
irregularly distributed, you might not find any B points in these cubes,
and then you'll need to look at the next layer of surrounding cubes,
etc. - the method starts to lose its shine a bit if this happens. I
guess delta should be chosen so that every cube contains at least one B
point.