От облака точек к поверхности (From point cloud to surface)

Автор: | 13.09.2020

Введение
Классификация алгоритмов перехода от облака точек к поверхности
Этапы перехода от облака точек к поверхности
Пакеты для 3D-визуализации сетей. Модуль PyVista
Примеры, демонстрирующие возможности PyVista
Форматы файлов облака точек
Визуализация данных из PLY и PCD файлов
Полезные ссылки

Введение

Облако точек (англ. point cloud) — набор вершин в трёхмерной системе координат. Эти вершины, как правило, определяются координатами X, Y и Z и, как правило, предназначены для представления внешней поверхности объекта.

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

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

Хотя облака точек могут быть непосредственно визуализированы и проверены, они, как правило, не используются напрямую в большинстве 3D-приложений, и поэтому, как правило, конвертируются в полигональную сетку, модели с NURBS—поверхностями или CAD-модели при помощи процесса, известного как «реконструкция поверхности» (англ. surface reconstruction).

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

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

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

Восстановление поверхности  — сложная задача:

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

Классификация алгоритмов перехода от облака точек к поверхности

Общая классификация алгоритмов проводится в соответствии с качеством (типом) исходных данных:

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

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

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

Другая классификация основана на способе представления поверхности:

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

Поверхность называется N-гладкой относительно заданной параметризации, если каждая из координатных функций x(u,v), y(u,v), z(u,v) имеют непрерывную производную до порядка N включительно.

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

Классификация по использованию данных:

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

Классификация алгоритмов по наличию дополнительной информации:

  • известен топологический тип поверхности (например, плоскость, цилиндр или сфера)

  • информация о структуре (чайник с носиком, ручкой и крышкой) и ориентации

Этапы перехода от облака точек к поверхности

Преобразование облака точек в поверхность состоит из четырех шагов:

Предварительная обработка:

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

Определение глобальной топологии поверхности объекта:

  • установление отношений между соседними частями поверхности
  • сортировка и рассмотрения возможных «ограничений» (например, определение структурных линий)
  • выделения особых свойств (например, кромок).

Создание многоугольной поверхности:

  • Это основная операция, которая преобразует набор точек в согласованную полигональную модель (сетку) с вершинами, ребрами и гранями.
  • Треугольная (или тетраэдрическая) сетка создается при требованиях определенного качества, например — ограничение на размер элемента сетки, отсутствие пересечения структурных линий и т. д.;
  • Хорошо известный метод строительства сетки — Делоне триангуляция (DT), который одновременно оптимизирует упомянутые выше показатели качества.На рисунке приведены диаграмма Вороного (слева) и триангуляция Делоне (справа) того же набора точек. Каждая область диаграммы Вороного состоит из части плоскости, ближайшей к этому узлу. Соединение узлов ячеек Вороного, имеющих общие границы образуют треугольники Делоне.  Критерий Делоне гарантирует, что никакая вершина не лежит внутри окружностей,  описанных вокруг любого из треугольников в сети.
  • Триангуляция выполняется в 2D или в 3D, в зависимости от геометрии входных данных. Триангуляция в 3D называется тетраэдризацией или
    тетраэдризация. Тетраэдризация — это разбиение входной домен в набор тетраэдров, которые встречаются только в общие грани (вершины, ребра или треугольники). Тетраэдризация результаты намного сложнее, чем двумерная триангуляция. Типы входных областей могут быть простыми многогранниками (сфера), непростыми многогранник (тор) или облаками точек.

Постобработка — операции редактирования  для уточнения и совершенствования многоугольной поверхности:

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

Пакеты для 3D-визуализации сетей. Модуль PyVista

В Python есть несколько вариантов 3D-визуализации сетей, среди которых несколько известных проектов:

  • Matplotlib (Hunter, 2007),
  • Mayavi (Ramachandran & Varoquaux, 2011),
  • the yt Project (Turk et al., 2010),
  • the Visualization Toolkit (VTK) (Schroeder, Lorensen, & Martin, 2006).

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

