Вычислить ограничивающий многоугольник альфа-формы из триангуляции Делоне

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

http://bl.ocks.org/gka/1552725

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

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

Ответ 1

Создайте граф, в котором узлы соответствуют треугольникам Делоне и в которых есть грань графика между двумя треугольниками тогда и только тогда, когда они разделяют две вершины.
Вычислить связанные компоненты графа.
Для каждого подключенного компонента найдите все узлы, которые имеют менее трех соседних узлов (то есть те, которые имеют степень 0, 1 или 2). Они соответствуют граничным треугольникам. Мы определяем ребра граничного треугольника, которые не разделяются с другим треугольником, как граничные ребра этого граничного треугольника.

В качестве примера я выделил эти граничные треугольники в вашем примере "вопросительный знак" триангуляции Делане:

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

Ответ 2

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

from scipy.spatial import Delaunay
import numpy as np


def alpha_shape(points, alpha, only_outer=True):
    """
    Compute the alpha shape (concave hull) of a set of points.
    :param points: np.array of shape (n,2) points.
    :param alpha: alpha value.
    :param only_outer: boolean value to specify if we keep only the outer border
    or also inner edges.
    :return: set of (i,j) pairs representing edges of the alpha-shape. (i,j) are
    the indices in the points array.
    """
    assert points.shape[0] > 3, "Need at least four points"

    def add_edge(edges, i, j):
        """
        Add a line between the i-th and j-th points,
        if not in the list already
        """
        if (i, j) in edges or (j, i) in edges:
            # already added
            assert (j, i) in edges, "Can't go twice over same directed edge right?"
            if only_outer:
                # if both neighboring triangles are in shape, it not a boundary edge
                edges.remove((j, i))
            return
        edges.add((i, j))

    tri = Delaunay(points)
    edges = set()
    # Loop over triangles:
    # ia, ib, ic = indices of corner points of the triangle
    for ia, ib, ic in tri.vertices:
        pa = points[ia]
        pb = points[ib]
        pc = points[ic]
        # Computing radius of triangle circumcircle
        # www.mathalino.com/reviewer/derivation-of-formulas/derivation-of-formula-for-radius-of-circumcircle
        a = np.sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
        b = np.sqrt((pb[0] - pc[0]) ** 2 + (pb[1] - pc[1]) ** 2)
        c = np.sqrt((pc[0] - pa[0]) ** 2 + (pc[1] - pa[1]) ** 2)
        s = (a + b + c) / 2.0
        area = np.sqrt(s * (s - a) * (s - b) * (s - c))
        circum_r = a * b * c / (4.0 * area)
        if circum_r < alpha:
            add_edge(edges, ia, ib)
            add_edge(edges, ib, ic)
            add_edge(edges, ic, ia)
    return edges

Если вы запустите его со следующим тестовым кодом, вы получите следующую цифру:

from matplotlib.pyplot import *

# Constructing the input point data
np.random.seed(0)
x = 3.0 * np.random.rand(2000)
y = 2.0 * np.random.rand(2000) - 1.0
inside = ((x ** 2 + y ** 2 > 1.0) & ((x - 3) ** 2 + y ** 2 > 1.0) & ((x - 1.5) ** 2 + y ** 2 > 0.09))
points = np.vstack([x[inside], y[inside]]).T

# Computing the alpha shape
edges = alpha_shape(points, alpha=0.25, only_outer=True)

# Plotting the output
figure()
axis('equal')
plot(points[:, 0], points[:, 1], '.')
for i, j in edges:
    plot(points[[i, j], 0], points[[i, j], 1])
show()

Ответ 3

Получается, что TopoJSON имеет алгоритм слияния, который выполняет именно эту задачу: https://github.com/mbostock/topojson/wiki/API-Reference#merge

Там даже пример показывает его в действии: http://bl.ocks.org/mbostock/9927735

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

Ответ 4

Основываясь на ответе @Timothy, я использовал следующий алгоритм для вычисления граничных колец триангуляции Делоне.

from matplotlib.tri import Triangulation
import numpy as np

def get_boundary_rings(x, y, elements):
    mpl_tri = Triangulation(x, y, elements)
    idxs = np.vstack(list(np.where(mpl_tri.neighbors == -1))).T
    unique_edges = list()
    for i, j in idxs:
        unique_edges.append((mpl_tri.triangles[i, j],
                             mpl_tri.triangles[i, (j+1) % 3]))
    unique_edges = np.asarray(unique_edges)
    ring_collection = list()
    initial_idx = 0
    for i in range(1, len(unique_edges)-1):
        if unique_edges[i-1, 1] != unique_edges[i, 0]:
            try:
                idx = np.where(
                    unique_edges[i-1, 1] == unique_edges[i:, 0])[0][0]
                unique_edges[[i, idx+i]] = unique_edges[[idx+i, i]]
            except IndexError:
                ring_collection.append(unique_edges[initial_idx:i, :])
                initial_idx = i
                continue
    # if there is just one ring, the exception is never reached,
    # so populate ring_collection before returning.
    if len(ring_collection) == 0:
        ring_collection.append(np.asarray(unique_edges))
    return ring_collection

Ответ 5

Немного пересмотреть ответ ханниэля для случая с 3 точками (тетраэдр).

def alpha_shape(points, alpha, only_outer=True):
    """
    Compute the alpha shape (concave hull) of a set of points.
    :param points: np.array of shape (n, 3) points.
    :param alpha: alpha value.
    :param only_outer: boolean value to specify if we keep only the outer border
    or also inner edges.
    :return: set of (i,j) pairs representing edges of the alpha-shape. (i,j) are
    the indices in the points array.
    """
    assert points.shape[0] > 3, "Need at least four points"

    def add_edge(edges, i, j):
        """
        Add a line between the i-th and j-th points,
        if not in the set already
        """
        if (i, j) in edges or (j, i) in edges:
            # already added
            if only_outer:
                # if both neighboring triangles are in shape, it not a boundary edge
                if (j, i) in edges:
                    edges.remove((j, i))
            return
        edges.add((i, j))

    tri = Delaunay(points)
    edges = set()
    # Loop over triangles:
    # ia, ib, ic, id = indices of corner points of the tetrahedron
    print(tri.vertices.shape)
    for ia, ib, ic, id in tri.vertices:
        pa = points[ia]
        pb = points[ib]
        pc = points[ic]
        pd = points[id]

        # Computing radius of tetrahedron Circumsphere
        # http://mathworld.wolfram.com/Circumsphere.html

        pa2 = np.dot(pa, pa)
        pb2 = np.dot(pb, pb)
        pc2 = np.dot(pc, pc)
        pd2 = np.dot(pd, pd)

        a = np.linalg.det(np.array([np.append(pa, 1), np.append(pb, 1), np.append(pc, 1), np.append(pd, 1)]))

        Dx = np.linalg.det(np.array([np.array([pa2, pa[1], pa[2], 1]),
                                     np.array([pb2, pb[1], pb[2], 1]),
                                     np.array([pc2, pc[1], pc[2], 1]),
                                     np.array([pd2, pd[1], pd[2], 1])]))

        Dy = - np.linalg.det(np.array([np.array([pa2, pa[0], pa[2], 1]),
                                       np.array([pb2, pb[0], pb[2], 1]),
                                       np.array([pc2, pc[0], pc[2], 1]),
                                       np.array([pd2, pd[0], pd[2], 1])]))

        Dz = np.linalg.det(np.array([np.array([pa2, pa[0], pa[1], 1]),
                                     np.array([pb2, pb[0], pb[1], 1]),
                                     np.array([pc2, pc[0], pc[1], 1]),
                                     np.array([pd2, pd[0], pd[1], 1])]))

        c = np.linalg.det(np.array([np.array([pa2, pa[0], pa[1], pa[2]]),
                                    np.array([pb2, pb[0], pb[1], pb[2]]),
                                    np.array([pc2, pc[0], pc[1], pc[2]]),
                                    np.array([pd2, pd[0], pd[1], pd[2]])]))

        circum_r = math.sqrt(math.pow(Dx, 2) + math.pow(Dy, 2) + math.pow(Dz, 2) - 4 * a * c) / (2 * abs(a))
        if circum_r < alpha:
            add_edge(edges, ia, ib)
            add_edge(edges, ib, ic)
            add_edge(edges, ic, id)
            add_edge(edges, id, ia)
            add_edge(edges, ia, ic)
            add_edge(edges, ib, id)
    return edges

Ответ 6

Альфа-формы определяются как триангуляция delaunay без границ, превышающих альфа. Сначала удалите все внутренние треугольники, а затем все ребра, превышающие альфа.