source: trunk/Packages/Common/UGeometry.pas

unit UGeometry;

{$mode delphi}

interface

uses
  Classes, SysUtils, Math;

type
  TPointArray = array of TPoint;

  { TLine }

  TLine = record
  private
    function GetDistance: Double;
    procedure SetDistance(AValue: Double);
  public
    P1: TPoint;
    P2: TPoint;
    function Create(const P1, P2: TPoint): TLine;
    function GetMiddle: TPoint;
    function GetAngle: Double;
    function GetSize: TPoint;
    function ToRect: TRect;
    function DotProduct: Double;
    procedure Rotate(const Angle: Double);
    class operator Equal(const A, B: TLine): Boolean;
    property Distance: Double read GetDistance write SetDistance;
  end;

  { TPolygon }

  TPolygon = record
    Points: TPointArray;
    function IsPointInside(const P: TPoint): Boolean;
    function Create(const Points: TPointArray): TPolygon; overload;
    function Create(const Rect: TRect): TPolygon; overload;
    procedure AddPoint(const P: TPoint);
    procedure Clear;
    procedure CutLine(const Vector: TLine; const PointInside: TPoint);
  end;

function Distance(const P1, P2: TPoint): Integer;
function Dot(const P1, P2: TPoint): Double;
function AddPoint(const P1, P2: TPoint): TPoint;
function SubPoint(const P1, P2: TPoint): TPoint;
function PointToLineDistance(const P, V, W: TPoint): Integer;
function ComparePoint(const P1, P2: TPoint): Boolean;
function RotatePoint(const Center, P: TPoint; Angle: Double): TPoint;
function RotatePoints(const Center: TPoint; P: TPointArray; Angle: Double): TPointArray;
function LineIntersect(const LineA, LineB: TLine; out Intersection: TPoint): Boolean;
function ArcTan2Point(const Point: TPoint): Float;
function ArcTanPoint(const Point: TPoint): Float;
function RectEquals(const A, B: TRect): Boolean;
function RectEnlarge(const Rect: TRect; Value: Integer): TRect;
function ShiftRect(const ARect: TRect; Delta: TPoint): TRect;
function PointsToRect(const P1, P2: TPoint): TRect;
function PointInRect(const P: TPoint; aRect: TRect): Boolean;
function HalfDistancePoint(const P1, P2: TPoint): TPoint;
function NormalizeAngle(const Angle: Double): Double;
function SubAngle(A1, A2: Double): Double;

implementation

function Distance(const P1, P2: TPoint): Integer;
begin
  Result := Trunc(Sqrt(Sqr(P2.X - P1.X) + Sqr(P2.Y - P1.Y)));
end;

function Dot(const P1, P2: TPoint): Double;
begin
  Result := P1.X * P2.X + P1.Y * P2.Y;
end;

function AddPoint(const P1, P2: TPoint): TPoint;
begin
  Result.X := P1.X + P2.X;
  Result.Y := P1.Y + P2.Y;
end;

function SubPoint(const P1, P2: TPoint): TPoint;
begin
  Result.X := P1.X - P2.X;
  Result.Y := P1.Y - P2.Y;
end;

function PointToLineDistance(const P, V, W: TPoint): Integer;
var
  l2, t: Double;
  tt: TPoint;
begin
  // Return minimum distance between line segment vw and point p
  L2 := Distance(V, W); // i.e. |w-v|^2 -  avoid a sqrt
  L2 := Power(l2, 2);
  if L2 = 0 then begin
    Result := Distance(P, V);   // v == w case
    Exit;
  end;
  // Consider the line extending the segment, parameterized as v + t (w - v).
  // We find projection of point p onto the line.
  // It falls where t = [(p-v) . (w-v)] / |w-v|^2
  T := Dot(SubPoint(P, V), SubPoint(W, V)) / L2;
  if T < 0 then begin
    Result := Distance(P, V);       // Beyond the 'v' end of the segment
    exit;
  end
  else if T > 1 then begin
    Result := Distance(P, W);  // Beyond the 'w' end of the segment
    Exit;
  end;
  TT.X := Trunc(V.X + T * (W.X - V.X));
  TT.Y := Trunc(V.Y + T * (W.Y - V.Y));
  Result := Distance(P, TT);
end;

function ComparePoint(const P1, P2: TPoint): Boolean;
begin
  Result := (P1.X = P2.X) and (P1.Y = P2.Y);
end;

function RotatePoint(const Center, P: TPoint; Angle: Double): TPoint;
var
  D: TPoint;
begin
  D := Point(P.X - Center.X, P.Y - Center.Y);
  Result := Point(Center.X + Round(D.X * Cos(Angle) - D.Y * Sin(Angle)),
    Center.Y + Round(D.X * Sin(Angle) + D.Y * Cos(Angle)));
end;

function RotatePoints(const Center: TPoint; P: TPointArray; Angle: Double): TPointArray;
var
  I: Integer;
