Реализация

{ Модуль для работы с геометрическими объектами }
{ Copyright (c) Антон Лапунов, 1998-2001 }
{ Версия для сайта ОНЗИ http://onzi.narod.ru }

Unit Geometry; {$N+,G+}

Interface

Type Real   = Extended;                       { Вещественное число }
     Point  = Record x,y : Real End;          { Точка }
     Line   = Record A,B,C : Real End;        { Прямая }
     Circle = Record O : Point; R : Real End; { Окружность }
     PClass = (_On,_Inside,_Outside,_Negative,_Positive);

Const Eps : Real = 1e-10;                     { "Почти" ноль }
      Inf = 1e+10;                            { "Бесконечность" }
      PInf : Point = (x:Inf; y:Inf);          { "Бесконечная" точка }
      PClassStr : Array [PClass] Of String =
        ('On','Inside','Outside','Negative','Positive');

Function Eq(Const a,b : Real) : Boolean;
Function Ne(Const a,b : Real) : Boolean;
Function Le(Const a,b : Real) : Boolean;
Function Lt(Const a,b : Real) : Boolean;
Function Ge(Const a,b : Real) : Boolean;
Function Gt(Const a,b : Real) : Boolean;
Function Btw(Const x,a,b : Real) : Boolean;
Function BtwS(Const x,a,b : Real) : Boolean;
Function Det(Const a,b,c,d : Real) : Real;
Function ArcTan2(Const y,x : Real) : Real;
Function PolarAngle(Const x,y : Real) : Real;
Function ADiff(Const a,b : Real) : Real;

Procedure PSet(Var P : Point; Const x,y : Real);
Procedure LSet(Var L : Line; Const A,B,C : Real);
Procedure CSet(Var C : Circle; Const x,y,R : Real);

Function PStr(Const P : Point; w,d : Integer) : String;
Function PBtw(Const P,P1,P2 : Point) : Boolean;
Function PBtwS(Const P,P1,P2 : Point) : Boolean;
Procedure PMoveTo(Const P1,P2 : Point; Const k : Real; Var P : Point);
Procedure PLProject(Const P : Point; Const L : Line; Var P1 : Point);
Procedure PLine(Const A,B : Real; Const P : Point; Var L : Line);
Procedure P2Line(Const P1,P2 : Point; Var L : Line);
Procedure PLPerp(Const P : Point; Const L : Line; Var L1 : Line);
Procedure SPerp(Const P1,P2 : Point; Var L : Line);
Procedure P3Circle(Const P1,P2,P3 : Point; Var C : Circle);

Function L2Point(Const L1,L2 : Line; Var P : Point) : Boolean;
Function CPPoints(Const C: Circle; Const P: Point; Var P1,P2: Point) : Byte;
Function CLPoints(Const C: Circle; Const L: Line; Var P1,P2: Point) : Byte;
Function C2Points(Const C1,C2 : Circle; Var P1,P2 : Point) : Byte;

Function PLClass(Const P : Point; Const L : Line) : PClass;
Function PCClass(Const P : Point; Const C : Circle) : PClass;
Function PSClass(Const P,P1,P2 : Point) : PClass;
Function PPolyClass(Const P : Point; N : Word; Const A : Array Of Point) : PClass;

Function P2Dist(Const P1,P2 : Point) : Real;
Function P2SqrDist(Const P1,P2 : Point) : Real;

Function SSquare(Const P1,P2 : Point) : Real;
Function P3Square(Const P1,P2,P3 : Point) : Real;
Function PolySquare(N : Word; Const A : Array Of Point) : Real;

Function PAngle(Const P : Point) : Real;
Function PPAngle(Const P1,P2 : Point) : Real;
Function PSAngle(Const P,P1,P2 : Point) : Real;

Procedure ConvexHull(N : Word; Const A : Array Of Point;
                     Var M : Word; Var Ind : Array Of Word);

Implementation

{ Проверяет два числа на равенство с точностью до Eps }
Function Eq(Const a,b : Real) : Boolean;
Begin
  Eq:=Abs(a-b)<=Eps;
End;

{ Проверяет два числа на неравенство с точностью до Eps }
Function Ne(Const a,b : Real) : Boolean;
Begin
  Ne:=Abs(a-b)>Eps;
End;

{ Проверяет два числа на "меньше или равно" с точностью до Eps }
Function Le(Const a,b : Real) : Boolean;
Begin
  Le:=(a<=b) Or Eq(a,b);
End;

{ Проверяет два числа на "меньше" с точностью до Eps }
Function Lt(Const a,b : Real) : Boolean;
Begin
  Lt:=(a<b) And Ne(a,b);
End;

{ Проверяет два числа на "больше или равно" с точностью до Eps }
Function Ge(Const a,b : Real) : Boolean;
Begin
  Ge:=(a>=b) Or Eq(a,b);
End;

{ Проверяет два числа на "больше" с точностью до Eps }
Function Gt(Const a,b : Real) : Boolean;
Begin
  Gt:=(a>b) And Ne(a,b);
End;

{ Проверяет принадлежность числа заданному диапазону с точностью до Eps }
Function Btw(Const x,a,b : Real) : Boolean;
Begin
  Btw:=Le(a,x) And Le(x,b) Or Ge(a,x) And Ge(x,b);
End;

{ Аналогично предыдущей функции, но проверка строгая }
Function BtwS(Const x,a,b : Real) : Boolean;
Begin
  BtwS:=(a<=x) And (x<=b) Or (a>=x) And (x>=b);
End;

{ Считает определитель 2x2 }
Function Det(Const a,b,c,d : Real) : Real;
Begin
  Det:=a*d-b*c;
End;

{ Вычисляет арктангенс — число из диапазона (-Pi,Pi] }
Function ArcTan2(Const y,x : Real) : Real; Assembler;
Asm
        fld     y
        fld     x
        fpatan
        fwait
End;

{ Находит полярный угол точки с координатами (x,y) }
{ Возвращает число из диапазона [0,2*Pi) }
{ ASSERT: хотя бы одно из чисел x и y отлично от нуля! }
Function PolarAngle(Const x,y : Real) : Real;
Var p : Real;
Begin
  p:=ArcTan2(y,x);
  If p<0 Then p:=p+2*Pi;
  PolarAngle:=p;
End;

{ Находит ориентированный угол от a до b }
{ Возвращает число из диапазона [0,2*Pi) }
Function ADiff(Const a,b : Real) : Real;
Begin
  If b>=a Then ADiff:=b-a Else ADiff:=b-a+2*Pi;
End;

{====================================================================}
{                                                                    }
{====================================================================}

{ Инициализирует объект типа Point }
Procedure PSet(Var P : Point; Const x,y : Real);
Begin
  P.x:=x; P.y:=y;
End;

{ Инициализирует объект типа Line }
Procedure LSet(Var L : Line; Const A,B,C : Real);
Begin
  L.A:=A; L.B:=B; L.C:=C;
End;

{ Инициализирует объект типа Circle }
Procedure CSet(Var C : Circle; Const x,y,R : Real);
Begin
  PSet(C.O,x,y);
  C.R:=R;
End;

{ Возвращает строковое представление точки }
Function PStr(Const P : Point; w,d : Integer) : String;
Var sx,sy : String;
Begin
  Str(P.x:w:d,sx);
  Str(P.y:w:d,sy);
  PStr:='('+sx+','+sy+')';
End;

{ Проверяет принадлежность точки замкнутому прямоугольнику, }
{ образованному двумя заданными точками, с точностью до Eps }
Function PBtw(Const P,P1,P2 : Point) : Boolean;
Begin
  PBtw:=Btw(P.x,P1.x,P2.x) And Btw(P.y,P1.y,P2.y);
End;

{ Аналогично предыдущей функции, но проверка строгая }
Function PBtwS(Const P,P1,P2 : Point) : Boolean;
Begin
  PBtwS:=BtwS(P.x,P1.x,P2.x) And BtwS(P.y,P1.y,P2.y);
End;

{ Проходит от точки P1 заданную часть отрезка [P1,P2] }
Procedure PMoveTo(Const P1,P2 : Point; Const k : Real; Var P : Point);
Begin
  P.x:=P1.x+(P2.x-P1.x)*k;
  P.y:=P1.y+(P2.y-P1.y)*k;
End;

{ Проектирует точку на прямую }
Procedure PLProject(Const P : Point; Const L : Line; Var P1 : Point);
Var k : Real;
Begin
  With P,L Do
  Begin
    k:=(A*x+B*y+C)/(Sqr(A)+Sqr(B));
    P1.x:=x-A*k;
    P1.y:=y-B*k;
  End;
End;

{ Проводит прямую с заданными коэффициентами A и B через точку P }
Procedure PLine(Const A,B : Real; Const P : Point; Var L : Line);
Begin
  L.A:=A; L.B:=B;
  L.C:=-A*P.x-B*P.y;
End;

{ Проводит прямую через две точки }
{ ASSERT: точки P1 и P2 различны! }
Procedure P2Line(Const P1,P2 : Point; Var L : Line);
Begin
  PLine(P1.y-P2.y,P2.x-P1.x,P1,L);
End;

{ Строит прямую, перпендикулярную заданной прямой и проходящую }
{ через указанную точку }
Procedure PLPerp(Const P : Point; Const L : Line; Var L1 : Line);
Begin
  PLine(L.B,-L.A,P,L1);
End;

{ Строит серединный перпендикуляр к отрезку }
{ ASSERT: точки P1 и P2 различны! }
Procedure SPerp(Const P1,P2 : Point; Var L : Line);
Var L1 : Line;
    P : Point;
Begin
  P2Line(P1,P2,L1);
  PMoveTo(P1,P2,1/2,P);
  PLPerp(P,L1,L);
End;

{ Проводит окружность через три точки }
{ ASSERT: точки P1, P2 и P3 различны! }
Procedure P3Circle(Const P1,P2,P3 : Point; Var C : Circle);
Var L1,L2 : Line;
Begin
  SPerp(P1,P2,L1);
  SPerp(P2,P3,L2);
  L2Point(L1,L2,C.O);
  C.R:=P2Dist(C.O,P1);
End;

{====================================================================}
{                            Пересечения                             }
{====================================================================}

{ Пересекает две прямые }
Function L2Point(Const L1,L2 : Line; Var P : Point) : Boolean;
Var D : Real;
Begin
  L2Point:=False;
  D:=Det(L1.A,L2.A,L1.B,L2.B);
  If Eq(D,0) Then P:=PInf
  Else Begin
    P.x:=-Det(L1.C,L2.C,L1.B,L2.B)/D;
    P.y:=-Det(L1.A,L2.A,L1.C,L2.C)/D;
    L2Point:=True;
  End;
End;

{ Пересекает окружность с прямой, проходящей через заданную точку  }
{ и перпендикулярной отрезку, соединяющему ее с центром окружности }
Function CPPoints(Const C : Circle; Const P: Point; Var P1,P2: Point) : Byte;
Var d,t,xx,yy : Real;
Begin
  P1:=PInf; P2:=PInf;
  With C Do
  Begin
    d:=P2Dist(O,P); CPPoints:=0;
    If Eq(d,R) Then CPPoints:=1
    Else If d<R Then CPPoints:=2
    Else Exit;

    t:=Sqr(R)-Sqr(d);
    If t<0 Then t:=0 Else t:=Sqrt(t)/d;
    xx:=(P.x-O.x)*t; yy:=(P.y-O.y)*t;
    P1.x:=P.x-yy; P1.y:=P.y+xx;
    P2.x:=P.x+yy; P2.y:=P.y-xx;
  End;
End;

{ Пересекает окружность и прямую }
Function CLPoints(Const C : Circle; Const L: Line; Var P1,P2: Point) : Byte;
Var P : Point;
Begin
  PLProject(C.O,L,P);
  CLPoints:=CPPoints(C,P,P1,P2);
End;

{ Пересекает две окружности }
Function C2Points(Const C1,C2 : Circle; Var P1,P2 : Point) : Byte;
Var P : Point;
    d2,t : Real;
Begin
  d2:=P2SqrDist(C1.O,C2.O);
  t:=(Sqr(C1.R)-Sqr(C2.R)+d2)/(2*d2);
  PMoveTo(C1.O,C2.O,t,P);
  C2Points:=CPPoints(C1,P,P1,P2);
End;

{====================================================================}
{                        Классификация точек                         }
{====================================================================}

{ Определяет расположение точки относительно прямой }
Function PLClass(Const P : Point; Const L : Line) : PClass;
Var d : Real;
Begin
  With P,L Do d:=A*x+B*y+C;
  If Eq(d,0) Then PLClass:=_On
  Else If d<0 Then PLClass:=_Negative
  Else PLClass:=_Positive;
End;

{ Определяет расположение точки относительно окружности }
Function PCClass(Const P : Point; Const C : Circle) : PClass;
Var d : Real;
Begin
  d:=P2Dist(C.O,P);
  If Eq(d,C.R) Then PCClass:=_On
  Else If d<C.R Then PCClass:=_Inside
  Else PCClass:=_Outside;
End;

{ Определяет расположение точки относительно отрезка }
{ ASSERT: точки P1 и P2 различны! }
Function PSClass(Const P,P1,P2 : Point) : PClass;
Var L: Line;
Begin
  P2Line(P1,P2,L);
  If (PLClass(P,L)=_On) And PBtw(P,P1,P2) Then PSClass:=_On
  Else PSClass:=_Outside;
End;

{ Определяет принадлежность точки многоугольнику, }
{ заданному массивом вершин A, N - число вершин   }
Function PPolyClass(Const P: Point; N: Word; Const A: Array Of Point): PClass;
Var i : Word;
    S : Real;
Begin
  PPolyClass:=_On;
  If PSClass(P,A[N-1],A[0])=_On Then Exit;
  S:=PSAngle(P,A[N-1],A[0]);
  For i:=0 To N-2 Do
  Begin
    If PSClass(P,A[i],A[i+1])=_On Then Exit;
    S:=S+PSAngle(P,A[i],A[i+1]);
  End;
  If Abs(S)<Pi Then PPolyClass:=_Outside { Точка P лежит снаружи }
  Else PPolyClass:=_Inside;              { Точка P лежит внутри }
End;

{====================================================================}
{                             Расстояния                             }
{====================================================================}

{ Считает расстояние между двумя точками }
Function P2Dist(Const P1,P2 : Point) : Real;
Begin
  P2Dist:=Sqrt(Sqr(P1.x-P2.x)+Sqr(P1.y-P2.y));
End;

{ Считает квадрат расстояния между двумя точками }
Function P2SqrDist(Const P1,P2 : Point) : Real;
Begin
  P2SqrDist:=Sqr(P1.x-P2.x)+Sqr(P1.y-P2.y);
End;

{====================================================================}
{                        Периметры и площади                         }
{====================================================================}

{ Вычисляет ориентированную площадь треугольника,  }
{ натянутого на отрезок [P1,P2] и начало координат }
Function SSquare(Const P1,P2 : Point) : Real;
Begin
  SSquare:=Det(P1.x,P2.x,P1.y,P2.y)/2;
End;

{ Вычисляет ориентированную площадь треугольника, }
{ натянутого на точки P1, P2 и P3 }
Function P3Square(Const P1,P2,P3 : Point) : Real;
Begin
  P3Square:=Det(P1.x-P3.x,P2.x-P3.x,P1.y-P3.y,P2.y-P3.y)/2;
End;

{ Вычисляет ориентированную площадь многоугольника, }
{ заданного массивом вершин A, N - число вершин }
Function PolySquare(N : Word; Const A : Array Of Point) : Real;
Var S : Real;
    i : Word;
Begin
  S:=SSquare(A[N-1],A[0]);
  For i:=0 To N-2 Do S:=S+SSquare(A[i],A[i+1]);
  PolySquare:=S;
End;

{====================================================================}
{                                Углы                                }
{====================================================================}

{ Определяет полярный угол точки }
Function PAngle(Const P : Point) : Real;
Begin
  PAngle:=PolarAngle(P.x,P.y);
End;

{ Определяет полярный угол точки P2 относительно точки P1 }
Function PPAngle(Const P1,P2 : Point) : Real;
Begin
  PPAngle:=PolarAngle(P2.x-P1.x,P2.y-P1.y);
End;

{ Определяет ориентированный угол, под которым виден отрезок }
{ из заданной точки }
Function PSAngle(Const P,P1,P2 : Point) : Real;
Var a : Real;
Begin
  a:=ADiff(PPAngle(P,P1),PPAngle(P,P2));
  If a<=Pi Then PSAngle:=a Else PSAngle:=a-2*Pi;
End;

{====================================================================}
{                         Выпуклая оболочка                          }
{====================================================================}

{ Строит выпуклую оболочку множества из N точек, заданного массивом A, }
{ занося в массив Ind номера точек из A, являющихся вершинами выпуклой }
{ оболочки, в порядке обхода против часовой стрелки, а в M - их количество }
Procedure ConvexHull(N : Word; Const A : Array Of Point;
                     Var M : Word; Var Ind : Array Of Word);
Var First,Cur,Min,i : Word;
    z : Real;
Begin
  { Метод Джарвиса }
  { Находим самую нижнюю из самых левых точек с минимальным номером }
  First:=0;
  For i:=1 To N-1 Do
    If (A[i].x<A[First].x) Or (A[i].x=A[First].x) And (A[i].y<A[First].y)
       Then First:=i;

  { "Заворачиваем подарок" }
  M:=0; Cur:=First;
  Repeat
    Ind[M]:=Cur; Inc(M); Min:=Cur;
    { Цикл в противоположном направлении и строгие неравенства! }
    For i:=N-1 DownTo 0 Do
    Begin
      If (A[i].x=A[Cur].x) And (A[i].y=A[Cur].y) Then Continue;
      z:=P3Square(A[Cur],A[i],A[Min]);
      If (Min=Cur) Or (z>0) Or (z=0) And PBtwS(A[Min],A[Cur],A[i]) Then Min:=i;
    End;
    Cur:=Min;
  Until Cur=First;
End;

End.