Всем привет! Сегодня мы хотели бы поделиться с вами нашим опытом анализа лидарных облаков. В заметке расскажем:

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

Заметка от партнера IT-центра МАИ и организатора магистерской программы “VR/AR & AI” — компании PHYGITALISM.

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

Инструменты для анализа лидарных данных

Итак, предположим что вы начинающий исследователь в области 3D ML и вам потребовалось работать с лидарными данными. Какие инструменты лучше для этого использовать? В чем особенность различных форматов данных? Какие подходы можно применить? В чем сложность обработки и анализа таких данных? Эти и многие другие вопросы встают перед исследователем. Попробуем начать с источников для изучения, инструментов и библиотек.

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

Из перечисленных в этом репозитории инструментов стоит отметить два основных, которыми пользуются многие исследователи и разработчики в области пространственного компьютерного зрения и 3D графики:

MeshLab (Open source, portable, and extensible system for the processing and editing of unstructured 3D triangular meshes).
Это быстрый, простой в освоении и использовании инструмент, поддерживающий большинство форматов 3D данных для облаков точек и полигонального меша. На официальном сайте можно найти много туториалов. В MeshLab есть инструменты для фильтрации, сегментации, сглаживания и многого другого. На хабре про MeshLab писали тут. MeshLab хорошо подходит для того, чтобы открыть не очень объемный скан или модель, посмотреть на нее, убрать шум или сегментировать на несколько частей, залатать дыры в меше и подобные операции. Для того чтобы работать с большими облаками, которые обычно записывают в .laz и .las формат, лучше подойдет следующий инструмент.
Cloud Compare (3D point cloud and mesh processing software Open Source Project)
Данный инструмент предназначен для визуализации и обработки больших облаков точек, он поддерживает большинство форматов, используемых для лидарных данных. Есть инструменты для редактирования облака точек, имеются модули из PCL, есть возможность обработки на GPU, для больших облаков можно генерировать октодеревья, для быстрой обработки.На YouTube есть много туториалов по этому инструменту, а на официальном сайте есть облака в .las формате для этих туториалов.
ROS and Webviz
В области беспилотного транспорта часто встречается ROS (The Robot Operating System) - фреймворк для написания программ в области робототехники. Для записи лидарных данных в ROS используют .bag формат, имеющий много особенностей и подводных камней. ROS и BAG файлы стоят отдельной большой заметки на Хабре (туториал текстовый и видео), и мы не будем подробно на этом останавливаться. Отметим только, что для визуализации BAG файлов и не только есть удобный браузерный визуализатор WebVIZ.

Библиотеки на Python и С++

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

1) Начнем со считывания данных и подготовки их к дальнейшему использованию. Если это .las файлы, можно воспользоваться Python библиотекой pylas.

Pylas не загружает в оперативную память сразу все облако, а создает специальный поток для чтения данных с диска (как memory map в numpy).

Пример чтения .las файла

import pylas

with pylas.open("some_big_file.laz") as f:
    for points in f.chunk_iterator(1_000_000):
        do_something_with(points)

С помощью методов pylas можно создать фильтр, который препроцессирует точки при первичном считывании.

2) Помимо .las формата, облако точек еще часто встречается в .pcd (PCD format) формате. Этот формат - представление данных в библиотеке PCL. Для работы с таким форматом в Python есть библиотека pypcd.

PyPCD в основном нужна для того, чтобы не возиться с кодом на C++ и чтобы считать или записать облака точек в PCL формате.

Пример чтения .pcd с помощью pypcd:

import pypcd
# also can read from file handles.
pc = pypcd.PointCloud.from_path('foo.pcd')
# pc.pc_data has the data as a structured array
# pc.fields, pc.count, etc have the metadata

# center the x field
pc.pc_data['x'] -= pc.pc_data['x'].mean()

# save as binary compressed
pc.save_pcd('bar.pcd', compression='binary_compressed')

Также pypcd совместим с rospy - Python оберткой над ROS.

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

3) PDAL (Point Data Abstraction Library) - С++ библиотека, альтернатива для LAStools, имеющая обертки на Java и Python. В ней реализованы несколько стандартных статистических и структурных фильтров. Подойдет для быстрого чтения и фильтрации больших облаков точек. В PDAL вы в виде конфигурационного файла прописываете пайплайн обработки данных, а после в несколько строк кода запускаете пайплайн.

