ENSEEIHT/projet-probleme-inverse-3D
1
0
Fork 0
Code Issues Pull requests Packages Projects Releases Wiki Activity

good enough

Browse source
This commit is contained in:
gdamms 2023-02-01 09:12:18 +01:00
parent f4e53a6e9e
commit 5cbb0d31ac
2 changed files with 124 additions and 61 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

51
src/draw.py
View file

@ -7,8 +7,8 @@ VOXEL_SIZE = 0.1
X_MIN, X_MAX = -5, 5 X_MIN, X_MAX = -5, 5
Y_MIN, Y_MAX = -5, 5 Y_MIN, Y_MAX = -5, 5
x_vals = np.arange(X_MIN, X_MAX, VOXEL_SIZE) x_vals = np.arange(X_MIN + VOXEL_SIZE / 2, X_MAX, VOXEL_SIZE)
y_vals = np.arange(Y_MIN, Y_MAX, VOXEL_SIZE) y_vals = np.arange(Y_MIN + VOXEL_SIZE / 2, Y_MAX, VOXEL_SIZE)
image_length = 1500 image_length = 1500
@ -23,15 +23,15 @@ Z = f(X - 2, Y) + f(X + 2, Y)
Z = (Z > 0.4).astype(np.float32) Z = (Z > 0.4).astype(np.float32)
Z *= np.random.rand(*Z.shape) Z *= np.random.rand(*Z.shape)
for i, x in enumerate(x_vals): for j, x in enumerate(x_vals):
for j, y in enumerate(y_vals): for i, y in enumerate(y_vals):
color = f"{hex(int(Z[j, i] * 255))[2:]}" color = f"{hex(int(Z[i, j] * 255))[2:]}"
if color == "0": if color == "0":
color = "#003" color = "#003"
else: else:
color = "#" + 3 * color color = "#" + 3 * color
plt.fill([x, x + 0.1, x + 0.1, x], plt.fill([x - VOXEL_SIZE/2, x + VOXEL_SIZE/2, x + VOXEL_SIZE/2, x - VOXEL_SIZE/2],
[y, y, y + 0.1, y + 0.1], color=color) [y - VOXEL_SIZE/2, y - VOXEL_SIZE/2, y + VOXEL_SIZE/2, y + VOXEL_SIZE/2], color=color)
nb_cams = 32 nb_cams = 32
cam_poses = np.array( cam_poses = np.array(
@ -49,8 +49,8 @@ cam2world_projs = np.array(
) )
for i in range(nb_cams): for i in range(nb_cams):
x = np.array([[0, 0, 1], [0.5, -0.2, 1], [0.5, 0.2, 1], [0, 0, 1]]).T x = np.array([[0, 0, 1], [0.5, -0.2, 1], [0.5, 0.2, 1], [0, 0, 1]])
x = cam2world_projs[i] @ x x = cam2world_projs[i] @ x.T
plt.plot(cam_poses[i][0], cam_poses[i][1], "ro") plt.plot(cam_poses[i][0], cam_poses[i][1], "ro")
plt.text(cam_poses[i][0], cam_poses[i][1], str(i)) plt.text(cam_poses[i][0], cam_poses[i][1], str(i))
plt.plot(x[0, :], x[1, :], "r-") plt.plot(x[0, :], x[1, :], "r-")
@ -79,16 +79,15 @@ for i in range(nb_cams):
for j in pixels_sort: for j in pixels_sort:
x, y = X.flatten()[j], Y.flatten()[j] x, y = X.flatten()[j], Y.flatten()[j]
color = f"{hex(int(Z.flatten()[j] * 255))[2:]}" color = int(Z.flatten()[j] * 255)
if color == "0": if color == 0:
continue continue
color = "#" + 3 * color
RT = np.linalg.inv(cam2world_projs[i])[:-1, :] RT = np.linalg.inv(cam2world_projs[i])[:-1, :]
px = np.array([[x, y, 1], px = np.array([[x - VOXEL_SIZE/2, y - VOXEL_SIZE/2, 1],
[x + 0.1, y, 1], [x + VOXEL_SIZE/2, y - VOXEL_SIZE/2, 1],
[x + 0.1, y + 0.1, 1], [x + VOXEL_SIZE/2, y + VOXEL_SIZE/2, 1],
[x, y + 0.1, 1]]) [x - VOXEL_SIZE/2, y + VOXEL_SIZE/2, 1]])
px = RT @ px.T px = RT @ px.T
px /= px[0, :] px /= px[0, :]
px += np.array([[0], [1.0]]) px += np.array([[0], [1.0]])
@ -98,8 +97,26 @@ for i in range(nb_cams):
x1 = px[1, :].max() x1 = px[1, :].max()
img[:, int(x0):int(x1), :] = np.array( img[:, int(x0):int(x1), :] = np.array(
3 * [int(Z.flatten()[j] * 255)], dtype=np.uint8) 3 * [color], dtype=np.uint8)
mask[:, int(x0):int(x1), :] = np.array([255, 255, 255], dtype=np.uint8) mask[:, int(x0):int(x1), :] = np.array([255, 255, 255], dtype=np.uint8)
# x, y = X_MIN + 27 * VOXEL_SIZE, Y_MIN + 66 * VOXEL_SIZE
# RT = np.linalg.inv(cam2world_projs[i])[:-1, :]
# px = np.array([[x + VOXEL_SIZE/2, y + VOXEL_SIZE/2, 1],
# [x, y, 1],
# [x + VOXEL_SIZE, y, 1],
# [x + VOXEL_SIZE, y + VOXEL_SIZE, 1],
# [x, y + VOXEL_SIZE, 1]])
# px = RT @ px.T
# px /= px[0, :]
# px += np.array([[0], [1.0]])
# px *= 0.5 * np.array([[1], [image_length]])
# img[48:52, int(px[1, 0]-1):int(px[1, 0]+1), :] = np.array([255, 0, 0], dtype=np.uint8)
# img[48:52, int(px[1, 1]-1):int(px[1, 1]+1), :] = np.array([255, 255, 0], dtype=np.uint8)
# img[48:52, int(px[1, 2]-1):int(px[1, 2]+1), :] = np.array([0, 255, 255], dtype=np.uint8)
# img[48:52, int(px[1, 3]-1):int(px[1, 3]+1), :] = np.array([255, 0, 255], dtype=np.uint8)
# img[48:52, int(px[1, 4]-1):int(px[1, 4]+1), :] = np.array([0, 255, 0], dtype=np.uint8)
cv2.imwrite(f"data/peanut/images/Image{i:04}.png", img) cv2.imwrite(f"data/peanut/images/Image{i:04}.png", img)
cv2.imwrite(f"data/peanut/masks/Image{i:04}.png", mask) cv2.imwrite(f"data/peanut/masks/Image{i:04}.png", mask)

132
src/sfs2d.py
View file

@ -8,14 +8,11 @@ from borders import update_border
from fvi import fast_voxel_intersect from fvi import fast_voxel_intersect
x_vals = np.linspace(-5, 5, 100)
y_vals = np.linspace(-5, 5, 100)
VOXEL_SIZE = 0.1 VOXEL_SIZE = 0.1
X_MIN, X_MAX = -5, 5 X_MIN, X_MAX = -5, 5
Y_MIN, Y_MAX = -5, 5 Y_MIN, Y_MAX = -5, 5
def plot_camera(cam2world_proj, name, color): 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 = np.array([[0, 0, 1], [0.5, -0.2, 1], [0.5, 0.2, 1], [0, 0, 1]]).T
points = cam2world_proj @ points points = cam2world_proj @ points
@ -26,12 +23,12 @@ def plot_camera(cam2world_proj, name, color):
def plot_voxels(voxels, color, alpha, background): def plot_voxels(voxels, color, alpha, background):
if background: if background:
plt.fill([X_MIN, X_MAX, X_MAX, X_MIN], [Y_MIN, Y_MIN, Y_MAX, Y_MAX], color="#000") 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): for i, j in np.argwhere(voxels):
x, y = X_MIN + i * VOXEL_SIZE, Y_MIN + j * VOXEL_SIZE 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) 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 nb_frame = 32
@ -40,8 +37,8 @@ points = np.array(
[ [
[x, y, 1.0] [x, y, 1.0]
for x, y in product( for x, y in product(
np.arange(X_MIN, X_MAX, VOXEL_SIZE), np.arange(X_MIN + VOXEL_SIZE / 2, X_MAX, VOXEL_SIZE),
np.arange(Y_MIN, Y_MAX, VOXEL_SIZE), np.arange(Y_MIN + VOXEL_SIZE / 2, Y_MAX, VOXEL_SIZE),
) )
] ]
) )
@ -54,7 +51,8 @@ cam_rots = np.linspace(np.pi, 3 * np.pi, nb_frame, endpoint=False)
cam2world_projs = np.array( cam2world_projs = np.array(
[ [
[[np.cos(theta), -np.sin(theta), cam_pose[0]], [[np.cos(theta), -np.sin(theta), cam_pose[0]],
[np.sin(theta), np.cos(theta), cam_pose[1]], [0, 0, 1]] [np.sin(theta), np.cos(theta), cam_pose[1]],
[0, 0, 1]]
for theta, cam_pose in zip(cam_rots, cam_poses) for theta, cam_pose in zip(cam_rots, cam_poses)
] ]
) )
@ -68,18 +66,16 @@ frames = []
for k in range(nb_frame): for k in range(nb_frame):
frame = cv2.imread( frame = cv2.imread(
f"data/peanut/masks/Image{k:04}.png", cv2.IMREAD_GRAYSCALE)[0, :] 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, :]) frames.append(cv2.imread(
# print(f"frame.shape: {frame.shape}") f"data/peanut/images/Image{k:04}.png", cv2.IMREAD_GRAYSCALE)[0, :])
RT = np.linalg.inv(cam2world_projs[k])[:-1, :] RT = np.linalg.inv(cam2world_projs[k])[:-1, :]
# a = np.array([0.0, 0.0, 1.0])
cam_points = RT @ points.T cam_points = RT @ points.T
cam_points /= cam_points[0] cam_points /= cam_points[0]
cam_points += np.array([[0], [1.0]]) cam_points += np.array([[0], [1.0]])
cam_points *= 0.5 * np.array([[1], [frame.shape[0]]]) cam_points *= 0.5 * np.array([[1], [frame.shape[0]]])
cam_points = np.round(cam_points).astype(np.int32) cam_points = np.round(cam_points).astype(np.int32)
# print(f"cam_points: {cam_points}")
visible = np.logical_and.reduce( visible = np.logical_and.reduce(
( (
@ -94,8 +90,10 @@ for k in range(nb_frame):
cam_points = cam_points[:, solid] cam_points = cam_points[:, solid]
points = points[solid, :] points = points[solid, :]
voxel = np.zeros((int((X_MAX-X_MIN)/VOXEL_SIZE), int((Y_MAX-Y_MIN)/VOXEL_SIZE))) voxel = np.zeros((int((X_MAX-X_MIN)/VOXEL_SIZE),
idx = np.floor_divide(points[:, :2] - np.array([X_MIN, Y_MIN]), VOXEL_SIZE).astype(int) 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 voxel[idx[:, 0], idx[:, 1]] = 1
@ -103,15 +101,14 @@ def f(x, y):
return np.exp(-((x**2) + y**2) / 3) return np.exp(-((x**2) + y**2) / 3)
x_vals = np.linspace(-5, 5, 100) x_vals = np.arange(X_MIN + VOXEL_SIZE / 2, X_MAX, VOXEL_SIZE)
y_vals = np.linspace(-5, 5, 100) y_vals = np.arange(Y_MIN + VOXEL_SIZE / 2, Y_MAX, VOXEL_SIZE)
X, Y = np.meshgrid(x_vals, y_vals) X, Y = np.meshgrid(x_vals, y_vals)
Z = f(X - 2, Y) + f(X + 2, Y) Z = f(X - 2, Y) + f(X + 2, Y)
Z = (Z > 0.4).astype(np.float32) Z = (Z > 0.4).astype(np.float32)
Z = np.swapaxes(Z, 0, 1) Z = Z.T
plt.figure() plt.figure()
plt.axis("equal") plt.axis("equal")
@ -119,37 +116,60 @@ plot_voxels(Z, "#fff", 1.0, True)
plot_voxels(voxel, "#00f", 0.5, False) plot_voxels(voxel, "#00f", 0.5, False)
for i in range(nb_frame): for i in range(nb_frame):
plot_camera(cam2world_projs[i], str(i), "#f00") plot_camera(cam2world_projs[i], str(i), "#f00")
# plt.show()
plt.savefig(f"init.png", dpi=300, bbox_inches="tight", pad_inches=0, transparent=True)
plt.close()
origin = np.array([X_MIN, Y_MIN]) origin = np.array([X_MIN, Y_MIN])
step = np.array([VOXEL_SIZE, VOXEL_SIZE]) step = np.array([VOXEL_SIZE, VOXEL_SIZE])
shape = np.array([int((X_MAX-X_MIN)/VOXEL_SIZE), int((Y_MAX-Y_MIN)/VOXEL_SIZE)]) shape = np.array([int((X_MAX-X_MIN)/VOXEL_SIZE),
int((Y_MAX-Y_MIN)/VOXEL_SIZE)])
# x, y = 2.45, -1.55
# i, j = np.floor_divide(np.array([x, y]) - origin, step).astype(int)
# print(i, j)
# exit()
for _ in range(5): for p in range(11):
border = update_border(voxel) border = update_border(voxel)
voxel_ = voxel.copy()
plt.figure() plt.figure()
plt.axis("equal") plt.axis("equal")
plot_voxels(voxel, "#fff", 1.0, True) plot_voxels(Z, "#fff", 1.0, True)
# plot_voxels(border, "#00f", 0.5, False) # plot_voxels(border, "#00f", 0.5, False)
for i in range(nb_frame): for i in range(nb_frame):
plot_camera(cam2world_projs[i], str(i), "#f00") plot_camera(cam2world_projs[i], str(i), "#f00")
# plt.show()
for i, j in track(np.argwhere(border)): for i, j in track(np.argwhere(border)):
x, y = X_MIN + i * VOXEL_SIZE, Y_MIN + j * VOXEL_SIZE # if i != 74 or j != 34:
start = np.array([x, y]) + 0.5 * step # continue
x, y = X_MIN + (i + 0.5) * VOXEL_SIZE, Y_MIN + (j + 0.5) * VOXEL_SIZE
# plt.fill([x, x + VOXEL_SIZE, x + VOXEL_SIZE, x], [y, y, y + VOXEL_SIZE, y + VOXEL_SIZE], "#0f0", alpha=0.7) start = np.array([x, y])
# plt.plot([start[0]], [start[1]], "o", color="#f00")
# plt.fill([x, x + VOXEL_SIZE, x + VOXEL_SIZE, x],
# [y, y, y + VOXEL_SIZE, y + VOXEL_SIZE], "#0f0", alpha=0.7)
values = [] values = []
# a b c d all four corners
a = np.array([x, y])
# plt.plot([a[0]], [a[1]], "o", color="#ff0")
b = np.array([x + VOXEL_SIZE, y])
# plt.plot([b[0]], [b[1]], "o", color="#0ff")
c = np.array([x + VOXEL_SIZE, y + VOXEL_SIZE])
# plt.plot([c[0]], [c[1]], "o", color="#f0f")
d = np.array([x, y + VOXEL_SIZE])
# plt.plot([d[0]], [d[1]], "o", color="#0f0")
for k in range(nb_frame): for k in range(nb_frame):
# plot_camera(cam2world_projs[k], str(k), "#0f0") # plot_camera(cam2world_projs[k], str(k), "#0f0")
end = cam_poses[k] end = cam_poses[k]
_, _, voxels_intersected = fast_voxel_intersect(start, end, origin, step, shape) _, _, voxels_intersected = fast_voxel_intersect(
start, end, origin, step, shape)
voxels_intersected = np.array(voxels_intersected, dtype=np.int32) voxels_intersected = np.array(voxels_intersected, dtype=np.int32)
# for voxel_id in voxels_intersected: # for voxel_id in voxels_intersected:
@ -175,29 +195,49 @@ for _ in range(5):
if len(voxels_intersected) == 0: if len(voxels_intersected) == 0:
visible = True visible = True
else: else:
visible = voxel[voxels_intersected[:, 0], voxels_intersected[:, 1]].sum() == 0 visible = voxel[voxels_intersected[:, 0],
voxels_intersected[:, 1]].sum() == 0
if visible: if visible:
RT = np.linalg.inv(cam2world_projs[k])[:-1, :] RT = np.linalg.inv(cam2world_projs[k])[:-1, :]
proj = RT @ np.array([start[0], start[1], 1.0]) # concatenate start with a b c and d
proj /= proj[0] points = np.concatenate([start[np.newaxis, :],
proj += np.array([0, 1.0]) a[np.newaxis, :],
proj *= 0.5 * np.array([1, frames[k].shape[0]]) b[np.newaxis, :],
c[np.newaxis, :],
d[np.newaxis, :]], axis=0)
# homogenous coordinates
points = np.concatenate([points, np.ones((5, 1))], axis=1)
proj = RT @ points.T
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) proj = np.round(proj).astype(np.int32)
values.append(frames[k][proj[1]]) values.append(frames[k][proj[1, 0]])
# fig = plt.gcf()
# plt.figure() # plt.figure()
# plt.imshow(np.repeat(frames[k][np.newaxis, :], 100, axis=0), cmap="gray") # plt.imshow(
# plt.plot(proj[1], 50, "ro") # np.repeat(frames[k][np.newaxis, :], 100, axis=0), cmap="gray")
# plt.show() # plt.plot(proj[1, 0], 50, "o", color="#f00")
# plt.plot(proj[1, 1], 50, "o", color="#ff0")
# plt.plot(proj[1, 2], 50, "o", color="#0ff")
# plt.plot(proj[1, 3], 50, "o", color="#f0f")
# plt.plot(proj[1, 4], 50, "o", color="#0f0")
# plt.figure(fig)
# else: # else:
# plot_camera(cam2world_projs[k], str(k), "#f00") # plot_camera(cam2world_projs[k], str(k), "#f00")
if np.std(values) < 2.0: # print(values)
voxel[i, j] = 1 # plt.show()
occurences = [[x, values.count(x)] for x in set(values)]
occurences = sorted(occurences, key=lambda x: x[1], reverse=True)
if occurences[0][1] >= len(values) * 0.68:
voxel_[i, j] = 1
color = "#0f0" color = "#0f0"
else: else:
voxel[i, j] = 0 voxel_[i, j] = 0
color = "#f00" color = "#f00"
plt.fill( plt.fill(
[ [
@ -216,9 +256,15 @@ for _ in range(5):
alpha=0.3, alpha=0.3,
) )
voxel = voxel_
plt.savefig(f"evol{p}.png", dpi=300, bbox_inches="tight", pad_inches=0, transparent=True)
plt.figure() plt.figure()
plt.axis("equal") plt.axis("equal")
plot_voxels(Z, "#fff", 1.0, True) plot_voxels(Z, "#fff", 1.0, True)
plot_voxels(voxel, "#00f", 0.5, False) plot_voxels(voxel, "#00f", 0.5, False)
plt.savefig(f"shape{p}.png", dpi=300, bbox_inches="tight", pad_inches=0, transparent=True)
plt.close("all")
plt.show() # plt.show()