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

Дано: у нас есть треугольник, нам известны только координаты его вершин. У нас есть точка, нам известны её координаты.

Что нужно узнать: нужно установить принадлежность точки треугольнику.

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

Метод сравнения площадей

Метод площадей
Рис. 1.

В данном методе сначала находятся площади 3-х треугольников, которые образует данная точка с каждой стороной треугольника. В нашем случае(рис. 1) это треугольники ABP, BCP, CAP и их площади s1, s2, s3 соответственно.

Затем находится площадь самого треугольника ABC.

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

Так как у нас известны только координаты точек, то все площади, находятся по формуле Герона, от обильности операций которой становится ясно, почему этот метод очень трудоёмкий.

Простейшая реализация алгоритма:

// проверка принадлежности точки треугольнику формулами Герона
function IsPointIn_Geron(aAx, aAy, aBx, aBy, aCx, aCy, aPx, aPy: single): boolean;
// funcs
  // площадь треугольника по точкам
  function Square(aAx, aAy, aBx, aBy, aCx, aCy: single): single;
  begin
    Result := abs(aBx*aCy - aCx*aBy - aAx*aCy + aCx*aAy + aAx*aBy - aBx*aAy);
  end;
var
  s : single;
begin
  s := Square(aPx, aPy, aAx, aAy, aBx, aBy) + Square(aPx, aPy, aBx, aBy, aCx, aCy) +
  Square(aPx, aPy, aCx, aCy, aAx, aAy);
  Result := abs(Square(aAx, aAy, aBx, aBy, aCx, aCy) - s) <= 0.01{погрешность};
end;

Атрибуты функции: aAx, aAy, aBx, aBy, aCx, aCy — координаты точек A, B, C треугольника; aPx, aPy — координаты точки, принадлежность которой надо определить.

Метод относительности

Метод относительности
Рис. 2.

Данный метод заключается в следующем. Сначала выбирается ориентация движения по вершинам треугольника(по часовой или против часовой стрелке). Я выбираю по часовой. На рисунке 2 выбранная ориентация движения(по часовой) показана стрелками. По данной ориентации проходим все стороны треугольника, рассматривая их как прямые, и рассчитываем по какую сторону от текущей прямой лежит наша точка. Не трудно догадаться, что если точка для всех прямых, при нашей ориентации, лежит с правой стороны, то значит точка принадлежит треугольнику, а если хоть для какой-то прямой она лежит с левой стороны, то значит условие принадлежности не выполняется.

На рисунке 2 продемонстрирована ситуация, когда точка только для одной прямой AB лежит по левую сторону, а значит не принадлежит треугольнику.

Реализация алгоритма:

// проверка прин. точки треуг. методом относительности положения
function IsPointIn_Relat(aAx, aAy, aBx, aBy, aCx, aCy, aPx, aPy: single): boolean;
// funcs
  function Q(ax, ay, bx, by, atx, aty: single): single;
  begin
    Result := atx * (by - ay) + aty * (ax - bx) + ay * bx - ax * by;
  end;
var
  q1, q2, q3 : single;
begin
  // выбираем определённую ориентацию по вершинам(чтоб было по порядку)
  
  // универсальный
  q1 := Q(aAx, aAy, aBx, aBy, aPx, aPy);
  q2 := Q(aBx, aBy, aCx, aCy, aPx, aPy);
  q3 := Q(aCx, aCy, aAx, aAy, aPx, aPy);
  Result := ((q1 >= 0) and (q2 >= 0) and (q3 >= 0)) or
    ((q1 < 0) and (q2 < 0) and (q3 < 0));
  //}
  {
  // для строгой ориентации по часовой
  Result := (Q(ftx1, fty1, ftx2, fty2, fpx, fpy) >= 0) and
    (Q(ftx2, fty2, ftx3, fty3, fpx, fpy) >= 0) and
    (Q(ftx3, fty3, ftx1, fty1, fpx, fpy) >= 0);
  //}
  {
  // для строгой ориентации против часовой
  Result := (Q(ftx1, fty1, ftx2, fty2, fpx, fpy) <= 0) and
    (Q(ftx2, fty2, ftx3, fty3, fpx, fpy) <= 0) and
    (Q(ftx3, fty3, ftx1, fty1, fpx, fpy) <= 0);
  //}
end;

Всё относительно!

Всё относительно
Рис. 3.

Тут надо кое что пояснить, весьма не маловажное, что может сыграть роль в оптимизации и выборе алгоритма. Обратите внимание, что в приведённом коде есть закомментированные блоки кода с комментариями «для строгой ориентации», в то время как рабочий код универсален — он предназначен для любой ориентации. Т.е. представленный код определит принадлежность точки для любого заданного треугольника. В моей тестирующей программе треугольники как раз таки строятся по random()-у координат вершин, а ориентация идёт по вершинам(A>B>C>A). Для рисунка 2 — это по часовой стрелки, но для рисунка 3 — это против часовой.

Так вот, в случае рисунка 3 точка должна лежать по левую сторону векторов, чтобы принадлежать треугольнику.

Вот тут и получается важный момент! Если вы уверены, что в вашем проекте все треугольники будут ориентированы по часовой стрелке(а т.е. вершина C будет всегда правее вектора AB), то вам можно закомментировать блок универсального решения и раскомментировать блок «для строгой ориентации по часовой» и данный алгоритм упрощается аж на 3 логических операции!

Векторный метод

Третий метод который я освещаю для меня самый интересный.

Идея его применения зарождается если взглянуть на треугольник как на половинку параллелограмма…

Данный метод я сначала проверил на бумаге. После всех оптимизаций формул, как всё сошлось, я реализовал его в коде, где он показал себя вполне успешным и результативным. Аж эффективнее 2-х предыдущих методов :]

Алгоритм такой:

1) одну вершину треугольника помещаем в координаты (0;0);

