Геометрия на плоскости , Алгоритмы и задачи

[Romtek]

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

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

---

Примечание: порядок ввода точек должен быть против часовой стрелки.

{
   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.

Сообщение отредактировано: Romtek - 02.03.05, 23:35

[tserega]

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

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

{ Построение выпуклой оболочки множества }
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.

Сообщение отредактировано: romtek - 23.05.04, 13:27

[tserega]

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

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

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.

Сообщение отредактировано: romtek - 23.05.04, 13:33

[tserega]

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

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

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.

Сообщение отредактировано: romtek - 23.05.04, 13:41

[tserega]

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

Из школьной программы известно, что через три точки, не лежащие на одной прямой, проходит окружность, притом только одна. Уравнение окружности с центром (x0, y0) и радиусом r выглядит так:
(x - x0)^2 + (y - y0)^2 = r^2.
Пусть известно, что окружность идет через точки с координатами (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
Получаем три уравнения с тремя неизвестными. Осталось аккуратно открыть скобки, произвести нехитрые действия (на самом деле, система сводиться к трем линейным уравнениям, которые решить можно как угодно, хотя бы методом Крамера, что и сделано в программе).

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;

Добавлено в 08.05.04, 14:34:

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

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

Даны три точки A(x1, y1), B(x2, y2) и C(x3, y3), через которые надо провести параболу.
Уравнение параболы y = ax^2 + bx + c.
Составим систему с неизвестными a, b и c:
y1 = a*x1^2 + b*x1 + c
y2 = a*x2^2 + b*x2 + c
y3 = a*x3^2 + b*x3 + c
которая эквивалентна данной:
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
Первые два уравения решаются методом Крамера через определители второго порядка. Коэфициент c находится из третьего уравнения.

{$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.

Сообщения были разделены в тему "Построение выпуклой оболочки"

Сообщение отредактировано: romtek - 23.05.04, 13:49