COLMAP已知相机内外参数重建稀疏/稠密模型

Posted 2020-11-23 li-minghao

COLMAP已知相机内外参数重建稀疏/稠密模型

Reconstruct sparse/dense model from known camera poses

参考官方Faq链接：https://colmap.github.io/faq.html#reconstruct-sparse-dense-model-from-known-camera-poses

1. 手动指定相机Pose和注册图像

在目录下手动新建cameras.txt, images.txt, 和 points3D.txt三个文本文件，目录示例：

+── %ProjectPath%/created/sparse
│   +── cameras.txt
│   +── images.txt
│   +── points3D.txt

将内参(camera intrinsics) 放入cameras.txt， 外参(camera extrinsics)放入 images.txt , points3D.txt 为空。具体格式如下：

cameras.txt 示例：

# Camera list with one line of data per camera:
# CAMERA_ID, MODEL, WIDTH, HEIGHT, PARAMS[fx,fy,cx,cy]
# Number of cameras: 2
1 PINHOLE 1280 720 771.904 771.896 639.993 360.001
2 PINHOLE 1280 720 771.899 771.898 639.999 360.001

COLMAP 设定了 SIMPLE_PINHOLE ，PINHOLE ， SIMPLE_RADIAL ，RADIAL， OPENCV ， FULL_OPENCV， SIMPLE_RADIAL_FISHEYE ，RADIAL_FISHEYE，OPENCV_FISHEYE，FOV，THIN_PRISM_FISHEYE 共11种相机模型，其中最常用的为PINHOLE，即小孔相机模型。一般正常拍摄的照片，不考虑畸变即为该模型。手动参数的情况下官方推荐OPENCV相机模型，相比PINHOLE，考虑了xy轴畸变。各模型区别具体可参考Camera Models 查看。请根据自己的已有相机参数选用合适的相机模型。

images.txt 示例：

# Image list with two lines of data per image:
# IMAGE_ID, QW, QX, QY, QZ, TX, TY, TZ, CAMERA_ID, NAME
# POINTS2D[] as (X, Y, POINT3D_ID)
# Number of images: 2, mean observations per image: 2
1 0.430115 0.411564 0.555504 -0.580543 10468.491287 380.313066 1720.465175 1 image001.jpg
# Make sure every other line is left empty
2 0.309712 0.337960 0.655221 -0.600456 10477.663284 446.4208 -1633.886712 2 image002.jpg

3 0.375916 0.401654 0.609703 -0.570633 10592.122754 263.672534 600.636247 3 image003.jpg

COLMAP 中多张图片可以使用一个相机模型，见 IMAGE_ID 所对应的CAMERA_ID， images.txt 中，每张图片由一行图像外参，一行2D点组成。在手动构建中由于没有2D点，第二行留空。QW, QX, QY, QZ 为四元数表示的相机旋转信息，TX, TY, TZ为平移向量。

接下来的命令会依赖目录机构，因此需先明确各级目录，本文中的命令参考以下目录结构：

E:CODECOLMAPMANUALProjectPath
├─created
│  └─sparse
│   +── cameras.txt
│   +── images.txt
│   +── points3D.txt
├─dense
├─images
│   +── images001.jpg
│   +── images002.jpg
│   +── ...
└─triangulated
    └─sparse

images/ 文件夹下放置需进行重建的图片，created/sparse/ 文件夹下放入上一步创建的 cameras.txt， images.txt 和 points3D.txt 。

2. 命令行执行稀疏重建

注意： 也可以不执行稀疏重建，直接执行MVS步骤生成深度图，可直接跳至 3.2?手动稠密重建 。

2.1. 抽取图像特征

在 ProjectPath 目录下运行命令：

colmap feature_extractor --database_path database.db --image_path images

目录下会自动生成一个 database.db 文件。终端输出示例：

==============================================================================
Feature extraction
==============================================================================
Processed file [1/25]
  Name:            image001.jpg
  Dimensions:      1280 x 720
  Camera:          #1 - SIMPLE_RADIAL
  Focal Length:    1536.00px
  Features:        10300
Processed file [2/25]
...
...
Elapsed time: 0.039 [minutes]

可以看出此时使用的相机模型是SIMPLE_RADIAL，这是没有关系的，因为此时只是抽取特征，官方文档的说法是此步之后才需要将 cameras.txt 中的相机内参复制到 database.db 中。

2.1.手动导入相机内参

有两种导入方式：

• 运行COLMAP 的GUI程序，选择 File->New Project->Database->Open 打开我们刚刚建立的 database.db ，Images->Select 选择存放图片的目录。点击 Save 保存。选择 Processing->Database management ，点击 Cameras 查看相机模型。这里就需要一行一行选中已有的相机模型进行修改。选中后点击 Set model ，先选中相机的model类型，如 PINHOLE ，或 OPENCV，注意不同相机的 params 个数可能不一样：

• 使用Python脚本更改 database.db ，官方在colmap源码 scriptspythondatabase.py 提供了一个fake的创建database脚本。但我们的目的其实是 update 原来的相机模型，这里我写了一个python脚本，其读取cameras.txt 中的数据，更新当前目录下的 database.db 中的相机模型 。代码如下：

  # This script is based on an original implementation by True Price.
  # Created by liminghao
  import sys
  import numpy as np
  import sqlite3
  
  IS_PYTHON3 = sys.version_info[0] >= 3
  
  def array_to_blob(array):
      if IS_PYTHON3:
          return array.tostring()
      else:
          return np.getbuffer(array)
  
  def blob_to_array(blob, dtype, shape=(-1,)):
      if IS_PYTHON3:
          return np.fromstring(blob, dtype=dtype).reshape(*shape)
      else:
          return np.frombuffer(blob, dtype=dtype).reshape(*shape)
  
  class COLMAPDatabase(sqlite3.Connection):
  
      @staticmethod
      def connect(database_path):
          return sqlite3.connect(database_path, factory=COLMAPDatabase)
  
      def __init__(self, *args, **kwargs):
          super(COLMAPDatabase, self).__init__(*args, **kwargs)
  
          self.create_tables = lambda: self.executescript(CREATE_ALL)
          self.create_cameras_table =             lambda: self.executescript(CREATE_CAMERAS_TABLE)
          self.create_descriptors_table =             lambda: self.executescript(CREATE_DESCRIPTORS_TABLE)
          self.create_images_table =             lambda: self.executescript(CREATE_IMAGES_TABLE)
          self.create_two_view_geometries_table =             lambda: self.executescript(CREATE_TWO_VIEW_GEOMETRIES_TABLE)
          self.create_keypoints_table =             lambda: self.executescript(CREATE_KEYPOINTS_TABLE)
          self.create_matches_table =             lambda: self.executescript(CREATE_MATCHES_TABLE)
          self.create_name_index = lambda: self.executescript(CREATE_NAME_INDEX)
  
      def update_camera(self, model, width, height, params, camera_id):
          params = np.asarray(params, np.float64)
          cursor = self.execute(
              "UPDATE cameras SET model=?, width=?, height=?, params=?, prior_focal_length=True WHERE camera_id=?",
              (model, width, height, array_to_blob(params),camera_id))
          return cursor.lastrowid
  
  def camTodatabase(txtfile):
      import os
      import argparse
  
      camModelDict = {'SIMPLE_PINHOLE': 0,
                      'PINHOLE': 1,
                      'SIMPLE_RADIAL': 2,
                      'RADIAL': 3,
                      'OPENCV': 4,
                      'FULL_OPENCV': 5,
                      'SIMPLE_RADIAL_FISHEYE': 6,
                      'RADIAL_FISHEYE': 7,
                      'OPENCV_FISHEYE': 8,
                      'FOV': 9,
                      'THIN_PRISM_FISHEYE': 10}
      parser = argparse.ArgumentParser()
      parser.add_argument("--database_path", default="database.db")
      args = parser.parse_args()
      if os.path.exists(args.database_path)==False:
          print("ERROR: database path dosen't exist -- please check database.db.")
          return
      # Open the database.
      db = COLMAPDatabase.connect(args.database_path)
  
      idList=list()
      modelList=list()
      widthList=list()
      heightList=list()
      paramsList=list()
      # Update real cameras from .txt
      with open(txtfile, "r") as cam:
          lines = cam.readlines()
          for i in range(0,len(lines),1):
              if lines[i][0]!='#':
                  strLists = lines[i].split()
                  cameraId=int(strLists[0])
                  cameraModel=camModelDict[strLists[1]] #SelectCameraModel
                  width=int(strLists[2])
                  height=int(strLists[3])
                  paramstr=np.array(strLists[4:12])
                  params = paramstr.astype(np.float64)
                  idList.append(cameraId)
                  modelList.append(cameraModel)
                  widthList.append(width)
                  heightList.append(height)
                  paramsList.append(params)
                  camera_id = db.update_camera(cameraModel, width, height, params, cameraId)
  
      # Commit the data to the file.
      db.commit()
      # Read and check cameras.
      rows = db.execute("SELECT * FROM cameras")
      for i in range(0,len(idList),1):
          camera_id, model, width, height, params, prior = next(rows)
          params = blob_to_array(params, np.float64)
          assert camera_id == idList[i]
          assert model == modelList[i] and width == widthList[i] and height == heightList[i]
          assert np.allclose(params, paramsList[i])
  
      # Close database.db.
      db.close()
  
  if __name__ == "__main__":
      camTodatabase("created/sparse/cameras.txt")
  

  注意，在 update_camera 函数中我将 prior_focal_length 设为 True，在数据库中值为1，表示信任预先给定的焦距，如有其他需求，可查看官方文档修改。

2.2. 特征匹配

在 ProjectPath 目录下运行命令：

colmap exhaustive_matcher --database_path database.db

终端输出示例：

==============================================================================
Exhaustive feature matching
==============================================================================

Matching block [1/1, 1/1] in 13.361s
Elapsed time: 0.225 [minutes]

2.3. 三角剖分

在 ProjectPath 目录下运行命令：

colmap point_triangulator --database_path database.db --image_path images --input_path created/sparse --output_path triangulated/sparse

终端输出示例：

==============================================================================
Loading database
==============================================================================

Loading cameras... 25 in 0.000s
Loading matches... 282 in 0.005s
Loading images... 25 in 0.011s (connected 25)
Building correspondence graph... in 0.048s (ignored 0)

Elapsed time: 0.001 [minutes]
==============================================================================
Triangulating image #1
==============================================================================
  => Image has 0 / 7503 points
  => Triangulated 3807 points
...
...
Bundle adjustment report
------------------------
    Residuals : 148710
   Parameters : 38385
   Iterations : 2
         Time : 0.153904 [s]
 Initial cost : 0.291198 [px]
   Final cost : 0.289322 [px]
  Termination : Convergence

  => Merged observations: 7
  => Completed observations: 0
  => Filtered observations: 11
  => Changed observations: 0.000242
==============================================================================
Extracting colors
==============================================================================

这里注意每张图片 Triangulated Points 数，如果太少，比如小于500，就说明输入的相机外参有问题，不同视图之间找不到共同可以观察到（Triangulate）的点，重建基本是失败的。关于 Bundle adjustment 是否收敛的问题，如果不收敛，一定失败了，如果收敛了，也不一定成功。我们可以仍通过运行COLMAP 的GUI程序, File->Import model 导入刚刚生成的在 triangulated/sparse 下的.bin格式的稀疏重建模型，观察重建效果。

3. 命令行执行稠密重建

3.1. 利用 三角剖分 生成的稀疏模型进行稠密重建

在 ProjectPath 目录下运行命令：

colmap image_undistorter --image_path images --input_path triangulated/sparse --output_path dense

继续运行命令：

colmap patch_match_stereo --workspace_path dense

这一步时间较长，与图片数量，分辨率和显卡能力都有关系，25张1280*720的图片在一张GTX 1070显卡上大概需要13分钟。如果需要生成三角网格，继续运行命令：

colmap stereo_fusion --workspace_path dense --output_path dense/fused.ply

