Message Boards Message Boards

Find optimal match of points in 3-Dimensional space?

Haobo Xia
Posted 4 years ago

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)
POSTED BY: Haobo Xia
4 Replies
Sort By:

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]}]

enter image description here

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
POSTED BY: Henrik Schachner
Haobo Xia
Haobo Xia
Posted 4 years ago

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]]

enter image description here
POSTED BY: Haobo Xia

If a heuristic method is acceptable thenResourceFunction["MaximizeOverPermutations"] can be used for this.
POSTED BY: Daniel Lichtblau
Haobo Xia
Haobo Xia
Posted 4 years ago

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]]

enter image description here
POSTED BY: Haobo Xia
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract