// created by florian berger (flockaroo) - 2018
// License Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.

// derived from  "medusas hairdo" https://www.shaderoo.org/?shader=IkKuv9
// using reinder's occlusion magic from  "Cubic space division #2"

const polygons = Polygons();
const Subdiv=Math.floor(Math.random()*2.999);
// the high subdivs were a bit edgy, but now more detail thanks to reinders optimizations
const RadialSegs=2+4*(3-Subdiv);
//const TangentialSegs=11;
//const RadialSegs=4*(3-Subdiv);
const TangentialSegs=6+6*(2-Subdiv);
const RandomDive=Math.random()>0.2;
const randQuat=getRandQuat();
const objQuat=getRandQuat();

Canvas.setpenopacity(1.);

const PI2 = Math.PI*2.0;

function getRandQuat()
{
    q=[Math.random(),Math.random(),Math.random(),Math.random()];
    ql=Math.sqrt(q[0]*q[0]+q[1]*q[1]+q[2]*q[2]+q[3]*q[3]);
    q[0]/=ql; q[1]/=ql; q[2]/=ql; q[3]/=ql;
    return q;
}

// Global code will be evaluated once.
const turtle = new Turtle();
const polygonList = [];
const quads = [];

function transformVecByQuat( v, q )
{
    //return v + 2.0 * cross( q.xyz, cross( q.xyz, v ) + q.w*v );
    return add3(v, scale3(cross( q, add3(cross( q, v ) , scale3(v,q[3]) )) ,2.0));
}

function mcos(x) {
    return Math.cos(x);
}

function msin(x) {
    return Math.sin(x);
}

function cos2(x) {
    return [Math.cos(x[0]),Math.cos(x[1])];
}

function sin2(x) {
    return [Math.sin(x[0]),Math.sin(x[1])];
}

function SC(x) {
    return [Math.sin(x),Math.cos(x)];
}

function add2(a,b) {
    return [a[0]+b[0],a[1]+b[1]];
}

function add3(a,b) {
    return [a[0]+b[0],a[1]+b[1],a[2]+b[2]];
}

function sub3(a,b) {
    return [a[0]-b[0],a[1]-b[1],a[2]-b[2]];
}

function dot3(a,b) {
    return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
}

function scale2(a,b) {
    return [a[0]*b,a[1]*b];
}

function scale3(a,b) {
    return [a[0]*b,a[1]*b,a[2]*b];
}

function scale4(a,b) {
    return [a[0]*b,a[1]*b,a[2]*b,a[3]*b];
}

function mul4(a,b) {
    return [a[0]*b[0],a[1]*b[1],a[2]*b[2],a[3]*b[3]];
}

function mymix(a,b,f) {
    return a*(1.0-f)+b*f;
}

function mymix22(a,b,f) {
    return [a[0]*(1.0-f[0])+b[0]*f[0],a[1]*(1.0-f[1])+b[1]*f[1]];
}

function mix3(a,b,f) {
    return add3(scale3(a,(1.0-f)),scale3(b,f));
}

function length2(a) {
    return Math.sqrt(a[0]*a[0]+a[1]*a[1]);
}

function length3(a) {
    return Math.sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);
}

function normalize3(a) {
    return scale3(a,1.0/length3(a));
}

function cross(a,b) {
    return [
        a[1]*b[2]-b[1]*a[2],
        a[2]*b[0]-b[2]*a[0],
        a[0]*b[1]-b[0]*a[1]
    ];
}

const G=(.5+Math.sqrt(5./4.));
//const PI2=(3.141592653*2.);
const PI=3.141592653;

function fract(a) {return a-Math.floor(a);}
function floor2(a) { return [Math.floor(a[0]),Math.floor(a[1])];}
function fract2(a) { return [fract(a[0]),fract(a[1])];}

// noise funcs by Morgan McGuire
https://www.shadertoy.com/view/4dS3Wd

