# -*- coding: utf-8 -*-

from __future__ import division, print_function

from itertools import permutations
import numpy as np

from collections import OrderedDict


class Triangulation(object):
    """Constrained delaunay triangulation

    Implementation based on [1]_.

    Parameters
    ----------
    pts : array
        Nx2 array of points.
    edges : array
        Nx2 array of edges (dtype=int).

    Notes
    -----
    * Delaunay legalization is not yet implemented. This produces a proper
      triangulation, but adding legalisation would produce fewer thin
      triangles.
    * The pts and edges arrays may be modified.

    References
    ----------
    .. [1] Domiter, V. and Žalik, B. Sweep‐line algorithm for constrained
       Delaunay triangulation


    """
    def __init__(self, pts, edges):
        self.pts = pts[:, :2].astype(np.float32)
        self.edges = edges
        if self.pts.ndim != 2 or self.pts.shape[1] != 2:
            raise TypeError('pts argument must be ndarray of shape (N, 2).')
        if self.edges.ndim != 2 or self.edges.shape[1] != 2:
            raise TypeError('edges argument must be ndarray of shape (N, 2).')

        # described in initialize()
        self._front = None
        self.tris = OrderedDict()
        self._edges_lookup = {}

    def _normalize(self):
        # Clean up data   (not discussed in original publication)

        # (i) Split intersecting edges. Every edge that intersects another
        #     edge or point is split. This extends self.pts and self.edges.
        self._split_intersecting_edges()

        # (ii) Merge identical points. If any two points are found to be equal,
        #      the second is removed and the edge table is updated accordingly.
        self._merge_duplicate_points()

        # (iii) Remove duplicate edges
        # TODO

    def _initialize(self):
        self._normalize()
        # Initialization (sec. 3.3)

        # sort points by y, then x
        flat_shape = self.pts.shape[0] * self.pts.shape[1]
        pts = self.pts.reshape(flat_shape).view([('x', np.float32),
                                                 ('y', np.float32)])
        order = np.argsort(pts, order=('y', 'x'))
        pts = pts[order]
        # update edges to match new point order
        invorder = np.argsort(order)
        self.edges = invorder[self.edges]
        self.pts = pts.view(np.float32).reshape(len(pts), 2)

        # make artificial points P-1 and P-2
        xmax = self.pts[:, 0].max()
        xmin = self.pts[:, 0].min()
        ymax = self.pts[:, 1].max()
        ymin = self.pts[:, 1].min()
        xa = (xmax-xmin) * 0.3
        ya = (ymax-ymin) * 0.3
        p1 = (xmin - xa, ymin - ya)
        p2 = (xmax + xa, ymin - ya)

        # prepend artificial points to point list
        newpts = np.empty((self.pts.shape[0]+2, 2), dtype=float)
        newpts[0] = p1
        newpts[1] = p2
        newpts[2:] = self.pts
        self.pts = newpts
        self.edges += 2

        # find topmost point in each edge
        self._tops = self.edges.max(axis=1)
        self._bottoms = self.edges.min(axis=1)

        # inintialize sweep front
        # values in this list are indexes into self.pts
        self._front = [0, 2, 1]

        # empty triangle list.
        # This will contain [(a, b, c), ...] where a,b,c are indexes into
        # self.pts
        self.tris = OrderedDict()

        # For each triangle, maps (a, b): c
        # This is used to look up the thrid point in a triangle, given any
        # edge. Since each edge has two triangles, they are independently
        # stored as (a, b): c and (b, a): d
        self._edges_lookup = {}

    def triangulate(self):
        """Do the triangulation."""
        self._initialize()

        pts = self.pts
        front = self._front

        # Begin sweep (sec. 3.4)
        for i in range(3, pts.shape[0]):
            pi = pts[i]

            # First, triangulate from front to new point
            # This applies to both "point events" (3.4.1)
            # and "edge events" (3.4.2).

            # get index along front that intersects pts[i]
            idx = 0
            while pts[front[idx+1], 0] <= pi[0]:
                idx += 1
            pl = pts[front[idx]]

            # "(i) middle case"
            if pi[0] > pl[0]:
                # Add a single triangle connecting pi,pl,pr
                self._add_tri(front[idx], front[idx+1], i)
                front.insert(idx+1, i)
            # "(ii) left case"
            else:
                # Add triangles connecting pi,pl,ps and pi,pl,pr
                self._add_tri(front[idx], front[idx+1], i)
                self._add_tri(front[idx-1], front[idx], i)
                front[idx] = i

            # Continue adding triangles to smooth out front
            # (heuristics shown in figs. 9, 10)
            for direction in -1, 1:
                while True:
                    # Find point connected to pi
                    ind0 = front.index(i)
                    ind1 = ind0 + direction
                    ind2 = ind1 + direction
                    if ind2 < 0 or ind2 >= len(front):
                        break

                    # measure angle made with front
                    p1 = pts[front[ind1]]
                    p2 = pts[front[ind2]]
                    err = np.geterr()
                    np.seterr(invalid='ignore')
                    try:
                        angle = np.arccos(self._cosine(pi, p1, p2))
                    finally:
                        np.seterr(**err)

                    # if angle is < pi/2, make new triangle
                    if angle > np.pi/2. or np.isnan(angle):
                        break

                    assert (i != front[ind1] and
                            front[ind1] != front[ind2] and
                            front[ind2] != i)
                    self._add_tri(i, front[ind1], front[ind2])
                    front.pop(ind1)

            # "edge event" (sec. 3.4.2)
            # remove any triangles cut by completed edges and re-fill
            # the holes.
            if i in self._tops:
                for j in self._bottoms[self._tops == i]:
                    # Make sure edge (j, i) is present in mesh
                    # because edge event may have created a new front list
                    self._edge_event(i, j)
                    front = self._front

        self._finalize()

        self.tris = np.array(list(self.tris.keys()), dtype=int)

    def _finalize(self):
        # Finalize (sec. 3.5)

        # (i) Add bordering triangles to fill hull
        front = list(OrderedDict.fromkeys(self._front))

        idx = len(front) - 2
        k = 1
        while k < idx-1:
            # if edges lie in counterclockwise direction, then signed area
            # is positive
            if self._iscounterclockwise(front[k], front[k+1], front[k+2]):
                self._add_tri(front[k], front[k+1], front[k+2])
                front.pop(k+1)
                idx -= 1
                continue
            k += 1

        # (ii) Remove all triangles not inside the hull
        #      (not described in article)

        tris = []  # triangles to check
        tri_state = {}  # 0 for outside, 1 for inside

        # find a starting triangle
        for t in self.tris:
            if 0 in t or 1 in t:
                tri_state[t] = 0
                tris.append(t)
                break

        while tris:
            next_tris = []
            for t in tris:
                v = tri_state[t]
                for i in (0, 1, 2):
                    edge = (t[i], t[(i + 1) % 3])
                    pt = t[(i + 2) % 3]
                    t2 = self._adjacent_tri(edge, pt)
                    if t2 is None:
                        continue
                    t2a = t2[1:3] + t2[0:1]
                    t2b = t2[2:3] + t2[0:2]
                    if t2 in tri_state or t2a in tri_state or t2b in tri_state:
                        continue
                    if self._is_constraining_edge(edge):
                        tri_state[t2] = 1 - v
                    else:
                        tri_state[t2] = v
                    next_tris.append(t2)
            tris = next_tris

        for t, v in tri_state.items():
            if v == 0:
                self._remove_tri(*t)

    def _edge_event(self, i, j):
        """Force edge (i, j) to be present in mesh.

        This works by removing intersected triangles and filling holes up to
        the cutting edge.
        """
        front_index = self._front.index(i)

        front = self._front

        # First just see whether this edge is already present
        # (this is not in the published algorithm)
        if (i, j) in self._edges_lookup or (j, i) in self._edges_lookup:
            return

        # traverse in two different modes:
        #  1. If cutting edge is below front, traverse through triangles. These
        #     must be removed and the resulting hole re-filled. (fig. 12)
        #  2. If cutting edge is above the front, then follow the front until
        #     crossing under again. (fig. 13)
        # We must be able to switch back and forth between these
        # modes (fig. 14)

        # Collect points that draw the open polygons on either side of the
        # cutting edge. Note that our use of 'upper' and 'lower' is not strict;
        # in some cases the two may be swapped.
        upper_polygon = [i]
        lower_polygon = [i]

        # Keep track of which section of the front must be replaced
        # and with what it should be replaced
        front_holes = []  # contains indexes for sections of front to remove

        next_tri = None   # next triangle to cut (already set if in mode 1)
        last_edge = None  # or last triangle edge crossed (if in mode 1)

        # Which direction to traverse front
        front_dir = 1 if self.pts[j][0] > self.pts[i][0] else -1

        # Initialize search state
        if self._edge_below_front((i, j), front_index):
            mode = 1  # follow triangles
            tri = self._find_cut_triangle((i, j))
            last_edge = self._edge_opposite_point(tri, i)
            next_tri = self._adjacent_tri(last_edge, i)
            assert next_tri is not None
            self._remove_tri(*tri)
            # todo: does this work? can we count on last_edge to be clockwise
            # around point i?
            lower_polygon.append(last_edge[1])
            upper_polygon.append(last_edge[0])
        else:
            mode = 2  # follow front

        # Loop until we reach point j
        while True:
            if mode == 1:
                # crossing from one triangle into another
                if j in next_tri:
                    # reached endpoint!
                    # update front / polygons
                    upper_polygon.append(j)
                    lower_polygon.append(j)
                    self._remove_tri(*next_tri)
                    break
                else:
                    # next triangle does not contain the end point; we will
                    # cut one of the two far edges.
                    tri_edges = self._edges_in_tri_except(next_tri, last_edge)

                    # select the edge that is cut
                    last_edge = self._intersected_edge(tri_edges, (i, j))
                    last_tri = next_tri
                    next_tri = self._adjacent_tri(last_edge, last_tri)
                    self._remove_tri(*last_tri)

                    # Crossing an edge adds one point to one of the polygons
                    if lower_polygon[-1] == last_edge[0]:
                        upper_polygon.append(last_edge[1])
                    elif lower_polygon[-1] == last_edge[1]:
                        upper_polygon.append(last_edge[0])
                    elif upper_polygon[-1] == last_edge[0]:
                        lower_polygon.append(last_edge[1])
                    elif upper_polygon[-1] == last_edge[1]:
                        lower_polygon.append(last_edge[0])
                    else:
                        raise RuntimeError("Something went wrong..")

                    # If we crossed the front, go to mode 2
                    x = self._edge_in_front(last_edge)
                    if x >= 0:  # crossing over front
                        mode = 2
                        next_tri = None

                        # where did we cross the front?
                        # nearest to new point
                        front_index = x + (1 if front_dir == -1 else 0)

                        # Select the correct polygon to be lower_polygon
                        # (because mode 2 requires this).
                        # We know that last_edge is in the front, and
                        # front[front_index] is the point _above_ the front.
                        # So if this point is currently the last element in
                        # lower_polygon, then the polys must be swapped.
                        if lower_polygon[-1] == front[front_index]:
                            tmp = lower_polygon, upper_polygon
                            upper_polygon, lower_polygon = tmp
                        else:
                            assert upper_polygon[-1] == front[front_index]

                    else:
                        assert next_tri is not None

            else:  # mode == 2
                # At each iteration, we require:
                #   * front_index is the starting index of the edge _preceding_
                #     the edge that will be handled in this iteration
                #   * lower_polygon is the polygon to which points should be
                #     added while traversing the front

                front_index += front_dir
                next_edge = (front[front_index], front[front_index+front_dir])

                assert front_index >= 0
                if front[front_index] == j:
                    # found endpoint!
                    lower_polygon.append(j)
                    upper_polygon.append(j)
                    break

                # Add point to lower_polygon.
                # The conditional is because there are cases where the
                # point was already added if we just crossed from mode 1.
                if lower_polygon[-1] != front[front_index]:
                    lower_polygon.append(front[front_index])

                front_holes.append(front_index)

                if self._edges_intersect((i, j), next_edge):
                    # crossing over front into triangle
                    mode = 1

                    last_edge = next_edge

                    # we are crossing the front, so this edge only has one
                    # triangle.
                    next_tri = self._tri_from_edge(last_edge)

                    upper_polygon.append(front[front_index+front_dir])

        # (iii) triangluate empty areas

        for polygon in [lower_polygon, upper_polygon]:
            dist = self._distances_from_line((i, j), polygon)
            while len(polygon) > 2:
                ind = np.argmax(dist)
                self._add_tri(polygon[ind], polygon[ind-1],
                              polygon[ind+1])
                polygon.pop(ind)
                dist.pop(ind)

        # update front by removing points in the holes (places where front
        # passes below the cut edge)
        front_holes.sort(reverse=True)
        for i in front_holes:
            front.pop(i)

    def _find_cut_triangle(self, edge):
        """
        Return the triangle that has edge[0] as one of its vertices and is
        bisected by edge.

        Return None if no triangle is found.
        """
        edges = []  # opposite edge for each triangle attached to edge[0]
        for tri in self.tris:
            if edge[0] in tri:
                edges.append(self._edge_opposite_point(tri, edge[0]))

        for oedge in edges:
            o1 = self._orientation(edge, oedge[0])
            o2 = self._orientation(edge, oedge[1])
            if o1 != o2:
                return (edge[0], oedge[0], oedge[1])

        return None

    def _edge_in_front(self, edge):
        """Return the index where *edge* appears in the current front.

        If the edge is not in the front, return -1
        """
        e = (list(edge), list(edge)[::-1])
        for i in range(len(self._front)-1):
            if self._front[i:i+2] in e:
                return i
        return -1

    def _edge_opposite_point(self, tri, i):
        """Given a triangle, return the edge that is opposite point i.

        Vertexes are returned in the same orientation as in tri.
        """
        ind = tri.index(i)
        return (tri[(ind+1) % 3], tri[(ind+2) % 3])

    def _adjacent_tri(self, edge, i):
        """Given a triangle formed by edge and i, return the triangle that shares
        edge. *i* may be either a point or the entire triangle.
        """
        if not np.isscalar(i):
            i = [x for x in i if x not in edge][0]

        try:
            pt1 = self._edges_lookup[edge]
            pt2 = self._edges_lookup[(edge[1], edge[0])]
        except KeyError:
            return None

        if pt1 == i:
            return (edge[1], edge[0], pt2)
        elif pt2 == i:
            return (edge[1], edge[0], pt1)
        else:
            raise RuntimeError("Edge %s and point %d do not form a triangle "
                               "in this mesh." % (edge, i))

    def _tri_from_edge(self, edge):
        """Return the only tri that contains *edge*.

        If two tris share this edge, raise an exception.
        """
        edge = tuple(edge)
        p1 = self._edges_lookup.get(edge, None)
        p2 = self._edges_lookup.get(edge[::-1], None)
        if p1 is None:
            if p2 is None:
                raise RuntimeError("No tris connected to edge %r" % (edge,))
            return edge + (p2,)
        elif p2 is None:
            return edge + (p1,)
        else:
            raise RuntimeError("Two triangles connected to edge %r" % (edge,))

    def _edges_in_tri_except(self, tri, edge):
        """Return the edges in *tri*, excluding *edge*."""
        edges = [(tri[i], tri[(i+1) % 3]) for i in range(3)]
        try:
            edges.remove(tuple(edge))
        except ValueError:
            edges.remove(tuple(edge[::-1]))
        return edges

    def _edge_below_front(self, edge, front_index):
        """Return True if *edge* is below the current front.

        One of the points in *edge* must be _on_ the front, at *front_index*.
        """
        f0 = self._front[front_index-1]
        f1 = self._front[front_index+1]
        return (self._orientation(edge, f0) > 0 and
                self._orientation(edge, f1) < 0)

    def _is_constraining_edge(self, edge):
        mask1 = self.edges == edge[0]
        mask2 = self.edges == edge[1]
        return (np.any(mask1[:, 0] & mask2[:, 1]) or
                np.any(mask2[:, 0] & mask1[:, 1]))

    def _intersected_edge(self, edges, cut_edge):
        """ Given a list of *edges*, return the first that is intersected by
        *cut_edge*.
        """
        for edge in edges:
            if self._edges_intersect(edge, cut_edge):
                return edge

    def _find_edge_intersections(self):
        """Return a dictionary containing, for each edge in self.edges, a list
        of the positions at which the edge should be split.
        """
        edges = self.pts[self.edges]
        cuts = {}  # { edge: [(intercept, point), ...], ... }
        for i in range(edges.shape[0]-1):
            # intersection of edge i onto all others
            int1 = self._intersect_edge_arrays(edges[i:i+1], edges[i+1:])
            # intersection of all edges onto edge i
            int2 = self._intersect_edge_arrays(edges[i+1:], edges[i:i+1])

            # select for pairs that intersect
            err = np.geterr()
            np.seterr(divide='ignore', invalid='ignore')
            try:
                mask1 = (int1 >= 0) & (int1 <= 1)
                mask2 = (int2 >= 0) & (int2 <= 1)
                mask3 = mask1 & mask2  # all intersections
            finally:
                np.seterr(**err)

            # compute points of intersection
            inds = np.argwhere(mask3)[:, 0]
            if len(inds) == 0:
                continue
            h = int2[inds][:, np.newaxis]
            pts = (edges[i, 0][np.newaxis, :] * (1.0 - h) +
                   edges[i, 1][np.newaxis, :] * h)

            # record for all edges the location of cut points
            edge_cuts = cuts.setdefault(i, [])
            for j, ind in enumerate(inds):
                if 0 < int2[ind] < 1:
                    edge_cuts.append((int2[ind], pts[j]))
                if 0 < int1[ind] < 1:
                    other_cuts = cuts.setdefault(ind+i+1, [])
                    other_cuts.append((int1[ind], pts[j]))

        # sort all cut lists by intercept, remove duplicates
        for k, v in cuts.items():
            v.sort(key=lambda x: x[0])
            for i in range(len(v)-2, -1, -1):
                if v[i][0] == v[i+1][0]:
                    v.pop(i+1)
        return cuts

    def _split_intersecting_edges(self):
        # we can do all intersections at once, but this has excessive memory
        # overhead.

        # measure intersection point between all pairs of edges
        all_cuts = self._find_edge_intersections()

        # cut edges at each intersection
        add_pts = []
        add_edges = []
        for edge, cuts in all_cuts.items():
            if len(cuts) == 0:
                continue

            # add new points
            pt_offset = self.pts.shape[0] + len(add_pts)
            new_pts = [x[1] for x in cuts]
            add_pts.extend(new_pts)

            # list of point indexes for all new edges
            pt_indexes = list(range(pt_offset, pt_offset + len(cuts)))
            pt_indexes.append(self.edges[edge, 1])

            # modify original edge
            self.edges[edge, 1] = pt_indexes[0]

            # add new edges
            new_edges = [[pt_indexes[i-1], pt_indexes[i]]
                         for i in range(1, len(pt_indexes))]
            add_edges.extend(new_edges)

        if add_pts:
            add_pts = np.array(add_pts, dtype=self.pts.dtype)
            self.pts = np.append(self.pts, add_pts, axis=0)
        if add_edges:
            add_edges = np.array(add_edges, dtype=self.edges.dtype)
            self.edges = np.append(self.edges, add_edges, axis=0)

    def _merge_duplicate_points(self):
        # generate a list of all pairs (i,j) of identical points
        dups = []
        for i in range(self.pts.shape[0]-1):
            test_pt = self.pts[i:i+1]
            comp_pts = self.pts[i+1:]
            eq = test_pt == comp_pts
            eq = eq[:, 0] & eq[:, 1]
            for j in np.argwhere(eq)[:, 0]:
                dups.append((i, i+1+j))

        dups_arr = np.array(dups)
        # remove duplicate points
        pt_mask = np.ones(self.pts.shape[0], dtype=bool)
        for i, inds in enumerate(dups_arr):
            # remove j from points
            # (note we pull the index from the original dups instead of
            # dups_arr because the indexes in pt_mask do not change)
            pt_mask[dups[i][1]] = False

            i, j = inds

            # rewrite edges to use i instead of j
            self.edges[self.edges == j] = i

            # decrement all point indexes > j
            self.edges[self.edges > j] -= 1
            dups_arr[dups_arr > j] -= 1

        self.pts = self.pts[pt_mask]

        # remove zero-length edges
        mask = self.edges[:, 0] != self.edges[:, 1]
        self.edges = self.edges[mask]

    def _distances_from_line(self, edge, points):
        # Distance of a set of points from a given line
        e1 = self.pts[edge[0]]
        e2 = self.pts[edge[1]]
        distances = []
        for i in points:
            p = self.pts[i]
            proj = self._projection(e1, p, e2)
            distances.append(((p - proj)**2).sum()**0.5)
        assert distances[0] == 0 and distances[-1] == 0
        return distances

    def _projection(self, a, b, c):
        """Return projection of (a,b) onto (a,c)
        Arguments are point locations, not indexes.
        """
        ab = b - a
        ac = c - a
        return a + ((ab*ac).sum() / (ac*ac).sum()) * ac

    def _cosine(self, A, B, C):
        # Cosine of angle ABC
        a = ((C - B)**2).sum()
        b = ((C - A)**2).sum()
        c = ((B - A)**2).sum()
        d = (a + c - b) / ((4 * a * c)**0.5)
        return d

    def _iscounterclockwise(self, a, b, c):
        # Check if the points lie in counter-clockwise order or not
        A = self.pts[a]
        B = self.pts[b]
        C = self.pts[c]
        return np.cross(B-A, C-B) > 0

    def _edges_intersect(self, edge1, edge2):
        """
        Return 1 if edges intersect completely (endpoints excluded)
        """
        h12 = self._intersect_edge_arrays(self.pts[np.array(edge1)],
                                          self.pts[np.array(edge2)])
        h21 = self._intersect_edge_arrays(self.pts[np.array(edge2)],
                                          self.pts[np.array(edge1)])
        err = np.geterr()
        np.seterr(divide='ignore', invalid='ignore')
        try:
            out = (0 < h12 < 1) and (0 < h21 < 1)
        finally:
            np.seterr(**err)
        return out

    def _intersect_edge_arrays(self, lines1, lines2):
        """Return the intercepts of all lines defined in *lines1* as they
        intersect all lines in *lines2*.

        Arguments are of shape (..., 2, 2), where axes are:

        0: number of lines
        1: two points per line
        2: x,y pair per point

        Lines are compared elementwise across the arrays (lines1[i] is compared
        against lines2[i]). If one of the arrays has N=1, then that line is
        compared against all lines in the other array.

        Returns an array of shape (N,) where each value indicates the intercept
        relative to the defined line segment. A value of 0 indicates
        intersection at the first endpoint, and a value of 1 indicates
        intersection at the second endpoint. Values between 1 and 0 are on the
        segment, whereas values outside 1 and 0 are off of the segment.
        """
        # vector for each line in lines1
        l1 = lines1[..., 1, :] - lines1[..., 0, :]
        # vector for each line in lines2
        l2 = lines2[..., 1, :] - lines2[..., 0, :]
        # vector between first point of each line
        diff = lines1[..., 0, :] - lines2[..., 0, :]

        p = l1.copy()[..., ::-1]  # vectors perpendicular to l1
        p[..., 0] *= -1

        f = (l2 * p).sum(axis=-1)  # l2 dot p
        # tempting, but bad idea!
        err = np.geterr()
        np.seterr(divide='ignore', invalid='ignore')
        try:
            h = (diff * p).sum(axis=-1) / f  # diff dot p / f
        finally:
            np.seterr(**err)

        return h

    def _orientation(self, edge, point):
        """ Returns +1 if edge[0]->point is clockwise from edge[0]->edge[1],
        -1 if counterclockwise, and 0 if parallel.
        """
        v1 = self.pts[point] - self.pts[edge[0]]
        v2 = self.pts[edge[1]] - self.pts[edge[0]]
        c = np.cross(v1, v2)  # positive if v1 is CW from v2
        return 1 if c > 0 else (-1 if c < 0 else 0)

    def _add_tri(self, a, b, c):
        # sanity check
        assert a != b and b != c and c != a

        # ignore flat tris
        pa = self.pts[a]
        pb = self.pts[b]
        pc = self.pts[c]
        if np.all(pa == pb) or np.all(pb == pc) or np.all(pc == pa):
            return

        # check this tri is unique
        for t in permutations((a, b, c)):
            if t in self.tris:
                raise Exception("Cannot add %s; already have %s" %
                                ((a, b, c), t))

        # TODO: should add to edges_lookup after legalization??
        if self._iscounterclockwise(a, b, c):
            assert (a, b) not in self._edges_lookup
            assert (b, c) not in self._edges_lookup
            assert (c, a) not in self._edges_lookup
            self._edges_lookup[(a, b)] = c
            self._edges_lookup[(b, c)] = a
            self._edges_lookup[(c, a)] = b
        else:
            assert (b, a) not in self._edges_lookup
            assert (c, b) not in self._edges_lookup
            assert (a, c) not in self._edges_lookup
            self._edges_lookup[(b, a)] = c
            self._edges_lookup[(c, b)] = a
            self._edges_lookup[(a, c)] = b

        tri = (a, b, c)

        self.tris[tri] = None

    def _remove_tri(self, a, b, c):
        for k in permutations((a, b, c)):
            if k in self.tris:
                break
        del self.tris[k]
        (a, b, c) = k

        if self._edges_lookup.get((a, b), -1) == c:
            del self._edges_lookup[(a, b)]
            del self._edges_lookup[(b, c)]
            del self._edges_lookup[(c, a)]
        elif self._edges_lookup.get((b, a), -1) == c:
            del self._edges_lookup[(b, a)]
            del self._edges_lookup[(a, c)]
            del self._edges_lookup[(c, b)]
        else:
            raise RuntimeError("Lost edges_lookup for tri (%d, %d, %d)" %
                               (a, b, c))

        return k


