http://mathling.com/sdf  library module

Imports

import module namespace sd2="http://mathling.com/sdf/2d"
       at "../sdf/sd2.xqy"

import module namespace sd3="http://mathling.com/sdf/3d"
       at "../sdf/sd3.xqy"

import module namespace cpoly="http://mathling.com/geometric/complex-polygon"
       at "../geo/complex-polygon.xqy"

import module namespace dist="http://mathling.com/type/distribution"
       at "../types/distributions.xqy"

import module namespace wrapper="http://mathling.com/type/wrapper"
       at "../types/wrapper.xqy"

import module namespace box="http://mathling.com/geometric/rectangle"
       at "../geo/rectangle.xqy"

import module namespace slot="http://mathling.com/type/slot"
       at "../types/slot.xqy"

import module namespace affine="http://mathling.com/geometric/affine"
       at "../geo/affine.xqy"

import module namespace rgb="http://mathling.com/colour/rgb"
       at "../colourspace/rgb.xqy"

import module namespace point="http://mathling.com/geometric/point"
       at "../geo/point.xqy"

import module namespace rand="http://mathling.com/core/random"
       at "../core/random.xqy"

import module namespace path="http://mathling.com/geometric/path"
       at "../geo/path.xqy"

import module namespace edge="http://mathling.com/geometric/edge"
       at "../geo/edge.xqy"

import module namespace xyz="http://mathling.com/colour/xyz"
       at "../colourspace/xyz.xqy"

import module namespace ellipse="http://mathling.com/geometric/ellipse"
       at "../geo/ellipse.xqy"

import module namespace ann="http://mathling.com/sdf/annotations"
       at "../sdf/annotations.xqy"

import module namespace solid="http://mathling.com/geometric/solid"
       at "../geo/solid.xqy"

import module namespace f="http://mathling.com/core/callable"
       at "../core/callable.xqy"

import module namespace gradient="http://mathling.com/svg/gradients"
       at "../svg/gradients.xqy"

import module namespace cs="http://mathling.com/colour/space"
       at "../colourspace/colour-space.xqy"

import module namespace util="http://mathling.com/core/utilities"
       at "../core/utilities.xqy"

import module namespace light="http://mathling.com/type/light"
       at "../types/light.xqy"

import module namespace v="http://mathling.com/core/vector"
       at "../core/vector.xqy"

import module namespace config="http://mathling.com/core/config"
       at "../core/config.xqy"

import module namespace errors="http://mathling.com/core/errors"
       at "../core/errors.xqy"

Functions

declare function this:trace(
  $f as item()
) as map(*)
{
  if (not(ann:is-sdf($f))) then errors:error("ML-BADARGS", ("f", $f)) else (),
  let $fn := f:function($f)
  return
  ann:f(util:function-name($f), 
    function ($point as xs:double*) as xs:double* {
      $fn($point)=>trace(util:quote($f)||"("||util:quote($point)||")")
    }
  )
}

declare function this:as-point(
  $f as item()
) as map(*)
{
  let $fn := f:function($f)
  return
  ann:f(util:function-name($f), 
    function ($point as xs:double*) as xs:double* {
      $fn(point:vector($point))=>point:pcoordinates()
    }
  )
}

declare function this:as-vector(
  $f as item()
) as map(*)
{
  let $fn := f:function($f)
  return
  ann:f(util:function-name($f),
    function ($point as map(xs:string,item()*)) as map(xs:string,item()*) {
      $fn(point:pcoordinates($point))=>point:vector()
    }
  )
}

declare function this:opConstant(
  $k as xs:double
) as map(*)
{
  ann:f("sdf:opConstant", function ($point as xs:double*) as xs:double* {
    $k
  })
}

declare function this:opMaterial(
  $base as item(),
  $material as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opMaterial",
    function ($point as xs:double*) as xs:double* {
      ($basef($point), $material)
    },
    ($base, $material)
  )
}

declare function this:opTranslate(
  $base as item(),
  $loc as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $matrix := affine:translation3(point:px($loc), point:py($loc), point:pz($loc))
  let $basef := f:function($base)
  return
    ann:f("sdf:opTranslate", function ($point as xs:double*) as xs:double* {
        $basef(affine:affine3($point, $matrix))
      },
      ($base, $loc)
    )
}

declare function this:opRotate(
  $base as item(), 
  $rotation as xs:double (: degrees :)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $matrix := affine:rotation2($rotation)
  let $basef := f:function($base)
  return
    ann:f("sdf:opRotate",
      function ($point as xs:double*) as xs:double* {
        $basef(affine:affine2($point, $matrix))
      },
      ($base, $rotation)
    )
}

declare function this:opRotate(
  $base as item(),
  $roll as xs:double, (: all degrees :)
  $pitch as xs:double,
  $yaw as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $matrix := affine:rotation3($roll, $pitch, $yaw)
  let $basef := f:function($base)
  return
    ann:f("sdf:opRotate",
      function ($point as xs:double*) as xs:double* {
        $basef(affine:affine3($point=>v:as-dimension(3), $matrix))
      },
      ($base, $roll, $pitch, $yaw)
    )
}

declare function this:opScale(
  $base as item(), 
  $scale as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opScale",
    function ($point as xs:double*) as xs:double* {
      $basef($point=>v:div($scale))=>v:times($scale)
    },
    ($base, $scale)
  )
}

declare function this:opReflect(
  $base as item(),
  $center as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $matrix := affine:reflection3($center)
  let $basef := f:function($base)
  return
    ann:f("sdf:opReflect",
      function ($point as xs:double*) as xs:double* {
        $basef(affine:affine3($point, $matrix))
      },
      ($base, $center)
    )
}

declare function this:opReflectXZ(
  $base as item(),
  $start as map(xs:string,item()*),
  $end as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $matrix := affine:reflection2($start, $end)
  let $basef := f:function($base)
  return
    ann:f("sdf:opReflectXZ",
      function ($point as xs:double*) as xs:double* {
        let $qxz := affine:affine2((v:px($point), v:pz($point)), $matrix)
        return $basef((v:px($qxz), v:py($point), v:pz($qxz)))
      },
      ($base, $start, $end)
    )
}

declare function this:opReflectXY(
  $base as item(),
  $start as map(xs:string,item()*),
  $end as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $matrix := affine:reflection2($start, $end)
  let $basef := f:function($base)
  return
    ann:f("sdf:opReflectXY",
      function ($point as xs:double*) as xs:double* { 
        let $qxy := affine:affine2((v:px($point), v:py($point)), $matrix)
        return $basef((v:px($qxy), v:py($qxy), v:pz($point)))
      },
      ($base, $start, $end)
    )
}

declare function this:opReflectYZ(
  $base as item(),
  $start as map(xs:string,item()*),
  $end as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $matrix := affine:reflection2($start, $end)
  let $basef := f:function($base)
  return
    ann:f("sdf:opReflectYZ",
      function ($point as xs:double*) as xs:double* {
        let $qyz := affine:affine2((v:py($point), v:pz($point)), $matrix)
        return $basef((v:px($point), v:py($qyz), v:pz($qyz)))
      },
      ($base, $start, $end)
    )
}

declare function this:opRepeat(
  $base as item(),
  $periods as map(xs:string,item()*) 
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $c := point:pcoordinates($periods)
  let $half-c := $c=>v:times(0.5)
  let $basef := f:function($base)
  return
    ann:f("sdf:opRepeat",
      function ($point as xs:double*) as xs:double* {
        let $q := $point=>v:add($half-c)=>v:modulo($c)=>v:sub($half-c)
        return $basef($q)
      },
      ($base, $periods)
    )
}

declare function this:opLimitedRepeat(
  $base as item(),
  $periods as map(xs:string,item()*),
  $lower-limits as map(xs:string,item()*),
  $upper-limits as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $c := point:pcoordinates($periods)
  let $lim-a := point:pcoordinates($lower-limits)
  let $lim-b := point:pcoordinates($upper-limits)
  let $basef := f:function($base)
  return
    ann:f("sdf:opLimitedRepeat",
      function ($point as xs:double*) as xs:double* {
        $basef(
          $point=>v:sub(
            v:clampv(v:round($point=>v:divide($c)), $lim-a, $lim-b)=>v:multiply($c)
          )
        )
      },
      ($base, $periods, $lower-limits, $upper-limits)
    )
}

declare function this:mod($x as xs:double, $k as xs:double) as xs:double
{
  $x - $k * floor($x div $k)
}

declare function this:vmod($v as xs:double*, $k as xs:double) as xs:double*
{
  v:map(function ($c as xs:double) as xs:double {this:mod($c, $k)}, $v)
}

declare function this:opRepeat2(
  $base as item(),
  $size as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  let $s2 := $size div 2
  return
    ann:f("sdf:opRepeat2",
      function ($point as xs:double*) as xs:double* {
        $basef($point=>v:plus($s2)=>this:vmod($size)=>v:minus($s2))
      },
      ($base, $size)
    )
}

declare function this:opMirrorRepeat2(
  $base as item(),
  $size as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  let $s2 := $size div 2
  return
    ann:f("sdf:opMirrorRepeat2",
      function ($point as xs:double*) as xs:double* {
        let $p2 := $point=>v:plus($s2)
        return
        $basef(
         v:multiply(
           ($p2=>this:vmod($size)=>v:minus($s2)),
           ($p2=>v:div($size)=>v:floor()=>this:vmod(2)=>v:times(2)=>v:minus(1))
			)
        )
      },
      ($base, $size)
    )
}

declare function this:opGridRepeat2(
  $base as item(),
  $size as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  let $s2 := $size div 2
  let $mirrorf := (: opMirrorRepeat2 with no base function :)
    function ($point as xs:double*) as xs:double* {
      let $p2 := $point=>v:plus($s2)
      return
       v:multiply(
         ($p2=>this:vmod($size)=>v:minus($s2)),
         ($p2=>v:div($size)=>v:floor()=>this:vmod(2)=>v:times(2)=>v:minus(1))
       )
    }
  return
    ann:f("sdf:opGridRepeat2",
      function ($point as xs:double*) as xs:double* {
        let $rm := $mirrorf($point)
        let $p := $rm=>v:minus($s2)
        let $p := if (v:px($p) > v:py($p)) then (v:py($p), v:px($p)) else $p
        return (
          $basef($p)
        )
      },
      ($base, $size)
    )
}

declare function this:opPolarRepeat2(
  $base as item(),
  $n as xs:integer
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  let $θ := 2 * math:pi() div $n
  let $halfθ := $θ div 2
  return
    ann:f("sdf:opPolarRepeat2",
      function ($point as xs:double*) as xs:double* {
        let $α := this:mod(math:atan2(v:py($point), v:px($point)) + $halfθ, $θ) - $halfθ
        let $p := (math:cos($α), math:sin($α))=>v:times(v:magnitude($point))
        return $basef($p)
      },
      ($base, $n)
    )
}

declare function this:opElongate(
  $base as item() ,
  $h as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $hv := point:pcoordinates($h)
  let $basef := f:function($base)
  return
    ann:f("sdf:opElongate",
      function ($point as xs:double*) as xs:double* {
        let $q := $point=>v:abs()=>v:sub($hv)
        return $basef($q=>v:at-max(0)) + min((max($q),0))
      },
      ($base, $h)
    )
}

declare function this:opRound(
  $base as item(),
  $r as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opRound",
    function ($point as xs:double*) as xs:double* {
      $basef($point) - $r
    },
    ($base, $r)
  )
}

declare function this:opOnion(
  $base as item(),
  $thickness as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $r := $thickness div 2
  let $basef := f:function($base)
  return
    ann:f("sdf:opOnion",
      function ($point as xs:double*) as xs:double* {
        abs($basef($point)) - $r
      },
      ($base, $thickness)
    )
}

declare function this:opManyOnion(
  $base as item(),
  $thickness as xs:double,
  $n as xs:integer
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  fold-left(
    1 to $n,
    $base,
    function ($f as item(), $i as xs:integer) as map(*) {
      this:opOnion($f, $thickness*math:pow(2, -($i - 1)))
    }
  )
}

declare function this:opExtrude(
  $base as item(),  (: 2D :) 
  $h as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opExtrude",
    function ($point as xs:double*) as xs:double* {
      let $d := $basef($point=>v:as-dimension(2))
      let $w := ( $d, abs(v:pz($point)) - $h )
      return min((max($w), 0)) + v:magnitude($w=>v:at-max(0))
    },
    ($base, $h)
  )
}

declare function this:opRevolve(
  $base as item(), (: 2D :)
  $o as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opRevolve",
    function ($point as xs:double*) as xs:double* {
      $basef( ( v:magnitude((v:px($point), v:pz($point))) - $o, v:py($point) ) )
    },
    ($base, $o)
  )
}

declare function this:opUnion(
  $a as item(),
  $b as item()
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opUnion",
    function ($point as xs:double*) as xs:double* {
      let $v1 := $af($point)
      let $v2 := $bf($point)
      return
        if (head($v1) < head($v2)) then $v1 else $v2
    },
    ($a, $b)
  )
}

declare function this:opUnion(
  $fs as item()*
) as map(*)
{
  for $f in $fs return if (not(ann:is-sdf($f))) then errors:error("ML-BADARGS", ("fs", $f)) else (),
  let $ffs := for $f in $fs return f:function($f)
  return
  ann:f("sdf:opUnion",
    function ($point as xs:double*) as xs:double* {
      let $ix :=
        util:min-index($ffs,
          function ($f as function(xs:double*) as xs:double*) as xs:double {
            head($f($point))
          }
        )
      return (
        $ffs[$ix]($point)
      )
    },
    $fs
  )
}

declare function this:min(
  $a as xs:double,
  $b as xs:double
) as xs:double
{
  min(($a,$b))
}

declare function this:sminQuad(
  $a as xs:double,
  $b as xs:double,
  $k as xs:double (: 0.1 :)
) as xs:double
{
  let $h := max(($k - abs($a - $b), 0)) div $k
  return min(($a,$b)) - $h*$h*$k*(1.0 div 4.0)
}

declare function this:sminCubic(
  $a as xs:double,
  $b as xs:double,
  $k as xs:double (: 0.1 :)
) as xs:double
{
  let $h := max(($k - abs($a - $b), 0)) div $k
  return min(($a,$b)) - $h*$h*$h*$k*(1.0 div 6.0)
}

declare function this:opUnion(
  $a as item(),
  $b as item(),
  $min as function(xs:double,xs:double) as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opUnion",
    function ($point as xs:double*) as xs:double* {
      let $v1 := $af($point)
      let $v2 := $bf($point)
      return (
        $min(head($v1), head($v2)),
        if (head($v1) < head($v2)) then tail($v1) else tail($v2)
      )
    },
    ($a, $b, $min)
  )
}

declare function this:opSubtract(
  $a as item(),
  $b as item()
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opSubtract",
    function ($point as xs:double*) as xs:double* {
      let $v1 := $af($point)
      let $v2 := $bf($point)
      return (
        max((head($v1), -head($v2))),
        if (head($v1) > -head($v2)) then tail($v1) else tail($v2)
      )
    },
    ($a, $b)
  )
}

declare function this:opIntersect(
  $a as item(),
  $b as item()
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opIntersect",
    function ($point as xs:double*) as xs:double* {
      let $v1 := $af($point)
      let $v2 := $bf($point)
      return (
        if (head($v1) > head($v2)) then $v1 else $v2
      )
    },
    ($a, $b)
  )
}

declare function this:opIntersect(
  $fs as item()*
) as map(*)
{
  for $f in $fs return if (not(ann:is-sdf($f))) then errors:error("ML-BADARGS", ("fs", $f)) else (),
  let $ffs := for $f in $fs return f:function($f)
  return
  ann:f("sdf:opIntersect",
    function ($point as xs:double*) as xs:double* {
      let $ix :=
        util:max-index($ffs,
          function ($f as function(xs:double*) as xs:double*) as xs:double {
            head($f($point))
          }
        )
      return (
        $ffs[$ix]($point)
      )
    },
    ($fs)
  )
}

declare function this:opRoundUnion(
  $a as item(), 
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opRoundUnion",
    function ($point as xs:double*) as xs:double* {
      let $v2 := $af($point)
      let $v1 := $bf($point)
      let $d1 := head($v1)
      let $d2 := head($v2)
      let $h := max(($k - abs($d1 - $d2), 0))
      return (
        min(($d1, $d2)) - $h*$h*0.25 div $k,
        if (head($v2) < head($v1)) then tail($v2) else tail($v1)
      )
    },
    ($a, $b, $k)
  )
}

declare function this:opRoundSubtract(
  $a as item(),
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opRoundSubtract",
    function ($point as xs:double*) as xs:double* {
      let $l := v:magnitude($point)
      let $v2 := $af($point)
      let $v1 := $bf($point)
      let $d1 := head($v1)
      let $d2 := head($v2)
      let $h :=  max(($k - abs(-$d1 - $d2),0))
      return (
        max((-$d1, $d2)) + $h*$h*0.25 div $k,
        if ($d2 > -$d1) then tail($v2) else tail($v1)
      )
    },
    ($a, $b, $k)
  )
}

declare function this:opRoundIntersect(
  $a as item(), 
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opRoundIntersect",
    function ($point as xs:double*) as xs:double* {
      let $v2 := $af($point)
      let $v1 := $bf($point)
      let $d1 := head($v1)
      let $d2 := head($v2)
      let $h := max(($k - abs($d1 - $d2), 0))
      return (
        max(($d1, $d2)) + $h*$h*0.25 div $k,
        if (head($v2) > head($v1)) then tail($v2) else tail($v1)
      )
    },
    ($a, $b, $k)
  )
}

declare function this:opSmoothUnion(
  $a as item(), 
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opSmoothUnion",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      let $h := util:clamp(0.5 + (0.5 * ($bd - $ad)) div $k, 0, 1)
      return (
        util:mix($bd, $ad, $h) - $k * $h * (1 - $h),
        if (head($av) < head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $k)
  )
}

declare function this:opSmoothSubtract(
  $a as item(),
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opSmoothSubtract",
    function ($point as xs:double*) as xs:double* {
      let $l := v:magnitude($point)
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      let $h :=  util:clamp(0.5 - (0.5 * ($ad + $bd)) div $k, 0, 1)
      return (
        util:mix($ad, -$bd, $h) + $k * $h * (1 - $h),
        if ($ad > -$bd) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $k)
  )
}

declare function this:opSmoothIntersect(
  $a as item(), 
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opSmoothIntersect",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      let $h := util:clamp(0.5 - (0.5 * ($bd - $ad)) div $k, 0, 1)
      return (
        util:mix($bd, $ad, $h) + $k * $h * (1 - $h),
        if (head($av) > head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $k)
  )
}

declare function this:opChamferUnion(
  $a as item(), 
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  let $sqrt_2_2 := math:sqrt(2) div 2
  return
  ann:f("sdf:opChamferUnion",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      return (
        min(($ad, $bd, ($ad - $k + $bd) * $sqrt_2_2)),
        if (head($av) > head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $k)
  )
}

declare function this:opChamferIntersect(
  $a as item(), 
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  let $sqrt_2_2 := math:sqrt(2) div 2
  return
  ann:f("sdf:opChamferIntersect",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      return (
        max(($ad, $bd, ($ad + $k + $bd) * $sqrt_2_2)),
        if (head($av) > head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $k)
  )
}

declare function this:opChamferSubtract(
  $a as item(), 
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  let $sqrt_2_2 := math:sqrt(2) div 2
  return
  ann:f("sdf:opChamferSubtract",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      return (
        max(($ad, -$bd, ($ad + $k - $bd) * $sqrt_2_2)),
        if (head($av) > head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $k)
  )
}

declare function this:opSteppedUnion(
  $a as item(), 
  $b as item(),
  $r as xs:double,
  $n as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opSteppedUnion",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      return (
        this:steps($ad, $bd, $r, $n),
        if (head($av) > head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $r, $n)
  )
}

declare function this:opSteppedIntersect(
  $a as item(), 
  $b as item(),
  $r as xs:double,
  $n as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opSteppedIntersection",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      return (
        -this:steps(-$ad, -$bd, $r, $n),
        if (head($av) > head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $r, $n)
  )
}

declare function this:opSteppedSubtract(
  $a as item(), 
  $b as item(),
  $r as xs:double,
  $n as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opSteppedSubtract",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      return (
        -this:steps(-$ad, $bd, $r, $n),
        if (head($av) > head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $r, $n)
  )
}

declare function this:opTwist(
  $base as item(), 
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opTwist",
    function ($point as xs:double*) as xs:double* {
      let $c := math:cos($k * v:py($point))
      let $s := math:sin($k * v:py($point))
      let $matrix := ($c, -$s, $s, $c)
      let $q := (
        $matrix[1] * v:px($point) + $matrix[2] * v:pz($point),
        $matrix[3] * v:px($point) + $matrix[4] * v:pz($point),
        v:py($point)
      )
      return $basef($q)
    },
    ($base, $k)
  )
}

declare function this:opBend(
  $base as item(), 
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opBend",
    function ($point as xs:double*) as xs:double* {
      let $c := math:cos($k * v:px($point))
      let $s := math:sin($k * v:px($point))
      let $matrix := ($c, -$s, $s, $c)
      let $q := (
        $matrix[1] * v:px($point) + $matrix[2] * v:py($point),
        $matrix[3] * v:px($point) + $matrix[4] * v:py($point),
        tail(tail($point))
      )
      return $basef($q)
    },
    ($base, $k)
  )
}

declare function this:opAdd(
  $a as item(),
  $b as item()
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opAdd",
    function ($point as xs:double*) as xs:double* {
      $af($point)=>v:add($bf($point))
    },
    ($a, $b)
  )
}

declare function this:opAddInto(
  $a as item(),
  $b as item()
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opAddInto",
    function ($point as xs:double*) as xs:double* {
      let $v1 := $af($point)
      let $v2 := $bf($point)
      return (head($v1) + head($v2), tail($v1))
    },
    ($a, $b)
  )
}

declare function this:opInvert(
  $base as item()
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opInvert",
    function ($point as xs:double*) as xs:double* {
      let $res := $basef($point)
      return (-head($res), tail($res))
    },
    $base
  )
}

declare function this:opRescale(
  $base as item(),
  $amplitude as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opRescale",
    function ($point as xs:double*) as xs:double* {
      let $res := $basef($point)
      return ($amplitude * head($res), tail($res))
    },
    ($base, $amplitude)
  )
}

declare function this:opModifyDistance(
  $base as item(),
  $f as item() 
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  if (not(ann:is-distance-function($f))) then errors:error("ML-BADARGS", ("f", $f)) else (),
  let $basef := f:function($base)
  let $ff := f:function($f)
  return
  ann:f("sdf:opModifyDistance",
    function ($point as xs:double*) as xs:double* {
      let $res := $basef($point)
      return ($ff(head($res)), tail($res))
    },
    ($base, $f)
  )
}

declare function this:modAlmostIdentity(
  $m as xs:double,
  $n as xs:double
) as item()
{
  let $a := 2.0*$n - $m
  let $b := 2.0*$m - 3.0*$n
  return  
    ann:f("sdf:modAlmostIdentity", function ($x as xs:double) as xs:double {
      if (abs($x) > $m) then $x
      else (
        let $t := $x div $m
        return (($a * $t + $b)*$t*$t + $n)
      )
    })
}

declare function this:modExpImpulse(
  $k as xs:double
) as item()
{
  ann:f("sdf:modExpImpulse", function ($x as xs:double) as xs:double {
    let $d := abs($x)
    let $sign := util:sign($x)
    return (
      let $h := $k * $d
      return ($sign * $h * math:exp(1.0 - $h), tail($d))
    )
  })
}

declare function this:modQuadImpulse(
  $k as xs:double
) as item()
{
  let $two-root-k := 2.0 * math:sqrt($k)
  return
    ann:f("sdf:modQuadImpulse", function ($x as xs:double) as xs:double {
      let $d := abs($x)
      let $sign := util:sign($x)
      return (
        $sign * $two-root-k * $d div (1.0 + $k*$d*$d)
      )
    })
}

declare function this:modExpSustainedImpulse(
  $k as xs:double,
  $f as xs:double
) as item()
{
  ann:f("sdf:modExpSustainedImpulse", function ($x as xs:double) as xs:double {
    let $d := abs($x)
    let $sign := util:sign($x)
    let $s := max(($d - $f, 0.0))
    return (
      $sign * min(( ($d*$d) div ($f*$f), 1 + (2.0 div $f)*$s*math:exp(-$k*$s)))
    )
  })
}

declare function this:modCubicPulse(
  $c as xs:double,
  $w as xs:double
) as item()
{
  ann:f("sdf:modCubicPulse", function ($x as xs:double) as xs:double {
    let $x := abs($x - $c)
    return (
      if ($x > $w) then 0.0
      else (
        let $x := $x div $w
        return (1.0 - $x*$x*(3.0 - 2.0*$x))
      )
    )
  })
}

declare function this:modExpStep(
  $k as xs:double,
  $n as xs:double
) as item()
{
  ann:f("sdf:modExpStep", function ($x as xs:double) as xs:double {
    let $d := abs($x)
    let $sign := if ($x=0) then 1 else util:sign($x)
    return (
      $sign * math:exp( -$k * math:pow($d, $n) )
    )
  })
}

declare function this:modGain(
  $k as xs:double
) as item()
{
  ann:f("sdf:modGain", function ($x as xs:double) as xs:double {
    let $d := abs($x)
    let $sign := util:sign($x)
    let $a := 0.5 * math:pow(2.0*(if ($d < 0.5) then $d else 1.0 - $d), $k)
    return (
      if ($d < 0.5) then $sign*$a else $sign*(1.0 - $a)
    )
  })
}

declare function this:modPowerCurve(
  $a as xs:double,
  $b as xs:double
) as item()
{
  let $k := math:pow($a + $b, $a + $b) div (math:pow($a, $a) * math:pow($b, $b))
  return
    ann:f("sdf:modPowerCurve", function ($x as xs:double) as xs:double {
      let $d := abs($x)
      let $sign := util:sign($x)
      return (
        $sign * $k * math:pow($d, $a) * math:pow(1.0 - $d, $b)
      )
    })
}

declare function this:modSinc(
  $k as xs:double,
  $scale as xs:double
) as item()
{
  ann:f("sdf:modSinc", function ($x as xs:double) as xs:double {
    let $d := abs($x)
    let $sign := util:sign($x)
    let $a := math:pi()*($k * $d - 1.0)
    return (
      if ($a = 0) then $sign * $scale
      else $sign * $scale * math:sin($a) div $a
    )
  })
}

declare function this:modSinc(
  $k as xs:double
) as item()
{
  this:modSinc($k, 1.0)
}

declare function this:renderTrace(
  $f as item()
) as map(*)
{
  if (not(ann:is-renderer($f))) then errors:error("ML-BADARGS", ("f", $f)) else (),
  let $fn := f:function($f)
  return
  ann:f(util:function-name($f), 
    function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      $fn($tm, $origin, $direction, $sdf)=>trace(util:quote($f)||"("||util:quote($tm)||")")
    }
  )
}

declare function this:renderPerformance(
  $gradient as xs:string
) as map(*)
{
  let $colours := gradient:gradient-points($gradient)
  return 
    ann:f("sdf:renderPerformance", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      (: 64 = max iterations :)
      point:pcoordinates(
        cs:rgb-to-xyz(cs:interpolate-colour($tm[3] div 64, $colours))
      )
    })
}

declare function this:renderDepth(
  $gradient as xs:string,
  $field as xs:double*
) as map(*)
{
  let $colours := gradient:gradient-points($gradient)
  let $n := count($colours)
  let $scale := $n div ($field[2] - $field[1])
  return 
    ann:f("sdf:renderDepth", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      let $adjusted-t := util:round(($tm[1] - $field[1]) * $scale) div $n
      return
        point:pcoordinates(
          cs:rgb-to-xyz(cs:interpolate-colour($adjusted-t, $colours))
        )
    })
}

declare function this:renderColourPick(
  $pick-table as map(xs:string,item()*)*
) as map(*)
{
  let $n := count($pick-table)
  let $xyz-pick-table := $pick-table!cs:rgb-to-xyz(.)
  return
    ann:f("sdf:renderColourPick", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      point:pcoordinates($xyz-pick-table[util:modix(util:round($tm[2]), $n)])
    })
}

declare function this:renderMaterial(
  $pick-table as item()*
) as map(*)
{
  for $pick in $pick-table return if (not(ann:is-renderer($pick))) then errors:error("ML-BADARGS", ("pick-table", $pick)) else (),
  let $pick-tablef := for $pick in $pick-table return f:function($pick)
  let $n := count($pick-table)
  return
    ann:f("sdf:renderMaterial",
      function (
        $tm as xs:double*,
        $origin as xs:double*,
        $direction as xs:double*,
        $sdf as item()
      ) as xs:double*
      {
        let $f := $pick-tablef[util:modix(util:round($tm[2]), $n)]
        return $f($tm, $origin, $direction, $sdf)
      },
      $pick-table
    )
}

declare function this:renderChoice(
  $t-pivot as xs:double,
  $background as item(),
  $foreground as item()
) as map(*)
{
  if (not(ann:is-renderer($background))) then errors:error("ML-BADARGS", ("background", $background)) else (),
  if (not(ann:is-renderer($foreground))) then errors:error("ML-BADARGS", ("foreground", $foreground)) else (),
  let $foregroundf := f:function($foreground)
  let $backgroundf := f:function($background)
  return
  ann:f("sdf:renderChoice",
    function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      if (head($tm) <= $t-pivot)
      then $backgroundf($tm, $origin, $direction, $sdf)
      else $foregroundf($tm, $origin, $direction, $sdf)
    },
    ($t-pivot, $background, $foreground)
  )
}

declare function this:renderConstant(
  $colour as map(xs:string,item()*)
) as map(*)
{
  let $colourv := point:pcoordinates(cs:rgb-to-xyz($colour))
  return
    ann:f("sdf:renderConstant", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      $colourv
    })
}

declare function this:renderRandom(
  $colours as map(xs:string,item()*)*
) as map(*)
{
  let $xyz-colours := $colours!cs:rgb-to-xyz(.)
  return
    ann:f("sdf:renderRandom", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      rand:select-random($xyz-colours)!point:pcoordinates(.)
    })
}

declare function this:renderGradient(
  $base-colour as map(xs:string,item()*),
  $strength as xs:double
) as map(*)
{
  let $colourv := point:pcoordinates(cs:rgb-to-xyz($base-colour))
  return
    ann:f("sdf:renderGradient", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      let $shade := head($tm) * $strength
      return $colourv=>v:sub(($shade, $shade, $shade))
    })
}

declare function this:renderXGradient(
  $base-colour as map(xs:string,item()*),
  $strength as xs:double
) as map(*)
{
  let $colourv := point:pcoordinates(cs:rgb-to-xyz($base-colour))
  return
    ann:f("sdf:renderXGradient", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      let $shade := v:px($direction) * $strength
      return $colourv=>v:sub(($shade, $shade, $shade))
    })
}

declare function this:renderYGradient(
  $base-colour as map(xs:string,item()*),
  $strength as xs:double
) as map(*)
{
  let $colourv := point:pcoordinates(cs:rgb-to-xyz($base-colour))
  return
    ann:f("sdf:renderYGradient", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      let $shade := v:py($direction) * $strength
      return $colourv=>v:sub(($shade, $shade, $shade))
    })
}

declare function this:renderZGradient(
  $base-colour as map(xs:string,item()*),
  $strength as xs:double
) as map(*)
{
  let $colourv := point:pcoordinates(cs:rgb-to-xyz($base-colour))
  return
    ann:f("sdf:renderZGradient", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      let $shade := v:pz($direction) * $strength
      return $colourv=>v:sub(($shade, $shade, $shade))
    })
}

declare function this:opRenderFuzz(
  $base-colour as item(),
  $field as xs:double*
) as map(*)
{
  if (not(ann:is-renderer($base-colour))) then errors:error("ML-BADARGS", ("base-colour", $base-colour)) else (),
  let $base-colourf := f:function($base-colour)
  return
  ann:f("sdf:opRenderFuzz",
    function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      let $col := $base-colourf($tm, $origin, $direction, $sdf)
      let $smooth := util:smoothstep($field[1], $field[2], head($tm))
      return $col=>v:times(1.0 - $smooth)  
    },
    ($base-colour, $field)
  )
}

declare function this:opRenderDiffuseLight(
  $surface-colour as item(),
  $light as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-renderer($surface-colour))) then errors:error("ML-BADARGS", ("surface-colour", $surface-colour)) else (),
  let $surface-colourf := f:function($surface-colour)
  let $light-dir := light:direction($light)
  let $light-colour := light:directional($light)
  let $intensity := light:intensity($light)
  let $ambient-intensity := light:ambient-intensity($light)
  let $ambient := light:ambient($light)=>v:times($ambient-intensity)
  return
    ann:f("sdf:opRenderDiffuseLight",
      function (
        $tm as xs:double*,
        $origin as xs:double*,
        $direction as xs:double*,
        $sdf as item()
      ) as xs:double*
      {
        let $surface := $surface-colourf($tm, $origin, $direction, $sdf)
        let $pos := this:linear($origin, $direction, head($tm))
        let $n := this:calc-normal($pos, $sdf)
        let $surface-brightness := max((v:dot($n, $light-dir), 0.0)) * $intensity
        let $directional := $light-colour=>v:times($surface-brightness)
        return
          $surface=>v:multiply( $directional=>v:add($ambient) )
      },
      ($surface-colour, $light)
    )
}

declare function this:opRenderHardShadow(
  $base-colour as item(),
  $light as map(xs:string,item()*),
  $field as xs:double*,
  $shadow-scale as xs:double,
  $shadow-strength as xs:double
) as map(*)
{
  if (not(ann:is-renderer($base-colour))) then errors:error("ML-BADARGS", ("base-colour", $base-colour)) else (),
  let $base-colourf := f:function($base-colour)
  let $light-dir := light:direction($light)
  return
    ann:f("sdf:opRenderHardShadow",
      function (
        $tm as xs:double*,
        $origin as xs:double*,
        $direction as xs:double*,
        $sdf as item()
      ) as xs:double*
      {
        let $col := $base-colourf($tm, $origin, $direction, $sdf)
        let $pos := this:linear($origin, $direction, head($tm))
        let $normal := this:calc-normal($pos, $sdf)
        let $shadow-origin := $pos=>v:add($normal=>v:times($shadow-scale))
        let $st := this:cast-ray($shadow-origin, $light-dir, $field, $sdf)[2]
        let $shadow := if ($st != -1.0) then 1.0 else 0.0
        return v:mix($col, $col=>v:times($shadow-strength), $shadow)
      },
      ($base-colour, $light, $field, $shadow-scale, $shadow-strength)
    )
}

declare function this:opRenderSoftShadow(
  $base-colour as item(),
  $light as map(xs:string,item()*),
  $field as xs:double*,
  $shadow-scale as xs:double,
  $shadow-strength as xs:double,
  $n-shadows as xs:integer,
  $shadow-falloff as xs:double
) as map(*)
{
  if (not(ann:is-renderer($base-colour))) then errors:error("ML-BADARGS", ("base-colour", $base-colour)) else (),
  let $base-colourf := f:function($base-colour)
  let $light-dir := light:direction($light)
  let $rand := dist:uniform(-1.0, 1.0 - $config:ε)
  return
    ann:f("sdf:opRenderSoftShadow",
      function (
        $tm as xs:double*,
        $origin as xs:double*,
        $direction as xs:double*,
        $sdf as item()
      ) as xs:double*
      {
        let $col := $base-colourf($tm, $origin, $direction, $sdf)
        let $pos := this:linear($origin, $direction, head($tm))
        let $normal := this:calc-normal($pos, $sdf)
        let $shadow-origin := $pos=>v:add($normal=>v:times($shadow-scale))
        let $shadow :=
          fold-left(
            1 to $n-shadows,
            0.0,
            function ($shadow as xs:double, $s as xs:integer) as xs:double {
              let $fall := $shadow-falloff * rand:randomize($rand)
              let $shadow-dir := $light-dir=>v:add(($fall, $fall, $fall))
              let $st := this:cast-ray($shadow-origin, $shadow-dir, $field, $sdf)[2]
              return (
                if ($st != -1.0) then $shadow + 1.0 else $shadow
              )
            }
          )
        return v:mix($col, $col=>v:times($shadow-strength), $shadow div $n-shadows)
      },
      ($base-colour, $light, $field, $shadow-scale, $shadow-strength, $n-shadows, $shadow-falloff)
    )
}

declare function this:opRenderComplexLighting(
  $base-colour as item(),
  $light as map(xs:string,item()*), (: light source :)
  $shadow-scale as xs:double, (: 0.01 :)
  $k as xs:double (: softness 32.0 :)
) as map(*)
{
  if (not(ann:is-renderer($base-colour))) then errors:error("ML-BADARGS", ("base-colour", $base-colour)) else (),
  let $base-colourf := f:function($base-colour)
  let $light-direction := light:direction($light)
  let $light-colour :=
    light:directional($light)=>
      v:times( light:intensity($light) )=>
      v:add( light:ambient($light)=>v:times( light:ambient-intensity($light) ) )
  return
    ann:f("sdf:opRenderComplexLighting",
      function (
        $tm as xs:double*,
        $origin as xs:double*,
        $direction as xs:double*,
        $sdf as item()
      ) as xs:double*
      {
        let $pos := this:linear($origin, $direction, head($tm))
        let $col := $base-colourf($tm, $origin, $direction, $sdf)
        let $normal := this:calc-normal($pos, $sdf)
        let $shade-origin := $direction=>v:add($normal=>v:times($shadow-scale))
        let $shade := head(
          util:while(
            function ($res as xs:double, $t as xs:double, $i as xs:integer, $working as xs:boolean) as xs:boolean { $i <= 16 and $working },
            function ($res as xs:double, $t as xs:double, $i as xs:integer, $working as xs:boolean) as item()* {
              let $h := head($sdf( this:linear($shade-origin, $light-direction, $t)))
              let $res := min(($res, $k*$h div $t))
              let $t := $t + util:clamp($h, 0.01, 1.0) (: fineness of calc :)
              return (
                $res, $t, $i + 1, not($h < 0.0001)
              )
            },
            1.0, 0.0, 1, true()
          )
        )
        let $shade :=
          util:clamp($shade, 0.0, 1.0) *
          util:clamp(v:dot($normal, $light-direction), 0.0, 1.0)
        let $occlusion := head(
          util:while(
            function ($res as xs:double, $sca as xs:double, $i as xs:integer) as xs:boolean { $i < 5 },
            function ($res as xs:double, $sca as xs:double, $i as xs:integer) as item()* {
              let $hr := 0.02 + 0.025 * xs:double($i) * xs:double($i)
              let $aopos := $pos=>v:add($normal=>v:times($hr))
              let $dd := head($sdf($aopos))
              return ($res + -($dd - $hr)*$sca, $sca * 0.995, $i + 1)
            },
            0.0, 1.0, 0
          )
        )
        let $occlusion :=
          (1.0 - util:clamp($occlusion, 0.0, 1.0))*util:clamp(v:py($normal), 0.0, 1.0)
        return
          ($col=>v:times($occlusion))=>             (: base + ambient shadows :)
            v:add($light-colour=>v:times($shade))=> (: lighting :)
            v:times(math:exp(-0.2*head($tm)))       (: softening :)
      },
      ($base-colour, $light, $shadow-scale, $k)
    )
}

