概要
画像レジストレーションを例に、ホモグラフィ行列を計算する方法を紹介します。
画像レジストレーション
画像レジストレーションとは、2つの画像を位置合わせすることであり、コンピュータビジョンや画像処理の分野で広く使用されています。
ホモグラフィ行列は、2つの平面の間の射影変換を表す行列であり、画像レジストレーションに使用されます。ホモグラフィ行列を使用することにより、2つの画像の座標系を変換することができ、2つの画像を正確に重ね合わせることができます。
ホモグラフィ行列は、最低4つの点の対応関係を使用して計算されます。つまり、1つの画像内の4つの点ともう1つの画像内の対応する4つの点を知っている必要があります。この情報を使用して、ホモグラフィ行列を計算し、2つの画像を正確に重ね合わせることができます。
ホモグラフィ行列の計算を解説します。ある対象物を2台のカメラで撮影した場合を考えます。それぞれのカメラで撮影された画像をAおよびBとします。画像間の対応する点の一つを$\boldsymbol{x}_A, \boldsymbol{x}_B$とします。
\boldsymbol{x}_A = \begin{pmatrix}
x_A \\ y_A
\end{pmatrix},
\boldsymbol{x}_B = \begin{pmatrix}
x_B \\ y_B
\end{pmatrix}.
射影変換は座標に1次元加えることで、すっきりとした式で書くことができます。座標値を1次元拡張し、値を$1$とします。
\tilde{\boldsymbol{x}}_A = \begin{pmatrix}
x_A \\ y_A \\ 1
\end{pmatrix},
\tilde{\boldsymbol{x}}_B = \begin{pmatrix}
x_B \\ y_B \\ 1
\end{pmatrix}.
ホモグラフィ行列を$(H)_{ij}=h_{ij}$とすると、射影変換は次のように書けます。
s \begin{pmatrix}
x_A \\ y_A \\ 1
\end{pmatrix}
= \begin{pmatrix}
h_{11} & h_{12} & h_{13} \\
h_{21} & h_{22} & h_{23} \\
h_{31} & h_{32} & 1
\end{pmatrix}
\begin{pmatrix}
x_B \\ y_B \\ 1
\end{pmatrix}
.
$s$は$h_{31},h_{32}$で表すことができるので、ホモグラフィ行列の自由度は8です。画像間の対応する点1つにつき方程式が2本立つので、最低4点あればホモグラフィ行列を求めることができます。4点より多くある場合は、最小二乗法等でホモグラフィ行列を求めることができます。本解説では、最小二乗法を使用してホモグラフィ行列を求めます。
式を展開して整理し直すと、次のように書き換えられます。
\begin{pmatrix}
x_B & y_B & 1 & 0 & 0 & 0 & -x_A x_B &- x_A y_B \\
0 & 0 & 0 & x_B & y_B & 1 & -y_A x_B &- y_A y_B
\end{pmatrix}
\begin{pmatrix}
h_{11} \\
h_{12} \\
h_{13} \\
h_{21} \\
h_{22} \\
h_{23} \\
h_{31} \\
h_{32}
\end{pmatrix}
=
\begin{pmatrix}
x_A \\
y_A \\
\end{pmatrix}
.
画像間の対応する点が$N$個ある場合、記号を整理して次のように書きます。
X = \begin{pmatrix}
x_B^{1} & y_B^{1} & 1 & 0 & 0 & 0 & -x_A^{1} x_B^{1} &- x_A^{1} y_B^{1} \\
0 & 0 & 0 & x_B^{1} & y_B^{1} & 1 & -y_A^{1} x_B^{1} &- y_A^{1} y_B^{1} \\
\cdots \\
\cdots \\
x_B^{N} & y_B^{N} & 1 & 0 & 0 & 0 & -x_A^{N} x_B^{N} &- x_A^{N} y_B^{N} \\
0 & 0 & 0 & x_B^{N} & y_B^{N} & 1 & -y_A^{N} x_B^{N} &- y_A^{N} y_B^{N}
\end{pmatrix}, \\
\boldsymbol{y} = \begin{pmatrix}
x_A^{1} \\
y_A^{1} \\
\cdots \\
\cdots \\
x_A^{N} \\
y_A^{N} \\
\end{pmatrix},
\boldsymbol{h} = \begin{pmatrix}
h_{11} \\
h_{12} \\
h_{13} \\
h_{21} \\
h_{22} \\
h_{23} \\
h_{31} \\
h_{32}
\end{pmatrix}
.
以上の記号を用いてホモグラフィ変換は次のようになります。
\begin{align}
X \boldsymbol{h} = \boldsymbol{y}
\end{align}
これは解析的に解くことができて、次のようになります。
\begin{align}
& X \boldsymbol{h} = \boldsymbol{y} \\
\Leftrightarrow & X^\top X \boldsymbol{h} = X^\top \boldsymbol{y} \\
\Leftrightarrow & \boldsymbol{h} = (X^\top X)^{-1} X^\top
\end{align}
$X^\top X$は$8\times8$行列なので、逆行列も比較的求めやすいです。
また、RANSAC回帰等を用いて、外れ値を除きながら回帰分析をしても良いです。画像レジストレーションの場合は、画像のディスクリプタを用いて座標を対応させます。これらの対応点には誤った対応点も含まれるため、外れ値に強い手法を使ったほうが良いです。
実装
データはOpenCVのリポジトリのデータをダウンロードしました。
"""
Copyright 2023 ground0state
Released under the MIT license
"""
from random import Random
from os.path import join
import cv2
import numpy as np
from sklearn.linear_model import LinearRegression
def get_n_trials(n, p=0.99, w=0.2):
nume = np.log(1-p)
deno = np.log(1-np.power(w, n))
return int(nume / deno) + 1
def get_sift_matches(img1, img2, n_matches=None):
sift = cv2.SIFT_create()
keypoints1, descriptors1 = sift.detectAndCompute(img1, None)
keypoints2, descriptors2 = sift.detectAndCompute(img2, None)
bf = cv2.BFMatcher()
matches = bf.match(descriptors1, descriptors2)
matches = sorted(matches, key=lambda x: x.distance)[:n_matches]
return matches, keypoints1, keypoints2
def _create_homography_matrix(src_point, dst_point):
x1, y1 = dst_point
x2, y2 = src_point
matrix = [[x2, y2, 1, 0, 0, 0, -x1*x2, -x1*y2],
[0, 0, 0, x2, y2, 1, -y1*x2, -y1*y2]]
return matrix, [x1, y1]
def _estimate_homography(
src_points, dst_points, max_iters=5000, reprojection_threshold=5, seed=0, verbose=0):
random = Random(seed)
X, y = zip(*[_create_homography_matrix(src_pt, dst_pt)
for src_pt, dst_pt in zip(src_points, dst_points)])
X = np.concatenate(X).astype(np.float32)
y = np.concatenate(y).astype(np.float32)
n_points = len(src_points)
n_sample = 4
magic_number = 15
n_trials = min(max_iters, get_n_trials(n_sample, w=magic_number/n_points))
if verbose >= 1:
print(" n_trials:", n_trials)
n_inlier = 0
for i_iter in range(n_trials):
sampled_idx = random.sample(range(n_points), n_sample)
sampled_X, sampled_y = zip(
*[_create_homography_matrix(src_points[i], dst_points[i]) for i in sampled_idx])
sampled_X = np.concatenate(sampled_X)
sampled_y = np.concatenate(sampled_y)
estimator = LinearRegression(
fit_intercept=False).fit(sampled_X, sampled_y)
y_pred = estimator.predict(X)
residual_distance = np.linalg.norm(
y_pred.reshape(-1, 2) - y.reshape(-1, 2), axis=1)
temp_inlier_mask = residual_distance <= reprojection_threshold
temp_n_inlier = temp_inlier_mask.sum()
if temp_n_inlier > n_inlier:
n_inlier = temp_n_inlier
inlier_mask = temp_inlier_mask
if verbose >= 2:
print(f" trial_id: {i_iter}, n_inlier_points: {n_inlier}")
if verbose >= 1:
print(f" n_inlier_points: {n_inlier}")
selected_X, selected_y = zip(*[_create_homography_matrix(
src_points[i], dst_points[i]) for i in range(len(inlier_mask)) if inlier_mask[i]])
selected_X = np.concatenate(selected_X)
selected_y = np.concatenate(selected_y)
estimator = LinearRegression(
fit_intercept=False).fit(selected_X, selected_y)
h = estimator.coef_.copy()
h = np.concatenate([h, [1]])
h = h.reshape(3, 3)
return h
def find_homography(
matches,
src_keypoints,
dst_keypoints,
max_iters=5000,
reprojection_threshold=5,
seed=0,
verbose=0
):
src_points, dst_points = [], []
for i_match in matches:
src_points.append(src_keypoints[i_match.trainIdx].pt)
dst_points.append(dst_keypoints[i_match.queryIdx].pt)
h = _estimate_homography(
src_points,
dst_points,
max_iters=max_iters,
reprojection_threshold=reprojection_threshold,
seed=seed,
verbose=verbose)
return h
def apply_homography(src_image, dst_image, h):
result = cv2.warpPerspective(src_image, h, dst_image.shape[:2][::-1])
return result
def blend_image(src_image, dst_image):
result = cv2.addWeighted(dst_image, 0.5, src_image, 0.5, 0)
return result
def find_homography_opencv(
matches,
src_keypoints,
dst_keypoints,
max_iters=5000,
reprojection_threshold=5.0,
):
dst_kpt = np.float32(
[dst_keypoints[m.queryIdx].pt for m in matches]).reshape(-1, 1, 2)
src_kpt = np.float32(
[src_keypoints[m.trainIdx].pt for m in matches]).reshape(-1, 1, 2)
h, mask = cv2.findHomography(
src_kpt, dst_kpt, cv2.RANSAC,
ransacReprojThreshold=reprojection_threshold,
maxIters=max_iters)
return h
def main():
# Parameters ---------------
seed = 0
data_dir = "../data"
src_image_name = "leuvenB.jpg"
dst_image_name = "leuvenA.jpg"
n_matches = 100
max_iters = 5000
verbose = 2
# --------------------------
dst_image = cv2.imread(
join(data_dir, dst_image_name), cv2.IMREAD_COLOR)
src_image = cv2.imread(
join(data_dir, src_image_name), cv2.IMREAD_COLOR)
dst_image_gray = cv2.cvtColor(dst_image, cv2.COLOR_BGR2GRAY)
src_image_gray = cv2.cvtColor(src_image, cv2.COLOR_BGR2GRAY)
matches, dst_keypoints, src_keypoints = get_sift_matches(
dst_image_gray, src_image_gray, n_matches)
h = find_homography(matches, src_keypoints, dst_keypoints,
max_iters=max_iters, seed=seed, verbose=verbose)
# OpenCVを使う場合はこっち
# h = find_homography_opencv(
# matches, src_keypoints, dst_keypoints,
# max_iters=max_iters)
result = apply_homography(src_image, dst_image, h)
blended = blend_image(result, dst_image)
cv2.imwrite("warp.jpg", result)
cv2.imwrite("blend.jpg", blended)
if __name__ == "__main__":
main()
dst_image
src_image
blended
Reference
中村恭之, 小枝正直, 上田悦子,『OpenCVによるコンピュータビジョン・機械学習入門』, 講談社サイエンティフィク, 2017.