function hash(n) { return fract(Math.sin(n) * 1.0e4); }
function hash2(p) { return fract(1.0e4 * Math.sin(17.0 * p[0] + p[1] * 0.1) * (0.1 + Math.abs(Math.sin(p[1] * 13.0 + p[0])))); }

function noise(x) {
    var i = Math.floor(x);
    var f = fract(x);
    var u = f * f * (3.0 - 2.0 * f);
    return mymix(hash(i), hash(i + 1.0), u);
}

function noise2(x) {
    var i = floor2(x);
    var f = fract2(x);
    var r1=mymix(hash2(add2(i,[0,0])), hash2(add2(i,[1,0])), f[0]);
    var r2=mymix(hash2(add2(i,[0,1])), hash2(add2(i,[1,1])), f[0]);
    return mymix(r1, r2, f[1]);
}

function getRand01Sph(pos)
{
    pos=transformVecByQuat( pos, randQuat );
    var res = [1024,1024];
    //var texc=((pos.xy*123.+pos.z)*res+.5)/res;
    var texc=[ pos[0]*123.+pos[2]+.5/res[0],
               pos[1]*123.+pos[2]+.5/res[1] ];
    var n=1.0-noise2(scale2(texc,256.));
    if(!RandomDive) return [0,0,0,0];
    return [n,n,n,n];
}

//const vec4 p0 = vec4( 1, G, -G ,0 )/length(vec2(1,G));

const pI = scale4([ 1, G, -G ,0 ],1.0/length2([1,G]));

/*vec3 icosaPosRaw[12] = vec3[] (
    -p0.xwz,  p0.xwy, -p0.xwy,  p0.xwz,
     p0.wyx, -p0.wzx,  p0.wzx, -p0.wyx,
     p0.yxw,  p0.zxw, -p0.zxw, -p0.yxw
);*/

const icosaPosRaw = [
    [-pI[0],-pI[3],-pI[2]/*.xwz*/],
    [ pI[0], pI[3], pI[1]/*.xwy*/],
    [-pI[0],-pI[3],-pI[1]/*.xwy*/],
    [ pI[0], pI[3], pI[2]/*.xwz*/],
    
    [ pI[3], pI[1], pI[0]/*.wyx*/],
    [-pI[3],-pI[2],-pI[0]/*.wzx*/],
    [ pI[3], pI[2], pI[0]/*.wzx*/],
    [-pI[3],-pI[1],-pI[0]/*.wyx*/],
     
    [ pI[1], pI[0], pI[3]/*.yxw*/],
    [ pI[2], pI[0], pI[3]/*.zxw*/],
    [-pI[2],-pI[0],-pI[3]/*.zxw*/],
    [-pI[1],-pI[0],-pI[3]/*.yxw*/]
];

const posIdx = [
0,  6, 1,
0, 11, 6,
1,  4, 0,
1,  8, 4,
1, 10, 8,
2,  5, 3,
2,  9, 5,
2, 11, 9,
3,  7, 2,
3, 10, 7,
4,  8, 5,
4,  9, 0,
5,  8, 3,
5,  9, 4,
6, 10, 1,
6, 11, 7,
7, 10, 6,
7, 11, 2,
8, 10, 3,
9, 11, 0
];

// get icosahedron triangle
function getIcosaTri(idx)
{
    var i1 = posIdx[(idx%20)*3+0];
    var i2 = posIdx[(idx%20)*3+1];
    var i3 = posIdx[(idx%20)*3+2];

    var p1=icosaPosRaw[i1];
    var p2=icosaPosRaw[i2];
    var p3=icosaPosRaw[i3];
    return [p1,p2,p3];
}


// subdivide 1 triangle into 4 triangles and give back closest triangle
function getTriSubDiv(idx, p1, p2, p3)
{
    var p4 = normalize3(add3(p1,p2));
    var p5 = normalize3(add3(p2,p3));
    var p6 = normalize3(add3(p3,p1));

    if     (idx==0) { p1=p1; p2=p4; p3=p6; }
    else if(idx==1) { p1=p6; p2=p5; p3=p3; }
    else if(idx==2) { p1=p6; p2=p4; p3=p5; }
    else if(idx==3) { p1=p4; p2=p2; p3=p5; }
    return [p1, p2, p3];
}


var triStripIndex = [0,1,2,1,3,2];

function mixSq3(a,b,f) { return mymix3(a,b,Math.cos(f*PI)*.5+.5); }

// cubic bezier curve in 2 dimensions
// equivalent to the function above, just in a more intuitive form
function bezierCurvePos(p, t)
{
    // combination of 2 quadric beziers
    var q=[]; 
    var r=[]; 
    for(var i=0;i<3;i++) q.push(mix3(p[i],p[i+1],t));
    for(var i=0;i<2;i++) r.push(mix3(q[i],q[i+1],t));
    return mix3(r[0],r[1],t);
} 
 
function tanCurve(p1, p2, t1, t2, x)
{
    var d1=Math.abs(dot3(t1,normalize3(sub3(p2,p1))));
    var d2=Math.abs(dot3(t2,normalize3(sub3(p2,p1))));
    var p=[];
    p.push(p1);
    p.push(add3(p1,scale3(t1,length3(sub3(p2,p1))*mymix(.7,.25,d1))));
    p.push(sub3(p2,scale3(t2,length3(sub3(p2,p1))*mymix(.7,.25,d2))));
    p.push(p2);
    return bezierCurvePos(p,x);
}

function geomTangentCurve(pos1, pos2, tan1, tan2, r1, r2, 
                          rSegNum, tSegNum, vIdx)
{
    var l = length3(sub3(pos1,pos2));
    l*=.4;
    var i=Math.floor(vIdx/3/2)%tSegNum;
    //{  // converted some loops into proper vertex index values
        var fact, fact2;
        fact=Math.max(0.,(i)/(tSegNum)); // force >=0 because of sqrt below
        fact2=Math.max(0.,(i+1)/(tSegNum)); // force >=0 because of sqrt below

        var ta = mix3(tan1,tan2,fact);
        var tn = mix3(tan1,tan2,fact2);

        var p1=tanCurve(pos1,pos2,tan1,tan2,fact);
        var p2=tanCurve(pos1,pos2,tan1,tan2,fact2);
        var p1d=tanCurve(pos1,pos2,tan1,tan2,fact+.001);
        var p2d=tanCurve(pos1,pos2,tan1,tan2,fact2+.001);

        var ta = normalize3(sub3(p1d,p1));
        var tn = normalize3(sub3(p2d,p2));

        var dph=PI*2./(rSegNum);
        //vec3 b1=normalize(vec3(ta.x,-ta.y,0));
        var b1=normalize3(cross(ta,p1));
        var b2=normalize3(cross(ta,b1));
        //vec3 b3=normalize(vec3(tn.x,-tn.y,0));
        var b3=normalize3(cross(tn,p2));
        var b4=normalize3(cross(tn,b3));
        var r_1 = mymix(r1,r2,fact);
        var r_2 = mymix(r1,r2,fact2);
        var j=Math.floor(vIdx/3/2/tSegNum)%rSegNum;
        //{
            var ph  = (j)*dph;
            var ph2 = ph+dph;
            var v1 = add3(p1,scale3(add3(scale3(b1,mcos(ph )),scale3(b2,msin(ph ))),r_1));
            var v2 = add3(p1,scale3(add3(scale3(b1,mcos(ph2)),scale3(b2,msin(ph2))),r_1));
            var v3 = add3(p2,scale3(add3(scale3(b3,mcos(ph )),scale3(b4,msin(ph ))),r_2));
            var v4 = add3(p2,scale3(add3(scale3(b3,mcos(ph2)),scale3(b4,msin(ph2))),r_2));
            var v = [v1,v2,v3,v4];
            var pos = v[triStripIndex[vIdx%6]];
            var normal = normalize3(cross(sub3(v[1],v[0]),sub3(v[2],v[0])));
        //}
    //}
    return [pos,normal]
}