declare function this:opRenderColourWaves(
  $base-colour as item(),
  $darkness as xs:double,  (: -3.0; more negative=lighter :)
  $frequency as xs:double, (: 75; higher=more bands :)
  $strength as xs:double (: 0.2 [0,1]; lower=gentler bands :)
) as map(*)
{
  if (not(ann:is-renderer($base-colour))) then errors:error("ML-BADARGS", ("base-colour", $base-colour)) else (),
  let $base-colourf := f:function($base-colour)
  return
  ann:f("sdf:opRenderColourWaves",
    function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      let $d := head($tm)
      return
        $base-colourf($tm, $origin, $direction, $sdf)=>
          v:times(1.0 - math:exp($darkness*abs($d)))=>
          v:times((1 - $strength) + $strength*math:cos($frequency*$d))
    },
    ($base-colour, $darkness, $frequency, $strength)
  )
}

declare function this:opRenderColourBlend(
  $base-colour as item(),
  $blend-colour as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-renderer($base-colour))) then errors:error("ML-BADARGS", ("base-colour", $base-colour)) else (),
  let $base-colourf := f:function($base-colour)
  let $blend := point:pcoordinates($blend-colour)
  return
    ann:f("sdf:opRenderColourBlend",
      function (
        $tm as xs:double*,
        $origin as xs:double*,
        $direction as xs:double*,
        $sdf as item()
      ) as xs:double*
      {
        let $col := $base-colourf($tm, $origin, $direction, $sdf)
        return v:mix($col, $blend, $tm[1])
      },
      ($base-colour, $blend-colour)
    )
}

declare function this:opRenderModifyColour(
  $base-colour as item(),
  $blend-colour as map(xs:string,item()*),
  $f as item()
) as map(*)
{
  if (not(ann:is-renderer($base-colour))) then errors:error("ML-BADARGS", ("base-colour", $base-colour)) else (),
  if (not(ann:is-distance-function($f))) then errors:error("ML-BADARGS", ("f", $f)) else (),
  let $base-colourf := f:function($base-colour)
  let $ff := f:function($f)
  let $blend := point:pcoordinates($blend-colour)
  return
    ann:f("sdf:opRenderModifyColour",
      function (
        $tm as xs:double*,
        $origin as xs:double*,
        $direction as xs:double*,
        $sdf as item()
      ) as xs:double*
      {
        let $col := $base-colourf($tm, $origin, $direction, $sdf)
        let $mod := $ff($tm[1])
        return v:mix($col, $blend, $mod)
      },
      ($base-colour, $blend-colour, $f)
    )
}

declare function this:camera(
  $position as map(xs:string,item()*),
  $target as map(xs:string,item()*),
  $perspective as xs:double (: Z scaling :)
) as map(xs:string,item()*)
{
  let $forward := point:normalize(point:sub($target, $position))
  let $right := point:normalize(point:cross(point:point(0,1,0), $forward))
  let $up := point:normalize(point:cross($forward, $right))
  return
    map {
      "kind": "camera",
      "position": point:pcoordinates($position),
      "view": point:pcoordinates($position),
      "target": point:pcoordinates($target),
      "perspective": $perspective,
      "forward": point:pcoordinates($forward)=>v:times($perspective),
      "right": point:pcoordinates($right),
      "up": point:pcoordinates($up)
    }
}

declare function this:camera(
  $angle as xs:double,
  $view as map(xs:string,item()*),
  $target as map(xs:string,item()*),
  $perspective as xs:double (: Z scaling :)
) as map(xs:string,item()*)
{
  let $radians := util:radians($angle)
  let $origin :=
    point:pcoordinates($view)=>
      v:multiply((math:cos($radians), 1.0, math:sin($radians)))
  let $forward := v:normalize(point:pcoordinates($target)=>v:sub($origin))
  let $right := v:normalize(v:cross($forward,(0,1,0)))
  let $up := v:normalize(v:cross($forward, $right))
  return
    map {
      "kind": "camera",
      "position": $origin,
      "angle": $angle,
      "view": point:pcoordinates($view),
      "target": point:pcoordinates($target),
      "perspective": $perspective,
      "forward": $forward=>v:times($perspective),
      "right": $right,
      "up": $up
    }
}

declare function this:castDepth(
  $field as xs:double*
) as map(*)
{
  ann:f("sdf:castDepth", function ($origin as xs:double*, $direction as xs:double*, $sdf as item()) as xs:double* {
    this:cast-ray($origin, $direction, $field, $sdf)
  })
}

declare function this:castDirect(
) as map(*)
{
  ann:f("sdf:castDirect", function ($origin as xs:double*, $direction as xs:double*, $sdf as item()) as xs:double* {
    $sdf($direction)
  })
}

declare function this:render(
  $uv as map(xs:string,item()*), (: normalized point :)
  $camera as map(xs:string,item()*)?,
  $caster as item(),
  $sdf as item(),
  $render as item()
) as map(xs:string,item()*) (: rgba :)
{
  (:
   : This is called a lot: you're better off with make-point-render() 
   :)
  if (not(ann:is-caster($caster))) then errors:error("ML-BADARGS", ("caster", $caster)) else (),
  if (not(ann:is-sdf($sdf))) then errors:error("ML-BADARGS", ("sdf", $sdf)) else (),
  if (not(ann:is-renderer($render))) then errors:error("ML-BADARGS", ("render", $render)) else (),
  let $casterf := f:function($caster)
  let $sdff := f:function($sdf)
  let $renderf := f:function($render)
  return (
    if (empty($camera)) then (
      (: Take each point verbatim; really only useful for 2D :)
      let $ray-origin := point:pcoordinates($point:ORIGIN3D)
      let $ray-direction := point:pcoordinates($uv)
      let $tm := $casterf($ray-origin, $ray-direction, $sdff)
      let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
      return cs:xyz-to-rgb(point:vector($col))
    ) else (
      let $ray-origin := $camera("position")
      let $ray-direction := this:camera-ray-direction($uv, $camera)
      let $tm := $casterf($ray-origin, $ray-direction, $sdff)
      let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
      return cs:xyz-to-rgb(point:vector($col))
    )
  )
}

declare function this:make-point-render(
  $camera as map(xs:string,item()*)?,
  $caster as item(),
  $sdf as item(),
  $render as item()
) as function(map(xs:string,item()*)) as map(xs:string,item()*)
{
  if (not(ann:is-caster($caster))) then errors:error("ML-BADARGS", ("caster", $caster)) else (),
  if (not(ann:is-sdf($sdf))) then errors:error("ML-BADARGS", ("sdf", $sdf)) else (),
  if (not(ann:is-renderer($render))) then errors:error("ML-BADARGS", ("render", $render)) else (),
  let $casterf := f:function($caster)
  let $sdff := f:function($sdf)
  let $renderf := f:function($render)
  return (
    if (empty($camera)) then (
      let $ray-origin := point:pcoordinates($point:ORIGIN3D) return
      function ($uv as map(xs:string,item()*)) as map(xs:string,item()*) {
        let $ray-direction := point:pcoordinates($uv)
        let $tm := $casterf($ray-origin, $ray-direction, $sdff)
        let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
        return cs:xyz-to-rgb(point:vector($col))
      }
    ) else (
      let $ray-origin := $camera("position") return
      function ($uv as map(xs:string,item()*)) as map(xs:string,item()*) {
        let $ray-direction := this:camera-ray-direction($uv, $camera)
        let $tm := $casterf($ray-origin, $ray-direction, $sdff)
        let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
        return cs:xyz-to-rgb(point:vector($col))
      }
    )
  )
}

declare function this:renderv(
  $uv as xs:double*, (: normalized point :)
  $camera as map(xs:string,item()*)?,
  $caster as item(),
  $sdf as item(),
  $render as item()
) as map(xs:string,item()*) (: rgba :)
{
  (:
   : This is called a lot: you're better off with make-vector-render() 
   :)
  if (not(ann:is-caster($caster))) then errors:error("ML-BADARGS", ("caster", $caster)) else (),
  if (not(ann:is-sdf($sdf))) then errors:error("ML-BADARGS", ("sdf", $sdf)) else (),
  if (not(ann:is-renderer($render))) then errors:error("ML-BADARGS", ("render", $render)) else (),
  let $casterf := f:function($caster)
  let $sdff := f:function($sdf)
  let $renderf := f:function($render)
  return (
    if (empty($camera)) then (
      (: Take each point verbatim; really only useful for 2D :)
      let $ray-origin := point:pcoordinates($point:ORIGIN3D)
      let $ray-direction := $uv
      let $tm := $casterf($ray-origin, $ray-direction, $sdff)
      let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
      return cs:xyz-to-rgb(point:vector($col))
    ) else (
      let $ray-origin := $camera("position")
      let $ray-direction := this:camera-ray-directionv($uv, $camera)
      let $tm := $casterf($ray-origin, $ray-direction, $sdff)
      let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
      return cs:xyz-to-rgb(point:vector($col))
    )
  )
}

declare function this:make-vector-render(
  $camera as map(xs:string,item()*)?,
  $caster as item(),
  $sdf as item(),
  $render as item()
) as function(xs:double*) as map(xs:string,item()*)
{
  if (not(ann:is-caster($caster))) then errors:error("ML-BADARGS", ("caster", $caster)) else (),
  if (not(ann:is-sdf($sdf))) then errors:error("ML-BADARGS", ("sdf", $sdf)) else (),
  if (not(ann:is-renderer($render))) then errors:error("ML-BADARGS", ("render", $render)) else (),
  let $casterf := f:function($caster)
  let $sdff := f:function($sdf)
  let $renderf := f:function($render)
  return (
    if (empty($camera)) then (
      let $ray-origin := point:pcoordinates($point:ORIGIN3D) return
      function ($uv as xs:double*) as map(xs:string,item()*) {
        let $ray-direction := $uv
        let $tm := $casterf($ray-origin, $ray-direction, $sdff)
        let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
        return cs:xyz-to-rgb(point:vector($col))
      }
    ) else (
      let $ray-origin := $camera("position") return
      function ($uv as xs:double*) as map(xs:string,item()*) {
        let $ray-direction := this:camera-ray-directionv($uv, $camera)
        let $tm := $casterf($ray-origin, $ray-direction, $sdff)
        let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
        return cs:xyz-to-rgb(point:vector($col))
      }
    )
  )
}

declare function this:camera-ray-direction(
  $uv as map(xs:string,item()*),
  $camera as map(xs:string,item()*)
) as xs:double*
{
  v:normalize(
    (
      ($camera("right"))=>v:times(point:px($uv))
    )=>v:add(
      ($camera("up"))=>v:times(point:py($uv))
    )=>v:add(
      $camera("forward")
    )
  )
}

declare function this:camera-ray-directionv(
  $uv as xs:double*,
  $camera as map(xs:string,item()*)
) as xs:double*
{
  v:normalize(
    (
      ($camera("right"))=>v:times(v:px($uv))
    )=>v:add(
      ($camera("up"))=>v:times(v:py($uv))
    )=>v:add(
      $camera("forward")
    )
  )
}

declare function this:cast-ray(
  $ray-origin as xs:double*,
  $ray-direction as xs:double*,
  $field as xs:double*, (: near field, far field :)
  $sdf as item()
) as xs:double* (: distance, material, iterations :)
{
  util:assert(count($ray-origin)=3, "bad ray origin"),
  (: util:assert(count($ray-direction)=3, "bad ray direction"), :)
  let $ray-direction := if (count($ray-direction)=2) then ($ray-direction,0) else $ray-direction
  let $ε := 0.0001
  let $start := $field[1]
  let $end := $field[2]
  let $max-iterations := 64
  let $data :=
    util:while(
      function ($working as xs:boolean, $t as xs:double, $m as xs:double, $i as xs:integer) as xs:boolean {
        $working and $i <= $max-iterations
      },
      function ($working as xs:boolean, $t as xs:double, $m as xs:double, $i as xs:integer) {
        let $pos := this:linear($ray-origin, $ray-direction, $t)
        let $result := $sdf($pos)
        let $distance := $result[1]
        let $material := ($result[2],0)[1]
        return (
          if ($distance < $ε) then (
            false(), $t, $m, $i + 1
          ) else if ($t + $distance > $end) then (
            false(), -1.0, -1.0, $i + 1
          ) else (
            true(), $t + $distance, $material, $i + 1
          )
        )
      },
      true(), $start, -1.0, 1
    )
  let $working := head($data)
  return
    if ($working) then (-1.0, -1.0, $data[4]) else tail($data)
}

declare function this:calc-normal(
  $point as xs:double*,
  $sdf as item()
) as xs:double*
{
  let $c := head($sdf($point))
  let $epsilon := 0.001
  (: Compare sdf in offsets from point :)
  return
    v:sub(
      v:normalize((
        head($sdf($point=>v:add(($epsilon, 0, 0)))),
        head($sdf($point=>v:add((0, $epsilon, 0)))),
        head($sdf($point=>v:add((0, 0, $epsilon))))
      )),
      ($c, $c, $c)
    )
}

declare function this:linear(
  $origin as xs:double*,
  $direction as xs:double*,
  $t as xs:double
) as xs:double*
{
  util:zip(
    function ($a as xs:double, $b as xs:double) as xs:double {$a + $t * $b},
      $origin,
      $direction
  )
}

declare function this:normalize-to-canvas(
  $pt as map(xs:string,item()*),
  $canvas as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $width := box:width($canvas)
  let $height := box:height($canvas)
  return (
    point:vector((
      2.0 * (point:px($pt) - ($width div 2)) div $width,
      2.0 * (point:py($pt) - ($height div 2)) div $height,
      point:pcoordinates($pt)[position() > 2]
    ))
  )
}

declare function this:normalize-to-canvasv(
  $pt as xs:double*,
  $canvas as map(xs:string,item()*)
) as xs:double*
{
  let $width := box:width($canvas)
  let $height := box:height($canvas)
  return (
    (
      2.0 * (v:px($pt) - ($width div 2)) div $width,
      2.0 * (v:py($pt) - ($height div 2)) div $height,
      ($pt)[position() > 2]
    )
  )
}

declare function this:normalize-to-canvas(
  $x as xs:double,
  $y as xs:double,
  $width as xs:double,
  $height as xs:double,
  $aspect as xs:double
) as map(xs:string,item()*)
{
  point:point(
    2.0 * ($x - ($width div 2)) div $width,
    2.0 * ($y - ($height div 2)) div ($aspect * $height)
  )
}

declare function this:normalize-to-canvasv(
  $x as xs:double,
  $y as xs:double,
  $width as xs:double,
  $height as xs:double,
  $aspect as xs:double
) as xs:double*
{
  (
    2.0 * ($x - ($width div 2)) div $width,
    2.0 * ($y - ($height div 2)) div ($aspect * $height)
  )
}

declare function this:scale-from-canvas(
  $d as xs:double,
  $canvas as map(xs:string,item()*)
) as xs:double
{
  if (empty($canvas)) then $d else
  let $width := $canvas("width")
  return (
    2.0 * $d div $width
  )
}

declare function this:scale-from-unit(
  $d as xs:double,
  $canvas as map(xs:string,item()*)
) as xs:double
{
  let $width := $canvas("width")
  return (
    0.5 * $d * $width
  )
}

declare function this:toSDF(
  $regions as map(xs:string,item()*)*,
  $options as map(xs:string,item()*)
) as map(*)
{
  this:toSDF($regions, (), $options)
}

declare function this:toSDF(
  $regions as map(xs:string,item()*)*,
  $canvas as map(xs:string,item()*)?,
  $options as map(xs:string,item()*)
) as map(*)
{
  let $options :=
    util:merge-into((
      map {
        "op": "union",
        "smoothing": 0.05,
        "interpolate": 0,
        "steps": 3,
        "radius": 1,
        "width": 1,
        "ud": false()
      },
      $options
    ))
  let $op := $options("op")
  let $k := $options("smoothing")
  let $n := $options("steps")
  let $translated-regions :=
    if (exists($canvas)) then (
      for $region in $regions return (
        this:regionToSDF($region, $op, $canvas, $options)
      )
    ) else (
      for $region in $regions return (
        this:regionToSDF($region, $op, $options)
      )
    )
  return (
    if ($op="union") then (
      this:opUnion($translated-regions)
    ) else if ($op="intersect") then (
      this:opIntersect($translated-regions)
    ) else (
      fold-left(
        tail($translated-regions), head($translated-regions),
        function ($sdf as map(*), $current as map(*)) as map(*) {
          if (empty($sdf)) then $current
          else switch($op)
          case "intersect" return this:opIntersect($sdf, $current)
          case "subtract" return this:opSubtract($sdf, $current)
          case "round-union" return this:opRoundUnion($sdf, $current, $k)
          case "round-intersect" return this:opRoundIntersect($sdf, $current, $k)
          case "round-subtract" return this:opRoundSubtract($sdf, $current, $k)
          case "smooth-union" return this:opSmoothUnion($sdf, $current, $k)
          case "smooth-intersect" return this:opSmoothIntersect($sdf, $current, $k)
          case "smooth-subtract" return this:opSmoothSubtract($sdf, $current, $k)
          case "chamfer-union" return this:opChamferUnion($sdf, $current, $k)
          case "chamfer-intersect" return this:opChamferIntersect($sdf, $current, $k)
          case "chamfer-subtract" return this:opChamferSubtract($sdf, $current, $k)
          case "stepped-union" return this:opSteppedUnion($sdf, $current, $k, $n)
          case "stepped-intersect" return this:opSteppedIntersect($sdf, $current, $k, $n)
          case "stepped-subtract" return this:opSteppedSubtract($sdf, $current, $k, $n)
          default (: union :) return this:opUnion($sdf, $current)
        }
      )
    )
  )
}

Original Source Code

