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

Ю.Л. Костюк, К.Г. Гульбин, С.В. Пешехонов

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

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

Лазерная локация - прикладная наука, изучающая вопросы использования лазерных локаторов (лидаров) для проведения воздушной топографо-геодезической съемки. При этом получаются массивы (облака) координат точек отражения объемом в многие миллионы при плотности до нескольких точек на 1 м2 [1]. Однако в настоящее время не известны такие методы обработки координат точек, которые могли бы автоматически строить полноценные модели объектов на поверхности Земли.

В статье предлагается вначале на множестве точек построить поверхностную триангуляцию. При этом нельзя непосредственно использовать плоскую триангуляцию в проекции на плоскость ХОУ, так как при наличии на поверхности зданий, сооружений, деревьев и других объектов после добавления к точкам координаты 1 треугольники, плоскости которых близки к вертикальной, будут неверно аппроксимировать эти объекты. Поэтому после построения плоской триангуляции предлагается выполнять локальное перестроение пространственных треугольников, уменьшающее длины их ребер. После этого можно использовать алгоритм, выделяющий сами объекты на поверхности Земли. Особенностью предлагаемого алгоритма является то, что в нем на первом этапе выделяются фрагменты плоскостей, аппроксимирующих группы смежных треугольников, а на втором этапе по смежным выделенным фрагментам плоскостей строятся более сложные аппроксимации в виде фрагментов сфер, цилиндров и др. Такой подход облегчает для человека-эксперта принятие решения о типе выделенного объекта.

Построение поверхностной триангуляции

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

1) те же облака точек, но разделенные по своей морфологической принадлежности (земля, растительность, поверхности водоемов, кровли зданий, провода ЛЭП и т.п.);

2) цифровые модели рельефа;

3) ортофотомозаика;

4) распознанные топографические объекты и их модели в трехмерном представлении.

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

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

Неоднозначность вертикальной координаты точки отражения от ее плоских координат может быть вызвана:

- наличием вертикальной плоскости;

- многоярусностью, когда под некоторым монолитным объектом находится пустота.

Такие топологии не присущи рельефу поверхности Земли и абсолютному большинству природных объектов, в этом случае можно использовать однозначную функцию, заданную на плоскости ХОУ. Что касается инженерных сооружений, большинство из них содержат вертикальные поверхности (такие, как стены), некоторые обладают многоярусностью (например, мост). В природе встречаются области, в которых наклон поверхности велик, но которая может быть однозначно задана на плоскости ХОУ (горы, овраги и т.п.). Таким образом, абсолютное большинство областей, где плоская триангуляция не соответствует пространственной, представляют собой группы треугольников, построен-

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

Рис. 1. Группа почти вертикальных треугольников в плоской триангуляции на фоне треугольников, соответствующих пологим участкам поверхности

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

Количество точек на единицу поверхности в первичных данных может быть различной, ее нижний предел -порядка нескольких точек на один квадратный метр [1]. Расстояние между соседними точками в плоской триангуляции прямо пропорционально плотности облака и обратно пропорционально косинусу угла наклона поверхности между этими точками. Это расстояние для областей малого уклона (менее 45°) не может быть меньше 25-30 см и, напротив, обычно мало (порядка относительной погрешности лазерно-локационного измерения) для областей со значительным уклоном. Таким образом, необходимым критерием вертикальности треугольника являются малые длины его ребер.

Две плоскости

относительно плоскости ХОУ, выражаемый соотношением (4), где а - меньший из углов ф1 и ф2, а у - некоторый параметр

а > 900 - у . (4)

Плоскость ХОУ описывается уравнением г = 0, а уравнение плоскости, в которой лежит треугольник, вычисляется по формуле [2]:

х - х0 у - у0 г - г0

х1 - х0 У1 - у0 г1 - г0 — 0,

х2 - х0 У2 - у0 г2 - г0

где х,, у, г{ - координаты 1-й точки треугольника.

Процедура проверки треугольника на вертикальность такова: рассматриваем длины ребер треугольника в ХОУ и, если они не превышают значения некоторого параметра е, рассматриваем вычисляемые по формуле

(3) углы наклона между плоскостью г = 0 и плоскостью треугольника, вычисляемой по формуле (5). В случае, если меньший из этих углов удовлетворяет критерию

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

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

Шаг 2. Пока есть непросмотренные треугольники, цикл:

2.1. Очередной непросмотренный треугольник помечается как просмотренный и проверяется на вертикальность.

2.2. Если треугольник почти вертикальный, то:

1) запускается процедура выделения группы смежных почти вертикальных треугольников;

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

Критерием того, что два треугольника лежат примерно в одной плоскости, является малый угол между их плоскостями

