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

Log in to post a comment.

// 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) {
            turtle.addTransform(EscherDroste);
        }
        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();
                p.addPoints(...cp);
                p.addOutline(i);
                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 = [];
        }
        addTransform(t) {
            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);
	    }
	    add(p) {
	        for (let i=0; i<this.children.length; i++) {
	            if (this.children[i].inside(p.aabb)) {
	                this.children[i].add(p);
	                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
		}
		addPoints(...points) {
		    // 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];
		}
		addSegments(...points) {
		    // add segments (each a pair of points)
		    points.forEach(p => this.dp.push(p));
		}
		addOutline() {
			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]);
			}
		}
		addHatching(a, d) {
			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);
			    if (addToVisList) polygonList.add(p);
			}
		}
	};
}

// 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 add3=(a,b)=>[a[0]+b[0],a[1]+b[1],a[2]+b[2]];
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;
}