到此，稠密重建结束。如果想可视化重建的深度图，法线等，可通过运行COLMAP 的GUI程序， Reconstruction->Dense reconstruction 进入稠密重建菜单，点击 Select 选择 dense 文件夹，就可以在应用程序界面选择每张图片对应的深度图和法线图。

3.2. 使用相机Pose和图片手动稠密重建

稀疏重建（三角剖分）对于稠密重建并不是必须的。在已知相机Pose的情况下可以直接进行稠密重建，首先在 ProjectPath 目录下运行命令：

colmap image_undistorter --image_path images --input_path created/sparse --output_path dense

在 dense/sparse/ 目录下的.bin文件的内容与 created/sparse 中我们手动建立的.txt文件内容相同。在 dense/stereo/ 目录下的 patch-match.cfg 规定了源图像进行块匹配的参考图像，默认格式如下：

image001.jpg
__auto__, 20
image002.jpg
__auto__, 20
...
__auto__, 20
image025.jpg
__auto__, 20

图像名字下一行规定了稠密重建时对应图像的参考图像， __auto__, 20 表示自动优先级前20张，但该选项只有进行了稀疏重建才可用。在图像数量少的情况下，可以使用 __all__ 指定所有图像为参考图像。或者，可以使用空格分隔得图像名字指定参考图像，示例：

image001.jpg
image002.jpg, image003.jpg, image004.jpg
image002.jpg
image001.jpg, image003.jpg, image004.jpg, image007.jpg

之后，我们需要指定最小深度 depth_min 和最大深度 depth_max 来执行mvs，具体数值根据场景来定：

colmap patch_match_stereo --workspace_path dense --PatchMatchStereo.depth_min 0.0 --PatchMatchStereo.depth_max 20.0

融合3D网格的命令与上一节一样，依然为：

colmap stereo_fusion --workspace_path dense --output_path dense/fused.ply

4. 其它

部分手动得到的相机内外参，以旋转矩阵+平移向量的方式实现，colmap需要转化为四元数，这里提供一部分matlab代码：

%旋转矩阵转四元数
% R为旋转矩阵 Quat为四元数向量
q0=0.5*sqrt(1+R(1,1)+R(2,2)+R(3,3));
q1=(R(3,2)-R(2,3))/(4*q0);
q2=(R(1,3)-R(3,1))/(4*q0);
q3=(R(2,1)-R(1,2))/(4*q0);
Quat=[q0 q1 q2 q3];

已有一组相机内参矩阵 K ，相机外参四元数： Quat ，平移向量： T ，写入至 cameras.txt 和 images.txt ：

cam_txt=fopen(cam_txt_path,'w');
image_txt=fopen(image_txt_path,'w');
fprintf(cam_txt,'#By liminghao
#CAMERA_ID, MODEL, WIDTH, HEIGHT, PARAMS[]
# Number of cameras: %d 
',numofimage);
fprintf(image_txt,'#By liminghao
#IMAGE_ID, QW, QX, QY, QZ, TX, TY, TZ, CAMERA_ID, NAME
# Number of images: %d 
',numofimage);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:numofimage
    fx=Ktemp{i}(1,1);
    fy=Ktemp{i}(2,2);
    cx=Ktemp{i}(1,3);
    cy=Ktemp{i}(2,3);
    fprintf(cam_txt,'%d PINHOLE %d %d %f %f %f %f
',i,Width,Height,fx,fy,cx,cy);
    fprintf(image_txt,'%d ',i);
    for qi=1:1:length(Quatvector{i})             
        fprintf(image_txt,'%f ',Quatvector{i}(qi));
    end
    for ti=1:1:length(Ttemp{i})             
        fprintf(image_txt,'%f ',Ttemp{i}(ti));
    end
    fprintf(image_txt,'%d image%03d.jpg

',i,i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fclose(cam_txt);
fclose(image_txt);

以上代码逻辑比较简单，用Python等其他语言也较容易实现，这里要注意的就是以上代码中设置相机模型为 PINHOLE ，使用者请根据自己情况再行修改。