Упрощение линий Рамера — Дугласа — Пекера

Над Алматы, в предгорьях и ущельях Заилийского Алатау, OpenStreetMap знает 1556 троп и грунтовок. Часть из них нарисована по GPS-трекам пешеходов — так и помечено в данных, — часть по спутниковым снимкам. Всего в них 58 758 точек: в среднем одна на 15 метров, а на самых подробных тропах — через каждые два-три метра, со всеми поворотами серпантина. Чтобы показать их на карте телефона или хранить в приложении, столько не нужно: на экране точки в паре метров друг от друга сливаются в одну. Двигайте ползунок допуска — все тропы упрощаются сразу, а выбранная показывает свой коридор: ни одна исходная точка не отходит от упрощённой линии дальше допуска.

При допуске пять метров от 58 758 точек остаётся 17 621 — в 3,3 раза меньше, а тропа «3 брата — Кумбель» из 450 точек превращается в 33. Если взять столько же точек, просто оставляя каждую пятнадцатую, линия отойдёт от тропы на 14,5 м, а у этого алгоритма — не больше чем на 4,9 м. Он выбирает, какие точки важны для формы, и тратит их там, где тропа поворачивает.

Идея

Оставим первую и последнюю точку и проведём между ними хорду. Найдём точку, которая дальше всех от хорды. Если она ближе допуска, все точки между концами можно выбросить: хорда и так проходит рядом с каждой. Если дальше — эта точка нужна; оставим её и сделаем то же для двух половин линии: от начала до неё и от неё до конца.

Это алгоритм Рамера — Дугласа — Пекера. Допуск ε — единственный параметр, и у него понятный смысл: насколько упрощённая линия может отойти от исходной. Расстояние здесь — от точки до отрезка хорды: до перпендикуляра, если он падает на отрезок, иначе до ближнего конца. В классическом описании берут расстояние до прямой, но на серпантинах это опасно: точка далеко за концом хорды может оказаться рядом с её продолжением и пропасть. Ниже — алгоритм на одной тропе, хорда за хордой.

Альпинистский маршрут 1Б на пик Школьник — 392 точки на 679 метрах, по точке на каждые 1,7 м. С допуском 5 м первая хорда соединяет начало и конец маршрута, а дальше серпантин раскладывается на всё более короткие хорды, пока каждая не окажется в пределах допуска от своих точек. Остаётся 24 точки.

Почему форма не портится

Гарантия односторонняя, но её хватает: вокруг упрощённой линии можно провести коридор шириной 2ε, и вся исходная линия в нём. Чего алгоритм не обещает — так это наименьшего числа точек. Он жадный сверху вниз: оставив точку, больше не пересматривает решение. Наименьшее число точек при том же допуске можно найти и точно: построить граф, где ребро i → j есть, если отрезок i–j проходит не дальше ε от всех точек между ними, и найти в нём кратчайший путь из первой точки в последнюю. Такие точные методы разбирали Хироси Имаи и Масао Ири в 1986 году, но даже их лучшие потомки тратят порядка n² действий — проверить приходится все пары точек. Рамер — Дуглас — Пекер на настоящих линиях намного быстрее, но может оставить больше точек, чем нужно. Не обещает он и наименьшего отклонения при заданном числе точек: при очень грубом допуске, когда от тропы остаются три-четыре точки, равномерная выборка иногда ложится ближе — в первой сцене это видно на допуске 50 м.

Код

Рекурсия здесь естественна, но на длинной линии она может уйти слишком глубоко, поэтому в листинге вместо неё стек хорд. Пример — короткий серпантин из тринадцати точек: программа печатает, что она решила на каждой хорде. Вторая часть — проход, который пригодится ниже: один запуск с нулевым допуском даёт каждой точке её порог.

А это первая сцена в цифрах: все тропы сразу, допуски от метра до пятидесяти, длина до и после, и сравнение с «каждой k-й точкой».

Ловушка: серпантин становится короче

Упрощённая линия срезает углы. Для формы на карте это и нужно, а для длины — нет: хорда всегда короче дуги, которую она заменяет. На всей сети троп потеря небольшая: 895,5 км превращаются в 884,1 км при допуске 5 м, на 1,3% меньше, и в 835,1 км при 50 м. Но на серпантине каждый поворот — это несколько метров, которые хорда срезает целиком. Маршрут на Школьник длиной 679 м после упрощения с допуском 5 м — 564 м, на 17% короче, а с допуском 20 м — 496 м. Приложение, которое посчитает длину похода по упрощённому треку, пообещает на четверть меньше ходьбы, чем будет на самом деле.

