d3.js 从Delaunay三角剖分计算alpha形状的边界多边形

rdrgkggo  于 2023-11-19  发布在  其他
关注(0)|答案(8)|浏览(264)

给定平面上的一组点,对于给定的正数alpha,alpha形状的概念通过找到Delaunay三角剖分并删除任何至少一条边的长度超过alpha的三角形来定义。下面是一个使用d3的例子:
http://bl.ocks.org/gka/1552725
问题是,当有成千上万个点时,简单地绘制所有内部三角形对于交互式可视化来说太慢了,所以我想只找到边界多边形。这并不简单,因为正如你从那个例子中看到的,有时可能有两个这样的多边形。
为了简化,假设已经执行了一些聚类,以保证每个三角剖分都有一个唯一的边界多边形。找到这个边界多边形的最佳方法是什么?特别是,边缘必须一致地排序,并且它必须支持“洞”的可能性(想想环面或甜甜圈形状-这在GeoJSON中可以表达)。

zphenhs4

zphenhs41#

下面是一些Python代码,可以满足您的需要。我修改了here的alpha形状(凹船体)计算,使其不插入内边(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 is 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.simplices:
        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

字符串
如果你运行下面的测试代码,你会得到this figure

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()

waxmsbnn

waxmsbnn2#

1.创建一个图,其中节点对应于Delaunay三角形,并且当且仅当两个三角形共享两个顶点时,两个三角形之间存在图边。
1.计算图的连通分量。
1.对于每个连通分支,找到所有相邻节点少于3个的节点(即度为0、1或2的节点)。这些节点对应于边界三角形。我们将边界三角形的边定义为该边界三角形的边界边
作为一个例子,我在你的“问号”Delaunay三角剖分示例中突出显示了这些边界三角形:


的数据
根据定义,每个边界三角形最多与另外两个边界三角形相邻。边界三角形的边界边形成循环。您可以简单地遍历这些循环来确定边界的多边形形状。如果在实现中记住了这些,这也适用于带孔的多边形。

e4eetjau

e4eetjau3#

我知道这是一个延迟的答案,但这里张贴的方法没有为我工作的各种原因。
上面提到的包Alphashape总体上是不错的,但它的缺点是它使用Shapely的cascade_union和三角测量方法来给予最终的凹形船体,这对于大型数据集和低alpha值来说非常慢(高精度)。在这种情况下,Iddo Hanniel发布的代码(以及Harold的修订版)将在内部保留大量的边缘,唯一的方法是使用前面提到的cascade_union和Shapely的三角剖分。
我通常使用复杂的表单,下面的代码运行良好,速度很快(100K 2D点需要2秒):

import numpy as np
from shapely.geometry import MultiLineString
from shapely.ops import unary_union, polygonize
from scipy.spatial import Delaunay
from collections import Counter
import itertools

def concave_hull(coords, alpha):  # coords is a 2D numpy array

    # i removed the Qbb option from the scipy defaults.
    # it is much faster and equally precise without it.
    # unless your coords are integers.
    # see http://www.qhull.org/html/qh-optq.htm
    tri = Delaunay(coords, qhull_options="Qc Qz Q12").vertices

    ia, ib, ic = (
        tri[:, 0],
        tri[:, 1],
        tri[:, 2],
    )  # indices of each of the triangles' points
    pa, pb, pc = (
        coords[ia],
        coords[ib],
        coords[ic],
    )  # coordinates of each of the triangles' points

    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) * 0.5  # Semi-perimeter of triangle

    area = np.sqrt(
        s * (s - a) * (s - b) * (s - c)
    )  # Area of triangle by Heron's formula

    filter = (
        a * b * c / (4.0 * area) < 1.0 / alpha
    )  # Radius Filter based on alpha value

    # Filter the edges
    edges = tri[filter]

    # now a main difference with the aforementioned approaches is that we dont
    # use a Set() because this eliminates duplicate edges. in the list below
    # both (i, j) and (j, i) pairs are counted. The reasoning is that boundary
    # edges appear only once while interior edges twice
    edges = [
        tuple(sorted(combo)) for e in edges for combo in itertools.combinations(e, 2)
    ]

    count = Counter(edges)  # count occurrences of each edge

    # keep only edges that appear one time (concave hull edges)
    edges = [e for e, c in count.items() if c == 1]

    # these are the coordinates of the edges that comprise the concave hull
    edges = [(coords[e[0]], coords[e[1]]) for e in edges]

    # use this only if you need to return your hull points in "order" (i think
    # its CCW)
    ml = MultiLineString(edges)
    poly = polygonize(ml)
    hull = unary_union(list(poly))
    hull_vertices = hull.exterior.coords.xy

    return edges, hull_vertices

字符串

6qfn3psc

6qfn3psc4#

现在有一个非常容易使用的python包alphashape,可以通过pipconda安装。
main函数的输入与@Iddo Hanniel给出的答案类似,调整第二个位置参数将为您提供所需的轮廓。或者,您可以将seconda位置参数留空,函数将为您优化该参数,以提供最佳凹船体。请注意,如果让函数优化值,则计算时间会大大增加。

nxagd54h

nxagd54h5#

Hanniel's answer的轻微修改,用于3d点的情况(四面体)。

from scipy.spatial import Delaunay
import numpy as np
import math

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's 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

字符串

xggvc2p6

xggvc2p66#

TopoJSON有一个合并算法,它只执行以下任务:https://github.com/mbostock/topojson/wiki/API-Reference#merge
甚至有一个例子显示了它的作用:http://bl.ocks.org/mbostock/9927735
在我的例子中,生成TopoJSON数据对我来说很容易,这个库函数完美地完成了我的任务。

brgchamk

brgchamk7#

在@Timothy的回答的基础上,我使用以下算法来计算Delaunay三角剖分的边界环。

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

字符串

vwhgwdsa

vwhgwdsa8#

Alpha形状定义为没有超过alpha的边的Delaunay三角剖分。首先删除所有内部三角形,然后删除所有超过alpha的边。

相关问题