begin
  SetLength(Result, Length(P));
  for I := 0 to High(P) do
    Result[I] := RotatePoint(Center, P[I], Angle);
end;

function LineIntersect(const LineA, LineB: TLine; out Intersection: TPoint): Boolean;
Var
  LDetLineA, LDetLineB, LDetDivInv: Double;
  LDiffLA, LDiffLB: TPoint;
  D: Integer;
begin
  if ComparePoint(LineA.P1, LineA.P2) or ComparePoint(LineB.P1, LineB.P2) then begin
    Result := False;
    Exit;
  end;
  LDetLineA := LineA.P1.X * LineA.P2.Y - LineA.P1.Y * LineA.P2.X;
  LDetLineB := LineB.P1.X * LineB.P2.Y - LineB.P1.Y * LineB.P2.X;

  LDiffLA := SubPoint(LineA.P1, LineA.P2);
  LDiffLB := SubPoint(LineB.P1, LineB.P2);

  D := (LDiffLA.X * LDiffLB.Y) - (LDiffLA.Y * LDiffLB.X);
  if D = 0 then begin
    // Parallel lines without intersection
    Result := False;
    Exit;
  end;
  LDetDivInv := 1 / D;

  Intersection.X := Round(((LDetLineA * LDiffLB.X) - (LDiffLA.X * LDetLineB)) * LDetDivInv);
  Intersection.Y := Round(((LDetLineA * LDiffLB.Y) - (LDiffLA.Y * LDetLineB)) * LDetDivInv);
  Result := True;
end;

function ArcTan2Point(const Point: TPoint): Float;
begin
  Result := ArcTan2(Point.Y, Point.X);
end;

function ArcTanPoint(const Point: TPoint): Float;
begin
  if Point.Y = 0 then Result := Infinity
    else Result := ArcTan(Point.X / Point.Y);
end;

function RectEquals(const A, B: TRect): Boolean;
begin
  Result := (A.Left = B.Left) and (A.Top = B.Top) and
    (A.Right = B.Right) and (A.Bottom = B.Bottom);
end;

function RectEnlarge(const Rect: TRect; Value: Integer): TRect;
begin
  Result.Left := Rect.Left - Value;
  Result.Right := Rect.Right + Value;
  Result.Top := Rect.Top - Value;
  Result.Bottom := Rect.Bottom + Value;
end;

function ShiftRect(const ARect: TRect; Delta: TPoint): TRect;
begin
  Result := Rect(ARect.Left + Delta.X, ARect.Top + Delta.Y,
    ARect.Right + Delta.X, ARect.Bottom + Delta.Y);
end;

function PointsToRect(const P1, P2: TPoint): TRect;
begin
  if P1.X < P2.X then Result.Left := P1.X else Result.Left := P2.X;
  if P1.Y < P2.Y then Result.Top := P1.Y else Result.Top := P2.Y;
  if P1.X > P2.X then Result.Right := P1.X else Result.Right := P2.X;
  if P1.Y > P2.Y then Result.Bottom := P1.Y else Result.Bottom := P2.Y;
end;

function PointInRect(const P: TPoint; aRect: TRect): Boolean;
begin
  Result := (P.X >= aRect.Left) and (P.X <= aRect.Right) and
    (P.Y >= aRect.Top) and (P.Y <= aRect.Bottom);
end;

function HalfDistancePoint(const P1, P2: TPoint): TPoint;
begin
  Result := Point(P1.X + (P2.X - P1.X) div 2, P1.Y + (P2.Y - P1.Y) div 2)
end;

function NormalizeAngle(const Angle: Double): Double;
begin
  if Angle < 0 then Result := Angle + (Trunc(Angle / (2 * Pi)) + 1) * (2 * Pi)
  else if Angle > 2 * Pi then Result := Angle - Trunc(Angle / (2 * Pi)) * (2 * Pi)
  else Result := Angle;
end;

function SubAngle(A1, A2: Double): Double;
begin
  A1 := NormalizeAngle(A1);
  A2 := NormalizeAngle(A2);
  if A1 < A2 then Result := A1 + 2 * Pi - A2
    else Result := A1 - A2;
end;

{ TPolygon }

function TPolygon.IsPointInside(const P: TPoint): Boolean;
var
  I, J: Integer;
begin
  Result := False;
  J := High(Points);
  for I := Low(Points) to High(Points) do begin
    if ((Points[I].Y <= P.Y) and (P.Y < Points[J].Y)) or
      ((Points[J].Y <= P.Y) and (P.Y < Points[I].Y)) then
      begin
      if (P.X < (Points[J].X - Points[I].X) *
         (P.Y - Points[I].Y) /
         (Points[J].Y - Points[I].Y) + Points[I].X) then
          Result := not Result;
      end;
    J := I;
  end;
end;

function PtInPoly(const Points: array of TPoint; Pos: TPoint): Boolean;
var
  Count, K, J : Integer;
