OpenCVやってみる-22. エピポーラ幾何

今回は、OpenCV チュートリアルの以下のページを参考に進めてみたいと思います。

エピポーラ幾何の理論はチュートリアルサイトの説明を参照で。

今回のテーマ画像
特徴点検出、特徴量算出
マッチング
Fundamental matrix計算
Epipolar line描画
Fundamental matrixの計算について
- Fundamental matrixについて
ここまで
参考

今回のテーマ画像

今回は、春のパン祭り台紙の2視点分の画像でマッチング、Epipolar lineの描画をしてみたいと思います。
似たような特徴点が複数あるので難しそうと思いましたが、RANSAC等適用すればうまくマッチングできるのではと。

import cv2
img1 = cv2.imread('harupan_210402_1.jpg')
img2 = cv2.imread('harupan_210402_2.jpg')
img1= cv2.resize(img1, None, fx=0.2, fy=0.2, interpolation=cv2.INTER_AREA)
img2= cv2.resize(img2, None, fx=0.2, fy=0.2, interpolation=cv2.INTER_AREA)

img12 = np.hstack((img1, img2))
cv2.imshow('Image', img12)
cv2.waitKey(0)

画像サイズはいずれも806x605になっています。

>>> img1.shape
(806, 605, 3)
>>> img2.shape
(806, 605, 3)

特徴点検出、特徴量算出

SIFTを使いました。

sift = cv2.SIFT_create()
kp1, des1 = sift.detectAndCompute(img1, None)
kp2, des2 = sift.detectAndCompute(img2, None)

