AddFind the convex hull or a set of 2D points
Uses the Graham's scan method to find the convex hull of a set of 2D points.

  1. type
  2.   TPointArray = array of TPoint;
  3.  
  4. // return the boundary points of the convex hull of a set of points
  5. // over-writes the input array - so make a copy first
  6. function FindConvexHull(var APoints : TPointArray) : boolean;
  7. var
  8.   LAngles : array of Real;
  9.   Lindex, LMinY, LMaxX, LPivotIndex : integer;
  10.   LPivot : TPoint;
  11.   LBehind, LInfront : TPoint;
  12.   LRightTurn : boolean;
  13.   LVecPoint : TPoint;
  14. begin
  15.   Result := True;
  16.  
  17.   if length(Points) = 3 then Exit; // already a convex hull
  18.   if length(Points) < 3 then begin // not enough points
  19.     Result := False;
  20.     Exit;
  21.   end;
  22.  
  23.   // find pivot point, which is known to be on the hull
  24.   // point with lowest y - if there are multiple, point with highest x
  25.   LMinY := 10000;
  26.   LMaxX := 10000;
  27.   LPivotIndex := 0;
  28.   for Lindex := 0 to High(APoints) do begin
  29.     if APoints[Lindex].Y = LMinY then begin
  30.       if APoints[Lindex].X > LMaxX then begin
  31.         LMaxX := APoints[Lindex].X;
  32.         LPivotIndex := Lindex;
  33.       end;
  34.     end else if APoints[Lindex].Y < LMinY then begin
  35.       LMinY := APoints[Lindex].Y;
  36.       LMaxX := APoints[Lindex].X;
  37.       LPivotIndex := Lindex;
  38.     end;
  39.   end;
  40.   // put pivot into seperate variable and remove from array
  41.   LPivot := APoints[LPivotIndex];
  42.   APoints[LPivotIndex] := APoints[High(APoints)];
  43.   SetLength(APoints, High(APoints));
  44.  
  45.   // calculate angle to pivot for each point in the array
  46.   SetLength(LAngles, length(APoints));
  47.   for Lindex := 0 to High(APoints) do begin
  48.     LVecPoint.X := LPivot.X - APoints[Lindex].X; // point vector
  49.     LVecPoint.Y := LPivot.Y - APoints[Lindex].Y;
  50.     // reduce to a unit-vector - length 1
  51.     LAngles[Lindex] := LVecPoint.X / Hypot(LVecPoint.X, LVecPoint.Y);
  52.   end;
  53.  
  54.   // sort the points by angle
  55.   QuickSortAngle(APoints, LAngles, 0, High(APoints));
  56.  
  57.   // step through array to remove points that are not part of the convex hull
  58.   Lindex := 1;
  59.   Repeat
  60.     // assign points behind and infront of current point
  61.     if Lindex = 0 then LRightTurn := True else begin
  62.       LBehind := APoints[Lindex-1];
  63.       if Lindex = High(APoints) then LInfront := LPivot
  64.                                 else LInfront := APoints[Lindex + 1];
  65.  
  66.       // work out if we are making a right or left turn using vector product
  67.       if ((LBehind.X-APoints[Lindex].X)*(LInfront.Y-APoints[Lindex].Y))-
  68.          ((LInfront.X-APoints[Lindex].X)*(LBehind.Y-APoints[Lindex].Y)) < 0 then
  69.         LRightTurn := true else
  70.         LRightTurn := false;
  71.     end;
  72.  
  73.     if LRightTurn then begin // point is currently considered part of the hull
  74.       Inc(Lindex); // go to next point
  75.     end else begin // point is not part of the hull
  76.       // remove point from convex hull
  77.       if Lindex = High(APoints) then begin
  78.         SetLength(APoints, High(APoints));
  79.       end else begin
  80.         Move(APoints[Lindex+1], APoints[Lindex], (High(APoints) - Lindex) * SizeOf(TPoint)+1);
  81.         SetLength(APoints, High(APoints));
  82.       end;
  83.  
  84.       Dec(Lindex); // backtrack to previous point
  85.     end;
  86.   until Lindex = High(APoints);
  87.  
  88.   // add pivot back into points array
  89.   SetLength(APoints, length(APoints) + 1);
  90.   APoints[High(APoints)] := LPivot;
  91. end;
  92.  
  93. // sort an array of points by angle
  94. procedure QuickSortAngle(var A: TPointArray ; Angles : array of Real; iLo, iHi: Integer);
  95. var
  96.   Lo, Hi : Integer;
  97.   Mid : real;
  98.   TempPoint : TPoint;
  99.   TempAngle : Real;
  100. begin
  101.   Lo := iLo;
  102.   Hi := iHi;
  103.   Mid := Angles[(Lo + Hi) shr 1];
  104.   repeat
  105.     while Angles[Lo] < Mid do Inc(Lo);
  106.     while Angles[Hi] > Mid do Dec(Hi);
  107.     if Lo  Hi;
  108.   // perform quicksorts on subsections
  109.   if Hi > iLo then QuickSortAngle(A, Angles, iLo, Hi);
  110.   if Lo < iHi then QuickSortAngle(A, Angles, Lo, iHi);
  111. end;
  • Author:
  • URL:
    http://www.geocities.com/peter_bone_uk
  • Date added:
    03 March, 2006
  • Views:
    7578
Latest News
Submit News Form Past News
Latest Forum Entries