Feature/qhull

[pmaciel]

This PR adds a ConvexHull base class and an implementation ConvexHullN which is templated to the number of dimensions of calculating a convex hull on. This exists in the eckit::maths library as discussed, and it is disabled by default (enable with cmake -DENABLE_CONVEX_HULL=ON ...)

The convex hull calculation itself is implemented by Qhull, hidden behind the interfacing classes, but exposing enough interface to catch exceptions including those internal to Qhull, so that the returned message(s) is meaningful to client code. Tests (and minimal examples) are in test_convex_hull.cc.

[codecov-commenter]

Codecov Report

Attention: 12 lines in your changes are missing coverage. Please review.

Comparison is base (391e6f4) 63.02% compared to head (19701fe) 63.01%.
Report is 19 commits behind head on develop.

Files	Patch %	Lines
src/eckit/net/Connector.cc	0.00%	12 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           develop      #99      +/-   ##
===========================================
- Coverage    63.02%   63.01%   -0.01%     
===========================================
  Files          792      792              
  Lines        45351    45358       +7     
===========================================
+ Hits         28583    28584       +1     
- Misses       16768    16774       +6     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[tlmquintino]

I've swapped @simondsmart for @danovaro as reviewer given his experience with computational geometry.

[wdeconinck]

See below comments.
I'm wondering if the list_triangles() function should be there at all as it is a special case. Open to discuss.
Atlas already has a bespoke Qhull triangulation API in any case.

[wdeconinck]

list_triangles() only works in particular case I suppose. This should be documented.

[pmaciel]

Yes, but triangulation is rather the common usage in most cases of both convex hull and delaunay triangulations, that's why I've included it. The Triangle element is a "flat" as can be. I didn't go as far as to specify a Polygon element though, I think it's overkill (and there is intersection with the ::geo/::geometry parts too), but it would be correct though.

In random lists of points (in most geometry cases for our applications, but not necessarily the case for polytope (where many of the axes are aligned or of the same size), the result is a triangulation. The common examples using Qhull (that I found) are with "Qt" which enforces this, so, I think it is natural to have it as base functionality

[wdeconinck]

list_triangles() only works in particular case. This should be documented.

[pmaciel]

(as above)

[wdeconinck]

list_triangles() only works in particular case I suppose. This should be documented.

[pmaciel]

(as above)

[pmaciel]

See below comments. I'm wondering if the list_triangles() function should be there at all as it is a special case. Open to discuss. Atlas already has a bespoke Qhull triangulation API in any case.

I'm not sure of the triangulation, but it is rather the common than the special case (hence why I included it). It was in fact the more "testable" situation as in that simpler geometry checks are needed to check correctness.

[wdeconinck]

See below comments. I'm wondering if the list_triangles() function should be there at all as it is a special case. Open to discuss. Atlas already has a bespoke Qhull triangulation API in any case.

I'm not sure of the triangulation, but it is rather the common than the special case (hence why I included it). It was in fact the more "testable" situation as in that simpler geometry checks are needed to check correctness.

What I refer to is that currently list_facets() also returns the triangles, albeit as vector<vector<size_t>> instead of vector<array<size_t,3>>

I would suggest an additional API for list_facets() to avoid the memory fragmentation and the counts of dynamic allocations / deallocations of vector<vector>>. An approach that returns two or three vector:

• two vectors:
  • vector 1: contiguous array that contains all with vertices of the facets
  • vector 2: array with offsets where a facet starts in the above array, size == nb_facets
• three vectors:
  • the two vectors from above, AND a third vector (size == nb_facets) that shows the number of vertices for each facet

The third vector is for the case that you'd want some byte alignment where a new facet starts in vector 1, and the nb_facets cannot be inferred from the difference in offsets in vector 2.

An additional API could wrap these 2/3 vectors in a Facets return type with an inline accessor for a Facet from these vectors

[pmaciel]

Following the recommendation, I decided to follow for the 2-vector solution while removing the list_triangles in favour of a overloaded list_facets(size_t n). For simplicity and to make the return type compatible in the two list_facets methods, the argument n filters the factes for those of the expected size.

This is keeping in mind the considered case (of Polytope?). I agree that for mesh generation one could make the API even more specific, but in that case I would instead pass a single vector, addressable as n*i + j as a block connectivity (as we've used before).

@danovaro @jameshawkes @tlmquintino @wdeconinck what do you think?

[pmaciel]

I've just updated the interface to allow low-level, optimal access to connectivity using a 1-vector solution. So now I see this as a good interface for pure geometrical operations. New methods are:

• facets_n returns a map[facets_vertex_count] = facets_count
• facets returns a std::vector(facets_vertex_count * facets_count) which can be further moved or reshaped, as it is a block connectivity for facets of the same numebr of nodes

Test are modified accordingly.

[wdeconinck]

Thanks for addressing my recommendations.
For my understanding, Is the intention of facets_n map like this?

for( auto& [nb_vertices, nb_facets] : qhull.facets_n() ) {
  auto facets_with_nb_vertices = qhull.facets(nb_vertices);
  ASSERT( facets_with_nb_vertices.size() == nb_vertices * nb_facets );
}

So that for each type of facets you get a separate list?
How would one use this for the square pyramid case?

[wdeconinck]

An additional comment to @jameshawkes , I understood from @tlmquintino that you are currently using CGAL in C++ somewhere, but I can only see python code for polytope. Could you show the use case for this, as this is basically supposed to be used for polytope.

[pmaciel]

Thanks for addressing my recommendations. For my understanding, Is the intention of facets_n map like this?

for( auto& [nb_vertices, nb_facets] : qhull.facets_n() ) {
  auto facets_with_nb_vertices = qhull.facets(nb_vertices);
  ASSERT( facets_with_nb_vertices.size() == nb_vertices * nb_facets );
}

So that for each type of facets you get a separate list? How would one use this for the square pyramid case?

I've covered in the tests specifically the square pyramid case, so that one gets 4 triangles and 1 quadrilateral. I've extended the test to show this functionality at play, showing the usage of both the more generic list_facets() method, and the more specific pair of methods facets_n() and facets(size_t).

[wdeconinck]

Amazing @pmaciel! Thanks for extending the square pyramid test! This is good to merge a.f.a.i.c.

[jameshawkes]

An additional comment to @jameshawkes , I understood from @tlmquintino that you are currently using CGAL in C++ somewhere, but I can only see python code for polytope. Could you show the use case for this, as this is basically supposed to be used for polytope.

We don't currently use CGAL, but as we are building a port of polytope to C++ (currently in Python) we need a convex hull interface we can use. QHull is what we use in Python, via scipy.spatial.

You can see how we currently use it in Python here:
https://github.com/ecmwf/polytope/blob/d4746f36cb70fda7ae6c58d8d24112c48b9056c0/polytope/engine/hullslicer.py#L173C5-L187C56

The interface for passing args in and receiving a result is pretty standard. The main challenge is the error handling. You can see in that code we are switching based on the error description to detect some failure edge cases which are important to us (detecting if the input points are "flat" or if the input shape is lower dimension than the space it occupies). This interface through scipy isn't great.

I discussed a bit with @pmaciel about the need to handle this error-handling a bit better in our C++ wrapper, if possible.

[wdeconinck]

The exception handling is now exceptional! 😅

I think no-one else is probably interested in reviewing this but I'd give some time to the end of the week.
If nobody else reviewed by then, I (or you) can merge this.