Реализация
{ Модуль для работы с геометрическими объектами }
{ 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.