Contour from triangulation in Python I have a set of 2d vertices (x,y) and a set of triangles referring to the vertices' index (p1,p2,p3).

The vertices represent a closed polygon shape, however those are not in order.

When 'merging' all the triangles together, we get a the area of the desired polygon (from which we can conclude the contour).

I've manage to achieve this with the following method (multiple polygon union), however not sure about its efficiency. Any advice?

My database is 1M+ entries.

Thanks

Nb this question has been forwarded from here

from shapely.geometry import Polygon
from shapely.ops import cascaded_union

def get_region(vertices, triangles):
    region = Polygon()

    for i, tri in enumerate(triangles):
        a = vertices[tri[0]]
        b = vertices[tri[1]]
        c = vertices[tri[2]]
        p = Polygon ([a, b, c])

        try:
            region = cascaded_union([region, p])
        except ValueError:
            print ('value error at: ', i)

    return (list(region.exterior.coords))

region = get_region(vertices, triangles)




Vertices:

[(2075, 2811),
 (1992, 2790),
 (1992, 2759),
 (3457, 2527),
 (2459, 2146),
 (3460, 2492),
 (3463, 2435),
 (3465, 2321),
 (2312, 2209),
 (1999, 2622),
 (2178, 2283),
 (2423, 2154),
 (2443, 2150),
 (1998, 2613),
 (1987, 2555),
 (2002, 2383),
 (1994, 2398),
 (1999, 2387),
 (1993, 2606),
 (1983, 2520),
 (1982, 2486),
 (2080, 2345),
 (2045, 2362),
 (1992, 2403),
 (1985, 2434),
 (2063, 2355),
 (1995, 2392),
 (3458, 2206),
 (3025, 1983),
 (3207, 1991),
 (3459, 2200),
 (3461, 2191),
 (3468, 2157),
 (3310, 1997),
 (3470, 2101),
 (3385, 2009),
 (2832, 1976),
 (2568, 2058),
 (2648, 2005),
 (2507, 2101),
 (2461, 2142),
 (2461, 2138),
 (2702, 1975),
 (2721, 1976),
 (2680, 1979),
 (2692, 1976),
 (2713, 1975),
 (2774, 1973),
 (2896, 1975),
 (3463, 2063),
 (3469, 2074),
 (3457, 2049),
 (3439, 2025),
 (3458, 2044),
 (3451, 2029),
 (3458, 2039),
 (3456, 2033),
 (3424, 2021),
 (3460, 2056),
 (3465, 2179),
 (3459, 2210),
 (3461, 2219),
 (3466, 2249),
 (3437, 2662),
 (3434, 2680),
 (3418, 2754),
 (3405, 2787),
 (1999, 2631),
 (3445, 2610),
 (3399, 2804),
 (3398, 2811),
 (1995, 2812),
 (1994, 2705),
 (3433, 2689),
 (2071, 2350),
 (1999, 2630),
 (1999, 2631)]


Triangles:

[(0, 1, 2),
 (3, 4, 5),
 (6, 4, 7),
 (8, 9, 10),
 (11, 12, 8),
 (13, 14, 10),
 (15, 16, 17),
 (14, 13, 18),
 (19, 20, 21),
 (22, 16, 15),
 (23, 16, 22),
 (23, 22, 24),
 (25, 20, 24),
 (26, 17, 16),
 (14, 19, 21),
 (14, 21, 10),
 (13, 10, 9),
 (12, 4, 8),
 (27, 28, 29),
 (30, 29, 31),
 (32, 33, 34),
 (33, 35, 34),
 (36, 37, 38),
 (39, 4, 40),
 (39, 40, 41),
 (42, 43, 44),
 (43, 38, 44),
 (44, 45, 42),
 (43, 42, 46),
 (47, 38, 43),
 (36, 38, 47),
 (48, 37, 36),
 (39, 37, 48),
 (39, 48, 28),
 (39, 28, 4),
 (34, 49, 50),
 (51, 52, 53),
 (52, 54, 55),
 (54, 56, 55),
 (52, 55, 53),
 (57, 52, 51),
 (57, 51, 35),
 (51, 58, 35),
 (58, 49, 35),
 (35, 49, 34),
 (29, 33, 32),
 (29, 32, 59),
 (29, 59, 31),
 (27, 29, 30),
 (60, 28, 27),
 (61, 28, 60),
 (62, 28, 61),
 (7, 28, 62),
 (4, 28, 7),
 (5, 4, 6),
 (0, 8, 63),
 (9, 8, 0),
 (63, 64, 0),
 (65, 66, 0),
 (9, 0, 67),
 (68, 63, 8),
 (69, 0, 66),
 (0, 69, 70),
 (0, 71, 1),
 (72, 0, 2),
 (65, 0, 73),
 (73, 0, 64),
 (3, 68, 4),
 (4, 68, 8),
 (74, 21, 25),
 (22, 25, 24),
 (21, 20, 25),
 (9, 67, 75),
 (67, 76, 75),
 (72, 67, 0)]
The cascaded_union function can take more than one object, so you should be able to do:

from shapely.geometry import Polygon, MultiPolygon
from shapely.ops import cascaded_union

def get_region(vertices, triangles):
    regions = MultiPolygon([Polygon([vertices[i] for i in triangle])
                            for triangle in triangles])
    return list(cascaded_union(regions).exterior.coords)


When deleting triangle[72] from your example, which is not really a triangle and the reason for your try...except ValueError, this gives the same area for this function and yours, 1060307.0.

If you cannot delete those zero-size triangles manually, you can do it programmatically as well:

def get_region(vertices, triangles):
    polygons = (Polygon([vertices[i] for i in triangle])
                for triangle in triangles)
    regions = MultiPolygon([p for p in polygons if p.area &gt; 0])
    return list(cascaded_union(regions).exterior.coords)


I would leave getting the points of the exterior to the caller, or another function, but I left it here in order to not change the interface.

Note that cascaded_union may be superseded by unary_union, depending on the version you are using.

I am not sure why you use the datastructure you have. It would be easier to directly store the vertices for each triangle as a three-tuple, so that you could just do:

def get_region(triangles):
    polygons = (Polygon(triangle) for triangle in triangles)
    regions = MultiPolygon([p for p in polygons if p.area &gt; 0])
    return list(cascaded_union(regions).exterior.coords)