### Cubic space division #3

This turtle is a combination of "Cubic space division" by Escher and the Droste spiral of "Print Gallery" (also by Escher).

Note that I divide by e^(0.05*z) vs the normal divide by z to project the vertices to the screen. This way I make sure the projected grid is a Droste image.

#Droste #polygons #3D #Escher

```// Cubic space division #3. Created by Reinder Nijhoff 2021
// @reindernijhoff
//
// https://turtletoy.net/turtle/e7a276c605
//

const Perspective = 0; // min=0, max=2, step=1 (Droste perspective and Escher spiral, Droste perspective, Real perspective)

let viewProjectionMatrix;

const lw = .25;
const ll =  5;
const lz =  5;
const dim = 9;
const dimz = 21;

let polys, turtle;

const ls = Math.log(15.642631884188175);
const sr = 4 * (Math.PI**2) / ((ls**2) + 4 * (Math.PI**2));
const si = 2 * Math.PI * ls / ((ls**2) + 4 * (Math.PI**2));

let zoom = 0;

function EscherDroste(p) {
// http://www.ams.org/notices/200304/fea-escher.pdf
// h(w) = w^((2πi + log scale)/(2πi))
const pr = Math.log(Math.hypot(p[0]/100, p[1]/100));

let pi = Math.atan2(p[0], p[1]) + .5; //+(0.4/256.)*deformationScale;

let l = Math.exp(pr * sr - pi * si);
const a = pr * si + pi * sr;

// barrel distortion
l =  Math.pow(l, 0.88) * 100;

return [l * Math.cos( a ), l * Math.sin( a )];
}

function walk(i, t) {
zoom = t;
if (i == 0) {
turtle = new Tortoise();
if (Perspective == 0) {
}
setupCamera();
polys = new Polygons();
}
const d = dim;
let z = 23 -((i/d)|0) * (2*lz+1);
let x, y = 0;

x = -((i%d)) * (2*ll+1);
drawCube(turtle, x,y+ll+lw*2,z, lw,ll-lw*2,lw ,2, true);
x = ((i%d)+1) * (2*ll+1);
drawCube(turtle, x,y+ll+lw*2,z, lw,ll-lw*2,lw ,0, true);

x = 0;
y = -((i%d)) * (2*ll+1);
drawCube(turtle, x+ll+lw*2,y,z, ll-lw*2,lw,lw ,1, true);
y = ((i%d)+1) * (2*ll+1);
drawCube(turtle, x+ll+lw*2,y,z, ll-lw*2,lw,lw ,3, true);

for (let j=0; j<d; j++) {
z = 23 -(((i/d)|0) - 1) * (2*lz+1);

x = -((i%d)) * (2*ll+1);
y = ((j+1)) * (2*ll+1);

drawCube(turtle, x,y,z,    1,1,1    ,3);
drawCube(turtle, x-ll,y,z, ll,lw,lw ,3, true);
drawCube(turtle, x,y+ll,z, lw,ll,lw ,2, true);
drawCube(turtle, x,y,z-lz, lw,lw,lz ,3, false);

x = -(i%d) * (2*ll+1);
y = -(j) * (2*ll+1);

drawCube(turtle,   x,y,z,    1,1,1  ,2);
drawCube(turtle, x-ll,y,z, ll,lw,lw ,1, true);
drawCube(turtle, x,y-ll,z, lw,ll,lw ,2, true);
drawCube(turtle, x,y,z-lz, lw,lw,lz ,2, false);

x = ((i%d) + 1) * (2*ll+1);
y = -(j) * (2*ll+1);

drawCube(turtle,   x,y,z,    1,1,1  ,1);
drawCube(turtle, x+ll,y,z, ll,lw,lw ,1, true);
drawCube(turtle, x,y-ll,z, lw,ll,lw ,0, true);
drawCube(turtle, x,y,z-lz, lw,lw,lz ,1, false);

x = ((i%d) + 1) * (2*ll+1);
y = ((j) + 1) * (2*ll+1);

drawCube(turtle,   x,y,z,    1,1,1  ,0);
drawCube(turtle, x+ll,y,z, ll,lw,lw ,3, true);
drawCube(turtle, x,y+ll,z, lw,ll,lw ,0, true);
drawCube(turtle, x,y,z-lz, lw,lw,lz ,0, false);
}

return i < dim*dimz;
}

function drawCube(turtle, x, y, z, w, h, d, mode = 0, hideCap = false) {
if (40 - (dimz + mode - 3) * (2*lz+1) > z - h) return;

let vl = [
[x-w,y-h,z+d,1],
[x+w,y-h,z+d,1],
[x+w,y+h,z+d,1],
[x-w,y+h,z+d,1],
[x-w,y+h,z-d,1],
[x+w,y+h,z-d,1],
[x+w,y-h,z-d,1],
[x-w,y-h,z-d,1],
];

let il;

switch(mode) {
case 0: il = [0,1,2,3, 3,4,7,0, 0,7,6,1]; break;
case 1: il = [0,1,2,3, 5,2,3,4, 3,4,7,0]; break;
case 2: il = [0,1,2,3, 1,2,5,6, 5,2,3,4]; break;
case 3: il = [0,1,2,3, 0,7,6,1, 1,2,5,6]; break;
}
const e = hideCap ? 2 : 3;
for (let i=0; i<e; i++) {
let vis = true;
const cp = [];
for (let j=0; j<4; j++) {
const v = transform4(vl[il[j+i*4]], viewProjectionMatrix);
const d = Perspective == 2 ? 100/v[3] : Math.exp(-0.05*v[3]) * 25;
cp.push([v[0]*d , -v[1]*d]);
if (v[3] < (Perspective == 2 ? 3 : -5)) {
vis = false;
break;
}
}
if (vis) {
vis = true;
for (let j=0; j<cp.length; j++) {
if (Math.abs(cp[j][0]) > 800 || Math.abs(cp[j][1]) > 800) {
vis = false;
break;
}
}
if (vis) {
const p =  polys.create();
polys.draw(turtle, p);
}
}
}
}

function setupCamera() {
viewMatrix = lookAt4m([ll+.5,ll+.5,20 - zoom*11], [ll+.5,ll+.5,0], [0.1,1,0]);
// viewMatrix = lookAt4m([ll,ll,31 - zoom*11], [ll,ll,0], [0,1,0]);

projectionMatrix = perspective4m(.5, 1, 0);
viewProjectionMatrix = multiply4m(projectionMatrix, viewMatrix);
}

////////////////////////////////////////////////////////////////
// Tortoise utility code. Created by Reinder Nijhoff 2019
// https://turtletoy.net/turtle/102cbd7c4d
////////////////////////////////////////////////////////////////

function Tortoise(x, y) {
class Tortoise extends Turtle {
constructor(x, y) {
super(x, y);
this.ps = Array.isArray(x) ? [...x] : [x || 0, y || 0];
this.transforms = [];
}
this.transforms.push(t);
this.jump(this.ps);
return this;
}
applyTransforms(p) {
if (!this.transforms) return p;
let pt = [...p];
this.transforms.map(t => { pt = t(pt); });
return pt;
}
goto(x, y, d = 0) {
const p = Array.isArray(x) ? [...x] : [x, y];
const pt = this.applyTransforms(p);
if (d > 10) {
super.jump(p);
}
if (this.isdown() && (this.pt[0]-pt[0])**2 + (this.pt[1]-pt[1])**2 > 12) {
this.goto((this.ps[0]+p[0])/2, (this.ps[1]+p[1])/2, d+1);
this.goto(p[0], p[1], d+1);
} else {
super.goto(pt);
this.ps = p;
this.pt = pt;
}
}
position() { return this.ps; }
}
return new Tortoise(x,y);
}

////////////////////////////////////////////////////////////////
// Polygon Clipping utility code - Created by Reinder Nijhoff 2019
// https://turtletoy.net/turtle/a5befa1f8d
////////////////////////////////////////////////////////////////

function Polygons() {
const Node = class {
constructor(aabb, depth, maxDepth) {
this.aabb = aabb;
this.depth = depth;
this.polys = [];
this.children = [];

if (depth < maxDepth) {
const half = aabb[2]/2;
this.children.push(new Node([aabb[0]-half, aabb[1]-half, half, half], depth+1, maxDepth));
this.children.push(new Node([aabb[0]+half, aabb[1]-half, half, half], depth+1, maxDepth));
this.children.push(new Node([aabb[0]+half, aabb[1]+half, half, half], depth+1, maxDepth));
this.children.push(new Node([aabb[0]-half, aabb[1]+half, half, half], depth+1, maxDepth));
}
}
inside(aabb) {
return (Math.abs(this.aabb[0] - aabb[0]) < this.aabb[2] - aabb[2] &&
Math.abs(this.aabb[1] - aabb[1]) < this.aabb[3] - aabb[3]);
}
cover(aabb) {
return !(Math.abs(this.aabb[0] - aabb[0]) - (aabb[2] + this.aabb[2]) >= 0 &&
Math.abs(this.aabb[1] - aabb[1]) - (aabb[3] + this.aabb[3]) >= 0);
}
for (let i=0; i<this.children.length; i++) {
if (this.children[i].inside(p.aabb)) {
return;
}
}
this.polys.push(p);
}
boolean(p) {
for (let j = 0; j < this.polys.length; j++) {
if (!p.boolean(this.polys[j])) {
return false;
}
}
for (let i=0; i<this.children.length; i++) {
if (this.children[i].cover(p.aabb)) {
if (!this.children[i].boolean(p)) {
return false;
}
}
}
return true;
}
}
const Polygon = class {
constructor() {
this.cp = [];       // clip path: array of [x,y] pairs
this.dp = [];       // 2d lines [x0,y0],[x1,y1] to draw
this.aabb = [];     // AABB bounding box
}
// add point to clip path and update bounding box
let xmin = 1e5, xmax = -1e5, ymin = 1e5, ymax = -1e5;
(this.cp = [...this.cp, ...points]).forEach( p => {
xmin = Math.min(xmin, p[0]), xmax = Math.max(xmax, p[0]);
ymin = Math.min(ymin, p[1]), ymax = Math.max(ymax, p[1]);
});
this.aabb = [(xmin+xmax)/2, (ymin+ymax)/2, (xmax-xmin)/2, (ymax-ymin)/2];
}
// add segments (each a pair of points)
points.forEach(p => this.dp.push(p));
}
for (let i = 0, l = this.cp.length; i < l; i++) {
this.dp.push(this.cp[i], this.cp[(i + 1) % l]);
}
}
draw(t) {
for (let i = 0, l = this.dp.length; i < l; i+=2) {
t.jump(this.dp[i]), t.goto(this.dp[i + 1]);
}
}
const tp = new Polygon();
tp.cp.push([-1e5,-1e5],[1e5,-1e5],[1e5,1e5],[-1e5,1e5]);
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);
this.dp = [...this.dp, ...tp.dp];
}
inside(p) {
let int = 0; // find number of i ntersection points from p to far away
for (let i = 0, l = this.cp.length; i < l; i++) {
if (this.segment_intersect(p, [0.1, -1000], this.cp[i], this.cp[(i + 1) % l])) {
int++;
}
}
return int & 1; // if even your outside
}
boolean(p, diff = true) {
// bouding box optimization by ge1doot.
if (Math.abs(this.aabb[0] - p.aabb[0]) - (p.aabb[2] + this.aabb[2]) >= 0 &&
Math.abs(this.aabb[1] - p.aabb[1]) - (p.aabb[3] + this.aabb[3]) >= 0) return this.dp.length > 0;

// 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.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];
int.sort( (a,b) =>  (a[0] - ls0[0]) * cmpx + (a[1] - ls0[1]) * cmpy -
(b[0] - ls0[0]) * cmpx - (b[1] - ls0[1]) * cmpy);

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.001) {
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]);
}
}
}
}
}
return (this.dp = ndp).length > 0;
}
//port of http://paulbourke.net/geometry/pointlineplane/Helpers.cs
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;
}
};
const polygonList = new Node([0,0,200,200], 0, 4);
return {
list: polygonList,
create: () => new Polygon(),
draw: (turtle, p, addToVisList=true) => {
const vis = polygonList.boolean(p);
if (vis) {
p.draw(turtle);
}
}
};
}

// vec3 functions
const scale3=(a,b)=>[a[0]*b,a[1]*b,a[2]*b];
const len3=(a)=>Math.sqrt(dot3(a,a));
const normalize3=(a)=>scale3(a,1/len3(a));
const sub3=(a,b)=>[a[0]-b[0],a[1]-b[1],a[2]-b[2]];
const dot3=(a,b)=>a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
const cross3=(a,b)=>[a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0]]
// vec4 functions
const transform4=(a,b)=>{
const d=new Float32Array(4);
for(let c=0;4>c;c++)d[c]=b[c]*a[0]+b[c+4]*a[1]+b[c+8]*a[2]+b[c+12]*a[3];
return d;
}
// mat4 functions
const lookAt4m=(a,b,d)=>{ // pos, lookAt, up
const c=new Float32Array(16);
b=normalize3(sub3(a,b));
d=normalize3(cross3(d,b));
const e=normalize3(cross3(b,d));
c[0]=d[0];c[1]=e[0];c[2]=b[0];c[3]=0;
c[4]=d[1];c[5]=e[1];c[6]=b[1];c[7]=0;
c[8]=d[2];c[9]=e[2];c[10]=b[2];c[11]=0;
c[12]=-(d[0]*a[0]+d[1]*a[1]+d[2]*a[2]);
c[13]=-(e[0]*a[0]+e[1]*a[1]+e[2]*a[2]);
c[14]=-(b[0]*a[0]+b[1]*a[1]+b[2]*a[2]);
c[15]=1;
return c;
}
const multiply4m=(a,b)=>{
const d=new Float32Array(16);
for(let c=0;16>c;c+=4)
for(let e=0;4>e;e++)
d[c+e]=b[c+0]*a[0+e]+b[c+1]*a[4+e]+b[c+2]*a[8+e]+b[c+3]*a[12+e];
return d;
}
const perspective4m=(a,b,d)=>{ // fovy, aspect. near
const c=(new Float32Array(16)).fill(0,0);
c[5]=1/Math.tan(a/2);
c[0]=c[5]/b;
c[10]=c[11]=-1;
c[14]=-2*d;
return c;
}
```