1

I would like to create a scipy.spatial.Delaunay object from an array of points: P = [[P1x,P1y],[P2x,P2y],...,[PNx,PNy]] and the corresponding triangulation: T = [[T1a,T1b,T1c],[T2a,T2b,T2c],...,[TMa,TMb,TMc]] A similar question has been asked here: How to add simplices to a scipy Delaunay triangulation object . Unfortunately as xdze2 points out, for my case I do need it to be a valid Delaunay object so that I can use find_simplex().

In matlab I can do the following:

TR = triangulation(double(T),double(P)); %generating the triangulation object
[ID,B] = pointLocation(TR,X(:),Y(:)); %X and Y the coordinates of the points for which I want 
                                      %to find the simplex (ID) they lie in as well the 
                                      % barycentric coordinates in that simplex

Can I do something similar with scipy.spatial.Delaunay?

Currently I have written a bruteforce version that works, but it being brute force takes some time to complete: for the same number of queries, my algorithm needs about 180s compared to 18s for the matlab build in function (pointLocation). I would like to see if I can get a speedup by using the find_simplex() method.

The bruteforce method I am using takes advantage of the barycentric transformation as described in the Delaunay doc. For each simplex I compute the barycentric coordinates for each point and find the correct simplex by finding the simplex for which all the coordinates are positive.

for num_simplex in range(T.shape[0]):
    #compute barycentric coordinates for all points in P
    b = transform_container[num_simplex,:2].dot(np.transpose(P-transform_container[num_simplex,2]))
    coords = np.c_[np.transpose(b), 1 - b.sum(axis=0)]
    # find negative coordinates
    neg_vect = coords < 0
    #find positions where ALL coords are positive
    pos_vect = np.sum(neg_vect,axis=1) == 0
   

transform_container is the concatenation of the Delaunay object transform matrix for each simplex. This matrix has been generated by calling Delaunay on subsets of points. These subsets form a partition of the entire set of points P.

My issue is that I cannot call Delaunay on my entire set of points as there are specific nodes that are not connected. That is why I am working with subsets of points to artificially prevent some of the edges.

Thank you in advance for your help

  • What does it mean that some nodes “are not connected”? – Davis Herring Jun 10 '21 at 13:42
  • Certain nodes should not have edges between them. I can actually rephrase my question: Suppose I have N different Delaunay triangulations objects, how can I concatenate them into a single Delaunay object? – Pookelyduke Jun 11 '21 at 04:58
  • It's not clear whether `scipy.spatial.Delaunay` can usefully represent a triangulation that isn't actually Delaunay (_e.g._, because it is a concatenation). It does have an incremental mode, but that seems to preserve the Delaunay property. As such, I would say your question is really just "how can I efficiently find which of my triangles contains a 2D point?", for which the usual answers are _k_-_d_ trees or interval trees. – Davis Herring Jun 11 '21 at 06:36
  • Yes, you are right, that is the ultimate goal. I was wondering if I could use `find_simplex()` . Using your suggestion, I'll look into [scipy.spatial.KDTree](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.KDTree.html), maybe that'll get me started – Pookelyduke Jun 11 '21 at 07:57
  • update: I forgot to mention, I do also need the barycentric transformation that the Delaunay object carries, hence my initial question of concatenation – Pookelyduke Jun 11 '21 at 11:16
  • I wanted to continue my original question. David Herring's comment does not reply to the question posed. How can I generate a [scipy.spatial.Delaunay](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.Delaunay.html) object from a list of points and the corresponding array of edges? – Pookelyduke Jul 07 '21 at 06:26
  • Does anyone else have a suggestion? – Pookelyduke Sep 03 '21 at 05:35

0 Answers0