add scripts to generate poses

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.

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)))

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.)

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()

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)
+