begin
  Result := False;
  Count := Length(Points) ;
  J := Count - 1;
  for K := 0 to Count - 1 do begin
  if ((Points[K].Y <= Pos.Y) and (Pos.Y < Points[J].Y)) or
    ((Points[J].Y <= Pos.Y) and (Pos.Y < Points[K].Y)) then
    begin
    if (Pos.X < (Points[j].X - Points[K].X) *
       (Pos.Y - Points[K].Y) /
       (Points[j].Y - Points[K].Y) + Points[K].X) then
        Result := not Result;
    end;
    J := K;
  end;
end;


function TPolygon.Create(const Points: TPointArray): TPolygon;
var
  I: Integer;
begin
  SetLength(Result.Points, Length(Points));
  for I := 0 to Length(Points) - 1 do
    Result.Points[I] := Points[I];
end;

function TPolygon.Create(const Rect: TRect): TPolygon;
begin
  SetLength(Result.Points, 4);
  Result.Points[0] := Point(Rect.Left, Rect.Top);
  Result.Points[1] := Point(Rect.Right, Rect.Top);
  Result.Points[2] := Point(Rect.Right, Rect.Bottom);
  Result.Points[3] := Point(Rect.Left, Rect.Bottom);
end;

procedure TPolygon.AddPoint(const P: TPoint);
begin
  SetLength(Points, Length(Points) + 1);
  Points[Length(Points) - 1] := P;
end;

procedure TPolygon.Clear;
begin
  SetLength(Points, 0);
end;

procedure TPolygon.CutLine(const Vector: TLine; const PointInside: TPoint);
var
  I: Integer;
  PointsChecked: Integer;
  L1, L2: TLine;
  Intersection: TPoint;
  NewPoly: TPolygon;
  NewPolygonStarted: Boolean;
  Success: Boolean;
begin
  NewPoly.Clear;
  Success := False;
  NewPolygonStarted := False;
  I := 0;
  PointsChecked := 0;
  L1 := Vector;
  L1.Rotate(Pi / 2);
  if Length(Points) > 0 then
  while True do begin
    L2 := TLine.Create(Points[I], Points[(I + 1) mod Length(Points)]);
    if LineIntersect(L1, L2, Intersection) then
    if PointInRect(Intersection, L2.ToRect) then begin
      if not NewPolygonStarted then begin
        // Crossing line, start new polygon
        NewPoly.Clear;
        NewPoly.AddPoint(Intersection);
        NewPolygonStarted := True;
      end else begin
        // Crossing line, end polygon. If point NewPolygonStarted, the use polygon as result
        NewPoly.AddPoint(Points[I]);
        NewPoly.AddPoint(Intersection);
        if NewPoly.IsPointInside(PointInside) then begin
          Success := True;
          Break;
        end else begin
          NewPoly.Clear;
          NewPoly.AddPoint(Intersection);
          NewPolygonStarted := True;
        end;
      end;
    end else
      NewPoly.AddPoint(Points[I]);
    I := (I + 1) mod Length(Points);
    Inc(PointsChecked);
    if PointsChecked > 2 * Length(Points) then Break;
  end;
  if Success then Points := NewPoly.Points;
end;

{ TLine }

function TLine.GetDistance: Double;
begin
  Result := Sqrt(Sqr(P2.X - P1.X) + Sqr(P2.Y - P1.Y));
end;

procedure TLine.SetDistance(AValue: Double);
var
  Angle: Double;
begin
  Angle := GetAngle;
  P2 := Point(Round(P1.X + Cos(Angle) * AValue),
    Round(P1.Y + Sin(Angle) * AValue));
end;

function TLine.Create(const P1, P2: TPoint): TLine;
begin
  Result.P1 := P1;
  Result.P2 := P2;
end;

function TLine.GetMiddle: TPoint;
begin
  Result := Point(P1.X + (P2.X - P1.X) div 2, P1.Y + (P2.Y - P1.Y) div 2);
end;

function TLine.GetAngle: Double;
begin
  Result := ArcTan2(P2.Y - P1.Y, P2.X - P1.X);
end;

function TLine.GetSize: TPoint;
begin
  Result := Point(P2.X - P1.X, P2.Y - P1.Y);
end;

function TLine.ToRect: TRect;
begin
  if P1.X < P2.X then Result.Left := P1.X else Result.Left := P2.X;
  if P1.Y < P2.Y then Result.Top := P1.Y else Result.Top := P2.Y;
  if P1.X > P2.X then Result.Right := P1.X else Result.Right := P2.X;
  if P1.Y > P2.Y then Result.Bottom := P1.Y else Result.Bottom := P2.Y;
end;

function TLine.DotProduct: Double;
begin
  Result := P1.X * P2.X + P1.Y * P2.Y;
end;

procedure TLine.Rotate(const Angle: Double);
begin
  P2 := RotatePoint(P1, P2, Angle);
end;

class operator TLine.Equal(const A, B: TLine): Boolean;
begin
  Result := ComparePoint(A.P1, B.P1) and ComparePoint(A.P2, B.P2);
end;

end.