xquery version "3.1";
(:~
 : SDF support functions
 :
 : SDF modifiers 
 : After https://www.iquilezles.org/www/articles/distfunctions/distfunctions.htm
 : and https://iquilezles.org/www/articles/functions/functions.htm
 : Also some aspects from https://github.com/electricsquare/raymarching-workshop
 : See also https://www.alanzucconi.com/2016/07/01/surface-shading/ for more
 : information on surface lighting models
 : Some operators from https://github.com/thi-ng/umbrella
 :
 : Modifier functions, like SDF functions, return a function that 
 : takes a point and returns the SDF+material.
 :
 : Rendering functions also return functions which take the vector returned
 : by the ray caster (distance, material, iterations), the ray origin and 
 : direction, and the overall SDF function (needed for normal calculations etc.)
 : Rendering functions return colours in linear (XYZ) colour space even
 : though for the most part they take RGB colours as constant parameters.
 : The main render function will need to map back to RGB at the end.
 :
 : As of 202302 these functions do not return annotated functions, but wrapped
 : callables. You'll need to use callable:function() to get the actual function to use.
 : sdf:make-point-render() and sdf:make-vector-render() take care of this for you, 
 : but when using SDFs outside of ray casting contexts, you'll have to manage this.
 :
 : Copyright© Mary Holstege 2021-2023
 : CC-BY (https://creativecommons.org/licenses/by/4.0/)
 : @since October 2021
 : @custom:Status Active
 :)
module namespace this="http://mathling.com/sdf";

declare namespace art="http://mathling.com/art";
declare namespace map="http://www.w3.org/2005/xpath-functions/map";
declare namespace array="http://www.w3.org/2005/xpath-functions/array";
declare namespace math="http://www.w3.org/2005/xpath-functions/math";

import module namespace errors="http://mathling.com/core/errors"
       at "../core/errors.xqy";
import module namespace config="http://mathling.com/core/config"
       at "../core/config.xqy";
import module namespace rand="http://mathling.com/core/random"
       at "../core/random.xqy";
import module namespace dist="http://mathling.com/type/distribution"
       at "../types/distributions.xqy";
import module namespace wrapper="http://mathling.com/type/wrapper"
       at "../types/wrapper.xqy";
import module namespace slot="http://mathling.com/type/slot"
       at "../types/slot.xqy";
import module namespace light="http://mathling.com/type/light"
       at "../types/light.xqy";
import module namespace util="http://mathling.com/core/utilities"
       at "../core/utilities.xqy";
import module namespace f="http://mathling.com/core/callable"
       at "../core/callable.xqy";
import module namespace v="http://mathling.com/core/vector"
       at "../core/vector.xqy";
import module namespace affine="http://mathling.com/geometric/affine"
       at "../geo/affine.xqy";
import module namespace point="http://mathling.com/geometric/point"
       at "../geo/point.xqy";
import module namespace box="http://mathling.com/geometric/rectangle"
       at "../geo/rectangle.xqy";
import module namespace edge="http://mathling.com/geometric/edge"
       at "../geo/edge.xqy";
import module namespace path="http://mathling.com/geometric/path"
       at "../geo/path.xqy";
import module namespace cpoly="http://mathling.com/geometric/complex-polygon"
       at "../geo/complex-polygon.xqy";
import module namespace ellipse="http://mathling.com/geometric/ellipse"
       at "../geo/ellipse.xqy";
import module namespace solid="http://mathling.com/geometric/solid"
       at "../geo/solid.xqy";
import module namespace gradient="http://mathling.com/svg/gradients"
       at "../svg/gradients.xqy";
import module namespace rgb="http://mathling.com/colour/rgb"
       at "../colourspace/rgb.xqy";
import module namespace xyz="http://mathling.com/colour/xyz"
       at "../colourspace/xyz.xqy";
import module namespace cs="http://mathling.com/colour/space"
       at "../colourspace/colour-space.xqy";
import module namespace sd2="http://mathling.com/sdf/2d"
       at "../sdf/sd2.xqy";
import module namespace sd3="http://mathling.com/sdf/3d"
       at "../sdf/sd3.xqy";
import module namespace ann="http://mathling.com/sdf/annotations"
       at "../sdf/annotations.xqy";

(:======================================================================
 : SDF functions
 :======================================================================:)

(:~
 : trace()
 : Pass through to function, with tracing
 :)
declare function this:trace(
  $f as item()
) as map(*)
{
  if (not(ann:is-sdf($f))) then errors:error("ML-BADARGS", ("f", $f)) else (),
  let $fn := f:function($f)
  return
  ann:f(util:function-name($f), 
    function ($point as xs:double*) as xs:double* {
      $fn($point)=>trace(util:quote($f)||"("||util:quote($point)||")")
    }
  )
};

(:~
 : as-point()
 : Pass through to point-as-point function, converting from vector
 : $f is function(map(xs:string,item()*)) as map(xs:string,item()*) or a callable
 : wrapper of same
 : Result is callable function(xs:double*) as xs:double*
 :)
declare function this:as-point(
  $f as item()
) as map(*)
{
  let $fn := f:function($f)
  return
  ann:f(util:function-name($f), 
    function ($point as xs:double*) as xs:double* {
      $fn(point:vector($point))=>point:pcoordinates()
    }
  )
};

(:~
 : as-vector()
 : Pass through to make a point-as-point function, converting from vector
 : Can be used to convert these functions to mutation functions, for example
 : $f is function(xs:double*) as xs:double* or a f: wrapper of same
 : Result is callable function(map(xs:string,item()*)) as map(xs:string,item()*)
 :)
declare function this:as-vector(
  $f as item()
) as map(*)
{
  let $fn := f:function($f)
  return
  ann:f(util:function-name($f),
    function ($point as map(xs:string,item()*)) as map(xs:string,item()*) {
      $fn(point:pcoordinates($point))=>point:vector()
    }
  )
};

(:~
 : opConstant()
 : Return a constant value, for any point.
 : @param $k: the constant value
 :)
declare function this:opConstant(
  $k as xs:double
) as map(*)
{
  ann:f("sdf:opConstant", function ($point as xs:double*) as xs:double* {
    $k
  })
};

(:~
 : opMaterial()
 : Append a constant value as an additional coordinate.
 : @param $material: the constant value
 :)
declare function this:opMaterial(
  $base as item(),
  $material as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opMaterial",
    function ($point as xs:double*) as xs:double* {
      ($basef($point), $material)
    },
    ($base, $material)
  )
};

(:~
 : opTranslate()
 : Translate the point before applying the base function.
 : This means the sign is the opposite of what you expect: to move a box,
 : say, to center at (1,1) do opTranslate(-1,-1)
 :)
declare function this:opTranslate(
  $base as item(),
  $loc as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $matrix := affine:translation3(point:px($loc), point:py($loc), point:pz($loc))
  let $basef := f:function($base)
  return
    ann:f("sdf:opTranslate", function ($point as xs:double*) as xs:double* {
        $basef(affine:affine3($point, $matrix))
      },
      ($base, $loc)
    )
};

(:~
 : opRotate()
 : Rotate (2D) the point before applying the base function.
 :)
declare function this:opRotate(
  $base as item(), 
  $rotation as xs:double (: degrees :)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $matrix := affine:rotation2($rotation)
  let $basef := f:function($base)
  return
    ann:f("sdf:opRotate",
      function ($point as xs:double*) as xs:double* {
        $basef(affine:affine2($point, $matrix))
      },
      ($base, $rotation)
    )
};

(:~
 : opRotate()
 : Rotate (3D) the point before applying the base function.
 :)
declare function this:opRotate(
  $base as item(),
  $roll as xs:double, (: all degrees :)
  $pitch as xs:double,
  $yaw as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $matrix := affine:rotation3($roll, $pitch, $yaw)
  let $basef := f:function($base)
  return
    ann:f("sdf:opRotate",
      function ($point as xs:double*) as xs:double* {
        $basef(affine:affine3($point=>v:as-dimension(3), $matrix))
      },
      ($base, $roll, $pitch, $yaw)
    )
};

(:~
 : opScale()
 : Scale the point before applying the base function.
 :)
declare function this:opScale(
  $base as item(), 
  $scale as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opScale",
    function ($point as xs:double*) as xs:double* {
      $basef($point=>v:div($scale))=>v:times($scale)
    },
    ($base, $scale)
  )
};

(:~
 : opReflect()
 : Reflect the point across a center before applying the base function.
 :)
declare function this:opReflect(
  $base as item(),
  $center as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $matrix := affine:reflection3($center)
  let $basef := f:function($base)
  return
    ann:f("sdf:opReflect",
      function ($point as xs:double*) as xs:double* {
        $basef(affine:affine3($point, $matrix))
      },
      ($base, $center)
    )
};

(:~
 : opReflect()
 : Reflect the point X and Z coordinates across a line before applying the base function.
 :)
declare function this:opReflectXZ(
  $base as item(),
  $start as map(xs:string,item()*),
  $end as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $matrix := affine:reflection2($start, $end)
  let $basef := f:function($base)
  return
    ann:f("sdf:opReflectXZ",
      function ($point as xs:double*) as xs:double* {
        let $qxz := affine:affine2((v:px($point), v:pz($point)), $matrix)
        return $basef((v:px($qxz), v:py($point), v:pz($qxz)))
      },
      ($base, $start, $end)
    )
};

(:~
 : opReflect()
 : Reflect the point X and Y coordinates across a line before applying the base function.
 :)
declare function this:opReflectXY(
  $base as item(),
  $start as map(xs:string,item()*),
  $end as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $matrix := affine:reflection2($start, $end)
  let $basef := f:function($base)
  return
    ann:f("sdf:opReflectXY",
      function ($point as xs:double*) as xs:double* { 
        let $qxy := affine:affine2((v:px($point), v:py($point)), $matrix)
        return $basef((v:px($qxy), v:py($qxy), v:pz($point)))
      },
      ($base, $start, $end)
    )
};

(:~
 : opReflect()
 : Reflect the point Y and Z coordinates across a line before applying the base function.
 :)
declare function this:opReflectYZ(
  $base as item(),
  $start as map(xs:string,item()*),
  $end as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $matrix := affine:reflection2($start, $end)
  let $basef := f:function($base)
  return
    ann:f("sdf:opReflectYZ",
      function ($point as xs:double*) as xs:double* {
        let $qyz := affine:affine2((v:py($point), v:pz($point)), $matrix)
        return $basef((v:px($point), v:py($qyz), v:pz($qyz)))
      },
      ($base, $start, $end)
    )
};

(:~
 : opRepeat()
 : Repeat the base function with the given periods (corresponding coordinates of
 : periods point).
 :)
declare function this:opRepeat(
  $base as item(),
  $periods as map(xs:string,item()*) 
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $c := point:pcoordinates($periods)
  let $half-c := $c=>v:times(0.5)
  let $basef := f:function($base)
  return
    ann:f("sdf:opRepeat",
      function ($point as xs:double*) as xs:double* {
        let $q := $point=>v:add($half-c)=>v:modulo($c)=>v:sub($half-c)
        return $basef($q)
      },
      ($base, $periods)
    )
};

(:~
 : opLimitedRepeat()
 : Repeat the base function with the given periods (corresponding coordinates of
 : periods point), with repetitions limited by by lower-limits and upper-limits.
 : we get lower-limits.x to upper-limits.x by lower-limits.y to upper-limits.y
 :)
declare function this:opLimitedRepeat(
  $base as item(),
  $periods as map(xs:string,item()*),
  $lower-limits as map(xs:string,item()*),
  $upper-limits as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $c := point:pcoordinates($periods)
  let $lim-a := point:pcoordinates($lower-limits)
  let $lim-b := point:pcoordinates($upper-limits)
  let $basef := f:function($base)
  return
    ann:f("sdf:opLimitedRepeat",
      function ($point as xs:double*) as xs:double* {
        $basef(
          $point=>v:sub(
            v:clampv(v:round($point=>v:divide($c)), $lim-a, $lim-b)=>v:multiply($c)
          )
        )
      },
      ($base, $periods, $lower-limits, $upper-limits)
    )
};

(:~ OpenGL mod: different from XQuery for negative numbers :)
declare function this:mod($x as xs:double, $k as xs:double) as xs:double
{
  $x - $k * floor($x div $k)
};

declare function this:vmod($v as xs:double*, $k as xs:double) as xs:double*
{
  v:map(function ($c as xs:double) as xs:double {this:mod($c, $k)}, $v)
};

(:~
 : opRepeat2()
 :)
declare function this:opRepeat2(
  $base as item(),
  $size as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  let $s2 := $size div 2
  return
    ann:f("sdf:opRepeat2",
      function ($point as xs:double*) as xs:double* {
        $basef($point=>v:plus($s2)=>this:vmod($size)=>v:minus($s2))
      },
      ($base, $size)
    )
};

(:~
 : opMirrorRepeat2()
 :)
declare function this:opMirrorRepeat2(
  $base as item(),
  $size as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  let $s2 := $size div 2
  return
    ann:f("sdf:opMirrorRepeat2",
      function ($point as xs:double*) as xs:double* {
        let $p2 := $point=>v:plus($s2)
        return
        $basef(
         v:multiply(
           ($p2=>this:vmod($size)=>v:minus($s2)),
           ($p2=>v:div($size)=>v:floor()=>this:vmod(2)=>v:times(2)=>v:minus(1))
			)
        )
      },
      ($base, $size)
    )
};

(:~
 : opGridRepeat2()
 :)
declare function this:opGridRepeat2(
  $base as item(),
  $size as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  let $s2 := $size div 2
  let $mirrorf := (: opMirrorRepeat2 with no base function :)
    function ($point as xs:double*) as xs:double* {
      let $p2 := $point=>v:plus($s2)
      return
       v:multiply(
         ($p2=>this:vmod($size)=>v:minus($s2)),
         ($p2=>v:div($size)=>v:floor()=>this:vmod(2)=>v:times(2)=>v:minus(1))
       )
    }
  return
    ann:f("sdf:opGridRepeat2",
      function ($point as xs:double*) as xs:double* {
        let $rm := $mirrorf($point)
        let $p := $rm=>v:minus($s2)
        let $p := if (v:px($p) > v:py($p)) then (v:py($p), v:px($p)) else $p
        return (
          $basef($p)
        )
      },
      ($base, $size)
    )
};

(:~
 : opPolarRepeat2()
 :)
declare function this:opPolarRepeat2(
  $base as item(),
  $n as xs:integer
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  let $θ := 2 * math:pi() div $n
  let $halfθ := $θ div 2
  return
    ann:f("sdf:opPolarRepeat2",
      function ($point as xs:double*) as xs:double* {
        let $α := this:mod(math:atan2(v:py($point), v:px($point)) + $halfθ, $θ) - $halfθ
        let $p := (math:cos($α), math:sin($α))=>v:times(v:magnitude($point))
        return $basef($p)
      },
      ($base, $n)
    )
};

(:~
 : opElongate()
 : Extends a base shape by h in the given dimensions
 : e.g. h=(2,0) will extend all x dimensions outward by 2 (+/-)
 : @param $h gives the amount of elongation in each dimension
 :)
declare function this:opElongate(
  $base as item() ,
  $h as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $hv := point:pcoordinates($h)
  let $basef := f:function($base)
  return
    ann:f("sdf:opElongate",
      function ($point as xs:double*) as xs:double* {
        let $q := $point=>v:abs()=>v:sub($hv)
        return $basef($q=>v:at-max(0)) + min((max($q),0))
      },
      ($base, $h)
    )
};

(:~
 : opRound()
 : Round the shape, with radius of rounding r
 :)
declare function this:opRound(
  $base as item(),
  $r as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opRound",
    function ($point as xs:double*) as xs:double* {
      $basef($point) - $r
    },
    ($base, $r)
  )
};

(:~
 : opOnion()
 : Create onion ring of the basic shape; kind of a scooping out of a smaller
 : version of the outer shape.
 : To create multiple layers, iterate onion
 : @param $thickness is width of layer between inner and outer shape; core shape is
 : extended +/- by this
 :)
declare function this:opOnion(
  $base as item(),
  $thickness as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $r := $thickness div 2
  let $basef := f:function($base)
  return
    ann:f("sdf:opOnion",
      function ($point as xs:double*) as xs:double* {
        abs($basef($point)) - $r
      },
      ($base, $thickness)
    )
};

declare function this:opManyOnion(
  $base as item(),
  $thickness as xs:double,
  $n as xs:integer
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  fold-left(
    1 to $n,
    $base,
    function ($f as item(), $i as xs:integer) as map(*) {
      this:opOnion($f, $thickness*math:pow(2, -($i - 1)))
    }
  )
};

(:~
 : opExtrude()
 : Convert a 2D shape to 3D by extruding into the third dimension by h
 :)
declare function this:opExtrude(
  $base as item(),  (: 2D :) 
  $h as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opExtrude",
    function ($point as xs:double*) as xs:double* {
      let $d := $basef($point=>v:as-dimension(2))
      let $w := ( $d, abs(v:pz($point)) - $h )
      return min((max($w), 0)) + v:magnitude($w=>v:at-max(0))
    },
    ($base, $h)
  )
};

(:~
 : opRevolve()
 : Convert a 2D shape into a 3D one by revolving around y axis at $o distance
 :)
declare function this:opRevolve(
  $base as item(), (: 2D :)
  $o as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opRevolve",
    function ($point as xs:double*) as xs:double* {
      $basef( ( v:magnitude((v:px($point), v:pz($point))) - $o, v:py($point) ) )
    },
    ($base, $o)
  )
};

(:~
 : opUnion()
 : Union of two shapes.
 :)
declare function this:opUnion(
  $a as item(),
  $b as item()
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opUnion",
    function ($point as xs:double*) as xs:double* {
      let $v1 := $af($point)
      let $v2 := $bf($point)
      return
        if (head($v1) < head($v2)) then $v1 else $v2
    },
    ($a, $b)
  )
};

(:~
 : opUnion()
 : Union of multiple shapes.
 :)
declare function this:opUnion(
  $fs as item()*
) as map(*)
{
  for $f in $fs return if (not(ann:is-sdf($f))) then errors:error("ML-BADARGS", ("fs", $f)) else (),
  let $ffs := for $f in $fs return f:function($f)
  return
  ann:f("sdf:opUnion",
    function ($point as xs:double*) as xs:double* {
      let $ix :=
        util:min-index($ffs,
          function ($f as function(xs:double*) as xs:double*) as xs:double {
            head($f($point))
          }
        )
      return (
        $ffs[$ix]($point)
      )
    },
    $fs
  )
};

declare function this:min(
  $a as xs:double,
  $b as xs:double
) as xs:double
{
  min(($a,$b))
};

declare function this:sminQuad(
  $a as xs:double,
  $b as xs:double,
  $k as xs:double (: 0.1 :)
) as xs:double
{
  let $h := max(($k - abs($a - $b), 0)) div $k
  return min(($a,$b)) - $h*$h*$k*(1.0 div 4.0)
};

declare function this:sminCubic(
  $a as xs:double,
  $b as xs:double,
  $k as xs:double (: 0.1 :)
) as xs:double
{
  let $h := max(($k - abs($a - $b), 0)) div $k
  return min(($a,$b)) - $h*$h*$h*$k*(1.0 div 6.0)
};

declare function this:opUnion(
  $a as item(),
  $b as item(),
  $min as function(xs:double,xs:double) as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opUnion",
    function ($point as xs:double*) as xs:double* {
      let $v1 := $af($point)
      let $v2 := $bf($point)
      return (
        $min(head($v1), head($v2)),
        if (head($v1) < head($v2)) then tail($v1) else tail($v2)
      )
    },
    ($a, $b, $min)
  )
};