Пример фильтрации облака точек:

import json
import pdal

def get_pipeline_filtering(path_to_input: str, path_to_output: str, filter:str):
    pipeline = [
            path_to_input,
            {
                "type":"filters." + filter
            },
            {
                "type":"filters.range",
                "limits":"Classification[2:2]"
            },
            path_to_output
    ]

    return pipeline

example = './example.laz'

# Пример фильтрации облака точек c параметрами по умолчанию
method = 'pmf' # Доступные фильтры: https://pdal.io/stages/filters.html
pipeline = get_pipeline_filtering(example, 'result.laz', method)
pipeline = json.dumps(pipeline, indent=2)
pipeline = pdal.Pipeline(pipeline)

if pipeline.validate():
    print("Validate: True")
    pipeline.execute()

4) pyntcloud - простая в использовании Python 3 библиотека для интерактивной визуализации и обработки полигональных мешей и облаков точек. В ней есть инструменты для конвертации форматов, фильтрации и вычисления некоторых скалярных полей над облаками точек (например ближайших соседей, кривизн поверхностей, собственных векторов облака и пр.)

Пример визуализации облака точек, полученного из меша и визуализации кривизны для точек из облака:

from pyntcloud import PyntCloud

anky = PyntCloud.from_file("data/ankylosaurus_mesh.ply")
anky_cloud = anky.get_sample("mesh_random", n=100000, rgb=True, normals=True, as_PyntCloud=True)

anky_cloud.plot()

curvature = anky_cloud.add_scalar_field("curvature", ev=eigenvalues)
anky_cloud.plot(use_as_color=curvature, cmap="jet")

5) Open3D - C++ и Python библиотека, в которой содержатся функции для обработки облаков точек и полигональных мешей. В отличие от всех предыдущих библиотек, Open3D имеет наиболее широкий функционал: это касается и количества поддерживаемых форматов и количества алгоритмов и задач. К библиотеке есть обширная документация с примерами.

Рассмотрим небольшой пример обработки облака точек в формате .ply (данное облако находится в репозитории проекта).

Загрузим и визуализируем исходное облако точек:

pcd = o3d.io.read_point_cloud("../../test_data/fragment.ply")
print(pcd)
print(np.asarray(pcd.points))
o3d.visualization.draw_geometries([pcd],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024])

Вывод:

PointCloud with 196133 points.
[[0.65234375 0.84686458 2.37890625]
 [0.65234375 0.83984375 2.38430572]
 [0.66737998 0.83984375 2.37890625]
 ...
 [2.00839925 2.39453125 1.88671875]
 [2.00390625 2.39488506 1.88671875]
 [2.00390625 2.39453125 1.88793314]]

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

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

downpcd = pcd.voxel_down_sample(voxel_size=0.05)
o3d.visualization.draw_geometries([downpcd],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024])

Вывод:

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

downpcd.estimate_normals(
    search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.1, max_nn=30))
o3d.visualization.draw_geometries([downpcd],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024],
                                  point_show_normal=True)

Вывод:

Практический пример: анализируем дорожные знаки

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

Цель - в потоке лидарных данных обнаружить дорожные знаки и измерить их линейные размеры.

Задача глобально разбивается на несколько этапов:

Анализ данных, пришедших на сервер;
Предобработка данных глобального облака точек проезда;
Предобработка локального облака точек;
Детектирование знака;
Определение пространственных характеристик.

Входные данные, которыми мы располагаем:

Облако точек в форматах .las, .laz глобального скана (проезд автомобиля по целой улице или больше);
Информация о распознанных знаках на маршруте (shapefile с координатами, классами и дополнительными атрибутами знака) - т.е. до того, как анализировать облако точек (3D данные) мы полагаем, что были проанализированы 2D данные проезда (панорамные снимки);
Метаданные (трек положения лидара в пространстве и временные метки);
2D видео проезда;
Размеченные 2D данные (json файл с координатами bounding boxes интересующих объектов);
Панорамные снимки проезда.

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

