triangulation

# -*- coding: utf8 -*- 
"""
Constrained delaunay implementation based on 

Domiter, V. and Žalik, B.
Sweep‐line algorithm for constrained Delaunay triangulation

(this implementation is not complete)
"""

import numpy as np

pts = [(0, 0),
       (10, 0),
       (10, 10),
       (20, 10),
       (20, 20),
       (10, 17),
       (9, 30),
       (5, 25),
       (6, 15),
       (10, 12),
       (0, 5)]

pts = np.array(pts, dtype=[('x', float), ('y', float)])
edges = np.array([(i, (i+1) % pts.shape[0]) for i in range(pts.shape[0])], dtype=int)

## Initialization (sec. 3.3)

# sort by y, then x
order = np.argsort(pts, order=('y', 'x'))  
pts = pts[order]

# update edges to match new point order
invorder = np.argsort(order)
edges = invorder[edges]

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

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

# find topmost point in each edge
tops = set(edges.max(axis=1))

# inintialize sweep front
front = [0, 2, 1]
tris = []


# visual debugging: draw edges, front, triangles
import pyqtgraph as pg
import time
app = pg.mkQApp()
win = pg.plot()
gpts = pts.view(float).reshape(len(pts), 2)
graph = pg.GraphItem(pos=gpts, adj=edges, pen={'width': 3, 'color':(0,100,0)})
win.addItem(graph)
front_line = pg.PlotCurveItem(pen={'width': 2, 'dash': [5, 5], 'color': 'y'})
win.addItem(front_line)
tri_shapes = []
def draw_state():
    front_pts = pts[np.array(front)]
    front_line.setData(front_pts['x'], front_pts['y'])
    for i in range(100):  # sleep ~1 sec, but keep ui responsive
        app.processEvents()
        time.sleep(0.01)
def draw_tri(tri):
    tpts = pts[np.array(tri)]
    path = pg.arrayToQPath(tpts['x'], tpts['y'])
    shape = pg.QtGui.QGraphicsPathItem(path)
    shape.setPen(pg.mkPen(255,255,255,100))
    shape.setBrush(pg.mkBrush(0,255,255,50))
    win.addItem(shape)
    tri_shapes.append(shape)
draw_state()


## Sweep (sec. 3.4)

def add_tri(f0, f1, p):
    # todo: legalize!
    global front, tris
    tris.append((front[f0], front[f1], p))
    front.insert(f1, p)
    
    # visual debugging
    draw_tri(tris[-1])

for i in range(3, len(pts)):
    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]
    l = 0
    while pts[front[l+1]]['x'] <= pi['x']:
        l += 1
    pl = pts[front[l]]
    pr = pts[front[l+1]]
    if pi['x'] > pl['x']:  # "(i) middle case"
        # Add a single triangle connecting pi,pl,pr
        add_tri(l, l+1, i)
    else: # "(ii) left case"
        ps = pts[l-1]
        # Add triangles connecting pi,pl,ps and pi,pl,pr
        add_tri(l, l+1, i)
        add_tri(l-1, l, i)
        front.pop(l+1)
        
    # todo: Continue adding triangles to smooth out front
    # (use heuristics shown in figs. 9, 10)
                    
    if i in tops: # this is an "edge event" (sec. 3.4.2)
        pass
        # (i) locate intersected triangles
        # (ii) remove intersected triangles
        # (iii) triangluate empty areas

    draw_state()

## Finalize (sec. 3.5)

# (i) Remove all triangles that include at least one artificial point
# (ii) Add bordering triangles to fill hull