(:~
 : opSubtract()
 : Difference of two shapes. a - b
 :)
declare function this:opSubtract(
  $a as item(),
  $b as item()
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opSubtract",
    function ($point as xs:double*) as xs:double* {
      let $v1 := $af($point)
      let $v2 := $bf($point)
      return (
        max((head($v1), -head($v2))),
        if (head($v1) > -head($v2)) then tail($v1) else tail($v2)
      )
    },
    ($a, $b)
  )
};

(:~
 : opIntersect()
 : Intersection of two shapes.
 :)
declare function this:opIntersect(
  $a as item(),
  $b as item()
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opIntersect",
    function ($point as xs:double*) as xs:double* {
      let $v1 := $af($point)
      let $v2 := $bf($point)
      return (
        if (head($v1) > head($v2)) then $v1 else $v2
      )
    },
    ($a, $b)
  )
};

(:~
 : opIntersect()
 : Intersection of multiple shapes.
 :)
declare function this:opIntersect(
  $fs as item()*
) as map(*)
{
  for $f in $fs return if (not(ann:is-sdf($f))) then errors:error("ML-BADARGS", ("fs", $f)) else (),
  let $ffs := for $f in $fs return f:function($f)
  return
  ann:f("sdf:opIntersect",
    function ($point as xs:double*) as xs:double* {
      let $ix :=
        util:max-index($ffs,
          function ($f as function(xs:double*) as xs:double*) as xs:double {
            head($f($point))
          }
        )
      return (
        $ffs[$ix]($point)
      )
    },
    ($fs)
  )
};

(:~
 : opRoundUnion()
 : Rounded union of two shapes
 :)
declare function this:opRoundUnion(
  $a as item(), 
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opRoundUnion",
    function ($point as xs:double*) as xs:double* {
      let $v2 := $af($point)
      let $v1 := $bf($point)
      let $d1 := head($v1)
      let $d2 := head($v2)
      let $h := max(($k - abs($d1 - $d2), 0))
      return (
        min(($d1, $d2)) - $h*$h*0.25 div $k,
        if (head($v2) < head($v1)) then tail($v2) else tail($v1)
      )
    },
    ($a, $b, $k)
  )
};

(:~
 : opRoundSubtract()
 : Rounded subtraction of two shapes a - b
 :)
declare function this:opRoundSubtract(
  $a as item(),
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opRoundSubtract",
    function ($point as xs:double*) as xs:double* {
      let $l := v:magnitude($point)
      let $v2 := $af($point)
      let $v1 := $bf($point)
      let $d1 := head($v1)
      let $d2 := head($v2)
      let $h :=  max(($k - abs(-$d1 - $d2),0))
      return (
        max((-$d1, $d2)) + $h*$h*0.25 div $k,
        if ($d2 > -$d1) then tail($v2) else tail($v1)
      )
    },
    ($a, $b, $k)
  )
};

(:~
 : opRoundIntersect()
 : Rounded intersection of two shapes
 :)
declare function this:opRoundIntersect(
  $a as item(), 
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opRoundIntersect",
    function ($point as xs:double*) as xs:double* {
      let $v2 := $af($point)
      let $v1 := $bf($point)
      let $d1 := head($v1)
      let $d2 := head($v2)
      let $h := max(($k - abs($d1 - $d2), 0))
      return (
        max(($d1, $d2)) + $h*$h*0.25 div $k,
        if (head($v2) > head($v1)) then tail($v2) else tail($v1)
      )
    },
    ($a, $b, $k)
  )
};

(:~
 : opSmoothUnion()
 : Smoothed union of two shapes
 :)
declare function this:opSmoothUnion(
  $a as item(), 
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opSmoothUnion",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      let $h := util:clamp(0.5 + (0.5 * ($bd - $ad)) div $k, 0, 1)
      return (
        util:mix($bd, $ad, $h) - $k * $h * (1 - $h),
        if (head($av) < head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $k)
  )
};

(:~
 : opSmoothSubtract()
 : Smoothed difference of two shapes
 :)
declare function this:opSmoothSubtract(
  $a as item(),
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opSmoothSubtract",
    function ($point as xs:double*) as xs:double* {
      let $l := v:magnitude($point)
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      let $h :=  util:clamp(0.5 - (0.5 * ($ad + $bd)) div $k, 0, 1)
      return (
        util:mix($ad, -$bd, $h) + $k * $h * (1 - $h),
        if ($ad > -$bd) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $k)
  )
};

(:~
 : opSmoothIntersect()
 : Smoothed intersection of two shapes
 :)
declare function this:opSmoothIntersect(
  $a as item(), 
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opSmoothIntersect",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      let $h := util:clamp(0.5 - (0.5 * ($bd - $ad)) div $k, 0, 1)
      return (
        util:mix($bd, $ad, $h) + $k * $h * (1 - $h),
        if (head($av) > head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $k)
  )
};

(:~
 : opChamferUnion()
 : Chamfered union of two shapes
 :)
declare function this:opChamferUnion(
  $a as item(), 
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  let $sqrt_2_2 := math:sqrt(2) div 2
  return
  ann:f("sdf:opChamferUnion",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      return (
        min(($ad, $bd, ($ad - $k + $bd) * $sqrt_2_2)),
        if (head($av) > head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $k)
  )
};

(:~
 : opChamferIntersect()
 : Chamfered intersection of two shapes
 :)
declare function this:opChamferIntersect(
  $a as item(), 
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  let $sqrt_2_2 := math:sqrt(2) div 2
  return
  ann:f("sdf:opChamferIntersect",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      return (
        max(($ad, $bd, ($ad + $k + $bd) * $sqrt_2_2)),
        if (head($av) > head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $k)
  )
};

(:~
 : opChamferSubtract()
 : Chamfered difference of two shapes
 :)
declare function this:opChamferSubtract(
  $a as item(), 
  $b as item(),
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  let $sqrt_2_2 := math:sqrt(2) div 2
  return
  ann:f("sdf:opChamferSubtract",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      return (
        max(($ad, -$bd, ($ad + $k - $bd) * $sqrt_2_2)),
        if (head($av) > head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $k)
  )
};

declare %private function this:steps(
  $a as xs:double,
  $b as xs:double,
  $r as xs:double,
  $n as xs:double
) as xs:double
{
  let $s := $r div $n
  let $u := $b - $r
  return min(($a, $b, 0.5 * ($u + $a + abs(($u - $a + $s) mod (2 * $s)))))
};

(:~
 : opSteppedUnion()
 : Stepped union of two shapes
 :)
declare function this:opSteppedUnion(
  $a as item(), 
  $b as item(),
  $r as xs:double,
  $n as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opSteppedUnion",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      return (
        this:steps($ad, $bd, $r, $n),
        if (head($av) > head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $r, $n)
  )
};

(:~
 : opSteppedIntersect()
 : Stepped intersection of two shapes
 :)
declare function this:opSteppedIntersect(
  $a as item(), 
  $b as item(),
  $r as xs:double,
  $n as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opSteppedIntersection",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      return (
        -this:steps(-$ad, -$bd, $r, $n),
        if (head($av) > head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $r, $n)
  )
};

(:~
 : opSteppedSubtract()
 : Stepped difference of two shapes
 :)
declare function this:opSteppedSubtract(
  $a as item(), 
  $b as item(),
  $r as xs:double,
  $n as xs:double
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opSteppedSubtract",
    function ($point as xs:double*) as xs:double* {
      let $av := $af($point)
      let $bv := $bf($point)
      let $ad := head($av)
      let $bd := head($bv)
      return (
        -this:steps(-$ad, $bd, $r, $n),
        if (head($av) > head($bv)) then tail($av) else tail($bv)
      )
    },
    ($a, $b, $r, $n)
  )
};

(:~
 : opTwist()
 : 3D twisting
 :)
declare function this:opTwist(
  $base as item(), 
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opTwist",
    function ($point as xs:double*) as xs:double* {
      let $c := math:cos($k * v:py($point))
      let $s := math:sin($k * v:py($point))
      let $matrix := ($c, -$s, $s, $c)
      let $q := (
        $matrix[1] * v:px($point) + $matrix[2] * v:pz($point),
        $matrix[3] * v:px($point) + $matrix[4] * v:pz($point),
        v:py($point)
      )
      return $basef($q)
    },
    ($base, $k)
  )
};

(:~
 : opBend()
 : Simple bend of the base
 :)
declare function this:opBend(
  $base as item(), 
  $k as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opBend",
    function ($point as xs:double*) as xs:double* {
      let $c := math:cos($k * v:px($point))
      let $s := math:sin($k * v:px($point))
      let $matrix := ($c, -$s, $s, $c)
      let $q := (
        $matrix[1] * v:px($point) + $matrix[2] * v:py($point),
        $matrix[3] * v:px($point) + $matrix[4] * v:py($point),
        tail(tail($point))
      )
      return $basef($q)
    },
    ($base, $k)
  )
};

(:~
 : opAdd()
 : Add the two distance functions.
 : e.g. to add in a wobble
 :)
declare function this:opAdd(
  $a as item(),
  $b as item()
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opAdd",
    function ($point as xs:double*) as xs:double* {
      $af($point)=>v:add($bf($point))
    },
    ($a, $b)
  )
};

(:~
 : opAddInto()
 : Add the second distance function to the first, but keep any materials or
 : other secondary values from the first.
 : e.g. to add in a wobble but keep the same material
 :)
declare function this:opAddInto(
  $a as item(),
  $b as item()
) as map(*)
{
  if (not(ann:is-sdf($a))) then errors:error("ML-BADARGS", ("a", $a)) else (),
  if (not(ann:is-sdf($b))) then errors:error("ML-BADARGS", ("b", $b)) else (),
  let $af := f:function($a)
  let $bf := f:function($b)
  return
  ann:f("sdf:opAddInto",
    function ($point as xs:double*) as xs:double* {
      let $v1 := $af($point)
      let $v2 := $bf($point)
      return (head($v1) + head($v2), tail($v1))
    },
    ($a, $b)
  )
};

(:~ 
 : opInvert()
 : Flip insides with outsides
 :)
declare function this:opInvert(
  $base as item()
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opInvert",
    function ($point as xs:double*) as xs:double* {
      let $res := $basef($point)
      return (-head($res), tail($res))
    },
    $base
  )
};

(:~ 
 : opRescale()
 : Rescale the output distance.
 :)
declare function this:opRescale(
  $base as item(),
  $amplitude as xs:double
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  let $basef := f:function($base)
  return
  ann:f("sdf:opRescale",
    function ($point as xs:double*) as xs:double* {
      let $res := $basef($point)
      return ($amplitude * head($res), tail($res))
    },
    ($base, $amplitude)
  )
};

(:~
 : opModifyDistance()
 : Apply a gain function to the raw distance value.
 :)
declare function this:opModifyDistance(
  $base as item(),
  $f as item() 
) as map(*)
{
  if (not(ann:is-sdf($base))) then errors:error("ML-BADARGS", ("base", $base)) else (),
  if (not(ann:is-distance-function($f))) then errors:error("ML-BADARGS", ("f", $f)) else (),
  let $basef := f:function($base)
  let $ff := f:function($f)
  return
  ann:f("sdf:opModifyDistance",
    function ($point as xs:double*) as xs:double* {
      let $res := $basef($point)
      return ($ff(head($res)), tail($res))
    },
    ($base, $f)
  )
};

(: 
 : Distance modifiers: 
 : NOTE: Some of these operating unit interval so be sure distance is scaled first 
 :)
 
(:~ 
 : modAlmostIdentity()
 : Modify value by returning the result except close to 0.
 :
 : @param $m: lower bound below which we modify value
 : @param $n: lower bound of output value
 :)
declare function this:modAlmostIdentity(
  $m as xs:double,
  $n as xs:double
) as item()
{
  let $a := 2.0*$n - $m
  let $b := 2.0*$m - 3.0*$n
  return  
    ann:f("sdf:modAlmostIdentity", function ($x as xs:double) as xs:double {
      if (abs($x) > $m) then $x
      else (
        let $t := $x div $m
        return (($a * $t + $b)*$t*$t + $n)
      )
    })
};

(:~ 
 : modExpImpulse()
 : Modify value fast growth plus slow decay.
 : Maximum value becomes 1 at d=1/k
 :
 : @param $k: stretching
 :)
declare function this:modExpImpulse(
  $k as xs:double
) as item()
{
  ann:f("sdf:modExpImpulse", function ($x as xs:double) as xs:double {
    let $d := abs($x)
    let $sign := util:sign($x)
    return (
      let $h := $k * $d
      return ($sign * $h * math:exp(1.0 - $h), tail($d))
    )
  })
};

(:~ 
 : modQuadImpulse()
 : Modify value fast growth plus slow decay.
 : Maximum value becomes 1 at d=sqrt(1/k)
 :
 : @param $k: stretching
 :)
declare function this:modQuadImpulse(
  $k as xs:double
) as item()
{
  let $two-root-k := 2.0 * math:sqrt($k)
  return
    ann:f("sdf:modQuadImpulse", function ($x as xs:double) as xs:double {
      let $d := abs($x)
      let $sign := util:sign($x)
      return (
        $sign * $two-root-k * $d div (1.0 + $k*$d*$d)
      )
    })
};

(:~ 
 : modExpSustainedImpulse()
 : Modify value fast growth plus slow decay, with control over width of
 : both attack and release.
 :
 : @param $k: attack
 : @param $f: release
 :)
declare function this:modExpSustainedImpulse(
  $k as xs:double,
  $f as xs:double
) as item()
{
  ann:f("sdf:modExpSustainedImpulse", function ($x as xs:double) as xs:double {
    let $d := abs($x)
    let $sign := util:sign($x)
    let $s := max(($d - $f, 0.0))
    return (
      $sign * min(( ($d*$d) div ($f*$f), 1 + (2.0 div $f)*$s*math:exp(-$k*$s)))
    )
  })
};

(:~ 
 : modCubicPulse()
 : Modify value zero with a spread curve in the middle; pseudo-Gaussian
 :
 : @param $c: center value
 : @param $w: width
 :)
declare function this:modCubicPulse(
  $c as xs:double,
  $w as xs:double
) as item()
{
  ann:f("sdf:modCubicPulse", function ($x as xs:double) as xs:double {
    let $x := abs($x - $c)
    return (
      if ($x > $w) then 0.0
      else (
        let $x := $x div $w
        return (1.0 - $x*$x*(3.0 - 2.0*$x))
      )
    )
  })
};

(:~ 
 : modExpStep()
 : Modify value exponential decay
 :
 : @param $k: decay rate
 : @param $n: sharpness of decay; higher n => more like a step 
 :)
declare function this:modExpStep(
  $k as xs:double,
  $n as xs:double
) as item()
{
  ann:f("sdf:modExpStep", function ($x as xs:double) as xs:double {
    let $d := abs($x)
    let $sign := if ($x=0) then 1 else util:sign($x)
    return (
      $sign * math:exp( -$k * math:pow($d, $n) )
    )
  })
};

(:~ 
 : modGain()
 : Modify value expanding sides and compressing center
 : k=1 is identity, k>1 is ess curve, k<1 is symmetric inverse (k=a vs k=1/a)
 :
 : @param $k: gain 
 :)
declare function this:modGain(
  $k as xs:double
) as item()
{
  ann:f("sdf:modGain", function ($x as xs:double) as xs:double {
    let $d := abs($x)
    let $sign := util:sign($x)
    let $a := 0.5 * math:pow(2.0*(if ($d < 0.5) then $d else 1.0 - $d), $k)
    return (
      if ($d < 0.5) then $sign*$a else $sign*(1.0 - $a)
    )
  })
};

(:~ 
 : modPowerCurve()
 : Modify value map [0,1] to [0,1] with corners mapped to 0 and a curve
 : in between.
 :
 : @param $a: curve shape constant
 : @param $b: curve shape constant
 :)
declare function this:modPowerCurve(
  $a as xs:double,
  $b as xs:double
) as item()
{
  let $k := math:pow($a + $b, $a + $b) div (math:pow($a, $a) * math:pow($b, $b))
  return
    ann:f("sdf:modPowerCurve", function ($x as xs:double) as xs:double {
      let $d := abs($x)
      let $sign := util:sign($x)
      return (
        $sign * $k * math:pow($d, $a) * math:pow(1.0 - $d, $b)
      )
    })
};

(:~ 
 : modSinc()
 : Modify value map with a phase shifted sinc curve (sine cardinal) 
 : sinc(x) = 1 if x=0, sin(x)/x otherwise
 :
 : @param $k: controls degree of phase shift
 : @param $scale: rescales result
 :)
declare function this:modSinc(
  $k as xs:double,
  $scale as xs:double
) as item()
{
  ann:f("sdf:modSinc", function ($x as xs:double) as xs:double {
    let $d := abs($x)
    let $sign := util:sign($x)
    let $a := math:pi()*($k * $d - 1.0)
    return (
      if ($a = 0) then $sign * $scale
      else $sign * $scale * math:sin($a) div $a
    )
  })
};

declare function this:modSinc(
  $k as xs:double
) as item()
{
  this:modSinc($k, 1.0)
};

(:======================================================================
 : Rendering functions
 :======================================================================:)

(:~
 : renderTrace()
 : Trace passthrough for rendering functions
 : f and result are both function(xs:double*, xs:double*, xs:double*, item()) as xs:double* or a callable of same
 : final argument is function(xs:double*) as xs:double* or callable of same
 :)
declare function this:renderTrace(
  $f as item()
) as map(*)
{
  if (not(ann:is-renderer($f))) then errors:error("ML-BADARGS", ("f", $f)) else (),
  let $fn := f:function($f)
  return
  ann:f(util:function-name($f), 
    function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      $fn($tm, $origin, $direction, $sdf)=>trace(util:quote($f)||"("||util:quote($tm)||")")
    }
  )
};

(:~
 : renderPerformance()
 : Render the scene by mapping iterations into the given gradient.
 :)
declare function this:renderPerformance(
  $gradient as xs:string
) as map(*)
{
  let $colours := gradient:gradient-points($gradient)
  return 
    ann:f("sdf:renderPerformance", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      (: 64 = max iterations :)
      point:pcoordinates(
        cs:rgb-to-xyz(cs:interpolate-colour($tm[3] div 64, $colours))
      )
    })
};

(:~
 : renderDepth()
 : Render the scene by mapping depth to the given gradient.
 : @param $field defines the range for scaling
 :)
declare function this:renderDepth(
  $gradient as xs:string,
  $field as xs:double*
) as map(*)
{
  let $colours := gradient:gradient-points($gradient)
  let $n := count($colours)
  let $scale := $n div ($field[2] - $field[1])
  return 
    ann:f("sdf:renderDepth", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      let $adjusted-t := util:round(($tm[1] - $field[1]) * $scale) div $n
      return
        point:pcoordinates(
          cs:rgb-to-xyz(cs:interpolate-colour($adjusted-t, $colours))
        )
    })
};

(:~
 : renderColourPick()
 : Map from a material index to a colour. Colours should be RGB; function
 : maps to XYZ.
 :)
