# 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!

• How many points do you have in how many dimensions? Knowing this would help with optimizing solutions. Jan 10, 2012 at 10:57
• @Szabolcs See my edit, but if you also have good idea for lower dimensions(1D~3D), please put it forward. Jan 10, 2012 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. Sep 15, 2012 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. Sep 15, 2012 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 = {#[], 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? Jan 10, 2012 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. Jan 10, 2012 at 12:58
• @kptnw My edit might address your situation more precisely. Jan 10, 2012 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}] Jan 12, 2012 at 12:03
• Actually, on version 7 this appears to be better: Join @@ f /@ b --- does your code do something special on version 8? Jan 12, 2012 at 12:06
• It does not unpack, the speed difference is minimal on my computer, but my approach uses less memory.
– user21
Jan 12, 2012 at 12:25
• I forgot about that. It looks like Join @@` was useful at least. Jan 12, 2012 at 12:47