はじめに

この記事は、NTTドコモR&D控え室 Advent Calendar 2020 7日目の記事です。
(こちらも是非ご覧ください！NTTドコモ R&D Advent Calendar 2020)

はじめまして(1年ぶりです)、NTTドコモのkitadeと申します。
普段は画像認識技術を中心に、農業用ロボットを作ったりドローン空撮画像を解析したり、新しいUIハードを試作したり。ハードウェアからソフトウェアまで幅広く雑食に取り扱っております。

昨今の事情から今年は在宅勤務が中心となったこともあり、バッチリ対策して外出した日はガッツリ見たものを記録して帰りたいと思うようになりました(記録オタク笑)。

この記事では複数の写真を使って３次元データを再構成する、フォトグラメトリ技術を取り扱います。既に市販ソフトウェアでは多くの高精度なフォトグラメトリソフトウェアが開発されていますが、中身を知る＆カスタマイズできるようにするために、OSSのみでやってみました。

車窓からその１：観覧車の車窓から

11月のよく晴れたある日、子供にせがまれ二人でみなとみらいの観覧車に載ることになりました。みなとみらいが一望できる素敵な観覧車なのですが、乗車料金それなりにするんだよなぁ‥大人も楽しみたいなぁ‥と考えていたら思いつきました。観覧車は、街を高いところからいろいろな視点から見ることができる、フォトグラメトリにピッタリなフィールドだと。

観覧車は静かに安定して回転するため、撮影用の台車としてもとても優秀です。車窓を電車の運転席に見立ててはしゃぐ子供を横目に、客車の側面からフォトグラメトリ用画像を収集することにしました。

データ収集の方法

至って単純で、なるべく遮蔽物のない窓にスマホを立て掛け、ビデオを撮る。それだけです。今回はiPhone11Proを使い、超広角レンズを選択して4K60fps撮影をしました。なお、後述しますがフォトグラメトリの基本は画像間の中で一致する特徴点を探し出し三次元的に繋げあわせていくことですので、多くの場合画像同士の重なりが多くなる方法、すなわち広角撮影が有利です。

使ったアルゴリズム

フォトグラメトリ界のOSSデファクト的な存在であるCOLMAPをベースに、特徴点抽出と画像ペア検索にDeeplearningベースのアルゴリズム(D2-net及びSuperGlue)を使いました。
といっても、Hierarchical-Localization(hloc)というライブラリを全面的に活用してパイプラインを組んでいます。

COLMAPについては「Structure-from-Motion Revisited」という論文でその理論が説明されています。下記が論文から抜粋した、パイプライン概要です。

ここではhlocのライブラリを使って、下記２つの組み合わせでFeature Extraction及びMatchingのプロセスを置き換えています。(オリジナルのCOLMAPではFeature ExtractionにSIFTを使っているようです。)それぞれの実装については、リンク先をご参照ください。

SuperPoints + SuperGlue
D2-net + NN

Geometric Verification以降はCOLMAPのコマンドラインインタフェースを使って進めます。

環境構築と実装

今回のパイプラインではCOLMAPの利用にあたりUbuntu20.04との組み合わせで処理が止まったりしてしまう現象が見られましたので、下記のDockerfileで環境をコンテナ化することにしました。もし試される際にはご参考にしてください。
手順は、COLMAP及びhlocの公式のものを踏襲しています。

FROM nvidia/cuda:10.2-cudnn7-devel-ubuntu18.04
RUN apt-get update -y
RUN apt-get install python3 python3-pip unzip wget -y
#install ceres-solver and colmap
RUN apt-get install \
    git \
    cmake \
    build-essential \
    libboost-program-options-dev \
    libboost-filesystem-dev \
    libboost-graph-dev \
    libboost-regex-dev \
    libboost-system-dev \
    libboost-test-dev \
    libeigen3-dev \
    libsuitesparse-dev \
    libfreeimage-dev \
    libgoogle-glog-dev \
    libgflags-dev \
    libglew-dev \
    qtbase5-dev \
    libqt5opengl5-dev \
    libcgal-dev -y

