Find the farthest pairs of two points in 2-dim

Question

I wonder if Mathematica can find the pairs of two points that are the farthest in a given list.

Let's say a list, XYlist, of points with (x,y) coordinates.

XYlist = Flatten[Table[{x, y}, {y, 0, 100}, {x, 0, 100}], 1]

We, humans, can easily find out that the following two pairs have the farthest distance in the XYlist by, for example, plotting the points:

{{0,0},{100,100}}
{{100,0},{0,100}}

ListPlot[XYlist, AspectRatio -> Automatic]

I, however, have not had any good idea of how to instruct Mathematica to find out the pairs from the list.

A close way I've got is using BoundingRegion.

br = BoundingRegion[XYlist, "MinDisk"]
(* Disk[{50, 50}, 50 Sqrt[2]] *)
max = 2*br[[2]]
(* 100 Sqrt[2] *)

It gives the minimum area Disk that includes all the points. In the above example, this method can correctly find out the farthest distance between the points of the list. However, it does not tell which pairs give the farthest distance, and furthermore, it does not work if the two points of the farthest pairs are not on a diameter of the returned Disk.

Another way is using Subsets and creating all the possible pairs at first, then checking the distance one by one:

ss = Subsets[XYlist, {2}]
nl = Parallelize[Norm /@ Subtract @@@ ss]
max = Max[nl]
pos = Position[nl, max]
pairs = Extract[ss, pos]
(* {{{0, 0}, {100, 100}}, {{100, 0}, {0, 100}}} *)

However, this method is very memory-consuming and slow, therefore it is not realistic if the number of points gets large.

Because BoundingRegion can find the minimum area Disk in a flush, which I think calculates fundamentally similar things, I can't help expecting there exist some very fast methods.

In the case that there are multiple pairs with the farthest distance, I want to get all the pairs.

Do you guys have any good ideas? Thanks in advance.

Answer 1

Let's say you have random points:

pts = RandomReal[{-100, 100}, {10000, 2}];

It would be useful to find the ConvexHullRegion and find distances between the points at the boundaries. This can be calculated in the following ways, as both outputs are identical. We will use bpts.

reg = ConvexHullRegion[pts]
bpts = PolygonCoordinates[reg] // Sort
cpts = Sequence @@@ MeshPrimitives[reg, 0] // Sort

farthest = 
 Union @@ MaximalBy[Tuples[bpts, 2], Apply@EuclideanDistance]

{{-99.1393, 98.7636}, {99.9366, -98.131}}

Visualization:

Region[reg
 , Epilog -> {Black, Point[pts]
   , AbsolutePointSize[8], Red
   , Point@farthest
   }
 ]

I hope you can find some usable pointers from this answer.

Answer 2

upd: fixed the error of farthestPointsDistance case. The final code is so verbose, whereas the result is definitely true.

---

You can implement your own functions closestPointsDistance and farthestPointsDistance.

farthestPointsDistance✅

Rotating Calipers Algorithm, $\text{average }\mathcal{O}(n \log n)$

Clear["Global`*"];

dist[p1_, p2_] := Norm[p2 - p1];


