====== Геометрия на плоскости ======

===== Определение принадлежности точки к многоугольнику =====

=== Алгоритм: ===

Точка лежит внутри многоугольника, если кол-во точек пересечения нечетно.
Все   отрезки кроме горизонтальных проверяются на пересечение с горизонтальным лучом, выходящим из проверяемой точки.\\ При попадании луча в вершину пересечение засчитывается только с теми отрезками, выходящими из вершины, для которых она является верхней.

== Примечание: ==
Порядок ввода точек должен быть **против** часовой стрелки. На границе значение функции не определено.

<code pascal>
{
   Code by VanDamM // [WRC]
}

Program Check_Point_In_Poly;
Uses Crt;

Type Point = Record   { тип точка }
      x, y : integer;
     End;

Var  PointXY : Point;                 { проверяемая точка }
     Poly    : array[0..24] of Point; { массив вершин многоугольника }
     C       : integer;               { кол-во вершин многоугольника }
     i, j    : integer;

Function Max( Num1, Num2 : integer ) : integer;
Begin
 If Num1>Num2 then Max:=Num1 else Max:=Num2;
End;

Function Min( Num1, Num2 : integer ) : integer;
Begin
 If Num1<Num2 then Min:=Num1 else Min:=Num2;
End;

Procedure EnterData; { Процедура ввода данных }
Begin
 Write('Enter poly`s vertex number: '); ReadLn ( C );
 For i:=0 to C-1 do
  begin
   Write('X[',i,']: '); ReadLn(Poly[i].x);
   Write('Y[',i,']: '); ReadLn(Poly[i].y);
  end;
   WriteLn;
   Write('Point X: '); ReadLn(PointXY.x);
   Write('Point Y: '); ReadLn(PointXY.y);
End;

Function PointInPoly ( A : Point; P : array of Point; N : integer) : integer;
Var Count : integer;
    T     : real;
Begin
  T:=0;
  Count:=0;
 For i:=0 to N-1 do
  begin
   j:=(i+1) mod N;
   If P[i].y = P[j].y then Continue;
   If (P[i].y > A.y) and (p[j].y > A.y) then Continue;
   If (P[i].y < A.y) and (p[j].y < A.y) then Continue;
   If Max(P[i].y, P[j].y) = A.y then
    Inc(Count)
   else
    If Min(P[i].y, P[j].y) = A.y then
     Continue
    else
     begin
      T := (A.y-P[i].y)/(P[j].y-P[i].y);
      If ((T>0) and (T<1)) and ((P[i].x + T*(P[j].x-P[i].x)) >= A.x) then
       Inc(Count);
     end;
  end;
 PointInPoly:= Count AND 1;
End;

Begin
 ClrScr;

 EnterData;

 writeln;

 If PointInPoly(PointXY, Poly, C) = 0 then
  Write('Answer: Point out of poly')
 else
  Write('Answer: Point in poly');

 readln;
End.
</code>


===== Построение выпуклой оболочки множества =====

//Выпуклая оболочка множества точек// - минимальный многоугольник, который содежит все точки данного множества.

=== Алгоритм решения задачи: ===
скажем, что начальная вершина V = V0 - самая нижняя (если таких несколько, то самая правая). Выберем произвольную ось (допустим, параллельную оси ОХ и идущую через первую вершину). Для каждой вершины T найдем углы между осью ОХ и отрезком VT. Выберем наименьший такой угол (углы считаются против часовой стрекли с "3 часов"). Выведем вершину T. Теперь присвоим V = T. Когда T = V0 - алгоритм заканчивается.

<code pascal>
{ Построение выпуклой оболочки множества }
type
  TPoint = record
    x,y: Integer;
  end;

const
  MaxN = 400;

var
  N: Integer;
  i: Integer;
  P: array [1..MaxN] of TPoint;     { Точки множества }
  Used: array [1..MaxN] of Boolean;

  First,Point,NewPoint: Integer;

  procedure WritePoint(i: Integer);
  begin
    Writeln('[',i:3,'] (',P[i].x, ';', P[i].y, ')');
  end;

  function Bigger(P,P1,P2: TPoint): Boolean;
  begin
    Bigger := (P1.x - P.x) * (P2.y - P.y) - (P1.y - P.y) * (P2.x - P.x) > 0;
  end;

begin
  assign(input,'input.txt'); reset(input);
  Readln(N);

  First:=1;
  for i:=1 to N do begin
    Readln(P[i].x, P[i].y);
    if (P[i].y < P[First].y) then First:=i else
    if (P[i].y = P[First].y) and (P[i].x < P[First].x) then First:=i;
  end;
  close(input);

  FillChar(Used,SizeOf(Used),False);

  Writeln('Выпуклая оболочка идет через вершины: ');

  Point:=First;

  repeat
    WritePoint(Point); Used[Point]:=True;
    NewPoint:=0;

    { Поиск следующей вершины }
    for i:=1 to N do
      if i <> Point then begin
        if (NewPoint = 0) or (Bigger(P[Point], P[i], P[NewPoint]))
          then NewPoint:=i;
      end;

    Point:=NewPoint;
  until Point = First;

