// Copyright 2016-2021 Doug Moen
// Licensed under the Apache License, version 2.0
// See accompanying file LICENSE or https://www.apache.org/licenses/LICENSE-2.0

/* Curv standard library, automatically loaded at startup. */
{

lib = file "lib";

null = #null;
deg = tau/360;
e = 2.71828182845904523536028747135266249775724709369995;
phi = sqrt 5 * .5 + .5;
X = 0;
Y = 1;
Z = 2;
T = 3;
X_axis = [1,0,0];
Y_axis = [0,1,0];
Z_axis = [0,0,1];
nl = char 0xA;   // newline character
tab = char 0x9;  // tab character
dol = #"$";  // dollar sign character
quot = #"""; // double quote character

mod[a,m] = a - m * floor(a/m);
rem[a,m] = a - m * trunc(a/m);
lerp[a,b,t] = a*(1-t) + b*t;
// Smooth Hermite interpolation between 0 and 1 when lo < x < hi.
// Used to construct a threshold function with a smooth transition.
// Results are undefined if lo >= hi.
// https://en.wikipedia.org/wiki/Smoothstep
smoothstep[lo,hi,x] =
    let t = clamp[(x - lo)/(hi - lo), 0, 1];
    in t*t*(3 - 2*t);
clamp[v,lo,hi] = min[max[v,lo],hi];

sec a = 1 / sin a;
csc a = 1 / cos a;
cot a = cos a / sin a;
atan2[y,x] = phase[x,y];

// lists
//String = {T: List Char, call: string};
encode = ucode; // backward compatibility to Curv 0.4
decode = char; // backward compatibility to Curv 0.4
strcat L = concat (map string L); // back compat to Curv 0.4
concat vv = [for (v in vv) for (i in v) i];
reverse v = v.[count(v)-1..0 by -1];
map f list = [for (x in list) f x];
filter p list = [for (x in list) if (p x) x];
reduce [zero, f] list =
    if (list == [])
        zero
    else
        do
            local r = list.[0];
            for (i in 1..<count list)
                r := f[r, list.[i]];
        in r;
//sum = reduce[0, [x,y]->x+y];
product = reduce[1, [x,y]->x*y];
sort a =
    if (a == [])
        []
    else
        let first = a.[0];
            rest = a.[1..<count a];
        in [...sort[for (e in rest) if (e < first) e],
            first,
            ...sort[for (e in rest) if (e >= first) e]];
contains [list, x] =
    do
        local i = 0;
        while (i < count list && list.[i] != x)
            i := i + 1;
    in i < count list;

// functions
id x = x;
is_fun x =
    do warning
        "'is_fun' is deprecated. Use is_func or is_primitive_func instead.";
    in is_primitive_func x;
into f rest first = f [first, ...rest];

merge rs = {for (r in rs) ...r};

ensure pred (x :: pred) = x;

// complex numbers: [RE,IM]
RE=0;
IM=1;
cmul [z,w] = [z.[RE]*w.[RE] - z.[IM]*w.[IM], z.[IM]*w.[RE] + z.[RE]*w.[IM]];
csqr z = [ z.[RE]*z.[RE] - z.[IM]*z.[IM], 2*z.[RE]*z.[IM] ];

////////////////////
// Linear Algebra //
////////////////////
is_vec2 x = is_list x && count x == 2 && is_num(x.[0]) && is_num(x.[1]);
is_vec3 x = is_list x && count x == 3 && is_num(x.[0]) && is_num(x.[1]) && is_num(x.[2]);
indices a = 0..<count a;
cross[p,q] = [p.[Y]*q.[Z] - p.[Z]*q.[Y], p.[Z]*q.[X] - p.[X]*q.[Z], p.[X]*q.[Y] - p.[Y]*q.[X]];
idmatrix n = [for(i in 1..n) [for(j in 1..n) if(i==j) 1 else 0]];
transpose a =
    if (a == []) a
    else [for (i in indices(a.[0])) [for (j in indices a) a.[j,i]]];
normalize v = v / mag v;

// convert phase angle to unit vector
cis theta = [cos theta, sin theta];

// perp: Rotate a 2D point by 90 degrees CCW. Multiply a complex number by i.
// It's the 2D analog of the 3D vector cross product (Cross in Mathematica).
// dot[perp a, b] is the "perp-dot" product:
// see: 'The Pleasures of "Perp-Dot" Products', Graphics Gems IV.
perp [x,y] = [-y, x];

// angle between two vectors. Deprecated.
//angle[a,b] = acos(dot[a,b] / (mag a * mag b));

////////////
// SHAPES //
////////////
make_shape r =
    let shape = {
            is_2d : false,
            is_3d : false,
            bbox : [[-inf,-inf,-inf],[inf,inf,inf]],
            colour p : sRGB[.8, .8, .5],
            ... r,
        };
    in do
        assert (shape.is_2d || shape.is_3d);
        assert (is_bbox3(shape.bbox));
        assert (is_func(shape.colour));
        assert (defined(shape.dist) && is_func(shape.dist));
    in shape;

is_shape x =
    defined(x.is_2d) && defined(x.is_3d) && defined(x.bbox)
    && defined(x.dist) && defined(x.colour);

size shape = shape.bbox.[MAX] - shape.bbox.[MIN];

MIN = 0;
MAX = 1;

is_bbox2 = match [
    [_ :: is_vec2, _ :: is_vec2] -> true;
    _ -> false;
];
is_bbox3 = match [
    [_ :: is_vec3, _ :: is_vec3] -> true;
    _ -> false;
];

/////////////////////
// CIRCULAR SHAPES //
/////////////////////

circle =
    let circle_e d =
            let r = d / 2;
            in make_shape {
                dist[x,y,z,t] = mag[x,y] - r;
                bbox = [[-r,-r,0],[r,r,0]];
                is_2d = true;
            };
    in make_shape {
        ... circle_e 2,
        call : circle_e,
    };

ellipse[dx,dy] = circle 1 >> stretch[dx,dy];

sphere =
    let sphere_e d =
            let r = d / 2;
            in make_shape {
                dist[x,y,z,t] = mag[x,y,z] - r;
                bbox = [[-r,-r,-r],[r,r,r]];
                is_3d = true;
            };
    in make_shape {
        ... sphere_e 2,
        call : sphere_e,
    };

ellipsoid [dx,dy,dz] = sphere 1 >> stretch[dx,dy,dz];

cylinder =
    let mitred{d,h} = extrude.mitred h (circle d);
        exact{d,h} = extrude.exact h (circle d);
    in {
        ... exact{d:2,h:2};
        call : exact;
        mitred : mitred;
        exact : exact;
    };

// Cone: +Z is up, origin is in the middle of the base.
// TODO: Provide an interface capable of constructing an infinite cone.
cone =
    {
        call = exact;
        mitred{d, h: height} =
            let radius = d / 2;
            in make_shape {
                dist[x,y,z,t] =
                    let q = [mag[x,y], z];
                        tip = q - [0, height];
                        mantleDir = normalize[height, radius];
                        mantle = dot[tip, mantleDir];
                    in max[mantle, -q.[Y]];
                bbox = [[-radius,-radius,0],[radius,radius,height]];
                is_3d = true
            };
        // Euclidean distance field.
        // Credit: MERCURY (http://mercury.sexy/hg_sdf)
        exact{d, h: height} =
            let radius = d / 2;
            in make_shape {
                dist[x,y,z,t] =
                    let q = [mag[x,y], z];
                        apex = q - [0, height];
                        mantleDir = normalize[height, radius];
                        mantle = dot[apex, mantleDir];
                        projected = dot[apex, [mantleDir.[Y], -mantleDir.[X]]];
                        d = max[mantle, -q.[Y]];
                    in do
                        // distance to apex
                        if (q.[Y] > height && projected < 0)
                            d := max[d, mag(apex)];
                        // distance to base ring
                        if (q.[X] > radius && projected > mag[height, radius])
                            d := max[d, mag(q - [radius, 0])];
                    in d;
                bbox = [[-radius,-radius,0],[radius,radius,height]];
                is_3d = true
            };
    };

// credit: based on sdCappedCone by Inigo Quilez
capped_cone {h,top,bottom} =
  let
    hh = h / 2;
    r1 = bottom / 2;
    r2 = top / 2;
    k1 = [r2, hh];
    k2 = [r2-r1, h];
  in make_shape {
    dist p = do
      local q = [ mag(p.[[X,Y]]), p.[Z] ];
      local ca = [q.[X]-min[q.[X],if (q.[Y]<0) r1 else r2], abs(q.[Y])-hh];
      local cb = q - k1 + k2*clamp[ dot[k1-q,k2]/dot[k2,k2], 0, 1 ];
      local s = if (cb.[X]<0 && ca.[Y]<0) -1 else 1;
    in
      s*sqrt( min[dot[ca,ca],dot[cb,cb]] );
    bbox = let r = max[r1,r2] in [[-r,-r,-hh],[r,r,hh]];
    is_3d = true;
  };

// credit: Inigo Quilez
capsule{d, from: a, to: b} =
    let ba = b - a;
        r = d / 2;
    in make_shape {
        dist p =
            let pa = p.[[X,Y,Z]] - a;
                h = clamp[dot[pa,ba] / dot[ba,ba], 0, 1];
            in mag(pa - ba*h) - r;
        bbox = [min[a-r,b-r], max[a+r,b+r]];
        is_3d = true;
    };

stroke{d, from: a, to: b} =
    let r = d / 2;
        ba = b - a;
    in make_shape {
        dist p =
            let pa = p.[[X,Y]] - a;
                h = clamp[dot[pa,ba] / dot[ba,ba], 0, 1];
            in mag(pa - ba*h) - r;
        bbox = [[...min[a-r,b-r],0], [...max[a+r,b+r],0]];
        is_2d = true;
    };

polyline {
        d :: (a->is_num a && a >= 0),
        v :: (a->is_list a && count a > 0 && and(map is_vec2 a))
    } =
    let r = d / 2;
    in make_shape {
        dist p =
            do
                local p = p.[[X,Y]];
                local num = count v;
                local d = dot[p-v.[0], p-v.[0]];
                local s = 1;
                local j = 0;
                for (i in 1..<num) (
                    local e = v.[j] - v.[i];
                    local w = p - v.[i];
                    local b = w - e*clamp[dot[w,e]/dot[e,e], 0, 1];
                    d := min[d, dot[b,b]];
                    j := i;
                );
            in s * sqrt d - r;
        bbox = [[...min v-r,0], [...max v+r,0]];
        is_2d = true;
    };

// d1: major diameter, from center of tube, through origin, to centre of tube.
// d2: minor diameter, of the tube.
torus {major:d1, minor:d2} = perimeter_extrude (circle d1) (circle d2);

// from OpenSCAD
revolve shape = perimeter_extrude (circle 0) shape;

// A generalized torus. Sweep a 2D shape around the perimeter of another
// 2D shape.
perimeter_extrude perimeter cross_section =
    make_shape {
        dist[x,y,z,t] = cross_section.dist[ perimeter.dist[x,y,0,t], z, 0, t ];
        bbox = [
            [ perimeter.bbox.[MIN,X] - cross_section.bbox.[MAX,X],
              perimeter.bbox.[MIN,Y] - cross_section.bbox.[MAX,X],
              cross_section.bbox.[MIN,Y] ],
            [ perimeter.bbox.[MAX,X] + cross_section.bbox.[MAX,X],
              perimeter.bbox.[MAX,Y] + cross_section.bbox.[MAX,X],
              cross_section.bbox.[MAX,Y] ],
        ];
        is_3d = true;
    };

///////////////
// POLYTOPES //
///////////////

half_space =
    let impl[a,n] = // a is point on boundary, n is unit normal
        make_shape {
            dist[x,y,z,_] = dot[-n, a-[x,y,z]];
            bbox = [ [ if (n.[Y]==0 && n.[Z]==0 && n.[X]<0) a.[X] else -inf,
                       if (n.[X]==0 && n.[Z]==0 && n.[Y]<0) a.[Y] else -inf,
                       if (n.[X]==0 && n.[Y]==0 && n.[Z]<0) a.[Z] else -inf, ],
                     [ if (n.[Y]==0 && n.[Z]==0 && n.[X]>0) a.[X] else inf,
                       if (n.[X]==0 && n.[Z]==0 && n.[Y]>0) a.[Y] else inf,
                       if (n.[X]==0 && n.[Y]==0 && n.[Z]>0) a.[Z] else inf, ] ];
            is_3d = true;
        };
    in match [
        {d, normal: n} -> impl[d*n,n];
        [a, b, c] -> impl[a, normalize(cross[b-a,c-b])];
        {at: a, normal: n} -> impl[a,n];
    ];

half_plane =
    let impl[a,n] = // a is point on boundary, n is unit normal
        make_shape {
            dist[x,y,_,_] = dot[-n, a-[x,y]];
            bbox = [ [ if (n.[Y]==0 && n.[X]<0) a.[X] else -inf,
                       if (n.[X]==0 && n.[Y]<0) a.[Y] else -inf,
                       0, ],
                     [ if (n.[Y]==0 && n.[X]>0) a.[X] else inf,
                       if (n.[X]==0 && n.[Y]>0) a.[Y] else inf,
                       0, ], ];
            is_2d = true;
        };
    in match [
        {d, normal: n} -> impl[d*n,n];
        [a, b] -> impl[a, normalize(perp(b-a))];
        {at: a, normal: n} -> impl[a,n];
    ];

square =
    let
        square_m d =
            let r = d / 2;
            in make_shape {
                dist[x,y,z,t] = max(abs[x,y] - r);
                bbox = [[-r,-r,0], [r,r,0]];
                is_2d = true;
            };
        square_e d = rect.exact[d,d];
    in {
        ... square_e 2,
        call : square_e,
        circumratio : sqrt 2,
        exact : {
            ... square_e 2,
            call : square_e,
        },
        mitred : {
            ... square_m 2,
            call : square_m,
        },
    };

// rect: an axis aligned rectangle.
rect = {
    call = exact;
    mitred x =
        x >> match [
        b :: is_vec2 ->
            let r = b / 2;
            in make_shape {
                dist[x,y,z,t] = max(abs[x,y] - r);
                bbox = [[-r.[X],-r.[Y],0],[r.[X],r.[Y],0]];
                is_2d = true;
            };
        b :: is_bbox2 ->
            let xmin = b.[MIN,X];
                ymin = b.[MIN,Y];
                xmax = b.[MAX,X];
                ymax = b.[MAX,Y];
            in intersection [
                if (xmin > -inf) make_shape {
                    dist[x,y,z,t] = xmin - x;
                    bbox = [[xmin, -inf, 0], [inf, inf, 0]];
                    is_2d = true;
                },
                if (ymin > -inf) make_shape {
                    dist[x,y,z,t] = ymin - y;
                    bbox = [[-inf, ymin, 0], [inf, inf, 0]];
                    is_2d = true;
                },
                if (xmax < inf) make_shape {
                    dist[x,y,z,t] = x - xmax;
                    bbox = [[-inf, -inf, 0], [xmax, inf, 0]];
                    is_2d = true;
                },
                if (ymax < inf) make_shape {
                    dist[x,y,z,t] = y - ymax;
                    bbox = [[-inf, -inf, 0], [inf, ymax, 0]];
                    is_2d = true;
                },
            ];
        {xmin=-inf, ymin=-inf, xmax=inf, ymax=inf} ->
            mitred[[xmin,ymin],[xmax,ymax]];
        ];
    exact x =
        x >> match [
        b :: is_vec2 ->
            let r = b / 2;
            in make_shape {
                dist[x,y,z,t] =
                    let d = abs[x,y] - r;
                    in min[max d, 0] + mag(max[d,0]);
                bbox = [[-r.[X],-r.[Y],0],[r.[X],r.[Y],0]];
                is_2d = true;
            };
        b :: is_bbox2 ->
            let lo = b.[MIN];
                hi = b.[MAX];
                outsideness[a,amin,amax] =
                    if (a < amin) amin - a
                    else if (a > amax) a - amax
                    else 0;
            in make_shape {
                dist[x,y,z,t] =
                    let ind = max[max(lo-[x,y]), max([x,y]-hi)];
                        outd = mag[outsideness[x,lo.[X],hi.[X]],
                                   outsideness[y,lo.[Y],hi.[Y]]];
                    in if (ind <= 0) ind else outd;
                bbox = [[...lo,0], [...hi,0]];
                is_2d = true;
            };
        {xmin=-inf, ymin=-inf, xmax=inf, ymax=inf} ->
            exact[[xmin,ymin],[xmax,ymax]];
        ];
};

// regular_polygon: n is # of sides, d is incircle diameter
regular_polygon n =
    let angle = tau / n;
        mitred d = half_plane{d: d/2, normal: [0,-1]}
                   >> repeat_radial n;
    in {
        ... mitred 2,
        call: mitred,
        mitred: mitred,
        circumratio : 1 / cos(angle/2),
    };

convex_polygon pts =
    do warning "'convex_polygon' is deprecated. Use 'polygon' instead.";
    in polygon pts;

// Finite polygon. May be nonconvex. Edges may cross, and the polygon may
// self intersect. Uses the "crossing number" method to determine which
// points are inside/outside the polygon. It doesn't matter if the vertexes
// are listed in clockwise or counterclockwise order.
// (https://web.archive.org/web/20200117192710/http://geomalgorithms.com/a03-_inclusion.html)
// Exact distance field.
// From https://www.shadertoy.com/view/wdBXRW
// which is Copyright 2019 Inigo Quilez (The MIT License)
polygon v = make_shape {
  dist p =
    do
      local p = p.[[X,Y]];
      local num = count v;
      local d = dot[p-v.[0], p-v.[0]];
      local s = 1;
      local j = num-1;
      for (i in 0..<num) (
        local e = v.[j] - v.[i];
        local w = p - v.[i];
        local b = w - e*clamp[dot[w,e]/dot[e,e], 0, 1];
        d := min[d, dot[b,b]];
        local cond = [p.[Y] >= v.[i,Y], p.[Y] < v.[j,Y], e.[X]*w.[Y] > e.[Y]*w.[X]];
        if (and cond || and(not cond)) s := -s;
        j := i;
      );
    in s * sqrt d;
  bbox = [[...min v,0], [...max v,0]];
  is_2d = true;
};

// box: an axis aligned cuboid.
box =
    let mitred a = a >> match [
            b :: is_vec3 ->
                let r = b / 2;
                in make_shape {
                    dist[x,y,z,t] = max(abs[x,y,z] - r);
                    bbox = [-r, r];
                    is_3d = true;
                };
            b :: is_bbox3 ->
                let lo = b.[MIN];
                    hi = b.[MAX];
                in make_shape {
                    dist[x,y,z,t] = max[max(lo-[x,y,z]), max([x,y,z]-hi)];
                    bbox = b;
                    is_3d = true;
                };
            {xmin=-inf, ymin=-inf, zmin=-inf, xmax=inf, ymax=inf, zmax=inf} ->
                mitred[[xmin,ymin,zmin],[xmax,ymax,zmax]];
            ];
        exact a = a >> match [
            b :: is_vec3 ->
                let r = b / 2;
                in make_shape {
                    dist[x,y,z,t] =
                        let d = abs[x,y,z] - r;
                        in min[max d, 0] + mag(max[d,0]);
                    bbox = [-r, r];
                    is_3d = true;
                };
            b :: is_bbox3 ->
                let lo = b.[MIN];
                    hi = b.[MAX];
                    outsideness[a,amin,amax] =
                        if (a < amin) amin - a
                        else if (a > amax) a - amax
                        else 0;
                in make_shape {
                    dist[x,y,z,t] =
                        let ind = max[max(lo-[x,y,z]), max([x,y,z]-hi)];
                            outd = mag[outsideness[x,lo.[X],hi.[X]],
                                       outsideness[y,lo.[Y],hi.[Y]],
                                       outsideness[z,lo.[Z],hi.[Z]]];
                        in if (ind <= 0) ind else outd;
                    bbox = b;
                    is_3d = true;
                };
            {xmin=-inf, ymin=-inf, zmin=-inf, xmax=inf, ymax=inf, zmax=inf} ->
                exact[[xmin,ymin,zmin],[xmax,ymax,zmax]];
            ];
    in {
        call : exact;
        mitred : mitred;
        exact : exact;
        ... exact [2,2,2];
    };

prism n =
    let prism_m {d,h} = regular_polygon n .mitred d >> extrude.mitred h;
    in {
        ... prism_m {d: 2, h: 2},
        call : prism_m,
        mitred : {
            ... prism_m {d: 2, h: 2},
            call : prism_m,
        },
    };

// TODO: pyramid[n,d,h]

tetrahedron =
    let tetrahedron_m d =
            let r = d / 2;
                // vertices of a tetrahedron
                v1 = [1,1,1];
                v2 = [-1,-1,1];
                v3 = [1,-1,-1];
                v4 = [-1,1,-1];

                // faces: [1,2,3] [1,4,2] [2,4,3] [1,3,4]
                // vertices in CCW order viewed from outside (right hand rule)

                // face normals
                n1 = normalize(cross[v1-v3,v2-v3]);
                n2 = normalize(cross[v1-v2,v4-v2]);
                n3 = normalize(cross[v2-v3,v4-v3]);
                n4 = normalize(cross[v4-v3,v1-v3]);
            in make_shape {
                dist p :
                    max [
                        dot[p.[[X,Y,Z]], n1],
                        dot[p.[[X,Y,Z]], n2],
                        dot[p.[[X,Y,Z]], n3],
                        dot[p.[[X,Y,Z]], n4],
                    ] - r,
                is_3d : true,
                bbox : sqrt(3) * [[-r,-r,-r], [r,r,r]],
            };
    in {
        ... tetrahedron_m 2,
        call : tetrahedron_m,
        mitred : {
            ... tetrahedron_m 2,
            call : tetrahedron_m,
        },
        // radius of circumscribed sphere / radius of inscribed sphere
        circumratio : (sqrt(3/2)/2) / (1/(2*sqrt 6)),
    };

cube =
    let
        cube_m d =
            let r = d / 2;
            in make_shape {
                dist[x,y,z,t] = max(abs[x,y,z] - r);
                bbox = [[-r, -r, -r], [r,r,r]];
                is_3d = true;
            };
        cube_e d = box.exact[d,d,d];
    in {
        ... cube_e 2,
        call: cube_e,
        circumratio: sqrt 3,
        exact: {
            ... cube_e 2,
            call: cube_e
        },
        mitred: {
            ... cube_m 2,
            call: cube_m
        },
    };

octahedron =
    let
        octahedron_m d =
            let r = d / 2;
                v0 = normalize[1,1,1];
                v1 = normalize[-1,1,1];
                v2 = normalize[1,-1,1];
                v3 = normalize[1,1,-1];
                // radius of circumscribed sphere / radius of inscribed sphere
                circumratio = (1/2*sqrt 2) / (1/6*sqrt 6);
            in make_shape {
                dist p :
                    max [
                        abs(dot[p.[[X,Y,Z]], v0]),
                        abs(dot[p.[[X,Y,Z]], v1]),
                        abs(dot[p.[[X,Y,Z]], v2]),
                        abs(dot[p.[[X,Y,Z]], v3]),
                    ] - r,
                is_3d : true,
                bbox : circumratio * [[-r,-r,-r], [r,r,r]],
            };
    in {
        ... octahedron_m 2,
        call : octahedron_m,
        circumratio : (1/2*sqrt 2) / (1/6*sqrt 6),
        mitred: {
            ... octahedron_m 2,
            call: octahedron_m
        },
    };

dodecahedron =
    let
        v0 = normalize[0, phi, 1];
        v1 = normalize[0, -phi, 1];
        v2 = normalize[1, 0, phi];
        v3 = normalize[-1, 0, phi];
        v4 = normalize[phi, 1, 0];
        v5 = normalize[-phi, 1, 0];
        midratio = (.25*(3+sqrt(5))) / (1/20*sqrt(250+110*sqrt(5)));
        circumratio = (.25*(sqrt 3 + sqrt 15)) / (1/20*sqrt(250+110*sqrt(5)));
        dodecahedron_m d =
            let r = d / 2;
            in make_shape {
                dist p :
                    max [
                        abs(dot[p.[[X,Y,Z]], v0]),
                        abs(dot[p.[[X,Y,Z]], v1]),
                        abs(dot[p.[[X,Y,Z]], v2]),
                        abs(dot[p.[[X,Y,Z]], v3]),
                        abs(dot[p.[[X,Y,Z]], v4]),
                        abs(dot[p.[[X,Y,Z]], v5]),
                    ] - r,
                is_3d : true,
                bbox : midratio * [[-r,-r,-r],[r,r,r]],
            };
    in {
        ... dodecahedron_m 2,
        call : dodecahedron_m,
        circumratio : (.25*(sqrt 3 + sqrt 15)) / (1/20*sqrt(250+110*sqrt(5))),
        mitred: {
            ... dodecahedron_m 2,
            call: dodecahedron_m
        },
    };

icosahedron =
    let
        v0 = normalize[1,1,1];
        v1 = normalize[-1,1,1];
        v2 = normalize[1,-1,1];
        v3 = normalize[1,1,-1];
        v4 = normalize[0, +1, phi+1];
        v5 = normalize[0, -1, phi+1];
        v6 = normalize[+phi+1, 0, 1];
        v7 = normalize[-phi-1, 0, 1];
        v8 = normalize[+1, phi+1, 0];
        v9 = normalize[-1, phi+1, 0];
        // radius of midsphere / radius of insphere
        midratio = (phi/2) / (1/12*(3*sqrt(3)+sqrt(15)));
        circumratio = (.25*sqrt(10+2*sqrt 5)) / (1/12*(3*sqrt(3)+sqrt(15)));
        icosahedron_m d =
            let r = d / 2;
            in make_shape {
                dist p :
                    max [
                        abs(dot[p.[[X,Y,Z]], v0]),
                        abs(dot[p.[[X,Y,Z]], v1]),
                        abs(dot[p.[[X,Y,Z]], v2]),
                        abs(dot[p.[[X,Y,Z]], v3]),
                        abs(dot[p.[[X,Y,Z]], v4]),
                        abs(dot[p.[[X,Y,Z]], v5]),
                        abs(dot[p.[[X,Y,Z]], v6]),
                        abs(dot[p.[[X,Y,Z]], v7]),
                        abs(dot[p.[[X,Y,Z]], v8]),
                        abs(dot[p.[[X,Y,Z]], v9]),
                    ] - r,
                is_3d : true,
                bbox : midratio * [[-r,-r,-r], [r,r,r]],
            };
    in {
        ... icosahedron_m 2,
        call : icosahedron_m,
        circumratio : (.25*sqrt(10+2*sqrt 5)) / (1/12*(3*sqrt(3)+sqrt(15))),
        mitred: {
            ... icosahedron_m 2,
            call: icosahedron_m
        },
    };

/////////////////////////////////////////////
// BOOLEAN / SET THEORETIC SHAPE OPERATORS //
/////////////////////////////////////////////

nothing =
    make_shape {
        dist p : inf,
        bbox : [[inf,inf,inf],[-inf,-inf,-inf]],
        is_2d : true,
        is_3d : true,
    };
everything =
    make_shape {
        dist p : -inf,
        bbox : [[-inf,-inf,-inf],[inf,inf,inf]],
        is_2d : true,
        is_3d : true,
    };
complement s =
    make_shape {
        dist p : -s.dist p,
        bbox : [[-inf,-inf,-inf],[inf,inf,inf]],
        is_2d : s.is_2d,
        is_3d : s.is_3d,
    };

// Approximate union that produces a mitred SDF inside. Fast.
// When unioning a list of coloured shapes, we paint the shapes from first to
// last order: the last shape is painted on top of its predecessors.
union list = reduce[nothing, _union2] list;
_union2 [s1,s2] =
    make_shape {
        dist p : min[s1.dist p, s2.dist p],
        colour p :
            let d1 = s1.dist p;
                d2 = s2.dist p;
            in if (d2 <= 0 || d2 <= d1) s2.colour p else s1.colour p,
        bbox : [min[s1.bbox.[MIN], s2.bbox.[MIN]], max[s1.bbox.[MAX], s2.bbox.[MAX]]],
        is_2d : s1.is_2d && s2.is_2d,
        is_3d : s1.is_3d && s2.is_3d,
    };

intersection list = reduce[everything, _intersection2] list;
_intersection2 [s1,s2] =
    make_shape {
        dist p : max[s1.dist p, s2.dist p],
        colour : s1.colour,
        bbox : [max[s1.bbox.[MIN], s2.bbox.[MIN]], min[s1.bbox.[MAX], s2.bbox.[MAX]]],
        is_2d : s1.is_2d && s2.is_2d,
        is_3d : s1.is_3d && s2.is_3d,
    };

difference[s1,s2] = _intersection2[s1, complement s2];
symmetric_difference shapes = difference[union shapes, intersection shapes];

row = let
    row_impl[d,shapes] =
        let widths = [for (s in shapes) size s .[X]];
            gap = if (d==null) max widths / 4 else d;
            width = sum widths + (count shapes - 1)*gap;
            xdelta = [];
            x = -width / 2;
            i = 0;
        in do 
            while (i < count shapes) (
                xdelta := concat[xdelta, [x - shapes.[i].bbox.[MIN,X]]];
                x := x + gap + size(shapes.[i]).[X];
                i := i + 1;
            );
        in union[for (i in 0..<count shapes) shapes.[i] >> move[xdelta.[i],0,0]];
in match [
    d :: is_num -> shapes -> row_impl[d,shapes];
    shapes -> row_impl[null, shapes];
];

smooth k = {
    union [s1,s2] =
        make_shape {
            dist p : smooth_min[s1.dist p, s2.dist p, k],
            colour p :
                let d1 = s1.dist p;
                    d2 = s2.dist p;
                in if (d2 <= 0 || d2 <= d1) s2.colour p else s1.colour p,
            bbox : [ min[s1.bbox.[MIN], s2.bbox.[MIN]] - k/4,
                     max[s1.bbox.[MAX], s2.bbox.[MAX]] + k/4 ],
            is_2d : s1.is_2d && s2.is_2d,
            is_3d : s1.is_3d && s2.is_3d,
        };

    intersection [s1,s2] =
        make_shape {
            dist p : smooth_max[s1.dist p, s2.dist p, k],
            colour : s1.colour,
            bbox : [ max[s1.bbox.[MIN], s2.bbox.[MIN]],
                     min[s1.bbox.[MAX], s2.bbox.[MAX]] ],
            is_2d : s1.is_2d && s2.is_2d,
            is_3d : s1.is_3d && s2.is_3d,
        };

    difference [s1,s2] = intersection [s1, complement s2];
};

// polynomial smooth min, from IQ
// http://www.iquilezles.org/www/articles/smin/smin.htm
smooth_min[a,b,k] =
    if (a == inf)
        // TODO: work around bug in 'max' that exists on ubuntu/nvidia GTX 1050
        // Bug is: max[-inf,x] returns NaN.
        b
    else
        let h = clamp[ 0.5+0.5*(b-a)/k, 0, 1 ];
        in lerp[ b, a, h ] - k*h*(1.0-h);
smooth_max[a,b,k] = -smooth_min[-a,-b,k];

chamfer r = {
    union [s1,s2] =
        make_shape {
            dist p =
                let a = s1.dist p;
                    b = s2.dist p;
                in chamfer_min[a,b,r];
            colour p =
                let d1 = s1.dist p;
                    d2 = s2.dist p;
                in if (d2 <= 0 || d2 <= d1) s2.colour p else s1.colour p;
            bbox = [ min[s1.bbox.[MIN], s2.bbox.[MIN]],
                     max[s1.bbox.[MAX], s2.bbox.[MAX]] ];
            is_2d = s1.is_2d && s2.is_2d;
            is_3d = s1.is_3d && s2.is_3d;
        };

    intersection [s1,s2] =
        make_shape {
            dist p =
                let a = s1.dist p;
                    b = s2.dist p;
                in chamfer_max[a,b,r];
            colour = s1.colour;
            bbox = [ max[s1.bbox.[MIN], s2.bbox.[MIN]],
                     min[s1.bbox.[MAX], s2.bbox.[MAX]] ];
            is_2d = s1.is_2d && s2.is_2d;
            is_3d = s1.is_3d && s2.is_3d;
        };

    difference [s1,s2] = intersection [s1, complement s2];

    chamfer_min[a, b, r] =
        let e = max[r - abs(a - b), 0];
        in min[a, b] - e*.5;
    chamfer_max[a,b,k] = -chamfer_min[-a,-b,k];
};

////////////////////////////
// AFFINE TRANSFORMATIONS //
////////////////////////////

move = translate;
translate _delta shape =
    let delta = if (count _delta == 2) [_delta.[X],_delta.[Y],0] else _delta;
    in do
        if (shape.is_2d) assert(delta.[Z] == 0);
    in make_shape {
        // TODO: dist[x,y,z,t] : shape.dist[...[x,y,z]-delta, t],
        // TODO: dist p : shape.dist[...p.[[X,Y,Z]]-delta, p.[T]],
        dist p : shape.dist[p.[X]-delta.[X], p.[Y]-delta.[Y], p.[Z]-delta.[Z], p.[T]],
        colour p : shape.colour[p.[X]-delta.[X],p.[Y]-delta.[Y],p.[Z]-delta.[Z],p.[T]],
        bbox : [shape.bbox.[MIN]+delta, shape.bbox.[MAX]+delta],
        is_2d : shape.is_2d,
        is_3d : shape.is_3d,
    };

// Isotropic scale operation, preserves structure of the distance field.
scale (s :: is_num) shape =
    make_shape {
        // TODO: dist p : shape.dist[...p.[[X,Y,Z]] / s, p.[T]] * s,
        dist[x,y,z,t] : shape.dist[x/s, y/s, z/s, t] * s,
        colour[x,y,z,t] : shape.colour[x/s, y/s, z/s, t],
        bbox : [s*shape.bbox.[MIN], s*shape.bbox.[MAX]],
        is_2d : shape.is_2d,
        is_3d : shape.is_3d,
    };

// Anisotropic scale operation, result is approximate distance field.
stretch _s shape =
    let s = if (is_num _s) [_s,1,1]
            else if (is_vec2 _s) [..._s, 1]
            else do assert(is_vec3(_s) && shape.is_3d) in _s;
    in make_shape {
        // TODO: dist p : shape.dist[...p.[[X,Y,Z]] / s, p.[T]] * min(s),
        dist[x,y,z,t] : shape.dist[x/s.[X], y/s.[Y], z/s.[Z], t] * min(s),
        colour[x,y,z,t] : shape.colour[x/s.[X], y/s.[Y], z/s.[Z], t],
        bbox : [s*shape.bbox.[MIN], s*shape.bbox.[MAX]],
        is_2d : shape.is_2d,
        is_3d : shape.is_3d,
    };

warp_domain_xy warp shape =
    let bv = map (warp.function) [
            [shape.bbox.[MIN,X], shape.bbox.[MIN,Y]],
            [shape.bbox.[MAX,X], shape.bbox.[MIN,Y]],
            [shape.bbox.[MAX,X], shape.bbox.[MAX,Y]],
            [shape.bbox.[MIN,X], shape.bbox.[MAX,Y]],
        ];
    in make_shape {
        dist[x,y,z,t] =
            let p = warp.inverse[x,y];
            in warp.fix_distance[x,y,shape.dist[p.[X], p.[Y], z, t]];
        colour[x,y,z,t] =
            let p = warp.inverse[x,y];
            in shape.colour[p.[X], p.[Y], z, t];
        bbox = [
            [ ...min bv, shape.bbox.[MIN,Z] ],
            [ ...max bv, shape.bbox.[MAX,Z] ],
        ];
        is_2d = shape.is_2d;
        is_3d = shape.is_3d;
    };

// rotate angle shape -- 2D rotation around the Z axis
// rotate {angle=,axis=} shape -- 3D rotation around specified axis
rotate = match [
    a :: is_num -> shape ->
        let // 2D rotation around the Z axis
            rot2[a,p] = cmul[p, cis a];
            b = if (shape.bbox.[MIN,X] == -inf || shape.bbox.[MIN,Y] == -inf
                    || shape.bbox.[MAX,X] == inf || shape.bbox.[MAX,Y] == inf)
                [[-inf,-inf], [inf,inf]]
            else
                let bv = map (p->rot2[a,p]) [
                        [shape.bbox.[MIN,X], shape.bbox.[MIN,Y]],
                        [shape.bbox.[MAX,X], shape.bbox.[MIN,Y]],
                        [shape.bbox.[MAX,X], shape.bbox.[MAX,Y]],
                        [shape.bbox.[MIN,X], shape.bbox.[MAX,Y]],
                    ];
                in [min bv, max bv];
        in make_shape {
            dist[x,y,z,t] :
                let xy = rot2[-a,[x,y]];
                in shape.dist[xy.[X],xy.[Y],z,t],
            colour[x,y,z,t] :
                let xy = rot2[-a,[x,y]];
                in shape.colour[xy.[X],xy.[Y],z,t],
            bbox : [
                [ b.[MIN,X], b.[MIN,Y], shape.bbox.[MIN,Z] ],
                [ b.[MAX,X], b.[MAX,Y], shape.bbox.[MAX,Z] ],
            ],
            is_2d : shape.is_2d,
            is_3d : shape.is_3d,
        };
    {angle, axis: a} -> shape ->
        do
            assert(shape.is_3d);
        in let
            axis = normalize a;
            rot3[a,axis,p] =
                    p * cos(a)
                    - cross[axis, p] * sin(a)
                    + axis * dot[axis, p * (1 - cos a)];
            b = if (shape.bbox.[MIN,X] == -inf || shape.bbox.[MIN,Y] == -inf
                    || shape.bbox.[MIN,Z] == -inf
                    || shape.bbox.[MAX,X] == inf || shape.bbox.[MAX,Y] == inf
                    || shape.bbox.[MAX,Z] == inf
                )
                    shape.bbox
                else
                    let bb = shape.bbox;
                        bv = map (p->rot3[-angle,axis,p]) [
                            [bb.[MIN,X], bb.[MIN,Y], bb.[MIN,Z]],
                            [bb.[MAX,X], bb.[MIN,Y], bb.[MIN,Z]],
                            [bb.[MAX,X], bb.[MAX,Y], bb.[MIN,Z]],
                            [bb.[MIN,X], bb.[MAX,Y], bb.[MIN,Z]],
                            [bb.[MIN,X], bb.[MIN,Y], bb.[MAX,Z]],
                            [bb.[MAX,X], bb.[MIN,Y], bb.[MAX,Z]],
                            [bb.[MAX,X], bb.[MAX,Y], bb.[MAX,Z]],
                            [bb.[MIN,X], bb.[MAX,Y], bb.[MAX,Z]],
                        ];
                    in [ min bv, max bv ];
        in make_shape {
            //TODO: dist p : shape.dist[...rot3[angle,axis,p.[[X,Y,Z]]],p.[T]],
            dist[x,y,z,t] :
                let q = rot3[angle, axis, [x,y,z]];
                in shape.dist[q.[X],q.[Y],q.[Z],t],
            colour[x,y,z,t] :
                let q = rot3[angle, axis, [x,y,z]];
                in shape.colour[q.[X],q.[Y],q.[Z],t],
            bbox : b,
            is_3d : true,
        };
    ];

reflect_x shape =
    make_shape {
        dist[x,y,z,t] = shape.dist[-x,y,z,t];
        colour[x,y,z,t] = shape.colour[-x,y,z,t];
        bbox = [
            [-shape.bbox.[MAX,X], shape.bbox.[MIN,Y], shape.bbox.[MIN,Z]],
            [-shape.bbox.[MIN,X], shape.bbox.[MAX,Y], shape.bbox.[MAX,Z]],
        ];
        is_2d = shape.is_2d;
        is_3d = shape.is_3d;
    };
reflect_y shape =
    make_shape {
        dist[x,y,z,t] = shape.dist[x,-y,z,t];
        colour[x,y,z,t] = shape.colour[x,-y,z,t];
        bbox = [
            [shape.bbox.[MIN,X], -shape.bbox.[MAX,Y], shape.bbox.[MIN,Z]],
            [shape.bbox.[MAX,X], -shape.bbox.[MIN,Y], shape.bbox.[MAX,Z]],
        ];
        is_2d = shape.is_2d;
        is_3d = shape.is_3d;
    };
reflect_z shape =
    do assert(is_shape shape && shape.is_3d)
    in make_shape {
        dist[x,y,z,t] = shape.dist[x,y,-z,t];
        colour[x,y,z,t] = shape.colour[x,y,-z,t];
        bbox = [
            [shape.bbox.[MIN,X], shape.bbox.[MIN,Y], -shape.bbox.[MAX,Z]],
            [shape.bbox.[MAX,X], shape.bbox.[MAX,Y], -shape.bbox.[MIN,Z]],
        ];
        is_2d = shape.is_2d;
        is_3d = shape.is_3d;
    };
reflect_xy shape = // reflect shape across plane X=Y
    make_shape {
        dist[x,y,z,t] = shape.dist[y,x,z,t];
        colour[x,y,z,t] = shape.colour[y,x,z,t];
        bbox = [
            [shape.bbox.[MIN,Y], shape.bbox.[MIN,X], shape.bbox.[MIN,Z]],
            [shape.bbox.[MAX,Y], shape.bbox.[MAX,X], shape.bbox.[MAX,Z]],
        ];
        is_2d = shape.is_2d;
        is_3d = shape.is_3d;
    };
reflect_xz shape = // reflect shape across plane X=Z
    do assert(is_shape shape && shape.is_3d)
    in make_shape {
        dist[x,y,z,t] = shape.dist[z,y,x,t];
        colour[x,y,z,t] = shape.colour[z,y,x,t];
        bbox = [
            [shape.bbox.[MIN,Z], shape.bbox.[MIN,Y], shape.bbox.[MIN,X]],
            [shape.bbox.[MAX,Z], shape.bbox.[MAX,Y], shape.bbox.[MAX,X]],
        ];
        is_2d = shape.is_2d;
        is_3d = shape.is_3d;
    };
reflect_yz shape = // reflect shape across plane Y=Z
    do assert(is_shape shape && shape.is_3d)
    in make_shape {
        dist[x,y,z,t] = shape.dist[x,z,y,t];
        colour[x,y,z,t] = shape.colour[x,z,y,t];
        bbox = [
            [shape.bbox.[MIN,X], shape.bbox.[MIN,Z], shape.bbox.[MIN,Y]],
            [shape.bbox.[MAX,X], shape.bbox.[MAX,Z], shape.bbox.[MAX,Y]],
        ];
        is_2d = shape.is_2d;
        is_3d = shape.is_3d;
    };

reflect normal shape =
    if (normal == X_axis) reflect_x shape
    else if (normal == Y_axis) reflect_y shape
    else if (normal == Z_axis) reflect_z shape
    else error "reflect: normal vector $(normal) is not supported";

shear_x kx shape =
    warp_domain_xy {
        function[x,y] = [x + kx*y, y];
        inverse[x,y] = [x - kx*y, y];
        fix_distance[x,y,d] = d; // TODO
    } shape;

// TODO: In `taper`, the Lipschitz compensation code seems like overkill.
local_taper_x {range:[y0,y1], scale:[kx0, kx1]} shape =
    do assert(kx0 > 0 && kx1 > 0);
       assert(y0 < y1);
    in make_shape {
        dist[x,y,z,t] =
            (if (y <= y0)
                shape.dist[x/kx0,y,z,t] * kx0
            else if (y >= y1)
                shape.dist[x/kx1,y,z,t] * kx1
            else
                let kx = lerp[kx0,kx1,(y-y0)/(y1-y0)];
                in shape.dist[x/kx,y,z,t] * kx
            ) * min[kx0,kx1];
        colour[x,y,z,t] =
            if (y <= y0)
                shape.colour[x/kx0,y,z,t]
            else if (y >= y1)
                shape.colour[x/kx1,y,z,t]
            else
                let kx = lerp[kx0,kx1,(y-y0)/(y1-y0)];
                in shape.colour[x/kx,y,z,t];
        bbox = shape.bbox * [[max[kx0,kx1],1,1],[max[kx0,kx1],1,1]];
        is_2d = shape.is_2d;
        is_3d = shape.is_3d;
    };
local_taper_xy {range:[z0,z1], scale:[[kx0, ky0], [kx1, ky1]]} shape =
    do assert(kx0 > 0 && ky0 > 0 && kx1 > 0 && ky1 > 0);
       assert(z0 < z1);
    in make_shape {
        dist[x,y,z,t] =
            (if (z <= z0)
                shape.dist[x/kx0,y/ky0,z,t] * min[kx0,ky0]
            else if (z >= z1)
                shape.dist[x/kx1,y/ky1,z,t] * min[kx1,ky1]
            else
                let kx = lerp[kx0,kx1,(z-z0)/(z1-z0)];
                    ky = lerp[ky0,ky1,(z-z0)/(z1-z0)];
                in shape.dist[x/kx,y/ky,z,t] * min[kx,ky]
            ) * min[kx0,kx1,ky0,ky1];
        colour[x,y,z,t] =
            if (z <= z0)
                shape.colour[x/kx0,y/ky0,z,t]
            else if (z >= z1)
                shape.colour[x/kx1,y/ky1,z,t]
            else
                let kx = lerp[kx0,kx1,(z-z0)/(z1-z0)];
                    ky = lerp[ky0,ky1,(z-z0)/(z1-z0)];
                in shape.colour[x/kx,y/ky,z,t];
        bbox =
            let kx = max[kx0,kx1];
                ky = max[ky0,ky1];
            in shape.bbox * [[kx,ky,1],[kx,ky,1]];
        is_2d = shape.is_2d;
        is_3d = shape.is_3d;
    };

// Classic 'swirl' image transformation.
// From: Python skimage.transform.swirl
// http://scikit-image.org/docs/dev/auto_examples/transform/plot_swirl.html
swirl{strength, d} shape =
    let r = log 2 * d/10;
        f[x,y,z,t] =
            let m = mag[x,y];
                phi = strength * e^(-m/r) + phase[x,y];
                xy = m * cis phi;
            in [xy.[X],xy.[Y],z,t];
    in make_shape {
        dist p = shape.dist(f p);
        colour p = shape.colour(f p);
        bbox = shape.bbox;
        is_2d = shape.is_2d;
        is_3d = shape.is_3d;
    };

at p t shape = shape >> move(-p) >> t >> move p;

// The slice operators convert a 3D shape to a 2D shape
// by extracting a planar cross section.
slice = slice_xy;
slice_xy shape =
    make_shape {
        dist p : shape.dist[p.[X],p.[Y],0,p.[T]],
        colour p : shape.colour[p.[X],p.[Y],0,p.[T]],
        bbox : [
            [shape.bbox.[MIN,X], shape.bbox.[MIN,Y], 0],
            [shape.bbox.[MAX,X], shape.bbox.[MAX,Y], 0],
        ],
        is_2d : true,
    };
slice_xz shape =
    make_shape {
        dist p : shape.dist[p.[X],0,p.[Y],p.[T]],
        colour p : shape.colour[p.[X],0,p.[Y],p.[T]],
        bbox : [
            [shape.bbox.[MIN,X], shape.bbox.[MIN,Z], 0],
            [shape.bbox.[MAX,X], shape.bbox.[MAX,Z], 0],
        ],
        is_2d : true,
    };
slice_yz shape =
    make_shape {
        dist p : shape.dist[0,p.[X],p.[Y],p.[T]],
        colour p : shape.colour[0,p.[X],p.[Y],p.[T]],
        bbox : [
            [shape.bbox.[MIN,Y], shape.bbox.[MIN,Z], 0],
            [shape.bbox.[MAX,Y], shape.bbox.[MAX,Z], 0],
        ],
        is_2d : true,
    };

// 'extrude' converts a 2D shape to a 3D shape, linearly extruding it along the
// positive and negative Z axis by a distance of 'h'. The result is centered
// on the XY plane with a height of 2*h. Similar to Autocad extrude and
// OpenSCAD linear_extrude.
//
// The .mitred variant is cheaper, and has a mitred distance field.
// The .exact variant is the same shape with a Euclidean distance field.
extrude = {
    call = exact;
    mitred d shape =
        let h = d / 2;
        in make_shape {
            dist p : max[abs(p.[Z]) - h, shape.dist[p.[X],p.[Y],0,p.[T]]],
            colour p : shape.colour[p.[X],p.[Y],0,p.[T]],
            bbox : [
                [shape.bbox.[MIN,X], shape.bbox.[MIN,Y], -h],
                [shape.bbox.[MAX,X], shape.bbox.[MAX,Y], +h],
            ],
            is_3d : true,
        };
    exact d shape =
        let h = d / 2;
        in make_shape {
            dist p :
                let dz = abs(p.[Z]) - h;
                    dxy = shape.dist[p.[X],p.[Y],0,p.[T]];
                in mag(max[[dz,dxy],0]) + min[max[dz,dxy], 0],
            colour p : shape.colour[p.[X],p.[Y],0,p.[T]],
            bbox : [
                [shape.bbox.[MIN,X], shape.bbox.[MIN,Y], -h],
                [shape.bbox.[MAX,X], shape.bbox.[MAX,Y], +h],
            ],
            is_3d : true,
        };
};

morph t [shape1, shape2] =
    make_shape {
        dist p : lerp[shape1.dist p, shape2.dist p, t],
        colour p : lerp[shape1.colour p, shape2.colour p, t],
        bbox : lerp[shape1.bbox, shape2.bbox, t],
        is_2d : shape1.is_2d && shape2.is_2d,
        is_3d : shape1.is_3d && shape2.is_3d,
    };

loft d [shape1, shape2] =
    let r = d / 2;
    in make_shape {
        dist[x,y,z,t] =
            max[abs z - r,
                lerp[shape1.dist[x,y,0,t], shape2.dist[x,y,0,t],
                    (clamp[z, -r, r] + r) / d]];
        is_3d = true;
    };

////////////////
// REPETITION //
////////////////

repeat_x d shape =
    let r = d / 2;
    in make_shape {
        dist[x,y,z,t] : shape.dist[mod[x + r, d] - r, y, z, t],
        colour[x,y,z,t] : shape.colour[mod[x + r, d] - r, y, z, t],
        bbox : [
            [-inf, shape.bbox.[MIN,Y], shape.bbox.[MIN,Z]],
            [+inf, shape.bbox.[MAX,Y], shape.bbox.[MAX,Z]],
        ],
        is_2d : shape.is_2d,
        is_3d : shape.is_3d,
    };

repeat_xy d shape =
    let r = d / 2;
    in make_shape {
        dist[x,y,z,t] =
            shape.dist[
                mod[x + r.[X], 2*r.[X]] - r.[X],
                mod[y + r.[Y], 2*r.[Y]] - r.[Y],
                z, t];
        colour[x,y,z,t] =
            shape.colour[
                mod[x + r.[X], 2*r.[X]] - r.[X],
                mod[y + r.[Y], 2*r.[Y]] - r.[Y],
                z, t];
        bbox = [
            [-inf, -inf, shape.bbox.[MIN,Z]],
            [+inf, +inf, shape.bbox.[MAX,Z]],
        ];
        is_2d = shape.is_2d;
        is_3d = shape.is_3d;
    };

repeat_xyz d shape =
    let r = d / 2;
    in make_shape {
        dist [x,y,z,t] =
            shape.dist[
                mod[x + r.[X], 2*r.[X]] - r.[X],
                mod[y + r.[Y], 2*r.[Y]] - r.[Y],
                mod[z + r.[Z], 2*r.[Z]] - r.[Z],
                t];
        colour [x,y,z,t] =
            shape.colour[
                mod[x + r.[X], 2*r.[X]] - r.[X],
                mod[y + r.[Y], 2*r.[Y]] - r.[Y],
                mod[z + r.[Z], 2*r.[Z]] - r.[Z],
                t];
        bbox = [
            [-inf, -inf, -inf],
            [+inf, +inf, +inf],
        ];
        is_2d = shape.is_2d;
        is_3d = shape.is_3d;
    };

repeat_finite d l shape =
    make_shape {
        dist [x,y,z,t] =
            shape.dist[
                x - d.[X]*clamp[round(x/d.[X]), 0, l.[X] - 1],
                y - d.[Y]*clamp[round(y/d.[Y]), 0, l.[Y] - 1],
                z - d.[Z]*clamp[round(z/d.[Z]), 0, l.[Z] - 1],
                t];
        colour [x,y,z,t] =
            shape.colour[
                x - d.[X]*clamp[round(x/d.[X]), 0, l.[X] - 1],
                y - d.[Y]*clamp[round(y/d.[Y]), 0, l.[Y] - 1],
                z - d.[Z]*clamp[round(z/d.[Z]), 0, l.[Z] - 1],
                t];
        bbox = [
            shape.bbox.[MIN],
            shape.bbox.[MIN] + d*l,
        ];
        is_2d = shape.is_2d;
        is_3d = shape.is_3d;
    };

repeat_mirror_x shape =
    make_shape {
        dist[x,y,z,t] : shape.dist[abs x, y, z, t],
        colour[x,y,z,t] : shape.colour[abs x, y, z, t],
        bbox : [
            [-shape.bbox.[MAX,X], shape.bbox.[MIN,Y], shape.bbox.[MIN,Z]],
            shape.bbox.[MAX],
        ],
        is_2d : shape.is_2d,
        is_3d : shape.is_3d,
    };

repeat_radial reps shape =
    let angle = tau/reps;
        inradius = -shape.bbox.[MIN,Y];
        circumradius = inradius / cos(angle/2);
        max_x =
            if (mod[reps,4]==0)
                inradius
            else if (mod[reps,2]==0)
                circumradius
            else
                circumradius * sin((floor(reps/4)+.5)*angle);
        ashift = tau/4 + angle/2;
        f[x,y,z,t] =
            let a = phase[x,y] + ashift;
                r = mag[x,y];
                a2 = mod[a,angle] - ashift;
                xy = cis(a2) * r;
            in [xy.[X], xy.[Y], z, t];
    in make_shape {
        dist p = shape.dist(f p);
        colour p = shape.colour(f p);
        // bbox is exact for regular_polygon, approximate for general case.
        bbox = [
            [ -max_x,
              -inradius,
              shape.bbox.[MIN,Z] ],
            [ max_x,
              if (mod[reps,2]==0) inradius else circumradius,
              shape.bbox.[MAX,Z] ]
        ];
        is_2d = shape.is_2d;
        is_3d = shape.is_3d;
    };

////////////
// COLOUR //
////////////

colour c shape = c >> match [
    c :: is_vec3 ->
        make_shape {
            ... shape,
            colour p : c,
        };
    [ifield, cmap] ->
        make_shape {
            ... shape,
            colour : compose[ifield, cmap],
        };
    c :: is_shape ->
        make_shape {
            ... shape,
            colour : c.colour,
        };
    c :: is_func ->
        make_shape {
            ... shape,
            colour : c,
        };
];

show_colour c = colour c everything;

// These conversions are not currently used
// Correct gamma encoding/compression, more accurate but slower than c^(1/2.2)
linear2sRGB c = if (c <= 0.0031308) 12.92*c else ( 1.055*c^(1/2.4) - 0.055 );
// Correct gamma decoding/expansion, more accurate but slower than c^2.2
sRGB2linear c = if (c <= 0.04045) c/12.92 else ((c+0.055)/1.055)^2.4;

sRGB = {
    call c = c^2.2;  // Faster approximate gamma decoding to convert to linearRGB
    hue h = HSV[h,1,1];
    grey i = call[i,i,i];
    // HSV to RGB adapted from https://www.shadertoy.com/view/MsS3Wc
    // which is Copyright 2014 Inigo Quilez (The MIT License)
    HSV [h,s,v] =
        let m = abs(mod[h + [1,2/3,1/3], 1] * 6 - 3);
        in call(v * lerp[1, clamp[m - 1, 0, 1], s]);
};

// OKLab to linear RGB
// https://bottosson.github.io/posts/oklab
// Uses updated matrices published on 2021-01-25
OKLab [L,a,b] =
    let toLMS = [
        [1.0         ,  1.0         ,  1.0         ],  // *L
        [0.3963377774, -0.1055613458, -0.0894841775],  // *a
        [0.2158037573, -0.0638541728, -1.2914855480]]; // *b
    LMStoCONE = [
        [4.0767416621, -1.2684380046, -0.0041960863],  // *l
       [-3.3077115913,  2.6097574011, -0.7034186147],  // *m
        [0.2309699292, -0.3413193965,  1.7076147010]]; // *s
    lms = dot[[L,a,b], toLMS];
    in clamp[dot[lms^3, LMStoCONE], 0, 1];
OKLCH [L,C,H] = OKLab[L, C*cos(H*tau), C*sin(H*tau)];
OKhue h = OKLCH [0.70,0.30,h];

gradient_cmap =
    {
        call = OKLab;
        linearRGB [col1, col2, i] = lerp[col1, col2, i];
        // Optimized linear-RGB color mix in OKLab space
        // Adapted from https://www.shadertoy.com/view/ttcyRS
        // Copyright © 2020 Inigo Quilez
        // The MIT License
        OKLab [col1, col2, i] =
            // https://bottosson.github.io/posts/oklab
            // Uses updated matrices published on 2021-01-25
            let CONEtoLMS = [
                [0.4122214708,  0.2119034982,  0.0883024619],  // *r = l
                [0.5363325363,  0.6806995451,  0.2817188376],  // *g = m
                [0.0514459929,  0.1073969566,  0.6299787005]]; // *b = s
            LMStoCONE = [
                [4.0767416621, -1.2684380046, -0.0041960863],  // *l
               [-3.3077115913,  2.6097574011, -0.7034186147],  // *m
                [0.2309699292, -0.3413193965,  1.7076147010]]; // *s
            // RGB to cone
            lmsA = dot[col1, CONEtoLMS]^(1/3);
            lmsB = dot[col2, CONEtoLMS]^(1/3);
            lms = lerp[lmsA, lmsB, i];
            in
            dot[lms^3, LMStoCONE];
    };

webRGB[r,g,b] = sRGB[r/255, g/255, b/255];

red = [1, 0, 0];
orange = sRGB[1, .5, 0];
yellow = [1, 1, 0];
chartreuse = sRGB.hue(1/4);
green = [0, 1, 0];
spring_green = sRGB.hue(5/12);
cyan = [0, 1, 1];
azure = sRGB.hue(7/12);
blue = [0, 0, 1];
indigo = sRGB.hue(3/4);
magenta = [1, 0, 1];
rose = sRGB.hue(11/12);

white = [1, 1, 1];
grey = sRGB[.5,.5,.5];
black = [0, 0, 0];

texture = colour;
make_texture x = colour x everything;

i_linear d [x,_,_,_] = x/d `mod` 1;
i_radial n [x,y,_,_] = (phase[x,y]/tau-.25)*n `mod` 1;
i_concentric d [x,y,_,_] = mag[x,y]/d `mod` 1;
i_gyroid [x,y,z,t] = (cos(x)*sin(y) + cos(y)*sin(z) + cos(z)*sin(x) + 1.5)/3;
i_animate period ifield [x,y,z,t] = ifield[x,y,z,t] + t/period `mod` 1;
show_ifield ifield = make_texture[ifield, sRGB.grey];

show_cmap f = union[
    rect[4.05,1.05] >> colour grey,
    rect[4,1] >> colour ([x,_,_,_]->f((x+2)/4))];

///////////////////////////////////////////////////////////////
// RESEARCH -- shapes and operators with bad distance fields //
///////////////////////////////////////////////////////////////

// This is the infinite connected 3D volume enclosed by the gyroid surface.
// Use `shell 0` to obtain the surface itself.
// https://en.wikipedia.org/wiki/Gyroid
// Lipschitz constant == 4/3. Apply 'lipschitz 1.33' after extracting isosurfaces
gyroid =
    make_shape {
        dist[x,y,z,t] = cos(x)*sin(y) + cos(y)*sin(z) + cos(z)*sin(x);
        is_3d = true;
    };

// Twist 3D space into a helix around the Z axis.
// ``tr`` is the twist rate: angle of twist (counterclockwise)
// per unit increase along the Z axis.
// Bad distance field.
twist tr shape =
    let b = shape.bbox;
        r = max(map mag [
                [b.[MIN,X], b.[MIN,Y]],
                [b.[MAX,X], b.[MIN,Y]],
                [b.[MIN,X], b.[MAX,Y]],
                [b.[MAX,X], b.[MAX,Y]],
            ]);
    in make_shape {
        dist p :
            let xy = cmul[p.[[X,Y]], cis(p.[Z]*-tr)];
            in shape.dist[xy.[X], xy.[Y], p.[Z], p.[T]],
        colour p :
            let xy = cmul[p.[[X,Y]], cis(p.[Z]*-tr)];
            in shape.colour[xy.[X], xy.[Y], p.[Z], p.[T]],
        bbox : [[-r,-r,b.[MIN,Z]],[r,r,b.[MAX,Z]]],
        is_3d : true,
    };

// Bend a 2D/3D shape upward/towards +Y, through the angle `a`,
// around a circle/cylinder of diameter `d` centred on the origin/Z axis.
// The angle defaults to `tau`/`360*deg`.
// TODO: bbox is larger than necessary.
bend {d = null, angle: a = tau} shape =
    let width = shape.bbox.[MAX,X] - shape.bbox.[MIN,X];
        rx = width / 2 / (a/tau);
        ry = if (d==null) width/tau else d/2;
        offset = [
            (shape.bbox.[MAX,X] + shape.bbox.[MIN,X]) / 2,
            shape.bbox.[MAX,Y] + ry,
        ];
        f[x,y,z,t] =
            let nx = (mod[phase[x,y]/pi+1.5,2]-1)*rx+offset.[X];
                ny = offset.[Y] - mag[x,y];
            in [nx,ny,z,t];
        ymax = shape.bbox.[MAX,Y] - shape.bbox.[MIN,Y] + ry;
    in make_shape {
        dist p = shape.dist(f p);
        colour p = shape.colour(f p);
        bbox = [
            [-ymax, -ymax, shape.bbox.[MIN,Z]],
            [+ymax, +ymax, shape.bbox.[MAX,Z]],
        ];
        is_2d = shape.is_2d;
        is_3d = shape.is_3d;
    };

///////////////////////////
// LOW LEVEL & DEBUGGING //
///////////////////////////
distance_field shape =
    make_shape {
        dist : everything.dist,
        colour p :
            let d = shape.dist p;
                c = if (d >= 0) ceil(d)-d else floor(d)-d;
            in [max[-c,0], bit(abs(d)==inf), max[c,0]],
        is_2d : true,
    };

// Select an isosurface of the shape's distance field.
// r > 0: inflate the shape, blow it up like a balloon.
// r == 0: no effect.
// r < 0: deflate the shape.
//
// If the distance field is Euclidean, then positive r yields
// the Minkowski sum of a sphere of radius r, aka rounded offset (CAD),
// aka dilation (Mathematical Morphology).
//
// If the distance field is mitred, the result is a mitred offset (CAD).
//
// WARNING: Two shapes can have identical boundaries, but different distance
// fields, so this is sensitive to the implementation of the shape.
// Useful for debugging, to visualize a distance field.
offset r shape =
    make_shape {
        dist p : shape.dist p - r,
        colour p : shape.colour p,
        bbox : [shape.bbox.[MIN]-r, shape.bbox.[MAX]+r],
        is_2d : shape.is_2d,
        is_3d : shape.is_3d,
    };

// Construct a shell of thickness d,
// whose boundaries are +/- d/2 from the original surface or perimeter.
// WARNING: this is implemented using isosurfaces, so the results are dependent
// on the distance field, not just on the boundary of the shape.
shell d shape =
    let r = d / 2;
    in make_shape {
        dist p : abs(shape.dist p) - r,
        colour : shape.colour,
        bbox : [shape.bbox.[MIN] - r, shape.bbox.[MAX] + r],
        is_2d : shape.is_2d,
        is_3d : shape.is_3d,
    };

pancake d shape = shape >> extrude.exact 0 >> offset (d/2);

// Super kludgy function to scale a shape's distance field,
// to compensate for a Lipschitz constant that is > 1.
// So 'k' is the Lipschitz constant.
// Try k=2 to fix a 3D shape that doesn't render properly.
lipschitz k shape =
    make_shape {
        dist p : shape.dist p / k,
        colour : shape.colour,
        bbox : shape.bbox,
        is_2d : shape.is_2d,
        is_3d : shape.is_3d,
    };

// Colour the XY plane with black where the distance field gradient <= 1 (okay)
// or with rainbow colour where it is > 1 and <= n (bad)
// or with white where it is > n (worst).
// Debug strategy: find the smallest value of n that doesn't produce white,
// and that's your Lipschitz constant. Try n=2 as an initial value.
show_gradient [j,k] shape =
    let xy = make_shape {
        dist[x,y,z,t] : abs(z),
        colour p :
            let eps = 0.01;
                d = shape.dist[p.[X],p.[Y],0,p.[T]];
                up = shape.dist[p.[X],p.[Y]+eps,0,p.[T]];
                down = shape.dist[p.[X],p.[Y]-eps,0,p.[T]];
                left = shape.dist[p.[X]-eps,p.[Y],0,p.[T]];
                right = shape.dist[p.[X]+eps,p.[Y],0,p.[T]];
                g = max[abs(d-up),abs(d-down),abs(d-left),abs(d-right)] / eps;
            in if (g <= j)
                [0,0,0]
            else if (g > k)
                [1,1,1]
            else
                sRGB.hue((g-j)/(k-j)),
        is_2d : true,
        is_3d : true,
    }
    in (if (shape.is_2d) xy else union[shape, xy])
        >> set_bbox (shape.bbox);

show_dist shape =
    let xy = make_shape {
        dist[x,y,z,t] : abs(z),
        colour p :
            let eps = 0.01;
                d = shape.dist p;
                up = shape.dist[p.[X],p.[Y]+eps,0,p.[T]];
                down = shape.dist[p.[X],p.[Y]-eps,0,p.[T]];
                left = shape.dist[p.[X]-eps,p.[Y],0,p.[T]];
                right = shape.dist[p.[X]+eps,p.[Y],0,p.[T]];
                g = max[abs(d-up),abs(d-down),abs(d-left),abs(d-right)] / eps;
                c = if (d >= 0) ceil(d)-d else floor(d)-d;
            in if (d == d)
                if (d == inf)
                    [0,1,1]
                else if (d == -inf)
                    sRGB[0,.5,.5]
                else
                    sRGB[clamp[g-1,0,1], max[-c,0], max[c,0]]
            else
                [1,1,1],
        bbox : [[-inf,-inf,0],[inf,inf,0]],
        is_2d : true,
        is_3d : shape.is_3d,
    };
    in (if (shape.is_2d) xy else union[shape, xy])
        >> set_bbox (shape.bbox);

show_axes shape =
    let eps = .05;
    in
        (if (shape.is_2d)
            union [
                shape,
                rect[eps, inf] >> colour green,
                rect[inf, eps] >> colour red,
            ]
        else
            union [
                shape,
                box[inf, eps, eps] >> colour red,
                box[eps, inf, eps] >> colour green,
                box[eps, eps, inf] >> colour blue,
            ])
        >> set_bbox (shape.bbox);

set_bbox (bbox :: is_bbox3) shape = {...shape; bbox: bbox};

show_bbox shape =
    let bb =
        if (shape.is_2d)
            rect[shape.bbox.[MIN,[X,Y]], shape.bbox.[MAX,[X,Y]]]
        else
            box[shape.bbox.[MIN], shape.bbox.[MAX]];
    in union[
        bb >> colour (sRGB.grey .2),
        shape `difference` bb >> colour red,
        shape >> colour green `intersection` bb,
    ];

/// VALUE PICKERS ///
slider[lo,hi] =
  { call = is_num;
    picker = {slider:[lo,hi]};
  };
int_slider[lo,hi] =
  { call = is_num;
    picker = {int_slider:[lo,hi]};
  };
scale_picker =
  { call = is_num;
    picker = #scale_picker;
  };
checkbox =
  { call = is_bool;
    picker = #checkbox;
  };
colour_picker =
  { call = is_vec3;
    picker = #colour_picker;
  };

}
