Содержание к диссертации
Введение
Глава 1. Построение триангуляции делоне 12
1.1. Введение 12
1.1.1. Диаграмма Вороного и триангуляция Делоне 14
1.1.2. Свойства триангуляции Делоне 16
1.1.3. Структуры данных для представления триангуляции 19
1.1.4. Алгоритмы триангуляции Делоне 20
1.2. Алгоритм триангуляции Делоне на основе предобработки набора исходных точек 26
1.3. Практическая реализация алгоритма 37
1.4. Вычислительный эксперимент 41
1.5. Выводы 45
Глава 2. Приближенное вычисление оптимальной триангуляции 48
2.1. Введение 48
2.2. Локальные приближенные алгоритмы построения оптимальной триангуляции 51
2.3. Экспериментальное исследование локальных перестроений на плоскости 55
2.4. Экспериментальное исследование локальных перестроений на однозначной поверхности 62
2.5. Выводы 69
Глава 3. Построение сечений однозначной поверхности 70
3.1. Введение 70
3.2. Построение коридора для изолиний
3.3. Построение ломаной минимальной длины в коридоре 82
3.4. Сглаживание ломаной минимальной длины кривыми Безье 89
3.5. Сглаживание ломаной минимальной длины локальными сплайнами...91
3.6. Построение визуально гладкой коридорной кривой 96
3.7. Выводы 101
Глава 4. Гладкая интерполяция однозначной поверхности 102
4.1. Введение 102
4.2. Расчет нормальных векторов в узлах сетки 111
4.3. Визуально гладкая интерполяция поверхности 117
4.4. Выводы 123
Глава 5. Реализация и использование алгоритмов построения цифровых моделей рельефа 124
5.1. Введение 124
5.2. Исходные данные и их топологическая корректность 127
5.3. Триангуляция с ограничениями 131
5.4. Быстрая проверка топологической корректности и правильности задания высот 139
5.5. Комплекс программ построения цифровых моделей рельефа
5.6. Преимущества использования предложенных алгоритмов в системах моделирования рельефа 147
5.7. Выводы 150
Заключение 152
Литература
- Свойства триангуляции Делоне
- Экспериментальное исследование локальных перестроений на плоскости
- Сглаживание ломаной минимальной длины кривыми Безье
- Визуально гладкая интерполяция поверхности
Свойства триангуляции Делоне
Триангуляция S получается путем соединения точек из S прямолинейными непересекающимися отрезками так, что любая грань, лежащая внутри выпуклой оболочки S, является треугольником. Поэтому триангуляцией называют получаемую сетку непересекающихся треугольников. При этом граница треугольной сетки совпадает с выпуклой оболочкой, вершинами треугольников являются исходные точки из набора S, а любое ребро будет либо общей границей двух смежных треугольников (внутреннее), либо одним из граничных ребер треугольной сетки. Для упрощения изложения будем далее называть триангуляцией и получаемую треугольную сетку, и сам процесс ее построения.
Если исходный набор включает п вершин, причем g из них принадлежат выпуклой оболочке, то, в соответствии с формулой Эйлера, любая триангуляция данного набора будет содержать ровно 3(п-1)-g ребер и 2{п-\)-g треугольников.
Задача триангуляции набора точек на плоскости является неоднозначной, поэтому необходимо вводить дополнительные критерии для сравнения разных систем треугольников. Триангуляция нужна, прежде всего, для структуризации исходного нерегулярного набора. Точки В{,В2,...,Вк, которые соединяются ребрами сетки с точкой А, определяют ее окрестность и должны быть близкими соседями А. Тогда естественным критерием будет минимум суммы длин ребер триангуляции.
Определение 1.5. Оптимальной или минимальной взвешенной называют триангуляцию, у которой сумма длин всех ребер минимальна.
В работе [111] доказана iVP-полнота задачи оптимальной триангуляции многосвязного полигона, поэтому есть все основания полагать, что и общая задача построения оптимальной триангуляции является не менее сложной. Следовательно, любой алгоритм ее решения будет экспоненциальным по трудоемкости и абсолютно непригодным для обработки реальных наборов точек. Для нахождения субоптимального решения теоретически можно использовать жадный алгоритм триангуляции, трудоемкость которого составляет 0(п2 logn) для п исходных точек [39], однако реально значение п ограничено несколькими тысячами из-за необходимости хранить информацию о Сп потенциальных ребрах треугольников.
Другим недостатком оптимальной триангуляции является возможное вырождение треугольников в линию в случае, когда 3 и более близкорасположенных точек лежат практически на одной прямой (рие. 1.1). В большинстве приложений требуется учитывать не только расстояния между парами точек, но и форму получаемых треугольников, причем желательно, чтобы треугольники были по возможности близки к равносторонним.
Две триангуляции одного набора точек: (а) - треугольники оптимальной триангуляции вырождаются в линию, (б) - треугольники локально наиболее близки к равносторонним
Определение 1.6. Ячейкой Вороного для точки А из заданного набора точек на плоскости называется множество точек плоскости, для которых А является ближайшей точкой из этого набора [39, 140]. Точка А называется центром ячейки.
Ячейка Вороного всегда будет выпуклым многоугольником: замкнутым, если ее центр находится внутри выпуклой оболочки набора, и открытым, если ее центр принадлежит выпуклой оболочке.
Определение 1.7. Диаграммой (мозаикой) Вороного для набора из п точек называется разбиение плоскости на п соответствующих ячеек Вороного [39,140].
В некоторых источниках диаграмму Вороного называют также разбиением Тиссена или Дирихле. Очевидно, что это разбиение единственно.
Границами ячеек Вороного являются ломаные линии. Внутренние точки их ребер лежат на одинаковом расстоянии от центров пар ячеек, а узловые точки одинаково удалены от центров трех или более ячеек. Если каждая из узловых точек является общей граничной точкой ровно трех ячеек, то для диаграммы Вороного можно единственным образом построить двойственную структуру - триангуляцию Делоне. Вершинами этой триангуляции будут исходные точки, а ребра соединят центры ячеек, имеющих общие границы (рис. 1.2а). Если в исходном наборе найдется к 3 точек, лежащих на одной окружности, внутри которой нет других исходных точек, то центром данной окружности будет общая граничная точка к соответствующих ячеек Вороного.
В триангуляцию Делоне потенциально можно включить Ск ребер, соединяющих все пары из данных к точек, но реально войдут: к граничных ребер выпуклого -угольника, вписанного в окружность; любые к-Ъ взаимно непересекающиеся из С\-к оставшихся ребер (рис. 1.26). Определение 1.8. Триангуляцией Делоне называется триангуляция, удовлетворяющая условию Делоне: внутрь окружности, описанной вокруг любого ее треугольника, не попадает никакая из вершин других треугольников [14, 39,44,105,111,129,140].
Экспериментальное исследование локальных перестроений на плоскости
Генерация нерегулярных сеток - вспомогательная операция при математическом моделировании и компьютерном представлении сложных объектов -в настоящее время фактически выделилась в особую междисциплинарную область. Это связано не только с высокой трудоемкостью построения сеток, но и с разнообразием предъявляемых к ним требований, в результате чего сетку часто приходится проектировать под конкретную задачу [34, 94].
Важнейшей областью применения сеток является математическое моделирование на основе метода конечных элементов (МКЭ) [36, 38, 49]. Для повышения точности решения на основе МКЭ обычно используют адаптивные сетки, построение которых включает расчет множества исходных точек, получение начальной сетки и итеративный процесс ее сгущения. На каждой итерации на основе алгебраических, дифференциальных или вариационных методов выделяется множество новых вершин сетки - точек с максимальными значениями градиентов [11, 34, 37, 94]. Процесс построения адаптивной сетки является очень трудоемким и может занимать до 80% общего времени решения задачи [94].
Если адаптивная сетка должна быть треугольной и строится на плоскости или однозначной поверхности, то как при получении начальной сетки, так и при последующем ее сгущении необходимо использовать одну из систематических триангуляции, конечный вид которых (с точностью до отдельных локальных участков) не зависит от порядка задания исходных точек. Систематическими являются 3 триангуляции - оптимальная, жадная и триангуляция Делоне [111, 115, 140].
Оптимальной называют плоскую триангуляцию с минимальной суммой длин ребер (минимальным весом) [39, 111, 115, 140]. При использовании МКЭ или приближении однозначной функции от двух аргументов такая триангуляция может быть предпочтительнее обычно применяемой триангуляции Делоне. Так, в [10] доказано неравенство, определяющее максимальную ошибку приближения гладкой двумерной функции сплайном первой степени на треугольнике. В этом неравенстве величина максимальной ошибки пропорциональна квадрату длины самой длинной стороны треугольника. Поэтому приближение, построенное на оптимальной триангуляции, для значительной части треугольников даст меньшую максимальную погрешность, чем любая другая триангуляция. Кроме того, триангуляцию Делоне не удается обобщить даже на такой простой трехмерный случай, как однозначная поверхность, хотя подобное обобщение оптимальной триангуляции очевидно.
Однако, в отличие от триангуляции Делоне, которую можно построить за время 0(n\ogn) в наихудшем и 0{п) в среднем [39, 44, 105], оптимальную триангуляцию можно вычислить практически лишь за экспоненциальное время. В [111] приведены веские соображения в пользу того утверждения, что задача получения оптимальной триангуляции является ЫР-пошой. С другой стороны, в [98] дан такой пример расположения точек, что вес построенной по ним триангуляция Делоне оказывается в 0{п) раз больше веса оптимальной.
Для получения приближенного решения можно использовать жадный алгоритм [39, 111, 115, 140], который строит триангуляцию следующим образом. Вначале вычисляются веса ребер (расстояния) между всеми п(п-\)12 парами точек. Затем ребра упорядочиваются по возрастанию весов, и первое ребро включается в триангуляцию. Далее в триангуляцию последовательно добавляются те ребра, которые не пересекаются с ранее включенными.
К сожалению, жадная триангуляция вычисляется за большое время 0(п logn), а в процессе работы требуется хранить информацию об п{п -1)/2 потенциальных ребрах. Это реально не позволит обработать больше, чем несколько тысяч точек. Кроме того, жадный алгоритм также не может гарантировать удовлетворительный результат: в [115] приведен пример, когда жадная і/ триангуляция имеет вес в 0{п/ъ) раз больший, чем вес оптимальной триангуляции. Однако примеры размещения точек в [98, 115] весьма искусственны, и в [111] доказано, что для равномерного с некоторыми дополнительными ограничениями распределения точек отношения весов жадной к оптимальной и триангуляции Делоне к оптимальной в худшем случае не превышают 0(\ogn). Ребра жадной триангуляции, расположенные вблизи границы, могут иметь большую длину, поэтому при работе жадного алгоритма приходится проводить проверки пересечения 0(п ) потенциальных ребер с уже включенными в триангуляцию ребрами. Однако большинство потенциальных ребер соединяет пары внутренних точек, которые располагаются на большом расстоянии друг от друга. Такие ребра не попадут не только в жадную триангуляцию, но и в триангуляцию Делоне.
Очевидно, что сумму длин внутренних ребер любого множества смежных треугольников в оптимальной триангуляции уменьшить невозможно. Но для произвольной исходной триангуляции локальные перестроения смежных треугольников могут привести к снижению суммарного веса. Если получаемый вес окажется близким к весу жадной триангуляции, то вместо жадного алгоритма можно использовать алгоритм локальных перестроений, имеющий трудоемкость О(п) и, следовательно, пригодный для обработки сотен тысяч точек. В качестве исходной триангуляции выгодно использовать триангуляцию Делоне, которая: не зависит от порядка задания исходных точек; может быть построена за время О(п); является хорошим начальным приближением оптимальной триангуляции.
В данной главе предлагаются простые локальные алгоритмы, строящие триангуляцию с весом, близким к весу жадной триангуляции. Рассматривается также обобщение этих алгоритмов для триангуляции точек, расположенных на однозначной поверхности z = F(x,y). Результаты, полученные автором, представлены в работах [22, 24].
Сглаживание ломаной минимальной длины кривыми Безье
Будем считать, что границы коридора являются ломаными, которые лежат в той же горизонтальной плоскости, что и искомая изолиния. Для выделения коридора необходима вспомогательная треугольная или прямоугольная сетка. Естественно предположить, что границы коридора и искомая линия должны пересекать одни и те же ребра, поэтому границы можно отслеживать точно так же, как изолинию. На каждом проверяемом ребре АВ необходимо найти отрезок LR, образованный точками, через которые в принципе может проходить изолиния при использовании различных методов интерполяции. Тогда L и R - это пара противолежащих вершин левой и правой границ коридора (рис. 3.2).
Рассмотрим задачу расчета узлов изолинии на треугольной сетке. Пусть плоскость z - zQ пересекает ребро АВ пространственного треугольника. Если предполагается, что поверхность гладкая, а линейная интерполяция не обеспечивает необходимой точности, то АВ нужно заменить таким сегментом кривой, что: данная кривая лежит в вертикальной плоскости, проходящей через АВ; наклоны в точках А и В непрерывны и определяются положением касательной плоскости; функция, описывающая кривую, является однозначной на АВ.
Пусть в вершине А задан нормальный вектор касательной плоскости Ntan (вопросы расчета нормалей касательных рассматриваются в главе 4 настоящей работы). Вектор Nver - {ув -уА,хА -хв,0} - это нормальный вектор вертикальной плоскости, проходящей через АВ. Направляющим вектором касательной прямой (линии пересечения плоскостей Ntan и Nver) в точке А является NА = Ntan х Nver. Отметим, что в силу однозначности поверхности в сетке нет треугольников, у которых все три вершины лежат в одной вертикальной плоскости, поэтому вектор Nver не нулевой, a 7V имеет положи -77 тельную Z-координату. Поэтому вектор NА не может быть вертикальным, т.е. его Х- и 7-координаты не могут одновременно равняться нулю. Пусть в конечных точках ребра ЛВ вычислены направляющие вектора касательных NA = {xA,yA,zA} и NB ={xB,yB,zB}. Для упрощения дальнейших выводов сдвинем систему координат в точку {хА,уА,0} и повернем вокруг оси OZTaK, чтобы точка {хв,ув,0} оказалась на положительной полуоси ОХ. В новой системе координат точки и вектора имеют координаты: A = {0,0,zA}, B = {Xl,0,zB}, NA={x0,0,zA}, NB={x[,0,zB}, где x, = yl(xB-xA)2+(yB-yA)2 ХО= (ХА)2 +(Ул)2 i = (хв)2 + (Ув)2 » причем х0 0,Xj 0. Далее координата Y не используется, и все построения производятся в вертикальной плоскости XOZ.
Отношения zA/x0 и zB\xx определяют наклоны в точках А и В, а по двум наклонам и двум значениям функции zA и zB можно построить отрезок непараметрического кубического сплайна (3.2) или в параметрическом виде Х-составляющая скорости движения точки вдоль данной кривой х (0 = i = const, поэтому общая скорость в точке будет тем выше, чем больше в ней угол наклона, и выпуклость кривой будет зависеть от положения отрезка АВ (рис. 3.3).
«Равномерно выпуклую» кривую можно построить с помощью локальных параметрических кубических сплайнов при условии равенства длин векторов производной на концах отрезка [18]. Для этого необходимо сначала нормировать вектора производных в А и В на результате получены векто ра {kAx 0,kAz A} и [квхх, kBzB}. Увеличивая или уменьшая значения всех четырех координат в равное число раз, можно получать более или менее выпуклые кривые, при этом в А и В сохраняются наклоны, т.е. отношения z /х (рис. 3.4).
Функция z(x) должна быть однозначной на АВ. Для этого необходимо, чтобы условие х (t) О выполнялось при любом t, 0 / 1. Дифференцируя (3.3), получим х (t) = 3a3t2 + 2a2t + ax. Составим систему уравнений для значений x(t), х (t) при t = 0 и t-l с учетом нормировки и, решив ее, вычислим коэффициент аъ = кАх0 + квхх -2х1. Если длину векторов производных увеличить в d раз, то аъ = d(kAxQ + квхх) - 2хх. Но х0 0, хх 0, а х (?) - парабола, поэтому х (t) 0 для любого t, если аъ 0. Следовательно, максимальная допустимая длина модулей производных в А и В составляет "max = 2 1 l\liAXQ + ВХ\ )
Изменяя длину векторов производных в точках А и В от 0 до dmax, можно получать различные однозначные по X кривые для интерполяции поверх -79 ности на ребре АВ. Сегмент на вертикальной плоскости, внутри которого лежат все эти кривые, ограничивают линии, соответствующие значениям d = О (отрезок АВ) и d = dmax (самая выпуклая кривая). Аналогичный по назначению сегмент можно ограничить также отрезком АВ и кривой (3.1). Точки пересечения L и R граничных линий сегмента и плоскости z = z0 определяют отрезок, который пересекают всевозможные линии уровня ZQ , поэтому будем считать L и R противолежащими вершинами двух границ коридора (рис. 3.5а). В более сложных случаях взаимного расположения границ сегмента и секущей плоскости используем следующие правила.
1. Если плоскость z = z0 пересекает граничную кривую, но не пересекает АВ, то будем считать, что пересечение вообще отсутствует.
2. Если граничная кривая имеет точку перегиба и немонотонна, то она может иметь 3 точки пересечения с плоскостью. В этом случае выберем в качестве L и R две наиболее удаленные друг от друга точки из множества точек пересечения (рис. 3.56).
3. Если АВ пересекают две плоскости z = z0 и z = zx, то соответствующие им отрезки L0R0 И L]Rl могут иметь общие точки, например, точки L0 и Rx лежат внутри LXR0. Чтобы коридоры разных изолиний не пересекались, принудительно разделим отрезки, заменяя L0 и Rx точкой Р - центром отрезка L0RX (рис. 3.5в).
Предложенный подход имеет два недостатка: сложность нахождения точек пересечения кривой и плоскости; сложность выделения нескольких отрезков при пересечении одного ребра сетки несколькими плоскостями, особенно на участках перегиба. Кроме того, практика показывает, что для получения приемлемого по качеству набора изолиний процесс их построения должен быть не полностью автоматизированным, а интерактивным.
Визуально гладкая интерполяция поверхности
Современные геоинформационные системы (ГИС) ориентированы, прежде всего, на обработку планарных объектов. Но во многих случаях объекты исследования характеризуются не только привязкой на плоскости, но и значениями пространственно распределенных параметров (например, температуры, уровня загрязнения или просто высот точек). Если каждой точке на плоскости XOY соответствует только одно значение параметра, то такая зависимость описывается некоторой функцией z = F(x,y) и может быть представлена трехмерной однозначной поверхностью. Наиболее сложной по форме и количеству дополнительно задаваемых ограничений, а также наиболее широко используемой поверхностью является земной рельеф. Поэтому лучшим способом проверки работоспособности алгоритмов, представленных в предыдущих главах, будет их использование для моделирования рельефа.
Цифровая модель рельефа - это цифровое представление рельефа с помощью прямоугольной или треугольной сетки, в узлах которой заданы высотные отметки (отметки глубин) [9, 25, 30, 35, 77, 103, 142].
Цифровая модель позволяет получать производные данные для анализа рельефа: вычислять углы наклона и экспозиции склонов, выделять зоны видимости, выпуклости и вогнутости, строить изолинии и профили поперечных сечений, вычислять объемы, расстояния и площади на поверхности, выделять структурные линии рельефа, проводить интерполяцию высот, визуализацию и другие операции [9, 17, 30, 35, 77, 114, 142]. Кроме того, цифровая модель рельефа необходима для построения ортофотопланов местности на основе космо- и аэроснимков [31, 35].
Известные системы построения и обработки поверхностей (в том числе, рельефа) чаще всего реализуются как отдельные компоненты многофункциональных ГИС или систем обработки данных дистанционного зондирования. В качестве примеров можно привести следующие системы.
1. Универсальная ГИС Arclnfo фирмы ESRI (США). Обработка поверхностей, в том числе рельефа, осуществляется с помощью компонент GRID (прямоугольная сетка) и TIN (треугольная сетка), при этом поддерживается переход от одного типа модели к другому. Одна из первых развитых систем построения моделей рельефа, которая во многом стала прототипом для остальных. Модели, построенные с помощью Arclnfo, можно использовать в широко известной системе ArcView GIS, ориентированной на интерактивный анализ пространственных объектов. Дополнительные модули Spatial Analyst и 3D Analyst расширяют возможности визуализации и обработки трехмерных объектов [65, 66].
2. Система обработки данных дистанционного зондирования ERDAS IMAGINE фирмы ERDAS (США). Наиболее мощная из систем дешифрирования и цифрового картографирования с использованием космо- и аэроснимков. Модули цифровой фотограмметрии OrthoMAX и OrthoBASE позволяют, кроме прочего, создать по стереопаре, обработать и отобразить цифровую модель рельефа на основе как регулярной, так и нерегулярной системы точек [86].
3. Полнофункциональная геоинформационная система Maplnfo фирмы Mapping Information Systems (США), широко используемая в России. Для разных ее версий было разработано несколько модулей обработки поверхностей - Vertical Mapper фирмы Northwood Geoscience (Канада), «Поверхность» фирмы ЭСТИ-МАП (г. Москва), SurfMapper (лаборатория ГИС
Томского политехнического университета, г. Томск). Данные модули поддерживают обработку нерегулярных наборов точек и переход к регулярным с помощью целого ряда методов локальной интерполяции, позволяют решать ряд задач анализа рельефа [17, 61, 137].
4. Пакет трехмерного моделирования поверхностей Surfer фирмы Golden Software (США). На основе различных методов интерполяции позволяет создавать равномерную сеть для отображения неравномерно распределенных данных и аналитических функций с учетом разломов и структурных линий. Используется для построения контуров, векторов, затененных изображений рельефа и проволочных моделей поверхности. Применяется в научных исследованиях [134].
5. Система обработки карт и изображений TNTmips компании Microimages (США). Используется для анализа, обработки и автоматизированного дешифрирования материалов дистанционного зондирования, фотограмметрической обработки изображений, цифровой картографии, геофизических и геологических приложений [135].
6. Интегрированная геоинформационная система ГрафИн, разработанная в НПО «Сибгеоинформатика» (г. Томск). Развитая подсистема моделирования рельефа позволяет использовать регулярную и треугольную сетки, строить разрезы и изолинии различных типов, вычислять расстояния, площади и объемы земляных работ, отображать поверхность [9, 46].
Компоненты для работы с поверхностями в этих системах функционально во многом близки. Анализ их возможностей позволил выделить множество задач, которые приходится решать при построении и обработке цифровых моделей рельефа. Полный перечень этих задач приведен в п. 5.6, где исследуется возможность их решения на основе алгоритмов, предложенных в настоящей работе.