Исследование входных данных

Для начала визуализируем имеющиеся данные и попробуем понять, какую полезную информацию можно из этого извлечь. Для каждого проезда лидарного комплекса у нас имелись пары .laz и .las файлов. В них отличается информация по типам источника: Intensity, Colour, Point source. Визуализируем с помощью встроенных возможностей Cloud Compare эти данные.

1. Intensity при одинаковом значении Exposure

2. Colour

В .laz, видимо, совсем нет информации т. к. всё чёрное, в .las наоборот информация есть, например, многие знаки ярко выделены. Такое ощущение, что они были светоотражающие и их подсветили фарами (это наблюдение можно использовать и пытаться искать знаки по признаку интенсивности отраженного света).

3. Point source

Можно сделать вывод, что данные получены с разных источников или в разное время суток. В данном случае Colour в .las является достаточно хорошим признаком для выделения знаков.

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

Рис 7. Трек лидара на карте. Для каждой точке геотрека нам известна соответствующая точка в глобальном облаке.

Теперь, используя информацию, полученную от анализа 2D данных (предварительные координаты и классы дорожных знаков, полученные с помощью методов детекции в 2D данных и проецировании точек на основе параметров камеры), проанализируем предварительное расположение объектов интереса.

Кстати, на тему проецирование точек из 2D в 3D на основе параметров камеры есть вот такая замечательная заметка на Medium, которая может быть полезна всем интересующимся исследователям в данной предметной области.

Рис 8. Пример разметки объектов интереса в облаке точек из датасета WAYMO.

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

Теперь у нас есть два способа разбиения глобального облака на локальные:

Брать сферические или кубические подоблака рядом с положением объектов из shape файлов;
Брать сферические подоблака на основе GPS трека.

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

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

Рис 9. Пример входного сферического изображения.

Рис 10. Пример полученного плоского изображения из участка сферического.

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

Разделим облако на сегменты и для каждого сегмента подсчитаем статистику плотности:

Рис 11. Гистограмма распределения количества точек в выделенных объектах (всего 529 объектов).

Всего в нашем распоряжении оказалось 529 объектов. Гистограмма распределения количества точек и типовые объекты приведены на рис.11-12. Среднее количество вершин в объекте 283. Более 75% объектов содержат менее 330-и точек. Помимо самих знаков в таких объектах присутствуют посторонние объекты. Если рассматривать плотность и количество точек только для знаков, то они оказываются еще меньшими.

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

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

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

Рис 13. Применение SOR фильтрации (изображение слева - до фильтрации, изображение справа - после фильтрации).

Для фильтрации воспользовались возможностями ранее упомянутой библиотеки PDAL.

Рассмотрим два базовых типа фильтра: статистические (наиболее типичный представитель - SOR [1] / Рис.13) и поиск опорных поверхностей (наиболее типичный представитель - CSF [2] / Рис. 14).

В зависимости от качества данных, фильтры могли как улучшить ситуацию (убрать посторонние точки), так и ухудшить (убрать точки из объекта интереса).

Основные выводы:

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

Рис 14. Применение CSF фильтра для выделения посторонних вершин, относящихся к опорной поверхности.

Данные для самостоятельных экспериментов

Мы не можем поделиться с вами исходными данными, в силу конфиденциальности проекта, но всем, кто хочет поэкспериментировать со схожими данными, мы рекомендуем обратить внимание на датасет WAYMO, который содержит знаки следующих типов: stop signs, yield signs, speed limit signs, warning signs, guide and recreation signs.

Сегментация знаков внутри под облака

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

Для сегментации были выбраны две модели: DGCNN [3] и PoinNet++ [4]. Кстати, про модели сегментации облаков точек мы писали на хабре вот в этой заметке. DGCNN выбран, как наиболее стабильная и качественная модель, основанная на графовом представлении. PoinNet++ выбрана для сравнения, как классическая модель, основанная на статистических признаках.

Для сравнения были выбраны следующие метрики:

F1-Score - агрегированный критерий качества, учитывающий точность и полноту.
Mean IoU - метрика, показывающая насколько у двух объектов (эталонного (ground true) и текущего) совпадает внутренний “объем”.

