96f4f461bf87 — Leonard Ritter a month ago
* voxel BVH
2 files changed, 480 insertions(+), 42 deletions(-)

M testing/voxel_bvh.sc
A => testing/voxel_bvh_bf.sc
M testing/voxel_bvh.sc +102 -42
@@ 13,7 13,10 @@ fn df (p)
     op := p
     p :=
         versor-rotate
-            versor (vec3 0 0 1) (radians 45.0)
+            versor (vec3 0 0 1)
+                radians 45.0
+                #radians 11.0
+                #radians 0.0
             p
     #sdOr
         sdSub

          
@@ 39,7 42,7 @@ fn df (p)
 
 let MAPDIM = 64
 
-local distmap : (Array f32)
+global distmap : (Array f32)
 'resize distmap (MAPDIM * MAPDIM) 0.0
 
 fn index (x y)

          
@@ 64,6 67,10 @@ case (bounds)
     let w h = (x2 - x1) (y2 - y1)
     w * h
 
+fn rectsize (bounds)
+    let x1 y1 x2 y2 = (unpack bounds)
+    vec2 (x2 - x1) (y2 - y1)
+
 fn diagvolume (bounds)
     let x1 y1 x2 y2 = (unpack bounds)
     let w h = (x2 - x1) (y2 - y1)

          
@@ 90,6 97,16 @@ for x y in (dim MAPDIM MAPDIM)
     idx := (index x y)
     distmap @ idx = d
 
+fn printrect (bounds)
+    let x1 y1 x2 y2 = (unpack bounds)
+    for y in (range y1 y2)
+        for x in (range x1 x2)
+            let d = (distmap @ (index x y))
+            io-write!
+                if (d <= 0) "Oo"
+                else "._"
+        io-write! "\n"
+
 for y in (range MAPDIM)
     for x in (range MAPDIM)
         let d = (distmap @ (index x y))

          
@@ 108,15 125,49 @@ for x y in (dim MAPDIM MAPDIM)
     s11 := (g x y)
     idx := (index x y)
     deltamap @ idx =
-        do
-            if (s11 == 1)
+        do  # count concave corners
+            if (s11 == 0)
+                s0 := (g (x + 1) y)
+                s1 := (g x (y + 1))
+                s2 := (g (x - 1) y)
+                s3 := (g x (y - 1))
+                +
+                    s0 * s1
+                    s1 * s2
+                    s2 * s3
+                    s3 * s0
+            # count convex corners
+
+            #if (s11 == 0)
+                s0 := (g (x + 1) y) ^ 1
+                s1 := (g x (y + 1)) ^ 1
+                s2 := (g (x - 1) y) ^ 1
+                s3 := (g x (y - 1)) ^ 1
+                +
+                    s0 * s1
+                    s1 * s2
+                    s2 * s3
+                    s3 * s0
+            else 0
+
+#for x y in (dim MAPDIM MAPDIM)
+    inline g (x y)
+        ? ((distmap @ (index x y)) <= 0) 1 0
+    #s00 := (g (x - 1) (y - 1))
+    s11 := (g x y)
+    idx := (index x y)
+    deltamap @ idx +=
+        do  # count lines
+            if (s11 == 0)
+                s0 := (g (x + 1) y)
+                s1 := (g x (y + 1))
+                s2 := (g (x - 1) y)
+                s3 := (g x (y - 1))
                 +
                     va-map
-                        inline (x) (abs (s11 - x))
-                        g (x - 1) y
-                        g (x + 1) y
-                        g x (y - 1)
-                        g x (y + 1)
+                        inline (x)
+                            abs (s11 - x)
+                        \ s0 s1 s2 s3
             else 0
 
 for y in (range MAPDIM)

          
@@ 228,15 279,15 @@ write-string outp
         <rect x="0" y="0" width="800" height="800" fill="white"/>
         <g transform="scale(12.5,12.5) translate(0,0)">
 
-fn descend (level sat sat2 outp bounds numblocks maxlevel)
+fn descend (id level sat sat2 outp bounds numblocks maxlevel)
     let dv = (volume sat bounds)
     if (dv == 0) # empty
         return;
-    #if (level > 3)
+    #if (level > 4)
         return;
     let v = (volume bounds)
     let trim = (v - dv)
-    #print "level" level "bounds" bounds "trim" trim "volume" dv
+    print id "level" level "bounds" bounds "trim" trim "volume" dv
     let x1 y1 x2 y2 = (unpack bounds)
     let w h = (x2 - x1) (y2 - y1)
     let m = 0.0

          
