From cf1f2d76a91e8395289a8d4701ecc0f99a359abc Mon Sep 17 00:00:00 2001 From: gdamms Date: Tue, 31 Jan 2023 14:51:14 +0100 Subject: [PATCH] =?UTF-8?q?le=202d=20preque=20r=C3=A9ussi=20=3F?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/borders.py | 28 +++++-- src/draw.py | 28 ++++++- src/sfs2d.py | 224 +++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 271 insertions(+), 9 deletions(-) create mode 100644 src/sfs2d.py diff --git a/src/borders.py b/src/borders.py index ed4d32a..8df1ba6 100644 --- a/src/borders.py +++ b/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 idx is None: + if len(voxel_values.shape) == 3: 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__": diff --git a/src/draw.py b/src/draw.py index bcde56a..f115aed 100644 --- a/src/draw.py +++ b/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() diff --git a/src/sfs2d.py b/src/sfs2d.py new file mode 100644 index 0000000..9e1b083 --- /dev/null +++ b/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()