function calcAngle(v1, v2)
{
    return Math.acos(dot3(v1,v2)/length3(v1)/length3(v2));
}

// distance to 2 torus segments in a triangle
// each torus segment spans from the middle of one side to the middle of another side
function geomTruchet(p1, p2, p3, dz, rSegNum, tSegNum, trNum, 
                     radius, idx )
{

    if (radius<0.0) radius=.45*dz;
    var d = 10000.0;
    var rnd =getRand01Sph(add3(add3(p1,p2),p3))[0];
    var rnd2=getRand01Sph(add3(add3(p1,p2),p3))[1];
    // random rotation of torus-start-edges
    if      (rnd>.75) { var d=p1; p1=p2; p2=d; }
    else if (rnd>.50) { var d=p1; p1=p3; p3=d; }
    else if (rnd>.25) { var d=p2; p2=p3; p3=d; }
    
    // FIXME: why is this necessary - very seldom actually!?
    if(dot3(cross(sub3(p2,p1),sub3(p3,p1)),p1)>0.0) {
        var dummy;
        dummy=p2; p2=p3; p3=dummy;
    }
    
    var tubeNum=rSegNum*tSegNum*2*3;
    var i=Math.floor(idx/(tubeNum))%trNum;
    {
        var pos1, pos2, tan1, tan2;

        var R0 = .23;
        if (i==0) { 
            pos1=mix3(p1,p2,R0); pos2=mix3(p1,p3,R0);
            tan1=scale3(normalize3(add3(cross(pos1,sub3(p2,p1)),scale3(pos1,.25))),-1.); 
            tan2=scale3(normalize3(add3(cross(pos2,sub3(p3,p1)),scale3(pos2,.25))),-1.); 
        }
        if (i==1) { 
            pos1=mix3(p1,p2,1.-R0); pos2=mix3(p1,p3,1.-R0); 
            tan1=scale3(normalize3(add3(cross(pos1,sub3(p2,p1)),scale3(pos1,-.25))),-1.); 
            tan2=scale3(normalize3(add3(cross(pos2,sub3(p3,p1)),scale3(pos2,-.25))),-1.);
        }
        if (i==2) { 
            pos1=mix3(p2,p3,R0); pos2=mix3(p3,p2,R0); 
            tan1=scale3(normalize3(add3(cross(pos1,sub3(p3,p2)),scale3(pos1,.25))),-1.); 
            tan2=scale3(normalize3(add3(cross(pos2,sub3(p2,p3)),scale3(pos2,.25))),-1.); 
        }
        var pn = geomTangentCurve(pos1,pos2,tan1,tan2,radius,radius,rSegNum,tSegNum,
                                  idx%tubeNum);
    
        return pn;
    }
}

// final shape
function geom_medusa(rNum, tNum, subdiv, idx)
{
    var p1,p2,p3;

    var icosaFaceNum = 20;
    var subDivNum = 4;
    
    var trNum = 3; // tubes per truchet segemnt
    var truchetNum=rNum*tNum*2*3*trNum; // 2 triangles * 3 vertices * trNum tubes
    
    //for(int i1=0;i1<icosaFaceNum;i1++)
    var idiv=truchetNum; for(var i=0;i<subdiv;i++) idiv*=subDivNum;
    var pi = getIcosaTri(Math.floor(idx/idiv));
    var p_subDivNum_i = 1;
    for(var i=0;i<subdiv;i++)
    {
        idiv=Math.floor(idiv/subDivNum);
        var isub = Math.floor(idx/idiv)%subDivNum;
        pi = getTriSubDiv(isub,pi[0],pi[1],pi[2]);
        p_subDivNum_i*=subDivNum;
    }
    var pn = geomTruchet(pi[0],pi[1],pi[2],0.12/(1+subdiv),rNum,tNum,trNum,-1.,idx%truchetNum);
    
    return pn;
}


function medusaTri(idx)
{
    var pos = [0,0,0];
    var normal = [0,0,0];
    var pn;
    pn=geom_medusa(RadialSegs,TangentialSegs,Subdiv,idx);
    //pn=geomTangentCurve([-2,0,0], [0,0,2], [1,0,0], [0,0,-1], .1, .1, 10, 10, idx);
    pos=pn[0];
    //pn = getIcosaTri(Math.floor(idx/3));
    //pos=pn[idx%3];
    return pos;
}


function rotX(ph,v) {
    return [ v[0],v[1]*mcos(ph)+v[2]*msin(ph), v[2]*mcos(ph)-v[1]*msin(ph) ];
}

function rotY(ph,v) {
    return [ v[0]*mcos(ph)+v[2]*msin(ph), v[1], v[2]*mcos(ph)-v[0]*msin(ph) ];
}

function project(p)
{
    p[2]+=180;
    return [p[0]/p[2]*180.,p[1]/p[2]*180.,p[2]];
}

function insertQuad(p0,p1,p2,p3)
{
    var z = p0[2]+p1[2]+p2[2]+p3[2];
    var idx=0;
    for(idx=0;idx<quads.length && quads[idx+8]<z;idx+=9);
    
    // hmm, why is the one below not working... !?
    //for(var i=0;i<quads.length;i+=9) {
    //    if(quads[i+8]>z) { idx=i; break; }
    //}
    quads.splice(idx, 0, p0[0], p0[1], p1[0], p1[1], p2[0], p2[1], p3[0], p3[1], z);
}

function walk(i) {
    var num = RadialSegs*TangentialSegs*3*20*Math.pow(4,Subdiv);
    if(i==0){
        for(let j=0;j<num;j++) {
            var p0=medusaTri(j*6);
            var p1=medusaTri(j*6+1);
            var p2=medusaTri(j*6+2);
            var p3=medusaTri(j*6+4);
            p0=scale3(p0,70.0);
            p1=scale3(p1,70.0);
            p2=scale3(p2,70.0);
            p3=scale3(p3,70.0);
            p0=transformVecByQuat(p0,objQuat);
            p1=transformVecByQuat(p1,objQuat);
            p2=transformVecByQuat(p2,objQuat);
            p3=transformVecByQuat(p3,objQuat);
            p0=project(p0);
            p1=project(p1);
            p2=project(p2);
            p3=project(p3);
    
            if(cross(sub3(p1,p0),sub3(p2,p0))[2]<0.0)
            {
                insertQuad(p0,p2,p3,p1);
            }
        }
    }
    
    var p0=[quads[i*9+0],quads[i*9+1]];
    var p1=[quads[i*9+2],quads[i*9+3]];
    var p2=[quads[i*9+4],quads[i*9+5]];
    var p3=[quads[i*9+6],quads[i*9+7]];
    const p = polygons.create();
    p.addPoints([p0[0], p0[1]], [p1[0], p1[1]], [p2[0], p2[1]], [p3[0], p3[1]]);
    p.addSegments([p0[0], p0[1]], [p1[0], p1[1]], [p2[0], p2[1]], [p3[0], p3[1]]);
    polygons.draw(turtle, p);
    /*turtle.penup();
    turtle.goto(p0);
    turtle.pendown();
    turtle.goto(p1);
    turtle.goto(p2);
    turtle.goto(p0);*/
    /*turtle.goto(p3);
    turtle.goto(p2);*/
    return i <= num/1.8;
}





////////////////////////////////////////////////////////////////
// reinder's occlusion code parts from "Cubic space division #2"
// Optimizations and code clean-up by ge1doot
////////////////////////////////////////////////////////////////