RUN apt-get install libcgal-qt5-dev -y
RUN apt-get install libatlas-base-dev libsuitesparse-dev -y
RUN git clone https://ceres-solver.googlesource.com/ceres-solver
WORKDIR /ceres-solver
RUN git checkout $(git describe --tags)
RUN mkdir build
WORKDIR /ceres-solver/build
RUN cmake .. -DBUILD_TESTING=OFF -DBUILD_EXAMPLES=OFF
RUN make -j
RUN make install
WORKDIR /

RUN git clone https://github.com/colmap/colmap.git
WORKDIR /colmap
RUN git checkout dev
RUN mkdir build
WORKDIR /colmap/build
RUN cmake ..
RUN make -j
RUN make install
WORKDIR /

#install requirements for Hierarchical-Localization
RUN pip3 install scikit-build
RUN pip3 install torch numpy opencv-python tqdm matplotlib scipy h5py
RUN pip3 install jupyterlab notebook
RUN pip3 install git+https://github.com/mihaidusmanu/pycolmap

WORKDIR /

実行方法の詳細は省略しますが、hlocのサンプルpipeline_SfM.ipynbをベースに下記の処理を進めて行きます。

特徴点抽出(hloc.extract_features)
画像ペア計算(hloc.match_features)
COLMAP形式のワークスペースへ変換(hloc.reconstruction)
疎なカメラ姿勢・パラメータ推定(colmap mapper)
画像歪み補正(colmap image_undistorter)
密なカメラ姿勢・パラメータ推定(colmap patch_match_stereo)
3D点群算出(colmap stereo_fusion)

7.まで行くと、.ply形式で3Dの点群が出力されます。

結果の比較

観覧車車窓から撮影した4K60fps動画を、ffmpegで静止画に切り出してパイプラインに投入します。下記のコマンドでframestepをN=300,600,1200,2400,4800,9600、すなわち5,10,20,40,80,160秒に1枚動画切り出したデータセットを作成し、それぞれSuperPoints及びD2-netアルゴリズムを使ってパイプラインを実行しました。

ffmpeg -i 入力動画ファイル -vf framestep=N -q:v 1 IMG_XXXX_%05d.jpg

各切り出し間隔、アルゴリズムごとにマッチングに利用できた画像数とパイプラインの処理の結果得られたポイント数(点群数)を下のグラフにまとめました。

画像数がそこそこ(10-80秒間隔)の場合はSuperPoints,D2-netともに大きな差はありませんが、画像数が多かったり少なかったりする場合はSuperPointsのロバスト性が垣間見えます。
例えば、5秒・160秒間隔の画像データセットではD2-netを使ったパイプラインは3D点群算出に失敗しています。
なお、算出されたポイント数はデータセット中の画像数と相関がありそうです。ただし、データセットに含まれる画像数が多ければ多いほどよいかというとそうでもなく、ノイズが増えたり計算時間が増大するデメリットがあるため、画像間で重複する面積を考慮しながら調整する必要があります。

SuperPointsパイプラインを使って得られた点群を、同じくパイプラインから計算されたカメラ姿勢(赤い四角錐)とともに各条件ごとに表示してみました。

カメラ姿勢がきれいに観覧車の形になって表示されていますね(右下の方)。
撮影間隔が短くなるほどみなとみらいのビル群がしっかり描画されていることがわかります。先述のグラフには記載しませんでしたが、1秒間隔のデータセットをパイプラインに通した結果は点群が濃く描画されている一方でノイズが増えていることがわかります。ちなみに、1秒間隔のパターンでは処理に数日かかっています(他のパターンは長くて数時間、Xeon64スレッド、QuadroRTX8000環境)。

メッシュを作ってみる

続いて、作成された点群を使ってメッシュをデータを作ってみます。
- hloc+COLMAPによるフォトグラメトリ(.ply点群出力)
- MeshLabを使った点群のトリミング
- open3Dによるメッシュ作成

open3Dを使ったメッシュ作成についても少し解説しておきます。

まず、COLMAPのパイプラインで出力された.ply形式の点群データをopen3dを使って読み込みます。可視化のUIもあって便利です。
ちなみに点群のダウンサンプルも可能で、あまりに点群が細かすぎる状況の場合はシュリンクして取り扱うこともできます。

open3d_mesher.py

import numpy as np
import open3d as o3d
from open3d import io

