0

My question is simple, but the solution might be very tricky. I have a collection of triangles, and I want to find their union. The triangles are given in the standard way: by a list of coordinate points (2 coordinates for each point) and a list of connections, where each line is a triangle given by the indices of its vertices:

points:

[ 15.02716923  81.72425842]
[ 21.42242702  79.91459549]
[ 24.87068939  79.0222168 ]
[ 29.25767326  77.96657562]
[ 34.07667923  76.65890503]

triangles:

7 8 9
8 18 20
8 20 10
8 10 9
9 10 11
9 11 108

I can very efficiently find the union using cascaded_union:

import shapely.geometry as geometry
from shapely.ops import cascaded_union, polygonize

polys = [geometry.Polygon([[points[point, 0], points[point, 1]] for point in triangle]) for triangle in triangles]
result = cascaded_union (polys)

The problem is, this gives the coordinates of the resulting polygon. Whereas, I would like the circumference of result as indices in the initial points array. So far I have not found a way to do this. One way may be to write my own combining function which spits out the point indices instead of the coordinates themselves.

John
  • 465
  • 1
  • 6
  • 13
  • Does this mean that you want the convex hull? – Bill Bell Mar 29 '17 at 13:21
  • @BillBell No, I need the concave hull. However, I do not need to find alpha shapes, since I know the connection triangles already. So I only need to merge them and produce the resulting shape. – John Mar 29 '17 at 16:26
  • what do you expect to get from the points of the exterior ring of the union that are not in your list? (this can easily happen). – eguaio Mar 29 '17 at 18:06
  • @eguaio Not sure what your question is exactly, but here is what I need: at the moment, my code above gives the exterior ring as a series of coordinates: (x1 y1, x2 y2, ..., xn yn). Instead, I want it in terms of point indices, where the indices are those which were used to denote the triangles. So the exterior ring would be given by (p1, p2, ..., pn), where pi would be an index in the `points` array, and point pi would have coordinates xi, yi. However, `cascaded_union` can only return a series of coordinates, and not a series of point indices (or point objects, for that matter). – John Mar 29 '17 at 23:30
  • I will illustrate what I mean with an example. If you make the union of the two triangles [(0,0),(1,1),(1,0)] and [(0,0), (0,1),(1,0)] the polygon you get is given by the coordinates: [(0,0),(0,1),(sqrt(2)/2,sqrt(2)/2),(1,1),(1,0)]. The point (sqrt(2)/2, sqrt(2)/2) is not a point in any of the original triangles (because it is a new point formed by two intersecting lines). The question is what index do you expect to get in those points that are not from the original triangles. I can imagine examples where there is a sequence of many not-original points. – eguaio Mar 30 '17 at 11:44
  • @eguaio No, in fact this should not happen with my triangles, as they are part of a triangle grid. – John Mar 30 '17 at 12:57
  • 1
    Do you have a definition of 'triangle grid' that can be used to generate a set of random triangles? – Bill Bell Mar 30 '17 at 13:09
  • I don't think `shapely` has any functionality to cover this specific case. Besides, in that conditions you can build an ad-hoc algorithm as you suggested that will beat in performance the shapely union, because no geometry is really needed. The only condition that I think could simplify the algorithm is if the triangles in the list are unique and no overlapping other than in the edges (which is true if my interpretation of 'grid' is correct). – eguaio Mar 30 '17 at 17:45

1 Answers1

0

If the union always return polygons with points from your list, you can just lookup the indexes. Build a dictionary where a tuple of coordinates point to the index you use (named coord_to_pointer in the code), then you can get the indexes easily:

polys = [geometry.Polygon([[points[point, 0], points[point, 1]] 
                          for point in triangle]) for triangle in triangles]
result = cascaded_union (polys)
geom_indexes = [coord_to_pointer[p] for p in list(result.exterior.coords)]
eguaio
  • 3,754
  • 1
  • 24
  • 38