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

Question

Есть код, который проверяет принадлежит ли точка многоугольнику

function IsInPoly(const ARegion: array of TPoint; const ASearchPoint: TPoint): Boolean;
var
  LPrevPoint, LCurPoint: TPoint;
  Li: Integer;
  DX, DY, StepX, StepY: Integer;
  Step: Double;
  Rotation: Double;
begin
  LPrevPoint := ARegion[0];
  Rotation := 0;
  for Li := 1 to Length(ARegion) - 1 do begin
    LCurPoint := ARegion[Li];
    DX := LCurPoint.X - LPrevPoint.X;
    DY := LCurPoint.Y - LPrevPoint.Y;
    StepX := LPrevPoint.X - ASearchPoint.X;
    StepY := LPrevPoint.Y - ASearchPoint.Y;
    if (StepX <> 0) or (StepY <> 0) then begin
      Step := (DX * StepY - DY * StepX) / (StepX * StepX + StepY * StepY);
      Rotation := Rotation + Step;
    end;
    LPrevPoint := LCurPoint;
  end;
  Result := (Abs(Rotation) > 5);
end;

Помогите идентифицировать этот алгоритм. Я не вижу, чтобы что-то подобное было описано в тех статьях, которые гуглятся. В итоге я не понимаю, работает ли код в принципе?

Тут, как я понял (DX, DY) - это вектор от предыдущей точки к следующей, а (StepX, StepY) от искомой к предыдущей. А что такое Step и Rotation?

Answer 1

Нет, алгоритм неправильный.

Контрпример — квадрат со углами (1, 0), (0, −1), (−1, 0), (0, 1), проверяемая точка (0, 0).

В этом случае (StepX * StepX + StepY * StepY) всегда равно 1, и мы суммируем площади треугольников OAᵢAᵢ₊₁, которые складываются в 4, что меньше 5.

(Алгоритму, судя по всему, нужно дублировать начальную точку в конце.)

Насколько я понимаю, автор алгоритма хотел сложить ориентированные углы AᵢXAᵢ₊₁ для всех i, эта сумма должна давать теоретически 0 для точки снаружи и 2π для точки внутри. Правильная реализация должна по идее выглядеть так:

DX1 := LCurPoint.X - ASearchPoint.X;
DY1 := LCurPoint.Y - ASearchPoint.Y;
DX2 := LPrevPoint.X - ASearchPoint.X;
DY2 := LPrevPoint.Y - ASearchPoint.Y;
S := DX1 * DY2 - DX2 * DY1; // ориентированная площадь, произведение сторон на синус
P := DX1 * DX2 + DY1 * DY2; // скалярное произведение, произведение сторон на косинус
if (S <> 0) or (P <> 0) then begin
  Angle := ArcTan2(S, P);
  Rotation := Rotation + Angle;
end;

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

Ещё популярные алгоритмы:

Вычислять последовательно векторы сторон AᵢAᵢ₊₁ многоугольника, и брать векторное произведение AᵢAᵢ₊₁ на XAᵢ (X — тестируемая точка). Если у всех этих векторных произведений одинаковый знак (геометрически: точка X на пути вдоль каждой стороны всё время слева (или справа), то точка внутри. Алгоритм работает только для выпуклых многоугольников. (Возможно, его можно обобщить для произвольных многоугольников, подсчитав чётность смен знака, доказать с лёта не могу.)
Провести через точку X произвольную прямую (например, горизонтальную), и подсчитать точки пересечения с каждым из отрезков, составляющих стороны. Если слева (или справа) от точки X чётное количество точек пересечения (например, ноль), то мы снаружи, или нечётное, то внутри. Чтобы не разбираться со случаями, когда прямая проходит через вершину, можно выбрать не обязательно горизонтальную прямую, любую не совпадающую ни с одной из прямых XAᵢ (их конечное число). Это решение в духе теоремы о жордановой кривой, работает для любого несамопересекающегося многоугольника.

БЛОГ НА HUSL

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

Ответы (1 шт):