declare function this:renderColourPick(
  $pick-table as map(xs:string,item()*)*
) as map(*)
{
  let $n := count($pick-table)
  let $xyz-pick-table := $pick-table!cs:rgb-to-xyz(.)
  return
    ann:f("sdf:renderColourPick", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      point:pcoordinates($xyz-pick-table[util:modix(util:round($tm[2]), $n)])
    })
};

(:~
 : renderMaterial()
 : Map from a material index to a colouring function.
 :)
declare function this:renderMaterial(
  $pick-table as item()*
) as map(*)
{
  for $pick in $pick-table return if (not(ann:is-renderer($pick))) then errors:error("ML-BADARGS", ("pick-table", $pick)) else (),
  let $pick-tablef := for $pick in $pick-table return f:function($pick)
  let $n := count($pick-table)
  return
    ann:f("sdf:renderMaterial",
      function (
        $tm as xs:double*,
        $origin as xs:double*,
        $direction as xs:double*,
        $sdf as item()
      ) as xs:double*
      {
        let $f := $pick-tablef[util:modix(util:round($tm[2]), $n)]
        return $f($tm, $origin, $direction, $sdf)
      },
      $pick-table
    )
};

(:~
 : renderChoice()
 : Render everything deeper than (or equal to) a cutoff depth using one
 : renderer; everything nearer with a different one.
 :)
declare function this:renderChoice(
  $t-pivot as xs:double,
  $background as item(),
  $foreground as item()
) as map(*)
{
  if (not(ann:is-renderer($background))) then errors:error("ML-BADARGS", ("background", $background)) else (),
  if (not(ann:is-renderer($foreground))) then errors:error("ML-BADARGS", ("foreground", $foreground)) else (),
  let $foregroundf := f:function($foreground)
  let $backgroundf := f:function($background)
  return
  ann:f("sdf:renderChoice",
    function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      if (head($tm) <= $t-pivot)
      then $backgroundf($tm, $origin, $direction, $sdf)
      else $foregroundf($tm, $origin, $direction, $sdf)
    },
    ($t-pivot, $background, $foreground)
  )
};

(:~ 
 : renderConstant()
 : Return a constant colour. Colour should be given in RGB; function maps  
 : to XYZ.
 :)
declare function this:renderConstant(
  $colour as map(xs:string,item()*)
) as map(*)
{
  let $colourv := point:pcoordinates(cs:rgb-to-xyz($colour))
  return
    ann:f("sdf:renderConstant", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      $colourv
    })
};

(:~ 
 : renderRandom()
 : Return a constant colour. Colour should be given in RGB; function maps  
 : to XYZ.
 :)
declare function this:renderRandom(
  $colours as map(xs:string,item()*)*
) as map(*)
{
  let $xyz-colours := $colours!cs:rgb-to-xyz(.)
  return
    ann:f("sdf:renderRandom", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      rand:select-random($xyz-colours)!point:pcoordinates(.)
    })
};

(:~
 : renderGradient()
 : Apply a shading gradient based on t
 : Base colour should be given in RGB, function maps to XYZ.
 :)
declare function this:renderGradient(
  $base-colour as map(xs:string,item()*),
  $strength as xs:double
) as map(*)
{
  let $colourv := point:pcoordinates(cs:rgb-to-xyz($base-colour))
  return
    ann:f("sdf:renderGradient", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      let $shade := head($tm) * $strength
      return $colourv=>v:sub(($shade, $shade, $shade))
    })
};

(:~
 : renderXGradient()
 : Apply a shading gradient based on x coordinate
 : Base colour should be given in RGB, function maps to XYZ.
 :)
declare function this:renderXGradient(
  $base-colour as map(xs:string,item()*),
  $strength as xs:double
) as map(*)
{
  let $colourv := point:pcoordinates(cs:rgb-to-xyz($base-colour))
  return
    ann:f("sdf:renderXGradient", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      let $shade := v:px($direction) * $strength
      return $colourv=>v:sub(($shade, $shade, $shade))
    })
};

(:~
 : renderYGradient()
 : Apply a shading gradient based on y coordinate
 : Base colour should be given in RGB, function maps to XYZ.
 :)
declare function this:renderYGradient(
  $base-colour as map(xs:string,item()*),
  $strength as xs:double
) as map(*)
{
  let $colourv := point:pcoordinates(cs:rgb-to-xyz($base-colour))
  return
    ann:f("sdf:renderYGradient", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      let $shade := v:py($direction) * $strength
      return $colourv=>v:sub(($shade, $shade, $shade))
    })
};


(:~
 : renderZGradient()
 : Apply a shading gradient based on y coordinate
 : Base colour should be given in RGB, function maps to XYZ.
 :)
declare function this:renderZGradient(
  $base-colour as map(xs:string,item()*),
  $strength as xs:double
) as map(*)
{
  let $colourv := point:pcoordinates(cs:rgb-to-xyz($base-colour))
  return
    ann:f("sdf:renderZGradient", function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      let $shade := v:pz($direction) * $strength
      return $colourv=>v:sub(($shade, $shade, $shade))
    })
};

(:~
 : opRenderFuzz()
 : Fuzz out the depth of field
 :)
declare function this:opRenderFuzz(
  $base-colour as item(),
  $field as xs:double*
) as map(*)
{
  if (not(ann:is-renderer($base-colour))) then errors:error("ML-BADARGS", ("base-colour", $base-colour)) else (),
  let $base-colourf := f:function($base-colour)
  return
  ann:f("sdf:opRenderFuzz",
    function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      let $col := $base-colourf($tm, $origin, $direction, $sdf)
      let $smooth := util:smoothstep($field[1], $field[2], head($tm))
      return $col=>v:times(1.0 - $smooth)  
    },
    ($base-colour, $field)
  )
};

(:~
 : opRenderDiffuseLight()
 : Diffuse lighting.
 :
 : Lambert reflection: Id = L·N C Il 
 : Id = surface brightness, L = direction of light, N = surface normal
 : C = colour, Il = incoming light intensity
 :)
declare function this:opRenderDiffuseLight(
  $surface-colour as item(),
  $light as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-renderer($surface-colour))) then errors:error("ML-BADARGS", ("surface-colour", $surface-colour)) else (),
  let $surface-colourf := f:function($surface-colour)
  let $light-dir := light:direction($light)
  let $light-colour := light:directional($light)
  let $intensity := light:intensity($light)
  let $ambient-intensity := light:ambient-intensity($light)
  let $ambient := light:ambient($light)=>v:times($ambient-intensity)
  return
    ann:f("sdf:opRenderDiffuseLight",
      function (
        $tm as xs:double*,
        $origin as xs:double*,
        $direction as xs:double*,
        $sdf as item()
      ) as xs:double*
      {
        let $surface := $surface-colourf($tm, $origin, $direction, $sdf)
        let $pos := this:linear($origin, $direction, head($tm))
        let $n := this:calc-normal($pos, $sdf)
        let $surface-brightness := max((v:dot($n, $light-dir), 0.0)) * $intensity
        let $directional := $light-colour=>v:times($surface-brightness)
        return
          $surface=>v:multiply( $directional=>v:add($ambient) )
      },
      ($surface-colour, $light)
    )
};

(:~
 : opRenderHardShadow()
 : Modify colours to provide hard shadows.
 :)
declare function this:opRenderHardShadow(
  $base-colour as item(),
  $light as map(xs:string,item()*),
  $field as xs:double*,
  $shadow-scale as xs:double,
  $shadow-strength as xs:double
) as map(*)
{
  if (not(ann:is-renderer($base-colour))) then errors:error("ML-BADARGS", ("base-colour", $base-colour)) else (),
  let $base-colourf := f:function($base-colour)
  let $light-dir := light:direction($light)
  return
    ann:f("sdf:opRenderHardShadow",
      function (
        $tm as xs:double*,
        $origin as xs:double*,
        $direction as xs:double*,
        $sdf as item()
      ) as xs:double*
      {
        let $col := $base-colourf($tm, $origin, $direction, $sdf)
        let $pos := this:linear($origin, $direction, head($tm))
        let $normal := this:calc-normal($pos, $sdf)
        let $shadow-origin := $pos=>v:add($normal=>v:times($shadow-scale))
        let $st := this:cast-ray($shadow-origin, $light-dir, $field, $sdf)[2]
        let $shadow := if ($st != -1.0) then 1.0 else 0.0
        return v:mix($col, $col=>v:times($shadow-strength), $shadow)
      },
      ($base-colour, $light, $field, $shadow-scale, $shadow-strength)
    )
};

(:~
 : opRenderSoftShadow()
 : Modify colours to provide soft shadows.
 :)
declare function this:opRenderSoftShadow(
  $base-colour as item(),
  $light as map(xs:string,item()*),
  $field as xs:double*,
  $shadow-scale as xs:double,
  $shadow-strength as xs:double,
  $n-shadows as xs:integer,
  $shadow-falloff as xs:double
) as map(*)
{
  if (not(ann:is-renderer($base-colour))) then errors:error("ML-BADARGS", ("base-colour", $base-colour)) else (),
  let $base-colourf := f:function($base-colour)
  let $light-dir := light:direction($light)
  let $rand := dist:uniform(-1.0, 1.0 - $config:ε)
  return
    ann:f("sdf:opRenderSoftShadow",
      function (
        $tm as xs:double*,
        $origin as xs:double*,
        $direction as xs:double*,
        $sdf as item()
      ) as xs:double*
      {
        let $col := $base-colourf($tm, $origin, $direction, $sdf)
        let $pos := this:linear($origin, $direction, head($tm))
        let $normal := this:calc-normal($pos, $sdf)
        let $shadow-origin := $pos=>v:add($normal=>v:times($shadow-scale))
        let $shadow :=
          fold-left(
            1 to $n-shadows,
            0.0,
            function ($shadow as xs:double, $s as xs:integer) as xs:double {
              let $fall := $shadow-falloff * rand:randomize($rand)
              let $shadow-dir := $light-dir=>v:add(($fall, $fall, $fall))
              let $st := this:cast-ray($shadow-origin, $shadow-dir, $field, $sdf)[2]
              return (
                if ($st != -1.0) then $shadow + 1.0 else $shadow
              )
            }
          )
        return v:mix($col, $col=>v:times($shadow-strength), $shadow div $n-shadows)
      },
      ($base-colour, $light, $field, $shadow-scale, $shadow-strength, $n-shadows, $shadow-falloff)
    )
};

(:~
 : opRenderComplexLighting()
 : Modify colours to provide soft shadows in diffuse light.
 :)
declare function this:opRenderComplexLighting(
  $base-colour as item(),
  $light as map(xs:string,item()*), (: light source :)
  $shadow-scale as xs:double, (: 0.01 :)
  $k as xs:double (: softness 32.0 :)
) as map(*)
{
  if (not(ann:is-renderer($base-colour))) then errors:error("ML-BADARGS", ("base-colour", $base-colour)) else (),
  let $base-colourf := f:function($base-colour)
  let $light-direction := light:direction($light)
  let $light-colour :=
    light:directional($light)=>
      v:times( light:intensity($light) )=>
      v:add( light:ambient($light)=>v:times( light:ambient-intensity($light) ) )
  return
    ann:f("sdf:opRenderComplexLighting",
      function (
        $tm as xs:double*,
        $origin as xs:double*,
        $direction as xs:double*,
        $sdf as item()
      ) as xs:double*
      {
        let $pos := this:linear($origin, $direction, head($tm))
        let $col := $base-colourf($tm, $origin, $direction, $sdf)
        let $normal := this:calc-normal($pos, $sdf)
        let $shade-origin := $direction=>v:add($normal=>v:times($shadow-scale))
        let $shade := head(
          util:while(
            function ($res as xs:double, $t as xs:double, $i as xs:integer, $working as xs:boolean) as xs:boolean { $i <= 16 and $working },
            function ($res as xs:double, $t as xs:double, $i as xs:integer, $working as xs:boolean) as item()* {
              let $h := head($sdf( this:linear($shade-origin, $light-direction, $t)))
              let $res := min(($res, $k*$h div $t))
              let $t := $t + util:clamp($h, 0.01, 1.0) (: fineness of calc :)
              return (
                $res, $t, $i + 1, not($h < 0.0001)
              )
            },
            1.0, 0.0, 1, true()
          )
        )
        let $shade :=
          util:clamp($shade, 0.0, 1.0) *
          util:clamp(v:dot($normal, $light-direction), 0.0, 1.0)
        let $occlusion := head(
          util:while(
            function ($res as xs:double, $sca as xs:double, $i as xs:integer) as xs:boolean { $i < 5 },
            function ($res as xs:double, $sca as xs:double, $i as xs:integer) as item()* {
              let $hr := 0.02 + 0.025 * xs:double($i) * xs:double($i)
              let $aopos := $pos=>v:add($normal=>v:times($hr))
              let $dd := head($sdf($aopos))
              return ($res + -($dd - $hr)*$sca, $sca * 0.995, $i + 1)
            },
            0.0, 1.0, 0
          )
        )
        let $occlusion :=
          (1.0 - util:clamp($occlusion, 0.0, 1.0))*util:clamp(v:py($normal), 0.0, 1.0)
        return
          ($col=>v:times($occlusion))=>             (: base + ambient shadows :)
            v:add($light-colour=>v:times($shade))=> (: lighting :)
            v:times(math:exp(-0.2*head($tm)))       (: softening :)
      },
      ($base-colour, $light, $shadow-scale, $k)
    )
};

(:~
 : opRenderColourWaves()
 : Waves of colour bands
 :)
declare function this:opRenderColourWaves(
  $base-colour as item(),
  $darkness as xs:double,  (: -3.0; more negative=lighter :)
  $frequency as xs:double, (: 75; higher=more bands :)
  $strength as xs:double (: 0.2 [0,1]; lower=gentler bands :)
) as map(*)
{
  if (not(ann:is-renderer($base-colour))) then errors:error("ML-BADARGS", ("base-colour", $base-colour)) else (),
  let $base-colourf := f:function($base-colour)
  return
  ann:f("sdf:opRenderColourWaves",
    function (
      $tm as xs:double*,
      $origin as xs:double*,
      $direction as xs:double*,
      $sdf as item()
    ) as xs:double*
    {
      let $d := head($tm)
      return
        $base-colourf($tm, $origin, $direction, $sdf)=>
          v:times(1.0 - math:exp($darkness*abs($d)))=>
          v:times((1 - $strength) + $strength*math:cos($frequency*$d))
    },
    ($base-colour, $darkness, $frequency, $strength)
  )
};

(:~
 : opRenderColourBlend()
 : Blend a colour into the source colouring based on the distance
 :)
declare function this:opRenderColourBlend(
  $base-colour as item(),
  $blend-colour as map(xs:string,item()*)
) as map(*)
{
  if (not(ann:is-renderer($base-colour))) then errors:error("ML-BADARGS", ("base-colour", $base-colour)) else (),
  let $base-colourf := f:function($base-colour)
  let $blend := point:pcoordinates($blend-colour)
  return
    ann:f("sdf:opRenderColourBlend",
      function (
        $tm as xs:double*,
        $origin as xs:double*,
        $direction as xs:double*,
        $sdf as item()
      ) as xs:double*
      {
        let $col := $base-colourf($tm, $origin, $direction, $sdf)
        return v:mix($col, $blend, $tm[1])
      },
      ($base-colour, $blend-colour)
    )
};

(:~
 : opRenderModifyColour()
 : Blend a colour into the source colouring based on a function of the distance
 :)
declare function this:opRenderModifyColour(
  $base-colour as item(),
  $blend-colour as map(xs:string,item()*),
  $f as item()
) as map(*)
{
  if (not(ann:is-renderer($base-colour))) then errors:error("ML-BADARGS", ("base-colour", $base-colour)) else (),
  if (not(ann:is-distance-function($f))) then errors:error("ML-BADARGS", ("f", $f)) else (),
  let $base-colourf := f:function($base-colour)
  let $ff := f:function($f)
  let $blend := point:pcoordinates($blend-colour)
  return
    ann:f("sdf:opRenderModifyColour",
      function (
        $tm as xs:double*,
        $origin as xs:double*,
        $direction as xs:double*,
        $sdf as item()
      ) as xs:double*
      {
        let $col := $base-colourf($tm, $origin, $direction, $sdf)
        let $mod := $ff($tm[1])
        return v:mix($col, $blend, $mod)
      },
      ($base-colour, $blend-colour, $f)
    )
};

(:======================================================================
 : Ray casting main drivers
 :======================================================================:)

(:~
 : camera()
 : Define the camera view: its position and target
 :)
declare function this:camera(
  $position as map(xs:string,item()*),
  $target as map(xs:string,item()*),
  $perspective as xs:double (: Z scaling :)
) as map(xs:string,item()*)
{
  let $forward := point:normalize(point:sub($target, $position))
  let $right := point:normalize(point:cross(point:point(0,1,0), $forward))
  let $up := point:normalize(point:cross($forward, $right))
  return
    map {
      "kind": "camera",
      "position": point:pcoordinates($position),
      "view": point:pcoordinates($position),
      "target": point:pcoordinates($target),
      "perspective": $perspective,
      "forward": point:pcoordinates($forward)=>v:times($perspective),
      "right": point:pcoordinates($right),
      "up": point:pcoordinates($up)
    }
};

(:~
 : camera()
 : Define camera POV based on view angles
 :)
declare function this:camera(
  $angle as xs:double,
  $view as map(xs:string,item()*),
  $target as map(xs:string,item()*),
  $perspective as xs:double (: Z scaling :)
) as map(xs:string,item()*)
{
  let $radians := util:radians($angle)
  let $origin :=
    point:pcoordinates($view)=>
      v:multiply((math:cos($radians), 1.0, math:sin($radians)))
  let $forward := v:normalize(point:pcoordinates($target)=>v:sub($origin))
  let $right := v:normalize(v:cross($forward,(0,1,0)))
  let $up := v:normalize(v:cross($forward, $right))
  return
    map {
      "kind": "camera",
      "position": $origin,
      "angle": $angle,
      "view": point:pcoordinates($view),
      "target": point:pcoordinates($target),
      "perspective": $perspective,
      "forward": $forward=>v:times($perspective),
      "right": $right,
      "up": $up
    }
};

(:~
 : castDepth()
 : A caster that calculates the depth relative to a particular ray
 :)
declare function this:castDepth(
  $field as xs:double*
) as map(*)
{
  ann:f("sdf:castDepth", function ($origin as xs:double*, $direction as xs:double*, $sdf as item()) as xs:double* {
    this:cast-ray($origin, $direction, $field, $sdf)
  })
};

(:~
 : castDirect()
 : A caster that just returns the SDF unaltered
 :)
declare function this:castDirect(
) as map(*)
{
  ann:f("sdf:castDirect", function ($origin as xs:double*, $direction as xs:double*, $sdf as item()) as xs:double* {
    $sdf($direction)
  })
};

(:~ 
 : render()
 : Compute rendering for a normalized point (scaled to canvas)
 :)
