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