計算視覺——基礎矩陣和極點極線

• 基礎矩陣

在計算機視覺中，基礎矩陣（Fundamental matrix）F是一個3×3的矩陣，表達了立體像對的像點之間的對應關系。在對極幾何中，對於立體像對中的一對同名點，它們的齊次化圖像坐標分別為p與 p'，

表示一條必定經過p'的直線（極線）。這意味着立體像對的所有同名點對都滿足： 

 

 

F矩陣中蘊含了立體像對的兩幅圖像在拍攝時相互之間的空間幾何關系（外參數）以及相機檢校參數（內參數），包括旋轉、位移、像主點坐標和焦距。因為F矩陣的秩為2，並且可以自由縮放（尺度化），所以只需7對同名點即可估算出F的值。
基礎矩陣這一概念由Q. T. Luong在他那篇很有影響力的博士畢業論文中提出。

[1] Faugeras則是在1992年發表的著作

[2] 中以上面的關系式給出了F矩陣的定義。盡管Longuet-Higgins提出的本質矩陣也滿足類似的關系式，但本質矩陣中並不蘊含相機檢校參數。

本質矩陣與基礎矩陣之間的關系可由下式表達：

 

其中K和K'分別為兩個相機的內參數矩陣。

• 極點和極線
•  極點

每一個極點之處，增益衰減-3db,並移相-45度。極點之后每十倍頻，增益下降20db.零點與極點相反；每一個零點之處，增益增加3db,並移相45度。零點之后，每十倍頻，增益增加20db。

閉環增益A0：a/1+ab=1/b(當a很大時），其中a為開環增益，b為反饋因子，可以理解為反饋量和輸出量的比值，當開環增益趨近於無窮大時，閉環增益就是反饋因子的倒數。

環路增益：T=a*b

對運放來說：閉環增益（1/b)的傳遞函數的零點是環路增益(ab) 傳遞函數的極點；閉環增益的傳遞函數的極點是環路增益傳遞函數的零點；而我們在反饋的時候，是希望在相位下降到180度之前，環路增益大於一，所以我們需要消除一個環路增益函數的極點（即閉環增益零點），以免發生震盪。

• 極線

在數學中， 極線通常是一個適用於圓錐曲線的概念，如果圓錐曲線的切於A,B兩點的切線相交於P點，那么P點稱為直線AB關於該曲線的極點(pole)，直線AB稱為P點的極線(polar)。

• 極線的幾何定義

上面定義僅適用於P點在此圓錐曲線外部的情況。實際上，在P點在圓錐曲線內部的時候同樣可以定義極線，這時我們可以認為極線是過P點做此圓錐曲線兩條虛切線切點的連線.特別的，如果這個圓錐曲線是一個圓，我們同樣有圓的極線和極點的概念。

• 極線的幾何性質

對於圓錐曲線,兩個點的切線的交點的極線即這兩點的連線。 [2] 此外，過不在圓錐曲線上任意一點做兩條和此曲線相交的直線得出四個點，那么這四個點確定的四邊形的對角線交點在該點的極線。我們也可以把這個性質作為圓錐曲線的極線的定義。 [3]
而當一個動點移動到曲線上，那么它的極線就退化為過這點的切線， 所以，極點和極線的思想實際上是曲線上點和過該點切線的思想的一般化。

• 求解圖像之間的基礎矩陣

 1.歸一化8點算法

 基本矩陣是由下述方程定義

 其中x相對於x′是兩幅圖像的任意一對匹配點。由於每一組點的匹配提供了計算FF系數的一個線性方程，當給定至少77個點（3×33×3的齊次矩陣減去一個尺度，以及一個秩為22的約束），方程就可以計算出未知的FF。

我們記點的坐標為x=(x,y,1)^T，x'=(x',y',1)^T,x=(x,y,1)^ T，x′=(x′,y′,1) ^T 

則對應的方程為 

 

 

 

 展開后有

 

 把矩陣FFF寫成列向量的形式，則有：

 

 給定n組點的集合，我們有如下方程：

 

 如果存在確定（非零）解，則系數矩陣A的秩最多是8。由於F是齊次矩陣，所以如果矩陣A的秩為8，則在差一個尺度因子的情況下解是唯一的。可以直接用線性算法解得。

 最終解為

•  生成基礎矩陣

 

from PIL import Image
from numpy import *
from pylab import *
import numpy as np

import PCV.geometry.camera as camera
import PCV.geometry.homography as homography
import PCV.geometry.sfm as sfm
import PCV.localdescriptors.sift as sift

im1 = array(Image.open('002.jpg'))
sift.process_image('002.jpg', 'im1.sift')

im2 = array(Image.open('001.jpg'))
sift.process_image('001.jpg', 'im2.sift')

l1, d1 = sift.read_features_from_file('im1.sift')
l2, d2 = sift.read_features_from_file('im2.sift')

matches = sift.match_twosided(d1, d2)

ndx = matches.nonzero()[0]
x1 = homography.make_homog(l1[ndx, :2].T)
ndx2 = [int(matches[i]) for i in ndx]
x2 = homography.make_homog(l2[ndx2, :2].T)

d1n = d1[ndx]
d2n = d2[ndx2]
x1n = x1.copy()
x2n = x2.copy()

figure(figsize=(16,16))
sift.plot_matches(im1, im2, l1, l2, matches, True)
show()

def F_from_ransac(x1, x2, model, maxiter=5000, match_threshold=1e-6):
    """ Robust estimation of a fundamental matrix F from point
    correspondences using RANSAC (ransac.py from
    http://www.scipy.org/Cookbook/RANSAC).

    input: x1, x2 (3*n arrays) points in hom. coordinates. """

    import PCV.tools.ransac as ransac
    data = np.vstack((x1, x2))
    d = 10 # 20 is the original
    F, ransac_data = ransac.ransac(data.T, model,
                                   8, maxiter, match_threshold, d, return_all=True)
    return F, ransac_data['inliers']

model = sfm.RansacModel()
F, inliers = F_from_ransac(x1n, x2n, model, maxiter=5000, match_threshold=1e-5)
print F

P1 = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0]])
P2 = sfm.compute_P_from_fundamental(F)

X = sfm.triangulate(x1n[:, inliers], x2n[:, inliers], P1, P2)


cam1 = camera.Camera(P1)
cam2 = camera.Camera(P2)
x1p = cam1.project(X)
x2p = cam2.project(X)

figure(figsize=(16, 16))
imj = sift.appendimages(im1, im2)
imj = vstack((imj, imj))

imshow(imj)

cols1 = im1.shape[1]
rows1 = im1.shape[0]
for i in range(len(x1p[0])):
    if (0<= x1p[0][i]<cols1) and (0<= x2p[0][i]<cols1) and (0<=x1p[1][i]<rows1) and (0<=x2p[1][i]<rows1):
        plot([x1p[0][i], x2p[0][i]+cols1],[x1p[1][i], x2p[1][i]],'c')
axis('off')
show()

d1p = d1n[inliers]
d2p = d2n[inliers]

# Read features
im3 = array(Image.open('003.jpg'))
sift.process_image('003.jpg', 'im3.sift')
l3, d3 = sift.read_features_from_file('im3.sift')

matches13 = sift.match_twosided(d1p, d3)

ndx_13 = matches13.nonzero()[0]
x1_13 = homography.make_homog(x1p[:, ndx_13])
ndx2_13 = [int(matches13[i]) for i in ndx_13]
x3_13 = homography.make_homog(l3[ndx2_13, :2].T)

figure(figsize=(16, 16))
imj = sift.appendimages(im1, im3)
imj = vstack((imj, imj))

imshow(imj)

cols1 = im1.shape[1]
rows1 = im1.shape[0]
for i in range(len(x1_13[0])):
    if (0<= x1_13[0][i]<cols1) and (0<= x3_13[0][i]<cols1) and (0<=x1_13[1][i]<rows1) and (0<=x3_13[1][i]<rows1):
        plot([x1_13[0][i], x3_13[0][i]+cols1],[x1_13[1][i], x3_13[1][i]],'c')
