Merge branch 'bordel_laurent' of git.inpt.fr:tocard-inc/enseeiht/projet-be into bordel_laurent

src/fvi.py

@@ -1,7 +1,7 @@
 import numpy as np
 
 
-def fast_voxel_intersect(start, end, origin, step, shape) -> tuple[list, list, list]:
+def fast_voxel_intersect(start_, end_, origin_, step_, shape_) -> tuple[list, list, list]:
     """Compute the voxels intersected by a line segment.
 
     Args:
@@ -17,14 +17,14 @@ def fast_voxel_intersect(start, end, origin, step, shape) -> tuple[list, list, l
     """
 
     # Convert to numpy arrays
-    start = np.asarray(start)
-    end = np.asarray(end)
-    origin = np.asarray(origin)
-    step = np.asarray(step)
+    start_ = np.asarray(start_)
+    end_ = np.asarray(end_)
+    origin_ = np.asarray(origin_)
+    step_ = np.asarray(step_)
 
     # Translate line segment to voxel grid
-    start = start - origin
-    end = end - origin
+    start = (start_ - origin_) / step_
+    end = (end_ - origin_) / step_
 
     # Initialize list of intersected voxels
     intersections = []
@@ -45,17 +45,18 @@ def fast_voxel_intersect(start, end, origin, step, shape) -> tuple[list, list, l
 
     # Initialize current position to start
     position = start.copy()
-    # print("position: ", position)
+
+    # Initialize early exit
+    is_in = False
 
     # Main loop
     while True:
         # Compute the distance to the next boundaries
-        next_boundaries = np.divide(position + step * direction_signs, step)
+        next_boundaries = position + direction_signs
         errored = np.abs(np.round(next_boundaries) - next_boundaries) < 1e-12
         next_boundaries[errored] = np.round(next_boundaries[errored])
-        distances = (
-            (1 - is_negative) * np.floor(next_boundaries) + is_negative * np.ceil(next_boundaries)
-        ) * step - position
+        distances = ((1 - is_negative) * np.floor(next_boundaries) +
+                     is_negative * np.ceil(next_boundaries)) - position
 
         # Determine the nearest boundary to be reached
         boundary_distances = np.abs(distances / direction)
@@ -69,29 +70,24 @@ def fast_voxel_intersect(start, end, origin, step, shape) -> tuple[list, list, l
 
         # Update position
         position = position + clothest_boundary_distance * direction
-        # print("position_update: ", position)
 
         # Correct position to be on boundary
-        position[clothest_boundary] = (
-            round(position[clothest_boundary] / step[clothest_boundary]) * step[clothest_boundary]
-        )
+        position[clothest_boundary] = round(position[clothest_boundary])
 
         # Get corresponding voxel
-        on_boundary = np.mod(position, step) == 0
-        voxel = np.floor_divide(position, step) * step - is_negative * on_boundary * step
+        on_boundary = np.mod(position, np.ones_like(position)) == 0
+        voxel = np.floor(position) - is_negative * on_boundary
 
         # Add voxel to list
-        idx = np.floor_divide(voxel, step).astype(int)
-        if np.any(idx < 0) or np.any(idx >= shape):
+        idx = np.floor(voxel).astype(int)
+        if np.any(idx < 0) or np.any(idx >= shape_):
+            if is_in:
+                break
             continue
-        # print(f"voxel: {voxel}, step: {step}, idx: {idx}")
-        intersections.append(position + origin)
-        voxels.append(voxel + origin)
+        intersections.append(position * step_ + origin_)
+        voxels.append(voxel * step_ + origin_)
         voxels_idx.append(idx)
-
-    # print(f"intersections: {intersections}")
-    # print(f"voxels: {voxels}")
-    # print(f"voxels_idx: {voxels_idx}")
+        is_in = True
 
     return intersections, voxels, voxels_idx
 
@@ -132,13 +128,15 @@ if __name__ == "__main__":
             )
 
         # Plot line segment
-        plt.plot([start[0], end[0]], [start[1], end[1]], "k-")
-        plt.plot(start[0], start[1], "go")
-        plt.plot(end[0], end[1], "ro")
+        plt.plot([start[0], end[0]], [start[1], end[1]], 'k-')
 
         # Plot intersection points
         for pos in positions:
-            plt.plot(pos[0], pos[1], "bo")
+            plt.plot(pos[0], pos[1], 'bo')
+        
+        # Plot start and end points
+        plt.plot(start[0], start[1], 'go')
+        plt.plot(end[0], end[1], 'ro')
 
         # Plot voxel grid
         plt.axis("equal")
@@ -157,8 +155,8 @@ if __name__ == "__main__":
             update_figure()
 
     # Define voxel grid
-    origin = np.array([-5.0, -5.0])
-    step = np.array([0.7, 0.7])
+    origin = np.array([-5., -5.])
+    step = np.array([1.4, 1.4])
     shape = (10, 10)
 
     # Define segment

src/levelset.py

@@ -0,0 +1,39 @@
+import mcubes
+import numpy as np
+import matplotlib.pyplot as plt
+import imageio.v2 as imageio
+import perlin_noise
+
+
+# X, Y, Z = np.mgrid[:100, :100, :100]
+# V = np.sqrt((X - 50)**2 + (Y - 50)**2 + (Z - 50)**2)
+
+# for i in range(5, 45):
+#     vertices, triangles = mcubes.marching_cubes(V, i)
+#     mcubes.export_obj(vertices, triangles, f"cube_{i}.obj")
+
+noise = perlin_noise.PerlinNoise(octaves=6, seed=1)
+
+X, Y = np.mgrid[:100, :100]
+V = [[10 * noise([x/100, y/100]) + np.sqrt((x-50)**2 + (y-50)**2) for y in range(100)] for x in range(100)]
+V = np.array(V)
+V = (V - V.min()) / (V.max() - V.min())
+
+frame_list = []
+for i in np.linspace(0.05, 0.55, 100):
+    plt.clf()
+    plt.subplot(1, 2, 1)
+    plt.imshow(V, cmap="gray")
+    plt.contour(V, [i], colors="r")
+    plt.plot([0, 0, 100, 100, 0], [0, 100, 100, 0, 0], "k-")
+    plt.axis('off')
+
+    plt.subplot(1, 2, 2)
+    plt.imshow(V > i, cmap="gray")
+    plt.plot([0, 0, 100, 100, 0], [0, 100, 100, 0, 0], "k-")
+    plt.axis('off')
+
+    plt.savefig(f"/tmp/frame.png")
+    frame_list.append(imageio.imread(f"/tmp/frame.png"))
+
+imageio.mimsave('picture.gif', frame_list + frame_list[::-1], fps=60)

src/main.py

@@ -12,8 +12,8 @@ from fvi import fast_voxel_intersect
 
 VOXEL_SIZE = 5e-3
 X_MIN, X_MAX = 0.7, 1.3
-Y_MIN, Y_MAX = -0.1, 0.1
-Z_MIN, Z_MAX = -0.1, 0.1
+Y_MIN, Y_MAX = -0.05, 0.05
+Z_MIN, Z_MAX = -0.05, 0.05
 
 nb_frame = 24
 
@@ -42,7 +42,7 @@ for k in range(nb_frame):
         matrices = pickle.load(file)
         proj_mat = matrices["P"]
         proj_mats.append(proj_mat)
-        position = matrices["RT"][:, 3]
+        position = -np.linalg.inv(matrices["RT"][:3,:3]) @ matrices["RT"][:3,3]
         positions.append(position)
 
     cam_points = proj_mat @ points.T
@@ -85,14 +85,7 @@ for k in range(nb_frame):
 #     cv2.imshow('Frame', frame)
 #     cv2.waitKey(0)
 
-
-voxel = np.zeros(
-    (
-        int((X_MAX - X_MIN) / VOXEL_SIZE + 1),
-        int((Y_MAX - Y_MIN) / VOXEL_SIZE + 1),
-        int((Z_MAX - Z_MIN) / VOXEL_SIZE + 1),
-    )
-)
+voxel = np.zeros((int((X_MAX-X_MIN)/VOXEL_SIZE), int((Y_MAX-Y_MIN)/VOXEL_SIZE), int((Z_MAX-Z_MIN)/VOXEL_SIZE)))
 idx = np.floor_divide(points[:, :3] - np.array([X_MIN, Y_MIN, Z_MIN]), VOXEL_SIZE).astype(int)
 voxel[idx[:, 0], idx[:, 1], idx[:, 2]] = 1
 
@@ -128,9 +121,11 @@ for idx in track(np.argwhere(border)):
 
         # si le rayon ne traverse aucun autre voxel entre le centre du voxel (idx) et le centre de la caméra
         _, _, voxels_intersected = fast_voxel_intersect(start, end, origin, step, shape)
-        print(f"ntm: {voxels_intersected}")
         voxels_intersected = np.array(voxels_intersected, dtype=np.int32)
-        visible = voxel[voxels_intersected].sum() == 0
+        if len(voxels_intersected) == 0:
+            visible = True
+        else:
+            visible = voxel[voxels_intersected[:, 0], voxels_intersected[:, 1], voxels_intersected[:, 2]].sum() == 0
 
         if visible:
             proj = proj_mats[i] @ np.array([start[0], start[1], start[2], 1.0])
@@ -140,8 +135,7 @@ for idx in track(np.argwhere(border)):
 
     # calcule écartype des valeurs
     std = np.std(values)
-    # print(std)
-    # print(values)
+    print(f"std: {std}, values: {values}, nb_values: {len(values)}")
 
     # changer le levelset en fonction de l'écartype
     if std < 2: