context:
I wrote a function to "uniformly" distribute points on the surface of a sphere/globe... I put a point at each pole, and then divided the rest of the east north space into latitude bands where the width/height of the latitude band is proportional to the number of points in it and the number of latitude bands and points in each latitude band is chosen to make make the east-west spacing roughly equal to the north-south spacing, and then I rotate each band independently to maximize the minimum longitude difference of any two points in consecutive latitude bands, and then I do a single spherical CVT (Centroidal Voronoi Tessellation) step to slightly improve the spacing of points. To do the CVT I use convhulln to get the delaunay tessellation of the surface of the sphere and compute the 1 Voronoi cell vertex per triangle, fractional area of each triangle that belongs to each of its nodes along with the centroid of the fractional areas, which I can add up (weighted by the fractional areas) to get the updated point positions. (the first CVT step improves the spacing then it decreases) anyway I'm getting a min/max point area ratio of about .91 which isn't a minimum energy/maxi-min design but it is quite close and it takes a second or less to generate. Using the point areas as weights in an integration also makes it a quite accurate at integrating a function over a sphere (getting between 3 and 5 sig figs of accuracy for O(1000) points, when the standard deviations equals the mean for the test function with an analytical solution).
That's context... put another way I DON'T have a rectangular grid of points on an east-west north-south map.
I would like to unwrap the function from the sphere onto an east north map, and I figured out the logic for a piece wise (triangle wise) constant unwrapping/flattening of sphere (which was the most complicated part of all of this, much more complicated than the point generation) and this meant that edges of the triangles would come out as curves on the map.
But really it would look a lot better (and be simpler) if I took an TP grid of lat long points on the map, translated them to xyz (on the surface of the sphere) points, identified the triangles on the sphere that contained those xyz points, did barycentric interpolation there, and then used those values with the "original" Tensor Product lat long coordinates to display the map.
The part I'm missing is a tsearchn equivalent for identifying which triangles on the sphere contain the xyz points; I am not a qhull expert, so what I would write would be no-where nearly as fast as the tsearchn function is. I was hoping that someone out there had already solved this.
Back when tsearchn was slow (took O(3) minutes for a few thousand points in O(6) dimensions), I wrote a replacement that took about a minute by using a space-filling curve as a short cut to restrict the number of triangles that I needed to check, and then using a simplex/neighbor based search, but shortly thereafter MATHWORKS released a update based on qhull that managed it in a few seconds, so like I said a qhull based (or similar) approach is going to be about 2 orders of magnitude faster than what I would write using my SFC+simplex approach.