def _triangulate_python(vertices_2d, segments):
    segments = segments.reshape(len(segments) // 2, 2)
    T = Triangulation(vertices_2d, segments)
    T.triangulate()
    vertices_2d = T.pts
    triangles = T.tris.ravel()
    return vertices_2d, triangles


def _triangulate_cpp(vertices_2d, segments):
    import triangle
    T = triangle.triangulate({'vertices': vertices_2d,
                              'segments': segments}, "p")
    vertices_2d = T["vertices"]
    triangles = T["triangles"]
    return vertices_2d, triangles


def triangulate(vertices):
    """Triangulate a set of vertices

    Parameters
    ----------
    vertices : array-like
        The vertices.

    Returns
    -------
    vertices : array-like
        The vertices.
    tringles : array-like
        The triangles.
    """
    n = len(vertices)
    vertices = np.asarray(vertices)
    zmean = vertices[:, 2].mean()
    vertices_2d = vertices[:, :2]
    segments = np.repeat(np.arange(n + 1), 2)[1:-1]
    segments[-2:] = n - 1, 0

    try:
        import triangle  # noqa: F401
    except (ImportError, AssertionError):
        vertices_2d, triangles = _triangulate_python(vertices_2d, segments)
    else:
        segments_2d = segments.reshape((-1, 2))
        vertices_2d, triangles = _triangulate_cpp(vertices_2d, segments_2d)

    vertices = np.empty((len(vertices_2d), 3))
    vertices[:, :2] = vertices_2d
    vertices[:, 2] = zmean
    return vertices, triangles