Calling C++ function that takes and returns Eigen Matrix from Julia

New to Julia
,
1

What is the best approach to call a C++ function which accepts and returns an Eigen Matrix? The function has signature

Eigen::Matrix3Xi makeconvexhull(Eigen::Matrix3Xd points)

Here, points is a 3 x N matrix of double and the function returns a M x 3 matrix of integers where every row represents a triangle (i.e. vertex ids of the triangle corresponding to coordinates in points).

  1. I found Polyhedra.jl which seems to provide a native Julia code to make convex hulls out of point clouds but I could not understand how to use it.
  2. I already have a C++ code that uses CGAL and Eigen and is specialized for my use case and thus very performant. So I would like to use it preferrably.
  3. My current approach is to use ccall which requires me to write a wrapper function in my C++ library that gives a C-style interface to my function.
extern "C"{
        struct Triangles{
               int* firsttriangle;
               size_t numtri;
        };

        Triangles wrapper(size_t N, double* points){
             // Use the pointer and size information to make an Eigen matrix
            //  and forward to the existing function
                ...
            // Take the output of my existing function and convert it to struct Triangles and return
            Trianlges triangles;
            triangles.firsttriangle = < address of first element of triangle matrix >;
            triangles.numtri = < number of triangles >;
            return triangles;
        } 
}

Is there a better approach that this?
4. In this approach, I have to do a malloc() in my wrapper function to allocate memory for the output matrix. How do I free this memory from my Julia code once it is not required any more?

2

My main comment is perhaps not what you were expecting, but that I suspect your best bet would be to reach out to the developers of Polyhedra.jl to clarify the things you need. If you volunteer to contribute to the documentation in exchange for some free consultation, I suspect it will be well received. Personally, I always find it difficult to write good documentation when I’m also the one writing the code, because I know too much about the internals to be able to approach it from the perspective of a newbie. The main developer, @blegat, is a very active & capable contributor.

To your specific questions,

6 Likes
3

Personally, I always find it difficult to write good documentation when I’m also the one writing the code, because I know too much about the internals to be able to approach it from the perspective of a newbie.

Exactly, I would be interested to hear what is confusing in the documentation. In your case, it seems you are interested to find the V-representation of the facets of a 3D polytope given its V-representation.
One way to do it is to use incidence information:
https://juliapolyhedra.github.io/Polyhedra.jl/stable/polyhedron/#Incidence-1
Given an halfspace of the H-representation, all the points incident to the halfspace give the V-representation of the facet exposed by the halfpsace.
The incidence information is currently only implemented by CDDLib but it should be possible to implemented also by QHull (by the way, this is the list of V-representation of the facets).

2 Likes
4

Thank you @tim.holy and @blegat for your time and answers.

  1. I am sorry that I did not express properly. I could not understand the documentation because I am not familiar with H-representation, V-representation terminology and I did not spend enough time trying to understand it. I was looking for a function that would take points as input and give triangulation as output.
  2. Qhull has such a function. I have used it before through Python. My CGAL implementation is faster in my measurements (I need to do several millions of triangulation for small point clouds) and I have found some corner cases in which Qhull gives fewer facets because it uses floating point arithmetic compared to CGAL’s exact predicates. That’s why I did not use Qhull.
  3. I came across CxxWrap.jl but I could not find any information if it can accept Eigen matrices as arguments. Perhaps I should reach out the developer.
  4. I found that I can decide in advance how much memory to allocate so I settled on passing a pre-allocated matrix to my C-function.

Thanks!

5

Here is an example to obtain a triangulation given a 3 x N Matrix using Polyhedra

using Polyhedra, CDDLib
A = randn(3, N)

# use floating point for the computation
P = polyhedron(vrep(A'), CDDLib.Library()) # polyhedra expects each point to be a *row*
triang = incidentpoints.(P, eachindex(halfspaces(P)))

# using exact arithmetic
Q = polyhedron(vrep(A'), CDDLib.Library(:exact))
triang_exact = incidentpoints.(Q, eachindex(halfspaces(Q)))
4 Likes