Sources
Delphi Russian Knowledge Base
DRKB - это самая большая и удобная в использовании база знаний по Дельфи в рунете, составленная Виталием Невзоровым

Find the convex hull of 2D points

01.01.2007
uses
 
Math;
 
type
 
TPointArray = array of TPoint;
 
 
TPointFloat = record
    X
: Real;
    Y
: Real;
 
end;
 
 
// return the boundary points of the convex hull of a set of points using Grahams scan
 
// over-writes the input array - so make a copy first
function FindConvexHull(var APoints: TPointArray): Boolean;
var
 
LAngles: array of Real;
 
Lindex, LMinY, LMaxX, LPivotIndex: Integer;
 
LPivot: TPoint;
 
LBehind, LInfront: TPoint;
 
LRightTurn: Boolean;
 
LVecPoint: TPointFloat;
begin
 
Result := True;
 
 
if Length(Points) = 3 then Exit; // already a convex hull
 
if Length(Points) < 3 then
 
begin // not enough points
   
Result := False;
   
Exit;
 
end;
 
 
// find pivot point, which is known to be on the hull
 
// point with lowest y - if there are multiple, point with highest x
 
LMinY := 1000;
 
LMaxX := 1000;
 
LPivotIndex := 0;
 
for Lindex := 0 to High(APoints) do
 
begin
   
if APoints[Lindex].Y = LMinY then
   
begin
     
if APoints[Lindex].X > LMaxX then
     
begin
       
LMaxX := APoints[Lindex].X;
       
LPivotIndex := Lindex;
     
end;
   
end
   
else if APoints[Lindex].Y < LMinY then
   
begin
     
LMinY := APoints[Lindex].Y;
     
LMaxX := APoints[Lindex].X;
     
LPivotIndex := Lindex;
   
end;
 
end;
 
// put pivot into seperate variable and remove from array
 
LPivot := APoints[LPivotIndex];
 
APoints[LPivotIndex] := APoints[High(APoints)];
 
SetLength(APoints, High(APoints));
 
 
// calculate angle to pivot for each point in the array
 
// quicker to calculate dot product of point with a horizontal comparison vector
 
SetLength(LAngles, Length(APoints));
 
for Lindex := 0 to High(APoints) do
 
begin
   
LVecPoint.X := LPivot.X - APoints[Lindex].X; // point vector
   
LVecPoint.Y := LPivot.Y - APoints[Lindex].Y;
   
// reduce to a unit-vector - length 1
   
LAngles[Lindex] := LVecPoint.X / Hypot(LVecPoint.X, LVecPoint.Y);
 
end;
 
 
// sort the points by angle
 
QuickSortAngle(APoints, LAngles, 0, High(APoints));
 
 
// step through array to remove points that are not part of the convex hull
 
Lindex := 1;
  repeat
   
// assign points behind and infront of current point
   
if Lindex = 0 then LRightTurn := True
   
else
   
begin
     
LBehind := APoints[Lindex - 1];
     
if Lindex = High(APoints) then LInfront := LPivot
     
else
       
LInfront := APoints[Lindex + 1];
 
     
// work out if we are making a right or left turn using vector product
     
if ((LBehind.X - APoints[Lindex].X) * (LInfront.Y - APoints[Lindex].Y)) -
       
((LInfront.X - APoints[Lindex].X) * (LBehind.Y - APoints[Lindex].Y)) < 0 then
       
LRightTurn := True
     
else
       
LRightTurn := False;
   
end;
 
   
if LRightTurn then
   
begin // point is currently considered part of the hull
     
Inc(Lindex); // go to next point
   
end
   
else
   
begin // point is not part of the hull
     
// remove point from convex hull
     
if Lindex = High(APoints) then
     
begin
       
SetLength(APoints, High(APoints));
     
end
     
else
     
begin
       
Move(APoints[Lindex + 1], APoints[Lindex],
       
(High(APoints) - Lindex) * SizeOf(TPoint) + 1);
       
SetLength(APoints, High(APoints));
     
end;
 
     
Dec(Lindex); // backtrack to previous point
   
end;
 
until Lindex = High(APoints);
 
 
// add pivot back into points array
 
SetLength(APoints, Length(APoints) + 1);
 
APoints[High(APoints)] := LPivot;
end;
 
// sort an array of points by angle
procedure
QuickSortAngle(var A: TPointArray; Angles: array of Real; iLo, iHi: Integer);
var
 
Lo, Hi: Integer;
 
Mid: Real;
 
TempPoint: TPoint;
 
TempAngle: Real;
begin
 
Lo  := iLo;
 
Hi  := iHi;
 
Mid := Angles[(Lo + Hi) div 2];
  repeat
   
while Angles[Lo] < Mid do Inc(Lo);
   
while Angles[Hi] > Mid do Dec(Hi);
   
if Lo <= Hi then
   
begin
     
// swap points
     
TempPoint := A[Lo];
      A
[Lo] := A[Hi];
      A
[Hi] := TempPoint;
     
// swap angles
     
TempAngle := Angles[Lo];
     
Angles[Lo] := Angles[Hi];
     
Angles[Hi] := TempAngle;
     
Inc(Lo);
     
Dec(Hi);
   
end;
 
until Lo > Hi;
 
// perform quicksorts on subsections
 
if Hi > iLo then QuickSortAngle(A, Angles, iLo, Hi);
 
if Lo < iHi then QuickSortAngle(A, Angles, Lo, iHi);
end;

Взято с сайта https://www.swissdelphicenter.ch/en/tipsindex.php