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;