BSP Tree

A Binary Space Partitioning Tree (BSP Tree) baseline implementation adapted to turtletoy

- resolves all z-index conflicts by splitting polygons
- tree can be traversed in perfect front-to-back order, then using Reinder's polygon's occlusion/clipping method

Created by ge1doot on 2019/2/12
153
4

Log in to post a comment.

// You can find the Turtle API reference here: https://turtletoy.net/syntax
Canvas.setpenopacity(1);
const turtle = new Turtle();
turtle.penup();
/* ===================================================
*  ------------ BSP Tree  -------------
* adapted from C++ BSP Tree Demo
* Copyright (c)1994-97 Bretton Wade
* http://www.gamers.org/dhs/helpdocs/bsp_faq.html
* =================================================== */

const polygons = Polygons();

///////////////////////// vec 4 /////////////////////////
const v4 = {
	multiplyScalar (v, s) {
		return [
			v[0] * s,
			v[1] * s,
			v[2] * s,
			v[3] * s
		];
	},
	subtract (v, p) {
		return [
			v[0] - p[0],
			v[1] - p[1],
			v[2] - p[2],
			v[3] - p[3]
		];
	},
	add (v, p) {
		return [
			v[0] + p[0],
			v[1] + p[1],
			v[2] + p[2],
			v[3] + p[3]
		];
	},
	multiplyMatrix (v, m) {
		const px = v[0], py = v[1], pz = v[2], pw = v[3];
		return [
			px * m[0] + py * m[4] + pz * m[8] +  pw * m[12],
			px * m[1] + py * m[5] + pz * m[9] +  pw * m[13],
			px * m[2] + py * m[6] + pz * m[10] + pw * m[14],
			px * m[3] + py * m[7] + pz * m[11] + pw * m[15]
		];
	},
	crossProduct (v, p) {
		return [
			v[1] * p[2] - v[2] * p[1], 
			v[2] * p[0] - v[0] * p[2],
			v[0] * p[1] - v[1] * p[0],
			0
		];
	},
	dotProduct (v, p) {
		return (
			v[0] * p[0] +
			v[1] * p[1] +
			v[2] * p[2] +
			v[3] * p[3]
		);
	},
	normalize (v) {
		const n = Math.sqrt(
			v[0] * v[0] +
			v[1] * v[1] +
			v[2] * v[2] +
			v[3] * v[3]
		);
		return [
			v[0] / n,
			v[1] / n,
			v[2] / n,
			v[3] / n
		];
	}
};
///////////////////////// mat 4 ///////////////////////
const m4 = {
	identity() {
		return [
			1, 0, 0, 0,
			0, 1, 0, 0,
			0, 0, 1, 0,
			0, 0, 0, 1
		];
	},
	multiply (m, n) {
		const
			m1XX = m[0], m1XY = m[1], m1XZ = m[2], m1XW = m[3],
			m1YX = m[4], m1YY = m[5], m1YZ = m[6], m1YW = m[7],
			m1ZX = m[8], m1ZY = m[9], m1ZZ = m[10], m1ZW = m[11],
			m1WX = m[12], m1WY = m[13], m1WZ = m[14], m1WW = m[15],
			m2XX = n[0], m2XY = n[1], m2XZ = n[2], m2XW = n[3],
			m2YX = n[4], m2YY = n[5], m2YZ = n[6], m2YW = n[7],
			m2ZX = n[8], m2ZY = n[9], m2ZZ = n[10], m2ZW = n[11],
			m2WX = n[12], m2WY = n[13], m2WZ = n[14], m2WW = n[15];
		return [
			m1XX * m2XX + m1XY * m2YX + m1XZ * m2ZX + m1XW * m2WX,
			m1XX * m2XY + m1XY * m2YY + m1XZ * m2ZY + m1XW * m2WY,
			m1XX * m2XZ + m1XY * m2YZ + m1XZ * m2ZZ + m1XW * m2WZ,
			m1XX * m2XW + m1XY * m2YW + m1XZ * m2ZW + m1XW * m2WW,
			m1YX * m2XX + m1YY * m2YX + m1YZ * m2ZX + m1YW * m2WX,
			m1YX * m2XY + m1YY * m2YY + m1YZ * m2ZY + m1YW * m2WY,
			m1YX * m2XZ + m1YY * m2YZ + m1YZ * m2ZZ + m1YW * m2WZ,
			m1YX * m2XW + m1YY * m2YW + m1YZ * m2ZW + m1YW * m2WW,
			m1ZX * m2XX + m1ZY * m2YX + m1ZZ * m2ZX + m1ZW * m2WX,
			m1ZX * m2XY + m1ZY * m2YY + m1ZZ * m2ZY + m1ZW * m2WY,
			m1ZX * m2XZ + m1ZY * m2YZ + m1ZZ * m2ZZ + m1ZW * m2WZ,
			m1ZX * m2XW + m1ZY * m2YW + m1ZZ * m2ZW + m1ZW * m2WW,
			m1WX * m2XX + m1WY * m2YX + m1WZ * m2ZX + m1WW * m2WX,
			m1WX * m2XY + m1WY * m2YY + m1WZ * m2ZY + m1WW * m2WY,
			m1WX * m2XZ + m1WY * m2YZ + m1WZ * m2ZZ + m1WW * m2WZ,
			m1WX * m2XW + m1WY * m2YW + m1WZ * m2ZW + m1WW * m2WW
		];
	},
	inverse (m) {
		const
			a00 = m[0],  a01 = m[1],  a02 = m[2],  a03 = m[3],
			a10 = m[4],  a11 = m[5],  a12 = m[6],  a13 = m[7],
			a20 = m[8],  a21 = m[9],  a22 = m[10], a23 = m[11],
			a30 = m[12], a31 = m[13], a32 = m[14], a33 = m[15],
			b00 = a00 * a11 - a01 * a10,
			b01 = a00 * a12 - a02 * a10,
			b02 = a00 * a13 - a03 * a10,
			b03 = a01 * a12 - a02 * a11,
			b04 = a01 * a13 - a03 * a11,
			b05 = a02 * a13 - a03 * a12,
			b06 = a20 * a31 - a21 * a30,
			b07 = a20 * a32 - a22 * a30,
			b08 = a20 * a33 - a23 * a30,
			b09 = a21 * a32 - a22 * a31,
			b10 = a21 * a33 - a23 * a31,
			b11 = a22 * a33 - a23 * a32;
		let det = b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06;
		if (!det) return this; 
		det = 1.0 / det;
		return [
			(a11 * b11 - a12 * b10 + a13 * b09) * det,
			(a02 * b10 - a01 * b11 - a03 * b09) * det,
			(a31 * b05 - a32 * b04 + a33 * b03) * det,
			(a22 * b04 - a21 * b05 - a23 * b03) * det,
			(a12 * b08 - a10 * b11 - a13 * b07) * det,
			(a00 * b11 - a02 * b08 + a03 * b07) * det,
			(a32 * b02 - a30 * b05 - a33 * b01) * det,
			(a20 * b05 - a22 * b02 + a23 * b01) * det,
			(a10 * b10 - a11 * b08 + a13 * b06) * det,
			(a01 * b08 - a00 * b10 - a03 * b06) * det,
			(a30 * b04 - a31 * b02 + a33 * b00) * det,
			(a21 * b02 - a20 * b04 - a23 * b00) * det,
			(a11 * b07 - a10 * b09 - a12 * b06) * det,
			(a00 * b09 - a01 * b07 + a02 * b06) * det,
			(a31 * b01 - a30 * b03 - a32 * b00) * det,
			(a20 * b03 - a21 * b01 + a22 * b00) * det
		];
	},
	rotateXYZ (angleX, angleY, angleZ) {
		const cw = Math.cos(angleX),
			sw = Math.sin(angleX),
			cy = Math.cos(angleY),
			sy = Math.sin(angleY),
			ck = angleZ ? Math.cos(angleZ) : 1,
			sk = angleZ ? Math.sin(angleZ) : 0;
		return [
			cy*ck, cw*sk+sw*sy*ck, sw*sk-cw*sy*ck, 0,
			-cy*sk, cw*ck-sw*sy*sk, sw*ck+cw*sy*sk, 0,
			sy, -sw*cy, cw*cy, 0,
			0, 0, 0, 1
		];
	},
	scale (m, x, y, z) {
		return this.multiply(m, 
			[
				x, 0, 0, 0,
				0, y, 0, 0,
				0, 0, z, 0,
				0, 0, 0, 1
			]
		);
	},
	translate (m, x, y, z) {
		return this.multiply(m,
			[
				1, 0, 0, 0,
				0, 1, 0, 0,
				0, 0, 1, 0,
				x, y, z, 1
			]
		);
	},
	rotateX (m, angle) {
		angle = angle * (Math.PI / 180);
		const c = Math.cos(angle),
			s = Math.sin(angle);
		return this.multiply(m, 
			[
				1, 0, 0, 0,
				0, c, s, 0,
				0,-s, c, 0,
				0, 0, 0, 1
			]
		);
	},
	rotateY (m, angle) {
		angle = angle * (Math.PI / 180);
		const c = Math.cos(angle),
			s = Math.sin(angle);
		return this.multiply(m, 
			[
				c, 0,-s, 0,
				0, 1, 0, 0,
				s, 0, c, 0,
				0, 0, 0, 1
			]
		);
	},
	rotateZ (m, angle) {
		angle = angle * (Math.PI / 180);
		const c = Math.cos(angle),
			s = Math.sin(angle);
		return this.multiply(m, 
			[
				c,-s, 0, 0,
				s, c, 0, 0,
				0, 0, 1, 0,
				0, 0, 0, 1
			]
		);
	}
};
////////////////////////// chaining ////////////////////////////
const mat4 = () => new Mat4();
const Mat4 = class {
	constructor () {
		this.m = m4.identity();
	}
	scale (x, y, z) {
		this.m = m4.scale(this.m, x, y, z);
		return this;
	}
	translate(x, y, z) {
		this.m = m4.translate(this.m, x, y, z);
		return this;
	}
	rotateX(a) {
		this.m = m4.rotateX(this.m, a);
		return this;
	}
	rotateY(a) {
		this.m = m4.rotateY(this.m, a);
		return this;
	}
	rotateZ(a) {
		this.m = m4.rotateZ(this.m, a);
		return this;
	}
};
/////////////////////// Geometry ////////////////////////
const cube = (transform) => {
	const p = [
		v4.multiplyMatrix([ 1,  1,  1, 1], transform),
		v4.multiplyMatrix([-1,  1,  1, 1], transform),
		v4.multiplyMatrix([-1, -1,  1, 1], transform),
		v4.multiplyMatrix([ 1, -1,  1, 1], transform),
		v4.multiplyMatrix([ 1,  1, -1, 1], transform),
		v4.multiplyMatrix([-1,  1, -1, 1], transform),
		v4.multiplyMatrix([-1, -1, -1, 1], transform),
		v4.multiplyMatrix([ 1, -1, -1, 1], transform)
	];
	return [
		polygon.createIndexed(p, [0, 1, 2, 3]),
		polygon.createIndexed(p, [7, 6, 5, 4]),
		polygon.createIndexed(p, [0, 3, 7, 4]),
		polygon.createIndexed(p, [0, 4, 5, 1]),
		polygon.createIndexed(p, [5, 6, 2, 1]),
		polygon.createIndexed(p, [3, 2, 6, 7])
	];
};
const pyramid = (transform, nosplit) => {
	const p = [
		v4.multiplyMatrix([ 1,  0,  1, 1], transform),
		v4.multiplyMatrix([-1,  0,  1, 1], transform),
		v4.multiplyMatrix([-1,  0, -1, 1], transform),
		v4.multiplyMatrix([ 1,  0, -1, 1], transform),
		v4.multiplyMatrix([ 0,  2,  0, 1], transform)
	];
	return [
		polygon.createIndexed(p, [0, 1, 2, 3]),
		polygon.createIndexed(p, [4, 1, 0]),
		polygon.createIndexed(p, [4, 0, 3]),
		polygon.createIndexed(p, [4, 2, 1]),
		polygon.createIndexed(p, [4, 3, 2])
	];
};
const diamond = (transform, nosplit) => {
	const p = [
		v4.multiplyMatrix([ 1,  0,  1, 1], transform),
		v4.multiplyMatrix([-1,  0,  1, 1], transform),
		v4.multiplyMatrix([-1,  0, -1, 1], transform),
		v4.multiplyMatrix([ 1,  0, -1, 1], transform),
		v4.multiplyMatrix([ 0,  1,  0, 1], transform),
		v4.multiplyMatrix([ 0, -1,  0, 1], transform)
	];
	return [
		polygon.createIndexed(p, [4, 1, 0]),
		polygon.createIndexed(p, [4, 0, 3]),
		polygon.createIndexed(p, [4, 2, 1]),
		polygon.createIndexed(p, [4, 3, 2]),
		polygon.createIndexed(p, [0, 1, 5]),
		polygon.createIndexed(p, [3, 0, 5]),
		polygon.createIndexed(p, [1, 2, 5]),
		polygon.createIndexed(p, [2, 3, 5])
	];
};
/////////////////////////// Polygon ///////////////////////////
const polygon = {
	createIndexed (p, list) {
		const points = [];
		for (const index of list) {
			points.push(p[index]);
		}
		points.plane = this.definePlane(points);
		return points;
	},
	create (p) {
		p.plane = this.definePlane(p);
		return p;
	},
	definePlane (poly) {
		let sum = [0,0,0,0];
		for (
			var i = 0, last = poly.length - 1; 
			i < poly.length; 
			last = i, i++
		) {
			const A = poly[last]; // assumes counter-clockwise points order
			const B = poly[i];
			sum[0] += (A[1] - B[1]) * (A[2] + B[2]);
			sum[1] += (A[2] - B[2]) * (A[0] + B[0]);
			sum[2] += (A[0] - B[0]) * (A[1] + B[1]);
		}
		sum = v4.normalize(sum);
		const p = poly[0];
		return [
			sum[0], 
			sum[1], 
			sum[2], 
			- (
				sum[0] * p[0] +
				sum[1] * p[1] +
				sum[2] * p[2] +
				sum[3] * p[3]
			)
		];
	},
	draw (camera, poly) {
		const p = polygons.create();
		for (const p0 of poly) {
			if (p0.frame !== camera.frame) {
				const p2 = v4.multiplyMatrix(p0, camera.viewing);
				p0.frame = camera.frame;
				p0.x = (p2[0] / p2[3]) *  camera.aspect;
				p0.y = (p2[1] / p2[3]) * -camera.aspect;
			}
			p.addPoints([p0.x, p0.y]);
		}
		p.addOutline(0);
		polygons.draw(turtle, p);
	},
	split (poly, plane) {
		const outpts = [], inpts = [];
		let outp, inp, out_c = 0, in_c = 0, poly_class = HC_ON;
		let ptA = poly[poly.length - 1], sideA = v4.dotProduct(ptA, plane);
		for (const ptB of poly) {
			const sideB = v4.dotProduct(ptB, plane);
			if (sideB > EPSILON) {
				if (poly_class === HC_ON) poly_class = HC_OUT;
				else if (poly_class !== HC_OUT) poly_class = HC_SPANNING;
				if (sideA < -EPSILON) {
					const v = v4.subtract(ptB, ptA);
					outpts[out_c++] = inpts[in_c++] = v4.add(ptA, v4.multiplyScalar(v, -v4.dotProduct(ptA, plane) / v4.dotProduct(v, plane)));
					poly_class = HC_SPANNING;
				}
				outpts[out_c++] = ptB;
			} else if (sideB < -EPSILON) {
				if (poly_class === HC_ON) poly_class = HC_IN;
				else if (poly_class !== HC_IN) poly_class = HC_SPANNING;
				if (sideA > EPSILON) {
					const v = v4.subtract(ptB, ptA);
					outpts[out_c++] = inpts[in_c++] = v4.add(ptA, v4.multiplyScalar(v, -v4.dotProduct(ptA, plane) / v4.dotProduct(v, plane)));
					poly_class = HC_SPANNING;
				}
				inpts[in_c++] = ptB;
			} else outpts[out_c++] = inpts[in_c++] = ptB;
			ptA = ptB;
			sideA = sideB;
		}
		switch (poly_class) {
			case HC_OUT: 
				outp = poly;
				break;
			case HC_IN:
				inp = poly;
				break;
			case HC_SPANNING:
				outp = this.create(outpts);
				inp = this.create(inpts);
				break;
		}
		return [outp, inp, poly_class];
	}
}
////////////////////////// BSP Tree class //////////////////////////////
const BspTree = class {
	constructor () {
		this.node = null;
	}
	insert (list, keep, cur = keep) {
		if (list.length === 0) return;
		if (this.node) {
			this.node.insert(list, keep);
		} else {
			if ((cur === keep) || (keep === HC_SPANNING)) {
				this.node = new BspNode(list.pop());
				if (list.length) this.node.insert(list, HC_SPANNING);
			}
		}
	}
	pushPoly (poly, result, keep, cur) {
		if (this.node !== null) this.node.pushPoly (poly, result, keep);
		else if (cur === keep) result.push(poly);
	}
	pushList (list, result, keep, cur) {
		if (list.length === 0) return;
		if (this.node !== null) result = this.node.pushList (list, result, keep);
		else if (cur === keep) result.push.apply(result, list);
	}
	reduce  () {
		if (this.node !== null) this.node.reduce();
	}
	draw (camera) {
		if (this.node !== null) this.node.draw(camera);
	}
}
////////////////////////// Bsp Node Class ////////////////////////////////	
const BspNode = class {
	constructor (poly) {
		this.plane = poly.plane;
		this.on    = [poly];
		this.in    = new BspTree();
		this.out   = new BspTree();
	}
	insert (list, keep) {
		const inside = [], outside = [];
		for (const poly of list) {
			const [outp, inp, sgn] = polygon.split(poly, this.plane);
			if (sgn === HC_ON) this.on.push(poly);
			else {
				if ((sgn === HC_IN) || (sgn === HC_SPANNING)) inside.push(inp);
				if ((sgn === HC_OUT) || (sgn === HC_SPANNING)) outside.push(outp);
			}
		}
		if (inside.length !== 0) this.in.insert(inside, keep, HC_IN);
		if (outside.length !== 0) this.out.insert(outside, keep, HC_OUT);
	}
	pushPoly (poly, result, keep) {
		const [outp, inp, sgn] = polygon.split(poly, this.plane);
		if (sgn === HC_ON) result.push(poly);
		else {
			if ((sgn === HC_IN) || (sgn === HC_SPANNING)) this.in.pushPoly (inp, result, keep, HC_IN);
			if ((sgn === HC_OUT) || (sgn === HC_SPANNING)) this.out.pushPoly (outp, result, keep, HC_OUT);
		}
	}
	pushList (list, result, keep) {
		const inside = [], outside = [];
		for (const poly of list) {
			const [outp, inp, sgn] = polygon.split(poly, this.plane);
			if (sgn === HC_ON) result.push(poly);
			else {
				if ((sgn === HC_IN) || (sgn === HC_SPANNING)) inside.push(inp);
				if ((sgn === HC_OUT) || (sgn === HC_SPANNING)) outside.push(outp);
			}
		}
		if (inside.length !== 0) this.in.pushList (inside, result, keep, HC_IN);
		if (outside.length !== 0) this.out.pushList (outside, result, keep, HC_OUT);
	}
	reduce () {	
		const results = [], boundary = [];
		for (const poly of this.on) {
			if (Math.abs(poly.plane[3] + this.plane[3]) > EPSILON) {
				this.in.pushPoly (poly, results, HC_IN, HC_IN);
				this.out.pushList (results, boundary, HC_OUT, HC_OUT);
			} else {
				this.out.pushPoly (poly, results, HC_OUT, HC_OUT);
				this.in.pushList (results, boundary, HC_IN, HC_IN);
			}
		}
		this.on = boundary;
		this.in.reduce();
		this.out.reduce();
	}
	// front to back tree traversal
	draw  (camera) { 
		var sgn = v4.dotProduct(camera.eye, this.plane);
		if (sgn >= 0) {
			this.out.draw(camera);
			for (const on of this.on) polygon.draw(camera, on);
			this.in.draw(camera);
		} else {
			this.in.draw(camera);
			this.out.draw(camera);
		}
	}
}
////////////////// camera /////////////////////
const camera = {
	frame: 1,
	aspect: 100,
	drawScene (world, rx, ry, rz) {
		this.frame++;
		const transformation = m4.rotateXYZ(rx, ry, rz);
		this.viewing = m4.multiply(transformation, this.view);
		this.eye = v4.multiplyMatrix(this.posEye, m4.inverse(transformation)); 
		world.draw(camera);
	},
	look (e, to, zoom) {
		this.posEye = [e[0],e[1],e[2],1];
		const vpn = v4.normalize(v4.subtract(this.posEye, [to[0],to[1],to[2],1]));
		const u = v4.crossProduct([0, 1, 0], vpn);
		const v = v4.crossProduct(vpn, u);
		const vrp = v4.add(this.posEye, v4.multiplyScalar(vpn, -zoom));
		this.view = m4.multiply(this.viewMatrix(u, v, vpn, vrp), this.perspective(zoom));
	},
	perspective (distance) {
		return [
			1, 0, 0, 0,
			0, 1, 0, 0,
			0, 0, 1,-1 / distance,
			0, 0, 0, 1
		];
	},
	viewMatrix (u, v, n, r) {
		return [
			u[0], v[0], n[0], 0,
			u[1], v[1], n[1], 0,
			u[2], v[2], n[2], 0,
			-(r[0] * u[0] + r[1] * u[1] + r[2] * u[2] + r[3] * u[3]),
			-(r[0] * v[0] + r[1] * v[1] + r[2] * v[2] + r[3] * v[3]),
			-(r[0] * n[0] + r[1] * n[1] + r[2] * n[2] + r[3] * n[3]), 1
		];
	}
}
		
////////////////////////////////////////////////////////////////////////////////

const HC_IN = -1, HC_ON = 0, HC_OUT = 1, HC_SPANNING = 2, INFINITY = 100000, EPSILON  = 1 / 100000;

const world = new BspTree();
world.insert(cube(mat4().scale(-5,-2.5,-5).m), 1);// outside walls
world.insert(cube(mat4().m), 1); // insert a basic cube
world.insert(cube(mat4().scale(-1.5,-0.875,-1.5).translate(0.625,0,0).m), -1); // cut out the cube to make a big "C"
world.insert(cube(mat4().scale(0.15,0.4,0.4).translate(-0.95,0,0).rotateX(45).m), 1); // window frame
world.insert(cube(mat4().scale(-0.3,-0.3,-0.3).translate(-1,0,0).rotateX(45).m), -1); // window hole
world.insert(cube(mat4().scale(0.15,0.4,0.4).translate(0.8,0,0).rotateX(45).m), 1); // window frame
world.insert(cube(mat4().scale(0.0625,1.8,0.0625).translate(0.8,0,0).m), 1); // support
world.insert(cube(mat4().scale(-1.5,-0.3,-0.3).translate(0.8,0,0).rotateX(45).m), -1); // window hole
world.insert(diamond(mat4().rotateZ(90).scale(4,0.15,0.15).rotateX(45).m), 1); // long diamond
world.insert(pyramid(mat4().scale(0.5,0.5,0.5).translate(0,-1.25,0).rotateZ(180).m), 1); // pyramid
world.insert(pyramid(mat4().scale(0.5,0.5,0.5).translate(0,-1.25,0).m), 1); // pyramid
world.insert(cube(mat4().scale(0.5,0.05,0.5).translate(0,-1.4,0).m), 1);// slab
world.insert(cube(mat4().scale(0.51,0.05,0.51).translate(0,-1.6,0).m), 1); // slab
world.insert(cube(mat4().scale(0.5,0.05,0.5).translate(0,1.4,0).m), 1); // slab
world.insert(cube(mat4().scale(0.51,0.05,0.51).translate(0,1.6,0).m), 1); // slab
world.reduce(); // extrude negative parts
camera.look([0,0,8], [0,0,0], 3);
camera.drawScene(world, 0.2, Math.random() * 2 * Math.PI, 0);


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

function Polygons() {
	const polygonList = [];
	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();
		}
		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]);
			}
		}
		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];
				t.penup();
				t.goto(d0);
				t.pendown();
				t.goto(d1);
			}
		}
		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
			const p1 = [0.1, -1000];
			let int = 0;
			for (let i = 0, l = this.cp.length; i < l; i++) {
				if ( (p[0]-this.cp[i][0])**2 +  (p[1]-this.cp[i][1])**2 <= 0.001) return false;
				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
			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 {
		create() {
			return new Polygon();
		},
		draw(t, 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(t);
				polygonList.push(p);
			}
		}
	};
}