@@ 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>
@@ 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"
+
+;