diff --git a/apps/calibration/align_colmap_ground.py b/apps/calibration/align_colmap_ground.py new file mode 100644 index 0000000..2c54a57 --- /dev/null +++ b/apps/calibration/align_colmap_ground.py @@ -0,0 +1,232 @@ +# 这个脚本用于对colmap的相机标定结果，寻找地面与场景中心 +# 方法： +# 1. 使用棋盘格 +# 2. 估计点云里的地面 +import os +from os.path import join +from easymocap.annotator.file_utils import save_json +from easymocap.mytools.debug_utils import myerror, run_cmd, mywarn, log +from easymocap.mytools.camera_utils import read_cameras, write_camera +from easymocap.mytools import read_json +from easymocap.mytools import batch_triangulate, projectN3, Undistort +import numpy as np +import cv2 +from apps.calibration.calib_extri import solvePnP + +def guess_ground(pcdname): + pcd = o3d.io.read_point_cloud(pcdname) + +def compute_rel(R_src, T_src, R_tgt, T_tgt): + R_rel = R_src.T @ R_tgt + T_rel = R_src.T @ (T_tgt - T_src) + return R_rel, T_rel + +def triangulate(cameras, areas): + Ps, k2ds = [], [] + for cam, _, k2d, k3d in areas: + k2d = Undistort.points(k2d, cameras[cam]['K'], cameras[cam]['dist']) + P = cameras[cam]['K'] @ np.hstack([cameras[cam]['R'], cameras[cam]['T']]) + Ps.append(P) + k2ds.append(k2d) + Ps = np.stack(Ps) + k2ds = np.stack(k2ds) + k3d = batch_triangulate(k2ds, Ps) + return k3d + +def best_fit_transform(A, B): + ''' + Calculates the least-squares best-fit transform that maps corresponding points A to B in m spatial dimensions + Input: + A: Nxm numpy array of corresponding points + B: Nxm numpy array of corresponding points + Returns: + T: (m+1)x(m+1) homogeneous transformation matrix that maps A on to B + R: mxm rotation matrix + t: mx1 translation vector + ''' + + assert A.shape == B.shape + + # get number of dimensions + m = A.shape[1] + + # translate points to their centroids + centroid_A = np.mean(A, axis=0) + centroid_B = np.mean(B, axis=0) + AA = A - centroid_A + BB = B - centroid_B + + # rotation matrix + H = np.dot(AA.T, BB) + U, S, Vt = np.linalg.svd(H) + R = np.dot(Vt.T, U.T) + + # special reflection case + if np.linalg.det(R) < 0: + Vt[m-1,:] *= -1 + R = np.dot(Vt.T, U.T) + + # translation + t = centroid_B.T - np.dot(R,centroid_A.T) + + return R, t + +def align_by_chessboard(cameras, path): + camnames = sorted(os.listdir(join(path, 'chessboard'))) + areas = [] + for ic, cam in enumerate(camnames): + imagename = join(path, 'images', cam, '000000.jpg') + chessname = join(path, 'chessboard', cam, '000000.json') + data = read_json(chessname) + k3d = np.array(data['keypoints3d'], dtype=np.float32) + k2d = np.array(data['keypoints2d'], dtype=np.float32) + # TODO + # pattern = (11, 8) + if 'pattern' in data.keys(): + pattern = data['pattern'] + else: + pattern = None + img = cv2.imread(imagename) + if args.scale2d is not None: + k2d[:, :2] *= args.scale2d + img = cv2.resize(img, None, fx=args.scale2d, fy=args.scale2d) + if args.origin is not None: + cameras[args.prefix+cam] = cameras.pop(args.origin+cam.replace('VID_', '0000')) + cam = args.prefix + cam + if cam not in cameras.keys(): + myerror('camera {} not found in {}'.format(cam, cameras.keys())) + continue + cameras[cam]['shape'] = img.shape[:2] + if k2d[:, -1].sum() < 1: + continue + # calculate the area of the chessboard + mask = np.zeros_like(img[:, :, 0]) + k2d_int = np.round(k2d[:, :2]).astype(int) + if pattern is not None: + cv2.fillPoly(mask, [k2d_int[[0, pattern[0]-1, -1, -pattern[0]]]], 1) + else: + cv2.fillPoly(mask, [k2d_int[[0, 1, 2, 3, 0]]], 1) + area = mask.sum() + print(cam, area) + areas.append([cam, area, k2d, k3d]) + areas.sort(key=lambda x: -x[1]) + best_cam, area, k2d, k3d = areas[0] + # 先解决尺度问题 + ref_point_id = np.linalg.norm(k3d - k3d[:1], axis=-1).argmax() + k3d_pre = triangulate(cameras, areas) + length_gt = np.linalg.norm(k3d[0, :3] - k3d[ref_point_id, :3]) + length = np.linalg.norm(k3d_pre[0, :3] - k3d_pre[ref_point_id, :3]) + log('gt diag={:.3f}, est diag={:.3f}, scale={:.3f}'.format(length_gt, length, length_gt/length)) + scale_colmap = length_gt / length + for cam, camera in cameras.items(): + camera['T'] *= scale_colmap + k3d_pre = triangulate(cameras, areas) + length = np.linalg.norm(k3d_pre[0, :3] - k3d_pre[-1, :3]) + log('gt diag={:.3f}, est diag={:.3f}, scale={:.3f}'.format(length_gt, length, length_gt/length)) + # 计算相机相对于棋盘格的RT + if False: + for cam, _, k2d, k3d in areas: + K, dist = cameras[cam]['K'], cameras[cam]['dist'] + R, T = cameras[cam]['R'], cameras[cam]['T'] + err, rvec, tvec, kpts_repro = solvePnP(k3d, k2d, K, dist, flag=cv2.SOLVEPNP_ITERATIVE) + # 不同视角的计算的相对变换应该是一致的 + R_tgt = cv2.Rodrigues(rvec)[0] + T_tgt = tvec.reshape(3, 1) + R_rel, T_rel = compute_rel(R, T, R_tgt, T_tgt) + break + else: + # 使用估计的棋盘格坐标与实际的棋盘格坐标 + X = k3d_pre[:, :3] + X_gt = k3d[:, :3] + R_rel, T_rel = best_fit_transform(X_gt, X) + # 从棋盘格坐标系映射到colmap坐标系 + T_rel = T_rel.reshape(3, 1) + centers = [] + for cam, camera in cameras.items(): + camera.pop('Rvec') + R_old, T_old = camera['R'], camera['T'] + R_new = R_old @ R_rel + T_new = T_old + R_old @ T_rel + camera['R'] = R_new + camera['T'] = T_new + center = - camera['R'].T @ camera['T'] + centers.append(center) + print('{}: ({:6.3f}, {:.3f}, {:.3f})'.format(cam, *np.round(center.T[0], 3))) + # 使用棋盘格估计一下尺度 + k3d_pre = triangulate(cameras, areas) + length = np.linalg.norm(k3d_pre[0, :3] - k3d_pre[ref_point_id, :3]) + log('{} {} {}'.format(length_gt, length, length_gt/length)) + log(k3d_pre) + transform = np.eye(4) + transform[:3, :3] = R_rel + transform[:3, 3:] = T_rel + return cameras, scale_colmap, np.linalg.inv(transform) + +# for 3D points X, +# in origin world: \Pi(RX + T) = x +# in new world: \Pi(R'Y+T') = x +# , where X = R_pY + T_p + +if __name__ == '__main__': + import argparse + parser = argparse.ArgumentParser() + parser.add_argument('path', type=str) + parser.add_argument('out', type=str) + parser.add_argument('--plane_by_chessboard', type=str, default=None) + parser.add_argument('--plane_by_point', type=str, default=None) + parser.add_argument('--num', type=int, default=1) + parser.add_argument('--scale2d', type=float, default=None) + parser.add_argument('--prefix', type=str, default='') + parser.add_argument('--origin', type=str, default=None) + parser.add_argument('--guess_plane', action='store_true') + parser.add_argument('--noshow', action='store_true') + parser.add_argument('--debug', action='store_true') + args = parser.parse_args() + + if not os.path.exists(join(args.path, 'intri.yml')): + assert os.path.exists(join(args.path, 'cameras.bin')), os.listdir(args.path) + cmd = f'python3 apps/calibration/read_colmap.py {args.path} .bin' + run_cmd(cmd) + + cameras = read_cameras(args.path) + if args.plane_by_point is not None: + # 读入点云 + import ipdb; ipdb.set_trace() + if args.plane_by_chessboard is not None: + cameras, scale, transform = align_by_chessboard(cameras, args.plane_by_chessboard) + if not args.noshow: + import open3d as o3d + sparse_name = join(args.path, 'sparse.ply') + dense_name = join(args.path, '..', '..', 'dense', 'fused.ply') + if os.path.exists(dense_name): + pcd = o3d.io.read_point_cloud(dense_name) + else: + pcd = o3d.io.read_point_cloud(sparse_name) + save_json(join(args.out, 'transform.json'), {'scale': scale, 'transform': transform.tolist()}) + points = np.asarray(pcd.points) + # TODO: read correspondence of points3D and points2D + points_new = (scale*points) @ transform[:3, :3].T + transform[:3, 3:].T + pcd.points = o3d.utility.Vector3dVector(points_new) + o3d.io.write_point_cloud(join(args.out, 'sparse_aligned.ply'), pcd) + grids = [] + center = o3d.geometry.TriangleMesh.create_coordinate_frame( + size=1, origin=[0, 0, 0]) + grids.append(center) + for cam, camera in cameras.items(): + center = - camera['R'].T @ camera['T'] + center = o3d.geometry.TriangleMesh.create_coordinate_frame( + size=0.5, origin=[center[0, 0], center[1, 0], center[2, 0]]) + if cam.startswith(args.prefix): + center.paint_uniform_color([1, 0, 1]) + center.rotate(camera['R'].T) + grids.append(center) + o3d.visualization.draw_geometries([pcd] + grids) + write_camera(cameras, args.out) + if args.prefix is not None: + cameras_ = {} + for key, camera in cameras.items(): + if args.prefix not in key: + continue + cameras_[key.replace(args.prefix, '')] = camera + os.makedirs(join(args.out, args.prefix), exist_ok=True) + write_camera(cameras_, join(args.out, args.prefix)) \ No newline at end of file diff --git a/apps/calibration/read_colmap.py b/apps/calibration/read_colmap.py new file mode 100644 index 0000000..11bed28 --- /dev/null +++ b/apps/calibration/read_colmap.py @@ -0,0 +1,532 @@ +# Copyright (c) 2018, ETH Zurich and UNC Chapel Hill. +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# * Redistributions of source code must retain the above copyright +# notice, this list of conditions and the following disclaimer. +# +# * Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# +# * Neither the name of ETH Zurich and UNC Chapel Hill nor the names of +# its contributors may be used to endorse or promote products derived +# from this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +# Author: Johannes L. Schoenberger (jsch-at-demuc-dot-de) + +import os +import sys +import collections +import numpy as np +import struct +import cv2 + +CameraModel = collections.namedtuple( + "CameraModel", ["model_id", "model_name", "num_params"]) +Camera = collections.namedtuple( + "Camera", ["id", "model", "width", "height", "params"]) +BaseImage = collections.namedtuple( + "Image", ["id", "qvec", "tvec", "camera_id", "name", "xys", "point3D_ids"]) +Point3D = collections.namedtuple( + "Point3D", ["id", "xyz", "rgb", "error", "image_ids", "point2D_idxs"]) + +class Image(BaseImage): + def qvec2rotmat(self): + return qvec2rotmat(self.qvec) + + +CAMERA_MODELS = { + CameraModel(model_id=0, model_name="SIMPLE_PINHOLE", num_params=3), + CameraModel(model_id=1, model_name="PINHOLE", num_params=4), + CameraModel(model_id=2, model_name="SIMPLE_RADIAL", num_params=4), + CameraModel(model_id=3, model_name="RADIAL", num_params=5), + CameraModel(model_id=4, model_name="OPENCV", num_params=8), + CameraModel(model_id=5, model_name="OPENCV_FISHEYE", num_params=8), + CameraModel(model_id=6, model_name="FULL_OPENCV", num_params=12), + CameraModel(model_id=7, model_name="FOV", num_params=5), + CameraModel(model_id=8, model_name="SIMPLE_RADIAL_FISHEYE", num_params=4), + CameraModel(model_id=9, model_name="RADIAL_FISHEYE", num_params=5), + CameraModel(model_id=10, model_name="THIN_PRISM_FISHEYE", num_params=12) +} +CAMERA_MODEL_IDS = dict([(camera_model.model_id, camera_model) \ + for camera_model in CAMERA_MODELS]) +CAMERA_MODEL_NAMES = dict([(camera_model.model_name, camera_model) + for camera_model in CAMERA_MODELS]) + +def read_next_bytes(fid, num_bytes, format_char_sequence, endian_character="<"): + """Read and unpack the next bytes from a binary file. + :param fid: + :param num_bytes: Sum of combination of {2, 4, 8}, e.g. 2, 6, 16, 30, etc. + :param format_char_sequence: List of {c, e, f, d, h, H, i, I, l, L, q, Q}. + :param endian_character: Any of {@, =, <, >, !} + :return: Tuple of read and unpacked values. + """ + data = fid.read(num_bytes) + return struct.unpack(endian_character + format_char_sequence, data) + + +def read_cameras_text(path): + """ + see: src/base/reconstruction.cc + void Reconstruction::WriteCamerasText(const std::string& path) + void Reconstruction::ReadCamerasText(const std::string& path) + """ + cameras = {} + with open(path, "r") as fid: + while True: + line = fid.readline() + if not line: + break + line = line.strip() + if len(line) > 0 and line[0] != "#": + elems = line.split() + camera_id = int(elems[0]) + model = elems[1] + width = int(elems[2]) + height = int(elems[3]) + params = np.array(tuple(map(float, elems[4:]))) + cameras[camera_id] = Camera(id=camera_id, model=model, + width=width, height=height, + params=params) + return cameras + + +def read_cameras_binary(path_to_model_file): + """ + see: src/base/reconstruction.cc + void Reconstruction::WriteCamerasBinary(const std::string& path) + void Reconstruction::ReadCamerasBinary(const std::string& path) + """ + cameras = {} + with open(path_to_model_file, "rb") as fid: + num_cameras = read_next_bytes(fid, 8, "Q")[0] + for camera_line_index in range(num_cameras): + camera_properties = read_next_bytes( + fid, num_bytes=24, format_char_sequence="iiQQ") + camera_id = camera_properties[0] + model_id = camera_properties[1] + model_name = CAMERA_MODEL_IDS[camera_properties[1]].model_name + width = camera_properties[2] + height = camera_properties[3] + num_params = CAMERA_MODEL_IDS[model_id].num_params + params = read_next_bytes(fid, num_bytes=8*num_params, + format_char_sequence="d"*num_params) + cameras[camera_id] = Camera(id=camera_id, + model=model_name, + width=width, + height=height, + params=np.array(params)) + assert len(cameras) == num_cameras + return cameras + + +def read_images_text(path): + """ + see: src/base/reconstruction.cc + void Reconstruction::ReadImagesText(const std::string& path) + void Reconstruction::WriteImagesText(const std::string& path) + """ + images = {} + with open(path, "r") as fid: + while True: + line = fid.readline() + if not line: + break + line = line.strip() + if len(line) > 0 and line[0] != "#": + elems = line.split() + image_id = int(elems[0]) + qvec = np.array(tuple(map(float, elems[1:5]))) + tvec = np.array(tuple(map(float, elems[5:8]))) + camera_id = int(elems[8]) + image_name = elems[9] + elems = fid.readline().split() + xys = np.column_stack([tuple(map(float, elems[0::3])), + tuple(map(float, elems[1::3]))]) + point3D_ids = np.array(tuple(map(int, elems[2::3]))) + images[image_id] = Image( + id=image_id, qvec=qvec, tvec=tvec, + camera_id=camera_id, name=image_name, + xys=xys, point3D_ids=point3D_ids) + return images + + +def read_images_binary(path_to_model_file): + """ + see: src/base/reconstruction.cc + void Reconstruction::ReadImagesBinary(const std::string& path) + void Reconstruction::WriteImagesBinary(const std::string& path) + """ + images = {} + with open(path_to_model_file, "rb") as fid: + num_reg_images = read_next_bytes(fid, 8, "Q")[0] + for image_index in range(num_reg_images): + binary_image_properties = read_next_bytes( + fid, num_bytes=64, format_char_sequence="idddddddi") + image_id = binary_image_properties[0] + qvec = np.array(binary_image_properties[1:5]) + tvec = np.array(binary_image_properties[5:8]) + camera_id = binary_image_properties[8] + image_name = "" + current_char = read_next_bytes(fid, 1, "c")[0] + while current_char != b"\x00": # look for the ASCII 0 entry + image_name += current_char.decode("utf-8") + current_char = read_next_bytes(fid, 1, "c")[0] + num_points2D = read_next_bytes(fid, num_bytes=8, + format_char_sequence="Q")[0] + x_y_id_s = read_next_bytes(fid, num_bytes=24*num_points2D, + format_char_sequence="ddq"*num_points2D) + xys = np.column_stack([tuple(map(float, x_y_id_s[0::3])), + tuple(map(float, x_y_id_s[1::3]))]) + point3D_ids = np.array(tuple(map(int, x_y_id_s[2::3]))) + images[image_id] = Image( + id=image_id, qvec=qvec, tvec=tvec, + camera_id=camera_id, name=image_name, + xys=xys, point3D_ids=point3D_ids) + return images + + +def read_points3D_text(path): + """ + see: src/base/reconstruction.cc + void Reconstruction::ReadPoints3DText(const std::string& path) + void Reconstruction::WritePoints3DText(const std::string& path) + """ + points3D = {} + with open(path, "r") as fid: + while True: + line = fid.readline() + if not line: + break + line = line.strip() + if len(line) > 0 and line[0] != "#": + elems = line.split() + point3D_id = int(elems[0]) + xyz = np.array(tuple(map(float, elems[1:4]))) + rgb = np.array(tuple(map(int, elems[4:7]))) + error = float(elems[7]) + image_ids = np.array(tuple(map(int, elems[8::2]))) + point2D_idxs = np.array(tuple(map(int, elems[9::2]))) + points3D[point3D_id] = Point3D(id=point3D_id, xyz=xyz, rgb=rgb, + error=error, image_ids=image_ids, + point2D_idxs=point2D_idxs) + return points3D + + +def read_points3d_binary(path_to_model_file): + """ + see: src/base/reconstruction.cc + void Reconstruction::ReadPoints3DBinary(const std::string& path) + void Reconstruction::WritePoints3DBinary(const std::string& path) + """ + points3D = {} + with open(path_to_model_file, "rb") as fid: + num_points = read_next_bytes(fid, 8, "Q")[0] + for point_line_index in range(num_points): + binary_point_line_properties = read_next_bytes( + fid, num_bytes=43, format_char_sequence="QdddBBBd") + point3D_id = binary_point_line_properties[0] + xyz = np.array(binary_point_line_properties[1:4]) + rgb = np.array(binary_point_line_properties[4:7]) + error = np.array(binary_point_line_properties[7]) + track_length = read_next_bytes( + fid, num_bytes=8, format_char_sequence="Q")[0] + track_elems = read_next_bytes( + fid, num_bytes=8*track_length, + format_char_sequence="ii"*track_length) + image_ids = np.array(tuple(map(int, track_elems[0::2]))) + point2D_idxs = np.array(tuple(map(int, track_elems[1::2]))) + points3D[point3D_id] = Point3D( + id=point3D_id, xyz=xyz, rgb=rgb, + error=error, image_ids=image_ids, + point2D_idxs=point2D_idxs) + return points3D + + +def read_model(path, ext): + if ext == ".txt": + cameras = read_cameras_text(os.path.join(path, "cameras" + ext)) + images = read_images_text(os.path.join(path, "images" + ext)) + points3D = read_points3D_text(os.path.join(path, "points3D") + ext) + else: + cameras = read_cameras_binary(os.path.join(path, "cameras" + ext)) + images = read_images_binary(os.path.join(path, "images" + ext)) + points3D = read_points3d_binary(os.path.join(path, "points3D") + ext) + return cameras, images, points3D + + +def qvec2rotmat(qvec): + return np.array([ + [1 - 2 * qvec[2]**2 - 2 * qvec[3]**2, + 2 * qvec[1] * qvec[2] - 2 * qvec[0] * qvec[3], + 2 * qvec[3] * qvec[1] + 2 * qvec[0] * qvec[2]], + [2 * qvec[1] * qvec[2] + 2 * qvec[0] * qvec[3], + 1 - 2 * qvec[1]**2 - 2 * qvec[3]**2, + 2 * qvec[2] * qvec[3] - 2 * qvec[0] * qvec[1]], + [2 * qvec[3] * qvec[1] - 2 * qvec[0] * qvec[2], + 2 * qvec[2] * qvec[3] + 2 * qvec[0] * qvec[1], + 1 - 2 * qvec[1]**2 - 2 * qvec[2]**2]]) + + +def rotmat2qvec(R): + Rxx, Ryx, Rzx, Rxy, Ryy, Rzy, Rxz, Ryz, Rzz = R.flat + K = np.array([ + [Rxx - Ryy - Rzz, 0, 0, 0], + [Ryx + Rxy, Ryy - Rxx - Rzz, 0, 0], + [Rzx + Rxz, Rzy + Ryz, Rzz - Rxx - Ryy, 0], + [Ryz - Rzy, Rzx - Rxz, Rxy - Ryx, Rxx + Ryy + Rzz]]) / 3.0 + eigvals, eigvecs = np.linalg.eigh(K) + qvec = eigvecs[[3, 0, 1, 2], np.argmax(eigvals)] + if qvec[0] < 0: + qvec *= -1 + return qvec + + +def write_cameras_text(cameras, path): + """ + see: src/base/reconstruction.cc + void Reconstruction::WriteCamerasText(const std::string& path) + void Reconstruction::ReadCamerasText(const std::string& path) + """ + HEADER = '# Camera list with one line of data per camera:\n' + '# CAMERA_ID, MODEL, WIDTH, HEIGHT, PARAMS[]\n' + '# Number of cameras: {}\n'.format(len(cameras)) + with open(path, "w") as fid: + fid.write(HEADER) + for _, cam in cameras.items(): + to_write = [cam.id, cam.model, cam.width, cam.height, *cam.params] + line = " ".join([str(elem) for elem in to_write]) + fid.write(line + "\n") + +def write_next_bytes(fid, data, format_char_sequence, endian_character="<"): + """pack and write to a binary file. + :param fid: + :param data: data to send, if multiple elements are sent at the same time, + they should be encapsuled either in a list or a tuple + :param format_char_sequence: List of {c, e, f, d, h, H, i, I, l, L, q, Q}. + should be the same length as the data list or tuple + :param endian_character: Any of {@, =, <, >, !} + """ + if isinstance(data, (list, tuple)): + bytes = struct.pack(endian_character + format_char_sequence, *data) + else: + bytes = struct.pack(endian_character + format_char_sequence, data) + fid.write(bytes) + +def write_cameras_binary(cameras, path_to_model_file): + """ + see: src/base/reconstruction.cc + void Reconstruction::WriteCamerasBinary(const std::string& path) + void Reconstruction::ReadCamerasBinary(const std::string& path) + """ + with open(path_to_model_file, "wb") as fid: + write_next_bytes(fid, len(cameras), "Q") + for _, cam in cameras.items(): + model_id = CAMERA_MODEL_NAMES[cam.model].model_id + camera_properties = [cam.id, + model_id, + cam.width, + cam.height] + write_next_bytes(fid, camera_properties, "iiQQ") + for p in cam.params: + write_next_bytes(fid, float(p), "d") + return cameras + + +def write_images_binary(images, path_to_model_file): + """ + see: src/base/reconstruction.cc + void Reconstruction::ReadImagesBinary(const std::string& path) + void Reconstruction::WriteImagesBinary(const std::string& path) + """ + with open(path_to_model_file, "wb") as fid: + write_next_bytes(fid, len(images), "Q") + for _, img in images.items(): + write_next_bytes(fid, img.id, "i") + write_next_bytes(fid, img.qvec.tolist(), "dddd") + write_next_bytes(fid, img.tvec.tolist(), "ddd") + write_next_bytes(fid, img.camera_id, "i") + for char in img.name: + write_next_bytes(fid, char.encode("utf-8"), "c") + write_next_bytes(fid, b"\x00", "c") + write_next_bytes(fid, len(img.point3D_ids), "Q") + for xy, p3d_id in zip(img.xys, img.point3D_ids): + write_next_bytes(fid, [*xy, p3d_id], "ddq") + +def write_images_text(images, path): + """ + see: src/base/reconstruction.cc + void Reconstruction::ReadImagesText(const std::string& path) + void Reconstruction::WriteImagesText(const std::string& path) + """ + if len(images) == 0: + mean_observations = 0 + else: + mean_observations = sum((len(img.point3D_ids) for _, img in images.items()))/len(images) + HEADER = '# Image list with two lines of data per image:\n' + '# IMAGE_ID, QW, QX, QY, QZ, TX, TY, TZ, CAMERA_ID, NAME\n' + '# POINTS2D[] as (X, Y, POINT3D_ID)\n' + '# Number of images: {}, mean observations per image: {}\n'.format(len(images), mean_observations) + + with open(path, "w") as fid: + fid.write(HEADER) + for _, img in images.items(): + image_header = [img.id, *img.qvec, *img.tvec, img.camera_id, img.name] + first_line = " ".join(map(str, image_header)) + fid.write(first_line + "\n") + + points_strings = [] + for xy, point3D_id in zip(img.xys, img.point3D_ids): + points_strings.append(" ".join(map(str, [*xy, point3D_id]))) + fid.write(" ".join(points_strings) + "\n") + +def write_points3D_text(points3D, path): + """ + see: src/base/reconstruction.cc + void Reconstruction::ReadPoints3DText(const std::string& path) + void Reconstruction::WritePoints3DText(const std::string& path) + """ + if len(points3D) == 0: + mean_track_length = 0 + else: + mean_track_length = sum((len(pt.image_ids) for _, pt in points3D.items()))/len(points3D) + HEADER = '# 3D point list with one line of data per point:\n' + '# POINT3D_ID, X, Y, Z, R, G, B, ERROR, TRACK[] as (IMAGE_ID, POINT2D_IDX)\n' + '# Number of points: {}, mean track length: {}\n'.format(len(points3D), mean_track_length) + + with open(path, "w") as fid: + fid.write(HEADER) + for _, pt in points3D.items(): + point_header = [pt.id, *pt.xyz, *pt.rgb, pt.error] + fid.write(" ".join(map(str, point_header)) + " ") + track_strings = [] + for image_id, point2D in zip(pt.image_ids, pt.point2D_idxs): + track_strings.append(" ".join(map(str, [image_id, point2D]))) + fid.write(" ".join(track_strings) + "\n") + + +def write_points3d_binary(points3D, path_to_model_file): + """ + see: src/base/reconstruction.cc + void Reconstruction::ReadPoints3DBinary(const std::string& path) + void Reconstruction::WritePoints3DBinary(const std::string& path) + """ + with open(path_to_model_file, "wb") as fid: + write_next_bytes(fid, len(points3D), "Q") + for _, pt in points3D.items(): + write_next_bytes(fid, pt.id, "Q") + write_next_bytes(fid, pt.xyz.tolist(), "ddd") + write_next_bytes(fid, pt.rgb.tolist(), "BBB") + write_next_bytes(fid, pt.error, "d") + track_length = pt.image_ids.shape[0] + write_next_bytes(fid, track_length, "Q") + for image_id, point2D_id in zip(pt.image_ids, pt.point2D_idxs): + write_next_bytes(fid, [image_id, point2D_id], "ii") + + +class FileStorage(object): + def __init__(self, filename, isWrite=False): + version = cv2.__version__ + self.version = int(version.split('.')[0]) + if isWrite: + self.fs = cv2.FileStorage(filename, cv2.FILE_STORAGE_WRITE) + else: + self.fs = cv2.FileStorage(filename, cv2.FILE_STORAGE_READ) + + def __del__(self): + cv2.FileStorage.release(self.fs) + + def write(self, key, value, dt='mat'): + if dt == 'mat': + cv2.FileStorage.write(self.fs, key, value) + elif dt == 'list': + if self.version == 4: # 4.4 + # self.fs.write(key, '[') + # for elem in value: + # self.fs.write('none', elem) + # self.fs.write('none', ']') + # import ipdb; ipdb.set_trace() + self.fs.startWriteStruct(key, cv2.FileNode_SEQ) + for elem in value: + self.fs.write('', elem) + self.fs.endWriteStruct() + else: # 3.4 + self.fs.write(key, '[') + for elem in value: + self.fs.write('none', elem) + self.fs.write('none', ']') + +def main(): + if len(sys.argv) != 3: + print("Usage: python read_model.py path/to/model/folder [.txt,.bin]") + return + + cameras, images, points3D = read_model(path=sys.argv[1], ext=sys.argv[2]) + import cv2 + cameras_out = {} + for key in cameras.keys(): + p = cameras[key].params + if cameras[key].model == 'SIMPLE_RADIAL': + f, cx, cy, k = p + K = np.array([f, 0, cx, 0, f, cy, 0, 0, 1]).reshape(3, 3) + dist = np.array([[k, 0, 0, 0, 0]]) + elif cameras[key].model == 'PINHOLE': + fx, fy, cx, cy = p + K = np.array([fx, 0, cx, 0, fy, cy, 0, 0, 1]).reshape(3, 3) + dist = np.array([[0, 0, 0, 0, 0]]) + else: + K = np.array([[p[0], 0, p[2], 0, p[1], p[3], 0, 0, 1]]).reshape(3, 3) + dist = np.array([[p[4], p[5], p[6], p[7], 0.]]) + cameras_out[key] = {'K': K, 'dist': dist} + mapkey = {} + cameras_new = {} + for key, val in images.items(): + cam = cameras_out[val.camera_id].copy() + t = val.tvec.reshape(3, 1) + R = qvec2rotmat(val.qvec) + cam['R'] = R + cam['Rvec'] = cv2.Rodrigues(R)[0] + cam['T'] = t + # mapkey[val.name.split('.')[0]] = val.camera_id + + cameras_new[val.name.split('.')[0]] = cam + # cameras_new[val.name.split('.')[0].split('/')[0]] = cam + keys = sorted(list(cameras_new.keys())) + cameras_new = {key:cameras_new[key] for key in keys} + print("num_cameras: {}/{}".format(len(cameras), len(cameras_new))) + print("num_images:", len(images)) + print("num_points3D:", len(points3D)) + if len(points3D) > 0: + keys = list(points3D.keys()) + xyz = np.stack([points3D[k].xyz for k in keys]) + rgb = np.stack([points3D[k].rgb for k in keys]) + import open3d as o3d + pcd = o3d.geometry.PointCloud() + pcd.points = o3d.utility.Vector3dVector(xyz) + pcd.colors = o3d.utility.Vector3dVector(rgb/255.) + from os.path import join + pcdname = join(sys.argv[1], 'sparse.ply') + o3d.io.write_point_cloud(pcdname, pcd) + if False: + o3d.visualization.draw_geometries([pcd]) + from easymocap.mytools.camera_utils import write_camera + write_camera(cameras_new, sys.argv[1]) + + +if __name__ == "__main__": + main()