voxels

@@ -0,0 +1,4 @@

+*.swp

+*.log

+venv

+tags

@@ -4,7 +4,7 @@ run:

 	python main.py

 dbg:

-	python -m pdb -c 'break debug_mode' main.py

+	python -m pdb -c 'break debug_mode ;; cont' main.py

 prof:

 	python -m cProfile main.py >$(PROF)

@@ -1,26 +1,12 @@

 from libs import *

-# https://iquilezles.org/articles/intersectors/

-def intersect_tri(ro, rd, tv):

-    e1 = tv[1] - tv[0]

-    e2 = tv[2] - tv[0]

-    to = ro - tv[0]

-    n = np.cross(e1, e2)

-    q = np.cross(to, rd)

-    d = 1.0 / np.dot(rd, n)

-    u = d * np.dot(-q, e2)

-    v = d * np.dot(q, e1)

-    t = d * np.dot(-n, to)

-    if u > 0.0 and v > 0.0 and u + v < 1.0 and t > 1e-3:

-        return t

-

 class Camera:

     def __init__(self, **kwargs):

         self.pos = vec([0, 0, 0])

         self.theta = 0

         self.phi = 0

         self.fovy = kwargs.get("fovy", 90)

-        self.aspect = 1

+        self.viewport = (640, 480)

     def looking(self):

         return vec([

@@ -32,6 +18,9 @@ class Camera:

     def up(self):

         return vec([0, 1, 0])

+    def aspect(self):

+        return self.viewport[0] / self.viewport[1]

+

     def view_mat(self):

         f = normalize(self.looking())

         s = normalize(np.cross(f, self.up()))

@@ -47,7 +36,7 @@ class Camera:

         f = 1.0 / np.tan(np.radians(self.fovy) / 2)

         depth = near - far

         return vec([

-            [f / self.aspect, 0, 0, 0],

+            [f / self.aspect(), 0, 0, 0],

             [0, f, 0, 0],

             [0, 0, (near + far) / depth, (2 * near * far) / depth],

             [0, 0, -1, 0]

@@ -56,7 +45,7 @@ class Camera:

     def mvp(self, ent):

         if ent.hud:

             am = np.identity(4, dtype=FLOAT)

-            am[0, 0] = 1 / self.aspect

+            am[0, 0] = 1 / self.aspect()

             return ent.translate() @ am @ ent.orient

         else:

             return self.persp_mat() @ self.view_mat() @ ent.translate() @ ent.orient

@@ -72,16 +61,7 @@ class Camera:

         l = self.looking()

         d, c, n = None, None, None

         for ent in scene:

-            if ent.hud:

-                continue

-            elif ent.draw_mode == GL_TRIANGLE_STRIP:

-                tl = ((i, ent.verts[i:i + 3]) for i in range(len(ent.verts) - 2))

-            elif ent.draw_mode == GL_TRIANGLES:

-                tl = ((i, ent.verts[3 * i:3 * (i + 1)]) for i in range(len(ent.verts) // 3))

-            else:

-                raise NotImplementedError

-            for (i, t) in tl:

-                r = intersect_tri(self.pos, l, np.dot(t, ent.orient[:3, :3].T) + ent.pos)

+            r, i = ent.raycast(self.pos, l)

             if r and (d is None or r < d):

                 d, c, n = r, ent, i

         return c, n

@@ -1,5 +1,19 @@

 from libs import *

+# https://iquilezles.org/articles/intersectors/

+def intersect_tri(ro, rd, tv):

+    e1 = tv[1] - tv[0]

+    e2 = tv[2] - tv[0]

+    to = ro - tv[0]

+    n = np.cross(e1, e2)

+    q = np.cross(to, rd)

+    d = 1.0 / np.dot(rd, n)

+    u = d * np.dot(-q, e2)

+    v = d * np.dot(q, e1)

+    t = d * np.dot(-n, to)

+    if u > 0.0 and v > 0.0 and u + v < 1.0 and t > 1e-3:

+        return t

+

 class Entity:

     def __init__(self, verts, col, pos, **kwargs):

         self.verts = vec(verts)

@@ -53,3 +67,19 @@ class Entity:

     def click(self, button, face):

         self.onclick(self, button, face)

+

+    def raycast(self, ro, rd):

+        if self.hud or self.draw_mode == GL_LINES:

+            return

+        elif self.draw_mode == GL_TRIANGLE_STRIP:

+            tl = ((i, self.verts[i:i + 3]) for i in range(len(self.verts) - 2))

+        elif self.draw_mode == GL_TRIANGLES:

+            tl = ((i, self.verts[3 * i:3 * (i + 1)]) for i in range(len(self.verts) // 3))

+        else:

+            raise NotImplementedError

+        d, n = None, None

+        for (i, t) in tl:

+            r = intersect_tri(ro, rd, np.dot(t, self.orient[:3, :3].T) + self.pos)

+            if r and (d is None or r < d):

+                d, n = r, i

+        return d, n

@@ -23,3 +23,8 @@ def roll_vec3(l, k):

             l[:] = [l[2], l[0], l[1]]

         case 2:

             l[:] = [l[1], l[2], l[0]]

+

+def axis2offset(a, x=1):

+    o = [0, 0, x * np.sign(a)]

+    roll_vec3(o, abs(a))

+    return o

@@ -1,6 +1,8 @@

 from OpenGL.GL import shaders

 from libs import *

 from entity import Entity

+from mapgen import MapGen

+from mesh import SquareMesh

 from player import Player

 from region import WorldRegion

@@ -24,28 +26,27 @@ void main() {

 }

 """

-FPS = 240

+FPS = 60

-def debug_mode():

+def debug_mode(scene, cam):

     pass

-def click_mesh(self, button, ind):

-    a, p = self.faces[ind]

+def click_mesh(self, button, coords):

     match button:

         case 1:

-            self.del_voxel(ind)

+            self.set_voxel(coords[1])

         case 2:

-            print(f"face {a} of {p}")

+            print(f"coords {coords}")

         case 3:

-            self.add_voxel(p + np.sign(a) * np.roll([0, 0, 1], abs(a)))

+            offset = axis2offset(coords[0])

+            self.set_voxel(coords[1] + offset, 2)

     if button % 2 == 1:

-        self.remesh()

         self.update_attr(0)

         self.update_attr(1)

-def build_scene():

+def build_scene(mg):

     box = WorldRegion()

-    print(f"map generated with {box.gen_perlin(np.random.rand(2, 2, 2, 3))} voxels")

+    print(f"map generated with {mg.world_natural(box)} voxels")

     box.remesh()

     box.onclick = click_mesh

     return [box]

@@ -58,13 +59,14 @@ def main():

         shaders.compileShader(VS, GL_VERTEX_SHADER),

         shaders.compileShader(FS, GL_FRAGMENT_SHADER)

     )

+    loc_mvp = glGetUniformLocation(sp, "MVP")

     glUseProgram(sp)

     glEnable(GL_DEPTH_TEST)

-    loc_mvp = glGetUniformLocation(sp, "MVP")

     glClearColor(0.1, 0.1, 0.1, 1.0)

-    scene = build_scene()

+    mg = MapGen()

+    scene = build_scene(mg)

     clock = pygame.time.Clock()

-    cam = Player()

+    cam = Player(speed=10/FPS)

     cam.pos = [8, 16, 8]

     cam.set_focus(True)

     while True:

@@ -73,16 +75,15 @@ def main():

                 pygame.quit()

                 return

             elif ev.type == pygame.KEYUP and ev.unicode == "d":

-                debug_mode()

+                debug_mode(scene, cam)

             cam.handle_event(ev, scene)

         cam.update_pos()

         glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

-        perf = clock.get_fps() / FPS

-        cam.hud["fps"].col = GREEN if 0.9 < perf < 1.1 else YELLOW if perf > 0.5 else RED

-        cam.hud["fps"].update_attr(1)

-        for ent in itertools.chain(scene, cam.hud.values()):

+        for ent in scene:

             cam.draw(ent, loc_mvp)

+        perf = clock.get_fps() / FPS

+        cam.draw_hud(loc_mvp, scene[0], perf)

         pygame.display.flip()

-        clock.tick(240)

+        clock.tick(FPS)

 main()

@@ -0,0 +1,57 @@

+from libs import *

+

+def smerp(ab, w):

+    assert all([abs(x) < 2 for x in ab])

+    for i in range(w.shape[0]):

+        ab = [a + w[i] ** 2 * (3 - 2 * w[i]) * (b - a) for (a, b) in zip(ab[::2], ab[1::2])]

+    return ab[0]

+

+class MapGen:

+    def __init__(self):

+        self.grads = [{} for _ in range(4)]

+

+    def perlin(self, name, pos, period):

+        dim = pos.shape[0]

+        gt = self.grads[dim].setdefault(name, {})

+        corner, fract = np.divmod(pos, period)

+        fract /= period

+        # https://stackoverflow.com/a/58827205

+        corners = np.arange(1 << dim)[:, None] >> np.arange(dim) & 1

+        dots = []

+        for c in corners:

+            ci = tuple(np.add(c, corner, dtype=int, casting="unsafe"))

+            dots.append(np.dot(fract - c, gt.setdefault(ci, np.random.rand(dim))))

+        return smerp(dots, fract)

+

+    def world_natural(self, wr):

+        n = 0

+        for i in range(wr.r):

+            for j in range(wr.r):

+                h = 8 + 4 * self.perlin("height", np.add(wr.pos[[0, 2]], (i, j)), np.array([10, 5]))

+                for k in range(int(h)):

+                    wr.grid[i, k, j] = 1

+                    n += 1

+        return n

+

+    def world_min(self, wr):

+        wr.add_voxel((0, 0, 0))

+        return 1

+

+    def world_grid(self, wr):

+        n = 0

+        for i in range(0, wr.r, 2):

+            for j in range(0, wr.r, 2):

+                for k in range(0, wr.r, 2):

+                    wr.grid[i, j, k] = 1

+                    n += 1

+        return n

+

+    def world_checker(self, wr):

+        n = 0

+        for i in range(wr.r):

+            for j in range(wr.r):

+                for k in range(wr.r):

+                    if not (i + j + k) % 2:

+                        wr.grid[i, j, k] = 1

+                        n += 1

+        return n

@@ -1,21 +1,14 @@

 from libs import *

 import entity

-def cube_col(axis):

-    return vec([

-        [0.0, 0.0, 0.0, 1.0],

-        [1.0, 0.0, 0.0, 1.0],

-        [0.0, 1.0, 0.0, 1.0],

-        [0.0, 0.0, 1.0, 1.0]

-    ][abs(axis)]) + 0.4 * np.sign(axis, dtype=FLOAT)

-

 class SquareMesh(entity.Entity):

-    def __init__(self):

+    def __init__(self, **kwargs):

         super().__init__(

             np.empty((0, 3), dtype=FLOAT),

             [0.0, 1.0, 0.0, 1.0],

             [0.0, 0.0, 0.0],

-            mode=GL_TRIANGLES

+            mode=GL_TRIANGLES,

+            **kwargs

         )

         self.faces = []

@@ -35,5 +28,31 @@ class SquareMesh(entity.Entity):

         self.verts = np.delete(self.verts, np.s_[3 * ind:3 * (ind + 2)], 0)

         del self.faces[ind:ind + 2]

+    def cube_col(self, fd, axis, corner):

+        if not isinstance(fd, int):

+            fd = self.grid[*fd]

+        match fd:

+            case 0:

+                return vec([0.0, 0.0, 0.0, 0.0])

+            case 1:

+                return vec([*np.add(0.5, axis2offset(axis, 0.5)), 1.0])

+            case 2:

+                top = vec([0.5, 0.5, 1.0, 1.0])

+                bottom = vec([0.0, 0.0, 0.5, 1.0])

+                match axis:

+                    case 3 | -3:

+                        return top if corner < 2 else bottom

+                    case 1 | -1:

+                        return top if not corner % 2 else bottom

+                    case 2:

+                        return top

+                    case -2:

+                        return bottom

+

     def col_list(self):

-        return np.repeat([cube_col(x[0]) for x in self.faces], 3, 0)

+        cl = np.empty((3 * len(self.faces), 4), dtype=FLOAT)

+        for s in range(len(self.faces) // 2):

+            axis, fd = self.faces[2 * s]

+            for (t, c) in enumerate([0, 1, 2, 3, 2, 1] if axis > 0 else [1, 2, 3, 2, 1, 0]):

+                cl[6 * s + t] = self.cube_col(fd, axis, c)

+        return cl

@@ -1,12 +1,13 @@

 from libs import *

 import camera

 import entity

+import mesh

-ARROWS = [pygame.K_UP, pygame.K_DOWN, pygame.K_LEFT, pygame.K_RIGHT, pygame.K_SPACE, pygame.K_LSHIFT]

+ARROWS = [pygame.K_RIGHT, pygame.K_LEFT, pygame.K_SPACE, pygame.K_LSHIFT, pygame.K_UP, pygame.K_DOWN]

 def reshape(w, h):

     glViewport(0, 0, w, h)

-    return w / h

+    return (w, h)

 class Player(camera.Camera):

     def __init__(self, **kwargs):

@@ -15,6 +16,8 @@ class Player(camera.Camera):

         self.rv = vec([0, 0, 0])

         self.speed = kwargs.get("speed", 0.05)

         self.hud = {

+            "skybox": mesh.SquareMesh(),

+            "shl": mesh.SquareMesh(),

             "crosshair": entity.Entity(

                 [[-0.03, 0.0, 0.0], [0.03, 0.0, 0.0], [0.0, -0.03, 0.0], [0.0, 0.03, 0.0]],

                 [1.0, 1.0, 1.0, 1.0], [0.0, 0.0, 0.0], mode=GL_LINES, hud=True

@@ -24,7 +27,14 @@ class Player(camera.Camera):

                 YELLOW, [-1.0, 1.0, 0.0], mode=GL_TRIANGLE_STRIP, hud=True

             ),

         }

-        self.viewport = (640, 480)

+        offset = [0, 0, 20]

+        for d in range(3):

+            roll_vec3(offset, 1)

+            for s in range(2):

+                self.hud["skybox"].add_sq((d + 1) * (-1) ** s, 40, (-1) ** s * vec(offset), 2)

+        self.hud["shl"].draw_mode = GL_LINE_STRIP

+        self.hud["shl"].add_sq(2, 0.96, [0, 0, 0], 2)

+        self.hud["shl"].verts[:] = self.hud["shl"].verts[[0, 1, 3, 2, 0, 3]]

     def set_focus(self, f): 

         pygame.event.set_grab(f)

@@ -40,11 +50,7 @@ class Player(camera.Camera):

         self.phi = min(np.pi / 2, max(-np.pi / 2, self.phi + self.sensitivity * dy))

     def update_motion(self, pk):

-        self.rv = vec([

-            pk[pygame.K_RIGHT] - pk[pygame.K_LEFT],

-            pk[pygame.K_SPACE] - pk[pygame.K_LSHIFT],

-            pk[pygame.K_UP] - pk[pygame.K_DOWN]

-        ])

+        self.rv = vec([pk[a] - pk[b] for (a, b) in zip(ARROWS[::2], ARROWS[1::2])])

     def vdof(self):

         l = self.looking()

@@ -62,17 +68,26 @@ class Player(camera.Camera):

                     self.set_focus(False)

                 elif ev.key in ARROWS:

                     self.update_motion(pygame.key.get_pressed())

-            case pygame.MOUSEMOTION:

-                if self.get_focus():

+            case pygame.MOUSEMOTION if self.get_focus():

                 self.mouse_move(ev.rel[0], -ev.rel[1])

-            case pygame.MOUSEBUTTONUP:

-                if ev.button == 1:

+            case pygame.MOUSEBUTTONUP if ev.button == 1:

                 self.set_focus(True)

-            case pygame.MOUSEBUTTONDOWN:

-                if self.get_focus():

+            case pygame.MOUSEBUTTONDOWN if self.get_focus():

                 t, n = self.target(scene)

                 if t:

                     t.click(ev.button, n)

             case pygame.WINDOWSIZECHANGED | pygame.WINDOWRESIZED:

-                self.aspect = reshape(ev.x, ev.y)

-                self.viewport = (ev.x, ev.y)

+                self.viewport = reshape(ev.x, ev.y)

+

+    def draw_hud(self, loc_mvp, wr, perf):

+        self.hud["skybox"].pos = self.pos

+        c, ap = self.target([wr])

+        if c:

+            self.hud["shl"].pos = ap[1] + axis2offset(ap[0], 0.52)

+            self.hud["shl"].orient[:3, :3] = np.sign(ap[0]) * np.roll(np.identity(3), abs(ap[0]) + 1, 0)

+        else:

+            self.hud["shl"].pos = vec([0, 0, 0])

+        self.hud["fps"].col = GREEN if 0.9 < perf < 1.1 else YELLOW if perf > 0.5 else RED

+        self.hud["fps"].update_attr(1)

+        for ent in self.hud.values():

+            self.draw(ent, loc_mvp)

@@ -1,23 +1,42 @@

 from libs import *

+import heapq

 import mesh

-def smoothstep(x):

-    return x ** 2 * (3 - 2 * x)

-

-def interp(ab, t):

-    return ab[0] + smoothstep(t) * (ab[1] - ab[0])

+def intersect_box(ro, rd, c, s):

+    m = 1 / rd

+    n = m * (ro - c)

+    k = np.abs(m) * s

+    tn = np.max(-n - k)

+    tf = np.min(-n + k)

+    if tf > tn and tf > 0:

+        return tn, tf

+    else:

+        return None, None

 class WorldRegion(mesh.SquareMesh):

     def __init__(self, size=16):

         super().__init__()

         self.grid = np.zeros((size, size, size), dtype=VOXEL_TYPE)

         self.r = size

+        self.dirty = 0

-    def add_voxel(self, pos, vt=1):

+    def set_voxel(self, pos, vt=0):

         self.grid[*pos] = vt

-

-    def del_voxel(self, ind):

-        self.grid[*self.faces[ind][1]] = 0

+        bp = bool(vt)

+        for d in range(3):

+            for s in range(2):

+                if pos[d] == self.r - 1 and s == 0 or pos[d] == 0 and s == 1:

+                    continue

+                a = (d + 1) * (-1) ** s

+                o = np.array(axis2offset(a), dtype=np.int64)

+                bc = bool(self.grid[*(pos + o)])

+                if bp ^ bc:

+                    a = (bp - bc) * (d + 1)

+                    self.add_sq(a, 1, pos + o / 2, pos + (0 if bp else o))

+        self.dirty += 1

+        if self.dirty > self.r:

+            self.dirty = 0

+            self.remesh()

     def remesh(self):

         offset = vec([0.5, 0.0, 0.0])

@@ -25,8 +44,7 @@ class WorldRegion(mesh.SquareMesh):

         prev, curr = 0, 0

         vl = []

         for d in range(3):

-            offset = np.array([1, 0, 0])

-            roll_vec3(offset, d)

+            offset = np.array(axis2offset(d + 1))

             for i in range(self.r):

                 for j in range(self.r):

                     for k in range(self.r + 1):

@@ -50,37 +68,19 @@ class WorldRegion(mesh.SquareMesh):

                                 self.faces.append(((bp - bc) * (d + 1), p))

         self.verts = vec(vl)

-    def prof_min(self):

-        self.add_voxel((0, 0, 0))

-        self.remesh()

-

-    def prof_med(self):

-        for i in range(self.r // 2):

-            for j in range(self.r // 2):

-                for k in range(self.r // 2):

-                    self.add_voxel((2 * i, 2 * j, 2 * k))

-        self.remesh()

-

-    def prof_max(self):

-        for i in range(self.r):

-            for j in range(self.r):

-                for k in range(self.r):

-                    if not (i + j + k) % 2:

-                        self.add_voxel((i, j, k))

-        self.remesh()

-

-    def gen_perlin(self, grads):

-        corners = np.array([(i, j, k) for i in range(2) for j in range(2) for k in range(2)])

-        n = 0

-        for i in range(self.r):

-            for j in range(self.r):

-                for k in range(self.r):

-                    w = np.array((i, j, k)) / self.r

-                    d = [np.dot(grads[*c], w) for c in corners]

-                    d = [interp(d[2 * k:2 * (k + 1)], w[2]) for k in range(4)]

-                    d = [interp(d[:2], w[1]), interp(d[2:4], w[1])]

-                    p = interp(d, w[0])

-                    if p > 2 * w[1] - 0.5:

-                        self.add_voxel((i, j, k))

-                        n += 1

-        return n

+    def raycast(self, ro, rd):

+        vr = np.repeat(self.r, 3)

+        tn, tf = intersect_box(ro, rd, self.pos + vr / 2 - 0.5, vr / 2)

+        if tn is not None:

+            if tn < 0:

+                tn = 0

+            ids = [[(d, i) for d in ((np.arange(self.r) + 0.5) * np.sign(rd[i]) - np.mod(ro[i] + tn * rd[i], 1)) / rd[i]] for i in range(3)]

+            for (te, i) in heapq.merge(*ids):

+                ts = tn + te

+                p = np.astype(np.round(ro + (ts + 1e-2) * rd - self.pos), np.int16)

+                if any(p < 0) or any(p >= self.r):

+                    break

+                elif self.grid[*p]:

+                    sa = np.astype((i + 1) * -np.sign(rd[i]), np.int8)

+                    return ts, (sa, p)

+        return None, None