Iteratively polygonize MultiLineString

I have an implementation of alpha-shape below. My goal is to compute the polygons for a list of alphas but do it as fast as possible.

LineProfiler showed that almost all of the time is spent in list(polygonize(m)). Since each alpha (in descending order) would add edges I figured I could polygonize in an iterative fashion to save time.

def alpha_shape(points, alphas):
    """
    Compute the alpha shape (concave hull) of a set
    of points.
    @param points: Iterable container of points.
    @param alpha: alpha value to influence the
        gooeyness of the border. Smaller numbers
        don't fall inward as much as larger numbers.
        Too large, and you lose everything!
    """

    if len(points) < 4:
        # When you have a triangle, there is no sense
        # in computing an alpha shape.
        return geometry.MultiPoint(list(points)).convex_hull
    coords = np.array([point.coords[0] for point in points])
    tri = Delaunay(coords)
    og_triangles = coords[tri.vertices]
    a = ((og_triangles[:, 0, 0] - og_triangles[:, 1, 0]) ** 2 + (og_triangles[:, 0, 1] - og_triangles[:, 1, 1]) ** 2) ** 0.5
    b = ((og_triangles[:, 1, 0] - og_triangles[:, 2, 0]) ** 2 + (og_triangles[:, 1, 1] - og_triangles[:, 2, 1]) ** 2) ** 0.5
    c = ((og_triangles[:, 2, 0] - og_triangles[:, 0, 0]) ** 2 + (og_triangles[:, 2, 1] - og_triangles[:, 0, 1]) ** 2) ** 0.5
    s = (a + b + c) / 2.0
    areas = (s*(s-a)*(s-b)*(s-c)) ** 0.5
    circums = a * b * c / (4.0 * areas)

    alpha_polygons = {}
    alphas.sort(reverse=True)
    for alpha in alphas:
        filtered = og_triangles[circums < (1.0 / alpha)]
        edge1 = filtered[:, (0, 1)]
        edge2 = filtered[:, (1, 2)]
        edge3 = filtered[:, (2, 0)]
        edge_points = np.unique(np.concatenate((edge1, edge2, edge3)), axis=0).tolist()
        m = geometry.MultiLineString(edge_points)

        triangles = list(polygonize(m))
        polygons = cascaded_union(triangles)
        if type(polygons) == shapely.geometry.Polygon:
            # In case it's just one polygon, this puts it in a list anyway
            alpha_polygons[alpha] = [polygons]
        else:
            alpha_polygons[alpha] = polygons
    return alpha_polygons

If I do something naive like the following it generates a few more polygons than if you do it all at once together.

list(polygonize(m[:10000])) + list(polygonize(m[10000:]))

I think what we’re basically trying to do is polygonize a set of lines, and then add lines from there but end up with the same result as if we had polygonized all the lines at once. Is this possible?