# Efficient way of finding shortest distance between two sets of points in mathematica

I have two sets of 3D points, say

a = RandomReal[1, {10, 3}];
b = RandomReal[1, {10, 3}];


I wanna find the first N pairs that have shortest distance between the two sets, my current approach is to form a NearestFunction from the first list and map it to the second(say, N=3):

f = Nearest[a];
c = Flatten[f/@b, 1];


Then get the result:

N=3;
Transpose[{c[[#]], b[[#]]}]&@ Ordering[Norm/@(c-b),N]

Output={{1st pair},{2nd pair},{3rd pair}}


Question:What is the more efficient way of finding shortest distance pairs of two sets in mathematica?

-edit-

In my actual case, I want to determine the intersection of two shapes(manifolds), but the analytic form of these shapes are not available. I can only compute points on these shapes through numerical process(NDSolve). So I want to use enough many points to approximate these shapes, and use Intersection-like operation to find points that lay approximately on both shapes. The number of points is very large(eg. 10^7), and the points may in high dimensional space(eg. 4D).

Excuse me for my crippled English. I am greatly appreciate any suggestions, thank you very much!

-

## migrated from stackoverflow.comSep 15 '12 at 0:19

This question came from our site for professional and enthusiast programmers.

How many points do you have in how many dimensions? Knowing this would help with optimizing solutions. –  Szabolcs Jan 10 '12 at 10:57
@Szabolcs See my edit, but if you also have good idea for lower dimensions(1D~3D), please put it forward. –  kptnw Jan 10 '12 at 11:48
@belisarius Interesting that this was migrated 8 months after it was asked and answered. Is there a reason? Perhaps a concerted effort to more stackoverflow questions over here? I mean, I got more! There's actually just a couple that might make sense. –  Mark McClure Sep 15 '12 at 2:09
@MarkMcClure I was in need of something like this and found this question and your excellent answer. After consulting it in chat with other users, I asked for migration because I think it is better to have it here as a reference (and surely it is more visible this way). Now, thinking again about it, perhaps I should have asked you if you agree. Too late for that, sorry. Re: concerted efforts. Not AFAIK, and surely not a massive movement, but from time to time I ask for migration of some good material that I think could be enhanced or expanded here. –  belisarius Sep 15 '12 at 3:50

One strength of InterpolatingFunctions is that they can be used just about like any other function. Thus, a more or less analytic approach may be possible. It's hard to say for sure, without seeing your example, but here's a contrived example.

pts1 = Table[{x, y, Sin[x] + Cos[y + 1]},
{x,-10,10}, {y,-10,10}];
f1 = Interpolation[Flatten[pts1, 1]];
pts2 = Table[{x, y, x^2 + y^2 + x + 2 y + 1},
{x,-10,10}, {y,-10,10}];
f2 = Interpolation[Flatten[pts2, 1]];
NMinimize[f2[x, y] - f1[x, y], {x, y}]


Edit

Here's an example that is, perhaps, more in line with your problem. We start with two periodically perturbed spheres - one centered at the origin and one shifted along the y-axis.

pts1 = Flatten[Table[{{p, t}, (1 + 0.1 Sin[6t]Sin[5p])*
{Sin[p] Cos[t], Sin[p] Sin[t], Cos[p]}},
{p, 0, Pi, Pi/12}, {t, 0, 2 Pi, Pi/12}], 1];
pts2 = pts1 /. {x_,y_,z_} -> {x,y+0.75,z};


We can interpolate these points and then plot the surfaces parametrically.

f1 = Interpolation[pts1];
f2 = Interpolation[pts2];
surfaces = ParametricPlot3D[{f1[p, t], f2[p, t]},
{p, 0, Pi}, {t, 0, 2 Pi}, MeshStyle -> None,
PlotStyle -> {
Directive[Yellow, Specularity[White, 40]],
Directive[Red, Specularity[White, 40]]},
ViewPoint -> {-3, -0.7, 1}]


The question is, where do these two surfaces intersect? The answer is provided by the equation f1[p1,t1]==f2[p2,t2]. Since f1 and f2 are both three dimensional vectors, this represents three equations in four unknowns. We can reduce the number of variables by fixing one and solving for the other three. For example, here is a point that lies on both perturbed spheres and at angle of Pi/4 from the positive z-axis.

Clear[t1, t2, p2];
{t1, t2, p2} = {t1, t2, p2} /. FindRoot[f1[Pi/4, t1] == f2[p2, t2],
{t1, 2.7, 2.8}, {t2, 3.5, 3.6}, {p2, Pi/4, Pi/4 + 0.1}];
{f1[Pi/4, t1], f2[p2, t2]}


Let's parametrize this process in terms of p1.

Clear[t1, t2, p2];
results = Table[{{p1, t1}, {p2, t2}} /.
FindRoot[f1[p1, t1] == f2[p2, t2],
{t1, 2.7, 2.8}, {t2, 3.5, 3.6},
{p2, p1, p1 + 0.1}], {p1, 0.5, 2.6, 0.1}]; // Quiet


We've suppressed some errors so let's check the results.

test[{pt1_, pt2_}] := Norm[f1 @@ pt1 - f2 @@ pt2];
{#, test[#]} & /@ results


Looks pretty good. Now, with these points, we can parametrize the intersection using an Interpolation, just as we parametrized the surfaces.

intersection = {#[[1]], f1 @@ #} & /@ First /@ results;
p = Interpolation[intersection]


Here's a visualization of the result.

intersectionPic = ParametricPlot3D[p[t], {t, 0.5, 2.8},
PlotStyle -> Directive[Thickness[0.02], Black]];
Show[{surfaces, intersectionPic}]


-
This idea is interesting. May be the last line in your code should be NMinimize[f2[x,y]-f1[x,y]// Abs, {x, y}]? But:First, you only get one point,the actual intersection may have one dimensional structure. Second, your pts1&pts2 lay in two dimensional surfaces, what if I want find the intersection of two 3D solids? Can InterpolatingFunction still apply? –  kptnw Jan 10 '12 at 12:53
@kptnw I don't need Abs in my example since I know that f2[x,y]>f1[x,y]. As I understand your problem, you might need to use FindRoot after parametrizing one variable. Again, it's hard to say for sure without more information on your setup. –  Mark McClure Jan 10 '12 at 12:58
@kptnw My edit might address your situation more precisely. –  Mark McClure Jan 10 '12 at 19:58

I think your approach is good, here is a small improvement:

c = Flatten[
DeveloperToPackedArray@Table[f[b[[i]], 1], {i, Length[b]}], 1];


edit: shorter code:

c = Join @@ Table[f[b[[i]], 1], {i, Length[b]}]

-
This appears to be faster, and much cleaner: c = Join @@ Table[f[i, 1], {i, b}] –  Mr.Wizard Jan 12 '12 at 12:03
Actually, on version 7 this appears to be better: Join @@ f /@ b --- does your code do something special on version 8? –  Mr.Wizard Jan 12 '12 at 12:06
It does not unpack, the speed difference is minimal on my computer, but my approach uses less memory. –  user21 Jan 12 '12 at 12:25
I forgot about that. It looks like Join @@` was useful at least. –  Mr.Wizard Jan 12 '12 at 12:47