/Triangulate.hx

package;

import ogl.GLM;

typedef P2 = {x:Float,y:Float};
typedef PV = {x:Float,y:Float,u:Float,cnt:Int};
typedef TValues = {e0:Float,e1:Float,e2:Float};
///
/// Triangulate.triangulate(fs:Array<P2>):ts:Array<Int>
///     ts[i,i+1,i+2]  (i%3) triangle indices
///
/// Triangulate.divide(fs:Array<P2>,ts:Array<Int>,vs:Array<TValues>):{fs:Array<PV>,ts:Array<Int>}
///     each triangle turned into 7, with 6 extra vertices.
///

class Triangulate {
    static inline var epsilon = 1e-12;

    static inline function h(a:Float,b:Float) return 0.5*(a+b);
    public static function divide(fs:Array<P2>,ts:Array<Int>,vs:Array<TValues>,skipTs=false):{fs:Array<PV>,ts:Null<Array<Int>>} {
        var nf = fs.length;
        var nt = Std.int(ts.length/3);

        // fs[i] -> ofs[i]
        // ed[i] -> ofs[#fs + i] for 1 vertex
        var ofs:Array<PV> = [];
        for (f in fs) ofs.push({x:f.x,y:f.y,u:0.0,cnt:0});

        function gen(a:P2,b:P2) return
            { x:h(a.x,b.x),
              y:h(a.y,b.y),
              u:0.0,
              cnt:0
            };

        // Generate edge vertices
        // (i,j) -> (i*#fs + j)
        var ed = new Map<Int,Int>();
        for (i in 0...nt) {
            var t = i*3;
            var t0 = ts[t+0];
            var t1 = ts[t+1];
            var t2 = ts[t+2];
            function insert(a, b) {
                var ind = a*nf + b;
                if (!ed.exists(ind)) {
                    var jnd = ed[ind] = ofs.length;
                    ofs.push(gen(fs[a],fs[b]));
                }
            }
            insert(t1,t2);
            insert(t2,t0);
            insert(t0,t1);
        }

        // And non-normalised colours
        // And new triangles.
        var ots:Array<Int> = if (skipTs) null else [];
        for (i in 0...Std.int(ts.length/3)) {
            var val = vs[i];
            var t = i*3;
            var t0 = ts[t+0];
            var t1 = ts[t+1];
            var t2 = ts[t+2];
            var u0 = ed[t1*nf+t2];
            var u1 = ed[t2*nf+t0];
            var u2 = ed[t0*nf+t1];

            var ind = ofs.length;
            var f0 = ofs[t0]; f0.cnt++; f0.u += val.e1+val.e2-val.e0;
            var f1 = ofs[t1]; f1.cnt++; f1.u += val.e2+val.e0-val.e1;
            var f2 = ofs[t2]; f2.cnt++; f2.u += val.e0+val.e1-val.e2;
            var e0 = ofs[ed[t1*nf+t2]]; e0.cnt++; e0.u += val.e0;
            var e1 = ofs[ed[t2*nf+t0]]; e1.cnt++; e1.u += val.e1;
            var e2 = ofs[ed[t0*nf+t1]]; e2.cnt++; e2.u += val.e2;

            if (!skipTs)
                ots = ots.concat([
                    t0,u2,u1, t1,u0,u2, t2,u1,u0,
                    u0,u1,u2
                ]);
        }

        return {
            fs: ofs,
            ts: ots
        };
    }

    // compute circumcircle of (c1,c2,c3) and whether point is contained within it.
    public static inline function circumscribed(point:P2, c1:P2, c2:P2, c3:P2, circle:Vec3) {
        var ax = c1.x;
        var ay = c1.y;
        var bx = c2.x - ax;
        var by = c2.y - ay;
        var cx = c3.x - ax;
        var cy = c3.y - ay;
        var px = point.x - ax;
        var py = point.y - ay;

        var d = bx*cy - by*cx;
        if (d*d < epsilon) return false;
        else {
            var id = 1/(2*d);
            var bl = bx*bx + by*by;
            var cl = cx*cx + cy*cy;
            var x = (cy*bl - by*cl)*id;
            var y = (bx*cl - cx*bl)*id;
            var cl = x*x + y*y;

            circle.x = ax + x;
            circle.y = ay + y;
            circle.z = Math.sqrt(cl);

            px -= x;
            py -= y;
            return (px*px + py*py) < cl;
        }
    }
    public static inline function lex(a:P2, b:P2) {
        return (a.x < b.x) || (a.x == b.x && a.y < b.y);
    }
    public static inline function sort<T>(xs:Array<T>, n:Int, lt:T->T->Bool) {
        for (i in 1...n) {
            var v = xs[i];
            var pos = i;
            while (pos > 0 && lt(v, xs[pos-1])) {
                xs[pos] = xs[pos-1];
                pos--;
            }
            xs[pos] = v;
        }
    }