axis('off')
show()

P3 = sfm.compute_P(x3_13, X[:, ndx_13])

print P1
print P2
print P3

基礎矩陣

 

 p1,p2,p3相機矩陣

 

 通過修改如圖所示參數，可以修改為7點法或10點法

• 極線代碼

 

import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt


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 = img1.shape
    img1 = cv.cvtColor(img1, cv.COLOR_GRAY2BGR)
    img2 = cv.cvtColor(img2, cv.COLOR_GRAY2BGR)
    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 = cv.line(img1, (x0, y0), (x1, y1), color, 1)
        img1 = cv.circle(img1, tuple(pt1), 5, color, -1)
        img2 = cv.circle(img2, tuple(pt2), 5, color, -1)
    return img1, img2


img1 = cv.imread('D:/Study/untitled1/333/5.jpg', 0)  # queryimage # left image
img2 = cv.imread('D:/Study/untitled1/333/6.jpg', 0)  # trainimage # right image
sift = cv.xfeatures2d.SIFT_create()
# find the keypoints and descriptors with SIFT
kp1, des1 = sift.detectAndCompute(img1, None)
kp2, des2 = sift.detectAndCompute(img2, None)
# FLANN parameters
FLANN_INDEX_KDTREE = 1
index_params = dict(algorithm=FLANN_INDEX_KDTREE, trees=5)
search_params = dict(checks=50)
flann = cv.FlannBasedMatcher(index_params, search_params)
matches = flann.knnMatch(des1, des2, k=2)
good = []
pts1 = []
pts2 = []
# ratio test as per Lowe's paper
for i, (m, n) in enumerate(matches):
    if m.distance < 0.8 * n.distance:
        good.append(m)
        pts2.append(kp2[m.trainIdx].pt)
        pts1.append(kp1[m.queryIdx].pt)

pts1 = np.int32(pts1)
pts2 = np.int32(pts2)
F, mask = cv.findFundamentalMat(pts1, pts2, cv.FM_LMEDS)
# We select only inlier points
pts1 = pts1[mask.ravel() == 1]
pts2 = pts2[mask.ravel() == 1]

# Find epilines corresponding to points in right image (second image) and
# drawing its lines on left image
lines1 = cv.computeCorrespondEpilines(pts2.reshape(-1, 1, 2), 2, F)
lines1 = lines1.reshape(-1, 3)
img5, img6 = drawlines(img1, img2, lines1, pts1, pts2)
# Find epilines corresponding to points in left image (first image) and
# drawing its lines on right image
lines2 = cv.computeCorrespondEpilines(pts1.reshape(-1, 1, 2), 1, F)
lines2 = lines2.reshape(-1, 3)
img3, img4 = drawlines(img2, img1, lines2, pts2, pts1)
plt.subplot(121), plt.imshow(img5)
plt.subplot(122), plt.imshow(img3)
plt.show()

 

•  實驗結果及分析

sift特征匹配：

 

 

 

 極點極線

 

 

 

 

 

 

 

•  實驗總結

 1. 在計算基礎矩陣時多次報錯或失敗，原因經了解主要是圖像壓縮之后會丟失很多匹配特征點，因此在計算基礎矩陣是一定要控制好圖像壓縮大小

 2. 由於兩幅圖像在匹配的時候有不少錯誤的匹配，所以計算的基本矩陣有較大的誤差。左右視圖中都可以發現很多的點在對極線附近但並沒有完全落在對極線上。我們可以觀察到左右視圖的對極線都響應地匯聚到一點，那點就是極點。這一對匹配中極點都落在圖像內，當然也有些情況對極線會落在圖像外

 3.  8點算法是計算基本矩陣的最簡單的方法。為了提高解的穩定性和精度，往往會對輸入點集的坐標先進行歸一化處理。

 4.  計算 基礎矩陣的前提是要先得到其sift特征匹配結果，通過sift特征匹配結果得到F基礎矩陣，因此要確保sift尋找出更多的特征匹配點。