pcd = o3d.io.read_point_cloud(INPUT_PLY_FILE)
#pcd = pcd.voxel_down_sample(voxel_size=0.05)
print("point cloud readed.")

o3d.visualization.draw_geometries([pcd])

最初に試すのはBall-Pivoting Algorithm（BPA）アルゴリズムです。BPAアルゴリズムでは、点群の表面にボールを転がす要領でメッシュを作り上げて行きます。ボールの大きさを決めるにあたって、前計算を行っています。シンプルで昔からあるアルゴリズムではあるのですが、調整が難しく（仕組み上穴ができやすい）、また計算時間が結構長いです。

#BPA
pcd.estimate_normals()
print("estimated normals.")
distances = pcd.compute_nearest_neighbor_distance()
print("compute nearest neighbor distance.")
avg_dist = np.mean(distances)
radius = 1.5 * avg_dist   

mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_ball_pivoting(pcd,o3d.utility.DoubleVector([radius, radius * 2]))
print("BPA mesh created. Simplifying...")
dec_mesh = mesh.simplify_quadric_decimation(100000)
dec_mesh.remove_degenerate_triangles()
dec_mesh.remove_duplicated_triangles()
dec_mesh.remove_duplicated_vertices()
dec_mesh.remove_non_manifold_edges()
print("Writing BPA mesh..")
o3d.io.write_triangle_mesh(output_path+"bpa_mesh.ply", dec_mesh)

次に試すのは、Poissonアルゴリズムです。このアルゴリズムは、点群に布をバサッとかぶせる要領でメッシュを作っていきます。depthパラメータでメッシュの解像度を高めることが可能です。ノイズはある程度ケアしてあげる必要がありますのでここではMeshlabを使ってトリミングを行いました。

#Poisson
poisson_mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(pcd, depth=11,width=0, linear_fit=False)[0]
print("Poisson mesh created.")
bbox = pcd.get_axis_aligned_bounding_box()
p_mesh_crop = poisson_mesh.crop(bbox)
print("Writing poisson mesh...")
o3d.io.write_triangle_mesh(output_path+"p_mesh_c.ply", p_mesh_crop)

下記、Poissonアルゴリズムでメッシュ化された観覧車の車窓から見える風景のメッシュデータです。

アルゴリズムの概要にもあった通り、バサっと布をかけた感じがみて取れると思います(Poissonアルゴリズム出力そのまま、編集をサボっているとも言える笑)。
実際には、もう少しノイズ処理などを工夫して綺麗なメッシュデータが作れると思います。

車窓からその2：飛行機の車窓から

続いて目をつけたのが、飛行機の窓から見える風景です。
離着陸時に見える街の風景をiPhone11Proで撮影し、フォトグラメトリしてみました。使ったアルゴリズムはSuperPoints＋SuperGlueです。
下の図がパイプラインによる処理結果の点群です。羽田空港着陸時のデータなのですが、きれいな一直線上に推定されたカメラ姿勢が表示されているのがわかります。平面が平面として取れているのも素晴らしいですね。

お気づきの方もいらっしゃるかもしれませんがこれ、新飛行経路(南風時C滑走路着陸)の経路っぽいですね。中心付近に、大井車両基地がうっすら見えていたりします。目視したときはドクターイエローまで見えたのですが、点群データにはうまく表現できていないのが残念です。とはいえ、オルソとは違ってビルなどの立体構造まで記録できるのは面白いですね。

終わりに

スマホ1台で撮影できるデータを使って、いろいろな窓から・いろいろなシーンをフォトグラメトリしてみました。高価なセンサーを使わずとも、意外と高精細な点群が得られたかと思います。
本当は生成したメッシュデータをVR空間上に配置したり、フルカラー3Dプリンタで印刷したりできれば面白かったのですが、またの機会にでも。

ちなみに、今回の記事には載せませんでしたが、この時期街中で多くみられるクリスマスツリーもこの方法で3Dモデル化できます。(社内に毎年設置されるクリスマスツリーをフォトグラメトリするつもりで準備していたら、在宅勤務体制の影響かまさかの今年は設置されず。。)
お外に出られた際は是非、試してみてはいかがでしょうか。

都会の車窓からフォトグラメトリ 〜LiDAR無しスマホでもできる！3Dスキャン〜