4

Main problem

In Python, I have triangulated the surface of a (unit) sphere using an icosahedral mesh. I have a list simplices of tuples containing the indices of the three vertices of each triangle, and I have a two lists describing the coordinates (in radians) of each vertex: its latitude and longitude.

For about a million points, I want to determine which triangle each point is in. I am looking for an efficient algorithm that returns the list indices of each triangle (indices corresponding to list simplices).

I am willing to sacrifice memory over efficiency, so I am fine with constructing a tree or using some lookup method.

Caveats

The triangles are of roughly equal size, but not exactly, so I suspect that a simple nearest-neighbor KDTree implementation is not exact.

Extra information

The icosahedral mesh has been obtained using the stripy package. It projects the vertices of the icosahedron onto the unit sphere and subsequently bisects the triangles, so that each edge is split in half, or conversely, each triangle is split in four. stripy has a built-in method for calculating the triangle a point is contained in, but for a mesh refinement of 6 (i.e. 6 bisections) and about a million points, this takes hours. I suspect that this method does not make use of a tree/lookup method and I hope there is a method that significantly improves on this.

falidoro
  • 399
  • 3
  • 14
  • I am not familiar with this but I think you might find some algorithms that are already implemented in some libraries: Qhull, CGAL, ..etc. For example, take a look at the [FAQ of Qhull](http://www.qhull.org/html/qh-faq.htm). [CGAL](https://www.cgal.org/), [binding-swing-cgal](https://github.com/CGAL/cgal-swig-bindings/wiki). Scipy has [spatial module](https://docs.scipy.org/doc/scipy/reference/spatial.html) that uses the `Qhull` library. – s.ouchene Oct 24 '19 at 20:44
  • Do you already have a means of testing containment in a single spherical triangle? – Davis Herring Oct 25 '19 at 02:12
  • In principle, `scipy` can determine the simplex a point is closest too. However, my Delauney triangulation is made on a sphere using `stripy`, and `scipy`'s Delauney triangulation function would return tetrahedrons when used in 3 dimensions. – falidoro Oct 25 '19 at 09:42

1 Answers1

1
  1. Compute a latitude/longitude bounding box for each triangle. Remember that the largest-magnitude latitudes may be on an edge (easily found by considering the normal to the great circle including each edge) or (if the pole is enclosed) in the interior.
  2. Divide all triangles that cross the periodic longitude boundary in two—or, to be cheap, just their bounding boxes.
  3. Build an extended-object k-d tree over the triangles (and triangle pieces from above). This still uses only the latitude/longitude values.
  4. Run the obvious recursive, conservative containment search to find candidate triangles. (It doesn’t matter which piece of the divided triangles you find.)
  5. Test carefully for triangle inclusion: for each triangle side, judge which hemisphere (defined by the great circle containing the segment) contains the query point in a fashion (probably just a cross product on the three-dimensional vectors) that doesn’t depend on the order in which the vertices are presented and never produces “on the dividing line”. Then every point is guaranteed to be in exactly one triangle.
Davis Herring
  • 36,443
  • 4
  • 48
  • 76