declare function this:render(
  $uv as map(xs:string,item()*), (: normalized point :)
  $camera as map(xs:string,item()*)?,
  $caster as item(),
  $sdf as item(),
  $render as item()
) as map(xs:string,item()*) (: rgba :)
{
  (:
   : This is called a lot: you're better off with make-point-render() 
   :)
  if (not(ann:is-caster($caster))) then errors:error("ML-BADARGS", ("caster", $caster)) else (),
  if (not(ann:is-sdf($sdf))) then errors:error("ML-BADARGS", ("sdf", $sdf)) else (),
  if (not(ann:is-renderer($render))) then errors:error("ML-BADARGS", ("render", $render)) else (),
  let $casterf := f:function($caster)
  let $sdff := f:function($sdf)
  let $renderf := f:function($render)
  return (
    if (empty($camera)) then (
      (: Take each point verbatim; really only useful for 2D :)
      let $ray-origin := point:pcoordinates($point:ORIGIN3D)
      let $ray-direction := point:pcoordinates($uv)
      let $tm := $casterf($ray-origin, $ray-direction, $sdff)
      let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
      return cs:xyz-to-rgb(point:vector($col))
    ) else (
      let $ray-origin := $camera("position")
      let $ray-direction := this:camera-ray-direction($uv, $camera)
      let $tm := $casterf($ray-origin, $ray-direction, $sdff)
      let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
      return cs:xyz-to-rgb(point:vector($col))
    )
  )
};

declare function this:make-point-render(
  $camera as map(xs:string,item()*)?,
  $caster as item(),
  $sdf as item(),
  $render as item()
) as function(map(xs:string,item()*)) as map(xs:string,item()*)
{
  if (not(ann:is-caster($caster))) then errors:error("ML-BADARGS", ("caster", $caster)) else (),
  if (not(ann:is-sdf($sdf))) then errors:error("ML-BADARGS", ("sdf", $sdf)) else (),
  if (not(ann:is-renderer($render))) then errors:error("ML-BADARGS", ("render", $render)) else (),
  let $casterf := f:function($caster)
  let $sdff := f:function($sdf)
  let $renderf := f:function($render)
  return (
    if (empty($camera)) then (
      let $ray-origin := point:pcoordinates($point:ORIGIN3D) return
      function ($uv as map(xs:string,item()*)) as map(xs:string,item()*) {
        let $ray-direction := point:pcoordinates($uv)
        let $tm := $casterf($ray-origin, $ray-direction, $sdff)
        let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
        return cs:xyz-to-rgb(point:vector($col))
      }
    ) else (
      let $ray-origin := $camera("position") return
      function ($uv as map(xs:string,item()*)) as map(xs:string,item()*) {
        let $ray-direction := this:camera-ray-direction($uv, $camera)
        let $tm := $casterf($ray-origin, $ray-direction, $sdff)
        let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
        return cs:xyz-to-rgb(point:vector($col))
      }
    )
  )
};


declare function this:renderv(
  $uv as xs:double*, (: normalized point :)
  $camera as map(xs:string,item()*)?,
  $caster as item(),
  $sdf as item(),
  $render as item()
) as map(xs:string,item()*) (: rgba :)
{
  (:
   : This is called a lot: you're better off with make-vector-render() 
   :)
  if (not(ann:is-caster($caster))) then errors:error("ML-BADARGS", ("caster", $caster)) else (),
  if (not(ann:is-sdf($sdf))) then errors:error("ML-BADARGS", ("sdf", $sdf)) else (),
  if (not(ann:is-renderer($render))) then errors:error("ML-BADARGS", ("render", $render)) else (),
  let $casterf := f:function($caster)
  let $sdff := f:function($sdf)
  let $renderf := f:function($render)
  return (
    if (empty($camera)) then (
      (: Take each point verbatim; really only useful for 2D :)
      let $ray-origin := point:pcoordinates($point:ORIGIN3D)
      let $ray-direction := $uv
      let $tm := $casterf($ray-origin, $ray-direction, $sdff)
      let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
      return cs:xyz-to-rgb(point:vector($col))
    ) else (
      let $ray-origin := $camera("position")
      let $ray-direction := this:camera-ray-directionv($uv, $camera)
      let $tm := $casterf($ray-origin, $ray-direction, $sdff)
      let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
      return cs:xyz-to-rgb(point:vector($col))
    )
  )
};

declare function this:make-vector-render(
  $camera as map(xs:string,item()*)?,
  $caster as item(),
  $sdf as item(),
  $render as item()
) as function(xs:double*) as map(xs:string,item()*)
{
  if (not(ann:is-caster($caster))) then errors:error("ML-BADARGS", ("caster", $caster)) else (),
  if (not(ann:is-sdf($sdf))) then errors:error("ML-BADARGS", ("sdf", $sdf)) else (),
  if (not(ann:is-renderer($render))) then errors:error("ML-BADARGS", ("render", $render)) else (),
  let $casterf := f:function($caster)
  let $sdff := f:function($sdf)
  let $renderf := f:function($render)
  return (
    if (empty($camera)) then (
      let $ray-origin := point:pcoordinates($point:ORIGIN3D) return
      function ($uv as xs:double*) as map(xs:string,item()*) {
        let $ray-direction := $uv
        let $tm := $casterf($ray-origin, $ray-direction, $sdff)
        let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
        return cs:xyz-to-rgb(point:vector($col))
      }
    ) else (
      let $ray-origin := $camera("position") return
      function ($uv as xs:double*) as map(xs:string,item()*) {
        let $ray-direction := this:camera-ray-directionv($uv, $camera)
        let $tm := $casterf($ray-origin, $ray-direction, $sdff)
        let $col := $renderf($tm, $ray-origin, $ray-direction, $sdff)
        return cs:xyz-to-rgb(point:vector($col))
      }
    )
  )
};

declare function this:camera-ray-direction(
  $uv as map(xs:string,item()*),
  $camera as map(xs:string,item()*)
) as xs:double*
{
  v:normalize(
    (
      ($camera("right"))=>v:times(point:px($uv))
    )=>v:add(
      ($camera("up"))=>v:times(point:py($uv))
    )=>v:add(
      $camera("forward")
    )
  )
};

declare function this:camera-ray-directionv(
  $uv as xs:double*,
  $camera as map(xs:string,item()*)
) as xs:double*
{
  v:normalize(
    (
      ($camera("right"))=>v:times(v:px($uv))
    )=>v:add(
      ($camera("up"))=>v:times(v:py($uv))
    )=>v:add(
      $camera("forward")
    )
  )
};

(:~
 : cast-ray()
 : Main ray-casting function
 :)
declare function this:cast-ray(
  $ray-origin as xs:double*,
  $ray-direction as xs:double*,
  $field as xs:double*, (: near field, far field :)
  $sdf as item()
) as xs:double* (: distance, material, iterations :)
{
  util:assert(count($ray-origin)=3, "bad ray origin"),
  (: util:assert(count($ray-direction)=3, "bad ray direction"), :)
  let $ray-direction := if (count($ray-direction)=2) then ($ray-direction,0) else $ray-direction
  let $ε := 0.0001
  let $start := $field[1]
  let $end := $field[2]
  let $max-iterations := 64
  let $data :=
    util:while(
      function ($working as xs:boolean, $t as xs:double, $m as xs:double, $i as xs:integer) as xs:boolean {
        $working and $i <= $max-iterations
      },
      function ($working as xs:boolean, $t as xs:double, $m as xs:double, $i as xs:integer) {
        let $pos := this:linear($ray-origin, $ray-direction, $t)
        let $result := $sdf($pos)
        let $distance := $result[1]
        let $material := ($result[2],0)[1]
        return (
          if ($distance < $ε) then (
            false(), $t, $m, $i + 1
          ) else if ($t + $distance > $end) then (
            false(), -1.0, -1.0, $i + 1
          ) else (
            true(), $t + $distance, $material, $i + 1
          )
        )
      },
      true(), $start, -1.0, 1
    )
  let $working := head($data)
  return
    if ($working) then (-1.0, -1.0, $data[4]) else tail($data)
};

declare function this:calc-normal(
  $point as xs:double*,
  $sdf as item()
) as xs:double*
{
  let $c := head($sdf($point))
  let $epsilon := 0.001
  (: Compare sdf in offsets from point :)
  return
    v:sub(
      v:normalize((
        head($sdf($point=>v:add(($epsilon, 0, 0)))),
        head($sdf($point=>v:add((0, $epsilon, 0)))),
        head($sdf($point=>v:add((0, 0, $epsilon))))
      )),
      ($c, $c, $c)
    )
};

declare function this:linear(
  $origin as xs:double*,
  $direction as xs:double*,
  $t as xs:double
) as xs:double*
{
  util:zip(
    function ($a as xs:double, $b as xs:double) as xs:double {$a + $t * $b},
      $origin,
      $direction
  )
};

(:~ 
 : normalize-to-canvas()
 : Map point from canvas space to unit box. Need to shift and scale.
 : Canvas: [0:w, 0:h] Unit box: [-1:1, -1:1]
 :)
declare function this:normalize-to-canvas(
  $pt as map(xs:string,item()*),
  $canvas as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $width := box:width($canvas)
  let $height := box:height($canvas)
  return (
    point:vector((
      2.0 * (point:px($pt) - ($width div 2)) div $width,
      2.0 * (point:py($pt) - ($height div 2)) div $height,
      point:pcoordinates($pt)[position() > 2]
    ))
  )
};

declare function this:normalize-to-canvasv(
  $pt as xs:double*,
  $canvas as map(xs:string,item()*)
) as xs:double*
{
  let $width := box:width($canvas)
  let $height := box:height($canvas)
  return (
    (
      2.0 * (v:px($pt) - ($width div 2)) div $width,
      2.0 * (v:py($pt) - ($height div 2)) div $height,
      ($pt)[position() > 2]
    )
  )
};

(:~ 
 : normalize-to-canvas()
 : Map point from canvas space to unit box. Need to shift and scale.
 : Canvas: [0:w, 0:h] Unit box: [-1:1, -1:1]
 : A version defined for efficiency when we are scanning through a whole
 : matrix of x/y values.
 :)
declare function this:normalize-to-canvas(
  $x as xs:double,
  $y as xs:double,
  $width as xs:double,
  $height as xs:double,
  $aspect as xs:double
) as map(xs:string,item()*)
{
  point:point(
    2.0 * ($x - ($width div 2)) div $width,
    2.0 * ($y - ($height div 2)) div ($aspect * $height)
  )
};

declare function this:normalize-to-canvasv(
  $x as xs:double,
  $y as xs:double,
  $width as xs:double,
  $height as xs:double,
  $aspect as xs:double
) as xs:double*
{
  (
    2.0 * ($x - ($width div 2)) div $width,
    2.0 * ($y - ($height div 2)) div ($aspect * $height)
  )
};

(:~
 : scale-from-canvas()
 : Scale a distance metric from canvas space to distance space. Scaling only.
 : (Scales by width.)
 : Canvas distance: [-w/2:w/2] Unit distance: [-1:1] 
 :)
declare function this:scale-from-canvas(
  $d as xs:double,
  $canvas as map(xs:string,item()*)
) as xs:double
{
  if (empty($canvas)) then $d else
  let $width := $canvas("width")
  return (
    2.0 * $d div $width
  )
};

(:~
 : scale-from-unit()
 : Scale a distance metric from distance space to canvas space. Scaling only.
 : (Scales by width.)
 : Unit distance: [-1:1] Canvas distance: [-w/2:w/2] 
 :)
declare function this:scale-from-unit(
  $d as xs:double,
  $canvas as map(xs:string,item()*)
) as xs:double
{
  let $width := $canvas("width")
  return (
    0.5 * $d * $width
  )
};

declare %private function this:regionToSDF(
  $region as map(xs:string,item()*),
  $op as xs:string,
  $canvas as map(xs:string,item()*)?,
  $options as map(xs:string,item()*)
) as map(*)
{
  switch (util:kind($region))
  case "point" return (
    let $size := this:scale-from-canvas($options("radius"),$canvas)
    let $pt := this:normalize-to-canvas($region, $canvas)=>point:minus()
    return this:opTranslate(sd2:sdCircle($size), $pt)
  )
  case "box" return (
    if (point:same(box:min-point($region), box:max-point($region)))
    then this:regionToSDF(box:min-point($region), $op, $canvas, $options)
    else (
      let $origin := this:normalize-to-canvas($point:ORIGIN, $canvas)
      let $min := this:normalize-to-canvas(box:min-point($region), $canvas)
      let $max := this:normalize-to-canvas(box:max-point($region), $canvas)
      let $sides := point:sub($max, $min)
      let $offset := $sides=>point:times(0.5)=>point:add($min)=>point:minus()=>point:add($origin)
      return
        this:opTranslate(
          sd2:sdBox(point:px($sides), point:py($sides)),
          $offset
        )
    )
  )
  case "block" return (
    let $origin := this:normalize-to-canvas($point:ORIGIN, $canvas)
    let $center := this:normalize-to-canvas(box:center($region), $canvas)
    let $offset := $center=>point:minus()=>point:add($origin)
    return (
      this:opTranslate(
        sd3:sdBox(box:width($region), box:height($region), box:depth($region)),
        $offset
      )
    )
  )
  case "space" return (
    if (point:same(box:min-point($region), box:max-point($region)))
    then this:regionToSDF(box:min-point($region), $op, $canvas, $options)
    else (
      let $origin := this:normalize-to-canvas($point:ORIGIN, $canvas)
      let $min := this:normalize-to-canvas(box:min-point($region), $canvas)
      let $max := this:normalize-to-canvas(box:max-point($region), $canvas)
      let $sides := point:sub($max, $min)
      let $offset := $sides=>point:times(0.5)=>point:add($min)=>point:minus()=>point:add($origin)
      return
        this:opTranslate(
          sd2:sdBox(point:px($sides), point:py($sides)),
          $offset
        )
    )
  )
  case "ellipse" return (
    if (ellipse:is-circle($region)) then (
      let $center := this:normalize-to-canvas(ellipse:center($region), $canvas)
      let $r := this:scale-from-canvas(ellipse:radius($region), $canvas)
      return
        this:opTranslate(
          sd2:sdCircle($r),
          $center=>point:minus()
        )
    ) else (
      let $center := this:normalize-to-canvas(ellipse:center($region), $canvas)
      let $a := this:scale-from-canvas(ellipse:rx($region), $canvas)
      let $b := this:scale-from-canvas(ellipse:ry($region), $canvas)
      let $rotate := ellipse:rotation($region)
      return
        this:opTranslate(
          if ($rotate = 0) 
          then sd2:sdEllipse($a, $b)
          else this:opRotate(sd2:sdEllipse($a, $b), $rotate),
          $center=>point:minus()
        ) 
    )
  )
  case "edge" return (
    if (point:same(edge:start($region), edge:end($region)))
    then this:regionToSDF(edge:start($region), $op, $canvas, $options)
    else (
      let $start := this:normalize-to-canvas(edge:start($region), $canvas)
      let $end := this:normalize-to-canvas(edge:end($region), $canvas)
      return sd2:sdSegment($start, $end)
    )
  )
  case "quad" return (
    if (point:same(edge:start($region), edge:end($region)))
    then this:regionToSDF(edge:start($region), $op, $canvas, $options)
    else (
      let $start := this:normalize-to-canvas(edge:start($region), $canvas)
      let $end := this:normalize-to-canvas(edge:end($region), $canvas)
      let $control := this:normalize-to-canvas(edge:controls($region)[1], $canvas)
      return
        if ($options("ud"))
        then sd2:udQuad(edge:quad($start, $end, $control))
        else sd2:sdQuad(edge:quad($start, $end, $control))
    )
  )
  case "cubic" return ( (: Approximate with interpolation :)
    if (point:same(edge:start($region), edge:end($region)))
    then this:regionToSDF(edge:start($region), $op, $canvas, $options)
    else (
      let $interpolate := $options("interpolate")
      let $n := if ($interpolate > 0) then $interpolate else 4
      let $edges := edge:to-edges(edge:interpolate($n, $region))
      return
        if ($op=("union","smooth-union","round-union","stepped-union","chamfer-union"))
        then this:toSDF($edges, $canvas, $options)
        else (
          this:toSDF($edges, $canvas,
            util:merge-into($options, map {"op":"union"})
        )
      )
    )
  )
  case "arc" return (
    if (point:same(edge:start($region), edge:end($region)))
    then this:regionToSDF(edge:start($region), $op, $canvas, $options)
    else (
      let $arc-width := this:scale-from-canvas($options("width"), $canvas)
      let $arc := edge:as-angle-arc($region)
      let $angles := ($arc=>edge:arc-angles())!util:remap-degrees(.)
      let $large := edge:arc-large($arc)
      let $flipped := edge:arc-flipped($arc)
      let $delta := $angles[2] - $angles[1]
      let $swap :=
        ($large and
          ($delta < -180 or ($delta > 0 and $delta < 180))
        ) or
        (not($large) and
          (($delta < 0 and abs($delta) < 180) or ($delta > 180))
        )
      let $adjust-from :=
        ($large and $delta > 0 and $delta < 180) or
        (not($large) and $delta > 180)
      let $adjust-to :=
        ($large and $delta < 0 and abs($delta) <= 180) or
        (not($large) and $delta <= -180)
      let $adjusted-from :=
        if ($adjust-from) then $angles[1] + 360 else $angles[1]
      let $adjusted-to :=
        if ($adjust-to) then $angles[2] + 360 else $angles[2]
      let $from-angle :=
        if ($swap) then $adjusted-to else $adjusted-from
      let $to-angle := 
        if ($swap) then $adjusted-from else $adjusted-to
      let $center := this:normalize-to-canvas($arc=>edge:arc-center(), $canvas)
      let $radius := this:scale-from-canvas($arc=>edge:arc-radius(), $canvas)
      (: 
       : If we run from A to B, orientation is 180 from angle in 
       : the middle of the range
       : orientation = 180 + (A + B)/2 
       : aperture = (B - A)
       : However, orientation of degrees is flipped, so we need to negate orientation
       :)
      let $orientation := util:remap-degrees(-(180 + ($to-angle + $from-angle) div 2.0))
      let $aperture := util:remap-degrees(($to-angle - $from-angle))
      (: center and radius already normalized :)
      let $from := point:destination($center, $from-angle, $radius)
      let $to := point:destination($center, $to-angle, $radius)
      return (
        this:opTranslate(
          if ( ($swap and not($flipped)) or (not($swap) and $flipped) ) then (
            this:opReflectXY(
              sd2:sdArc($orientation, $aperture, $radius, $arc-width),
              $from,
              $to
            )
          ) else (
            sd2:sdArc($orientation, $aperture, $radius, $arc-width)
          ),
          $center=>point:minus()
        )
      )
    )
  )
  case "path" return (
    if ($op=("union","smooth-union","round-union","chamfer-union","stepped-union"))
    then this:toSDF(path:edges($region), $canvas, $options)
    else (
      this:toSDF(path:edges($region), $canvas,
        util:merge-into($options, map {"op":"union"})
      )
    )
  )
  case "polygon" return (
    let $interpolate := $options("interpolate") return
    if ($interpolate=0) then (
      sd2:sdPolygon(
        for $pt in path:vertices($region)
        return this:normalize-to-canvas($pt, $canvas)
      )
    ) else if ($interpolate > 0) then (
      sd2:sdPolygon(
        for $pt in edge:interpolate($interpolate, path:edges($region))
        return this:normalize-to-canvas($pt, $canvas)
      )
    ) else (
      if ($op=("union","smooth-union","round-union","chamfer-union","stepped-union"))
      then this:toSDF(path:edges($region), $canvas, $options)
      else (
        this:toSDF(path:edges($region), $canvas,
          util:merge-into($options, map{"op":"union"})
        )
      )
    )
  )
  case "complex-polygon" return (
    let $inner-options := 
      if ($op=("union","smooth-union","round-union","chamfer-union","stepped-union"))
      then $options
      else util:merge-into($options, map{"op":"union"})
    let $outer := this:toSDF(cpoly:outer($region), $canvas, $options)
    let $inners := this:toSDF(cpoly:inners($region), $canvas, $inner-options)
    let $k := $options("smoothing")
    let $n := $options("steps")
    return (
      switch ($op)
      case "subtract" return this:opSubtract($outer, $inners)
      case "smooth-subtract" return this:opSmoothSubtract($outer, $inners, $k)
      case "round-subtract" return this:opRoundSubtract($outer, $inners, $k)
      case "chamfer-subtract" return this:opChamferSubtract($outer, $inners, $k)
      case "stepped-subtract" return this:opSteppedSubtract($outer, $inners, $k, $n)
      default return this:opSubtract($outer, $inners)
    )
  )
  case "sphere" return (
    let $center := this:normalize-to-canvas(solid:center($region), $canvas)
    let $r := this:scale-from-canvas(solid:radius($region), $canvas)
    let $origin := this:normalize-to-canvas($point:ORIGIN, $canvas)
    let $offset := $center=>point:minus()=>point:add($origin)
    return
      this:opTranslate(
        sd3:sdSphere($r),
        $offset
      )
  )
  case "ellipsoid" return (
    let $center := this:normalize-to-canvas(solid:center($region), $canvas)
    let $ra := this:scale-from-canvas(solid:rx($region), $canvas)
    let $rb := this:scale-from-canvas(solid:ry($region), $canvas)
    let $rc := this:scale-from-canvas(solid:rz($region), $canvas)
    let $origin := this:normalize-to-canvas($point:ORIGIN, $canvas)
    let $offset := $center=>point:minus()=>point:add($origin)
    return
      this:opTranslate(
        sd3:sdEllipsoid($ra, $rb, $rc),
        $offset
      )
  )
  case "tetrahedron" return (
    let $center := this:normalize-to-canvas(solid:center($region), $canvas)
    let $size := this:scale-from-canvas(solid:scale($region), $canvas)
    let $origin := this:normalize-to-canvas($point:ORIGIN, $canvas)
    let $offset := $center=>point:minus()=>point:add($origin)
    return             
      this:opTranslate(
        sd3:sdPyramid($size * math:sqrt(3) div 2, $size),
        $offset
      )
  )
  case "cube" return (
    let $center := this:normalize-to-canvas(solid:center($region), $canvas)
    let $size := this:scale-from-canvas(solid:scale($region), $canvas)
    let $origin := this:normalize-to-canvas($point:ORIGIN, $canvas)
    let $offset := $center=>point:minus()=>point:add($origin)
    return
      this:opTranslate(
        sd3:sdBox($size, $size, $size),
        $offset
      )
  )
  case "octahedron" return (
    let $center := this:normalize-to-canvas(solid:center($region), $canvas)
    let $size := this:scale-from-canvas(solid:scale($region), $canvas)
    let $origin := this:normalize-to-canvas($point:ORIGIN, $canvas)
    let $offset := $center=>point:minus()=>point:add($origin)
    return
      this:opTranslate(
        sd3:sdOctahedron($size),
        $offset
      )
  )
  case "face" return (
    let $vertices :=
      for $v in solid:vertices($region)[position() < last()]
      return this:normalize-to-canvas($v, $canvas)
    return (
      if (count($vertices) = 3) then (
        sd3:udTriangle($vertices[1], $vertices[2], $vertices[3])
      ) else (
        let $triangles :=
          for $i in 2 to count($vertices) - 1
          return sd3:udTriangle($vertices[1], $vertices[$i], $vertices[$i + 1])
        let $k := $options("smoothing")
        let $n := $options("steps")
        return (
          if ($op="smooth-union") then (
            fold-left(tail($triangles), head($triangles),
              function ($sdf as map(*), $triangle as map(*)) as map(*) {
                this:opSmoothUnion($sdf, $triangle, $k)
              }
            )
          ) else if ($op="round-union") then (
            fold-left(tail($triangles), head($triangles),
              function ($sdf as map(*), $triangle as map(*)) as map(*) {
                this:opRoundUnion($sdf, $triangle, $k)
              }
            )
          ) else if ($op="chamfer-union") then (
            fold-left(tail($triangles), head($triangles),
              function ($sdf as map(*), $triangle as map(*)) as map(*) {
                this:opChamferUnion($sdf, $triangle, $k)
              }
            )
          ) else if ($op="stepped-union") then (
            fold-left(tail($triangles), head($triangles),
              function ($sdf as map(*), $triangle as map(*)) as map(*) {
                this:opSteppedUnion($sdf, $triangle, $k, $n)
              }
            )
          ) else (
            this:opUnion($triangles)
          )
        )
      )
    )
  )
  case "polyhedron" return (
    if ($op=("union","smooth-union","round-union","chamfer-union","stepped-union"))
    then this:toSDF(solid:faces($region), $canvas, $options)
    else (
      this:toSDF(solid:faces($region), $canvas,
        util:merge-into($options, map{"op":"union"})
      )
    )
  )
  case "wrapper" return (
    this:toSDF(wrapper:body($region), $canvas, $options)
  )
  case "slot" return (
    (: XYZZY: NOT accounting for transforms :)
    this:toSDF(slot:body($region), $canvas, $options)
  )
  default return errors:error("GEOM-BADREGION", ($region, "regionToSDF"))
};

