...die konvexe Hülle von 2D Punkten finden?

Autor: Peter Bone
Homepage: http://www.geocities.com/peter_bone_uk

Kategorie: Mathematik

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;

 

printed from
www.swissdelphicenter.ch
developers knowledge base