    public static function triangulate(points:Array<P2>, ?count:Null<Int>):Array<Int> {
        var n = if (count == null) points.length else count;
        if (n < 3) return [];

        var triangles = [];
        var complete = [];
        var edges = [];

        sort(points, n, lex);

        var minx = points[0].x;
        var miny = points[0].y;
        var maxx = minx;
        var maxy = miny;
        for (i in 1...n) {
            var p = points[i];
            if (p.x < minx) minx = p.x;
            if (p.y < miny) miny = p.y;
            if (p.x > maxx) maxx = p.x;
            if (p.y > maxy) maxy = p.y;
        }
        var dx = maxx - minx;
        var dy = maxy - miny;
        var dmax = if (dx > dy) dx else dy;
        var midx = (maxx + minx)*0.5;
        var midy = (maxy + miny)*0.5;

        var ntri = 1;
        triangles[0] = points.length;
        triangles[1] = points.length+1;
        triangles[2] = points.length+2;
        complete[0] = false;

        points.push({x:midx-2*dmax, y:midy-dmax  });
        points.push({x:midx,        y:midy+2*dmax});
        points.push({x:midx+2*dmax, y:midy-dmax  });

        var circle:Vec3 = [0,0,0];
        for (i in 0...n) {
            var p = points[i];
            var nedge = 0;
            var j = 0;
            while (j < ntri) {
                if (complete[j]) {
                    j++;
                    continue;
                }

                var tri = j*3;
                var inside = circumscribed(
                    p, points[triangles[tri+0]], points[triangles[tri+1]], points[triangles[tri+2]], circle
                );

                if (circle.x + circle.z < p.x) complete[j] = true;

                if (inside) {
                    edges[(nedge<<1)+0] = triangles[tri+0];
                    edges[(nedge<<1)+1] = triangles[tri+1];
                    edges[(nedge<<1)+2] = triangles[tri+1];
                    edges[(nedge<<1)+3] = triangles[tri+2];
                    edges[(nedge<<1)+4] = triangles[tri+2];
                    edges[(nedge<<1)+5] = triangles[tri+0];
                    nedge += 3;
                    var ktri = (ntri-1)*3;
                    triangles[tri+0] = triangles[ktri+0];
                    triangles[tri+1] = triangles[ktri+1];
                    triangles[tri+2] = triangles[ktri+2];
                    complete[j] = complete[ntri-1];
                    ntri--;
                }
                else j++;
            }

            for (j in 0...nedge-1) {
                var je = j<<1;
                for (k in j+1...nedge) {
                    var ke = k<<1;
                    if (edges[je] == edges[ke+1] && edges[je+1] == edges[ke]) {
                        edges[je] = edges[je+1] = -1;
                        edges[ke] = edges[ke+1] = -1;
                    }
                }
            }

            for (j in 0...nedge) {
                var je = j<<1;
                if (edges[je] == -1 || edges[je+1] == -1) continue;
                var tri = ntri*3;
                triangles[tri+0] = edges[je];
                triangles[tri+1] = edges[je+1];
                triangles[tri+2] = i;
                complete[ntri] = false;
                ntri++;
            }
        }

        var i = 0;
        while (i < ntri) {
            var tri = i*3;
            if (triangles[tri+0] >= n
             || triangles[tri+1] >= n
             || triangles[tri+2] >= n) {
                var ktri = (ntri-1)*3;
                triangles[tri+0] = triangles[ktri+0];
                triangles[tri+1] = triangles[ktri+1];
                triangles[tri+2] = triangles[ktri+2];
                ntri--;
            }
            else i++;
        }

        while (triangles.length > ntri*3) triangles.pop();

        points.pop();
        points.pop();
        points.pop();

        return triangles;
    }
}