9
$\begingroup$

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.

$\endgroup$
10
  • 1
    $\begingroup$ The brute force method scales as $\mathcal{O}(n^2)$, where $n$ is the number of points. As I understand, even a faster algorithm is desired. In principle, it is possible to speed up by divide and conquer method. Find out the smallest rectangle containing all points and divide it into $10\times10$ cells. Find a pair of furthest non-empty cells, this operation is constant in the number of points $\mathcal{O}(1)$. Divide the 2 cells further and repeat the procedure until only 1 point in each cell remains. The overall scaling should be $\mathcal{O}(\log(n))$ $\endgroup$
    – yarchik
    Aug 2, 2023 at 9:07
  • 1
    $\begingroup$ The algorithm can be more complicated if there are many pairs of furthest points, e.g regular $n$-gon. What do you expect in this case? $\endgroup$
    – yarchik
    Aug 2, 2023 at 9:15
  • $\begingroup$ @yarchik Yes. The actual my case could contain multiple pairs with the exact same distance that is the farthest. $\endgroup$ Aug 2, 2023 at 9:19
  • $\begingroup$ @yarchik, I believe this algorithm fails. Consider the 3x3 cell case for ease of demonstration. if you have a point in the bottom left cell, and then a point at the far edge of the top-center cell, and a point at the closest corner of the top-right cell, the algorithm will decide to focus on bottom-left + top-right, but in reality the top-center point is most distant from the bottom-left point. $\endgroup$
    – Kaia
    Aug 3, 2023 at 16:29
  • 1
    $\begingroup$ Interestingly enough one can measure the maximum distance in v13.3 quite efficiently (seems to scale roughly linearly), but not acquire information on actual points involved: RandomReal[{-100, 100}, {1000000, 2}] // Point // RegionFarthestDistance[#, #] & $\endgroup$
    – kirma
    Aug 29, 2023 at 6:19

5 Answers 5

9
$\begingroup$

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

enter image description here

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

$\endgroup$
5
  • $\begingroup$ Updated to include the definition for reg. Sorry about the earlier omission. $\endgroup$
    – Syed
    Aug 2, 2023 at 6:29
  • $\begingroup$ Is this a property of ConvexHullthat it contains the points with maximal distance? $\endgroup$ Aug 2, 2023 at 9:30
  • $\begingroup$ Similar post on SO and this one. I can't say that I know the rigorous answer @UlrichNeumann. $\endgroup$
    – Syed
    Aug 2, 2023 at 9:36
  • 3
    $\begingroup$ Note that the Union may have unintended behavior. If there are multiple pairs with the farthest distance and they share a point, the information is lost that which two points are paired. $\endgroup$ Aug 2, 2023 at 10:05
  • 2
    $\begingroup$ I think using Subsets instead of Tuples may help. Subsets doesn't create duplicate pairs with the same points but in the opposite ordering. Therefore, deleting duplicated entries using Union will not be needed. $\endgroup$ Aug 2, 2023 at 10:11
12
$\begingroup$

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]

enter image description here

enter image description here

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

enter image description here

$\endgroup$
1
  • $\begingroup$ Is it possible to beat the $\mathcal{O}(n^2)$ scaling in the worst case scenario? $\endgroup$
    – yarchik
    Aug 2, 2023 at 8:55
7
$\begingroup$
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]]]
$\endgroup$
2
  • $\begingroup$ I don't get what your Transpose is for. I thought the last two lines of your answer can be replaced with the following: Part[coordinates, #] & /@ Position[matrix, max] // Column. Is this correct? $\endgroup$ Aug 2, 2023 at 9:51
  • $\begingroup$ @TaikiBessho, take a look at DistanceMatrix docs. Max returns indices of matrix, transpose is used to relate them to list of points, sure, a better approach can be used with part or extract $\endgroup$
    – I.M.
    Aug 2, 2023 at 10:44
6
$\begingroup$

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

enter image description here

Hope it helps!

$\endgroup$
2
  • 2
    $\begingroup$ Nearest tends to be O(n log n) or so, and here it is applied n times, so that's worse than O(n^2) in toto. $\endgroup$ Aug 3, 2023 at 19:33
  • 2
    $\begingroup$ @DanielLichtblau That is of course correct, but QP only asked for a simple solution $\endgroup$ Aug 4, 2023 at 6:25
5
$\begingroup$

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}}} *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.