Voronoi tess / ccall with Fortran

[Robert DJ]

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

[Tom Breloff]

I don't have a complete answer, but wanted to make you aware of this package: https://github.com/Voxel8/Voronoi.jl

[Lutfullah Tomak]

I think you already know Fortran functions take arguments by reference/pointer. If there is not a typo in ccall some of the arguments are not Ptr but IntXX in the example. I suggest you use Ref instead of Ptr and just pass any variable since it will be converted automatically to a pointer. You can use Ref if you are in julia 0.4+. Warning suggest you use unsafe_convert since casting an integer a pointer is unsafe(I think).
Reference to ccall and Ref
http://docs.julialang.org/en/release-0.4/manual/calling-c-and-fortran-code/

[Robert DJ]

Actually, I wasn't aware of that package -- thanks for the link! It looks like the output from this package would also need some work, though.

[Robert DJ]

I'm no next to nothing about Fortran -- I was mostly inspired by the package https://github.com/simonster/GLMNet.jl/

I'll look into Ptr/Ref !

[Robert DJ]

I ended up wrapping the Fortran code from the R package. If it's of interest I'm sharing the code:

https://github.com/robertdj/Deldir.jl

[Tony Kelman]

Note that if you're redistributing GPL Fortran code, you should include the text of the GPL in your license file (with a note, as you have in the readme, that says it's applicable to the Fortran code).

[Tom Breloff]

Thanks Robert... it is of interest :)

I'm curious... is this really allowed?

The Julia code in this package is MIT licensed and the Fortran code is licensed under GPL.

Can I call this from another package that is fully MIT (and without mentioning the license), or does that GPL notice need to be attached to everything up the dependency chain?

Also... if possible it would be preferable to remove the dependency on DataFrames (but I understand that's a tough decision to make after the fact).

Thanks for the contribution.

Tom

[Stefan Karpinski]

Yes, it is completely allowed. As long as you do not violate the GPL license of the code you are creating a derived work with – and the releasing your code under the MIT license satisfies the GPL, so it's fine. What you cannot do is release a derived work that is or uses proprietary code since that would mean that you are not abiding by the terms of the GPL and your right to use or redistribute the GPL code would be nullified.

[Tom Breloff]

Just to be super clear... if I was to use this package from within Plots, would I need to change anything about my license file?

[Stefan Karpinski]

Not legally, but it's best to be clear that the combined work is GPL.