ru.algorithms

 
 - RU.ALGORITHMS ----------------------------------------------------------------
 From : Michael Sedov                        2:5015/185.2   04 Apr 2002  18:32:37
 To : All
 Subject : [2] Faq по Геометрии.
 -------------------------------------------------------------------------------- 
 

   +  Константы и типы.
   1. Площадь многоугольника.
   2. Проверка вхождения точки в треугольник.
   3. Точка пересечения перпендикуляра и прямой.
   4. Задача ь 3 с использованием параметрического задания прямой.
   5. Расстояние от точки до прямой.
   6. Принадлежность точки отрезку.
   7. Расстояние от точки до прямой (прямая задана через Ax + By + C = 0).
   8. Принадлежность точки прямой.
   9. Получение уравнения прямой.
  10. Расстояние от точки до отрезка. 
  11. Проверка выпуклости многоугольника.
  12. Проверка вхождения точки в многоугольник.
 
 const
      maxn = 100;  { Максимальное количество вершин в многоугольнике }
      eps  = 1e-9; { Малое число выбирается в зависимости от необходимой 
                     точности ответа, подробнее в 1.3 }
 
 { Используемые типы }
 type
     float = double; { Для простоты замены используемого типа в зависимости 
                       от необходимой точности ограничения на память. }
     tpoint    = record x, y : float; end;
     tpiece    = record { Отрезок, задается кому как нравится }
       case byte of
       0: (a: array [1..2] of tpoint;);
       1: (p1, p2: tpoint);
       2: (x1, y1, x2, y2: float);
     end;
     tpoly = array[0..maxn] of tpoint;
       { Первый способ задания полигона, количество вершин
         храниться в отдельной переменной, лучше использовать
         когда в задаче рассматривается 1-2 полигона }
     tpolygon = record
       n : integer;
       p : tpoly;
       { Второй способ задания полигона, используется когда
         полигонов много}
     end;
     ttriangle = array[1..3] of tpoint;
     tline = record a,b,c: float; end; { Для Ax + By + C = 0 }
 
 { *Внимание*:
  Все функции написаны в расчете на то что выполнено "замыкание"
 многоугольника, 
  т.е. p[0] := p[n]; p[n+1] := p[1]; }
 
 { _Функция возвращает площадь многоугольника._
   Используется свойство векторного произведения : векторное произведение равно
  ориентированной прощади  параллелограмма, построенного  на векторах, входящих
  в векторное прозведение. Складывать площадь будем из треугольников (одна вер-
  шина - (0, 0), другие - две соседние из массива tpoly). }
 
 function square(p: tpoly; n: integer): float;
 var
   i: integer;
   s: float;
 begin
   s:= 0.0;
   for i:= 1 to n do
     s:= s + (p[i].x*p[i-1].y-p[i-1].x*p[i].y);
   square:= abs(s)*0.5;
 end;
 
 { _Проверка вхождения точки в треугольник._
   Опять же векторное произведение. Проверяем лежит ли точка по одну сторону от
   всех сторон треугольника.}
 
 function intriangle(q: tpoint; t: ttriangle): boolean;
 var v: float;
 begin
   v:= (q.x-t[1].x)*(t[2].y-t[1].y)-(q.y-t[1].y)*(t[2].x-t[1].x);
   intriangle:= false;
   if((q.x-t[2].x)*(t[3].y-t[2].y)-(q.y-t[2].y)*(t[3].x-t[2].x))*v>-eps then
     if((q.x-t[3].x)*(t[1].y-t[3].y)-(q.y-t[3].y)*(t[1].x-t[3].x))*v>-eps then
     intriangle:= true;
 end;
 
 { _Точка пересечения перпендикуляра опущенного из точки, и прямой_
   _к которой он опущен(прямая задана точками через которые проходит)._
   Голая математика, с использованием свойств в.п. и с.п.}
 
 procedure crsort(q: tpoint; p: tpiece; var w: tpoint);
 var a,b,c,d,e,f,g: float;
 begin
   a:= p.x2-p.x1; b:= p.y2-p.y1; c:= a*q.x+b*q.y;
   d:= b;         e:= -a;        f:= d*p.x1+e*p.y1;
   g:= a*e-b*d;
   w.x:= (c*e-b*f)/g;
   w.y:= (a*f-c*d)/g;
 end;
 
 { Тоже самое только с использованием _параметрического_ задания прямой (может 
   понадобиться для определения попала точка на отрезок, или нет). }
 
 procedure crsortp(q: tpoint; p: tpiece; var t: float; var w: tpoint);
 begin
   t:= -((p.x1-q.x)*(p.x2-p.x1)+(p.y1-q.y)*(p.y2-p.y1))/
         (sqr(p.x2-p.x1)+sqr(p.y2-p.y1));
   w.x:= p.x1*(1.0-t)+p.x2*t;
   w.y:= p.y1*(1.0-t)+p.y2*t;
 end;
 
 { _Расстояние от точки до прямой (прямая задана точками через отрезок, оба 
   конца которого лежат на данной прямой)._
   Используем свойства в.п.}
 
 function point2piece(q: tpoint; p: tpiece): float;
 begin
   point2piece:= abs((q.x-p.x1)*(p.y2-p.y1)-(q.y-p.y1)*(p.x2-p.x1))/
                     sqrt(sqr(p.x2-p.x1)+sqr(p.y2-p.y1));
 end;
 
 { _Принадлежность точки отрезку._
   Считаем расстояние от данной точки до отрезка, и сравниваем с eps.
   1) Квадрат расстояния меньше sqr(eps)
   2) Hаходится между перпендикулярными прямыми, проведенными через концы 
      отрезка.
   Проверяется через с.п.}
 
 function belong2piece(q: tpoint; p: tpiece): boolean;
 var t: float;
 begin
   t:= sqr((q.x-p.x1)*(p.y2-p.y1)-(q.y-p.y1)*(p.x2-p.x1))/
            (sqr(p.x2-p.x1)+sqr(p.y2-p.y1));
   belong2piece:= false;
 
   if t>sqr(eps) then exit;
 
   if((q.x-p.x1)*(p.x2-p.x1)+(q.y-p.y1)*(p.y2-p.y1))*
     ((q.x-p.x2)*(p.x2-p.x1)+(q.y-p.y2)*(p.y2-p.y1))<eps then
       belong2piece:= true;
 end;
 
 { _Расстояние от точки до прямой (прямая задана через Ax + By + C = 0)._
   "Голая теория". Без комментариев.}
 
 function point2line(q: tpoint; s: tline): float;
 begin
   point2line:= abs(s.a*q.x+s.b*q.y+s.c)/sqrt(sqr(s.a)+sqr(s.b));
 end;
 
 { _Принадлежность прямой._
   Считаем квадрат расстояния от точки до прямой и сравниваем с sqr(eps).
   Попутно используем умножение вместо деления. }
 
 function belong2line(q: tpoint; s: tline): boolean;
 begin
   belong2line:= false;
   if sqr(s.a*q.x+s.b*q.y+s.c) < sqr(eps)*(sqr(s.a)+sqr(s.b)) then
     belong2line:= true;
 end;
 
 { _Получение уравнения прямой_ через координаты отрезка, концы которого лежат 
  на этой прямой.
     Задача решена теоретически.  При решение использовался такой способ: запи-
  сывалось параметрическое уравнение прямой и, путём исключения параметра, была
  получена формула. Если часто используется поиск расстояния от точки до
 прямой, 
  или прямую надо двигать на паралленое расстояние, то можно это уравнение нор-
  мировать. Т. е. разделить на sqrt(A^2+B^2). Тогда при нахождении  расстояния 
  не надо делить на эту величину.}
 
 procedure piece2line(p: tpiece; var s: tline);
 begin
   if (abs(p.x1-p.x2) < eps) and (abs(p.y1-p.y2) < eps) then
     writeln('За последствия не ручаюсь!')
   else begin
     s.a:= p.y2-p.y1;
     s.b:= p.x1-p.x2;
     s.c:= -s.a*p.y1-p.x1*s.b;
   end;
 end;
 
 { _Расстояние от точки до отрезка._
     Задача особенна тем, что  рассояние до отрезка не всегда равно расстоянию
  до прямой, содержашей отрезок. Hадо просто немного переделать уже написанный 
  belong2piece}
 
 function distance2piece(q: tpoint; p: tpiece): float;
 var t,w: float;
 begin
   if((q.x-p.x1)*(p.x2-p.x1)+(q.y-p.y1)*(p.y2-p.y1))*
     ((q.x-p.x2)*(p.x2-p.x1)+(q.y-p.y2)*(p.y2-p.y1))>-eps then
   begin
     t:= sqr(q.x-p.x1)+sqr(q.y-p.y1);
     w:= sqr(q.x-p.x2)+sqr(q.y-p.y2);
     if w<t then t:= w;
   end else
     t:= sqr((q.x-p.x1)*(p.y2-p.y1)-(q.y-p.y1)*(p.x2-p.x1))/
              (sqr(p.x2-p.x1)+sqr(p.y2-p.y1));
   distance2piece:= sqrt(t);
 end;
 
 { _Проверка выпуклости многоугольника._
   "Классическая задача". Попарно берётся два соседних вектора и векторное
   произведение должно иметь один и тот же знак для любой вершины.}
 
 function isconvex(p: tpoly; n: integer): boolean;
 var
   t: float;
   i: integer;
 begin
   isconvex:= true;
   t:= (p[1].x-p[n].x)*(p[2].y-p[1].y)-(p[2].x-p[1].x)*(p[1].y-p[n].y);
   for i:= 1 to n-1 do
   begin
     if abs(t)<eps then
       t:= (p[i+1].x-p[i].x)*(p[i+2].y-p[i+1].y)-
           (p[i+2].x-p[i+1].x)*(p[i+2].y-p[i].y) else
     if t*((p[i+1].x-p[i].x)*(p[i+2].y-p[i+1].y)-
           (p[i+2].x-p[i+1].x)*(p[i+2].y-p[i].y))<-eps then
       begin isconvex:= false; break end;
   end;
 end;
 
 { _Проверка вхождения точки в многоугольник._
   Считается количество пересечений луча из этой точки и многоугольника.}
 
 function inpoly(p: tpoint; q: tpoly; n: integer): shortint;
 var ans: boolean;
   function crs(a1,a2: tpoint): integer;
   var x: float;
   begin
     crs:= 0;
     if abs(a1.y-a2.y) < eps then
     begin
       if(abs(p.y-a1.y) < eps) and ((p.x-a1.x)*(p.x-a2.x)<0.0)then ans:= false;
       exit;
     end;
     if((a1.y-p.y)*(a2.y-p.y)>0.0)then exit;
     x:= a2.x-(a2.y-p.y)/(a2.y-a1.y)*(a2.x-a1.x);
     if abs(x-p.x)<eps then ans:= false else
     if(x<p.x)then
     begin
       crs:= 1;
       if(abs(a1.y-p.y)<eps)and(a1.y<a2.y)then crs:= 0 else
       if(abs(a2.y-p.y)<eps)and(a2.y<a1.y)then crs:= 0;
     end;
   end;
 var i,c: integer;
 begin
   c:= 0; ans:= true;
   for i:= 1 to n do
   begin
     inc(c,crs(q[i],q[i+1]));
     if not ans then break;
   end;
   if not ans then inpoly:= -1 else inpoly:= c and 1;
 end;
 --- WP/95 Rel 1.78E (215.0) Reg.
  * Origin: [2] Faq по Геометрии. (2:5015/185.2)
 
 

Вернуться к списку тем, сортированных по: возрастание даты  уменьшение даты  тема  автор