Отсюда правило: длину, набор высоты и время считают по исходному треку, а упрощают только для рисования и хранения формы. Вторая известная беда — топология: упрощённая линия может пересечь сама себя или соседнюю, хотя исходные не пересекались. Мы проверили тропы этой статьи: ни одна из них не пересекла сама себя ни при 5, ни при 50 метрах — расстояние до отрезка здесь осторожнее расстояния до прямой. Но для границ участков, где пересечение — ошибка в данных, используют варианты алгоритма, которые это проверяют.

Все масштабы одним проходом

Карта показывает тропы на разных масштабах, и каждому нужен свой допуск — примерно в пиксель. Запускать алгоритм заново на каждый масштаб не обязательно. Дерево делений хорд от допуска не зависит: самая далёкая от хорды точка одна и та же при любом ε. Меняется только то, где деление останавливается. Поэтому один проход с нулевым допуском записывает каждой точке её расстояние до хорды, а при допуске ε остаются точки, у которых это расстояние больше ε у них самих и у всех их предков по дереву.

Слова про предков важны. Расстояния вниз по дереву не обязаны убывать: в примере из листинга точка 5 отстоит от своей хорды на 11,70, а её «родитель», точка 2, — всего на 7,15. При допуске 8 точку 2 выбросят, хорда 0–12 останется целой, и до точки 5 дело не дойдёт. Поэтому порог точки — наименьшее из расстояний на пути от корня: min(11,70; 7,15) = 7,15. Листинг проверяет, что такой фильтр даёт ровно то же, что полный запуск, — на всех допусках. Первая сцена так и работает: ползунок только сравнивает пороги с допуском. Похожий приём есть в TopoJSON: каждой точке заранее записывают вес — площадь треугольника с соседями, как в алгоритме Висвалингам — Уайатта, — а упрощение под масштаб выбрасывает точки с весом ниже порога.

Сколько это стоит

Каждая хорда — проход по точкам между её концами. Если самая далёкая точка делит линию примерно пополам, уровней деления около log₂ n, и на каждом в сумме O(n) работы — всего O(n log n). Так ведут себя настоящие тропы: во второй сцене видно, как быстро укорачиваются хорды. Но если самая далёкая точка всякий раз стоит рядом с концом хорды, деление отщипывает по одной точке, и выходит 1 + 2 + … + n ≈ n²/2. Такая линия есть: зигзаг, сходящийся к своему концу, где каждая следующая вершина чуть ближе к оси, чем предыдущая.

На тысяче точек зигзаг обходится в десятки раз дороже тропы, и разрыв растёт вместе с n. Джон Хершбергер и Джек Сноуинк в 1992 году показали, как сделать и худший случай O(n log n): самая далёкая от хорды точка всегда лежит на выпуклой оболочке куска линии, и если поддерживать оболочки кусков при делении, искать её можно быстро. Памяти нужно O(n): сами точки, флажки «оставить» и стек хорд, в котором не больше n записей.

Где работает

• ГИС и картография. Функция ST_Simplify в PostGIS, simplify в Shapely и DouglasPeuckerSimplifier в библиотеке JTS упрощают геометрию этим алгоритмом.
• GPS-треки. Часы и телефоны часто пишут трек пробежки или похода точкой в секунду, а для карты и хранения хватает точек на поворотах.
• Этот каталог. Формы улиц на карте Алматы в сценах с Дейкстрой и A* упрощены им с допуском 3 м — это видно в скрипте подготовки данных.
• Роботы. Лазерный дальномер даёт облако точек стен, и тот же приём — «разделить в самой далёкой точке» — выделяет из него прямые отрезки.

История

Алгоритм придумали дважды и почти одновременно. Урс Рамер опубликовал его в 1972 году в журнале Computer Graphics and Image Processing — как способ приближать плоские кривые ломаными. Годом позже Дэвид Дуглас и Томас Пекер описали в картографическом журнале тот же приём для оцифрованных линий на картах; их статья называлась «Алгоритмы сокращения числа точек, нужных, чтобы представить оцифрованную линию или её карикатуру». В литературе по распознаванию образов тот же метод известен как iterative end-point fit — подгонка по концам с повторением.

Задачи

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