А1 х + В-[ у + С-^ г + Dl — 0 х + В2 у + С 2 г + Е2 — 0

(1)

(2)

образуют четыре двухгранных угла, попарно равных и вычисляемых по формуле [2]:

008 ф1 2 — ±

Л А + В В + С С

^Л2 + в2 + с2 ^Л2 + в22 + с22 '

(3)

р <е,

(6)

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

где в - меньший из двух углов ф! и ф2 , полученных по формуле (3), в качестве параметров плоскостей в которой подставлены параметры плоскостей треугольников, вычисленные по формуле (5); е - некоторый произвольный параметр, соответствующий максимальному углу между плоскостями, при котором треугольники все еще считаются лежащими примерно в одной плоскости.

Координаты треугольников каждой подгруппы преобразуем в координаты плоскости эталонного треугольника. При замене системы координат ХУ2 новой

и

системой ХнУн1н с тем же началом старые координаты точки (х, у, г) выражаются через новые (хн, ун, гн) по формуле

где

А =

(х, y, г) = (хн, ун, гн)А ,

( СО$,(їн , Г) С08(]н, Г) соь(кн, Г)Л

С08(Гн , ]) С0Є(]н , ]) С0Б(к н , ])

соб(/н , к) С08(Ун, к) С08(кн, к)

(7)

і,],к - единичные векторы старой системы координат, а ін, ]н, кн - единичные векторы новой системы координат [2]. В свою очередь, новые координаты точки выражаются через старые:

(хн , Ун , гн ) = У, г)А '.

(8)

Косинус угла между двумя векторами а1 — (х1, у1, г1) и а2 — (г2, у2, г2), лежащими в одной плоскости, находится через их скалярное произведение:

С0Є ф =

Х1 Х2 + Уі У 2 + г1 г2

Vх!2 + у2 + г2 -\/х2 + у2 + г2

(9)

Осями базовой системы координат являются векторы (1,0,0), (0,1,0) и (0,0,1), а в качестве осей новой системы координат примем вектор, образованный одной из сторон треугольника, вектор нормали к поверхности и вектор, дополняющий их до прямоугольной системы координат. Этот третий вектор является вектором нормали по отношению к плоскости, образованной первыми двумя. Пусть заданы два вектора а1 = (х1, у1, г1) и а2 = (г2, у2, г2), лежащие в одной плоскости. Они перпендикулярны вектору нормали „ = (хп,уп, г„), т.е. соответствующие скалярные произведения равны нулю:

Х1 Х„ + у у„ + г1 г„ = 0,

Х2 Х„ + у2 у„ + г2 г„ = 0.

(10)

Умножим верхнее уравнение на х2 и вычтем из него нижнее, умноженное на х1. Получим

(у х2 - у2 Х1 ) у„ = ( г2 Х1 - г1 Х2 ) г„ Х2 Х„ + у2 у„ + г2 г„ = °.

(11)

Приняв гп = у1 х2 - у2 х1, получим уп = г2 х1 - г1 х2. Подставив теперь гп и уп во второе из уравнений системы (11), найдем хп:

у 2 (г 2 Х1 - г1 Х2) + г2(у1 Х2 - у 2 Х1)

. (12)

Воспользовавшись вышеприведенными выкладками, можно вычислить вектор нормали к плоскости тре-

угольника, взяв в качестве векторов а1 и а2 его стороны, а затем вычислить третий осевой вектор новой системы координат.

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

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

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

Выделение фрагментов плоскостей, аппроксимирующих группы смежных треугольников пространственной триангуляции

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

Алгоритм выделения фрагментов плоскостей (граней).

Шаг 1. В цикле все треугольники поверхностной триангуляции помечаются как непросмотренные.

Шаг 2. Цикл, пока есть непросмотренные треугольники.

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

Х

2

2.2. По выбранному множеству треугольников (точнее, их вершин, т.е. точек) строится аппроксимирующая плоскость, которая задаётся уравнением Ах + Ву + С2 + Б = 0. Аппроксимация производится в смысле метода наименьших квадратов [5] (сумма квадратов расстояний от точек до полученной плоскости должна быть минимальной).

2.3. Аналогично тому, как это делалось на шаге 2.1, начиная от первого треугольника, в ширину выбираются смежные треугольники. Выбор треугольника происходит в том случае, если все расстояния от его трех вершин до аппроксимирующей плоскости не превосходят некоторой величины 5 (параметра в программе). Если треугольник выбран, то рассматриваются его соседи и т.д. Все выбранные треугольники помечаются как просмотренные.

2.4. Если треугольников, которые были выбраны на шаге 2.3, достаточно для выделения грани (третий параметр в программе), то такой фрагмент выделяется, а треугольники, его составляющие, помечаются принадлежащими этому объекту.

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

1

кД/ 2\ 0 /\ / /3 V/

\ ^ /\ 7 > / \/ ]\1 \

\ 6/\ 0 /\ /

Д/ 2\ /3 \/

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

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

