def find_local_boundary(tri, triangles):
    """Find and return the outside edges of a collection of natural neighbor triangles.

    There is no guarantee that this boundary is convex, so ConvexHull is not
    sufficient in some situations.

    Parameters
    ----------
    tri: `scipy.spatial.Delaunay`
        A Delaunay Triangulation
    triangles: (N, ) array
        List of natural neighbor triangles.

    Returns
    -------
    edges: list
        List of vertex codes that form outer edges of
        a group of natural neighbor triangles.

    """
    edges = []
    for triangle in triangles:
        for i in range(3):
            pt1 = tri.simplices[triangle][i]
            pt2 = tri.simplices[triangle][(i + 1) % 3]
            if (pt1, pt2) in edges:
                edges.remove((pt1, pt2))
            elif (pt2, pt1) in edges:
                edges.remove((pt2, pt1))
            else:
                edges.append((pt1, pt2))
    return edges
