c1eadf26b4b5 — Leonard Ritter tip a month ago
* inverse spherical functions
4 files changed, 191 insertions(+), 43 deletions(-)

M lib/tukan/normal.sc
A => lib/tukan/spherical.sc
M testing/test_subjective.sc
M testing/test_subjective_mt.sc
M lib/tukan/normal.sc +2 -1
@@ 6,7 6,8 @@ using import glm
 fn oct_wrap (v)
     *
         1.0 - (abs v.yx)
-        vec2
+        sign v
+        #vec2
             ? (v.x >= 0.0) 1.0 -1.0
             ? (v.y >= 0.0) 1.0 -1.0
 

          
A => lib/tukan/spherical.sc +110 -0
@@ 0,0 1,110 @@ 
+
+using import glm
+
+FIX_EQUAL_AREA := true
+
+# everything spherical coordinates
+
+inline oct_wrap (v)
+    (1.0 - (abs v.yx)) * (sign v)
+
+# vec2 -> vec3
+fn octahedral-surface (uv)
+    z := 1.0 - (abs uv.x) - (abs uv.y)
+    vec3
+        ? (z >= 0.0) uv (oct_wrap uv)
+        z
+
+# n is a point on the L1 surface
+# vec3 -> vec2
+inline octahedral-coordinates (n)
+    #n = L1normalize(n); # not necessary if point is already normalized
+    ? (n.z >= 0.0) n.xy (oct_wrap n.xy)
+
+# spherical coordinates from points on L1 surface
+# vec3 -> vec2
+fn L1-spherical (o)
+    # L1 to L2 arclength
+    r := (abs o.x) + (abs o.y)
+    h := o.x * (? (r == 0.0) 0.0 (1.0 / r)) - 1.0
+    phi := (? (o.y >= 0.0) -h h) * (pi * 0.5)
+    # L1 to L2 radius
+    let theta =
+        static-if FIX_EQUAL_AREA
+            # inverse of https://en.wikipedia.org/wiki/Collignon_projection
+            z := 1.0 - (abs o.z)
+            (pi * 0.5) - (asin (1.0 - z * z)) * (sign o.z)
+        else
+            (1.0 - o.z) * (pi * 0.5)
+    vec2 phi theta
+
+# point on L1 surface from spherical coordinates
+fn spherical-L1 (uv)
+    # L2 to L1 radius
+    let z =
+    static-if FIX_EQUAL_AREA
+        # forward formula of https://en.wikipedia.org/wiki/Collignon_projection
+        theta := (pi * 0.5) - uv.y
+        (1.0 - (sqrt (1.0 - (sin (abs theta))))) * (sign theta)
+    else
+        1.0 - uv.y / (pi * 0.5)
+    # L2 to L1 arclength
+    a := uv.x / (pi * 0.25)
+    r := 1.0 - (abs z)
+    a := (mod (a + 4.0) 8.0) - 4.0
+    cos_x := 1.0 - (abs a) / 2.0
+    xy := r * (vec2 cos_x ((sign a) * (1.0 - (abs cos_x))))
+    vec3 xy z
+
+# uv is [-pi..pi]x[0..pi]
+# vec2 -> vec3
+fn spherical-surface (uv)
+    phi_s := (sin uv.x)
+    phi_c := (cos uv.x)
+    theta_s := (sin uv.y)
+    theta_c := (cos uv.y)
+    vec3
+        theta_s * phi_c
+        theta_s * phi_s
+        theta_c
+
+# n is a normal vector
+fn spherical-coordinates (n)
+    vec2 (atan n.y n.x) (acos n.z)
+
+# [0..1]x[0..1]+offset -> [-1..1]x[-1..1]
+# https://en.wikipedia.org/wiki/Peirce_quincuncial_projection
+# vec2 -> vec2
+fn tile-quincuncial (p)
+    p := (mod p 2.0)
+    q := (ivec2 p)
+    uv := p - (vec2 q)
+    (uv * 2.0 - 1.0) * (? (((q.x ^ q.y) & 1) != 0) -1.0 1.0)
+
+# [0..2]x[0..1]+offset -> [-1..1]x[-1..1]
+# https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection
+# vec2 -> vec2
+inline tile-guyou (p)
+    tile-quincuncial
+        (vec2 (p.x + p.y) (p.y - p.x - 1.0)) / 2.0
+
+inline L1normalize (n)
+    n / ((abs n.x) + (abs n.y) + (abs n.z))
+
+# [-1..1]x[-1..1] -> [0..2]x[0..1]
+fn inverse-guyou (vec2 o)
+    o := (vec2 (o.x - o.y + 1.0 / 2.0) (o.x + o.y - 1.0 / 2.0))
+    o := o + (vec2 -0.5 0.5)
+    let o =
+        if (o.y >= 1.0)
+            o = 2.0 - o
+        elseif (o.x < -1.0)
+            o.x = o.x + 4.0
+        elseif (o.y < -1.0)
+            o = (vec2 2.0 -2.0) - o
+    (o + 1.0) / 2.0
+
+do
+    let octahedral-surface L1-spherical spherical-surface spherical-coordinates
+    let spherical-L1 tile-guyou tile-quincuncial inverse-guyou octahedral-coordinates
+    locals;
  No newline at end of file

          
M testing/test_subjective.sc +19 -6
@@ 11,6 11,7 @@ using import tukan.gl
 using import tukan.color
 using import tukan.normal
 using import tukan.isosurface
+using import tukan.spherical
 using import tukan.ResourceGroup
 using import .testfragment
 

          
@@ 74,7 75,7 @@ fn unpack_normal_snorm (encN)
             ? (z >= 0.0) encN.xy (oct_wrap encN.xy)
             z
 
-N := 32
+N := 16
 
 fn persp1 (p)
     static-if true

          
@@ 85,7 86,19 @@ fn persp1 (p)
         #r := (r / (1.0 - (abs r))) #+ 1.0
         #r := (-(log2 (1.0 - r)))
         #r := (/ r)
-        (unpack_normal_snorm p.xy) * r
+
+        #p := (unpack_normal_snorm p.xy)
+        #p :=
+            spherical-surface
+                L1-spherical
+                    octahedral-surface p.xy
+        p :=
+            spherical-surface
+                L1-spherical
+                    octahedral-surface
+                        tile-guyou ((p.xy - 1.0) / 2.0)
+
+        p * r
     else
         p := (expand p)
         p := p / rcp_c

          
@@ 156,14 169,14 @@ fn draw (size pg-test frame)
             p * 4.0
 
     #wireframe;
-    for x y z in (dim N N N)
+    for x y z in (dim (2 * N) N N)
         let bias = 0.5
         #let bias = 0.0
         c := ((vec3 x y z) + bias) / N
         pc := c * 2 - 1
         #r := c.z * 0.8
-        #r := 0.5
-        r := 1.0
+        r := 0.5
+        #r := 1.0
         p0 := pc - d * r
         p1 := pc + d * r
         color c

          
@@ 191,7 204,7 @@ fn draw (size pg-test frame)
                 tf (vec3 p1.x p0.y p1.z)
                 tf (vec3 p1.x p1.y p0.z)
                 tf (vec3 p1.x p1.y p1.z)
-        static-if false
+        static-if true
             cuboid (unpack verts)
         elseif true
             for i in (range 6)

          
M testing/test_subjective_mt.sc +60 -36
@@ 41,18 41,19 @@ using import tukan.projection
 using import tukan.derivative
 using import tukan.isosurface
 using import tukan.hash
+using import tukan.spherical
 using import tukan.ResourceGroup
 using import .testfragment
 
 let RG = (ResourceGroup "RG")
 from (import tukan.math) let expmix
 
-SAMPLE_CAMERA_OFFSET := true
-PROJECT_FINAL_VERTEX := true
+SAMPLE_CAMERA_OFFSET := false
+PROJECT_FINAL_VERTEX := false
 VISUALIZE_IDS := false
 POST_TRANSFORM := false # true is worse
-OCCLUSION_CULLING := true
-FOG := true
+OCCLUSION_CULLING := false
+FOG := false
 
 # to reach fog density D at depth Z, FOG_RATE = -log2(1 - D)/Z
 FOG_RATE := 0.01 # 50% at 100 units

          
@@ 61,6 62,7 @@ let MAX_VERTICES = (20 * (1 << 20))
 
 let WORLD_SIZE = (uvec3 256)
 let WORLD_SCALE = (vec3 310.0)
+let MAX_WORLD_LOD = 8.0
 
 let CUBE_SIZE = (uvec3 128)
 let GROUP_SIZE = 4

          
@@ 243,9 245,10 @@ fn matmapf (p)
     sdmDist
         do
             S := ((length p) - 300.0)
-            p := p * 0.01
+            #p := p * 0.01
+            p := p * 0.2
             local d = 0.0
-            for i in (range 8)
+            for i in (range 1)
                 s := (exp2 (i as f32))
                 d += ((triquad-noise3 (p * s)) * 2.0 - 1.0) / s
             max S

          
@@ 303,7 306,14 @@ fn map_onion_radius (p)
         k := 0.01
         r := r / (k * (1.0 - r))
 
-    _ ((unpack_normal_snorm p.xy) * r) (r * 2.5)
+    p := (unpack_normal_snorm p.xy)
+    #p :=
+        spherical-surface
+            L1-spherical
+                octahedral-surface
+                    tile-guyou ((p.xy - 1.0) / 2.0)
+
+    _ (p * r) r
 
 fn map_onion (p)
     let p r = (map_onion_radius p)

          
@@ 373,17 383,21 @@ fn generate-world ()
     imageStore world-out ipos (vec4 density 0 0 0)
     ;
 
-inline mapf (p)
-    d := (textureLod smp-world ((p / WORLD_SCALE) * 0.5 + 0.5) 0) . r
-    1.0 - d * 2.0
+inline mapf (p lod)
+    d := (textureLod smp-world ((p / WORLD_SCALE) * 0.5 + 0.5) lod) . r
+    0.5 - d
     #d