2) две стороны, выходящие из этой вершины, представляем как вектора.

Таким образом из всего этого появляется система простых условий нахождения точки P между векторами b и c.(рис. 4)

Векторный метод
Рис. 4.

Смотрим код:

function IsPIn_Vector(aAx, aAy, aBx, aBy, aCx, aCy, aPx, aPy: single): boolean;
var
  Bx, By, Cx, Cy, Px, Py : single;
  m, l : single; // мю и лямбда
begin
  Result := False;
  // переносим треугольник точкой А в (0;0).
  Bx := aBx - aAx; By := aBy - aAy;
  Cx := aCx - aAx; Cy := aCy - aAy;
  Px := aPx - aAx; Py := aPy - aAy;
  //
  m := (Px*By - Bx*Py) / (Cx*By - Bx*Cy);
  if (m >= 0) and (m <= 1) then
  begin
    l := (Px - m*Cx) / Bx;
    if (l >= 0) and ((m + l) <= 1) then
      Result := True;
  end;
end;

По коду можно увидеть, что находятся новые координаты точек B и C, которые одновременно являются векторами b и c (рис. 4.). А новые координаты точки P являются вектором (Px, Py). Далее идёт формула, которую я предварительно свёл к общему виду и упростил.

Кол-во основных операций получается 13(+4). Совсем не плохо :]

Метод геометрического луча

Это достаточно известный метод, особенно когда определяется принадлежность точки многоугольникам. Часто данный метод называют «трассировка луча», хотя это не совсем правильно, т.к. трассировка луча — это расчёт хода световых лучей в 3D сцене.

Рис. 5. Метод геометрического луча.
Рис. 5. Метод геометрического луча.

Суть в том, что из данной точки пускается луч по какой-либо оси в каком-либо направлении.
Затем проверяются пересечения со сторонами многоугольника и ведётся подсчёт пересечений.
Таким образом если кол-во пересечений чётное, то значит точка лежит вне многоугольника, если же кол-во НЕчётное, то значит точка лежит внутри.

На рисунке 5 изображены две подопытные точки P и K, у луча из точки P одно пересечение со сторонами треугольника, таким образом точка P принадлежит фигуре, а точке K не повезло — у неё два пересечения.

Реализация алгоритма:

function IsPIn_RayCast(aAx, aAy, aBx, aBy, aCx, aCy, aPx, aPy: single): boolean;
// funcs
  // p1 - начало 1ого отрезка, p2 - конец 1ого отрезка, p3 - начало 2ого отрезка, p4 - конец 2ого отрезка
  function peresechenie(p1, p2, p3, p4: TPoint): boolean;
  var
    zn, ua, ub : single;
  begin // 25 операций
    zn := (p4.y - p3.y) * (p2.x - p1.x) - (p4.x - p3.x) * (p2.y - p1.y);
    ua := ((p4.x-p3.x)*(p1.y-p3.y) - (p4.y-p3.y)*(p1.x-p3.x)) / zn;
    ub := ((p2.x-p1.x)*(p1.y-p3.y) - (p2.y-p1.y)*(p1.x-p3.x)) / zn;
    // если 'ua' и 'ub' принадлежат [0,1] то отрезки пересекаются
    Result := (ua >= 0) and (ua <= 1) and (ub >= 0) and (ub <= 1);
  end;
var
  cros_cnt : integer;
  s1, s2 : single; // коэф.
begin
  cros_cnt := 0; // кол-во пересечений граней (лучом из точки)
  // луч пускаем по оси X вправо
  // 1-я сторона треугла
  s2 := (aPy - aAy) / (aBy - aAy);
  s1 := aAx - aPx + s2 * (aBx - aAx);
  if (s1 >= 0) and (s2 >= 0) and (s2 <= 1) then
    inc(cros_cnt);
  // 2-я сторона треугла
  s2 := (aPy - aBy) / (aCy - aBy);
  s1 := aBx - aPx + s2 * (aCx - aBx);
  if (s1 >= 0) and (s2 >= 0) and (s2 <= 1) then
    inc(cros_cnt);
  // 3-я сторона треугла
  if cros_cnt < 2 then
  begin
    s2 := (aPy - aCy) / (aAy - aCy);
    s1 := aCx - aPx + s2 * (aAx - aCx);
    if (s1 >= 0) and (s2 >= 0) and (s2 <= 1) then
      inc(cros_cnt);
  end;
  Result := cros_cnt = 1;
end;

Вроде нормально упростил(для треугольника имею ввиду), но что-то мне кажется, что может быть и по круче…

А так, получается примерно 30 операций.

Сравнение скоростей

Ну вот мы и подошли к самому интересному! Кто быстрее и сильнее!? :]

Я провёл тест со следующими параметрами(хотя всё зависит от процессора):

  • кол-во повторений алгоритма за 1 имитацию = 4 миллиона.
  • кол-во имитаций для каждого алгоритма = 1000.

Таблица результатов:

Рис. 6. Сравнение скоростей.
Рис. 6. Сравнение скоростей.

Ну что сказать, векторный метод хорош)

С ним конечно соперничает  второй способ относительности точки, но главное отличие в том, что для отн. точки необходима строгая ориентация сторон треугольника, а для векторного это не важно, поэтому он круче)

Так же можно скачать написанную в ходе экспериментов программу: prin_tochki_proga. Программа реализована на Delphi 2007.