function Polygons() {
	const polygonList = [];
	const linesDrawn = [];
	const Polygon = class {
		constructor() {
			this.cp = [];       // clip path: array of [x,y] pairs
			this.dp = [];       // 2d line to draw
			this.aabb = [];     // AABB bounding box
		}
		addPoints(...points) {
		    for (let i = 0; i < points.length; i++) this.cp.push(points[i]);
		    this.aabb = this.AABB();
		}
		addSegments(...points) {
		    for (let i = 0; i < points.length; i++) this.dp.push(points[i]);
		}
		addOutline(s = 0) {
			for (let i = s, l = this.cp.length; i < l; i++) {
				this.dp.push(this.cp[i], this.cp[(i + 1) % l]);
			}
		}
		createPoly(x, y, c, r, a) {
			this.cp.length = 0;
			for (let i = 0; i < c; i++) {
				this.cp.push([
					x + Math.sin(i * Math.PI * 2 / c + a) * r,
					y + Math.cos(i * Math.PI * 2 / c + a) * r
				]);
			}
			this.aabb = this.AABB();
		}
		draw(t) {
			if (this.dp.length === 0) return;
			for (let i = 0, l = this.dp.length; i < l; i+=2) {
				const d0 = this.dp[i];
				const d1 = this.dp[i + 1];
				const line_hash =
					Math.min(d0[0], d1[0]).toFixed(2) +
					"-" +
					Math.max(d0[0], d1[0]).toFixed(2) +
					"-" +
					Math.min(d0[1], d1[1]).toFixed(2) +
					"-" +
					Math.max(d0[1], d1[1]).toFixed(2);

				if (!linesDrawn[line_hash]) {
					t.penup();
					t.goto(d0);
					t.pendown();
					t.goto(d1);
					linesDrawn[line_hash] = true;
				}
			}
		}
		AABB() {
			let xmin = 2000;
			let xmax = -2000;
			let ymin = 2000;
			let ymax = -2000;
			for (let i = 0, l = this.cp.length; i < l; i++) {
				const x = this.cp[i][0];
				const y = this.cp[i][1];
				if (x < xmin) xmin = x;
				if (x > xmax) xmax = x;
				if (y < ymin) ymin = y;
				if (y > ymax) ymax = y;
			}
			// Bounding box: center x, center y, half w, half h
			return [
				(xmin + xmax) * 0.5,
				(ymin + ymax) * 0.5,
				(xmax - xmin) * 0.5,
				(ymax - ymin) * 0.5
			];
		}
		addHatching(a, d) {
			const tp = new Polygon();
			tp.cp.push(
			    [this.aabb[0] - this.aabb[2], this.aabb[1] - this.aabb[3]],
			    [this.aabb[0] + this.aabb[2], this.aabb[1] - this.aabb[3]],
			    [this.aabb[0] + this.aabb[2], this.aabb[1] + this.aabb[3]],
			    [this.aabb[0] - this.aabb[2], this.aabb[1] + this.aabb[3]]
			);
			const dx = Math.sin(a) * d, dy = Math.cos(a) * d;
			const cx = Math.sin(a) * 200, cy = Math.cos(a) * 200;
			for (let i = 0.5; i < 150 / d; i++) {
				tp.dp.push([dx * i + cy, dy * i - cx], [dx * i - cy, dy * i + cx]);
				tp.dp.push([-dx * i + cy, -dy * i - cx], [-dx * i - cy, -dy * i + cx]);
			}
			tp.boolean(this, false);
			for (let i = 0, l = tp.dp.length; i < l; i++) this.dp.push(tp.dp[i]);
		}
		inside(p) {
			// find number of i ntersection points from p to far away
			// if even your outside
			const p1 = [0.1, -1000];
			let int = 0;
			for (let i = 0, l = this.cp.length; i < l; i++) {
				if (
					this.vec2_find_segment_intersect(
						p,
						p1,
						this.cp[i],
						this.cp[(i + 1) % l]
					) !== false
				) {
					int++;
				}
			}
			return int & 1;
		}
		boolean(p, diff = true) {
			// polygon diff algorithm (narrow phase)
			const ndp = [];
			for (let i = 0, l = this.dp.length; i < l; i+=2) {
				const ls0 = this.dp[i];
				const ls1 = this.dp[i + 1];
				// find all intersections with clip path
				const int = [];
				for (let j = 0, cl = p.cp.length; j < cl; j++) {
					const pint = this.vec2_find_segment_intersect(
						ls0,
						ls1,
						p.cp[j],
						p.cp[(j + 1) % cl]
					);
					if (pint !== false) {
						int.push(pint);
					}
				}
				if (int.length === 0) {
					// 0 intersections, inside or outside?
					if (diff === !p.inside(ls0)) {
						ndp.push(ls0, ls1);
					}
				} else {
					int.push(ls0, ls1);
					// order intersection points on line ls.p1 to ls.p2
					const cmpx = ls1[0] - ls0[0];
					const cmpy = ls1[1] - ls0[1];
					for (let i = 0, len = int.length; i < len; i++) {
					    let j = i;
					    const item = int[j];
						for (
							const db = (item[0] - ls0[0]) * cmpx + (item[1] - ls0[1]) * cmpy;
							j > 0 && (int[j - 1][0] - ls0[0]) * cmpx + (int[j - 1][1] - ls0[1]) * cmpy < db;
							j--
						) int[j] = int[j - 1];
						int[j] = item;
					}
					for (let j = 0; j < int.length - 1; j++) {
						if (
							(int[j][0] - int[j + 1][0]) ** 2 + (int[j][1] - int[j + 1][1]) ** 2 >= 0.01
						) {
							if (
								diff ===
								!p.inside([
									(int[j][0] + int[j + 1][0]) / 2,
									(int[j][1] + int[j + 1][1]) / 2
								])
							) {
								ndp.push(int[j], int[j + 1]);
							}
						}
					}
				}
			}
			this.dp = ndp;
			return this.dp.length > 0;
		}
		//port of http://paulbourke.net/geometry/pointlineplane/Helpers.cs
		vec2_find_segment_intersect(l1p1, l1p2, l2p1, l2p2) {
			const d =
				(l2p2[1] - l2p1[1]) * (l1p2[0] - l1p1[0]) -
				(l2p2[0] - l2p1[0]) * (l1p2[1] - l1p1[1]);
			if (d === 0) return false;
			const n_a =
				(l2p2[0] - l2p1[0]) * (l1p1[1] - l2p1[1]) -
				(l2p2[1] - l2p1[1]) * (l1p1[0] - l2p1[0]);
			const n_b =
				(l1p2[0] - l1p1[0]) * (l1p1[1] - l2p1[1]) -
				(l1p2[1] - l1p1[1]) * (l1p1[0] - l2p1[0]);
			const ua = n_a / d;
			const ub = n_b / d;
			if (ua >= 0 && ua <= 1 && ub >= 0 && ub <= 1) {
				return [
					l1p1[0] + ua * (l1p2[0] - l1p1[0]),
					l1p1[1] + ua * (l1p2[1] - l1p1[1])
				];
			}
			return false;
		}
	};
	return {
		list() {
			return polygonList;
		},
		create() {
			return new Polygon();
		},
		draw(turtle, p) {
			let vis = true;
			for (let j = 0; j < polygonList.length; j++) {
				const p1 = polygonList[j];
				// AABB overlapping test - still O(N2) but very fast
				if (
					Math.abs(p1.aabb[0] - p.aabb[0]) - (p.aabb[2] + p1.aabb[2]) < 0 &&
					Math.abs(p1.aabb[1] - p.aabb[1]) - (p.aabb[3] + p1.aabb[3]) < 0
				) {
					if (p.boolean(p1) === false) {
						vis = false;
						break;
					}
				}
			}
			if (vis) {
				p.draw(turtle);
				polygonList.push(p);
			}
		}
	};
}