(*really naive, so slow sortPointsCounterClockwise*)
sortPointsCounterClockwise[pts_] := 
  Module[{center = Mean[pts], angles}, 
   angles = ArcTan @@ (# - center) & /@ pts;
   pts[[Ordering[angles]]]];

convexHull[pts_] := 
  sortPointsCounterClockwise@MeshCoordinates@ConvexHullMesh@pts;

TriangleArea[p1_, p2_, p3_] := 
  Abs[(p1[[1]]*(p2[[2]] - p3[[2]]) + p2[[1]]*(p3[[2]] - p1[[2]]) + 
      p3[[1]]*(p1[[2]] - p2[[2]]))/2];



(*Proposed by Preparata and Shamos in 1985 that avoided calculation of angles *)
GetAllAntiPodalPairs[p_List] := 
 Module[{n = Length[p], i0, i = 1, j, j0}, i0 = n;
  j = i + 1;
  While[TriangleArea[p[[i]], p[[Mod[i + 1, n, 1]]], 
     p[[Mod[j + 1, n, 1]]]] > 
    TriangleArea[p[[i]], p[[Mod[i + 1, n, 1]]], p[[Mod[j, n, 1]]]], 
   j = Mod[j + 1, n, 1]];
  j0 = j;
  Reap[While[i != j0, i = Mod[i + 1, n, 1];
     Sow[{i, j}];
     While[
      TriangleArea[p[[i]], p[[Mod[i + 1, n, 1]]], 
        p[[Mod[j + 1, n, 1]]]] > 
       TriangleArea[p[[i]], p[[Mod[i + 1, n, 1]]], p[[Mod[j, n, 1]]]],
       j = Mod[j + 1, n, 1];
      If[{i, j} != {j0, i0}, Sow[{i, j}]];];
     If[TriangleArea[p[[j]], p[[Mod[i + 1, n, 1]]], 
        p[[Mod[j + 1, n, 1]]]] == 
       TriangleArea[p[[i]], p[[Mod[i + 1, n, 1]]], p[[Mod[j, n, 1]]]],
       If[{i, j} != {j0, i0}, Sow[{i, Mod[j + 1, n, 1]}], 
        Sow[{Mod[i + 1, n, 1], j}]];];]][[2, 1]]]





pts = RandomReal[{-100, 100}, {50, 2}];
polygonPoints = sortPointsCounterClockwise@convexHull@pts;
antiPodalPairs = GetAllAntiPodalPairs@polygonPoints;
numberOfPolygonPoints = Length[polygonPoints];
calculateDistance =(((dist @@ {polygonPoints[[Mod[#1, numberOfPolygonPoints, 1]]], polygonPoints[[Mod[#2, numberOfPolygonPoints, 1]]]}) &) @@ # &);
calculateDistance /@ antiPodalPairs // Max

(* bruteForce to verify the result. *)
dist @@ # & /@ Subsets[pts, {2}] // Max



GraphicsGrid[
 Partition[
  Table[Graphics[{(*Draw the polygon*){EdgeForm[Directive[Black]], 
      FaceForm[None], Polygon[polygonPoints]},(*Draw the anti-
     podal pairs*){Red, PointSize[Large], 
      Point[{polygonPoints[[antiPodalPairs[[i, 1]]]], 
        polygonPoints[[antiPodalPairs[[i, 
            2]]]]}]},(*Draw lines connecting anti-
     podal pairs*){Blue, Dashed, 
      Line[{polygonPoints[[antiPodalPairs[[i, 1]]]], 
        polygonPoints[[antiPodalPairs[[i, 2]]]]}]}}], {i, 
    Length[antiPodalPairs]}], UpTo[2] (*Maximum 2 plots in a row*)], 
 Frame -> All]

maxPointsPairs = 
  MaximalBy[antiPodalPairs, calculateDistance ];

maxPoints = 
  Part[{polygonPoints[[#1]], polygonPoints[[#2]]} & @@@ 
    maxPointsPairs, 1];

(*Visualization*)ListPlot[pts, 
 Epilog -> {Red, PointSize[Large], Point[maxPoints], Green, 
   Line[maxPoints]}, AspectRatio -> 1]

closestPointsDistance✅

Divide-and-conquer Algorithm, $\text{average }\mathcal{O}(n \log n)$

euclideanDistance[p1_, p2_] := EuclideanDistance[p1, p2]

bruteForce[points_] := 
 Module[{minDistance = Infinity, dist, minPoints}, 
  Do[dist = euclideanDistance[points[[i]], points[[j]]];
   If[dist < minDistance, minDistance = dist; 
    minPoints = {points[[i]], points[[j]]}], {i, 
    Length[points] - 1}, {j, i + 1, Length[points]}];
  {minDistance, minPoints}]

stripClosest[strip_, d_, pair_] := 
 Module[{minDistance = d, dist, minPoints = pair}, 
  Do[dist = euclideanDistance[strip[[i]], strip[[j]]];
   If[dist < minDistance, minDistance = dist; 
    minPoints = {strip[[i]], strip[[j]]}], {i, Length[strip] - 1}, {j,
     i + 1, Min[i + 7, Length[strip]]}];
  {minDistance, minPoints}]

closestPointsDistance[points_] := 
 Module[{mid, midPoint, dl, dr, d, pair, strip, result}, 
  If[Length[points] <= 3, Return[bruteForce[points]]];
  mid = Ceiling[Length[points]/2];
  midPoint = points[[mid]];
  dl = closestPointsDistance[Take[points, mid]];
  dr = closestPointsDistance[Drop[points, mid]];
  {d, pair} = If[First[dl] < First[dr], dl, dr];
  strip = Select[points, Abs[#[[1]] - midPoint[[1]]] < d &];
  result = stripClosest[strip, d, pair];
  result]

minDistPoints[points_List] := 
 Module[{sortedPoints}, sortedPoints = SortBy[points, First];
  closestPointsDistance[sortedPoints]]

pts = RandomReal[{-100, 100}, {50, 2}];
{minDist, minPoints} = minDistPoints[pts];
Print[minDist];

(*Visualization*)
ListPlot[pts, 
 Epilog -> {Red, PointSize[Large], Point[minPoints], Green, 
   Line[minPoints]}, AspectRatio -> 1]

(* bruteForce to verify the result. *)
euclideanDistance @@ # & /@ Subsets[pts, {2}] // Min

Answer 3

XYlist = Flatten[Table[{x, y}, {y, 0, 100}, {x, 0, 100}], 1];
region = ConvexHullMesh[XYlist] ;
coordinates = MeshCoordinates[region] ;
matrix = DistanceMatrix[coordinates] ;
max = Max[matrix] ;

{i, j}= Transpose[Position[matrix, max] ] ;
Transpose[{coordinates[[i]], coordinates[[j]]}] // Column
(* or *)
Partition[coordinates[[Flatten[Position[matrix, max]]]], Last[Dimensions[coordinates]]]

Answer 4

Nearestgives a more direct way I think. Therefor we only use the last point found by pointwise Nearest-search

pts = RandomReal[{-100, 100}, {100 , 2}];
maxpair =Map[(np = Nearest[pts, #, All][[-1]]; {#, np,Sqrt[(np - #) . (np - #)]}) &, pts];

maxpair contains a list of point-pairs with maximal distance.

max=Select[maxpair, #[[3]] == Max[maxpair[[All, 3]]] &]
Graphics[{Point[pts], Red, Line[max[[All, {1, 2}]]]}]

Hope it helps!

Answer 5

This is a variation of @Syed's answer. Distance computation of convex hull boundary points is improved by using DistanceMatrix - which is very well optimised, even if half of the matrix is unnecessary. This optimisation is beneficial when there's lots of boundary points (try generating input with N@CirclePoints[1000] for instance). Both methods suffer from the fact that memory usage grows by the square of the number of the boundary points...

With[{bpts = PolygonCoordinates@ConvexHullMesh@pts}, 
 DistanceMatrix[bpts, DistanceFunction -> SquaredEuclideanDistance] //
   bpts[[FirstPosition[#, Max@Flatten@#, {2}]]] &]

To avoid this memory usage issue, one can attempt to implement a function which doesn't use memory unnecessarily for storage of matrices or tuple lists:

farthestPair[pts_] :=
 ParallelTable[
   pts[[i + 1 ;;]] - ConstantArray[pts[[i]], Length@pts - i] //
     Map@Norm //
    With[{pos = First@PositionLargest@#},
      {#[[pos]], pos + i}] &,
   {i, Length@pts - 1}, Method -> "CoarsestGrained"] //
  With[{pos = First@PositionLargest[#[[All, 1]]]},
    pts[[{pos, #[[pos, 2]]}]]] &

farthestPair@PolygonCoordinates@ConvexHullMesh@pts

It must be stated neither DistanceMatrix nor the above can't compete performance-wise with the Rotating Calipers algorithm since its complexity is $\text{average }\mathcal{O}(n \log n)$, and this is obviously $\mathcal{O}(n^2)$, but at least this implementation is easy to understand and could also trivially be generalised to three dimensions. Also it should be noted that the $n$ is the number of boundary points of the convex hull, which is typically small in comparison to all input points.

---

@138_Aspen's Rotating Calipers is much faster with large number of boundary points though, so I decided to write it to use a compiled function:

farthestPoints =
  With[{rotatingCalipers =
     FunctionCompile[
      Function[Typed[p, "PackedArray"::["Real64", 2]],
       Module[
        {area =
          Abs[#1[[1]] (#2[[2]] - #3[[2]]) +
              #2[[1]] (#3[[2]] - #1[[2]]) +
              #3[[1]] (#1[[2]] - #2[[2]])] &,
         res = {}},
        Module[{n = Length[p], i0, i = 1, j, j0}, i0 = n;
         j = i + 1;
         While[area[p[[i]], p[[Mod[i + 1, n, 1]]], 
            p[[Mod[j + 1, n, 1]]]] > 
           area[p[[i]], p[[Mod[i + 1, n, 1]]], p[[Mod[j, n, 1]]]], 
          j = Mod[j + 1, n, 1]];
         j0 = j;
         While[i != j0, i = Mod[i + 1, n, 1];
          AppendTo[res, {i, j}];
          While[area[p[[i]], p[[Mod[i + 1, n, 1]]], 
             p[[Mod[j + 1, n, 1]]]] > 
            area[p[[i]], p[[Mod[i + 1, n, 1]]], p[[Mod[j, n, 1]]]], 
           j = Mod[j + 1, n, 1];
           If[{i, j} != {j0, i0}, AppendTo[res, {i, j}]];];
          If[area[p[[j]], p[[Mod[i + 1, n, 1]]], 
             p[[Mod[j + 1, n, 1]]]] == 
            area[p[[i]], p[[Mod[i + 1, n, 1]]], p[[Mod[j, n, 1]]]], 
           If[{i, j} != {j0, i0}, 
             AppendTo[res, {i, Mod[j + 1, n, 1]}], 
             AppendTo[res, {Mod[i + 1, n, 1], j}]];];];
         res]]]]},
   
   With[{p = MeshPrimitives[ConvexHullMesh[#], 2][[1, 1]]},
     p[[#]] & /@ rotatingCalipers[p] // 
       TakeLargestBy[Apply@SquaredEuclideanDistance, 1] //
      First] &];

It finds farthest points in a 10000 point boundary set on a circle in half a second (DistanceMatrix takes 9 seconds on my laptop):

(Timing@*farthestPoints)@RandomSample@N@CirclePoints[10000]

(* {0.480575, {{-0.895386, -0.445292}, {0.895386, 0.445292}}} *)