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