séance castatrof

src/fvi.py

@@ -45,6 +45,7 @@ def fast_voxel_intersect(start, end, origin, step, shape) -> tuple[list, list, l
 
     # Initialize current position to start
     position = start.copy()
+    # print("position: ", position)
 
     # Main loop
     while True:
@@ -67,6 +68,7 @@ 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(
@@ -74,16 +76,21 @@ def fast_voxel_intersect(start, end, origin, step, shape) -> tuple[list, list, l
 
         # Get corresponding voxel
         on_boundary = np.mod(position, step) == 0
-        voxel = np.floor(position) - is_negative * on_boundary * step
+        voxel = np.floor_divide(position, step) * step - is_negative * on_boundary * step
 
         # Add voxel to list
         idx = np.floor_divide(voxel, step).astype(int)
         if np.any(idx < 0) or np.any(idx >= shape):
             continue
+        # print(f"voxel: {voxel}, step: {step}, idx: {idx}")
         intersections.append(position + origin)
         voxels.append(voxel + origin)
         voxels_idx.append(idx)
 
+    # print(f"intersections: {intersections}")
+    # print(f"voxels: {voxels}")
+    # print(f"voxels_idx: {voxels_idx}")
+
     return intersections, voxels, voxels_idx
 
 
@@ -137,7 +144,7 @@ if __name__ == '__main__':
 
     # Define voxel grid
     origin = np.array([-5., -5.])
-    step = np.array([1.0, 1.0])
+    step = np.array([0.7, 0.7])
     shape = (10, 10)
 
     # Define segment

src/intersec_line_voxel.py

@@ -9,7 +9,14 @@ def check_line_voxel(
     vx, vy, vz,
     c
 ):
-    """Check if a line intersects a voxel."""
+    """Check if a line intersects a voxel.
+    
+    Parameters:
+    - px, py, pz: voxel coner coordinates
+    - dx, dy, dz: line origin coordinates
+    - vx, vy, vz: line direction coordinates
+    - c: voxel size
+    """
 
     # Compute the intersection bounds
     kx1 = (px - dx) / vx

src/main.py

@@ -9,10 +9,10 @@ from rich.progress import track
 from borders import update_border
 from fvi import fast_voxel_intersect
 
-VOXEL_SIZE = 1e-1
-X_MIN, X_MAX = -2, 2
-Y_MIN, Y_MAX = -2, 2
-Z_MIN, Z_MAX = -1, 1
+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
 
 nb_frame = 24
 
@@ -21,15 +21,15 @@ points = np.array([[x, y, z, 1.0] for x, y, z in product(
     np.arange(Y_MIN, Y_MAX, VOXEL_SIZE),
     np.arange(Z_MIN, Z_MAX, VOXEL_SIZE))])
 
-mask = 255 * np.ones((3,))
+mask = 255
 
 positions = []
 proj_mats = []
 frames = []
 
 for k in range(nb_frame):
-    frame = cv2.imread(f'data/torus/masks/Image{k:04}.png')
-    frames.append(cv2.imread(f'data/torus/images/Image{k:04}.png'))
+    frame = cv2.imread(f'data/torus/masks/Image{k:04}.png', cv2.IMREAD_GRAYSCALE)
+    frames.append(cv2.imread(f'data/torus/images/Image{k:04}.png', cv2.IMREAD_GRAYSCALE))
 
     with open(f"data/torus/cameras/{k:04d}.pickle", 'rb') as file:
         matrices = pickle.load(file) 
@@ -46,7 +46,7 @@ for k in range(nb_frame):
     cam_points = cam_points[:,visible]
     points = points[visible,:]
 
-    solid = (frame[cam_points[1,:],cam_points[0,:]] == mask).all(axis=1)
+    solid = frame[cam_points[1,:],cam_points[0,:]] == mask
     cam_points = cam_points[:,solid]
     points = points[solid,:]
 
@@ -72,7 +72,7 @@ for k in range(nb_frame):
 #     cv2.waitKey(0)
 
 
-voxel = np.zeros((int((X_MAX-X_MIN)/VOXEL_SIZE), int((Y_MAX-Y_MIN)/VOXEL_SIZE), int((Z_MAX-Z_MIN)/VOXEL_SIZE)))
+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)))
 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
 
@@ -90,7 +90,6 @@ step = np.array([VOXEL_SIZE, VOXEL_SIZE, VOXEL_SIZE])
 shape = np.array([int((X_MAX-X_MIN)/VOXEL_SIZE), int((Y_MAX-Y_MIN)/VOXEL_SIZE), int((Z_MAX-Z_MIN)/VOXEL_SIZE)])
 
 for idx in track(np.argwhere(border)):
-    
     # coordonnées du centre du voxel
     start = np.array([
         X_MIN + (idx[0] + 0.5) * VOXEL_SIZE,
@@ -107,7 +106,12 @@ for idx in track(np.argwhere(border)):
         end = positions[i]
 
         # si le rayon ne traverse aucun autre voxel entre le centre du voxel (idx) et le centre de la caméra
-        if not fast_voxel_intersect(start, end, origin, step, shape):
+        _, _, 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 visible:
             proj = proj_mats[i] @ np.array([start[0], start[1], start[2], 1.0])
             proj /= proj[2]
             proj = np.round(proj).astype(np.int32)
@@ -115,6 +119,8 @@ for idx in track(np.argwhere(border)):
 
     # calcule écartype des valeurs
     std = np.std(values)
+    # print(std)
+    # print(values)
 
     # changer le levelset en fonction de l'écartype
     if std < 2:
@@ -124,5 +130,5 @@ for idx in track(np.argwhere(border)):
 
 
 
-vertices, triangles = mcubes.marching_cubes(border, 0)
+vertices, triangles = mcubes.marching_cubes(voxel, 0)
 mcubes.export_obj(vertices, triangles, "result.obj")