img1_kp = cv2.drawKeypoints(img1, kp1, None, flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
img2_kp = cv2.drawKeypoints(img2, kp2, None, flags=cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)

img12_kp = np.hstack((img1_kp, img2_kp))
cv2.imshow('Keypoints', img12_kp)
cv2.waitKey(0)

マッチング

ここでは総当たりマッチングを使います。

bf = cv2.BFMatcher()
matches = bf.match(des1, des2)
matches = sorted(matches, key = lambda x:x.distance)

img_matches = cv2.drawMatches(img1, kp1, img2, kp2, matches[:200], None, flags=cv2.DrawMatchesFlags_NOT_DRAW_SINGLE_POINTS)
cv2.imshow('Matches', img_matches)
cv2.waitKey(0)

いまいちな結果なので、ratio testも入れてみます。

matches_rt = bf.knnMatch(des1, des2, k=2)
good_matches = []
for m,n in matches_rt:
    if m.distance < 0.75*n.distance:
        good_matches.append(m)

good_matches = sorted(good_matches, key = lambda x:x.distance)
img_good_matches = cv2.drawMatches(img1, kp1, img2, kp2, good_matches[:200], None, flags=cv2.DrawMatchesFlags_NOT_DRAW_SINGLE_POINTS)
cv2.imshow('Good matches', img_good_matches)
cv2.waitKey(0)

これでもまだまだですが、この結果を使っていきます。

ちなみに、good_matchesには1169個のマッチが含まれていました。

>>> len(good_matches)
1169

Fundamental matrix計算

Fundamental matrixの計算をやってみます。
この過程でマッチング結果のinlier、outlierを判定してくれるので、それも見てみます。

pts1 = []
pts2 = []

for i,m in enumerate(good_matches):
    pts1.append(kp1[m.queryIdx].pt)
    pts2.append(kp2[m.trainIdx].pt)

import numpy as np
pts1 = np.int32(pts1)
pts2 = np.int32(pts2)
F, mask = cv2.findFundamentalMat(pts1, pts2, cv2.FM_LMEDS)

pts1_inliers = pts1[mask.ravel() == 1]
pts2_inliers = pts2[mask.ravel() == 1]

good_matches_inliers = [m for i,m in enumerate(good_matches[:200]) if mask[i,0] == 1]

img_good_matches_inliers = cv2.drawMatches(img1,kp1,img2,kp2,good_matches_inliers[:200],None,flags=cv2.DrawMatchesFlags_NOT_DRAW_SINGLE_POINTS)
cv2.imshow('Good matches(inliers)', img_good_matches_inliers)
cv2.waitKey(0)

しっかり合っているものだけ残っているようです。すごい！

ちなみに、マッチング点が1169点あったのに対して、inlierは866点になっていました。

>>> len(pts1_inliers)
866

Epipolar line描画

では計算したFundamental matrixを元にEpipolar lineを描画するとどうなるかというと、

def drawlines(img1,img2,lines,pts1,pts2):
    ''' img1 - image on which we draw the epilines for the points in img2
        lines - corresponding epilines '''
    r,c,ch = img1.shape
    for r,pt1,pt2 in zip(lines,pts1,pts2):
        color = tuple(np.random.randint(0,255,3).tolist())
        x0,y0 = map(int, [0, -r[2]/r[1] ])
        x1,y1 = map(int, [c, -(r[2]+r[0]*c)/r[1] ])
        img1 = cv2.line(img1, (x0,y0), (x1,y1), color,1)
        img1 = cv2.circle(img1,tuple(pt1),5,color,-1)
        img2 = cv2.circle(img2,tuple(pt2),5,color,-1)
    return img1,img2

pts1_inliers_200 = pts1_inliers[0:200,:]
pts2_inliers_200 = pts2_inliers[0:200,:]

lines1 = cv2.computeCorrespondEpilines(pts2_inliers_200.reshape(-1,1,2),2,F)
lines1 = lines1.reshape(-1,3)
img5,img6 = drawlines(img1.copy(),img2.copy(),lines1,pts1,pts2)
lines2 = cv2.computeCorrespondEpilines(pts1_inliers_200.reshape(-1,1,2),1,F)
lines2 = lines2.reshape(-1,3)
img3,img4 = drawlines(img2.copy(),img1.copy(),lines2,pts2,pts1)
img35 = np.hstack((img3,img5))
cv2.imshow('Epipolar lines', img35)
cv2.waitKey(0)

こんな感じになっちゃいました。思っていたのと違う…

試しにFundamental matrixの計算で使うマッチング点数を減らしてみたりしましたが、

F, mask = cv2.findFundamentalMat(pts1[0:200,:], pts2[0:200,:], cv2.FM_LMEDS)

pts1_inliers = pts1[mask.ravel() == 1]
pts2_inliers = pts2[mask.ravel() == 1]

pts1_inliers_200 = pts1_inliers[0:200,:]
pts2_inliers_200 = pts2_inliers[0:200,:]

lines1 = cv2.computeCorrespondEpilines(pts2_inliers_200.reshape(-1,1,2),2,F)
lines1 = lines1.reshape(-1,3)
img5,img6 = drawlines(img1.copy(),img2.copy(),lines1,pts1,pts2)
lines2 = cv2.computeCorrespondEpilines(pts1_inliers_200.reshape(-1,1,2),1,F)
lines2 = lines2.reshape(-1,3)
img3,img4 = drawlines(img2.copy(),img1.copy(),lines2,pts2,pts1)
img35 = np.hstack((img3,img5))
cv2.imshow('Epipolar lines', img35)
cv2.waitKey(0)

これまた変な結果に。

チュートリアルの最後のほうをよく見てみると、特徴点が1つの平面上にあると正しくFundamental matrixを計算できないという記載がありました。

Fundamental matrixの計算について

考察してみます。
以下のサイトのpdfを参考にしています。ケンブリッジ大学出版の教科書です。

Computer Vision Models

Amazonで購入も可能なようです。

Amazon | Computer Vision: Models, Learning, and Inference | Prince, Simon J. D. | Machine Vision

Fundamental matrixについて

Fundamental matrixがどういうものかというと、途中経過は端折りますが、

任意の実空間座標 $\boldsymbol{w}$ と、2視点からの画像上でこれに対応する点 $\boldsymbol{x_1}, \boldsymbol{x_2}$ (同次座標で)について、以下の式が成り立つ

$\boldsymbol{x_1} \boldsymbol{F} \boldsymbol{x_2} = 0 \tag{1}$

というものです。

$\boldsymbol{x_1} = \lbrack x_1 y_1 1 \rbrack ^T , \boldsymbol{x_2} = \lbrack x_2 y_2 1 \rbrack ^T , \boldsymbol{F} = \begin{bmatrix} f_{11} & f_{12} & f_{13} \\ f_{21} & f_{22} & f_{23} \\ f_{31} & f_{32} & f_{33} \end{bmatrix}$

というように要素を書き下すと、式(1)は

$x_2 x_1 f_{11} + x_2 y_1 f_{12} + x_2 f_{13} + y_2 x_1 f_{21} + y_2 y_1 f_{f22} + y_2 f_{23} + x_1 f_{31} + y_1 f_{32} + f_{33} = 0 \tag{2}$ $\lbrack x_2 x_1 \: x_2 y_1 \: x_2 \: y_2 x_1 \: y_2 y_1 \: y_2 \: x_1 \: y_1 \: 1 \rbrack \boldsymbol{f} = 0 \tag{3}$ ( $\boldsymbol{f} = \lbrack f_{11} \: f_{12} \: f_{13} \: f_{21} \: f_{22} \: f_{23} \: f_{31} \: f_{32} \: f_{33} \rbrack ^T$ )

と書けます。
$\boldsymbol{F}$ はスケールの任意性があるので、自由度が8となっていて、8個以上のマッチングペアがあれば一意に決まります。

ここで、マッチングペアの点に対する実空間座標が全て同じ平面上にあった場合を考えます。

実空間座標を $\boldsymbol{w}=\lbrack u \; v \; w \rbrack ^T$ とすると、画像上の点(同次座標) $\boldsymbol{x_1}, \boldsymbol{x_2}$ は $u, v, w$ を使って以下のように表されます。

$\lambda_1 \boldsymbol{x_1} = \boldsymbol{\Omega_1} \boldsymbol{w} + \boldsymbol{\tau_1}$
$\lambda_2 \boldsymbol{x_2} = \boldsymbol{\Omega_2} \boldsymbol{w} + \boldsymbol{\tau_2}$

要素ごとに書くと、

$\lambda_1 x_1 = \Omega_{1_{11}} u + \Omega_{1_{12}} v + \Omega_{1_{13}} w + \tau_{1_1}$
$\lambda_1 y_1 = \Omega_{1_{21}} u + \Omega_{1_{22}} v + \Omega_{1_{23}} w + \tau_{1_2}$
$\lambda_1 = \Omega_{1_{31}} u + \Omega_{1_{32}} v + \Omega_{1_{33}} w + \tau_{1_3}$
$\lambda_2 x_2 = \Omega_{2_{11}} u + \Omega_{2_{12}} v + \Omega_{2_{13}} w + \tau_{2_1}$
$\lambda_2 y_2 = \Omega_{2_{21}} u + \Omega_{2_{22}} v + \Omega_{2_{23}} w + \tau_{2_2}$
$\lambda_2 = \Omega_{2_{31}} u + \Omega_{2_{32}} v + \Omega_{2_{33}} w + \tau_{2_3}$

式(2)に $\lambda_1 \lambda_2$ を掛けると、全ての項が $\lambda_1 x_1, \lambda_2 x_2, \lambda_1, \lambda_2$ の積で表されるので、上の6式を代入してみると、式(2)における係数がいずれも $u, v, w, uv, uw, vw, u^2, v^2, w^2, 1$ の線型結合で表されることがわかります。なので、式(3)について8個の線形独立な係数ベクトルを用意することができ、 $\boldsymbol{f}$ を一意に(定数倍の自由度を除いて)決めることができます。