le 2d preque réussi ?

src/borders.py

@@ -2,11 +2,9 @@ import numpy as np
 import matplotlib.pyplot as plt
 
 
-def update_border(voxel_values, idx=None):
+def update_border(voxel_values):
 
-    voxel_values = voxel_values > 0
+    if len(voxel_values.shape) == 3:
-
-    if idx is None:
         x_m1 = voxel_values[1:, :, :]
         x_m1 = np.concatenate((x_m1, np.zeros((1, x_m1.shape[1], x_m1.shape[2]))), axis=0)
 
@@ -36,7 +34,27 @@ def update_border(voxel_values, idx=None):
             )
         )
 
-    # TODO: update only concidered voxels (idx)
+    elif len(voxel_values.shape) == 2:
+        x_m1 = voxel_values[1:, :]
+        x_m1 = np.concatenate((x_m1, np.zeros((1, x_m1.shape[1]))), axis=0)
+
+        x_p1 = voxel_values[:-1, :]
+        x_p1 = np.concatenate((np.zeros((1, x_p1.shape[1])), x_p1), axis=0)
+
+        y_m1 = voxel_values[:, 1:]
+        y_m1 = np.concatenate((y_m1, np.zeros((y_m1.shape[0], 1))), axis=1)
+
+        y_p1 = voxel_values[:, :-1]
+        y_p1 = np.concatenate((np.zeros((y_p1.shape[0], 1)), y_p1), axis=1)
+
+        return np.logical_or.reduce(
+            (
+                voxel_values != x_m1,
+                voxel_values != x_p1,
+                voxel_values != y_m1,
+                voxel_values != y_p1,
+            )
+        )
 
 
 if __name__ == "__main__":

src/draw.py

@@ -1,5 +1,6 @@
 import matplotlib.pyplot as plt
 import numpy as np
+import cv2
 
 
 def f(x, y):
@@ -17,6 +18,7 @@ Z = (Z > 0.4).astype(np.float32)
 Z *= np.random.rand(*Z.shape)
 
 
+start = np.array([-3.75, -0.25])
 for i, x in enumerate(x_vals):
     for j, y in enumerate(y_vals):
         color = f"{hex(int(Z[j, i] * 255))[2:]}"
@@ -25,7 +27,7 @@ for i, x in enumerate(x_vals):
         else:
             color = "#" + 3 * color
         plt.fill([x, x + 0.1, x + 0.1, x], [y, y, y + 0.1, y + 0.1], color=color)
