Hi,
I am looking for way to work with Voronoi tesselations. I have a set of 2D
points where the stored order is not important, but has to be preserved in
the final result.
I would like to compute the bounded Voronoi tesselation of the points
(i.e., the Voronoi tesselation intersected with a bounding box) and
- compute the corners of each Voronoi cell (which may be on the bounding
box).
- compute the area of each Voronoi cell; the vector of areas should be in
the same order as the Voronoi centers.
(If the corners of the cells are returned in a data structure that is
ordered like the Voronoi centers, the second point is of course obsolete).
The VoronoiDelaunay package computes the corners of each Voronoi cell, but
it would require post-processing to reorder and calculate the intersection
with a bounding box.
I've looked into two alternatives that motivate the second part of the
title:
- Qhull computes the corners of each Voronoi cell and as I recall from
Matlab they are ordered as I need. Unfortunately, I can't even figure out
how to compile Qhull and other posts around here suggest that Qhull is not
easy to work with...
- The R package deldir is a wrapper for a Fortran library and (in R) it
returns exactly what I need. I've written a function that calls the Fortran
library, but Julia segfaults and I can't figure out why.
If I can make it work easily, I think that calling an external library
would be a quick solution here and now. But maybe I have overlooked a
better way?
Best,
Robert
The easiest way to get the shared Fortran library from deldir is probably
to install the R package:
install.package("deldir")
The Julia function that calls the Fortran library gives a few warnings of
the type
WARNING: convert(::Type{Ptr}, ::Int64) methods should be converted to be
methods of unsafe_convert
Changing every integer and Int-type to Int32 doesn't remove the warning.
All the variables needed by the Fortran "master" function are copied from
the R script/wrapper without dwelving into their background. I don't know
what the "master" function returns; first I just want it not to crash.
# Mac
const libdeldir =
"/Library/Frameworks/R.framework/Versions/3.2/Resources/library/deldir/libs/deldir.so"
# Linux
const libdeldir = expanduser(
"~/R/x86_64-pc-linux-gnu-library/3.2/deldir/libs/deldir.so")
function deldir(x::Vector{Float64}, y::Vector{Float64};
sort::Int32=1, rw::Vector{Float64}=[0.0;1.0;0.0;1.0],
epsi::Float64=1e-9 )
@assert (num_points = length(x)) == length(y) "Coordinate vectors must be
of equal length"
@assert epsi >= eps(Float64)
# Dummy points
num_dum_points = 0
npd = num_points + num_dum_points
# The total number of points
ntot = npd + 4
X = [zeros(4); x; zeros(4)]
Y = [zeros(4); y; zeros(4)]
# Set up fixed dimensioning constants
ntdel = 4*npd
ntdir = 3*npd
# Set up dimensioning constants which might need to be increased:
madj = max( 20, ceil(Int32,3*sqrt(ntot)) )
tadj = (madj+1)*(ntot+4)
ndel = Int32( madj*(madj+1)/2 )
tdel = 6*ndel
ndir = ndel
tdir = 8*ndir
nadj = zeros(Int32, tadj)
ind = zeros(Int32, npd)
tx = zeros(Float64, npd)
ty = zeros(Float64, npd)
ilist = zeros(Int32, npd)
delsgs = zeros(Float64, tdel)
delsum = zeros(Float64, ntdel)
dirsgs = zeros(Float64, tdir)
dirsum = zeros(Float64, ntdir)
nerror = Int32[1]
ccall( (:master_, libdeldir), Void,
(Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Float64}, Ptr{Int32},
Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Float64}, Ptr{Float64},
Ptr{Int32}, Ptr{Float64}, Ptr{Float64}, Ptr{Int32}, Ptr{Float64},
Ptr{Float64}, Ptr{Int32}, Ptr{Float64}, Ptr{Int32}),
X, Y, &sort, rw, npd,
&ntot, nadj, &madj, ind, tx, ty,
ilist, &epsi, delsgs, ndel, delsum,
dirsgs, &ndir, dirsum, nerror
)
end
