3D SVG renderer for menger sponges

menger.py

import numpy as np
from svg import *

compose = lambda f: lambda g: lambda *a, **k: f(g(*a, **k))

PURPLE = '#6a1eb0'
ORANGE = '#ff824a'
BLUE   = '#75c1ff'

VP_WIDTH = 2.0
VP_HEIGHT = 2.0

VP_EYE = -np.array((11.90, 11, 27))
VP_DIR = -np.array((-2.00, -1, -3))
VP_FOC = 2.5

# I don't know numpy, there is probably a builtin thing for doing this
B3 = -VP_DIR;                B3 = B3 / np.linalg.norm(B3)
B2 = (0, 1, 0) - B3[2] * B3; B2 = B2 / np.linalg.norm(B2)
B1 = np.cross(B2, B3);       B1 = B1 / np.linalg.norm(B1)
B = np.matrix((B1, B2, B3)).transpose()
VIEWBASIS = B.I

def zcoordinate(u):
    x = VIEWBASIS * np.matrix(u - (VP_EYE + VP_DIR)).transpose()
    return np.ndarray.item(x[2]), np.ndarray.item(VP_FOC*x[0]/(1-x[2]), 0), np.ndarray.item(VP_FOC*x[1]/(1-x[2]), 0)

class View:
    def __init__(self, position, direction):
        self.pos = position
        self.dir = direction

        bz = -self.dir / np.linalg.norm(self.dir)
        by = (0, 1, 0) - bz[2] * bz
        by = by / np.linalg.norm(by)
        bx = np.cross(by, bz)

        self.basis = bx, by, bz
        self.invbasis = np.matrix(self.basis).transpose().I

    def global_to_view(self, v):
        v = v - (self.pos + self.dir)
        return self.invbasis * v.transpose()

# Basic header
svg_header()
svg_open(width=1024, height=1024,
         x_bias=-VP_WIDTH/2, y_bias=-VP_HEIGHT/2,
         xw=VP_WIDTH, yh=VP_HEIGHT)

class Poly:
    def __init__(self, *points, **opts):
        self.points = [np.array(p) for p in points]
        self.opts = opts
        centre = sum(self.points[:-1]) / len(points)
        zcoords = [zcoordinate(point) for point in self.points]
        self.projpoints = tuple((x, y) for _, x, y in zcoords)
        self.z = zcoordinate(centre)[0]

    def draw(self):
        svg_polygon(*self.projpoints, **self.opts)

    def __lt__(self, other):
        return self.z < other.z

class Cube:
    def __init__(self, x, y, z, w=1):
        self.faces = [
            Poly((x,y,z), (x+w,y,z), (x+w,y+w,z), (x,y+w,z), fill=PURPLE, stroke=PURPLE, stroke_width="0.001"),
            Poly((x,y,z), (x,y,z+w), (x,y+w,z+w), (x,y+w,z), fill=ORANGE, stroke=ORANGE, stroke_width="0.001"),
            Poly((x,y,z), (x+w,y,z), (x+w,y,z+w), (x,y,z+w), fill=BLUE, stroke=BLUE, stroke_width="0.001"),
            ]

    def projections(self):
        yield from self.faces

class Level:
    def __init__(self, x, y, z, n=1):
        if n == 1:
            f = Cube
        else:
            f = lambda *a, **k: Level(*a, **k, n=n-1)
        self.cubes = [
            # Back
            f(x+3**(n-1)*2, y+3**(n-1)*2, z+3**(n-1)*2), f(x+3**(n-1)*1, y+3**(n-1)*2, z+3**(n-1)*2), f(x+3**(n-1)*0, y+3**(n-1)*2, z+3**(n-1)*2),
            f(x+3**(n-1)*0, y+3**(n-1)*1, z+3**(n-1)*2), f(x+3**(n-1)*2, y+3**(n-1)*1, z+3**(n-1)*2),
            f(x+3**(n-1)*2, y+3**(n-1)*0, z+3**(n-1)*2), f(x+3**(n-1)*1, y+3**(n-1)*0, z+3**(n-1)*2), f(x+3**(n-1)*0, y+3**(n-1)*0, z+3**(n-1)*2),
            # Mid
            f(x+3**(n-1)*2, y+3**(n-1)*2, z+3**(n-1)*1), f(x+3**(n-1)*0, y+3**(n-1)*2, z+3**(n-1)*1),
            f(x+3**(n-1)*2, y+3**(n-1)*0, z+3**(n-1)*1), f(x+3**(n-1)*0, y+3**(n-1)*0, z+3**(n-1)*1),
            # Front
            f(x+3**(n-1)*2, y+3**(n-1)*2, z+3**(n-1)*0), f(x+3**(n-1)*1, y+3**(n-1)*2, z+3**(n-1)*0), f(x+3**(n-1)*0, y+3**(n-1)*2, z+3**(n-1)*0),
            f(x+3**(n-1)*0, y+3**(n-1)*1, z+3**(n-1)*0), f(x+3**(n-1)*2, y+3**(n-1)*1, z+3**(n-1)*0),
            f(x+3**(n-1)*2, y+3**(n-1)*0, z+3**(n-1)*0), f(x+3**(n-1)*1, y+3**(n-1)*0, z+3**(n-1)*0), f(x+3**(n-1)*0, y+3**(n-1)*0, z+3**(n-1)*0),
            ]

    def projections(self):
        for cube in self.cubes:
            yield from cube.projections()

def draw(obj):
    polygons = sorted(obj.projections())
    for polygon in polygons:
        polygon.draw()

draw(Level(0, 0, 0, n=2))

# Footer
svg_close()

menger.svg

svg.py

def xml(tag, _xml_tag_is_a_singleton=True, **options):
    s = f'<{tag}'
    kw_attrib = lambda x: x.replace('_', '-')
    if options:
        s += ' '
        s += ' '.join(f'{kw_attrib(k)}="{str(v)}"' for k, v in options.items())
    if _xml_tag_is_a_singleton:
        s += ' />'
    else:
        s += '>'
    print(s)

def xml_open(*args, **kwargs):
    xml(*args, **kwargs, _xml_tag_is_a_singleton=False)

def xml_close(tag):
    print(f'</{tag}>')

def svg_header():
    '<?xml version="1.0" encoding="UTF-8" standalone="no"?>'

def svg_open(width, height, x_bias, y_bias, xw, yh, **opts):
    vb = f'{x_bias} {y_bias} {xw} {yh}'
    ns = 'http://www.w3.org/2000/svg'
    xml_open('svg', width=width, height=height, viewBox=vb, xmlns=ns, **opts)

def svg_close():
    xml_close('svg')

def svg_poly(*points, **opts):
    point_str = ' '.join(f'{x},{y}' for x, y in points)
    xml('polyline', points=point_str, **opts)

def svg_polygon(*points, **opts):
    point_str = ' '.join(f'{x},{y}' for x, y in points)
    xml('polygon', points=point_str, **opts)

def svg_circle(point, radius, **opts):
    xml('circle', cx=point[0], cy=point[1], r=radius, **opts)

def svg_line(p1, p2, **opts):
    xml('line', x1=p1[0], y1=p1[1], x2=p2[0], y2=p2[1], **opts)

def svg_rect(p, width, height, x_radius=0, y_radius=0, **opts):
    xml('rect', x=p[0], y=p[1], width=width, height=height,
              rx=x_radius, ry=y_radius, **opts)