Find optimal match of points in 3-Dimensional space?

Questions
, , ,
1

There are two sets of points in 3-Dimensional space both contain n points. For example:

A={p1,p2,…pn} pi={pxi,pyi,pzi} (i from 1 to n)

B={q1,q2,…qn} qi={qxi,qyi,qzi} (i from 1 to n)

I want to find a match of points between A and B such that the sum of the distance between each pair of points has minimal value.

How to solve it in Mathematica?

Some days ago, I found the problem can be solved by Kuhn-Munkres algorithm. But it seems that there are no Built-in functions in Mathematica contain Kuhn-Munkres algorithm.

So I turn the problem into an equivalent problem which can be solved by DistanceMatrix and FindMinimumCostFlow.

It did work, but still too slow. So I try to find other ways.

Yesterday, I found Built-in functions: NearestNeighborGraph and Nearest.

If the points in A and B meet certain conditions than NearestNeighborGraph may give right answer in a very fast speed, but not always right or it need further processing.

So I wonder whether the problem can be solved in Mathematica in some clever ways? (Maybe use Nearest, NearestNeighborGraph or some other built-in functions)

2

If a heuristic method is acceptable thenResourceFunction["MaximizeOverPermutations"] can be used for this.

3

Hello Haobo Xia,

my idea here is to use FindShortestTour - and to prohibit any direct step between points of the same set:

aPts = RandomReal[{-1, 0}, {10, 3}];
bPts = RandomReal[{0, 1}, {10, 3}];
pts = Join[aPts, bPts];

sameListQ[a_, b_] := (MemberQ[aPts, a] && MemberQ[aPts, b]) || (MemberQ[bPts, a] && MemberQ[bPts, b])

{length, order} = 
  FindShortestTour[pts, 
   DistanceFunction -> (If[sameListQ[#1, #2], Infinity, EuclideanDistance[#1, #2]] &)];

Graphics3D[{Black, Dashed, Line[pts[[order]]], PointSize[Large], Red, Point[aPts], Green, Point[bPts]}]

twosets.png
twosets.png361×397 6.55 KB

Well, so much for this idea. The problem now seems to be that with this kind of DistanceFunction the function FindShortestTour does not work for more than about 20 points - at least on my system (“12.1.0 for Linux x86 (64-bit) (March 14, 2020)”). Does that help, anyway?

Regards – Henrik

4

Hello Henrik Schachner,

Thank you, it is a creative idea but it need further processing and although it work for large amount of points on my system (Windows 10), it still slower than FindMinimumCostFlow, and sometimes it does not give right answer for my problem.

Regards – Haobo Xia

(*My Function*)
ca[n_] := ca[n] = ConstantArray[0, {n, n}]
PointsMatch[s1_, s2_] := 
 If[s1 == s2, {IdentityMatrix[#], 0}, 
    Extract[FindMinimumCostFlow[
       ArrayFlatten[{{0, 1., 0, 0}, {0, 0, 
          DistanceMatrix[s1, s2, 
           DistanceFunction -> EuclideanDistance], 0}, {0, ca[#], 
          0, -1.}, {0, 0, 0, 0}}], 1, 2*# + 2, 
       "OptimumFlowData"][{"FlowMatrix", "CostValue"}], {{1, 
       Range[2, 1 + #], Range[# + 2, 2*# + 1]}, {2}}]] &@Length[s1]
aPts = RandomReal[{-10, 10}, {8, 3}];
bPts = RandomReal[{-10, 10}, {8, 3}];
pts = Join[aPts, bPts];
(*By MaximizeOverPermutations*)
AbsoluteTiming[
 MaximizeOverPermutations[-Plus @@ 
     MapThread[EuclideanDistance, {aPts, bPts[[#]]}] &, 8]]
(*By FindShortestTour*)
AbsoluteTiming[{length, order} = 
  FindShortestTour[pts, 
   DistanceFunction -> (If[sameListQ[#1, #2], Infinity, 
       EuclideanDistance[#1, #2]] &)];
 MinimalBy[{Plus @@ (EuclideanDistance@(Sequence @@ 
            pts[[#]]) & /@ #), #} & /@ (Partition[#, 2] & /@ {Rest[
       order], order}), First]]
(*By FindMinimumCostFlow*)
AbsoluteTiming[MatrixForm /@ PointsMatch[aPts, bPts]]

Mathematica0.jpg
Mathematica0.jpg1637×708 160 KB

5

Thank you,

It does work, but if it is guaranteed to be optimal, it has to enumerate n! cases.
So it is slower than FindMinimumCostFlow. Thanks anyway.

(*My Function*)
ca[n_] := ca[n] = ConstantArray[0, {n, n}]
PointsMatch[s1_, s2_] := 
 If[s1 == s2, {IdentityMatrix[#], 0}, 
    Extract[FindMinimumCostFlow[
       ArrayFlatten[{{0, 1., 0, 0}, {0, 0, 
          DistanceMatrix[s1, s2, 
           DistanceFunction -> EuclideanDistance], 0}, {0, ca[#], 
          0, -1.}, {0, 0, 0, 0}}], 1, 2*# + 2, 
       "OptimumFlowData"][{"FlowMatrix", "CostValue"}], {{1, 
       Range[2, 1 + #], Range[# + 2, 2*# + 1]}, {2}}]] &@Length[s1]
aPts = RandomReal[{-10, 10}, {8, 3}];
bPts = RandomReal[{-10, 10}, {8, 3}];
pts = Join[aPts, bPts];
(*By MaximizeOverPermutations*)
AbsoluteTiming[
 MaximizeOverPermutations[-Plus @@ 
     MapThread[EuclideanDistance, {aPts, bPts[[#]]}] &, 8]]
(*By FindShortestTour*)
AbsoluteTiming[{length, order} = 
  FindShortestTour[pts, 
   DistanceFunction -> (If[sameListQ[#1, #2], Infinity, 
       EuclideanDistance[#1, #2]] &)];
 MinimalBy[{Plus @@ (EuclideanDistance@(Sequence @@ 
            pts[[#]]) & /@ #), #} & /@ (Partition[#, 2] & /@ {Rest[
       order], order}), First]]
(*By FindMinimumCostFlow*)
AbsoluteTiming[MatrixForm /@ PointsMatch[aPts, bPts]]

Mathematica0.jpg
Mathematica0.jpg1637×708 160 KB