Предварительная обработка исходных данных для построения цифровой модели рельефа местности Текст научной статьи по специальности «Математика»

Ю.Л. Костюк, А.Л. Фукс

ПРЕДВАРИТЕЛЬНАЯ ОБРАБОТКА ИСХОДНЫХ ДАННЫХ ДЛЯ ПОСТРОЕНИЯ ЦИФРОВОЙ МОДЕЛИ РЕЛЬЕФА МЕСТНОСТИ

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

Важным объектом исследования современных геоин-формационных систем является земной рельеф. Как правило, исходные данные для построения цифровой модели рельефа (ЦМР) получаются с помощью методов дистанционного зондирования и представляют собой наборы высотных отметок и структурных линий, задающих дополнительные ограничения на форму рельефа. При этом координаты высотных отметок и узлов некоторых структурных линий (границ оврагов, обрывов, береговых линий) задаются с высокой точностью, поэтому целесообразно использовать цифровую модель на основе треугольной сетки, узлами которой являются указанные точки. Если для построения сетки используются алгоритмы триангуляции с ограничениями, то в получаемую модель просто и естественно включаются структурные линии.

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

В настоящей работе для обработки структурных линий рельефа предлагается использовать метод геометрической аппроксимации сечений однозначной поверхности [1, 2]. В соответствии с ним для каждой исходной линии необходимо выделить коридор - область допустимого расположения данной линии, а затем внутри коридора построить по возможности плавно изменяющуюся расчетную линию, которая и будет заменять исходную.

ВЫДЕЛЕНИЕ КОРИДОРА ДЛЯ СТРУКТУРНОЙ ЛИНИИ

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

Корректно устранить топологические нарушения можно только интерактивно. Далее мы будем предполагать, что эта работа уже проведена. Если существует допустимое пересечение двух структурных линий, то при проверке ограничений точка пересечения должна автоматически включаться в каждую из этих линий и, естественно, стать узлом триангуляции. В результате будет действительно получена треугольная сетка, для которой все отрезки исходных структурных линий совпадают с ребрами треугольников (рис. 1).

Рис. 1. Триангуляция со слабыми ограничениями на основе исходных изолиний

Такая сетка с заданными в ее узлах высотами используется в большинстве известных геоинформаци-онных систем в качестве ЦМР. Однако, как показано в [6], эта модель является недопустимо грубым приближением и требует дополнительного уточнения. В настоящей статье треугольная сетка применяется только как вспомогательная структура для выделения коридоров исходных линий.

При любом методе получения структурных линий для каждого типа линий (горизонтали, линии обрывов, береговые линии) можно определить два значения абсолютной погрешности: еХУ - точность задания координат узлов в плоскости ХОУ и ег - точность расчета высот в узлах. Рассмотрим участок структурной линии АВС и несколько примыкающих к нему треугольников сетки, расположенных справа по ходу движения вдоль линии (рис. 2).

Рис. 2. Расчет узловых точек границы коридора

Будем считать, что точка R на ребре AU принадлежит коридору линии ABC, если выполняются три условия:

V(ХЯ - xA ) + (уЯ - yA ) - е XV ,

\2Я - ZAІ — ^ , (1)

V(ХЯ - xA ) + (уЯ - yA ) — "2>/(Хи - XA ) + (уи - yA ) ,

т.е. расстояние между R и A как в плоскости ХОУ, так и по высоте не превышает величины погрешности, и, кроме того, точка R находится ближе к A, чем к и. Пусть - это наиболее удаленная от A точка на AU, для которой выполняются условия (1), тогда принадлежит границе коридора. Проходя по смежным треугольникам, расположенным вдоль линии ABC справа, можно последовательно вычислять узловые точки правой границы коридора Яш , Я^, ЯВУ,

ЯВШ. Узлы левой границы коридора рассчитываются аналогично путем просмотра треугольников, примыкающих к ABC слева. Пример коридора для изолинии приведен на рис. 3.

Рис. 3. Выделение коридора для изолинии на треугольной сетке

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