declare %private function this:regionToSDF(
  $region as map(xs:string,item()*),
  $op as xs:string,
  $options as map(xs:string,item()*)
) as map(*)
{
  switch (util:kind($region))
  case "point" return (
    let $size := $options("radius")
    let $pt := $region=>point:minus()
    return this:opTranslate(sd2:sdCircle($size), $pt)
  )
  case "box" return (
    if (point:same(box:min-point($region), box:max-point($region)))
    then this:regionToSDF(box:min-point($region), $op, $options)
    else (
      let $min := box:min-point($region)
      let $max := box:max-point($region)
      let $sides := point:sub($max, $min)
      let $offset := $sides=>point:times(0.5)=>point:add($min)=>point:minus()
      return
        this:opTranslate(
          sd2:sdBox(point:px($sides), point:py($sides)),
          $offset
        )
    )
  )
  case "block" return (
    let $center := box:center($region)
    let $offset := $center=>point:minus()
    return (
      this:opTranslate(
        sd3:sdBox(box:width($region), box:height($region), box:depth($region)),
        $offset
      )
    )
  )
  case "space" return (
    if (point:same(box:min-point($region), box:max-point($region)))
    then this:regionToSDF(box:min-point($region), $op, $options)
    else (
      let $min := box:min-point($region)
      let $max := box:max-point($region)
      let $sides := point:sub($max, $min)
      let $offset := $sides=>point:times(0.5)=>point:add($min)=>point:minus()
      return
        this:opTranslate(
          sd2:sdBox(point:px($sides), point:py($sides)),
          $offset
        )
    )
  )
  case "ellipse" return (
    if (ellipse:is-circle($region)) then (
      let $center := ellipse:center($region)
      let $r := ellipse:radius($region)
      return
        this:opTranslate(
          sd2:sdCircle($r),
          $center=>point:minus()
        )
    ) else (
      let $center := ellipse:center($region)
      let $a := ellipse:rx($region)
      let $b := ellipse:ry($region)
      let $rotate := ellipse:rotation($region)
      return
        this:opTranslate(
          if ($rotate = 0) 
          then sd2:sdEllipse($a, $b)
          else this:opRotate(sd2:sdEllipse($a, $b), $rotate),
          $center=>point:minus()
        ) 
    )
  )
  case "edge" return (
    if (point:same(edge:start($region), edge:end($region)))
    then this:regionToSDF(edge:start($region), $op, $options)
    else (
      let $start := edge:start($region)
      let $end := edge:end($region)
      return sd2:sdSegment($start, $end)
    )
  )
  case "quad" return (
    if (point:same(edge:start($region), edge:end($region)))
    then this:regionToSDF(edge:start($region), $op, $options)
    else (
      let $start := edge:start($region)
      let $end := edge:end($region)
      let $control := edge:controls($region)[1]
      return
        if ($options("ud"))
        then sd2:udQuad(edge:quad($start, $end, $control))
        else sd2:sdQuad(edge:quad($start, $end, $control))
    )
  )
  case "cubic" return ( (: Approximate with interpolation :)
    if (point:same(edge:start($region), edge:end($region)))
    then this:regionToSDF(edge:start($region), $op, $options)
    else (
      let $interpolate := $options("interpolate")
      let $n := if ($interpolate > 0) then $interpolate else 4
      let $edges := edge:to-edges(edge:interpolate($n, $region))
      return
        if ($op=("union","smooth-union","round-union","stepped-union","chamfer-union"))
        then this:toSDF($edges, $options)
        else (
          this:toSDF($edges, util:merge-into($options, map {"op":"union"})
        )
      )
    )
  )
  case "arc" return (
    if (point:same(edge:start($region), edge:end($region)))
    then this:regionToSDF(edge:start($region), $op, $options)
    else (
      let $arc-width := $options("width")
      let $arc := edge:as-angle-arc($region)
      let $angles := ($arc=>edge:arc-angles())!util:remap-degrees(.)
      let $large := edge:arc-large($arc)
      let $flipped := edge:arc-flipped($arc)
      let $delta := $angles[2] - $angles[1]
      let $swap :=
        ($large and
          ($delta < -180 or ($delta > 0 and $delta < 180))
        ) or
        (not($large) and
          (($delta < 0 and abs($delta) < 180) or ($delta > 180))
        )
      let $adjust-from :=
        ($large and $delta > 0 and $delta < 180) or
        (not($large) and $delta > 180)
      let $adjust-to :=
        ($large and $delta < 0 and abs($delta) <= 180) or
        (not($large) and $delta <= -180)
      let $adjusted-from :=
        if ($adjust-from) then $angles[1] + 360 else $angles[1]
      let $adjusted-to :=
        if ($adjust-to) then $angles[2] + 360 else $angles[2]
      let $from-angle :=
        if ($swap) then $adjusted-to else $adjusted-from
      let $to-angle := 
        if ($swap) then $adjusted-from else $adjusted-to
      let $center := $arc=>edge:arc-center()
      let $radius := $arc=>edge:arc-radius()
      (: 
       : If we run from A to B, orientation is 180 from angle in 
       : the middle of the range
       : orientation = 180 + (A + B)/2 
       : aperture = (B - A)
       : However, orientation of degrees is flipped, so we need to negate orientation
       :)
      let $orientation := util:remap-degrees(-(180 + ($to-angle + $from-angle) div 2.0))
      let $aperture := util:remap-degrees(($to-angle - $from-angle))
      (: center and radius already normalized :)
      let $from := point:destination($center, $from-angle, $radius)
      let $to := point:destination($center, $to-angle, $radius)
      return (
        this:opTranslate(
          if ( ($swap and not($flipped)) or (not($swap) and $flipped) ) then (
            this:opReflectXY(
              sd2:sdArc($orientation, $aperture, $radius, $arc-width),
              $from,
              $to
            )
          ) else (
            sd2:sdArc($orientation, $aperture, $radius, $arc-width)
          ),
          $center=>point:minus()
        )
      )
    )
  )
  case "path" return (
    if ($op=("union","smooth-union","round-union","chamfer-union","stepped-union"))
    then this:toSDF(path:edges($region), $options)
    else (
      this:toSDF(path:edges($region),
        util:merge-into($options, map {"op":"union"})
      )
    )
  )
  case "polygon" return (
    let $interpolate := $options("interpolate") return
    if ($interpolate=0) then (
      sd2:sdPolygon(path:vertices($region))
    ) else if ($interpolate > 0) then (
      sd2:sdPolygon(
        edge:interpolate($interpolate, path:edges($region))
      )
    ) else (
      if ($op=("union","smooth-union","round-union","chamfer-union","stepped-union"))
      then this:toSDF(path:edges($region), $options)
      else (
        this:toSDF(path:edges($region), 
          util:merge-into($options, map{"op":"union"})
        )
      )
    )
  )
  case "complex-polygon" return (
    let $inner-options := 
      if ($op=("union","smooth-union","round-union","chamfer-union","stepped-union"))
      then $options
      else util:merge-into($options, map{"op":"union"})
    let $outer := this:toSDF(cpoly:outer($region), $options)
    let $inners := this:toSDF(cpoly:inners($region), $inner-options)
    let $k := $options("smoothing")
    let $n := $options("steps")
    return (
      switch ($op)
      case "subtract" return this:opSubtract($outer, $inners)
      case "smooth-subtract" return this:opSmoothSubtract($outer, $inners, $k)
      case "round-subtract" return this:opRoundSubtract($outer, $inners, $k)
      case "chamfer-subtract" return this:opChamferSubtract($outer, $inners, $k)
      case "stepped-subtract" return this:opSteppedSubtract($outer, $inners, $k, $n)
      default return this:opSubtract($outer, $inners)
    )
  )
  case "sphere" return (
    let $center := solid:center($region)
    let $r := solid:radius($region)
    let $offset := $center=>point:minus()
    return
      this:opTranslate(
        sd3:sdSphere($r),
        $offset
      )
  )
  case "ellipsoid" return (
    let $center := solid:center($region)
    let $ra := solid:rx($region)
    let $rb := solid:ry($region)
    let $rc := solid:rz($region)
    let $offset := $center=>point:minus()
    return
      this:opTranslate(
        sd3:sdEllipsoid($ra, $rb, $rc),
        $offset
      )
  )
  case "tetrahedron" return (
    let $center := solid:center($region)
    let $size := solid:scale($region)
    let $offset := $center=>point:minus()
    return             
      this:opTranslate(
        sd3:sdPyramid($size * math:sqrt(3) div 2, $size),
        $offset
      )
  )
  case "cube" return (
    let $center := solid:center($region)
    let $size := solid:scale($region)
    let $offset := $center=>point:minus()
    return
      this:opTranslate(
        sd3:sdBox($size, $size, $size),
        $offset
      )
  )
  case "octahedron" return (
    let $center := solid:center($region)
    let $size := solid:scale($region)
    let $offset := $center=>point:minus()
    return
      this:opTranslate(
        sd3:sdOctahedron($size),
        $offset
      )
  )
  case "face" return (
    let $vertices := solid:vertices($region)[position() < last()]
    return (
      if (count($vertices) = 3) then (
        sd3:udTriangle($vertices[1], $vertices[2], $vertices[3])
      ) else (
        let $triangles :=
          for $i in 2 to count($vertices) - 1
          return sd3:udTriangle($vertices[1], $vertices[$i], $vertices[$i + 1])
        let $k := $options("smoothing")
        let $n := $options("steps")
        return (
          if ($op="smooth-union") then (
            fold-left(tail($triangles), head($triangles),
              function ($sdf as map(*), $triangle as map(*)) as map(*) {
                this:opSmoothUnion($sdf, $triangle, $k)
              }
            )
          ) else if ($op="round-union") then (
            fold-left(tail($triangles), head($triangles),
              function ($sdf as map(*), $triangle as map(*)) as map(*) {
                this:opRoundUnion($sdf, $triangle, $k)
              }
            )
          ) else if ($op="chamfer-union") then (
            fold-left(tail($triangles), head($triangles),
              function ($sdf as map(*), $triangle as map(*)) as map(*) {
                this:opChamferUnion($sdf, $triangle, $k)
              }
            )
          ) else if ($op="stepped-union") then (
            fold-left(tail($triangles), head($triangles),
              function ($sdf as map(*), $triangle as map(*)) as map(*) {
                this:opSteppedUnion($sdf, $triangle, $k, $n)
              }
            )
          ) else (
            this:opUnion($triangles)
          )
        )
      )
    )
  )
  case "polyhedron" return (
    if ($op=("union","smooth-union","round-union","chamfer-union","stepped-union"))
    then this:toSDF(solid:faces($region), $options)
    else (
      this:toSDF(solid:faces($region), 
        util:merge-into($options, map{"op":"union"})
      )
    )
  )
  case "wrapper" return (
    this:toSDF(wrapper:body($region), $options)
  )
  case "slot" return (
    (: XYZZY: NOT accounting for transforms :)
    this:toSDF(slot:body($region), $options)
  )
  default return errors:error("GEOM-BADREGION", ($region, "regionToSDF"))
};

declare function this:toSDF(
  $regions as map(xs:string,item()*)*,
  $options as map(xs:string,item()*)
) as map(*)
{
  this:toSDF($regions, (), $options)
};

declare function this:toSDF(
  $regions as map(xs:string,item()*)*,
  $canvas as map(xs:string,item()*)?,
  $options as map(xs:string,item()*)
) as map(*)
{
  let $options :=
    util:merge-into((
      map {
        "op": "union",
        "smoothing": 0.05,
        "interpolate": 0,
        "steps": 3,
        "radius": 1,
        "width": 1,
        "ud": false()
      },
      $options
    ))
  let $op := $options("op")
  let $k := $options("smoothing")
  let $n := $options("steps")
  let $translated-regions :=
    if (exists($canvas)) then (
      for $region in $regions return (
        this:regionToSDF($region, $op, $canvas, $options)
      )
    ) else (
      for $region in $regions return (
        this:regionToSDF($region, $op, $options)
      )
    )
  return (
    if ($op="union") then (
      this:opUnion($translated-regions)
    ) else if ($op="intersect") then (
      this:opIntersect($translated-regions)
    ) else (
      fold-left(
        tail($translated-regions), head($translated-regions),
        function ($sdf as map(*), $current as map(*)) as map(*) {
          if (empty($sdf)) then $current
          else switch($op)
          case "intersect" return this:opIntersect($sdf, $current)
          case "subtract" return this:opSubtract($sdf, $current)
          case "round-union" return this:opRoundUnion($sdf, $current, $k)
          case "round-intersect" return this:opRoundIntersect($sdf, $current, $k)
          case "round-subtract" return this:opRoundSubtract($sdf, $current, $k)
          case "smooth-union" return this:opSmoothUnion($sdf, $current, $k)
          case "smooth-intersect" return this:opSmoothIntersect($sdf, $current, $k)
          case "smooth-subtract" return this:opSmoothSubtract($sdf, $current, $k)
          case "chamfer-union" return this:opChamferUnion($sdf, $current, $k)
          case "chamfer-intersect" return this:opChamferIntersect($sdf, $current, $k)
          case "chamfer-subtract" return this:opChamferSubtract($sdf, $current, $k)
          case "stepped-union" return this:opSteppedUnion($sdf, $current, $k, $n)
          case "stepped-intersect" return this:opSteppedIntersect($sdf, $current, $k, $n)
          case "stepped-subtract" return this:opSteppedSubtract($sdf, $current, $k, $n)
          default (: union :) return this:opUnion($sdf, $current)
        }
      )
    )
  )
};