diff --git a/src/fvi.py b/src/fvi.py index 7ecd277..8ff1d6d 100644 --- a/src/fvi.py +++ b/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 diff --git a/src/levelset.py b/src/levelset.py new file mode 100644 index 0000000..3345ef2 --- /dev/null +++ b/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) diff --git a/src/main.py b/src/main.py index 9801901..69b6490 100644 --- a/src/main.py +++ b/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: