0

我想使用 Julia 计算一组点的 alpha 形状（甚至只是凹壳）。在其他问题中，他们通过使用包含给定点集的 Delaunay tesselation Boundary在 python 中解决了这个问题。

Julia 中的这个包可以获得 Delaunay 镶嵌 https://github.com/JuliaGeometry/VoronoiDelaunay.jl（虽然我不确定它是否针对 julia v0.7 进行了更新）。我想知道 julia v0.7 是否已经有一个实现可以得到 eh alpha 形状，甚至只是一组点的凹壳。

或者，有没有办法有效地调用 python (scipy.spatial.Delaunay) 来完成这项工作？

4

1 回答 1

2

VoronoiDelaunay.jl适用于 Julia 1.0 和 1.1。它也应该适用于 Julia 0.7。

VoronoiDelaunay.jl有一些数值限制，即(1.0+eps(), 2.0-eps())坐标，因此您可能需要重新缩放数据点。

要DelaunayTesselation使用您自己的点类型创建一个，请确保您的类型是 的子类型AbstractPoint2D，即<: AbstractPoint2D和 定义getx和gety方法。

我相信，下面的示例代码使用一组点找到了所谓的凹壳，DelaunayTesselation并绘制了结果。它在这个答案中基本上使用了相同的算法。您可以轻松地编辑代码以获得 alpha 形状。

我没有将一些代码片段包装到一个函数中。如果您需要高性能，请这样做。我===在检查点是否相等时使用，它实际上检查两个点是否是同一个对象（即内存中的地址）。如果你以某种方式最终得到一个破坏这部分的代码，你可以扩展==和使用它而不是===.

using Random, VoronoiDelaunay, Plots

import Base.==

struct MyEdge{T<:AbstractPoint2D}
    _a::T
    _b::T
end

==(e1::MyEdge{T}, e2::MyEdge{T}) where {T<:AbstractPoint2D} = ((e1._a === e2._a) && (e1._b === e2._b)) || ((e1._b === e2._a) && (e2._b === e1._a))

###==(p1::T, p2::T) where {T<:AbstractPoint2D} = (getx(p1) == getx(p2)) && (gety(p1) == gety(p2))

### Create a Delaunay tesselation from random points
tess = DelaunayTessellation2D(46)

for _ in 1:23
    push!(tess, Point2D(rand()+1, rand()+1))
end

edges = MyEdge[]

function add_edge!(edges, edge)
    i = findfirst(e -> e == edge, edges)

    if isnothing(i) # not found
        push!(edges, edge)
    else # found so not outer, remove this edge
        deleteat!(edges, i) 
    end
end

for trig in tess
    a, b, c = geta(trig), getb(trig), getc(trig)
    add_edge!(edges, MyEdge(b, c))
    add_edge!(edges, MyEdge(a, b))
    add_edge!(edges, MyEdge(a, c))
end

### PLOT
x, y = Float64[], Float64[] # outer edges
for edge in edges
    push!(x, getx(edge._a))
    push!(x, getx(edge._b))
    push!(x, NaN)
    push!(y, gety(edge._a))
    push!(y, gety(edge._b))
    push!(y, NaN)
end

xall, yall = getplotxy(delaunayedges(tess)) # all the edges

plot(xall, yall, color=:blue, fmt=:svg, size=(400,400))
plot!(x, y, color=:red, linewidth=3, opacity=0.5)
于 2019-03-11T22:14:46.613 回答