-
+# plt.plot(start[0], start[1], "ro")
 
 nb_cams = 32
 cam_poses = np.array(
@@ -50,10 +52,11 @@ for i in range(nb_cams):
 plt.xlim(-7, 7)
 plt.ylim(-7, 7)
 plt.axis("equal")
+plt.show()
 
 
 # draw 1d image of the scene for each camera
-for i in [0, 8, 16, 24]:
+for i in range(13, 21):
     plt.figure(f"cam {i}")
     # sort pixels by distance to camera
     cam_pose = cam_poses[i]
@@ -84,7 +87,24 @@ for i in [0, 8, 16, 24]:
     plt.xlim(-1, 1)
     plt.ylim(0, 0.2)
     plt.axis("off")
-    plt.savefig(f"data/peanut/images/Image{i:04}.png", dpi=300, bbox_inches="tight", pad_inches=0, transparent=True)
+    plt.savefig(f"/tmp/Image{i:04}.png", dpi=300, bbox_inches="tight", pad_inches=0, transparent=True)
+    # plt.close()
 
 
-plt.show()
+    RT = np.linalg.inv(cam2world_projs[i])[:-1, :]
+    proj = RT @ np.array([start[0], start[1], 1.0])
+    proj /= proj[0]
+    plt.plot(proj[1], 0.1, "ro")
+
+
+    plt.figure()
+    frame = cv2.imread(f"/tmp/Image{i:04}.png")
+    proj += np.array([0, 1.0])
+    proj *= 0.5 * np.array([1, frame.shape[1]])
+    proj = np.round(proj).astype(np.int32)
+    plt.imshow(frame, cmap="gray")
+    plt.plot(proj[1], frame.shape[0] / 2, "ro")
+
+    plt.show()
+
+# plt.show()

src/sfs2d.py

@@ -0,0 +1,224 @@
+import numpy as np
+import cv2
+import matplotlib.pyplot as plt
+from itertools import product
+from rich.progress import track
+
+from borders import update_border
+from fvi import fast_voxel_intersect
+
+
+
+x_vals = np.linspace(-5, 5, 100)
+y_vals = np.linspace(-5, 5, 100)
+
+VOXEL_SIZE = 0.1
+X_MIN, X_MAX = -5, 5
+Y_MIN, Y_MAX = -5, 5
+
+def plot_camera(cam2world_proj, name, color):
+    points = np.array([[0, 0, 1], [0.5, -0.2, 1], [0.5, 0.2, 1], [0, 0, 1]]).T
+    points = cam2world_proj @ points
+    plt.plot(points[0, 0], points[1, 0], color=color, marker="o")
+    plt.plot(points[0, :], points[1, :], color=color, linestyle="-")
+    plt.text(points[0, 0], points[1, 0], s=name, color=color)
+
+
+def plot_voxels(voxels, color, alpha, background):
+    if background:
+        plt.fill([X_MIN, X_MAX, X_MAX, X_MIN], [Y_MIN, Y_MIN, Y_MAX, Y_MAX], color="#000")
+    for i, j in np.argwhere(voxels):
+        x, y = X_MIN + i * VOXEL_SIZE, Y_MIN + j * VOXEL_SIZE
+        plt.fill([x, x + VOXEL_SIZE, x + VOXEL_SIZE, x], [y, y, y + VOXEL_SIZE, y + VOXEL_SIZE], color=color, alpha=alpha)
+
+
+
+
+nb_frame = 32
+
+points = np.array(
+    [
+        [x, y, 1.0]
+        for x, y in product(
+            np.arange(X_MIN, X_MAX, VOXEL_SIZE),
+            np.arange(Y_MIN, Y_MAX, VOXEL_SIZE),
+        )
+    ]
+)
+
+cam_poses = np.array(
+    [[6 * np.cos(theta), 6 * np.sin(theta)]
+     for theta in np.linspace(0, 2 * np.pi, nb_frame, endpoint=False)]
+)
+cam_rots = np.linspace(np.pi, 3 * np.pi, nb_frame, endpoint=False)
+cam2world_projs = np.array(
+    [
+        [[np.cos(theta), -np.sin(theta), cam_pose[0]],
+         [np.sin(theta), np.cos(theta), cam_pose[1]], [0, 0, 1]]
+        for theta, cam_pose in zip(cam_rots, cam_poses)
+    ]
+)
+
+mask = 255
+
+positions = []
+proj_mats = []
+frames = []
+
+for k in range(nb_frame):
+    frame = cv2.imread(
+        f"data/peanut/masks/Image{k:04}.png", cv2.IMREAD_GRAYSCALE)[0, :]
+    frames.append(cv2.imread(f"data/peanut/images/Image{k:04}.png", cv2.IMREAD_GRAYSCALE)[0, :])
+    # print(f"frame.shape: {frame.shape}")
+
+    RT = np.linalg.inv(cam2world_projs[k])[:-1, :]
+
+    # a = np.array([0.0, 0.0, 1.0])
+    cam_points = RT @ points.T
+    cam_points /= cam_points[0]
+    cam_points += np.array([[0], [1.0]])
+    cam_points *= 0.5 * np.array([[1], [frame.shape[0]]])
+    cam_points = np.round(cam_points).astype(np.int32)
+    # print(f"cam_points: {cam_points}")
+
+    visible = np.logical_and.reduce(
+        (
+            0 <= cam_points[1, :],
+            cam_points[1, :] < frame.shape[0],
+        )
+    )
+    cam_points = cam_points[:, visible]
+    points = points[visible, :]
+
+    solid = frame[cam_points[1, :]] == mask
+    cam_points = cam_points[:, solid]
+    points = points[solid, :]
+
+voxel = np.zeros((int((X_MAX-X_MIN)/VOXEL_SIZE), int((Y_MAX-Y_MIN)/VOXEL_SIZE)))
+idx = np.floor_divide(points[:, :2] - np.array([X_MIN, Y_MIN]), VOXEL_SIZE).astype(int)
+voxel[idx[:, 0], idx[:, 1]] = 1
+
+
+def f(x, y):
+    return np.exp(-((x**2) + y**2) / 3)
+
+
+x_vals = np.linspace(-5, 5, 100)
+y_vals = np.linspace(-5, 5, 100)
+
+
+X, Y = np.meshgrid(x_vals, y_vals)
+
+Z = f(X - 2, Y) + f(X + 2, Y)
+Z = (Z > 0.4).astype(np.float32)
+Z = np.swapaxes(Z, 0, 1)
+
+plt.figure()
+plt.axis("equal")
+plot_voxels(Z, "#fff", 1.0, True)
+plot_voxels(voxel, "#00f", 0.5, False)
+for i in range(nb_frame):
+    plot_camera(cam2world_projs[i], str(i), "#f00")
+
+
+origin = np.array([X_MIN, Y_MIN])
+step = np.array([VOXEL_SIZE, VOXEL_SIZE])
+shape = np.array([int((X_MAX-X_MIN)/VOXEL_SIZE), int((Y_MAX-Y_MIN)/VOXEL_SIZE)])
+
+
+for _ in range(5):
+    border = update_border(voxel)
+
+    plt.figure()
+    plt.axis("equal")
+    plot_voxels(voxel, "#fff", 1.0, True)
+    # plot_voxels(border, "#00f", 0.5, False)
+    for i in range(nb_frame):
+        plot_camera(cam2world_projs[i], str(i), "#f00")
+
+    for i, j in track(np.argwhere(border)):
+        x, y = X_MIN + i * VOXEL_SIZE, Y_MIN + j * VOXEL_SIZE
+        start = np.array([x, y]) + 0.5 * step
+
+        # plt.fill([x, x + VOXEL_SIZE, x + VOXEL_SIZE, x], [y, y, y + VOXEL_SIZE, y + VOXEL_SIZE], "#0f0", alpha=0.7)
+
+        values = []
+
+        for k in range(nb_frame):
+
+            # plot_camera(cam2world_projs[k], str(k), "#0f0")
+
+            end = cam_poses[k]
+            _, _, voxels_intersected = fast_voxel_intersect(start, end, origin, step, shape)
+            voxels_intersected = np.array(voxels_intersected, dtype=np.int32)
+
+            # for voxel_id in voxels_intersected:
+            #     if not voxel[voxel_id[0], voxel_id[1]]:
+            #         continue
+            #     plt.fill(
+            #         [
+            #             origin[0] + voxel_id[0] * step[0],
+            #             origin[0] + (voxel_id[0] + 1) * step[0],
+            #             origin[0] + (voxel_id[0] + 1) * step[0],
+            #             origin[0] + voxel_id[0] * step[0],
+            #         ],
+            #         [
+            #             origin[1] + voxel_id[1] * step[1],
+            #             origin[1] + voxel_id[1] * step[1],
+            #             origin[1] + (voxel_id[1] + 1) * step[1],
+            #             origin[1] + (voxel_id[1] + 1) * step[1],
+            #         ],
+            #         color="#f00",
+            #         alpha=0.3,
+            #     )
+
+            if len(voxels_intersected) == 0:
+                visible = True
+            else:
+                visible = voxel[voxels_intersected[:, 0], voxels_intersected[:, 1]].sum() == 0
+
+            if visible:
+                RT = np.linalg.inv(cam2world_projs[k])[:-1, :]
+                proj = RT @ np.array([start[0], start[1], 1.0])
+                proj /= proj[0]
+                proj += np.array([0, 1.0])
+                proj *= 0.5 * np.array([1, frames[k].shape[0]])
+                proj = np.round(proj).astype(np.int32)
+                values.append(frames[k][proj[1]])
+
+                # plt.figure()
+                # plt.imshow(np.repeat(frames[k][np.newaxis, :], 100, axis=0), cmap="gray")
+                # plt.plot(proj[1], 50, "ro")
+                # plt.show()
+            # else:
+            #     plot_camera(cam2world_projs[k], str(k), "#f00")
+
+        if np.std(values) < 2.0:
+            voxel[i, j] = 1
+            color = "#0f0"
+        else:
+            voxel[i, j] = 0
+            color = "#f00"
+        plt.fill(
+            [
+                origin[0] + i * step[0],
+                origin[0] + (i + 1) * step[0],
+                origin[0] + (i + 1) * step[0],
+                origin[0] + i * step[0],
+            ],
+            [
+                origin[1] + j * step[1],
+                origin[1] + j * step[1],
+                origin[1] + (j + 1) * step[1],
+                origin[1] + (j + 1) * step[1],
+            ],
+            color=color,
+            alpha=0.3,
+        )
+
+    plt.figure()
+    plt.axis("equal")
+    plot_voxels(Z, "#fff", 1.0, True)
+    plot_voxels(voxel, "#00f", 0.5, False)
+
+    plt.show()