Расчет линии, целиком лежащей внутри коридора, требует проверок взаимного расположения отдельных участков коридора и линии. Однако выделение нужного участка коридора, т.е. пары начальных и пары конечных узлов участка границ, будет весьма сложной задачей, так как границы коридора вычисляются независимо, а число их узлов может существенно различаться. Данную проблему можно решить путем использования кратных узлов границ. Приведем алгоритм модификации границ коридора для незамкнутой структурной линии, выходящей на границы треугольной сетки. Этот алгоритм фактически строит внутри коридора триангуляцию Делоне.

Алгоритм модификации границ коридора

Предполагается, что заданы наборы из т узлов ¡г левой и к узлов Ту правой границы коридора. Новые

границы Ь и Я формируются путем последовательного добавления к ним узлов исходных границ. На каждом шаге выбираются 2 текущие точки (¡г и Ту) и 2

следующие за ними (¡г+1 и т]+1). Далее для треугольника ¡¡ТуТ]-+1 и точки ¡г+1 проверяется условие Делоне (¡г+1 находится вне окружности, описанной около треугольника). При его выполнении текущими п-ми вершинами новых границ становятся ¡г и Ту+1, в противном случае треугольником Делоне будет ¡т^м , а в новые границы добавляются точки ¡м и т'у.

I := 1; у := 1; п := 1; Ь[1]:= ¡[1]; Я[1] := т[1];

пока (г < т) или (у < к ) цикл

п := п +1;

если г = т то начало

Ь[п] := ¡[/]; Я[п]:= т[у +1]; у := у +1;

конец

инес у = кто начало

Ь[п] := ¡[г +1]; Я[п] := т[у ]; г := г +1; конец

инес ¡1+1 вне описанной окружности IТуТу+1 то начало

Ь[п] := ¡[г]; Я[п]:= т[ у +1]; у := у +1; конец

иначе начало

Ь[п] := ¡[г +1]; Я[п] := т[у ]; г := г +1;

конец

кц

Конец алгоритма

Данный алгоритм имеет линейную трудоемкость О(к + т). Новые границы включают одинаковое число п узловых точек, п = т + к -1, причем для любых р и q, 1 - р — д — п , узлы Ьр,...,Ьд на левой границе и

Яр,...,Яд на правой соответствуют одному участку

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

ПОСТРОЕНИЕ ЛИНИИ ВНУТРИ КОРИДОРА

Расчетная структурная линия, которая будет заменять исходную, должна целиком лежать внутри коридора и изменяться по возможности плавно. Прежде всего, необходимо избавиться от лишних узловых точек и уменьшить осцилляцию линии. В [1] в качестве начального приближения расчетной линии предлагается использовать ломаную линию минимальной длины с заданными начальной и конечной точками, а также приводится эффективный алгоритм ее построения. Эта ломаная по форме подобна резиновой нити с двумя закрепленными концами, растянутой внутри

коридора, поэтому она имеет минимальное число осцилляций (рис. 4).

Рисунок иллюстрирует характерное свойство ломаной минимальной длины - на участках выпуклости (вогнутости) она прижимается к левой или правой границе коридора, а на участках перегиба переходит с одной границы на другую. С этим связан и основной ее недостаток: для замкнутых выпуклых линий она всегда совпадает с одной из границ.

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

- Построение коридора и ломаной минимальной длины в нем.

- Модификация коридора, при которой все вершины границ Ьг или Яг , совпадающие с узлами ломаной, сдвигаются к центру коридора (к центру соответствующих отрезков ЬгЯг).

- Построение ломаной минимальной длины в модифицированном коридоре.

В результате узлы полученной ломаной располагаются, в основном, в центре начального коридора, число осцилляций увеличивается незначительно, а порядок трудоемкости построения ломаной не увеличивается. Рис. 5 иллюстрирует процесс получения данной линии.

Рис. 5. Модификация коридора и перестроение ломаной минимальной длины

К сожалению, ломаная минимальной длины визуально очень далека от гладкой, поэтому не может использоваться непосредственно в качестве расчетной структурной линии. Однако по ее узловым точкам с помощью локальных параметрических кубических

сплайнов можно построить гладкую интерполирующую кривую, которая будет иметь минимальное количество точек перегиба и целиком располагаться в коридоре [1, 2].

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

В вычислительном отношении наиболее простой способ интерполяции - это замена отрезков ломаной минимальной длины кубическими кривыми Безье [3]:

р(*)=Е 3=0 р^ (*), 0—*—1,

где (*) = С\? (1 - * )3-г - кубические полиномы

Бернштейна, а Рг - вершины характеристического многоугольника.

Пусть Му, Му+1 - две соседние вершины минимальной ломаной. Кривую Безье на у-м отрезке определят четыре вершины многоугольника, причем Р0 = Му, Р3 = Му+1, а дополнительные точки Р1 и Р2

задают наклоны в Му и Му+1. Если на (у -1 )-м участке кривую Безье определяют точки Q0, Q1, Q2 и Q3, то, естественно, Q3 = Р0 = Му, а для обеспечения непрерывности наклона в Му точки Q2, Му и Р1

должны лежать на одной прямой.

Наклоны кривой Безье в каждой вершине минимальной ломаной можно вычислить, используя простые методы для локальных кубических сплайнов [4]. Затем на прямых, определяемых этими наклонами, необходимо так расположить точки Р1 и Р2 , чтобы кривая целиком попадала в участок коридора, соответствующий отрезку М уМ у+1 ломаной минимальной

длины (рис. 6).

Рис. 6. Расчет отрезка кривой Безье в коридоре: а - на выпуклом участке, б - на участке перегиба

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

ПОСТРОЕНИЕ ВИЗУАЛЬНО ГЛАДКОЙ КОРИДОРНОЙ КРИВОЙ

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

Для решения данной проблемы предлагается заменять кривую с п узлами ломаной с предельным числом вершин N > п , которая будет визуально настолько близка к гладкой линии, насколько это позволяет значение N.

Будем считать, что ломаная является визуально гладкой, если любые три ее последовательные вершины кажутся лежащими на одной прямой. Очевидно, что это свойство определяется масштабом изображения. Для любого масштаба 1: М и толщины изображаемых линий Д¡ (в единицах задания координат точек) можно вычислить ДЬ = Д¡ • М - реальное расстояние, соответствующее толщине линии на изображении. Если A, В, С - последовательно расположенные вершины ломаной и расстояние от В до отрезка AC меньше ДЬ , то на рисунке с масштабом 1: М данные точки будут визуально располагаться на одной прямой.

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

Назовем отклонением вершины р расстояние от

нее до отрезка Рг-1 Рм , соединяющего две соседние с Рг точки ломаной (отклонение начальной и конечной вершин незамкнутой линии равно нулю). Пусть Р-2, Р— , Р, Р+! и Р+2 - последовательно выбранные вершины ломаной, которая интерполирует заданную кривую. Если отклонение Р превышает ДЬ , то ломаную необходимо перестроить, добавив в нее либо дополнительную вершину Ц на отрезке Р^Р , либо Ь>2 на РР+1. В качестве Д следует выбрать наиболее удаленную от хорды Рг-1Рг точку на соответствующем отрезке кривой, Ь>2 выбирается аналогично

на РР+1 (рис. 7).

Включение Д приводит к изменению отклонений в точках Р-1 и ^ , а также к необходимости вычис-

ления отклонения В1 от РІ-1РІ. Аналогично, три значения отклонений нужно вычислить и при добавлении точки Ь>2. Включение новой точки должно максимально приближать ломаную к кривой, поэтому выбор Д или Ь>2 будет определяться минимумом значения максимального из трех отклонений.

Рис. 7. Включение новых точек при построении визуально гладкой ломаной

В последующем описании номера точек всегда будут определять не порядок их следования в ломаной, а порядок добавления вершин к ломаной. Начальные п точек - это вершины сглаживаемой ломаной минимальной длины. Добавляемые вершины будут иметь номера п +1,...,N .

Новую точку необходимо включать всегда в окрестности вершины с максимальным отклонением. Для упрощения выбора такой вершины предлагается по значениям отклонений точек строить пирамиду (дерево сортировки), как в алгоритме пирамидальной сортировки [5]. Пирамида для п значений представляет из себя массив длины п, в котором для любого элемента аг при г < п /2 выполняется: аг > а2г,

аг > а2г+1. Таким образом, элемент а1 является максимальным. Для построения такого массива и добавления новых элементов используется операция просеивания элемента в пирамиде.

