From 08f98f8f3a491e3d6a70d57ded749ad5d2a4d157 Mon Sep 17 00:00:00 2001 From: Kai-46 Date: Tue, 20 Oct 2020 22:54:43 -0400 Subject: [PATCH] add scripts to generate poses --- README.md | 9 + colmap_runner/extract_sfm.py | 133 ++++++++ colmap_runner/normalize_cam_dict.py | 57 ++++ colmap_runner/read_write_model.py | 483 ++++++++++++++++++++++++++++ colmap_runner/run_colmap.py | 166 ++++++++++ 5 files changed, 848 insertions(+) create mode 100644 colmap_runner/extract_sfm.py create mode 100644 colmap_runner/normalize_cam_dict.py create mode 100755 colmap_runner/read_write_model.py create mode 100644 colmap_runner/run_colmap.py diff --git a/README.md b/README.md index e907213..ae71921 100644 --- a/README.md +++ b/README.md @@ -57,3 +57,12 @@ Plese cite our work if you use the code. year = {2020}, } ``` + +## Generating poses with [COLMAP SfM](https://colmap.github.io/) +You can use the scripts inside 'colmap_runner/' to generate camera poses from images with COLMAP SfM. +* Specify 'img_dir' and 'out_dir' in 'colmap_runner/run_colmap.py'. +* Inside 'colmap_runner/', execute command 'python run_colmap.py'. +* After program finishes running, you would see the posed images in the folder 'out_dir/posed_images'. + * Distortion-free images are inside 'out_dir/posed_images/images'. + * Raw COLMAP poses is stored as a json file 'out_dir/posed_images/kai_cameras.json'. + * Normalized camderas is stored in 'out_dir/posed_images/kai_cameras_normalized.json'. See the 'Scene normalization method' in the 'Data' section. diff --git a/colmap_runner/extract_sfm.py b/colmap_runner/extract_sfm.py new file mode 100644 index 0000000..d04ba7d --- /dev/null +++ b/colmap_runner/extract_sfm.py @@ -0,0 +1,133 @@ +from read_write_model import read_model +import numpy as np +import json +import os +from pyquaternion import Quaternion +import trimesh + + +def parse_tracks(colmap_images, colmap_points3D): + all_tracks = [] # list of dicts; each dict represents a track + all_points = [] # list of all 3D points + view_keypoints = {} # dict of lists; each list represents the triangulated key points of a view + + + for point3D_id in colmap_points3D: + point3D = colmap_points3D[point3D_id] + image_ids = point3D.image_ids + point2D_idxs = point3D.point2D_idxs + + cur_track = {} + cur_track['xyz'] = (point3D.xyz[0], point3D.xyz[1], point3D.xyz[2]) + cur_track['err'] = point3D.error.item() + + cur_track_len = len(image_ids) + assert (cur_track_len == len(point2D_idxs)) + all_points.append(list(cur_track['xyz'] + (cur_track['err'], cur_track_len) + tuple(point3D.rgb))) + + pixels = [] + for i in range(cur_track_len): + image = colmap_images[image_ids[i]] + img_name = image.name + point2D_idx = point2D_idxs[i] + point2D = image.xys[point2D_idx] + assert (image.point3D_ids[point2D_idx] == point3D_id) + pixels.append((img_name, point2D[0], point2D[1])) + + if img_name not in view_keypoints: + view_keypoints[img_name] = [(point2D[0], point2D[1]) + cur_track['xyz'] + (cur_track_len, ), ] + else: + view_keypoints[img_name].append((point2D[0], point2D[1]) + cur_track['xyz'] + (cur_track_len, )) + + cur_track['pixels'] = sorted(pixels, key=lambda x: x[0]) # sort pixels by the img_name + all_tracks.append(cur_track) + + return all_tracks, all_points, view_keypoints + + +def parse_camera_dict(colmap_cameras, colmap_images): + camera_dict = {} + for image_id in colmap_images: + image = colmap_images[image_id] + + img_name = image.name + cam = colmap_cameras[image.camera_id] + + # print(cam) + assert(cam.model == 'PINHOLE') + + img_size = [cam.width, cam.height] + params = list(cam.params) + qvec = list(image.qvec) + tvec = list(image.tvec) + + # w, h, fx, fy, cx, cy, qvec, tvec + # camera_dict[img_name] = img_size + params + qvec + tvec + camera_dict[img_name] = {} + camera_dict[img_name]['img_size'] = img_size + + fx, fy, cx, cy = params + K = np.eye(4) + K[0, 0] = fx + K[1, 1] = fy + K[0, 2] = cx + K[1, 2] = cy + camera_dict[img_name]['K'] = list(K.flatten()) + + rot = Quaternion(qvec[0], qvec[1], qvec[2], qvec[3]).rotation_matrix + W2C = np.eye(4) + W2C[:3, :3] = rot + W2C[:3, 3] = np.array(tvec) + camera_dict[img_name]['W2C'] = list(W2C.flatten()) + + return camera_dict + + +def extract_all_to_dir(sparse_dir, out_dir, ext='.bin'): + if not os.path.exists(out_dir): + os.mkdir(out_dir) + + camera_dict_file = os.path.join(out_dir, 'kai_cameras.json') + xyz_file = os.path.join(out_dir, 'kai_points.txt') + track_file = os.path.join(out_dir, 'kai_tracks.json') + keypoints_file = os.path.join(out_dir, 'kai_keypoints.json') + + colmap_cameras, colmap_images, colmap_points3D = read_model(sparse_dir, ext) + + camera_dict = parse_camera_dict(colmap_cameras, colmap_images) + with open(camera_dict_file, 'w') as fp: + json.dump(camera_dict, fp, indent=2, sort_keys=True) + + all_tracks, all_points, view_keypoints = parse_tracks(colmap_images, colmap_points3D) + all_points = np.array(all_points) + np.savetxt(xyz_file, all_points, header='# format: x, y, z, reproj_err, track_len, color(RGB)', fmt='%.6f') + + mesh = trimesh.Trimesh(vertices=all_points[:, :3].astype(np.float32), + vertex_colors=all_points[:, -3:].astype(np.uint8)) + mesh.export(os.path.join(out_dir, 'kai_points.ply')) + + with open(track_file, 'w') as fp: + json.dump(all_tracks, fp) + + with open(keypoints_file, 'w') as fp: + json.dump(view_keypoints, fp) + + +if __name__ == '__main__': + mvs_dir = '/home/zhangka2/sg_render/run_mvs/scan114_train_5/colmap_mvs/mvs' + sparse_dir = os.path.join(mvs_dir, 'sparse') + out_dir = os.path.join(mvs_dir, 'sparse_inspect') + extract_all_to_dir(sparse_dir, out_dir) + + xyz_file = os.path.join(out_dir, 'kai_points.txt') + reproj_errs = np.loadtxt(xyz_file)[:, 3] + with open(os.path.join(out_dir, 'stats.txt'), 'w') as fp: + fp.write('reprojection errors (px) in SfM:\n') + fp.write(' percentile value\n') + for a in [50, 70, 90, 99]: + fp.write(' {} {:.3f}\n'.format(a, np.percentile(reproj_errs, a))) + + print('reprojection errors (px) in SfM:') + print(' percentile value') + for a in [50, 70, 90, 99]: + print(' {} {:.3f}'.format(a, np.percentile(reproj_errs, a))) \ No newline at end of file diff --git a/colmap_runner/normalize_cam_dict.py b/colmap_runner/normalize_cam_dict.py new file mode 100644 index 0000000..f32a04c --- /dev/null +++ b/colmap_runner/normalize_cam_dict.py @@ -0,0 +1,57 @@ +import numpy as np +import json +import copy + + +def get_tf_cams(cam_dict, target_radius=1.): + cam_centers = [] + for im_name in cam_dict: + W2C = np.array(cam_dict[im_name]['W2C']).reshape((4, 4)) + C2W = np.linalg.inv(W2C) + cam_centers.append(C2W[:3, 3:4]) + + def get_center_and_diag(cam_centers): + cam_centers = np.hstack(cam_centers) + avg_cam_center = np.mean(cam_centers, axis=1, keepdims=True) + center = avg_cam_center + dist = np.linalg.norm(cam_centers - center, axis=0, keepdims=True) + diagonal = np.max(dist) + return center.flatten(), diagonal + + center, diagonal = get_center_and_diag(cam_centers) + radius = diagonal / 2. * 1.1 + + translate = -center + scale = target_radius / radius + + return translate, scale + + +def normalize_cam_dict(in_cam_dict_file, out_cam_dict_file, target_radius=1.): + with open(in_cam_dict_file) as fp: + in_cam_dict = json.load(fp) + + translate, scale = get_tf_cams(in_cam_dict, target_radius=target_radius) + + def transform_pose(W2C, translate, scale): + C2W = np.linalg.inv(W2C) + cam_center = C2W[:3, 3] + cam_center = (cam_center + translate) * scale + C2W[:3, 3] = cam_center + return np.linalg.inv(C2W) + + out_cam_dict = copy.deepcopy(in_cam_dict) + for img_name in out_cam_dict: + W2C = np.array(out_cam_dict[img_name]['W2C']).reshape((4, 4)) + W2C = transform_pose(W2C, translate, scale) + assert(np.isclose(np.linalg.det(W2C[:3, :3]), 1.)) + out_cam_dict[img_name]['W2C'] = list(W2C.flatten()) + + with open(out_cam_dict_file, 'w') as fp: + json.dump(out_cam_dict, fp, indent=2, sort_keys=True) + + +if __name__ == '__main__': + in_cam_dict_file = '' + out_cam_dict_file = '' + normalize_cam_dict(in_cam_dict_file, out_cam_dict_file, target_radius=1.) diff --git a/colmap_runner/read_write_model.py b/colmap_runner/read_write_model.py new file mode 100755 index 0000000..6d195c7 --- /dev/null +++ b/colmap_runner/read_write_model.py @@ -0,0 +1,483 @@ +# 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 argparse + + +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 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 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 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_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 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 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_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 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 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") + + +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 write_model(cameras, images, points3D, path, ext): + if ext == ".txt": + write_cameras_text(cameras, os.path.join(path, "cameras" + ext)) + write_images_text(images, os.path.join(path, "images" + ext)) + write_points3D_text(points3D, os.path.join(path, "points3D") + ext) + else: + write_cameras_binary(cameras, os.path.join(path, "cameras" + ext)) + write_images_binary(images, os.path.join(path, "images" + ext)) + write_points3d_binary(points3D, 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 main(): + parser = argparse.ArgumentParser(description='Read and write COLMAP binary and text models') + parser.add_argument('input_model', help='path to input model folder') + parser.add_argument('input_format', choices=['.bin', '.txt'], + help='input model format') + parser.add_argument('--output_model', metavar='PATH', + help='path to output model folder') + parser.add_argument('--output_format', choices=['.bin', '.txt'], + help='outut model format', default='.txt') + args = parser.parse_args() + + cameras, images, points3D = read_model(path=args.input_model, ext=args.input_format) + + print("num_cameras:", len(cameras)) + print("num_images:", len(images)) + print("num_points3D:", len(points3D)) + + if args.output_model is not None: + write_model(cameras, images, points3D, path=args.output_model, ext=args.output_format) + + +if __name__ == "__main__": + main() diff --git a/colmap_runner/run_colmap.py b/colmap_runner/run_colmap.py new file mode 100644 index 0000000..ad93f03 --- /dev/null +++ b/colmap_runner/run_colmap.py @@ -0,0 +1,166 @@ +import os +import subprocess +from extract_sfm import extract_all_to_dir +from normalize_cam_dict import normalize_cam_dict + +######################################################################### +# Note: configure the colmap_bin to the colmap executable on your machine +######################################################################### + +def bash_run(cmd): + colmap_bin = '/home/zhangka2/code/colmap/build/__install__/bin/colmap' + cmd = colmap_bin + ' ' + cmd + print('\nRunning cmd: ', cmd) + + subprocess.check_call(['/bin/bash', '-c', cmd]) + + +gpu_index = '-1' + + +def run_sift_matching(img_dir, db_file, remove_exist=False): + print('Running sift matching...') + + if remove_exist and os.path.exists(db_file): + os.remove(db_file) # otherwise colmap will skip sift matching + + # feature extraction + # if there's no attached display, cannot use feature extractor with GPU + cmd = ' feature_extractor --database_path {} \ + --image_path {} \ + --ImageReader.single_camera 1 \ + --ImageReader.camera_model SIMPLE_RADIAL \ + --SiftExtraction.max_image_size 5000 \ + --SiftExtraction.estimate_affine_shape 0 \ + --SiftExtraction.domain_size_pooling 1 \ + --SiftExtraction.use_gpu 1 \ + --SiftExtraction.max_num_features 16384 \ + --SiftExtraction.gpu_index {}'.format(db_file, img_dir, gpu_index) + bash_run(cmd) + + # feature matching + cmd = ' exhaustive_matcher --database_path {} \ + --SiftMatching.guided_matching 1 \ + --SiftMatching.use_gpu 1 \ + --SiftMatching.max_num_matches 65536 \ + --SiftMatching.max_error 3 \ + --SiftMatching.gpu_index {}'.format(db_file, gpu_index) + + bash_run(cmd) + + +def run_sfm(img_dir, db_file, out_dir): + print('Running SfM...') + + cmd = ' mapper \ + --database_path {} \ + --image_path {} \ + --output_path {} \ + --Mapper.tri_min_angle 3.0 \ + --Mapper.filter_min_tri_angle 3.0'.format(db_file, img_dir, out_dir) + + bash_run(cmd) + + +def prepare_mvs(img_dir, sparse_dir, mvs_dir): + print('Preparing for MVS...') + + cmd = ' image_undistorter \ + --image_path {} \ + --input_path {} \ + --output_path {} \ + --output_type COLMAP \ + --max_image_size 2000'.format(img_dir, sparse_dir, mvs_dir) + + bash_run(cmd) + + +def run_photometric_mvs(mvs_dir, window_radius): + print('Running photometric MVS...') + + cmd = ' patch_match_stereo --workspace_path {} \ + --PatchMatchStereo.window_radius {} \ + --PatchMatchStereo.min_triangulation_angle 3.0 \ + --PatchMatchStereo.filter 1 \ + --PatchMatchStereo.geom_consistency 1 \ + --PatchMatchStereo.gpu_index={} \ + --PatchMatchStereo.num_samples 15 \ + --PatchMatchStereo.num_iterations 12'.format(mvs_dir, + window_radius, gpu_index) + + bash_run(cmd) + + +def run_fuse(mvs_dir, out_ply): + print('Running depth fusion...') + + cmd = ' stereo_fusion --workspace_path {} \ + --output_path {} \ + --input_type geometric'.format(mvs_dir, out_ply) + + bash_run(cmd) + + +def run_possion_mesher(in_ply, out_ply, trim): + print('Running possion mesher...') + + cmd = ' poisson_mesher \ + --input_path {} \ + --output_path {} \ + --PoissonMeshing.trim {}'.format(in_ply, out_ply, trim) + + bash_run(cmd) + + +def main(img_dir, out_dir, run_mvs=False): + os.makedirs(out_dir, exist_ok=True) + + #### run sfm + sfm_dir = os.path.join(out_dir, 'sfm') + os.makedirs(sfm_dir, exist_ok=True) + + img_dir_link = os.path.join(sfm_dir, 'images') + if os.path.exists(img_dir_link): + os.remove(img_dir_link) + os.symlink(img_dir, img_dir_link) + + db_file = os.path.join(sfm_dir, 'database.db') + run_sift_matching(img_dir, db_file, remove_exist=False) + sparse_dir = os.path.join(sfm_dir, 'sparse') + os.makedirs(sparse_dir, exist_ok=True) + run_sfm(img_dir, db_file, sparse_dir) + + # undistort images + mvs_dir = os.path.join(out_dir, 'mvs') + os.makedirs(mvs_dir, exist_ok=True) + prepare_mvs(img_dir, sparse_dir, mvs_dir) + + # extract camera parameters and undistorted images + os.makedirs(os.path.join(out_dir, 'posed_images'), exist_ok=True) + extract_all_to_dir(os.path.join(mvs_dir, 'sparse'), os.path.join(out_dir, 'posed_images')) + undistorted_img_dir = os.path.join(mvs_dir, 'images') + posed_img_dir_link = os.path.join(out_dir, 'posed_images/images') + if os.path.exists(posed_img_dir_link): + os.remove(posed_img_dir_link) + os.symlink(undistorted_img_dir, posed_img_dir_link) + # normalize average camera center to origin, and put all cameras inside the unit sphere + normalize_cam_dict(os.path.join(out_dir, 'posed_images/kai_cameras.json'), + os.path.join(out_dir, 'posed_images/kai_cameras_normalized.json')) + + if run_mvs: + # run mvs + run_photometric_mvs(mvs_dir, window_radius=7) + + out_ply = os.path.join(out_dir, 'mvs/fused.ply') + run_fuse(mvs_dir, out_ply) + + out_mesh_ply = os.path.join(out_dir, 'mvs/meshed_trim_3.ply') + run_possion_mesher(out_ply, out_mesh_ply, trim=3) + + +if __name__ == '__main__': + img_dir = '' + out_dir = '' + run_mvs = False + main(img_dir, out_dir, run_mvs=run_mvs) +