CC BY
174
159
УДК 004.93:550.8
К. С. Алсынбаев
АЛГОРИТМЫ ОПРЕДЕЛЕНИЯ ТЕЛ ОБЪЕМНЫХ ОБЪЕКТОВ В ТРЕХМЕРНОМ НЕРЕГУЛЯРНОМ ОБЛАКЕ ТОЧЕК
Описана методика распознавания форм тел в трехмерном облаке точек. На первом шаге проводится кластеризация облака с критерием максимального расстояния между точками для идентификации распознаваемых тел. Далее производится переход к вокселному представлению. Критерием принадлежности воксела телу считается существование тетраэдра с вершинами в точках облака и ограниченными по размеру ребрами, для которого тестируемый воксел — внутренний. Для оптимизации разработан и используется быстрый алгоритм заполнение вокселных тетраэдров. Работа входит в состав программного обеспечения обработки данных микросейсмического мониторинга.
A technique for recognition of shapes of bodies in a 3D point cloud is described. At the first step cloud clustering is conducted using the criterion of maximum distance between the points to identify the bodies being determined.
© К. С. Алсынбаев, 2015
Вестник Балтийского федерального университета им. И. Канта. 2015. Вып. 10. С. 159 — 165.
Then proceeding to voxel presentation is being done. Existence of a tetrahedron with its points in cloud points and size-limited sides for which the voxel being tested is internal is treated as a criterion of belongment of voxel to a body. A fast algorithm for voxel tetrahedron filling is developed and used for optimization. The work is a part of software for micro-seismic monitoring data processing.
Ключевые слова: облако точек, кластеризация, воксел, микросейсмический мониторинг, алгоритмы распознавания, алгоритм Брезенхейма.
Key words: a point cloud, Clustering, Voxel, micro-seismic monitoring, Recognition algorithms, Bresenham's line algorithm.
1. Контекст разработки алгоритмов
В рамках работ по созданию программного обеспечения обработки данных микросейсмического мониторинга [1] появилась потребность выделения объемных объектов в множестве микросейсмических событий, которое по сути представляет собой трехмерное нерегулярное облако параметризованных точек. Эти алгоритмы были необходимы для разработки функций оценки характеристик объемных объектов горной природы, а также в качестве промежуточного шага для сепарации специфичных объектов, например разломов и трещиноватостей, представляющих образования, близкие к двумерным поверхностям.
Решение задачи распознавания объемных объектов предполагает:
• определение локализации и формы пространства, которое соответствует близким в метрическом смысле точкам, включая возможность невыпуклости этой формы и наличия отверстий и полостей;
• порождение идентификаторов объектов, соотнесение вычисленных форм и присвоенных им идентификаторов;
• возможность вычисления объемов объектов и подготовки отчетов в виде списков идентификаторов объектов и их объемов;
• возможность граневой визуализации выделенного объекта.
2. Основной алгоритм формирования фигур распознаваемых тел
Наиболее известное решение задачи определения форм тел, представленных множеством точек, — построение тетраэдризации Делоне (трехмерной триангуляции Делоне) с заданной границей [2]. Существуют также алгоритмы поиска границ, например [4]. Полученные таким образом тетраэдры, составляющие граневую аппроксимацию, позволяют вычислить объем тела, организовать интерактивную ЗО-визуа-лизацию. В контексте данной задачи также эффективен подход, когда строится триангуляция Делоне для выпуклой оболочки, затем удаляются тетраэдры, у которых радиусы описанных окружностей превышают заданный [3]. Алгоритмы имеют вычислительную эффективность O(NlogN). Одной из проблем тетраэдризации является чувствительность к точности вычислений.
Другой подход основан на представлении трехмерного пространства в виде дискретной сетки, элемент которой по аналогии с пикселем называется «воксел». Иначе говоря, производится своеобразная про-
странственная растеризация трехмерной области. С каждым вокселом связывается определенное значение или совокупность значений. Искомое трехмерное тело составляется в виде множества вокселов, имеющих необходимые значения.
В случае, когда исходное облако точек является «плотным» и регулярным, например полученным в результате томографического обследования, каждый воксел соответствует одной точке, изначально с точкой соотнесено определенное значение, вокселная матрица образуются автоматически, и выделение отдельного тела (соответствующего, например, какому либо человеческому органу) в виде вокселной модели трудности не представляет. Если же исходные данные представляют разреженное нерегулярное облако точек, как в случае микросей- 161
смических данных, необходимо алгоритмическое решение вопроса о -
принадлежности или непринадлежности каждого воксела телу.
В данной работе в некотором смысле объединены оба подхода. Область пространства с точками предварительно дискретизируется трехмерной сеткой, то есть представляется в виде трехмерной матрицы вокселов. Каждому вокселу соответствует область пространства кубической формы. Часть вокселов соответствует точкам облака (соответствующий кубик содержит одну или несколько точек облака) и помечается признаком «кластерный воксел» (выбор термина будет понятен ниже), остальным вокселам присваивается значение «0». Для каждого воксела со значением «0» определяется, находится ли он внутри объекта. Для определения понятия «внутри» устанавливается значение О максимального возможно диаметра внутренней полости тел. Внутренние вокселы тетраэдров, вершинами которых являются точки облака, и длины ребер не превышают значение О, считаются внутренними для тела тела (объекта).
Чтобы решить проблему предварительной локализации каждого тела и присвоения ему идентификатора, еще до этапа дискретизации проводится кластеризация, при которой точки считаются принадлежащими одному кластеру, если их для любой точки кластера существует точка этого же кластера, расстояние между которыми не превышает заданного. Предварительная кластеризация проводится с заданным максимальным расстоянием, равным О. После кластеризации каждый кластер получает идентификатор, равный порядковому номеру кластера, и считается, что точки кластера представляют отдельное тело, возможно, несвязное. Вокселы, соответствующие после дискретизации точкам кластера, будем называть кластерными вокселами, а в программной реализации будем помечать соответствующим значением. Кластерные вокселы изначально считаются принадлежащими телу.
Согласно формулированному выше воксел считается принадлежащим телу, если найдется тетраэдр, вершинами которого являются кластерные вокселы, содержащий данный воксел. Суть алгоритма заключается в переборе всех вокселов и проверке для каждого из них, существует ли допустимый вокселный тетраэдр с вершинами в кластерных вокселах, содержащий текущий воксел.
Вокселное пространство проверяется маской кубической формы размера МММ, где М = 2-т1(Ц/Р + 1) + 1, с шагом, равным Ы(О/Р+1)+1, то есть примерно половине стороны маски. Для кластерных вокселов внутри маски строится тетраэдризация Делоне [2], удаляются тетра-
162
эдры с ребрами, превышающими D, и все вокселы внутри полученных тетраэдров приобретают специальное значение «принадлежит телу и не кластерный». Возможны случаи, когда воксел содержится внутри вырожденного тетраэдра, внутри треугольника или отрезка. В этом случае заполняется треугольник в трехмерном пространстве или проводится отрезок.
В качестве оптимизации используется следующая операция. Если тестирование определило, что воксел содержится внутри кластерного тетраэдра, то все внутренние точки этого тетраэдра помечаются как «принадлежащие телу, но не кластерные», с использованием быстрого алгоритма заполнения тетраэдра. Алгоритмы заполнения тетраэдра, треугольника и трехмерного отрезка используют принципы растровой графики, восходящие к алгоритму Брезенхейма, и при вычислениях в основных циклах не используется плавающая арифметика, операции умножения и деления. О реализации алгоритмов изложено ниже.
Полученная совокупность вокселов тела — приближенная относительно решения в вещественных значениях с точностью, равной размеру ребра куба, соответствующего вокселу. Поскольку каждый кластер обрабатывается отдельно, то порождаемая при этом вокселная модель является отдельным объектом.
В работе реализованы средства граневой визуализации, простого изображения кубиками (рис. 1), или полученную алгоритмом marching-cubes (реализована для данной задачи М. В. Козловым).
Рис. 1. Изображения оболочек трех распознанных тел. Показаны в вокселной модели гранями кубиков
Приведем оценку вычислительной эффективности алгоритма. Пусть N — число точек облака, D — максимальный диаметр полости тела, Р — размер стороны воксела, K — количество кластеров, R — среднее число точек облака в кластере, V — среднее количество вокселов в кластерах, L — среднее количество точек облака в области маски, C — средний линейный размер кластера, L — среднее количество точек облака (или кластерных вокселов) в пределах маски.
Пусть значок ~ означает «приблизительно пропорционально». В кластере необходимо ~ C3/ D3 перемещений маски. Для перебора всех тетраэдров в пределах маски нужно ~ L logL операций для заполнения тетраэдров в пределе маски ~ М3 операций. Таким образом, получаем оценку мощности вычислений ~ K(C3/ D3)(L logL + М3) (не будем использовать префикс «О», чтобы избежать некорректности его использования). Принимая R ~ N/K, М ~ D/P, L ~ RD3/C3 ~ ND3/(KC3), получаем оценку мощности вычислений
N(logN + 3logD - logK + 3logC) + KC3/P3 .
Таким образом, при постоянных C, P и K мощность алгоритма соответствует O(NlogN), причем такой рост обуславливается только ростом числа точек в окрестности тестируемой точки. Заметим, что при росте C уменьшается K. Время вычислений растет с увеличением C и резко растет с уменьшением P.
При значениях K = 300 и P = 0,5 м плотного облака в 480000 точек данных микросейсмического мониторинга длилась 2 часа 39 минут (без этапа кластеризации) на компьютере с процессором Intel® Xeon CPU W3565 3,19 GHz 12 ГБ.
Реализация позволяет производить поиск тел только для указанных кластеров, что обеспечивает возможность подбора значений P при повторных запусках и завершения алгоритма для всех кластеров в приемлемое время. Алгоритм легко распараллеливается, в том числе для кластерных суперкомпьютеров.
3. Трехмерные аналоги алгоритмов Брезенхейма
Как известно, алгоритмы Брезенхейма и их производные предназначены для порождения графических объектов в растровых множествах. Особенность их в том, что они не содержат операций деления и умножения, а в основном цикле вычисления производятся исключительно на основе операций целочисленных сложения, вычитания, сравнения, благодаря чему достигается высокая скорость выполнения графических функций.
Вокселная матрица — трехмерный аналог растровых данных, и для них возможно создание аналогичных функций.
В данной статье были разработаны алгоритмы для генерации трехмерных вокселных векторов, трехмерных многоугольников, тетраэдров.
В основе генерации заполненных тетраэдров — быстрый алгоритм порождения плоских многоугольников. Для каждого горизонтального слоя тетраэдра генерируется отдельный выпуклый плоский растровый многоугольник, соответствующий сечению — треугольнику или четырехугольнику, вершины которых вычисляются трехмерными алгоритмами Брезенхейма проведения линий. На рисунке 2, а показаны два горизонтальных слоя, создаваемых быстрыми алгоритмами, а ребра тетраэдра проведены (для демонстрации) трехмерными вокселными векторами. На рисунке 2б — соответствующий заполненный вокселный тетраэдр.
164
Рис. 2: а — элементы алгоритма создания заполненного тетраэдра. б — заполненный вокселный тетраэдр
4. Алгоритм кластеризации облака точек
Оптимизация алгоритма осуществляется на исключении из сравнения точек, находящихся на достаточно большом расстоянии от текущей проверяемой. Принцип схожего алгоритма для двумерного случая, называемый «кластеризация точек на регулярной сети», приведен, например в [5]. В рассматриваемой задаче он модифицирован для повышения вычислительной эффективности.
Алгоритм заключатся в следующем. Пусть задано максимальное расстояние между точками одного кластера Dmax и минимальный размер каждого кластера. Вычисляется максимальная дистанция для сравнения по каждой из координат, равная GRIDV = Dmax*2. Каждая точка (x, y, z) снабжается полем, задающим номер кластера типа «целое» и полем для индексации типа «целое». Создается массив SizeOfCluster, содержащий размеры каждого кластера. На предварительном этапе вычисляется покоординатные индексы iX, iY, iZ каждой точки по формуле iX = int((x-min(x))/GRIDV) и общий индекс-свертка на основе позиционного способа ind(x) = iZ*max(iX+1)*max(iY+1) + iY* max(iX+1) + iX. Очевидно, зная max(iX+1) и max(iY+1), из ind(x) легко восстановить iX, iY, iZ. Значения номеров кластеров для всех точек «сбрасываются» присвоением условного значения «-1».
В основном блоке производится обход всех точек во внешнем цикле, фиксируется очередная точка, если она не принадлежит никакому кластеру, порождается новый кластер и номер для нее, если точка уже в кластере, ее номер запоминается как текущий; производится обход всех точек во внутреннем цикле 2-го уровня, назовем эти точки тестируемыми. Проверяется, является ли тестируемая точка близкой к текущей на основе проверки условия abs(iX-iXtest)<=1 and abs(iY-iYtest)<=1 and abs(iZ-iZtest)<=1. Далее, если выполняется условие, что текущая и тестируемая точка расположены по метрике Евклида ближе, чем Dmax, то тестируемая точка либо присоединяется (если она не принадлежит ни к какому кластеру) к кластеру текущей точки, либо два кластера, текущей и тестируемой точки, объединяются. Для объединения кластеров необходим еще один внутренний цикл (3-го уровня) по просмотру номеров кластеров точек. В заключительном блоке проводится отсеивание мелких кластеров, оптимизация нумерации кластеров и массива их размерности.
Если использовать наиболее простой алгоритм, мы имеем три вложенных цикла по массиву точек (событий), вычислительная эффективность алгоритма равна О(^), где N — количество точек (событий). Возможно улучшение эффективности алгоритма до при организации кластеров в списочные структуры. Если провести упорядочивание по индексу-свертке i(x) и индексацию точек массива точек, то можно улучшить эффективность до О(N log2 (N)). Однако каждая оптимизация по времени требует дополнительной памяти: для организации списочной структуры — 4 целочисленных массива размером в количество обрабатываемых событий, упорядочивание по индексу-свертке — еше двух таких же массивов.
Эксперименты показали, что время выполнения кластеризации на тестовом массиве 7000 точек около одной секунды на компьютере с процессором Intel® Xeon CPU W3565 3,19 GHz 12 ГБ. Время кластеризации на тестовом массиве 480 000 точек составило около 3 часов.
После оптимизации алгоритма - хранение кластеров в списках и объединение их путем слияния списков (при этом при каждом объединении просмотр точек присоединяемого кластера и переназначение номера кластера) — время поиска кластеров для 480 000 точек составило 16 минут. Прогноз поиска кластеров для миллиона точек длится, таким образом, около 2 часов, что является приемлемым.
В нашей задаче целесообразен выбор Dmax равным D или немного меньшим D.
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (код проекта 14-07-00699).
Список литературы
1. Алсынбаев К. С., Козлов А. В. Средства распознавания и визуализации разломов и зон техногенной трешиноватости на основе обработки данных микросейсмического мониторинга // Вестник Балтийского федерального университета им. И. Канта. 2014. Вып. 4. С. 127 — 134.
2. Боровиков С. Н., Иванов И. Э., Крюков И. А. Построение тетраэдризации Делоне с ограничениями для тел с криволинейными границами // Журнал вычислительной математики и математической физики. 2005. Т. 45, № 8. С. 1407—1423.
3. Суков С. А. Методы генерации тетраэдральных сеток и их программные реализации // Препринты ИПМ им. М. В. Келдыша. 2015. № 23.
4. Baidurja Ray, Avi Lin, Jianfu Ma. Unconventional micro-seismicity based enhanced 3D SRV estimator using advanced parameter-free concave methodology // SEG Technical Program Expanded Abstracts. 2014. P. 2304-2308.
5. Кластеризация точек на регулярной сети. URL: http://habrahabr.ru/post/ 138185/ (дата обрашения: 19.08.2015).
Об авторе
Камил Салихович Алсынбаев — канд. техн. наук, доц., Балтийский федеральный университет им. И. Канта, Калининград.
E-mail: [email protected]
About author
Dr Kamil Alsynbaev — Ass. Prof., head of the laboratory, I. Kant Baltic Federal University, Kaliningrad.
E-mail: [email protected]