PyVista (ранее vtki) — это вспомогательный модуль для VTK, который использует другой подход к взаимодействию с VTK — через NumPy и прямой доступ к массиву. В отличие от VTK  пакет PyVista  предоставляет краткий, хорошо документированный интерфейс, демонстрирующий мощную систему визуализации VTK. Этот модуль можно использовать для научного построения презентаций и исследовательских работ, а также в качестве вспомогательного модуля для других модулей Python, имеющих отношение к сеточным структурам.

PyVista — это:

  • «VTK для людей» — это высокоуровневый API для Visualization Toolkit (VTK)
  • Сеточные структуры данных и методы фильтрации для пространственных наборов данных
  • Простое построение трехмерных графиков для больших и сложных геометрических объектов

Примеры, демонстрирующие возможности PyVista

Ниже приводятся несколько  быстрых примеров, демонстрирующих возможности PyVista.

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

import numpy as np
import pyvista as pv
# mesh points
vertices = np.array([[0, 0, 0],
                     [1, 0, 0],
                     [1, 1, 0],
                     [0, 1, 0],
                     [0.5, 0.5, -1]])
# mesh faces
faces = np.hstack([[4, 0, 1, 2, 3], # square
                   [3, 0, 1, 4], # triangle
                   [3, 1, 2, 4]]) # triangle
surf = pv.PolyData(vertices, faces)
# plot each face with a different color
surf.plot(scalars=np.arange(3), cpos=[-1, 1, 0.5])

Пример 2.

import pyvista as pv
import numpy as np
pts = np.random.rand(512*3).reshape(-1,3)
# Make vtkPolyData of the points array
point_cloud = pv.PolyData(pts)
point_cloud.plot(render_points_as_spheres=True, point_size=10)
# runs the delaunay 3D algorithm
mesh = point_cloud.delaunay_3d(alpha=0.25)
# Make a wireframe with some data
wires = mesh.compute_cell_sizes(length=False, area=False, volume=True).wireframe()
# Plot it
wires.plot(line_width=2)

Пример 3. Создание поверхности из набора точек с помощью триангуляции Делоне — с отверстием.

# sphinx_gallery_thumbnail_number = 2
import pyvista as pv
import numpy as np
# First, create some points for the surface.
# Define a simple Gaussian surface
n = 20
x = np.linspace(-200, 200, num=n) + np.random.uniform(-5, 5, size=n)
y = np.linspace(-200, 200, num=n) + np.random.uniform(-5, 5, size=n)
xx, yy = np.meshgrid(x, y)
A, b = 100, 100
zz = A * np.exp(-0.5 * ((xx / b) ** 2.0 + (yy / b) ** 2.0))

# Get the points as a 2D NumPy array (N by 3)
points = np.c_[xx.reshape(-1), yy.reshape(-1), zz.reshape(-1)]
points[0:5, :]

# simply pass the numpy points to the PolyData constructor
cloud = pv.PolyData(points)
cloud.plot(point_size=15)

# Now that we have a PyVista data structure of the points, we can perform a
# triangulation to turn those boring discrete points into a connected surface.
surf = cloud.delaunay_2d()
surf.plot(show_edges=True)

# Masked Triangulations
x = np.arange(10, dtype=float)
xx, yy, zz = np.meshgrid(x, x, [0])
points = np.column_stack((xx.ravel(order="F"),
                          yy.ravel(order="F"),
                          zz.ravel(order="F")))
# Perturb the points
points[:, 0] += np.random.rand(len(points)) * 0.3
points[:, 1] += np.random.rand(len(points)) * 0.3
# Create the point cloud mesh to triangulate from the coordinates
cloud = pv.PolyData(points)
cloud

# Run the triangulation on these points
surf = cloud.delaunay_2d()
surf.plot(cpos="xy", show_edges=True)

# Note that some of the outer edges are unconstrained and the triangulation
# added unwanted triangles. We cn mitigate that with the ``alpha`` parameter.
surf = cloud.delaunay_2d(alpha=1.0)
surf.plot(cpos="xy", show_edges=True)

