Удаление шума и деформация 3D‑облаков точек в Python с помощью графовых лапласианов / Хабр

Облако точек — это набор измерений в формате точек в 3D‑системе координат, где каждая точка представляет одно пространственное измерение на поверхности некоторого объекта. Вместе все точки описывают внешнюю поверхность объекта в пространстве. Такие облака точек можно получить с помощью 3D‑лазерных сканеров, LiDAR-сканеров или фотограмметрии; их часто используют для построения точных цифровых 3D‑представлений в строительстве и архитектуре, геодезии и картографии для создания цифровых двойников зданий или моделей рельефа и местности, а также в промышленности для создания 3D-моделей деталей (CAD), контроля качества продукции и анализа износа оборудования.

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

Ключевое отличие облаков точек от 3D‑сеток (mesh) в том, что облака точек — это несвязанные точки данных, а не поверхность, построенная из треугольников или полигонов (рисунок 1).

Рисунок 1: 3D‑модели Stanford Bunny: сетка (слева) и облако точек (справа)

Таким образом, в отличие от сеточных моделей, у облаков точек нет явной информации о связности (рёбра и грани), из‑за ч��го трудно вычислять на них производные [1]. Главная проблема — извлечь глобальную информацию из облака точек при отсутствии связности. Графовый лапласиан позволяет обойти эту проблему, задавая связность посредством построения графа (например, по методу ближайших соседей) и определения дискретного оператора, который аппроксимирует непрерывные геометрические свойства (кривизну, нормали поверхности и т.п.) на дискретных данных.

Графовый лапласиан широко примененяется в LiDAR-SLAM, навигации автономных систем, а также в задачах реконструкции поверхностей. Он эффективен при генерации сеток, триангуляции разрывов и восстановлении поврежденных 3D-сканов

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

План:

Облака точек и kNN‑графы: от kNN к разреженной матрице смежности
Вычисление графового лапласиана
Стабильное сглаживание + простое притяжение к ближайшему соседу
Визуализация деформации
Тестирование кода

Облака точек и kNN‑графы: от kNN к разреженной матрице смежности

Графовый лапласиан — это дифференциальный оператор, заданный на графе, построенном на точках (в отличие от сеточной поверхности или непрерывной области). Он задаёт дискретный оператор диффузии на графе, который можно использовать для сглаживания и регуляризации деформаций. Мы вычислим его и применим для трансформации облаков точек из репозитория Stanford 3D Scanning [3].

Разобьём задачу на 6 шагов:

1. Построить граф по точкам с помощью ближайших соседей (kNN).

2. Построить разреженную матрицу , вычислив веса $w_{ij}$ для пар соседей.

3. Вычислить матрицу степеней (degree matrix) , суммируя веса построчно.

4. Вычислить лапласиан (ненормированный) или $L_{norm} = \mathbf{I} - D^{-1/2} W D^{-1/2}$ (нормированный).

5. Деформировать форму 1 в форму 2, минимизируя энергию.

6. Визуализировать процесс деформации с помощью Plotly.

Для вычислений и визуализации мы будем использовать несколько модулей Python:

import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as sl
from sklearn.neighbors import NearestNeighbors
import plotly.graph_objects as go
import open3d as o3d
from tests.tests import *

Последний импорт — это файл tests.py, в котором живут модульные (unit) тесты. Его можно закомментировать, если нет желания производить проверки кода.

Пусть $P = \{p_i\}_i^n$ , где $p_i \in \mathbb{R}^3$ , — облако точек. Поскольку у облака точек нет встроенной связности, чтобы вычислить лапласиан, нужно «симулировать поверхность», определив, какие точки являются соседями. Это можно сделать, построив граф ближайших соседей (k‑NN). Для каждой точки в исходной форме найдём её ближайших соседей (например, ).

Соседи задают локальную топологию через матрицу весов $W \in \mathbb{R}^{n \times n}$ с элементами

$w_{ij}=\begin{cases} e^{-\frac{||p_i - p_j||^2}{2\sigma^2}} \quad \text{если $i, j$ — соседи}\\ 0 \quad \quad \quad \quad \text{иначе.} \end{cases}$

Таким образом, содержит информацию о смежности и о локальной структуре: более близким соседям соответствует больший вес. Так как мы ищем только $k \ll n$ ближайших точек, матрица будет разреженной.

Ещё одна матрица, которую потребуется вычислить на этом этапе, — матрица степеней $D$. Это диагональная матрица с ненулевыми элементами $d_{ii} = \sum_j w_{ij}$ , то есть суммой весов для точки .

Далее вычислим графовый лапласиан (ненормированный вариант), который фиксирует локальную форму, сравнивая точку со средним по её соседям. Мы также добавим опцию для сравнения его с нормированным лапласианом $L_{norm} = \mathbf{I} - D^{-1/2} W D^{-1/2}$ , который иногда бывает устойчивее при неравномерной дискретизации.

Вычисление графового лапласиана

Для реализации этих шагов в Python мы определим функцию build_knn_laplacian(). Она принимает облако точек в виде матрицы, число соседей и (опционально) параметр $\sigma$ для гауссовых весов, опцию симметризации матрицы , а также опцию нормирования лапласиана. Функция строит граф, находя соседей с помощью sklearn.neighbors.NearestNeighbors, вычисляет евклидовы расстояния между соседними точками и собирает разреженную матрицу , применяя гауссову функцию weights = np.exp(-(dists**2) / (2.0 sigma*2)) поэлементно к матрице расстояний между соседями. Затем мы определим диагональную матрицу , суммируя веса по строкам , и вычислим два варианта лапласиана — нормированный и ненормированный.

def build_knn_laplacian(P, k=20, sigma=None, symmetrize=True, normalize=False):
    
    n = P.shape[0]
    nn = NearestNeighbors(n_neighbors=k+1, algorithm="auto").fit(P)
    dists, inds = nn.kneighbors(P)

    # Исключить саму точку из списка соседей
    dists = dists[:, 1:]
    inds = inds[:, 1:]

    # назначить дефолтное значение sigma (медианное расстояние между соседями)
    if sigma is None:
        sigma = np.median(dists)
        sigma = max(sigma, 1e-8)

    # веса: исарользовать функцию Гаусса к расстояниям между соседями
    weights = np.exp(-(dists**2) / (2.0 * sigma**2))

    # построить разреженную матрицу W
    rows = np.repeat(np.arange(n), k)
    cols = inds.reshape(-1)
    vals = weights.reshape(-1)

    W = sp.coo_matrix((vals, (rows, cols)), shape=(n, n)).tocsr()

    if symmetrize:
        # применить max(W, W^T) для симметрии
        W = W.maximum(W.T)

    d = np.array(W.sum(axis=1)).reshape(-1)
    d = np.maximum(d, 1e-12)
    D = sp.diags(d, offsets=0, shape=(n, n), format="csr")
    
    if normalize:
        D_inv_sqrt = sp.diags(1.0 / np.sqrt(d))
        L = sp.eye(n, format="csr") - (D_inv_sqrt @ W @ D_inv_sqrt)
    else:
        L = D - W
    
    return L, d

Стабильное сглаживание + простое притяжение к ближайшему соседу

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

Для матрицы координат точек $P \in \mathbb{R}^{n \times 3}$ , когда применяется к каждой координатной компоненте, для каждой точки получаются лапласиановы координаты

$\delta_i = \sum_j w_{ij} (p_i - p_j),$

которые содержат информацию о «локальных деталях» или кривизне в окрестности точки.

Эту информацию можно использовать, чтобы преобразовать одну форму облака точек, скажем $P \in \mathbb{R}^{n \times 3}$ , в другую, скажем $T \in \mathbb{R}^{m \times 3}$ . Ключевая идея заключается в том, чтобы сохранить $\delta_P$ (локальные детали исходной формы), одновременно постепенно двигаясь к целевой форме. Обычно это делают через минимизацию энергии на каждом шаге времени $t = 1, \dots, \mathcal{T}$ :

$\min_{P_t} \big(L P_t - \delta(t)^2 + \lambda||P_t - T(t)||^2 \big), \quad \quad (*)$

где — положения точек в момент времени , $\delta(t) = (1-t) \delta_P + t \delta_T$ линейная интерполяция лапласиановых координат между исходной и целевой формой, $T(t) = (1-t)P + t T_{corr}$ — позиционная интерполяция к целевым точкам (с учётом соответствий), а $\lambda$ играет роль параметра «жёсткости»: большое значение $\lambda$ заставляет форму точнее совпадать с целью, но может «исказить» локальные детали; маленькое значение держит форму более «жёсткой» (ближе к исходной). Первое слагаемое в (*) отвечает за сохранение локальных дифференциалов (сглаживание/денойзинг), а второе — за притяжение к целевым положениям.

(*) — это стандартная задача наименьших квадратов с регуляризацией, у которой есть решение в закрытой форме. Нормальные уравнения выглядят так:

$(L^T L +\lambda I)P' = L^T \delta + \lambda T. \quad \quad (**)$

Пока мы не определили матрицу соответствий $T_{corr}$ и интерполированные лапласиановы координаты $\delta$ . $T_{corr}$ показывает, какая точка исходной формы «хочет» двигаться к какой точке целевой формы . Задать это соответствие можно разными способами. Если облака имеют одинаковый порядок и размер, можно использовать отображение 1:1. Но такие ситуации — скорее редкость, и часто применяют сопоставление по методу ближайшего соседа. Ещё один вариант — неригидная регистрация (CPD, варианты ICP), но чтобы не усложнять, мы снова используем kNN.

Вычисление $T_{corr}$ через проекцию на ближайшего соседа естественным образом поддерживает разные размерности матриц и , задавая отображение

$\phi:\{1,\dots, n\} \rightarrow \{1,\dots,m\},$

где для каждой точки мы выбираем ближайшую целевую точку: $\phi(i) = \text{argmin}_{j \in 1,\dots,m}||p_i - t_j||$ и определяем соответствующую целевую позицию как $t^{corr}_i = t_{\phi(i)}$ , так что $T_{corr} \in \mathbb{R}^{n \times 3}$ . при такой формулировке возможны и отображения «многие‑к‑одному», и «один‑ко‑многим».

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

def recompute_Tcorr(Pcur, T):
    nn = NearestNeighbors(n_neighbors=1, algorithm="auto").fit(T)
    _, idx = nn.kneighbors(Pcur)
    return T[idx[:, 0]]

Дифференциальное смещение во времени , $\delta(t)$ , — это простая линейная интерполяция между векторами смещения источника ( $\delta_P$ ) и цели ( $\delta_T$ ) для каждой точки облака.

Чтобы реализовать этот шаг в Python, мы определим функцию laplacian_morph_frames(), которая принимает исходную матрицу точек и целевую , параметр жёсткости $\lambda$ , число шагов (итераций), интервал обновления лапласиана и соответствий, а также (опционально) $\sigma$ для вычисления лапласиана.

После вычисления дифференциальных координат deltaP = L_current @ Pcur и deltaT = L_current @ Tcorr мы решим систему нормальных уравнений (**).

Для решения нормальных уравнений мы используем LU‑разложение (для удобного и устойчивого решения разреженной СЛАУ) через scipy.sparse.linalg.factorized() — прекрасно работает для разреженных матриц. Затем вызываем метод solve, который использует уже построенное разложение.

Мы будем пересчитывать лапласиан build_knn_laplacian(Pcur, k=k, sigma=sigma, symmetrize=True, normalize=False) каждые refresh_every шагов внутри временного цикла. По мере деформации Pcur мы обновляем и граф, чтобы «поверхностная связность», по которой идёт сглаживание, не оставалась замороженной в исходной геометрии. Если не обновлять, связи между точками будут соответствовать исходной форме, даже когда она деформируется в цель — это плохо согласуется с физическим смыслом. Обновление позволяет этому «псевдосеточному» соседству адаптироваться к новой форме.

def laplacian_morph_frames(
    P, T, k=20, lam=10.0, n_frames=30,
    refresh_every=5, sigma=None, normalize=True
):

    n = P.shape[0]
    I = sp.eye(n, format="csr")

    ts = np.linspace(0.0, 1.0, n_frames)
    frames = []

    Pcur = P.copy()

    # создать плейсхолдеры
    L_current = None
    solve_P = None
    deltaP = None
    deltaT = None
    Tcorr = None

    for i, t in enumerate(ts):
        if (i % refresh_every == 0) or (solve_P is None):
            # 1) обновить корреспо��денции
            Tcorr = recompute_Tcorr(Pcur, T)  # (n,3)

            # 2) обновить оператор
            L_current, _ = build_knn_laplacian(Pcur, k=k, sigma=sigma, symmetrize=True, normalize=normalize)

            # 3) обновить лапласианские координаты + факторизовать солвер
            deltaP = L_current @ Pcur
            deltaT = L_current @ Tcorr

            A = (L_current.T @ L_current) + lam * I
            solve_P = sl.factorized(A.tocsc())

        # интерполировать координаты
        delta_t = (1.0 - t) * deltaP + t * deltaT
        Tt = (1.0 - t) * P + t * Tcorr

        rhs = (L_current.T @ delta_t) + lam * Tt

        Pnext = np.empty_like(Pcur)
        for dim in range(3):
            Pnext[:, dim] = solve_P(rhs[:, dim])

        frames.append(Pnext)
        Pcur = Pnext

    return ts, frames

Визуализация деформации

Для визуализации процесса трансформации можно использовать 3D‑график Plotly со слайдером, который позволяет прокручивать кадры деформации как шаги времени. Облако точек будем рисовать через graphical_objects.Scatter3d() по координатам , и : инициализируем graphical_objects.Figure, задаём начальный кадр (frames[0]), определяем слайдер, добавляя остальные кадры от 1 до n_frames и описывая шаги анимации, затем обновляем layout, чтобы создать слайдер и кнопки управления (“Pause”, “Play”), и выставляем нужный ракурс камеры. Цвет каждой точки задаётся расстоянием от начала координат — это помогает отслеживать движение разных частей формы.

def plotly_slider(frames, title="Laplacian Morph", marker_size=2):

    # стабилизировать границы камеры во всех кадрах
    all_pts = np.vstack(frames)
    xmin, ymin, zmin = all_pts.min(axis=0)
    xmax, ymax, zmax = all_pts.max(axis=0)

    def dist_from_origin(X):
        return np.sqrt((X**2).sum(axis=1))

    # Исходный фрейм
    X0 = frames[0]
    c0 = dist_from_origin(X0)

    fig = go.Figure(
        data=[
            go.Scatter3d(
                x=X0[:, 0], y=X0[:, 1], z=X0[:, 2],
                mode="markers",
                marker=dict(
                    size=marker_size,
                    color=c0,
                    colorscale="Viridis",
                    showscale=True,
                    colorbar=dict(title="||x||")
                ),
                name="Deformed"
            )
        ]
    )

    # последующие кадры (последовательность для анимации)
    plotly_frames = []
    for i, Xt in enumerate(frames):
        ct = dist_from_origin(Xt)
        plotly_frames.append(
            go.Frame(
                data=[go.Scatter3d(
                    x=Xt[:, 0], y=Xt[:, 1], z=Xt[:, 2],
                    marker=dict(color=ct) 
                )],
                name=str(i)
            )
        )

    fig.frames = plotly_frames

    # слайдер - шаги, кнопки, ракурс камеры
    steps = [
        dict(
            method="animate",
            args=[[str(i)], dict(mode="immediate",
                                frame=dict(duration=0, redraw=True),
                                transition=dict(duration=0))],
            label=str(i),
        )
        for i in range(len(frames))
    ]

    fig.update_layout(
        title=title,
        sliders=[dict(active=0, currentvalue={"prefix": "frame: "}, pad={"t": 30}, steps=steps)],
        updatemenus=[dict(
            type="buttons",
            showactive=False,
            buttons=[
                dict(label="Play", method="animate",
                     args=[None, dict(frame=dict(duration=60, redraw=True),
                                      transition=dict(duration=0),
                                      fromcurrent=True)]),
                dict(label="Pause", method="animate",
                     args=[[None], dict(frame=dict(duration=0, redraw=False),
                                       mode="immediate")]),
            ],
            x=0.0, y=1.0
        )],
        scene=dict(                       # установить ракурс камеры
        xaxis=dict(autorange='reversed'), 
        camera=dict(
            up=dict(x=0, y=1, z=0),      
            eye=dict(x=0, y=0, z=2),      
            center=dict(x=0, y=0, z=0)
        ),
        dragmode='orbit'                
    )
    )
    
    fig.show()

Теперь сделаем обёртку над всеми функциями: эта функция принимает данные облаков точек (исходное и целевое), число соседей для kNN, параметр жёсткости, и число шагов. Затем она нормализует облака, находит соответствия между исходными и целевыми точками, вычисляет лапласиан, решает задачу минимизации энергии и визуализирует результат.

Мы будем использовать o3d.io.read_point_cloud() для чтения облака точек и затем преобразуем его в numpy‑массив типа float64. Во многих задачах перед вычислением графового лапласиана оправдана нормализация данных к единичной сфере. Хорошая практика в задачах деформации — нормализовать оба облака одним и тем же преобразованием, используя центр и масштаб исходной формы для каждого облака. Нормализация критична и для вычисления весов смежности при построении графа (предполагается примерно один и тот же масштаб координат для исходной и целевой формы), и для величины «притяжения к цели» (шаг/сила зависит от масштаба координат).

# загрузить данные
pcd_init = o3d.io.read_point_cloud("data/bunny.ply")
initial_points = np.asarray(pcd_init.points, dtype=np.float64)

pcd_tgt = o3d.io.read_point_cloud("data/dragon.ply")
target_points = np.asarray(pcd_tgt.points, dtype=np.float64)


def morph_ply_pointclouds(P, T, k=20, lam=10.0, n_frames=11, 
                          refresh_every=5, sigma=None, normalize=True):

    # нормализовать
    mu = P.mean(axis=0, keepdims=True)
    P0 = P - mu
    scale = np.sqrt((P0**2).sum(axis=1)).mean()

    P = P0 / (scale + 1e-12)
    T = (T - mu) / (scale + 1e-12)
    
    # вычислить деформации (фреймы для Plotly)
    ts, frames = laplacian_morph_frames(P, T, k=k, lam=lam, n_frames=n_frames,
    refresh_every=refresh_every, sigma=sigma, normalize=normalize)
    
    # визуализация
    plotly_slider(frames, title=f"Laplacian morph")

Теперь осталось лишь вызвать функцию morph_ply_pointclouds и наблюдать деформацию.

morph_ply_pointclouds(initial_points, target_points)

Бонус: unit-тестирование

Для деформации облаков точек через лапласиан можно использовать unit‑тесты, которые ловят неправильный шаблон разреженности, сломанную симметрию, разорванные графы, неустойчивые решения СЛАУ и ошибки в соответствии.

Мы пройдёмся по 6 полезным тестам, которые помогают поймать несколько вещей, которые могут незаметно сломаться.

1) Тесты построения лапласиана (форма + разреженность + свойства)

Тест 1 проверяет размерность лапласиана: должно быть $n \times n$ .

Тест 2 проверяет разреженность: ожидаемое число ненулевых элементов — порядка $n \cdot k$ после симметризации.

Тест 3a проверяет свойство нулевой суммы по строкам: для ненормированного графового лапласиана ожидаем: $L\mathbf{I}=0$ .

Тест 3b проверяет свойство нулевой суммы по строкам: для нормированного графового лапласиана $L=\mathbf{I} - D^{-1/2}W D^{-1/2}$ ожидаем: $L_{norm}(D^{1/2} \mathbf{I})=0$ .

2) Тест дифференциальных координат

Тест 4 проверяет инвариантность лапласиановых координат к переносу. Ожидаем

, то есть координаты не меняются, если всё облако точек сдвинуть на вектор . Мы тестируем это, генерируя случайный сдвиг и подтверждая, что $L(X+c) \approx LX$ .

3) Сквозные поведенческие тесты (на маленьких синтетических облаках)

Тест 5 проверяет трансформацию $P \rightarrow P$ . Если взять синтетическое облако (или данное) и задать , то каждый кадр (в том числе, последний) должен быть $\approx P$ (или очень близок к нему).

Тест 6 проверяет предсказуемое поведение деформации. Мы строим , масштабируя координату : $T[:,2] = 2 \cdot P[:,2]$ . Ожидаем, что расстояние монотонно убывает по .

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

Описанные 6 тестов собраны в tests/tests.py.

Полный код + данные + описание: https://github.com/yuliadm/pointcloud_laplacian

Ссылки

1. Liang, J., Lai, R., Wong, T.W. and Zhao, H., 2012, June. Geometric understanding of point clouds using Laplace-Beltrami operator. In 2012 IEEE conference on computer vision and pattern recognition (pp. 214-221). IEEE. https://ww3.math.ucla.edu/camreport/cam12-26.pdf

2. Alexiou, E., Ebrahimi, T., Bernardo, M.V., Pereira, M., Pinheiro, A., Cruz, L.A.D.S., Duarte, C., Dmitrovic, L.G., Dumic, E., Matkovics, D. and Skodras, A., 2018, May. Point cloud subjective evaluation methodology based on 2D rendering. In 2018 Tenth international conference on quality of multimedia experience (QoMEX) (pp. 1-6). IEEE.

https://www.epfl.ch/labs/mmspg/downloads/reconstructed-point-clouds-results/

3. Stanford 3D Scanning Repository: https://graphics.stanford.edu/data/3Dscanrep/