@@ 285,6 336,8 @@ fn descend (level sat sat2 outp bounds n
 
             let lsq = (sqvolume lbounds)
             let rsq = (sqvolume rbounds)
+            let lsa = (sqaspect lbounds)
+            let rsa = (sqaspect rbounds)
 
             let lp = (perimeter lbounds)
             let rp = (perimeter rbounds)

          
@@ 294,15 347,24 @@ fn descend (level sat sat2 outp bounds n
             #score := lsq + rsq + ((lb - lv) + (rb - rv)) ** 2
             #score := (sqrt ((lsq + rsq) as f32)) + ((lb - lv) ** 2 + (rb - rv) ** 2) as f32
             let lm rm = (lb - lv) (rb - rv)
-            #score := lm ** 2 + rm ** 2
+            #err_trim := lm ** 2 + rm ** 2
+            err_trim := ((lm / b) ** 2) + ((rm / b) ** 2)
+
+            print id err_trim
 
             let lbb rbb = (lp - lo) (rp - ro)
             score :=
-                min
-                    lm ** 2 + rm ** 2
-                    +
-                        (lo - ro) ** 2
-                        lbb ** 2 + rbb ** 2
+                do
+                    # prefer perfect solutions
+                    if ((lm == rm) & (lm == 0)) 0.0
+                    #elseif (((min lo ro) == 0) & ((max lo ro) == 4)) 1
+                    else
+                        #2 + (max (lo - 1) (ro - 1))
+                        ((2 + (max (lo // 2) (ro // 2))) as f32) #+ err_trim
+                        #2 + (max lo ro)
+                        #(lo - ro) ** 2
+                        #err_trim
+
             #score := lm ** 2 + rm ** 2
             #score := (lo - ro) ** 2 + lbb ** 2 + rbb ** 2
             #score := (lo - ro) ** 2

          
@@ 334,17 396,19 @@ fn descend (level sat sat2 outp bounds n
             testsplit x xl
                 testsplit x xr bestscore b1 b2
 
+    let w h = (x2 - x1) (y2 - y1)
     #let score b1 b2 =
-        if ((level & 1) == 0)
-            # horizontal split
-            split inf (ivec4 0) (ivec4 0) (y1 + 1) y2
-                inline (y)
-                    _ (ivec4 x1 y1 x2 y) (ivec4 x1 y x2 y2)
-        else
-            # vertical split
-            split inf (ivec4 0) (ivec4 0) (x1 + 1) x2
-                inline (x)
-                    _ (ivec4 x1 y1 x y2) (ivec4 x y1 x2 y2)
+        do
+            if (w < h) #((level & 1) == 1)
+                # horizontal split
+                split inf (ivec4 0) (ivec4 0) (y1 + 1) y2
+                    inline (y)
+                        _ (ivec4 x1 y1 x2 y) (ivec4 x1 y x2 y2)
+            else
+                # vertical split
+                split inf (ivec4 0) (ivec4 0) (x1 + 1) x2
+                    inline (x)
+                        _ (ivec4 x1 y1 x y2) (ivec4 x y1 x2 y2)
 
 
     let b1 b2 =

          
@@ 359,34 423,30 @@ fn descend (level sat sat2 outp bounds n
                 split inf (ivec4 0) (ivec4 0) (x1 + 1) x2
                     inline (x)
                         _ (ivec4 x1 y1 x y2) (ivec4 x y1 x2 y2)
+            print id hscore vscore
             if (hscore < vscore)
                 _ hb1 hb2
             elseif (hscore > vscore)
                 _ vb1 vb2
-            #else
-                if ((level & 1) == 0)
+            else
+                shb := (rectsize hb1) + (rectsize hb2)
+                svb := (rectsize vb1) + (rectsize vb2)
+                if (w < h)
                     _ hb1 hb2
-                else
+                elseif (w > h)
                     _ vb1 vb2
-            else
-                let hv1 = (sqaspect hb1)
-                let hv2 = (sqaspect hb2)
-                let vv1 = (sqaspect vb1)
-                let vv2 = (sqaspect vb2)
-                e0 := hv1 + hv2
-                e1 := vv1 + vv2
-                if (e0 < e1)
+                elseif ((level & 1) == 0)
                     _ hb1 hb2
                 else
                     _ vb1 vb2
 
     level := level + 1
-    this-function level sat sat2 outp b1 numblocks maxlevel
-    this-function level sat sat2 outp b2 numblocks maxlevel
+    this-function (id << 1) level sat sat2 outp b1 numblocks maxlevel
+    this-function ((id << 1) + 1) level sat sat2 outp b2 numblocks maxlevel
 
 local numblocks = 0
 local maxlevel = 0
-descend 0 sat sat2 (view outp) bounds numblocks maxlevel
+descend 0 0 sat sat2 (view outp) bounds numblocks maxlevel
 
 write-string outp
     """"</g>

          
A => testing/voxel_bvh_bf.sc +378 -0
@@ 0,0 1,378 @@ 
+
+using import glm
+using import itertools
+using import Array
+
+import ..lib.tukan.use
+using import tukan.sdf
+using import tukan.File
+using import tukan.color
+using import tukan.rotation
+
+fn df (p)
+    op := p
+    p :=
+        versor-rotate
+            versor (vec3 0 0 1) (radians 0.0)
+            p
+    sdOr
+        sdSub
+            sdSphere p 0.8
+            sdSphere (p - (vec3 0.3 -0.3 0)) 0.8
+        sdTorus (vec3 (p.x - 0.3) 0 (p.y + 0.3)) (vec2 0.4 0.1)
+    #sdOr
+        sdBox op (vec3 0.3)
+        sdSub
+            sdBox p (vec3 0.6)
+            sdSphere p 0.5
+    #sdSub
+        sdBox p (vec3 0.6)
+        sdBox op (vec3 0.4)
+    #fold (d = inf) for i x y in (enumerate (dim 4 4))
+        q := (((vec2 x y) / 3.0) * 2.0 - 1.0) * 0.8
+        sdOr d
+            sdSphere (p - (vec3 q 0)) (0.0625 + (i as f32 / 16.0) * 0.0625)
+    #sdSmoothOr
+        sdTorus ((p - (vec3 -0.5 0 0)) . xzy) (vec2 0.35 0.05)
+        sdTorus ((p - (vec3 0.5 0 0)) . xzy) (vec2 0.35 0.05)
+        0.6
+
+let MAPDIM = 64
+
+local distmap : (Array f32)
+'resize distmap (MAPDIM * MAPDIM) 0.0
+
+fn index (x y)
+    (mod y MAPDIM) * MAPDIM + (mod x MAPDIM)
+
+fn sat@ (sat x y)
+    if ((x < 0) | (y < 0)) 0
+    else (deref (sat @ (index x y)))
+
+fn... volume (sat bounds)
+    let x1 y1 x2 y2 =
+        va-map
+            inline (x) (x - 1)
+            unpack bounds
+    p00 := (sat@ sat x1 y1)
+    p10 := (sat@ sat x2 y1)
+    p01 := (sat@ sat x1 y2)
+    p11 := (sat@ sat x2 y2)
+    p00 + p11 - p10 - p01
+case (bounds)
+    let x1 y1 x2 y2 = (unpack bounds)
+    let w h = (x2 - x1) (y2 - y1)
+    w * h
+
+fn diagvolume (bounds)
+    let x1 y1 x2 y2 = (unpack bounds)
+    let w h = (x2 - x1) (y2 - y1)
+    abs (w - h)
+
+fn sqvolume (bounds)
+    let x1 y1 x2 y2 = (unpack bounds)
+    let w h = (x2 - x1) (y2 - y1)
+    w * w + h * h
+
+fn perimeter (bounds)
+    let x1 y1 x2 y2 = (unpack bounds)
+    let w h = (x2 - x1) (y2 - y1)
+    2 * (w + h)
+
+fn sqaspect (bounds)
+    let x1 y1 x2 y2 = (unpack bounds)
+    let w h = (x2 - x1) (y2 - y1)
+    (w - h) ** 2
+
+for x y in (dim MAPDIM MAPDIM)
+    p := ((vec2 x y) / (MAPDIM - 1)) * 2.0 - 1.0
+    d := (df (vec3 p 0.0))
+    idx := (index x y)
+    distmap @ idx = d
+
+for y in (range MAPDIM)
+    for x in (range MAPDIM)
+        let d = (distmap @ (index x y))
+        io-write!
+            if (d <= 0) "Oo"
+            else "._"
+    io-write! "\n"
+
+local deltamap : (Array i32)
+'resize deltamap (MAPDIM * MAPDIM) 0
+
+for x y in (dim MAPDIM MAPDIM)
+    inline g (x y)
+        ? ((distmap @ (index x y)) <= 0) 1 0
+    #s00 := (g (x - 1) (y - 1))
+    s11 := (g x y)
+    idx := (index x y)
+    deltamap @ idx =
+        do
+            if (s11 == 1)
+                +
+                    va-map
+                        inline (x) (abs (s11 - x))
+                        g (x - 1) y
+                        g (x + 1) y
+                        g x (y - 1)
+                        g x (y + 1)
+            else 0
+
+for y in (range MAPDIM)
+    for x in (range MAPDIM)
+        let s = (deltamap @ (index x y))
+        if (s < 16)
+            io-write! "0"
+        io-write! (hex s)
+    io-write! "\n"
+
+inline makesat (dst src qf)
+    for x y in (dim MAPDIM MAPDIM)
+        idx := (index x y)
+        d := (src @ idx)
+        let sum =
+            -
+                +
+                    qf d
+                    sat@ dst (x - 1) y
+                    sat@ dst x (y - 1)
+                sat@ dst (x - 1) (y - 1)
+        dst @ idx = sum
+
+# summed-area table - volume
+local sat : (Array i32)
+'resize sat (MAPDIM * MAPDIM) 0
+
+makesat sat distmap
+    inline (d)
+        ? (d <= 0) 1 0
+
+# summed-area table - perimeter
+local sat2 : (Array i32)
+'resize sat2 (MAPDIM * MAPDIM) 0
+
+makesat sat2 deltamap
+    inline (d) d
+
+for y in (range MAPDIM)
+    for x in (range MAPDIM)
+        let s = (sat @ (index x y))
+        if (s < 16)
+            io-write! "0"
+        io-write! (hex s)
+    io-write! "\n"
+
+################################################################################
+
+fn tighten (sat bounds)
+    let v = (volume sat bounds)
+    let x1 y1 x2 y2 = (unpack bounds)
+    if (v == 0)
+        return (ivec4 x1 y1 x1 y1)
+    #assert (y2 > y1)
+    #assert (x2 > x1)
+    inline optimize (x1 x2 rnd cmpf vecf)
+        loop (l r = x1 x2)
+            x := (l + r + rnd) // 2
+            if ((r - l) <= 1)
+                break x
+            nv := (volume sat (vecf x))
+            if (cmpf nv v)
+                repeat l x
+            else
+                repeat x r
+    let x1 =
+        optimize x1 x2 0 (_ <) (inline (x) (ivec4 x y1 x2 y2))
+    let x2 =
+        optimize x1 x2 1 (_ >=) (inline (x) (ivec4 x1 y1 x y2))
+    #assert (x1 < x2) (.. (repr x1) " < " (repr x2))
+    let y1 =
+        optimize y1 y2 0 (_ <) (inline (y) (ivec4 x1 y x2 y2))
+    let y2 =
+        optimize y1 y2 1 (_ >=) (inline (y) (ivec4 x1 y1 x2 y))
+    #assert (y1 < y2) (.. (repr y1) " < " (repr y2))
+    ivec4 x1 y1 x2 y2
+
+let rootbounds = (ivec4 0 0 MAPDIM MAPDIM)
+let rootvolume = (volume sat rootbounds)
+let bounds = (tighten sat rootbounds)
+assert ((volume sat bounds) == rootvolume)
+assert (bounds == (tighten sat bounds))
+
+print "volume:" rootvolume
+print "tightened:" bounds "trim:" ((volume bounds) - rootvolume)
+
+let outp =
+    try (File.open (.. module-dir "/voxel_bvh.svg") "w")
+    else
+        error "failed to open file"
+
+fn write-string (outp str)
+    'write outp (str as rawstring) (countof str)
+
+write-string outp
+    """"<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+        <!-- Created with Inkscape (http://www.inkscape.org/) -->
+        <svg
+        xmlns:dc="http://purl.org/dc/elements/1.1/"
+        xmlns:cc="http://creativecommons.org/ns#"
+        xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
+        xmlns:svg="http://www.w3.org/2000/svg"
+        xmlns="http://www.w3.org/2000/svg"
+        xmlns:xlink="http://www.w3.org/1999/xlink"
+        width="800"
+        height="800"
+        id="svg">
+
+        <rect x="0" y="0" width="800" height="800" fill="white"/>
+        <g transform="scale(12.5,12.5) translate(0,0)">
+
+using import Map
+local cache : (Map (tuple i32 i32 i32 i32 i32) (tuple i32 i32))
+
+fn descend (level cache sat sat2 outp bounds test?)
+    # max level, num blocks
+    returning i32 i32
+    key := (tupleof level (unpack bounds))
+    if test?
+        try
+            return (unpack (copy ('get cache key)))
+        else;
+    let _descend = this-function
+    let dv = (volume sat bounds)
+    if ((dv == 0) | (level > 12)) # empty
+        'set cache key (tupleof 1000 1000)
+        return 1000 1000 # bad
+    #if (level > 3)
+        return;
+    let v = (volume bounds)
+    let trim = (v - dv)
+    let x1 y1 x2 y2 = (unpack bounds)
+    let w h = (x2 - x1) (y2 - y1)
+    let m = 0.0
+    if (not test?)
+        print "level" level "bounds" bounds "trim" trim "volume" dv
+        #if (((level % 2) == 0) | (trim == 0))
+        #let m = (level as f32 * 0.01)
+        write-string outp
+            .. "<rect x=\""
+                tostring (x1 as f32 + m)
+                "\" y=\""
+                tostring (y1 as f32 + m)
+                "\" width=\""
+                tostring (w as f32 - (2 * m))
+                "\" height=\""
+                tostring (h as f32 - (2 * m))
+                "\" fill=\""
+                do
+                    color := (viridis (level as f32 / 9.0))
+                    let color =
+                        if (trim == 0) color
+                        else (mix color (vec3 1.0) 0.5)
+                    let r g b = (unpack (color * 255.0))
+                    .. "rgb(" (tostring r) "," (tostring g) "," (tostring b) ")"
+                "\" stroke=\"#000000\" stroke-width=\"0.05px\"/>\n"
+    if (trim == 0) # done
+        'set cache key (tupleof level 1)
+        return level 1
+
+    inline split (score b1 b2 r1 r2 splitf)
+        inline testsplit (x bestscore b1 b2)
+            let lbounds rbounds = (splitf x)
+            let lbounds = (tighten sat lbounds)
+            let rbounds = (tighten sat rbounds)
+
+            let sublevel = (level + 1)
+            let level0 numblocks0 = (_descend sublevel cache sat sat2 outp lbounds true)
+            let level1 numblocks1 = (_descend sublevel cache sat sat2 outp rbounds true)
+
+
+            score := (max level0 level1)
+
+            score as:= f32
+            if (score < bestscore)
+                _ score lbounds rbounds
+            else
+                _ bestscore b1 b2
+
+        c := (r1 + r2) // 2
+        w := (r2 - r1 + 1) // 2
+        #if (level < 1)
+            let x = (clamp c r1 (r2 - 1))
+            let lbounds rbounds = (splitf x)
+            let lbounds = (tighten sat lbounds)
+            let rbounds = (tighten sat rbounds)
+            return 0.0 lbounds rbounds
+
+        fold (bestscore b1 b2 = score b1 b2)
+            \ for x in (range 0 w)
+            xl := (clamp (c - x - 1) r1 (r2 - 1))
+            xr := (clamp (c + x) r1 (r2 - 1))
+            testsplit xl
+                testsplit xr bestscore b1 b2
+
+    let b1 b2 =
+        do
+            # horizontal split
+            let hscore hb1 hb2 =
+                split inf (ivec4 0) (ivec4 0) (y1 + 1) y2
+                    inline (y)
+                        _ (ivec4 x1 y1 x2 y) (ivec4 x1 y x2 y2)
+            # vertical split
+            let vscore vb1 vb2 =
+                split inf (ivec4 0) (ivec4 0) (x1 + 1) x2
+                    inline (x)
+                        _ (ivec4 x1 y1 x y2) (ivec4 x y1 x2 y2)
+            if (hscore < vscore)
+                _ hb1 hb2
+            elseif (hscore > vscore)
+                _ vb1 vb2
+            #else
+                if ((level & 1) == 0)
+                    _ hb1 hb2
+                else
+                    _ vb1 vb2
+            else
+                let hv1 = (sqaspect hb1)
+                let hv2 = (sqaspect hb2)
+                let vv1 = (sqaspect vb1)
+                let vv2 = (sqaspect vb2)
+                e0 := hv1 + hv2
+                e1 := vv1 + vv2
+                if (e0 < e1)
+                    _ hb1 hb2
+                else
+                    _ vb1 vb2
+
+    local newlevel = level
+    local numblocks = 0
+    level := level + 1
+
+    let level0 numblocks0 = (this-function level cache sat sat2 outp b1 test?)
+    let level1 numblocks1 = (this-function level cache sat sat2 outp b2 test?)
+
+    if (numblocks0 > 0)
+        newlevel = (max newlevel level0)
+        numblocks = (numblocks + numblocks0)
+    if (numblocks1 > 0)
+        newlevel = (max newlevel level1)
+        numblocks = (numblocks + numblocks1)
+    if (numblocks > 0)
+        # count this block
+        numblocks += 1
+
+    'set cache key (tupleof newlevel numblocks)
+    return newlevel numblocks
+
+let maxlevel numblocks =
+    descend 0 cache sat sat2 (view outp) bounds false
+
+write-string outp
+    """"</g>
+        </svg>
+drop outp
+
+print numblocks "blocks," (maxlevel + 1) "levels"
+
+;