# We could also add a polygon to ignore during the triangulation via the
# ``edge_source`` parameter.

# Define a polygonal hole with a clockwise polygon
ids = [22, 23, 24, 25, 35, 45, 44, 43, 42, 32]

# Create a polydata to store the boundary
polygon = pv.PolyData()
# Make sure it has the same points as the mesh being triangulated
polygon.points = points
# But only has faces in regions to ignore
polygon.faces = np.array([len(ids),] + ids)

surf = cloud.delaunay_2d(alpha=1.0, edge_source=polygon)

p = pv.Plotter()
p.add_mesh(surf, show_edges=True)
p.add_mesh(polygon, color="red", opacity=0.5)
p.show(cpos="xy")

Пример 4. Использование фильтра связности, чтобы удалить зашумленные изоповерхности

Изображенная на рисунке слева сетка очень зашумлена. Мы можем извлечь самые большие связанные isosurface в этой сетке с помощью func: `pyvista.DataSetFilters.connectivity` или с помощью фильтра: func: pyvista.DataSetFilters.extract_largest` (обе эквивалентны).

import pyvista as pv
from pyvista import examples
#######################################################################
# Load a dataset that has noisy isosurfaces
mesh = examples.download_pine_roots()
cpos = [(40.6018, -280.533, 47.0172),
        (40.6018, 37.2813, 50.1953),
        (0.0, 0.0, 1.0)]
# Plot the raw data
p = pv.Plotter()
p.add_mesh(mesh, color='#965434')
p.add_mesh(mesh.outline())
p.show(cpos=cpos)

#######################################################################
# Grab the largest connected volume present
largest = mesh.connectivity(largest=True)
# or: largest = mesh.extract_largest()
p = pv.Plotter()
p.add_mesh(largest, color='#965434')
p.add_mesh(mesh.outline())
p.camera_position = cpos
p.show()

Форматы файлов облака точек

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

Используемые типы файлов поддерживают как ASCII код,  передающий данные в виде строк текста, так и двоичные форматы. Каждая строка текста в файле ASCII представляет запись в виде пространственных координат (x, y, z). Дополнительная информация, такая как цвет и интенсивность, может быть включена в некоторые форматы. Файлы ASCII можно открывать в текстовых редакторах. Однако использование текста для передачи данных связано с расходами. Форматы двоичных файлов более компактны. Также требуется меньше времени для обработки и визуализации двоичных файлов, поскольку их можно проиндексировать, что позволяет читать их по частям, а не последовательно.

PCD (Point Cloud Data) — формат необработанных облаков точек

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

STL (от англ. stereolithography) — формат файла, широко используемый для хранения трёхмерных моделей объектов c целью использования в аддитивных технологиях. Информация об объекте хранится как список треугольных граней, которые описывают его поверхность, и их нормалей.

OBJ — это простой формат данных, который содержит только 3D геометрию, а именно, позицию каждой вершины, связь координат текстуры с вершиной, нормаль для каждой вершины, а также параметры, которые создают полигоны.

 

Визуализация данных из PLY и PCD файлов

По ссылке   PLY Files an ASCII Polygon Format   можно получить множество файлов PLY формата. Один из них (airplane.ply) использую в программе на Python для визуализации данных.

Создаю в редакторе «Блокнот» файл  airplane.ply   с текстовкой, скопированной со страницы по указанной выше ссылке.

import numpy as np
import open3d as o3d

def main():
 cloud = o3d.io.read_point_cloud("airplane.ply") # Read the point cloud
 print(cloud)
 o3d.visualization.draw_geometries([cloud]) # Visualize the point cloud

if __name__ == "__main__":
 main()

Программа чтения PCL может правильно прочитать и PCD файл. По ссылке   получаю файл формата .pcd. (chair.pcd) и читаю его в той же программе.

cloud = o3d.io.read_point_cloud("chair.pcd")

 

 

Полезные ссылки:

 

Автор: Николай Свирневский