Точки, полигоны и функторы

[thedeemon]

Долго колебался, стоит ли об этом писать, но после того, как знающий человек похвалил выбор алгоритма, решил все же рассказать. Интересная часть экстремистской задачи из прошедшего конкурса ПФП про вырезание московской области из России заключалась в том, чтобы определить для заданной точки, входит ли она в фигуру, заданную набором ломанных (набор невыпуклых многоугольников, в которых могут быть дырки из других многоугольников). Сперва меня терзали сомнения, можно ли использовать подходы из планиметрии в задаче, где координаты точек даны в терминах долготы и широты на совсем не плоской планете. Ведь если посмотреть на карту в районе полюса, становится хорошо заметно, что прямые в этой системе координат не такие уж прямые. Но потом забил на эти сомнения и стал решать планиметрически.

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



Для начала определим тип для отрезков:

type line = { x1 : float; y1: float; x2 : float; y2 : float } 

Когда структура целиком состоит из флоатов, они хранятся компактно, unboxed, что хорошо для производительности.

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

let intersects x y line = 
  if x < line.x1 || x >= line.x2 then 0 else 
  let ly = (x -. line.x1) /. (line.x2 -. line.x1) *. (line.y2 -. line.y1) +. line.y1 in 
  if ly >= y then 1 else 0

Функция возвращает не bool, а 0 или 1, чтобы результаты по всем отрезкам можно было просуммировать и узнать, четно их количество или нет. Граница московской области содержит около 1500 отрезков, и если проверять их все, то получается довольно медленно - за секунду получается проверить около 20000 точек, в то время как на карте России их около 10 миллионов. Кто-то на Common Lisp'е решал эту проблему путем кодогенерации через макросы, что уменьшает время одного шага, но не количество шагов. Я же решил переходом от линейного алгоритма к логарифмическому. Для этого при чтении отрезков фигуры их х-координаты собирались в массив, он сортировался, повторы выкидывались и получался набор интервалов, на которые проекции точек фигуры разбивают ось Х (см. картинку). Для каждого такого интервала запоминался набор невертикальных отрезков, проекции которых на Х с ним пересекаются. Эти данные помещались в сбалансированное дерево, реализованное с помощью входящего в стандартную библиотеку Окамла функтора Map. Функторами в Окамле называются параметризованные модули, по сути это функции над модулями (где модуль - набор типов и функций). Map требует на вход модуль, описывающий некий тип, для которого задана операция сравнения. Нужно было лишь задать такой тип и операцию для интервалов:

type interval = { left : float; right : float } 
module Interval = struct  
  type t = interval  
  let compare a b =  
    if a.right <= b.left then -1 else 
    if b.right <= a.left then 1 else 0 
end   
module M = Map.Make(Interval)

Применение функтора Map.Make к модулю Interval дает модуль M, который реализует деревья с интервалами в качестве ключей и произвольным типом в качестве значений. В моем случае это был массив из отрезков (точнее, фактически, ссылок на них). Поиск в таком дереве - операция логарифмической сложности. В результате функция проверки принадлежности точки фигуре выглядела так:

let is_inside figure x y = 
  try 
    let a = M.find { left = x; right = x } figure in 
    let n =  DynArray.enum a |> Enum.map (intersects x y) |> Enum.fold (+) 0 in 
    n mod 2 > 0 
  with Not_found -> false

Вот и вся реализация. Она уже работала со скоростью около 2 миллионов точек в секунду, т.е. вся карта России могла быть обработана за 5 секунд. Дальше все упиралось лишь в скорость чтения и разбора карты и записи результата.
По-моему, гораздо элегантней, чем кодогенерация или разбиение карты на клетки и заливка floodfill'ом, как у кого-то.

Comments

[enternet.livejournal.com]

Я когда эту задачу думал делать, этот алгоритм на алголисте подсмотрел. http://algolist.manual.ru/maths/geom/belong/poly2d.php

[justy-tylor.livejournal.com]

Да, идея сведения к одномерной сетке хорошая и простая.

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

[thedeemon.livejournal.com]

Вроде, полюса не входят ни в одну карту страны, так что проблемы особой нет.

[usovalx.livejournal.com]

Сам полюс правильно обработать довольно нетривиально, а всё остальное пространство карты можно рассматривать как плоскость, как описал автор. Hint -- поискать в сети карту земли американского образца. Они используют прямоугольную сетку.

[zevlg.livejournal.com]

ээээ, и это работает для впуклых многоугольников?

[thedeemon.livejournal.com]

Конечно.

[sastanin.livejournal.com]

Красивая идея с подсчётом пересечений. Спасибо.

[vp.livejournal.com]

красота. Спасибо! Запомним на будущее.

[soonts.livejournal.com]

Думается мне, твой алгоритм пойдёт по пизде в случаях, когда тестируемая точка является вершиной.

Если пофиксить функцию intersects шоб она возвращала 0 в случае когда точка на отрезке, из-за особенностей float-арифметики алгоритм пойдёт по пизде в случаях, когда тестируемая точка близка к точке пересечения сегментов ломанной, не являющейся вершиной ломанной.

Правильное решение для общего случая сука сложное и довольно медленное, что-то вроде этого только 2D..

[thedeemon.livejournal.com]

Насколько я помню из обсуждения задачи, на граничные случаи можно было забить.

[sleepy-drago.livejournal.com]

на алгоритм проблемы плавучки тень не бросают. Обычно или считают в integer с заданной точностью :) или забивают на эти случаи. Заниматься refinement'ом деревьев и split'ами в граничных случах можно всю жизнь :) и оно так и не будет работать в некоторых случаях на практике.