-inline matmapf (p)
-    sdmDist (mapf p)
+inline matmapf (p lod)
+    sdmDist (mapf p lod)
         sdMaterial
             vec4 0.5 0.3 1.0 1.0
-fn normalmapf (p r)
-    r := (0.25 / WORLD_SIZE.x) * WORLD_SCALE
-    - (sdNormalFast mapf p r)
+fn normalmapf (p lod)
+    r := (1.0 / WORLD_SIZE.x) * WORLD_SCALE
+    -
+        sdNormalFast
+            inline (p)
+                mapf p lod
+            \ p r
 
 fn generate-verts ()
     local_size GROUP_SIZE GROUP_SIZE GROUP_SIZE

          
@@ 408,21 422,28 @@ fn generate-verts ()
     flipmask := (? flipped? 1 0)
     flipsign := (? flipped? -1.0 1.0)
 
-    local cp : (array vec3 8)
+    inline ztolod (z)
+        z := z / WORLD_SCALE.x #* WORLD_SIZE.x
+        clamp (log2 z) 0.0 MAX_WORLD_LOD
+
+    local cp : (array vec4 8)
     local cd : (array f32 8)
     local mind = inf
     local maxd = -inf
     for i in (range 8)
         k := i ^ flipmask
         pos := (coord + ((vec3 ((uvec3 k (k >> 1:u32) (k >> 2:u32)) & 1:u32)) * d))
-        tpos := (map_vertex pos)
-        d := (mapf (map-translation tpos))
+        let tpos tr = (map_vertex_rlimit pos)
+        lod := (ztolod tr)
+        d := (mapf (map-translation tpos) lod)
         cp @ i =
-            do
-                static-if PROJECT_FINAL_VERTEX
-                    static-if POST_TRANSFORM pos
-                    else tpos
-                else pos
+            vec4
+                do
+                    static-if PROJECT_FINAL_VERTEX
+                        static-if POST_TRANSFORM pos
+                        else tpos
+                    else pos
+                lod
         cd @ i = d * flipsign
         mind = (min mind d)
         maxd = (max maxd d)

          
@@ 440,7 461,7 @@ fn generate-verts ()
 
         let idxs = (ivec4 0 k1 7 k3)
 
-        local p : (array vec3 4)
+        local p : (array vec4 4)
         p @ 0 = cp @ 0
         p @ 1 = cp @ k1
         p @ 2 = cp @ 7

          
@@ 465,19 486,19 @@ fn generate-verts ()
 
         if (c == 1:u32)
             # generate triangle
-            entries @ (ofs + 0) . xyz = (tf i.x i.y)
-            entries @ (ofs + 1) . xyz = (tf i.x i.z)
-            entries @ (ofs + 2) . xyz = (tf i.x i.w)
+            entries @ (ofs + 0) = (tf i.x i.y)
+            entries @ (ofs + 1) = (tf i.x i.z)
+            entries @ (ofs + 2) = (tf i.x i.w)
         elseif (c == 2:u32)
             p0 := (tf i.x i.z)
             p2 := (tf i.y i.w)
             # generate quad
-            entries @ (ofs + 0) . xyz = p0
-            entries @ (ofs + 1) . xyz = (tf i.x i.w)
-            entries @ (ofs + 2) . xyz = p2
-            entries @ (ofs + 3) . xyz = p2
-            entries @ (ofs + 4) . xyz = (tf i.y i.z)
-            entries @ (ofs + 5) . xyz = p0
+            entries @ (ofs + 0) = p0
+            entries @ (ofs + 1) = (tf i.x i.w)
+            entries @ (ofs + 2) = p2
+            entries @ (ofs + 3) = p2
+            entries @ (ofs + 4) = (tf i.y i.z)
+            entries @ (ofs + 5) = p0
 
     ;
 

          
@@ 487,20 508,22 @@ inout albedo : vec4
 inout matdata : vec4
 fn rasterize-vert ()
     let vertex-index = ((deref gl_VertexID) as u32)
-    let coord = (vec3 (vertex-in.entries @ vertex-index . xyz))
+    let vin = (deref (vertex-in.entries @ vertex-index))
+    let coord = (vec3 vin.xyz)
+    let lod = (vin.w as f32)
 
     let tcoord =
         static-if PROJECT_FINAL_VERTEX (map-translation coord)
         else (map_vertex coord)
 
-    let dist = (matmapf tcoord)
+    let dist = (matmapf tcoord lod)
     let material =
         dist.material
     #let dist = dist0
     #let material =
         'mix dist0.material dist1.material l
     let n =
-        normalmapf tcoord 0.001 # (r * 0.5)
+        normalmapf tcoord lod # (r * 0.5)
 
     # rotate it a little
     #embed

          
@@ 537,6 560,7 @@ fn rasterize-vert ()
             vec4 coord 1.0
 
     normal.out =
+        #(viridis (lod / MAX_WORLD_LOD)) * 2.0 - 1.0
         do
             static-if VISUALIZE_IDS ((vec3hash (vertex-index as f32)) * 2.0 - 1.0)
             else n