Выделение более сложных фрагментов поверхностей

После выделения плоскостных объектов имеются следующие данные:

1) начальный массив вершин (точек);

2) массив треугольников, представляющий поверхностную триангуляцию, при этом каждый треугольник содержит пометку - индекс плоскости, в которую он входит, или пометку о том, что данный треугольник не принадлежит ни одной из выделенных плоскостей;

3) массив плоских граней, в котором каждая плоскость содержит индекс первого треугольника, принадлежащего данному объекту.

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

Шаг 1. Строится граф топологической связности граней. Две грани (вершины графа) имеют ребро, если одной грани принадлежит треугольник Т1, смежный треугольнику Т2, принадлежащему другой грани.

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

Шаг 3. Для выбранного набора граней строится в первом приближении фигура, которую они ограничивают.

Шаг 4. По точкам, составляющим выбранные на шаге 2 грани, итеративно, отбрасывая точки, находящиеся достаточно далеко от текущего приближения фигуры, строится новое приближение фигуры по методу наименьших квадратов [5].

Шаг 5. Полученное приближение фигуры используется для определения принадлежности других граней данной фигуре. Для этого на графе топологической связности производится поиск в ширину. Если какая-либо грань оказывается принадлежащей данной фигуре, то она прикрепляется к данной фигуре и исключается из дальнейшего рассмотрения (для данного типа фигур). Считается, что некоторая грань принадлежит данной фигуре, если хотя бы половина её точек лежит на её поверхности (с некоторым допуском).

Шаг 6. Если на шаге 5 было выделено достаточно граней (параметром задаётся минимальное общее количество вершин в гранях, принадлежащих этой фигуре), то данная фигура выделяется, иначе все изменения, сделанные на шаге 5, отменяются.

Шаг 7. Если есть непросмотренные грани, то повторение действий с шага 2.

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

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

Любую сферу можно точно определить через четыре плоскости, которые её касаются. Для того чтобы определить сферу, решается система линейных уравнений, каждое из которых представляет собой расстояние от плоскости до центра сферы. Расстояние от некоторой плоскости, заданной уравнением Ах + Ву + С2 + Б = 0, до точки Р0 = (х0, у0,10) может быть вычислено по формуле

d —

Ax0 + Вуд + Czq + D

Ц2

+ B2 + C2

(13)

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

' A Хд + В Уд + C Zg + D | .

(14)

Для того чтобы узнать, какой знак будет иметь выражение после раскрытия модуля, примем во внимание, что центр сферы будет лежать по ту же сторону от плоскости, по которую лежат остальные грани, точнее, подавляющее множество ее точек. Обозначим знак перед радиусом г после раскрытия модуля в формуле (14) символом Е' и получим уравнение

А'Хд + В 'Уд + C 'Zg — E Г — —D'.

(15)

f a; b; C\ —e; 1 f x л л0 f— D’l Л

А2 B2 C2 —E2 Уо — D2

a; B'3 C ^3 — E3 Z0 — D3

1 A4 B'4 C 4 —E4 1гJ I— D4 J

Решение системы (16) даёт точку Р0(х0,у0,і0) -центр сферы, которую касаются четыре грани, а также ее радиус.

На шаге 4 из четырёх плоскостей, образующих кусок поверхности сферы, берутся все точки, и, используя приближение, полученное в результате решения системы (16) в качестве начального для метода наименьших квадратов, получается более точное решение. Метод наименьших квадратов для выделения таких фигур, как сфера, тор, конус и цилиндр, описан в [5].

На рис. 2 представлен пример работы алгоритма выделения фрагмента сферы. Различными цветами закрашены грани, точки которых образовали сферу.

Аналогично получаются оставшиеся три уравнения для остальных граней, в результате получается система

(16)

Рис. 3. Фрагмент сферы с выделенными на ней гранями

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

ЛИТЕРАТУРА

1. Данилин И.М., Медведев Е.М., Мельников С.Р. Лазерная локация земли и леса: Учебное пособие. Красноярск: Ин-т леса им. В.Н. Сукачева

СО РАН, 2005.

2. Выгодский М.Я. Справочник по высшей математике. М.: Наука, 1973.

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

во Том. ун-та, 1998. С. 61-66.

4. Костюк Ю.Л., Фукс А.Л. Приближенное вычисление оптимальной триангуляции // Межд. конф. «Дискретный анализ и исследование опера-

ций»: Матер. конф. Новосибирск: Изд-во ин-та математики, 2000. С. 152.

5. Shakarji C.M. Least-Squares Fitting Algorithms of the NIST Algorithm Testing System // Journal of Research of the National Institute of Standards

and Technology. 1998. Vol. 103, № 6. P. 633-635 [Электронный ресурс] - Режим доступа: http://nvl.nist.gov/pub/nistpubs/jres/103/6/j36sha.pdf, свободный.

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