Поскольку знаков в реальных данных было не очень много и сегментируемые формы достаточно “простые”, мы написали генератор синтетических данных, который позволил получить в большом количестве знаки разных форм, прикрепленные к столбу.

Рис 15. Пример объекта из синтетического датасета.

Ниже представлен график эволюции mean IoU на этапе обучения модели (размер датасета - 5 тыс., количество эпох - 100). Метрика рассчитывалась в конце каждой эпохи.

Рис 16. PoinNet++ mean IoU (слева), DGCNN mean IoU(справа).

Далее было произведено сравнение обученных моделей (инференс) на тестовом наборе данных из 2 тыс. моделей. Результаты представлены в табл. 1.

Таблица 1 “Метрики качества сегментации на синтетическом датасете”

	F1-score				Доля правильно классифицированных точек
	Плоскость	Основа знака	Треугольник	Круг	Доля правильно классифицированных точек
PoinNet++	0.545420	0.685618	0.700799	0.743857	0.668281
DGCNN	0.807792	0.888354	0.959071	0.847574	0.880837

Измерили время работы на видеокарте (GPU) и на процессоре (CPU):

Таблица 2 “Время сегментирования на этапе инференса”

	CPU, c	GPU, c
PoinNet++	0,38	0,01
DGCNN	1,28	0,04

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

На синтетических данных:

На реальных данных (предварительно отфильтрованных):

Здесь стоит отметить, что на самом деле данное облако точек не является круговым сегментом плоскости, а представляет собой заполнение цилиндрического объема.

В результате мы можем наблюдать следующую картину:

DGCNN справляется с задачей сегментации лучше, чем PointNet++;
Обученная на синтетических данных модель неудовлетворительно справляется с сегментацией реальных данных.

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

Рис 17. Найденная плоскость внутри облака точек.

Расчет пространственных характеристик

После того, как мы выделили точки, которые предположительно относятся исключительно только к объекту интереса, мы можем приступать к определению пространственных характеристик. Сама по себе задача определения различных пространственных характеристик облака точек не тривиальна. Так, например, в работе [5] авторы определяют множество точек внутри облака, которые соответствуют границе объекта. А в работе [6]авторы использовали механизм геометрического внимания и глубокие архитектуры для того, чтобы выделять острые края объектов.

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

Вторым подходом является использование готового образца (например гостированный знак в виде параметрической или полигональной модели) для “наложения его на имеющиеся облако точек”. При таком подходе можно в качестве пространственных характеристик облака использовать пространственные характеристики образца (с поправкой на масштабирование). Такое сопоставление можно производить с помощью ICP алгоритмов, реализация которых есть, например, в trimesh.

Рис 18. Результат наложение образца на облако методом ICP.

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

Рис 19. Измерение линейного размера дорожного знака в Cloud Compare.

Заключение

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

Источники

Rusu, Radu Bogdan, et al. “Towards 3D point cloud based object maps for household environments.” Robotics and Autonomous Systems 56.11 (2008): 927-941. [paper]
Zhang, Wuming, et al. “An easy-to-use airborne LiDAR data filtering method based on cloth simulation.” Remote Sensing 8.6 (2016): 501. [paper]
Phan, A.V., Le Nguyen, M., Nguyen, Y.L.H. and Bui, L.T., 2018. DGCNN: A convolutional neural network over large-scale labeled graphs. Neural Networks, 108, pp.533-543. [paper]
Qi, C.R., Yi, L., Su, H. and Guibas, L.J., 2017. Pointnet++: Deep hierarchical feature learning on point sets in a metric space. In Advances in neural information processing systems (pp. 5099-5108). [paper]
Mineo C., Pierce S. G., Summan R. Novel algorithms for 3D surface point cloud boundary detection and edge reconstruction //Journal of Computational Design and Engineering. – 2019. – Т. 6. – №. 1. – С. 81-91. [paper]
Matveev A. et al. Geometric Attention for Prediction of Differential Properties in 3D Point Clouds //IAPR Workshop on Artificial Neural Networks in Pattern Recognition. – Springer, Cham, 2020. – С. 113-124. [paper]

Осваиваем анализ лидарных данных и измеряем дорожные знаки