При построении визуально гладкой ломаной линии потребуются следующие структуры данных:

Пирамида, каждый элемент которой содержит отклонение (используемое при просеивании) и номер соответствующей точки. Пирамида включает от п (начальная) до N (максимальная) элементов.

Текущая ломаная. Для простоты модификации представляется двумя массивами, в которых для каждой точки хранятся номера предыдущей и следующей вершин ломаной.

Массив, определяющий положение вершин текущей ломаной в пирамиде. Модифицируется при просеивании элементов.

Массив вершин характеристических многоугольников Безье для всех п -1 отрезков гладкой кривой.

Массивы, определяющие положение вершин текущей ломаной на кривой (номер отрезка кривой, значение параметра).

Эти структуры используются в следующем алгоритме.

Алгоритм построения визуально гладкой ломаной

Шаг 1. Задание или определение максимального отклонения AL и максимального числа вершин N.

Шаг 2. Вычисление вершин характеристических многоугольников Безье для всех отрезков гладкой кривой в коридоре.

Шаг 3. Построение текущей (начальной) ломаной с n вершинами по узлам кривой в коридоре.

Шаг 4. Включение в текущую ломаную точек перегиба кривой.

Шаг 5. Расчет начальных отклонений и построение пирамиды.

Шаг 6. В цикле пока отклонение 1-го элемента в пирамиде больше AL и число вершин текущей ломаной меньше N необходимо:

выбрать вершину P из 1-го элемента пирамиды; определить для р предыдущую Pj и последующую Pk вершины в текущей ломаной;

найти на отрезках кривых PjPt и PiPk точки D1 и

D2, наиболее удаленные от хорд;

рассчитать новые отклонения для точек Pj, D1, P

(2 значения, соответствующие включению либо D1, либо D2), D2 и Pk ;

выбрать в качестве новой вершины точку D1 или D2 по критерию минимума максимального отклонения;

просеять в пирамиде точку P сверху вниз; добавить к текущей ломаной вершину D1 (D2) и просеять ее снизу вверх;

просеять точку Pj ( Pk ) вниз при уменьшении значения ее отклонения или вверх при увеличении.

Конец алгоритма.

Трудоемкость шагов 2-5 составляет O(n), а каждое просеивание на шаге 6 требует не более O(log N) операций. Таким образом, общая трудоемкость алгоритма не превышает O(N log N). Вычислительный эксперимент, проведенный на реальных данных, показал, что если n > 50, то для получения визуально гладкой ломаной достаточно задать N и 2n . На рис. 8 приводятся наборы исходных и визуально гладких расчетных изолиний одного участка местности.

Рис. 8. Исходные и визуально гладкие изолинии участка рельефа

ЗАКЛЮЧЕНИЕ

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

ЛИТЕРАТУРА

1. Костюк Ю.Л., Фукс А.Л. Построение и аппроксимация изолиний однозначной поверхности, заданной набором исходных точек // Геоинформатика: Теория и практика. Вып. 1. Томск: Изд-во Том. ун-та, 1998. С. 119-126.

2. Костюк Ю.Л., Фукс А.Л. Гладкая аппроксимация изолиний однозначной поверхности, заданной нерегулярным набором точек // Труды межд. научно-практ. конф. «Геоинформатика-2000». Томск: Изд-во Том. ун-та, 2000. С. 37-41.

3. Роджерс Д., Адамс Дж. Математические основы машинной графики: Пер. с англ. М.: Машиностроение, 1980. 240 с.

4. Костюк Ю.Л. Применение сплайнов для изображения линий в машинной графике // Автоматизация эксперимента и машинная графика. Томск: Изд-во Том. ун-та, 1977. С. 116-130.

5. Вирт Н. Алгоритмы и структуры данных: Пер. с англ. М.: Мир, 1989. 360 с.

6. Костюк Ю.Л., Фукс А.Л. Построение цифровой модели рельефа местности на основе структурных линий и высотных отметок // Вестник ТГУ. 2003. № 280. С. 286-289.

Статья представлена кафедрой теоретических основ информатики факультета информатики Томского государственного университета, поступила в научную редакцию 15 мая 2003 г.