Аннотация When constructing a digital terrain model from the available vector data, iso-lines (relief contour lines) are used as the main source of information. Also the maps may contain the point elevation marks. To obtain the relief model, the algorithms for construction of constrained Delaunay triangulations are used, which take the contour lines as the hard edges (constraints). The algorithms may produce the artefacts horizontal "steps consisting of triangles with all the three vertices lying on isolines of the same height. These "steps"appear, both on the ridges, and in the valleys when the neighboring isolines are located suficiently far apart from each other. The article considers the algorithms for elimination of the horisontal triangle artefacts. The proposed algorithm is based upon adding to the triangulation the edges of the approximation of the skeleton (mid-line) of the horizontal triangle groups. To calculate the heights of the added points, a greedy algorithm is used, and we prove its optimality according to the criterion of minimizing the maximum gradient of the added edges. Аннотация При построении цифровой модели рельефа по имеющимся векторным данным в качестве основного источника информации выступают горизонтали (изолинии рельефа). Также на картах могут присутствовать точечные отметки высот. Для получения модели рельефа при этом используются алгоритмы построения триангуляций с ограничениями, в которых отрезки горизонталей выступают в качестве жёстких рёбер (ограничений). В результате применения таких алгоритмов на триангуляции, построенной по изолиниям, появляются артефакты - горизонтальные «ступени», состоящие из треугольников, все вершины которых лежат на изолиниях одной высоты. Эти «ступени» появляются, как на гребнях, так и в распадках в том случае, когда соседние изолинии достаточно далеко отстоят друг от друга. В работе рассматриваются алгоритмы, предназначенные для устранения таких артефактов. Предлагаемый алгоритм основан на добавлении в триангуляцию рёбер аппроксимации скелета горизонтального участка. Для вычисления высот добавляемых точек используется жадный алгоритм, доказывается его оптимальность по критерию минимизации максимального уклона добавляемых рёбер.
При построении цифровой модели рельефа (ЦМР) по имеющимся векторным данным в качестве основного источника информации выступают горизонтали (изолинии рельефа). Также на картах для отдельных вершин могут присутствовать точечные отметки высот. Для получения модели рельефа при этом используются алгоритмы построения триангуляций с ограничениями [1],[3], в которых отрезки горизонталей выступают в качестве жёстких рёбер (ограничений). В результате применения таких алгоритмов на триангуляции, построенной по изолиниям, появляются артефакты - горизонтальные «ступени», состоящие из треугольников, все вершины которых лежат на изолиниях одной высоты. Эти «ступени» появляются, как на гребнях, так и в распадках в том случае, когда соседние изолинии достаточно далеко отстоят друг от друга (Рис. 1, 2, 3, 4).
Хотя во многих случаях для устранения горизонтального участка достаточно выполнить переброску рёбер без добавления новых точек, в общем случае необходимо учитывать, что такой участок может иметь сложную невыпуклую форму, поэтому предлагаемый алгоритм основан на построении аппроксимации скелета горизонтального участка по образующим его треугольникам с последующим вычислением высот в добавляемых точках — вершинах скелетных линий.
Для вычисления высот в добавляемых точках был предложен жадный алгоритм, который сначала находит на графе путь с максимальным уклоном между точками с известной высотой, а затем присваивает высоты точкам с ещё не известной высотой, находящимся на этом пути. Далее процесс повторяется. Доказано, что предложенный алгоритм является оптимальным, т.е. то, что он позволяет минимизировать максимальный уклон ребра. Результат применения алгоритма показан на Рис. 15.
Рис. 1. Артефакты триангуляции, построенной по изолиниям 2
Обзор имеющихся решений Хотя во всех популярных ГИС реализована операция построения ЦМР по изолиниям, готовых решений по устранению горизонтальных участков триангуляций они не предлагают. При этом на проблему горизонтальных участков в разное время обращали внимание различные исследователи и предлагали свои решения. Кратко рассмотрим основные принципы этих решений. Сделаем это не в хронологическом порядке, а по мере возрастания сложности предлагаемых решений.
В работе [4] упоминается проблема появления горизонтальных треугольников, для её решения предлагается добавлять вручную жёсткие рёбра для исправления фрагментов триангуляции, содержащих артефакты. Позже в работе [5] этого коллектива авторов предложены два алгоритма для автоматизации добавления таких ограничений. Алгоритмы рассматривают фрагменты триангуляции, состоящие из горизонтальных треугольников. Простой алгоритм добавляет жёсткие рёбра из вершины за пределами горизонтальной области находящейся рядом с её самым длинным граничным ребром. Добавляемые рёбра должны пересекать самое длинное граничное ребро, но не другие граничные рёбра. Данная эвристика не всегда сработает правильно (самым длинным может оказаться жёсткое ребро), кроме того, никак не контролируется уклон создаваемых треугольников, что может привести к созданию практически отвесных фрагментов с большим углом наклона. Более сложный алгоритм использует идею построения осевой линии из [6], с тем отличием, что в триангуляцию добавляются не просто точки, а отрезки осевой линии. При перемещении по осевой линии в точке ветвления выбирается направление в сторону более длинного ребра, либо в сторону с меньшим углом поворота (заметим, что оба подхода не гарантируют однозначного выбора при изменении начальной точки обхода). В статье не рассмотрен способ вычисления высот в добавляемых точках.
В работе [6] предлагается строить аппроксимации осевых линий (скелетов) горизонтальных областей и добавлять в триангуляцию соответствующие точки. Аппроксимация скелета строится по серединам нежёстких рёбер области, высоты в добавляемых точках вычисляются методом линейной интерполяции. При наличии ветвлений скелета в первую очередь выполняется интерполяция для точек самого длинного пути (что может привести к очень большим перепадам высоты в этих точках ветвления и появлению фрагментов с большой крутизной склонов).
В статье [7] вводится понятие сильных ограничений для триангуляции, представляющей рельеф, основным из которых является требование, чтобы все горизонтальные рёбра проходили по изолиниям или соединяли их конечные точки. Для коррекции также по сути предлагается использовать осевые линии горизонтальных участков. Алгоритм начинает работу с точек ветвления скелета - горизонтальных треугольников, не имеющих жёстких рёбер. Для каждого такого треугольника в триангуляцию добавляется его центральная точка с некоторым изменением высоты. Описан способ определения знака приращения высоты, но не рассматривается способ определения его величины. После обработки всех точек ветвления остаются лишь линейные горизонтальные участки, для конечных точек осевых линий которых высота известна. На этих участках предлагается использовать следующий подход: находим самое короткое ребро, пересекаемое осевой линией и задаём высоту в точке пересечения – середине этого ребра; если высоты в конечных точках осевой линии отличаются, то это – среднее между высотами в конечных точках; иначе (высоты совпадают) считаем точку седловой и вычисляем высоту, как взвешенное среднее высоты конечных точек и высоты граничных рёбер горизонтальной области с весами равными длине отрезка и длине средней линии соответственно; вычисленная точка добавляется в триангуляцию и процесс повторяется. Такой способ вычисления высот придаёт особое внимание наиболее узким участкам горизонтальных фрагментов. При наличии нескольких таких участков или равномерной ширине результат становится непредсказуемым. Кроме того, перепад высот в добавляемой точке также может оказаться значительным, т.к. алгоритм его никак не ограничивает.
Наиболее работоспособные алгоритмы устранения горизонтальных фрагментов рассмотрены в документации к программе DIGIMINE [8], предназначенной для обработки горногеологической информации. Реализован алгоритм триангуляции с соединением точек в пределах прямой видимости, который позволяет устранить часть горизонтальных участков, на которых отсутствуют изгибы, препятствующие проведению таких рёбер. Более сложный алгоритм триангуляции с добавлением дополнительных точек (тальвегов/водоразделов, седловин, вершин/впадин) добавляет в триангуляцию какую-то разновидность скелетных линий. Также имеется алгоритм комбинированной триангуляции, при которой производится и соединение точек в пределах прямой видимости и построение дополнительных точек. Другие подробности работы алгоритмов (способ выбора добавляемых точек и вычисления высот в них) не описаны. К сожалению, не удалось найти публикации с описанием реализованных алгоритмов.
Таким образом, все предложенные алгоритмы используют следующие подходы: перестроение простого горизонтального участка триангуляции (соединение точек в пределах прямой видимости) и добавление точек осевых линий (скелета). Для второго варианта могут различаться способы выбора добавляемых точек и вычисления высот в них.
Рассматриваемый в данной работе алгоритм также использует аппроксимацию осевых линий. Его особенностями являются улучшения, позволяющие получить более качественный результат, приближающийся к диаграмме Вороного отрезков границы, а также строгое обоснование используемого способа вычисления высот в добавляемых точках, как минимизирующего максимальный уклон добавляемых рёбер. 3
Алгоритм исправления артефактов Для исправления артефактов триангуляции будем использовать следующий подход: – Находим области, требующие изменения – горизонтальные участки с примыкающими к ним через нежёсткие рёбра треугольниками – Строим граф - аппроксимацию скелета для каждой такой области. – По конечным вершинам этого графа с известными высотами рассчитываем высоты в промежуточных вершинах. Выбранный способ расчёта определяет результирующий алгоритм.
– Добавляем в триангуляцию точки скелета с вычисленными высотами, а рёбра скелета, соединяющие эти точки, делаем жёсткими рёбрами в триангуляции с ограничениями. Рассмотрим основные этапы работы алгоритма и связанные с ними понятия более подробно. 3.1 Скелет области Рассмотрим множество точек плоскости . Проекция точки на множество - это подмножество точек , ближайших к : Обозначим границу области , как . Скелетом () области является множество точек плоскости, проекция которых на границу области содержит более одной точки: (, ) = ∈ | − | () = { | |(, )| > 1} (1) (2) Внутренним скелетом () области является пересечение скелета () с самой областью: () = () ∩ = { ∈ | |(, )| > 1} (3) Рис. 5. Скелет области 3.2 Аппроксимация скелета фрагмента триангуляции в отдельных треугольниках Скелет многоугольника является графом, состоящим из границ ячеек диаграммы Вороного для точек (вершин многоугольника) и отрезков прямых (рёбер многоугольника). Поэтому скелет состоит из фрагментов кривых, являющихся геометрическими местами точек, равноудалённых от: – двух точек (является отрезком прямой); – двух прямых (является отрезком прямой); – точки и прямой (является фрагментом параболы).
Построение этого скелета является вычислительно сложной задачей, при решении которой никак не используется уже имеющееся разбиение области на треугольники.
Имеющуюся триангуляцию можно использовать для построения аппроксимации внутреннего скелета области, которая должна для большинства триангуляций не слишком отличаться от реального внутреннего скелета. На рис.6 рассмотрены основные виды треугольников и предлагаемые фрагменты графа аппроксимации скелета для этих треугольников. Эти фрагменты строятся без использования информации о соседях треугольника и являются приближением фрагмента настоящего скелета для всех возможных видов и положений этих соседних треугольников. Будем различать граничные и внутренние (общие с другими треугольниками области) рёбра треугольников. Виды треугольников определяются числом граничных рёбер и формой треугольника.
На рис.6 рассмотрены следующие случаи: – угловой треугольник, два ребра которого являются граничными. Для настоящего скелета необходимо провести биссектрису угла между граничными рёбрами, которая, однако, ближе к третьей внутренней стороне может изменить направление или перейти в параболу (в зависимости от того, как продолжаются границы области). Из-за этого, а также для обеспечения совместимости с фрагментами аппроксимации для других видов треугольников вместо биссектрисы используем медиану; – боковой треугольник, имеющий одно граничное ребро. Истинный скелет области в таком треугольнике чаще всего будет содержать фрагмент параболы - геометрического места точек, равноудалённых от граничного ребра и вершины, общей для внутренних рёбер. Эта линия, однако, ближе к краям также может изменить направление. В качестве аппроксимации скелета выбираем отрезок, соединяющий середины внутренних сторон; – центральный треугольник, все рёбра которого являются внутренними. Скелет области для такого треугольника будет содержать фрагменты серединных перпендикуляров к сторонам треугольника, точкой пересечения которых является центр описанной вокруг треугольника окружности. У остроугольного треугольника этот центр попадает внутрь треугольника, а у тупоугольного оказывается снаружи. Два этих случая рассматриваем отдельно: в остроугольном треугольнике проводим серединные перпендикуляры, а в тупоугольном - соединяем середины коротких сторон с серединой длинной (ближайшей к центру описанной окружности точке треугольника).
Случай с тремя граничными рёбрами не рассматриваем, поскольку в этом случае вся область состоит из одного треугольника и не может быть скорректирована предлагаемыми алгоритмами.
Рис. 6. Аппроксимация скелета в треугольниках (жирными линиями обозначены внешние границы области) Для всех видов треугольников фрагменты графа аппроксимации скелета заканчиваются в серединах внутренних сторон, что позволяет корректно совмещать любые фрагменты между собой.
Пример результата применения этого алгоритма для некоторого фрагмента триангуляции показан на рис.8. Во многих местах заметно, что линии скелета оказываются сильно изломанными при чередовании направлений боковых треугольников, но это всё, что можно сделать без привлечения информации о других треугольниках. Аппроксимация скелета фрагмента триангуляции с учётом соседних треугольников Более качественный результат может быть получен, если вместо середины внутренней стороны области в некоторых случаях использовать другую точку этого отрезка, вычисленную с учётом вида и положения обоих треугольников, к которым относится ребро. Поскольку, как уже было сказано ранее, наиболее неудовлетворительный результат получается при чередовании направления боковых треугольников, скорректируем положение точки конца ребра графа в этом случае.
Рис. 7. Точка деления биссектрисой Геометрическим местом точек, равноудалённых от множеств внутренних точек граничных рёбер, является отрезок биссектрисы угла между прямыми, содержащими эти рёбра. В том случае, если такой отрезок биссектрисы пересекает внутреннее ребро, перенесём в эту точку конечную вершину ребра графа (7).
Это изменение затрагивает внутренние рёбра, разделяющие боковые треугольники, граничные рёбра которых касаются разных вершин этого внутреннего ребра в том случае, когда оба угла между внутренним ребром и граничными являются острыми.
Результат применения скорректированного алгоритма показан на рис.9. Рис.10 позволяет сравнить результаты двух алгоритмов. На рисунках заметно, что в большинстве случаев исправленный алгоритм позволил устранить лишние зигзаги на линиях аппроксимации скелета. 3.4
Алгоритмы расчёта высот После построения аппроксимации скелетов для корректируемых областей с примыкающими к ним треугольниками необходимо рассчитать высоты в добавленных вершинах графа. Вычисленные значения должны обеспечивать плавное изменение высоты при переходе по графу между точками с известными высотами. Было предложено два алгоритма, удовлетворяющие этому требованию: – Резиновые нити; – Минимизация максимального уклона.
Рассмотрим общие для обоих алгоритмов понятия и принципы их работы. 4
Линейная интерполяция высоты Рассмотрим некоторый путь, соединяющий две точки с известной высотой. Далее будем использовать следующие обозначения:
Рис. 8. Скелеты горизонтальных областей, полученные по серединам сторон Рис. 9. Скелеты горизонтальных областей, полученные алгоритмом, использующие точки пересечения сторон с биссектрисами ∙ ℎ - высота в точке ∙ , - некоторый путь по графу из точки в точку , для которых ℎ и ℎ – известны ∙ | | - длина пути ∙ ,, = |ℎ−ℎ | – уклон пути между точками и
|| ∙ ∈ - промежуточная точка на пути с неизвестной высотой ∙ = |,| / |,| – параметрическая координата точки на пути ,
При использовании линейной интерполяции высоты значение высоты в точке пути линейно зависит от её параметрической координаты: ℎ = ℎ = ℎ + (ℎ − ℎ) · (4) Рис. 10. Сравнение скелетов горизонтальных областей, полученных без использования (зелёные линии) и с использованием биссектрис (фиолетовые линии)
Рис. 11. Скелеты горизонтальных областей - середины (2) Докажем, что именно линейная интерполяция высоты минимизирует максимальный уклон участка пути. Лемма 1 (Об оптимальности линейной интерполяции) . Любое отклонение от линейной интерполяции высоты в промежуточных точках пути увеличивает максимальный уклон участка пути. Доказательство. Пусть для определённости ℎ ≤ ℎ и в точке высота ℎ = ℎ + , т.е. отклоняется от высоты при линейной интерполяции ℎ на величину . Тогда уклон участка пути от до будет иметь значение: ,, = (ℎ + − ℎ)/ |,| = (ℎ − ℎ)/ |,| + / |,| = ,, + / |,|
Рис. 12. Скелеты горизонтальных областей - биссектрисы (2) Поэтому, при > 0 получаем: ,, > ,,. Аналогично, при < 0 получим ,, > ,,. 5
Минимизация максимального уклона Из двух возможных вариантов задания высот для промежуточных точек графа предпочтём тот, у которого максимальный уклон ребра графа будет меньше. Таким образом, получаемый рельеф будет более плавным, исключающим резкие переходы.
Для того, чтобы найти высоты, минимизирующие максимальный уклон ребра, будем использовать следующий алгоритм, который относится к классу жадных (т.е. выбирающих на каждом шаге локально оптимальное решение) (Алгоритм 1).
Алгоритм 1: Жадный алгоритм минимизации максимального уклона 1 пока существуют вершины графа с ещё не известной высотой выполнять 2 Найти путь с максимальным уклоном между двумя вершинами с известной высотой, проходящий через вершины с неизвестной высотой; 3 Задать высоты промежуточным точкам пути линейной интерполяцией по длине пути;
Жадные алгоритмы не всегда позволяют найти оптимальный результат. Покажем, что в данном случае предложенный Алгоритм 1 позволяет найти высоты для промежуточных точек, обеспечивающие минимально возможный максимальный уклон ребра.
В дополнение к ранее введённым обозначениям далее будем использовать следующие: ∙ *, - кратчайший путь из в ∙ , = ,,*, = |ℎ−ℎ | – уклон кратчайшего пути из в
|*,| ∙ * = (,) , - максимальный уклон пути на графе.
Сначала покажем, что для найденного пути с максимальным уклоном * необходимо выполнить линейную интерполяцию (как это и делается в алгоритме), чтобы не превысить уклон этого пути на одном из его рёбер.
ℎ ℎ ? ? ? ? ℎ ? ℎ Рис. 13. Путь с максимальным уклоном Лемма 2 (О необходимости линейной интерполяции вдоль пути с максимальным уклоном, выполняемой в алгоритме 1). Выполнение линейной интерполяции высоты в промежуточных точках пути с максимальным уклоном * необходимо, чтобы не превысить этот уклон на одном из рёбер графа, при этом максимальный уклон ребра графа не может быть ниже * . Доказательство. Рассмотрим две вершины графа с известными высотами и , между которыми существует один или несколько путей, проходящих по вершинам с ещё не известными высотами. Наибольший уклон пути между двумя этими точками достигается на кратчайшем из этих путей *,. Из Леммы 1 следует, что при задании высот промежуточных точек любое отклонение от линейной интерполяции вдоль этого пути приведёт к тому, что уклон некоторых рёбер этого пути будет превышать уклон , = ,,*, пути в целом.
Таким образом, если не выполнить линейную интерполяцию для пути с максимальным уклоном, то для некоторого ребра этого пути (и для графа в целом) найденный максимальный уклон будет превышен. При линейной интерполяции все рёбра пути будут иметь уклон * и это - тот минимум уклона, который можно достичь для рёбер пути с максимальным уклоном. Значит, максимальный уклон всех рёбер графа не может быть меньше этого значения.
После задания высот для точек пути с максимальным уклоном для других путей, которые его пересекают, произойдёт отклонение от линейной интерполяции высоты. Можно было бы предположить, что в результате может появиться путь с уклоном больше * . Покажем, что это не так. Докажем, что такое отклонение не может привести к превышению значения * на рёбрах путей, затронутых заданием высот для точек пути с максимальным уклоном. Теорема 1 (Oб оптимальности жадного алгоритма) . Пусть (, ) = (,) ,, т.е.
, = * ≥ ,∀(, ) для которых ℎ и ℎ известны (5) Тогда линейная интерполяция высоты в промежуточные точки пути *, не приведёт к возрастанию максимального уклона для оставшихся рёбер. Доказательство. Рассмотрим путь с максимальным уклоном (, ) и некоторый путь (, ), который пересекается с ним в точке (Рис.13). После нахождения пути (, ) его промежуточным точкам, включая точку , присваиваются значения высот, вычисленные при помощи линейной интерполяции. Предположим, что в результате задания значения ℎ стало выполняться , > * = ,, т.е. уклон пути (, ) стал превышать ранее известный максимальный уклон. Пусть для определённости ℎ > ℎ и ℎ > ℎ. Тогда (из сделанного предположения): (6) Но тогда должно выполняться , > ,, действительно: что противоречит условию (5), утверждающему, что максимальный уклон достигается на пути , .
Таким образом, когда жадный алгоритм находит две вершины графа с известной высотой, вдоль кратчайшего пути между которыми достигается максимальный уклон и задаёт высоту в промежуточных точках с использованием линейной интерполяции, то на всех рёбрах этого пути устанавливается тот же найденный максимальный уклон. В силу Леммы 1 именно линейная интерполяция даёт такой результат, любые другие способы задания высот приведут к увеличению уклона на одном из участков. В результате граф делится на связные области из промежуточных точек с ещё не заданными высотами. Следующие шаги алгоритма постепенно находят высоты для этих внутренних точек, при этом максимальный уклон на оставшихся рёбрах может лишь снижаться. 6
Примеры работы алгоритма На рис.14 показан фрагмент триангуляции, полученной по изолинииям, содержащий большое количество артефактов - горизонтальных участков. Слева показано изображение фрагмента без выделения рёбер, в справа – рёбра выделены: жёсткие рёбра – чёрным цветом, остальные - серым. На рис.15 показаны результаты обработки этого фрагмента триангуляции методом минимизации максимального уклона ребра, также слева - без выделения рёбер, а справа - с выделением. На рисунке с выделением рёбер видны добавленные жёсткие рёбра, относящиеся к аппроксимациям скелетов горизонтальных областей.
На рис.16 показан фрагмент триангуляции, содержащей горизонтальные участки, соединяющие разные изолинии и, поэтому, имеющие более одной вершины с высотой, отличной от высоты изолинии, в составе скелетного графа (более одного негоризонтального треугольника, примыкающего к горизонтальной области через нежёсткие рёбра).
Рис.17 позволяет сравнить трёхмерные модели рельефа, полученные по изолиниям, до (слева) и после (справа) применения рассматриваемого алгоритма исправления. Видно, что в результате исправления триангуляции были устранены ступеньки, как на хребтах, так и в распадках рельефа.
Рис. 14. Артефакты триангуляции, построенной по изолиниям Рис. 15. Тот же фрагмент триангуляции после устранения артефактов Рис. 17. Полученная 3D-модель рельефа до и после исправления методом максимального уклона Предложен алгоритм исправления артефактов триангуляций, полученных по изолиниям рельефа и отметкам высот – горизонтальных участков, образованных точками изолиний одинаковой высоты. Алгоритм строит аппроксимацию скелета для каждого горизонтального участка, к которому примыкают негоризонтальные треугольники через нежёсткие рёбра. Для промежуточных точек скелета вычисляются высоты с использованием предложенного алгоритма минимизации максимального уклона ребра. Точки и рёбра скелетов областей включаются в триангуляцию, что позволяет устранить горизонтальные участки.
В дальнейшем алгоритм может быть дополнительно улучшен за счёт разработки методов коррекции аппроксимаций скелетов областей с удалением лишних вершин. Список литературы