3D Reconstruction (3次元再構成) についての概要とCOLMAPを用いた実践

3Dのデータをいくつか取り扱っている中で、3D Reconstruction という技術に当たりました。ここでは、大まかにどのような技術があるのかという意味での概要と、具体的な手法の1つとしてSfM・MVSを用いた実践を1つの記事としてまとめてみました。

## 3D Reconstruction (3次元再構成) とは
* コンピュータービジョンの一分野で2Dの画像やその他のセンサーデータなどを用いて現実世界の物体の3次元的な情報(位置・形状など)を取得する手法のこと
* ロボティクス・自動運転・医療画像診断・XR(AR/VR) など様々な分野で応用されている

## データの取得方法から見た手法の分類
大きくActiveな手法とPassiveな手法に分類される。

### Activeな手法
測定機器のほうからActive(能動的)に信号を送信し、その反射や干渉データを基に物体との距離を推定する方法。RGB画像に加えて深度マップが付与されたRGB-Dデータを取得することができる。

* レーザー
  LiDAR(Light Detection and Ranging)などに代表される、ビームを照射して反射時間(Time of Flight, ToF)を計測する
* 構造化光 (structured light)
  物体に対し、特定のパターン(ストライプ状など)の光を当て、その歪みを測定する方法

Passiveな方法と比べて精度がよかったり、撮影条件(照明・テクスチャ)の制約が少なかったりする一方、一般的には機器が高価であることが多い。
しかしながら、最近は例えば[iPhoneやiPadの一部機種で搭載されていたり](https://support.apple.com/ja-jp/102468#ipad)と、以前よりは気軽になってきているように思われる。

### Passiveな手法
Activeと違って、機器側から信号を送信せず、物体から反射した自然光を利用して(要するに普通の写真を使って)3次元情報を推定する方法。

* 単眼(1枚の画像)
  1枚のRGB画像から3次元の情報を幾何学的な情報だけで求めることはill-definedな問題(方程式の数が足りない)であり、付加的に情報を取得したり、推定したりする必要がある。

  従来的な方法として、ピント(カメラの焦点を変えて、ボケ具合を見る)や寸法が既知の物体などを利用する方法がある。

  <span></span>

* 複眼(複数の画像)
  人間を含む多くの動物は複数の目からの見え方の違いを利用して、物体との距離を認識している。複数のカメラを使ったり、カメラの位置を変えたりして得られた複数の視点の画像から幾何学的に物体とカメラの位置を計算して求める。

  <a href="https://en.wikipedia.org/wiki/File:Aufnahme_mit_zwei_Kameras.svg">
    <img src="https://upload.wikimedia.org/wikipedia/commons/thumb/7/7c/Aufnahme_mit_zwei_Kameras.svg/1920px-Aufnahme_mit_zwei_Kameras.svg.png" width="400"/>
  </a>

  Structure from Motion (SfM)・Multi-View Stereo (MVS) を用いた[実践](#colmap-%E3%82%92%E7%94%A8%E3%81%84%E3%81%9F%E5%AE%9F%E8%B7%B5)を後述する。

また、最近では機械学習的な方法を用いて推定する方法が多く出されている。この記事では深くは取り扱わないが、この分野の進展において近年急速に進展している領域である。有名な例としては[NeRF](https://arxiv.org/abs/2003.08934)や[3DGS](https://arxiv.org/abs/2308.04079)が挙げられる。

## Sparse・Dense な Reconstruction
出力されるデータの種類に応じて、Sparse(疎)・Dense(密)なReconstructionとして分類される。

* SparseなReconstruction
  画像の中からエッジ(辺)やコーナー(角)などの特徴点のみに基づいてそれぞれの点の対応関係からReconstruction を行うもの
* DenseなReconstruction
  画像全体を解析し、ピクセルごとに深度などの情報を得るもの

## Structure from Motion (SfM)
異なる視点から撮影された複数の画像から物体とカメラの3Dの位置を同時に推定するSparse Reconstructionの手法。

1. 特徴点の抽出
  各画像から特徴点を抽出する。[SIFT](https://en.wikipedia.org/wiki/Scale-invariant_feature_transform)(Scale-Invariant Feature Transform), [SURF](https://en.wikipedia.org/wiki/Speeded_up_robust_features)(Speeded-Up Robust Features), [ORB](https://en.wikipedia.org/wiki/Oriented_FAST_and_rotated_BRIEF)(Oriented FAST and Rotated BRIEF)などが知られている。
2. 特徴点のマッチング
  特徴点は多次元の特徴量で表現されている。2枚の画像の特徴量が似ている点をマッチングさせる。誤対応を取り除くためにRANSACなどの手法を用いる。
3. 3Dの位置推定
  画像間の対応点をもとに三角測量(triangulation)を利用して、特徴点の3D座標を計算する。

## Multi-View Stereo (MVS)
SfMで得られたカメラとsparseな特徴点情報をもとに、denseな3Dデータを得る方法。

1. denseな深度マップの計算
  2つ以上の画像の特徴点の対応関係をもとに、各ピクセルの深度を計算する。
2. 複数の画像の深度マップをもとにした点群の作成
  3次元空間上の点群データを作成する。

## ツール
3D reconstruction を行えるツールがフリーのものでいくつかある。
* [COLMAP](https://github.com/colmap/colmap)
* [OpenSfM](https://opensfm.org/)
* [VisualSfM](http://ccwu.me/vsfm/)
* [Meshroom](https://github.com/alicevision/Meshroom)
* [OpenMVG](https://github.com/openMVG/openMVG)

ただし、いくつかのツールは数年間大きな更新がないものもある。この記事では現在(2024年)も開発が続けられているCOLMAPを使用して3D reconstruction の実践を行った。

## [COLMAP](https://github.com/colmap/colmap) を用いた実践

### インストール

[Installation](https://colmap.github.io/install.html)を参照。Windows環境であれば、左記リンク先に記述されているように[こちら](https://demuc.de/colmap/)から Pre-release または Release バイナリをダウンロードする。`colmap.exe` が実行できるようにパスを通しておくとよい。以下、Windows環境を想定する。

### Automatic Reconstruction
`COLMAP.bat` を実行するとGUIが起動する。メニューバーの Reconstruction > Automatic reconstruction をクリックすると、フォルダーといくつかのパラメーターの指定より、SfM・MVSの各処理をつなげたパイプライン処理が自動で実行される。

`colmap automatic_reconstructor` により、コマンドベースでも実行できる。

### 個別の処理の実行
上述したパイプライン処理は個別に実行できる。それぞれの処理の内容と生成物を確認する。

準備として、作業用のフォルダーとその下に `images` フォルダーを作成し、 `images` フォルダー内に使用する画像すべてを置く(ファイル名は任意)。

```
/path/to/project
└── images/
    ├── image1.jpg
    ├── image2.jpg
    ├── ...
```

以下のコマンドは `/path/to/project` で実行することを想定する(`cd /path/to/project`)

1. 特徴点の抽出
    ```powershell
    colmap feature_extractor --database_path database.db --image_path images
    ```
    SQLiteのファイル `database.db` が作成され、画像ファイルとその特徴点の情報などが保存される。

    <span></span>

2. 特徴点のマッチング
    ```powershell
    colmap exhaustive_matcher --database_path database.db
    ```
    1.で抽出した特徴点のマッチングを行う。同じデータベース(`database.db`)にマッチングの情報が保存される。exhaustive_matcher 以外にもいくつかのmatcherがあり、撮影条件に応じて適切なものを選ぶとよい([参考](https://colmap.github.io/tutorial.html#feature-matching-and-geometric-verification))。

    <span></span>

3. Sparse 3D reconstruction
    ```powershell
    mkdir sparse
    colmap mapper --database_path database.db --image_path images --output_path sparse
    ```
    1.2.のデータをもとに SfM による Sparse 3D reconstruction を行い、カメラと特徴点の3D空間上の座標を計算する。
     
    * `cameras.bin`, `images.bin`, `points3D.bin` などが作成される。これらはバイナリファイルであるが、以下によりテキストファイル化することができる。
        ```powershell
        mkdir sparse_txt/0
        colmap model_converter --input_path sparse/0 --output_path sparse_txt/0 --output_type TXT
        ```

        * 出力例
        ```cameras.txt
        # Camera list with one line of data per camera:
        #   CAMERA_ID, MODEL, WIDTH, HEIGHT, PARAMS[]
        # Number of cameras: 1
        1 SIMPLE_RADIAL 4032 3024 3296.67251130364 2016 1512 0.017093516036745148
        ```

        ```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: 78, mean observations per image: 670.16666666666663
        65 0.55891783123726846 0.022128111277298114 0.78666459580772718 0.26130445522224305 -3.2328704066921916 0.028955970486720498 0.19995538039688623 1 IMG_0183.JPG
        2241.623779296875 13.273831367492676 -1 2241.416748046875 25.041868209838867 -1 2036.967529296875 89.324897766113281 -1 2035.1534423828125 96.087562561035156 -1 538.929443359375 117.58818054199219 -1 2008.7135009765625 152.4844970703125 -1 1902.857421875 155.88566589355469 -1 1999.121826171875 163.20806884765625 -1 
        ...
        ```

        ```points3D.txt
        # 3D point list with one line of data per point:
        #   POINT3D_ID, X, Y, Z, R, G, B, ERROR, TRACK[] as (IMAGE_ID, POINT2D_IDX)
        # Number of points: 13372, mean track length: 3.9091384983547712
        944 -3.0240020883344054 -0.37937594844595512 3.7383005881167892 207 206 207 2.1946228533193333 12 11507 77 8709 59 7156 14 9123 56 12883 57 9873 13 9854 67 5626 66 11246 79 16150
        1 -1.9330339681318982 -0.75438385790511264 4.6082875792386151 130 130 106 1.5793058294304585 31 307 11 5768 12 9835 79 12845 78 5009 77 4382 32 9851 30 5536 26 10988 13 9275 14 8693 56 12514 62 8781
        ...
        ```

    <span></span>
     
4. 画像の歪みの除去
    ```powershell
    mkdir dense
    colmap image_undistorter --image_path images --input_path sparse/0 --output_path dense
    ```
    MVS を行うために必要な歪みを除去した画像を作成する。
    * 画像は `dense\images` に保存される
    * `dense\sparse` に3.で作成されたファイルのコピーが保存される

    <span></span>

5. Dense 3D reconstruction
    ```powershell
    colmap patch_match_stereo --workspace_path dense --workspace_format COLMAP --PatchMatchStereo.geom_consistency true
    ```
    MVS による dense な 3D reconstruction を行う。`dense\stereo` 以下に
    * `depth_maps` ... 深度マップ
    * `normal_maps` ... 法線マップ
    
    のフォルダーが作成され、その中にbinファイルが作成される。

    ```
    /dense/stereo/
    ├── depth_maps
    │   ├── image1.jpg.geometric.bin
    │   ├── image1.jpg.photometric.bin
    │   ├── image2.jpg.geometric.bin
    │   ├── image2.jpg.photometric.bin
    │   ├── ...
    ├── normal_maps
    │   ├── image1.jpg.geometric.bin
    │   ├── image1.jpg.photometric.bin
    │   ├── image2.jpg.geometric.bin
    │   ├── image2.jpg.photometric.bin
    │   ├── ...
    ```

    このbinファイルの形式はやや特殊で、先頭に `width&height&channel&` (`width`、`height`、`channel` は整数)がテキスト形式で、その後にfloat32 `width` × `height` × `channel` の配列としてバイナリ形式のデータ部分が続く。
    
    * 深度マップを画像として表示する例(python):
        ```python
        import numpy as np
        import matplotlib.pyplot as plt
        
        filepath = '*****.bin'
        
        with open(filepath, 'rb') as f:
    
            # 3つ目の「&」までをヘッダーとして取得する 
            header = b''
            i = 0
            while i < 3:
                c = f.read(1)
                header += c
                
                if c == b'&':
                    i += 1
            
            header = header.decode('utf-8').strip("&")
            width, height, channel = map(int, header.split("&"))
        
            data = np.frombuffer(f.read(), dtype=np.float32)
        
        image = data.reshape((height, width, channel))
    
        plt.imshow(image) 
        plt.colorbar()
        plt.show()
        ```

    * depth map、 `***.geometric.bin` の出力例

        ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3627611/6e0c4203-7c68-e167-e5de-199353473fb4.png)
    
    <span></span>

6. 深度マップの融合
    ```powershell
    colmap stereo_fusion --workspace_path dense --workspace_format COLMAP --input_type geometric --output_path dense/fused.ply
    ```
    5.のデータを融合して、1つの3D空間上にある点群データとしてまとめる。`fused.ply` に点群データとしてデータが保存される。このplyファイルには各点の座標(`x`, `y`, `z`)、法線ベクトル(`nx`, `ny`, `nz`)、色(`red`, `green`, `blue`)の情報が格納されている。

    <span></span>

7. メッシュデータの作成
    ```powershell
    colmap poisson_mesher --input_path dense/fused.ply --output_path dense/meshed-poisson.ply
    ```
    and/or
    ```powershell
    colmap delaunay_mesher --input_path dense --output_path dense/meshed-delaunay.ply
    ```
    COLMAPには点群をメッシュ化するアルゴリズムとして上の2つがある。メッシュデータを含んだplyファイルが出力される。出力例は下の[実践 - 実行結果](#%E5%AE%9F%E8%A1%8C%E7%B5%90%E6%9E%9C)を参照。

### 画像データ
出力結果が画像データの質に依存する。[COLMAPのチュートリアル](https://colmap.github.io/tutorial.html#structure-from-motion)のガイドラインとして、

* 質感のある画像を撮影する
  * 白い壁や何もない机のような、完全に質感がない画像は避ける
  * 必要なら背景にポスターなどを追加して質感を増やす

<span></span>

* 照明条件を統一する
  * 明暗差が激しいシーン（逆光、影、窓越しの撮影など）は避ける
  * 光沢面の反射（スペキュラリティ）を避ける

<span></span>

* 画像に十分な視覚的重複を持たせる
  * 各オブジェクトが最低3枚以上の画像に映るようにする
  * 画像が多ければ多いほど良いが、必要以上に多すぎると処理が遅くなる可能性がある

<span></span>

* 異なる視点から撮影する
  * 同じ場所でカメラを回転させるだけではなく、撮影位置を少しずつ変える

が挙げられている。

### 実践

以下の条件で3Dオブジェクトの生成を試みた。

* 入力データ ... iPhone SE で撮影した 4032x3024 のJPG画像79枚(481MB)。

    以下はそのうちの2枚(弊社の受付)。

    <span></span>
    
    ![IMG_0150.JPG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3627611/c03cc7ce-de0a-8771-448f-0c54bec84073.jpeg)

    <span></span>

* マシンスペック
    |||
    |:---|:---|
    |OS|Windows 11 Pro 24H2|
    |CPU|CPU Intel(R) Core(TM) i9-9900K CPU @ 3.60GHz 8コア|
    |RAM|32GB|
    |GPU|NVIDIA GeForce RTX 2080 Ti|

#### 実行結果
* [個別の処理の実行](#個別の処理の実行) 1.～7.までの処理を約1時間で実行。うち55分は5. (Dense 3D reconstruction)にかかった。
* 各処理によって出力されたファイルのサイズ
  * `database.db` 133MB
  * `sparse` 21MB
  * `dense` 7.3 GB ... うち `depth_maps` 1.7GB, `normal_maps` 5.2 GB

* メッシュの表示のためにpythonのライブラリ[Open3D](https://www.open3d.org/)を用いた。
    ```python
    import open3d as o3d

    filepath = 'path/to/project/meshed-poisson.ply'
    
    mesh = o3d.io.read_triangle_mesh(filepath)
    o3d.visualization.draw_geometries([mesh])
    ```

    ![Animation.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3627611/98132d1b-816c-e3c6-4230-592ad8362ddb.gif)

    白抜きになっている部分はおおむね、撮影枚数が不足しているかlow-textureな部分であり、そうでない箇所についてある程度の質で再現できているように見える。

## 最後に
3D Reconstruction の方法は記事の前半にも記載したようにいくつかの方法がある。実践では比較的古典的な手法を用いたが、機械学習や生成AI的な手法は今後さらに成長すると予想される。また、LiDARなどによる手法ももっと気軽にできるようになるかもしれないが、さまざまなシーンや撮影条件に合わせてこれらの使い分けが必要になってくるように思われる。