1：はじめに
Image Matching Challenge 2025に参加したのですが、PyColmapのドキュメンテーションがあまりにも少ないので、備忘録もかねてPythonでの動かし方を説明したいと思います。
2：SFM(Structure From Motion)とは
SfM (Structure from Motion) とは、さまざまな角度や位置から撮影された複数枚の写真から、対象物の3次元形状を復元する技術です。一般的に、ドローンやカメラで撮影した画像を解析して、3D点群データを作成する際に用いられます。また、三次元形状は以下の手順を経て復元されます。
- 特徴点抽出
- 特徴点マッチング
- 画像の３次元復元
3.1：特徴点抽出
特徴点とは、画像の中で明るさの変化が顕著に現れる場所、例えば角やエッジなどを指します。これらの特徴点を抽出することで、異なる画像間で同じ物体上の点を特定し、その位置関係を解析することで3D形状を推定することが出来ます。
詳しくは以下のリンクを参照してください。
画像から特徴点を抽出すると以下のようになります。
3.2：特徴点マッチング
抽出した特徴点を対応している者同士でマッチングします。これによって対応関係がわかります。
先ほどの画像でマッチングを実行すると以下のようになります。
このように画像中の特徴点がすべてマッチングに使われるわけではなくて、ペアの画像同士の特徴点の中で一致しているものがあればマッチングさせることになります。
その他、細かい説明などは以下が参考になると思います。
4：Pythonによる実装
実行環境：Google Colab(CPUで実行できます)
モデル：ALIKED(特徴量抽出)、Lightglue(特徴量マッチング)、Colmap(３次元構造復元)
データセット：South building (Colmap公式 https://demuc.de/colmap/datasets/ より)
4.1：データやライブラリのimport
以下のコードを実行します。pycolmapのライブラリをGithubからクローンしていますが、パスを指定できるならcolmapディレクトリで実行する必要はありません。
from google.colab import drive
drive.mount('/content/drive')
%cd /content/drive/MyDrive/library
!git clone https://github.com/cvlab-epfl/disk.git
%cd /content/drive/MyDrive/library/disk/colmap
!pip install git+https://github.com/cvg/LightGlue.git
! pip install pycolmap
from colmap.database import *
import cv2
import torch
from lightglue import match_pair
from lightglue import ALIKED, LightGlue
from lightglue.utils import load_image, rbd
import os
import pycolmap
import matplotlib.pyplot as plt
from tqdm import tqdm
import numpy as np
import gc
import os
from tqdm import tqdm
from time import time
from fastprogress import progress_bar
import h5py
from IPython.display import clear_output
from collections import defaultdict
from copy import deepcopy
import concurrent.futures
from collections import Counter
import itertools
from pathlib import Path
import kornia as K
import kornia.feature as KF
from transformers import AutoImageProcessor, AutoModel
import torch.nn.functional as F
次に、Colmapで使用する関数を定義します。Colmapは画像から抽出した特徴点、マッチング情報をdatabase形式で保存しないと実行できません。また、database形式変換を公式がライブラリを用意してくれているのですが、変換元がh5ファイルに限定されているのと、h5ファイルからの具体的な変換の仕方を公式が一切書いてくれていないので自分で関数を再定義してカスタマイズする必要があります。
以下がその関数になります。
import os, argparse, h5py, warnings
import numpy as np
from tqdm import tqdm
from PIL import Image, ExifTags
from colmap.database import COLMAPDatabase, image_ids_to_pair_id
#Copyright [2020] [Michał Tyszkiewicz, Pascal Fua, Eduard Trulls]
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
def get_focal(image_path, err_on_default=False):
image = Image.open(image_path)
max_size = max(image.size)
exif = image.getexif()
focal = None
if exif is not None:
focal_35mm = None
# https://github.com/colmap/colmap/blob/d3a29e203ab69e91eda938d6e56e1c7339d62a99/src/util/bitmap.cc#L299
for tag, value in exif.items():
focal_35mm = None
if ExifTags.TAGS.get(tag, None) == 'FocalLengthIn35mmFilm':
focal_35mm = float(value)
break
if focal_35mm is not None:
focal = focal_35mm / 35. * max_size
if focal is None:
if err_on_default:
raise RuntimeError("Failed to find focal length")
# failed to find it in exif, use prior
FOCAL_PRIOR = 1.2
focal = FOCAL_PRIOR * max_size
return focal
def create_camera(db, image_path, camera_model):
image = Image.open(image_path)
width, height = image.size
focal = get_focal(image_path)
if camera_model == 'simple-pinhole':
model = 0 # simple pinhole
param_arr = np.array([focal, width / 2, height / 2])
if camera_model == 'pinhole':
model = 1 # pinhole
param_arr = np.array([focal, focal, width / 2, height / 2])
elif camera_model == 'simple-radial':
model = 2 # simple radial
param_arr = np.array([focal, width / 2, height / 2, 0.1])
elif camera_model == 'opencv':
model = 4 # opencv
param_arr = np.array([focal, focal, width / 2, height / 2, 0., 0., 0., 0.])
return db.add_camera(model, width, height, param_arr)
def add_keypoints(db, h5_path, image_path, img_ext, camera_model, single_camera = False):
keypoint_f = h5py.File(os.path.join(h5_path, 'keypoints.h5'), 'r')
camera_id = None
fname_to_id = {}
for filename in tqdm(list(keypoint_f.keys())):
keypoints = keypoint_f[filename][()]
fname_with_ext = filename# + img_ext
path = os.path.join(image_path, fname_with_ext)
if not os.path.isfile(path):
raise IOError(f'Invalid image path {path}')
if camera_id is None or not single_camera:
camera_id = create_camera(db, path, camera_model)
image_id = db.add_image(fname_with_ext, camera_id)
fname_to_id[filename] = image_id
db.add_keypoints(image_id, keypoints.squeeze(0))
return fname_to_id
def add_matches(db, h5_path, fname_to_id):
match_file = h5py.File(os.path.join(h5_path, 'matches.h5'), 'r')
added = set()
n_keys = len(match_file.keys())
n_total = (n_keys * (n_keys - 1)) // 2
with tqdm(total=n_total) as pbar:
for key_1 in match_file.keys():
group = match_file[key_1]
for key_2 in group.keys():
id_1 = fname_to_id[key_1]
id_2 = fname_to_id[key_2]
pair_id = image_ids_to_pair_id(id_1, id_2)
if pair_id in added:
warnings.warn(f'Pair {pair_id} ({id_1}, {id_2}) already added!')
continue
matches = group[key_2][()]
db.add_matches(id_1, id_2, matches)
added.add(pair_id)
pbar.update(1)
4.2：ディレクトリの指定
以下のようにディレクトリを指定します。この辺りは好みなので、好みに合わせて設定してください。
ROOT_PATH = '/content/drive/MyDrive/dataset/south-building'
IMAGE_PATH = '/content/drive/MyDrive/dataset/south-building/images/'
device = torch.device('cuda') if torch.cuda.is_available() else torch.device('cpu')
os.makedirs(f'{ROOT_PATH}/output', exist_ok=True)
OUT_PATH = f'{ROOT_PATH}/output'
DEBUG = True
4.3：少数の画像で試してみる
ここで、二枚の画像を使用して、特徴点抽出とマッチングをやってみます。
以下のコードで実行します。EXTRACTERは特徴点抽出を、MATCHERは特徴点マッチングを行います。
ALIKEDはn16が、lightglueはv0-1がダウンロードされると思います。
EXTRACTOR = ALIKED(max_num_keypoints=4096,detection_threshold=0.01).eval().to(device)
MATCHER = LightGlue(features='aliked').eval().to(device) # load the matcher
Image1 = load_image(f'{IMAGE_PATH}/P1180160.JPG').to(device)
Image2 = load_image(f'{IMAGE_PATH}/P1180157.JPG').to(device)
with torch.inference_mode():
feats0, feats1, matches01 = match_pair(EXTRACTOR, MATCHER, Image1, Image2)
from lightglue.utils import load_image, rbd
from lightglue import viz2d
kpts0, kpts1, matches = feats0["keypoints"], feats1["keypoints"], matches01["matches"]
m_kpts0, m_kpts1 = kpts0[matches[..., 0]], kpts1[matches[..., 1]]
axes = viz2d.plot_images([Image1, Image2])
viz2d.plot_matches(m_kpts0, m_kpts1, color="lime", lw=0.2)
viz2d.add_text(0, f'Stop after {matches01["stop"]} layers', fs=20)
kpc0, kpc1 = viz2d.cm_prune(matches01["prune0"]), viz2d.cm_prune(matches01["prune1"])
viz2d.plot_images([Image1, Image2])
viz2d.plot_keypoints([kpts0, kpts1], colors=[kpc0, kpc1], ps=10)
これを実行すると、以下の画像が生成されます。
4.4：Global featureを用いた類似度測定
Kaggleでは割とメジャーな手法ですが、すべての画像でマッチングをすると計算量が多くなるなどの課題があるので、画像の類似度を測定してマッチングに使う画像の数を減らします。
モデルはDINOv2を使用します。DinoV2はhuggin faceのfacebookディレクトリに入っています。これをロードするためにTransformersからAutoImageProcessorとAutoModelをImportしています。
ここで、マッチングに使う画像が一定数以下ならすべてのマッチを使用します。また、計算時間の関係でデータセット中のすべての画像は使っていません。以下の関数では画像のペアのindexリストを返します。
def get_global_desc(fnames, device = torch.device('cpu')):
processor = AutoImageProcessor.from_pretrained('facebook/dinov2-base')
model = AutoModel.from_pretrained("facebook/dinov2-base")
model = model.eval()
model = model.to(device)
global_descs_dinov2 = []
for i, img_fname_full in tqdm(enumerate(fnames),total= len(fnames)):
key = os.path.splitext(os.path.basename(img_fname_full))[0]
timg = load_image(img_fname_full)
with torch.inference_mode():
inputs = processor(images=timg, return_tensors="pt", do_rescale=False).to(device)
outputs = model(**inputs)
dino_mac = F.normalize(outputs.last_hidden_state[:,1:].max(dim=1)[0], dim=1, p=2)
global_descs_dinov2.append(dino_mac.detach().cpu())
global_descs_dinov2 = torch.cat(global_descs_dinov2, dim=0)
return global_descs_dinov2
def get_img_pairs_exhaustive(img_fnames):
index_pairs = []
for i in range(len(img_fnames)):
for j in range(i+1, len(img_fnames)):
index_pairs.append((i,j))
return index_pairs
def get_image_pairs_shortlist(fnames,
sim_th = 0.8, # should be strict
min_pairs = 20,
exhaustive_if_less = 15,
device=torch.device('cpu')):
num_imgs = len(fnames)
if num_imgs <= exhaustive_if_less:
return get_img_pairs_exhaustive(fnames)
descs = get_global_desc(fnames, device=device)
dm = torch.cdist(descs, descs, p=2).detach().cpu().numpy()
# removing half
mask = dm <= sim_th
total = 0
matching_list = []
ar = np.arange(num_imgs)
already_there_set = []
for st_idx in range(num_imgs-1):
mask_idx = mask[st_idx]
to_match = ar[mask_idx]
if len(to_match) < min_pairs:
to_match = np.argsort(dm[st_idx])[:min_pairs]
for idx in to_match:
if st_idx == idx:
continue
if dm[st_idx, idx] < 200:
matching_list.append(tuple(sorted((st_idx, idx.item()))))
total+=1
matching_list = sorted(list(set(matching_list)))
return matching_list
DebagをTrueにしていたら、マッチングが実行されるようにします。この結果120個の画像のペアが返されました。
if DEBUG:
image_list = list(Path(IMAGE_PATH).glob('*.JPG'))[::5]
print(image_list)
pairs = get_image_pairs_shortlist(image_list)
print(pairs)
print(len(pairs))
4.5：KeypointsとDescriptorsを求める
先ほど定義したEXTRACTERを使って特徴点を求めます。先ほどの参考資料に明記されていますが、
画像上で物体の特徴を決める上では、従来から2段階の方法をとることが多いです。 これは前回記事でも紹介したとおり以下のような流れで、1段階目を特徴点抽出(Detector)、2段階目(Descriptor)を特徴量記述といいます。 この2段階をまとめて特徴量抽出(Feature Extraction)と呼びます。
同様にDEBAGがTrueの時に実行されるようにします。
この時、それぞれの情報はkeypoints.h5とdescritor.h5に書き出されます。
def detect_keypoints(paths, feature_dir=OUT_PATH):
feature_dir_path = Path(feature_dir)
with h5py.File(f'{feature_dir}/keypoints.h5', mode="w") as f_keypoints, h5py.File(f'{feature_dir}/descriptors.h5', mode="w") as f_descriptors:
for path in tqdm(paths, desc='computing keypoints....'):
key = path.name
with torch.inference_mode():
image = load_image(path).to(device)
feats = EXTRACTOR.extract(image)
f_keypoints[key] = feats["keypoints"].cpu().numpy()
f_descriptors[key] = feats["descriptors"].cpu().numpy()
print('computing finished!!')
if DEBUG:
detect_keypoints(image_list)
4.6：マッチング情報を求める
matcherとKeypointsとDescriptorsを使ってマッチング情報を求めます。
これはmatches.h5に書き出されます。
def keypoint_distances(
paths: list[Path],
index_pairs: list[tuple[int, int]],
feature_dir: Path,
min_matches: int = 15,
verbose: bool = True,
device: torch.device = torch.device("cpu"),
):
matcher = KF.LightGlueMatcher("aliked").eval().to(device)
with h5py.File(f'{OUT_PATH}/keypoints.h5', mode="r") as f_keypoints,h5py.File(f'{OUT_PATH}/descriptors.h5', mode="r") as f_descriptors,h5py.File(f'{OUT_PATH}/matches.h5', mode="w") as f_matches:
for idx1, idx2 in tqdm(index_pairs, desc="Computing keypoing distances"):
key1, key2 = paths[idx1].name, paths[idx2].name
keypoints1 = torch.from_numpy(f_keypoints[key1][...]).to(device).squeeze(0)
keypoints2 = torch.from_numpy(f_keypoints[key2][...]).to(device).squeeze(0)
descriptors1 = torch.from_numpy(f_descriptors[key1][...]).to(device).squeeze(0)
descriptors2 = torch.from_numpy(f_descriptors[key2][...]).to(device).squeeze(0)
with torch.inference_mode():
distances, indices = matcher(
descriptors1,
descriptors2,
KF.laf_from_center_scale_ori(keypoints1[None]),
KF.laf_from_center_scale_ori(keypoints2[None]),
)
# We have matches to consider
n_matches = len(indices)
if n_matches:
if verbose:
print(f"{key1}-{key2}: {n_matches} matches")
# Store the matches in the group of one image
if n_matches >= min_matches:
group = f_matches.require_group(key1)
group.create_dataset(key2, data=indices.detach().cpu().numpy().reshape(-1, 2))
if DEBUG:
keypoint_distances(image_list, pairs, OUT_PATH, verbose=False)
2個の画像のマッチ情報を知りたいときは以下を実行します。
if DEBUG:
with h5py.File(f'{OUT_PATH}/keypoints.h5', mode="r") as f_keypoints, h5py.File(f'{OUT_PATH}/descriptors.h5', mode="r") as f_descriptors:
matcher_params = {
"width_confidence": -1,
"depth_confidence": -1,
"mp": True if 'cuda' in str(device) else False,
}
matcher = KF.LightGlueMatcher("aliked", matcher_params).eval().to(device)
idx1, idx2 = pairs[0]
key1, key2 = image_list[idx1].name, image_list[idx2].name
keypoints1 = torch.from_numpy(f_keypoints[key1][...]).to(device).squeeze(0)
keypoints2 = torch.from_numpy(f_keypoints[key2][...]).to(device).squeeze(0)
print("Keypoints:", keypoints1.shape, keypoints2.shape)
descriptors1 = torch.from_numpy(f_descriptors[key1][...]).to(device).squeeze(0)
descriptors2 = torch.from_numpy(f_descriptors[key2][...]).to(device).squeeze(0)
print("Descriptors:", descriptors1.shape, descriptors2.shape)
with torch.inference_mode():
distances, indices = matcher(
descriptors1,
descriptors2,
KF.laf_from_center_scale_ori(keypoints1[None]),
KF.laf_from_center_scale_ori(keypoints2[None]),
)
print(distances, indices)
4.7：Colmapを使った3次元構造推定
ここまでお疲れさまでした。これが最後のステップです。
COLMAPdatabaseでdatabaseを作成してpycolmap.incremental_mappingから３次元構造の推定を実行します。
def import_into_colmap(path = IMAGE_PATH ,feature_dit = OUT_PATH , detabase_path = f'{OUT_PATH}/colmap.db'):
db = COLMAPDatabase.connect(detabase_path)
db.create_tables()
single_camera = False
fname_to_id = add_keypoints(db ,OUT_PATH,IMAGE_PATH ,'','simple-pinhole',single_camera)
add_matches(db ,OUT_PATH , fname_to_id)
db.commit()
db.close()
if DEBUG:
import_into_colmap()
# This does RANSAC
pycolmap.match_exhaustive(f'{OUT_PATH}/colmap.db')
mapper_options = pycolmap.IncrementalPipelineOptions()
mapper_options.min_model_size = 3
mapper_options.max_num_models = 2
maps = pycolmap.incremental_mapping(
database_path=f'{OUT_PATH}/colmap.db',
image_path=IMAGE_PATH,
output_path=OUT_PATH,
options=mapper_options,
)
rotationなどの情報が欲しい場合は以下のコードを実行してください。
if DEBUG:
print(maps[0].summary())
for k, im in maps[0].images.items():
print("Rotation", im.cam_from_world.rotation.matrix(), "Translation:", im.cam_from_world.translation, sep="\n")
print()
4.8：3次元構造の図示
最後に推定した構造の図示をします。3次元に図示するためにpycolmap.Reconstructionを読み込んでplyファイルに書き換えて3次元で図示します。
result = pycolmap.Reconstruction(f'{OUT_PATH}/0')
result.export_PLY(f'{OUT_PATH}/0/rec.ply')
!pip install trimesh scipy pyglet
import trimesh
import matplotlib.pyplot as plt
# PLYファイルのパスを指定
file_path = f'{OUT_PATH}/0/rec.ply' # PLYファイルのパスを指定
# 点群データの読み込み
point_cloud = trimesh.load(file_path)
# 点群データの表示
fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(111, projection='3d')
# 点群の座標を抽出してプロット
ax.scatter(point_cloud.vertices[:, 0],
point_cloud.vertices[:, 1],
point_cloud.vertices[:, 2],
s=0.1) # s=0.1 で点の大きさを設定
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('3D Point Cloud')
plt.show()
これを実行すると以下のような画像が出力されます。
マウスを使って詳細な3次元構造を確認したい場合は以下のコードを実行してください。
import plotly.graph_objects as go
import trimesh
# PLYファイルのパスを指定
file_path = f'{OUT_PATH}/0/rec.ply'
# 点群データを読み込みます
try:
point_cloud = trimesh.load(file_path)
# 点群の座標と色を取得します
points = point_cloud.vertices
# colors = point_cloud.colors # 元の色は使用しない
# Plotlyの3D散布図を作成します
fig = go.Figure(data=[go.Scatter3d(
x=points[:, 0],
y=points[:, 1],
z=points[:, 2],
mode='markers',
marker=dict(
size=1, # 点のサイズを調整できます
color='blue', # 点の色を赤に変更
# colorscale='Viridis', # カラーマップは不要
opacity=0.8
)
)])
# レイアウトを設定します
fig.update_layout(
margin=dict(l=0, r=0, b=0, t=0),
scene=dict(
xaxis_title='X',
yaxis_title='Y',
zaxis_title='Z'
),
width=400, # 図の幅を小さくする例
height=400 # 図の高さを小さくする例
)
# プロットを表示します
fig.show()
except FileNotFoundError:
print(f"エラー: ファイル '{file_path}' が見つかりません。")
except Exception as e:
print(f"エラーが発生しました: {e}")
これを実行すると以下のようにマウスを使って三次元構造をいろいろな角度から見られるようになります。
データ数がそこまで多くないので正確には推定できていませんが、まずまずの仕上がりと思います。
さいごに
２度とやりたくない！！！ドキュメントがなさすぎる！！！
Discussion