end.</code>


===== Определить, что точка находится внутри треугольника =====

Определить, что точка находиться внутри треугольника можно намного проще.\\ Пусть даны координаты вершин треуголькика, т.е. координаты вершин A, B и C. Надо определить, находится ли точка D внутри или снаружи. Если S(ABC) = S(ABD) + S(ACD) + S(BCD), то точка лежит внутри (это видно чисто геометрически). В противном случае точка лежит вне треуголькика. 

== Вот реализация функции: ==
<code pascal>Const
   Epsilon = 1e-6;
Function InTriangle(X1,Y1,X2,Y2,X3,Y3,X,Y:Real):Boolean;
{    /
    / \
   b   c
 /      \
------a---
}
Function Dist(XX1,YY1,XX2,YY2:Real):Real;
Begin
 Dist:=Sqrt (Sqr(XX2 - XX1) + Sqr(YY2 - YY1));
End;

Function Square(A,B,C:Real):Real;
Var P:Real;
Begin
 P := (A + B + C) / 2;
 Square := Sqrt (P * (P - A) * (P - B) * (P - C));
End;

Var
   Ab,Ac,Bc:Real;
   S,S1:Real;
   Da,Db,Dc:Real;
Begin
 Ab := Dist(X1, Y1, X2, Y2);
 Bc := Dist(X2, Y2, X3, Y3);
 Ac := Dist(X1, Y1, X3, Y3);
 S := Square(Ab, Bc, Ac);
 Da := Dist(X1, Y1, X, Y);
 Db := Dist(X2, Y2, X, Y);
 Dc := Dist(X3, Y3, X, Y);
 S1 := Square(Ac, Da, Dc) + Square(Ab, Db, Da) + Square(Bc, Db, Dc);
 InTriangle := Abs(S - S1) < Epsilon;
End.</code>

===== Найти радиус вписанной окружности =====

Дан треугольник со сторонами A, B, C. Требуется найти радиус вписанной окружности. Вывести рисунок на экран в графическом режиме.
<code pascal>
uses crt,graph;

{Инициализация графики}
Procedure InitGr;
var
  grDriver: Integer;
  grMode: Integer;
begin
  grDriver := Detect;
  InitGraph(grDriver, grMode,'..\BGI'); { путь можно указывать относительно или абсолютно }
End;

{Вычисление площади треугольника}
Function Square(A, B, C : Extended) : Extended;
Var P : Extended;
Begin
 P:=(A + B + C) / 2;
 Square:=Sqrt(P * (P - A) * (P - B) * (P - C));
End;

{Вычисление радиуса треугольника}
Function Radius(A, B, C : Extended) : Extended;
{S = p*r => r = S/p}
Begin
 Radius:=2 * Square(A, B, C) / (A + B + C);
End;

{Рисование всей картинки: A, B, C - стороны треугольника, X, Y - координаты на экране, Scale - масштаб}
Procedure Draw(A, B, C : Extended; X, Y : Longint; Scale : Extended);
Var R : Extended;
    X1, Y1, X2, Y2, X3, Y3 : Extended;
    SinAlpha, CosAlpha : Extended;
    TgAlphaDiv2, SinAlphaDiv2, CosAlphaDiv2 : Extended;
    Xr, Yr : Extended;
Begin
 A := A * Scale; B := B * Scale; C := C * Scale;
 R := Radius(A, B, C);
 {Пусть одна из сторон будет нарисована горизонтальной}
 X1 := 0; Y1 := 0;
 X2 := C; Y2 := 0;
 {Найдем координаты третей вершины. Сначала найдем синус и косинус одного из углов треугольника}
 CosAlpha := (Sqr(B) + Sqr( C ) - Sqr(A)) / (2 * B * C);
 SinAlpha := Sqrt(1 - Sqr(CosAlpha));
 X3 := B * CosAlpha;
 Y3 := B * SinAlpha;
 {Чертим треугольник}
 Line(X + Round(X1), Y + Round(Y1), X + Round(X2), Y + Round(Y2));
 Line(X + Round(X1), Y + Round(Y1), X + Round(X3), Y + Round(Y3));
 Line(X + Round(X3), Y + Round(Y3), X + Round(X2), Y + Round(Y2));
 {Центр вписанной окружности лежит в точке пересечения биссектрис}
 CosAlphaDiv2 := Sqrt((1 + CosAlpha) / 2);
 SinAlphaDiv2 := Sqrt(1 - Sqr(CosAlphaDiv2));
 TgAlphaDiv2 := SinAlphaDiv2 / CosAlphaDiv2;
 {Найдем координаты центра окружности}
 Xr := R / TgAlphaDiv2;
 Yr := R;
 {Рисуем окружность}
 Circle(X + Round(Xr), Y + Round(Yr), Round(R));
End;

Begin
  InitGr;

  Draw(6, 5, 6, 100, 200, 30.0);

  repeat until keypressed;
  CloseGraph;
End.</code>


===== Нахождение центра описанной окружности =====

=== Алгоритм: ===
Из школьной программы известно, что через три точки, не лежащие на одной прямой, проходит окружность, притом только одна. Уравнение окружности с центром (x<sub>0</sub>, y<sub>0</sub>) и радиусом r выглядит так:
(x - x0)<sup>2</sup> + (y - y0)<sup>2</sup> = r<sup>2</sup>.
Пусть известно, что окружность идет через точки с координатами (x1, y1), (x2, y2) и (x3, y3). Это обозначает, что будут выполняться условия:
(x1 - x0)^2 + (y1 - y0)^2 = r^2
(x2 - x0)^2 + (y2 - y0)^2 = r^2
(x3 - x0)^2 + (y3 - y0)^2 = r^2
Получаем три уравнения с тремя неизвестными. Осталось аккуратно открыть скобки, произвести нехитрые действия (на самом деле, система сводится к трем линейным уравнениям, которые решить можно как угодно, хотя бы методом Крамера, что и сделано в программе).
<code pascal>
Type
     TPoint = Record {Тип точка}
      X, Y : Extended;
     End;

{Поиск центра описанной окружности. Если точки лежат на одной прямой, функция возвращает false, иначе true}

Function FindOuterRadius(A, B, C : TPoint; Var Rr : TPoint): Boolean;
Var M : Array[1..2,1..3] Of Extended;
    D, Dx, Dy : Extended;
Begin
 M[1, 1] := 2 * (A.X - B.X);
 M[1, 2] := 2 * (A.Y - B.Y);
 M[1, 3] := Sqr(A.X) +Sqr(A.Y) - (Sqr(B.X) + Sqr(B.Y));

 M[2, 1] := 2 * (B.X - C.X);
 M[2, 2] := 2 * (B.Y - C.Y);
 M[2, 3] := Sqr(B.X) +Sqr(B.Y) - (Sqr(C.X) + Sqr(C.Y));

 D := M[1, 1] * M[2, 2] - M[2, 1] * M[1, 2];
 Dx := M[1, 3] * M[2, 2] - M[2, 3] * M[1, 2];
 Dy := M[1, 1] * M[2, 3] - M[2, 1] * M[1, 3];

 If D <> 0 Then
  Begin
   Rr.X := Dx/D;
   Rr.Y := Dy/D;
   FindOuterRadius := True;
  End Else
  Begin
   Rr.X := 0;
   Rr.Y := 0;
   FindOuterRadius := False;
  End;
End;</code> 


===== Проведение параболы через три точки =====

Через три точки, не лежащие на одной прямой проходит парабола, причем только одна (если три точки лежат на одной прямой, то парабола вырождается в прямую линию).

=== Алгоритм: ===
Даны три точки A(x1, y1), B(x2, y2) и C(x3, y3), через которые надо провести параболу.
Уравнение параболы y = ax<sup>2</sup> + bx + c.
Составим систему с неизвестными a, b и c:
<code>y1 = a*x1^2 + b*x1 + c
y2 = a*x2^2 + b*x2 + c
y3 = a*x3^2 + b*x3 + c
</code>
которая эквивалентна данной:
<code>
a(x1^2 - x2^2) + b(x1 - x2) = y1 - y2
a(x2^2 - x3^2) + b(x2 - x3) = y2 - y3
c = y1 - a*x1^2 - b*x1
</code>
Первые два уравения решаются методом Крамера через определители второго порядка. Коэфициент c находится из третьего уравнения.
<code pascal>
{$N+,E+}
Type TPoint = Record
      X, Y : Extended;
     End;

Procedure Parabola(A, B, C : TPoint; Var KoefA, KoefB, KoefC : Extended);
Var M : Array[1..2, 1..3] Of Extended;
    D, Da, Db : Extended;
Begin
 M[1, 1] := Sqr(A.X) - Sqr(B.X);
 M[1, 2] := A.X - B.X;
 M[1, 3] := A.Y - B.Y;

 M[2, 1] := Sqr(B.X) - Sqr(C.X);
 M[2, 2] := B.X - C.X;
 M[2, 3] := B.Y - C.Y;

 D := M[1, 1] * M[2, 2] - M[1, 2] * M[2, 1];

 If D <> 0 Then
  Begin
   Da := M[1, 3] * M[2, 2] - M[2, 3] * M[1, 2];
   Db := M[1, 1] * M[2, 3] - M[1, 2] * M[1, 3];
   KoefA := Da / D;
   KoefB := Db / D;
   KoefC := A.Y - KoefA * Sqr(A.X) - KoefB * A.X;
  End
  Else
  Begin
   KoefA := 0; KoefB := 0; KoefC := 0;
  End;
End;

Var A, B, C : TPoint;
    KoefA, KoefB, KoefC : Extended;

Begin
 A.X := -1; A.Y := 1;
 B.X :=  0; B.Y := 0;
 C.X :=  1; C.Y := 2;
 Parabola(A, B, C, KoefA, KoefB, KoefC);
 Writeln('y = ', KoefA:10:10, 'x^2 + ', KoefB:10:10, 'x + ', KoefC:10:10);
End.
</code>