[OpenCV-Python] OpenCV 中攝像機標定和 3D 重構 部分 VII

部分 VII
攝像機標定和 3D 重構

OpenCV-Python 中文教程（搬運）目錄

42 攝像機標定

目標
　　• 學習攝像機畸變以及攝像機的內部參數和外部參數
　　• 學習找到這些參數，對畸變圖像進行修復

42.1 基礎
　　今天的低價單孔攝像機（照相機）會給圖像帶來很多畸變。畸變主要有兩種：徑向畸變和切想畸變。如下圖所示，用紅色直線將棋盤的兩個邊標注出來，但是你會發現棋盤的邊界並不和紅線重合。所有我們認為應該是直線的也都凸出來了。你可以通過訪問Distortion (optics)獲得更多相關細節。

　　　　
這種畸變可以通過下面的方程組進行糾正：
　　　　　　

於此相似，另外一個畸變是切向畸變，這是由於透鏡與成像平面不可能絕對平行造成的。這種畸變會造成圖像中的某些點看上去的位置會比我們認為的位置要近一些。它可以通過下列方程組進行校正：
　　　　　　
簡單來說，如果我們想對畸變的圖像進行校正就必須找到五個造成畸變的系數：
　　　　　　
除此之外，我們還需要再找到一些信息，比如攝像機的內部和外部參數。
內部參數是攝像機特異的。它包括的信息有焦距（f x ,f y ），光學中心（c x ,c y ）
等。這也被稱為攝像機矩陣。它完全取決於攝像機自身，只需要計算一次，以后就可以已知使用了。可以用下面的 3x3 的矩陣表示：
　　　　　　
外部參數與旋轉和變換向量相對應，它可以將 3D 點的坐標轉換到坐標系統中。
在 3D 相關應用中，必須要先校正這些畸變。為了找到這些參數，我們必須要提供一些包含明顯圖案模式的樣本圖片（比如說棋盤）。我們可以在上面找到一些特殊點（如棋盤的四個角點）。我們起到這些特殊點在圖片中的位置以及它們的真是位置。有了這些信息，我們就可以使用數學方法求解畸變系數。這就是整個故事的摘要了。為了得到更好的結果，我們至少需要 10 個這樣的圖案模式。
42.2 代碼
如上所述，我們至少需要 10 圖案模式來進行攝像機標定。OpenCV 自帶了一些棋盤圖像（/sample/cpp/left001.jpg--left14.jpg）, 所以我們可以使用它們。為了便於理解，我們可以認為僅有一張棋盤圖像。重要的是在進行攝像機標定時我們要輸入一組 3D 真實世界中的點以及與它們對應 2D 圖像中的點。2D 圖像的點可以在圖像中很容易的找到。（這些點在圖像中的位置是棋盤上兩個黑色方塊相互接觸的地方）

那么真實世界中的 3D 的點呢？這些圖像來源與靜態攝像機和棋盤不同的擺放位置和朝向。所以我們需要知道（X，Y，Z）的值。但是為了簡單，我們可以說棋盤在 XY 平面是靜止的，（所以 Z 總是等於 0）攝像機在圍着棋盤移動。這種假設讓我們只需要知道 X，Y 的值就可以了。現在為了求 X，Y 的值，我們只需要傳入這些點（0,0），（1,0），（2,0）...，它們代表了點的位置。在這個例子中，我們的結果的單位就是棋盤（單個）方塊的大小。但是如果我們知道單個方塊的大小（加入說 30mm），我們輸入的值就可以是（0,0），（30,0），（60,0）...，結果的單位就是 mm。（在本例中我們不知道方塊的大小，因為不是我們拍的，所以只能用前一種方法了）。
3D 點被稱為 對象點，2D 圖像點被稱為 圖像點。

42.2.1 設置
　　為了找到棋盤的圖案，我們要使用函數 cv2.findChessboardCorners()。我們還需要傳入圖案的類型，比如說 8x8 的格子或 5x5 的格子等。在本例中我們使用的恨死 7x8 的格子。（通常情況下棋盤都是 8x8 或者 7x7）。它會返回角點，如果得到圖像的話返回值類型（Retval）就會是 True。這些角點會按順序排列（從左到右，從上到下）。
其他：這個函數可能不會找出所有圖像中應有的圖案。所以一個好的方法是編寫代碼，啟動攝像機並在每一幀中檢查是否有應有的圖案。在我們獲得圖案之后我們要找到角點並把它們保存成一個列表。在讀取下一幀圖像之前要設置一定的間隔，這樣我們就有足夠的時間調整棋盤的方向。繼續這個過程直到我們得到足夠多好的圖案。就算是我們舉得這個例子，在所有的 14 幅圖像中也不知道有幾幅是好的。所以我們要讀取每一張圖像從其中找到好的能用的。
其他：除了使用棋盤之外，我們還可以使用環形格子，但是要使用函數cv2.findCirclesGrid() 來找圖案。據說使用環形格子只需要很少的圖像就可以了。
在找到這些角點之后我們可以使用函數 cv2.cornerSubPix() 增加准確度。我們使用函數 cv2.drawChessboardCorners() 繪制圖案。所有的這些步驟都被包含在下面的代碼中了：

import numpy as np
import cv2
import glob

# termination criteria
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)

# prepare object points, like (0,0,0), (1,0,0), (2,0,0) ....,(6,5,0)
objp = np.zeros((6*7,3), np.float32)
objp[:,:2] = np.mgrid[0:7,0:6].T.reshape(-1,2)

# Arrays to store object points and image points from all the images.
objpoints = [] # 3d point in real world space
imgpoints = [] # 2d points in image plane.

images = glob.glob('*.jpg')

for fname in images:
    img = cv2.imread(fname)
    gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

    # Find the chess board corners
    ret, corners = cv2.findChessboardCorners(gray, (7,6),None)

    # If found, add object points, image points (after refining them)
    if ret == True:
        objpoints.append(objp)

        corners2 = cv2.cornerSubPix(gray,corners,(11,11),(-1,-1),criteria)
        imgpoints.append(corners2)

        # Draw and display the corners
        img = cv2.drawChessboardCorners(img, (7,6), corners2,ret)
        cv2.imshow('img',img)
        cv2.waitKey(500)

cv2.destroyAllWindows()

一副圖像和被繪制在上邊的圖案：
　　　　

42.2.2 標定
　　在得到了這些對象點和圖像點之后，我們已經准備好來做攝像機標定了。
我們要使用的函數是 cv2.calibrateCamera()。它會返回攝像機矩陣，畸變系數，旋轉和變換向量等。

ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(objpoints, imgpoints, gray.shape[::-1],None,None)

 

42.2.3 畸變校正
　　現在我們找到我們想要的東西了，我們可以找到一幅圖像來對他進行校正了。OpenCV 提供了兩種方法，我們都學習一下。不過在那之前我們可以使用從函數 cv2.getOptimalNewCameraMatrix() 得到的自由縮放系數對攝像機矩陣進行優化。如果縮放系數 alpha = 0，返回的非畸變圖像會帶有最少量的不想要的像素。它甚至有可能在圖像角點去除一些像素。如果 alpha = 1，所有的像素都會被返回，還有一些黑圖像。它還會返回一個 ROI 圖像，我們可以用來對結果進行裁剪。
我們讀取一個新的圖像（left2.ipg）

img = cv2.imread('left12.jpg')
h,  w = img.shape[:2]
newcameramtx, roi=cv2.getOptimalNewCameraMatrix(mtx,dist,(w,h),1,(w,h))

使用 cv2.undistort() 這是最簡單的方法。只需使用這個函數和上邊得到的 ROI 對結果進行裁剪。

# undistort
dst = cv2.undistort(img, mtx, dist, None, newcameramtx)

# crop the image
x,y,w,h = roi
dst = dst[y:y+h, x:x+w]
cv2.imwrite('calibresult.png',dst)

使用 remapping 這應該屬於“曲線救國”了。首先我們要找到從畸變圖像到非畸變圖像的映射方程。再使用重映射方程。

# undistort
mapx,mapy = cv2.initUndistortRectifyMap(mtx,dist,None,newcameramtx,(w,h),5)
dst = cv2.remap(img,mapx,mapy,cv2.INTER_LINEAR)

# crop the image
x,y,w,h = roi
dst = dst[y:y+h, x:x+w]
cv2.imwrite('calibresult.png',dst)

這兩中方法給出的結果是相同的。結果如下所示：

　　　　

你會發現結果圖像中所有的邊界都變直了。
現在我們可以使用 Numpy 提供寫函數（np.savez，np.savetxt 等）
將攝像機矩陣和畸變系數保存以便以后使用。

42.3 反向投影誤差
　　我們可以利用反向投影誤差對我們找到的參數的准確性進行估計。得到的結果越接近 0 越好。有了內部參數，畸變參數和旋轉變換矩陣，我們就可以使用 cv2.projectPoints() 將對象點轉換到圖像點。然后就可以計算變換得到圖像與角點檢測算法的絕對差了。然后我們計算所有標定圖像的誤差平均值。

mean_error = 0
for i in xrange(len(objpoints)):
    imgpoints2, _ = cv2.projectPoints(objpoints[i], rvecs[i], tvecs[i], mtx, dist)
    error = cv2.norm(imgpoints[i],imgpoints2, cv2.NORM_L2)/len(imgpoints2)
    tot_error += error

print "total error: ", mean_error/len(objpoints)

 

43 姿勢估計

目標
　　• 本節我們要學習使用 calib3D 模塊在圖像中創建 3D 效果

43.1 基礎
　　在上一節的攝像機標定中，我們已經得到了攝像機矩陣，畸變系數等。有了這些信息我們就可以估計圖像中圖案的姿勢，比如目標對象是如何擺放，如何旋轉等。對一個平面對象來說，我們可以假設 Z=0，這樣問題就轉化成攝像機在空間中是如何擺放（然后拍攝）的。所以，如果我們知道對象在空間中的姿勢，我們就可以在圖像中繪制一些 2D 的線條來產生 3D 的效果。我們來看一下怎么做吧。
我們的問題是，在棋盤的第一個角點繪制 3D 坐標軸（X，Y，Z 軸）。X軸為藍色，Y 軸為綠色，Z 軸為紅色。在視覺效果上來看，Z 軸應該是垂直與棋盤平面的。
首先我們要加載前面結果中攝像機矩陣和畸變系數。

import cv2
import numpy as np
import glob

# Load previously saved data
with np.load('B.npz') as X:
    mtx, dist, _, _ = [X[i] for i in ('mtx','dist','rvecs','tvecs')]

現在我們來創建一個函數：draw，它的參數有棋盤上的角點（使用cv2.findChessboardCorners() 得到）和要繪制的 3D 坐標軸上的點。

def draw(img, corners, imgpts):
    corner = tuple(corners[0].ravel())
    img = cv2.line(img, corner, tuple(imgpts[0].ravel()), (255,0,0), 5)
    img = cv2.line(img, corner, tuple(imgpts[1].ravel()), (0,255,0), 5)
    img = cv2.line(img, corner, tuple(imgpts[2].ravel()), (0,0,255), 5)
    return img

和前面一樣，我們要設置終止條件，對象點（棋盤上的 3D 角點）和坐標軸點。3D 空間中的坐標軸點是為了繪制坐標軸。我們繪制的坐標軸的長度為3。所以 X 軸從（0,0,0）繪制到（3,0,0），Y 軸也是。Z 軸從（0,0,0）繪制到（0,0,-3）。負值表示它是朝着（垂直於）攝像機方向。

criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)
objp = np.zeros((6*7,3), np.float32)
objp[:,:2] = np.mgrid[0:7,0:6].T.reshape(-1,2)

axis = np.float32([[3,0,0], [0,3,0], [0,0,-3]]).reshape(-1,3)

很通常一樣我們需要加載圖像。搜尋 7x6 的格子，如果發現，我們就把它優化到亞像素級。然后使用函數:cv2.solvePnPRansac() 來計算旋轉和變換。但我們有了變換矩陣之后，我們就可以利用它們將這些坐標軸點映射到圖像平面中去。簡單來說，我們在圖像平面上找到了與 3D 空間中的點（3,0,0）,(0,3,0),(0,0,3) 相對應的點。然后我們就可以使用我們的函數 draw() 從圖像上的第一個角點開始繪制連接這些點的直線了。搞定！！！

for fname in glob.glob('left*.jpg'):
    img = cv2.imread(fname)
    gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
    ret, corners = cv2.findChessboardCorners(gray, (7,6),None)

    if ret == True:
        corners2 = cv2.cornerSubPix(gray,corners,(11,11),(-1,-1),criteria)

        # Find the rotation and translation vectors.
        rvecs, tvecs, inliers = cv2.solvePnPRansac(objp, corners2, mtx, dist)

        # project 3D points to image plane
        imgpts, jac = cv2.projectPoints(axis, rvecs, tvecs, mtx, dist)

        img = draw(img,corners2,imgpts)
        cv2.imshow('img',img)
        k = cv2.waitKey(0) & 0xff
        if k == 's':
            cv2.imwrite(fname[:6]+'.png', img)

cv2.destroyAllWindows()

結果如下，看到了嗎，每條坐標軸的長度都是 3 個格子的長度。

　　　　

43.1.1 渲染一個立方體
　　如果你想繪制一個立方體的話要對 draw() 函數進行如下修改：
修改后的 draw() 函數：

def draw(img, corners, imgpts):
    imgpts = np.int32(imgpts).reshape(-1,2)

    # draw ground floor in green
    img = cv2.drawContours(img, [imgpts[:4]],-1,(0,255,0),-3)

    # draw pillars in blue color
    for i,j in zip(range(4),range(4,8)):
        img = cv2.line(img, tuple(imgpts[i]), tuple(imgpts[j]),(255),3)

    # draw top layer in red color
    img = cv2.drawContours(img, [imgpts[4:]],-1,(0,0,255),3)

    return img

修改后的坐標軸點。它們是 3D 空間中的一個立方體的 8 個角點：

axis = np.float32([[0,0,0], [0,3,0], [3,3,0], [3,0,0],
                   [0,0,-3],[0,3,-3],[3,3,-3],[3,0,-3] ])

結果如下：

　　　　
如果你對計算機圖形學感興趣的話，為了增加圖像的真實性，你可以使用OpenGL 來渲染更復雜的圖形。（下一個目標）

44 對極幾何（Epipolar Geometry ）
目標
　　• 本節我們要學習多視角幾何基礎
　　• 學習什么是極點，極線，對極約束等

44.1 基本概念
　　在我們使用針孔相機時，我們會丟失大量重要的信心，比如說圖像的深度，或者說圖像上的點和攝像機的距離，因這是一個從 3D 到 2D 的轉換。因此一個重要的問題就產生了，使用這樣的攝像機我們能否計算除深度信息呢？答案就是使用多個相機。我們的眼睛就是這樣工作的，使用兩個攝像機（兩個眼睛），這被稱為立體視覺。我們來看看 OpenCV 在這方面給我們都提供了什么吧。
（《學習 OpenCV》一書有大量相關知識。）
在進入深度圖像之前，我們要先掌握一些多視角幾何的基本概念。在本節中我們要處理對極幾何。下圖為使用兩台攝像機同時對一個一個場景進行拍攝的示意圖。

　　　　　　
如果只是用一台攝像機我們不可能知道 3D 空間中的 X 點到圖像平面的距離，因為 OX 連線上的每個點投影到圖像平面上的點都是相同的。但是如果我們也考慮上右側圖像的話，直線 OX 上的點將投影到右側圖像上的不同位置。
所以根據這兩幅圖像，我們就可以使用三角測量計算出 3D 空間中的點到攝像機的距離（深度）。這就是整個思路。
直線 OX 上的不同點投射到右側圖像上形成的線 l
′
被稱為與 x 點對應的極
線
極
線。也就是說，我們可以在右側圖像中沿着這條極線找到 x 點。它可能在這條直線上某個位置（這意味着對兩幅圖像間匹配特征的二維搜索就轉變成了沿着極線的一維搜索。這不僅節省了大量的計算，還允許我們排除許多導致虛假匹配的點）。這被稱為 對極約束。與此相同，所有的點在其他圖像中都有與之對應的極線。平面 XOO' 被稱為 對極平面。
O 和 O' 是攝像機的中心。從上面的示意圖可以看出，右側攝像機的中心O' 投影到左側圖像平面的 e 點，這個點就被稱為 極點。極點就是攝像機中心連線與圖像平面的交點。因此點 e' 是左側攝像機的極點。有些情況下，我們可能不會在圖像中找到極點，它們可能落在了圖像之外（這說明這兩個攝像機不能拍攝到彼此）。
所有的極線都要經過極點。所以為了找到極點的位置，我們可以先找到多條極線，這些極線的交點就是極點。
本節我們的重點就是找到極線和極點。為了找到它們，我們還需要兩個元素， 本征矩陣（E ）和 基礎矩陣（F ）。本征矩陣包含了物理空間中兩個攝像機相關的旋轉和平移信息。如下圖所示（本圖來源自：學習 OpenCV）

　　　　　　
基礎矩陣 F 除了包含 E 的信息外還包含了兩個攝像機的內參數。由於 F包含了這些內參數，因此它可以它在像素坐標系將兩台攝像機關聯起來。（如果使用是校正之后的圖像並通過除以焦距進行了歸一化，F=E）。簡單來說，基礎矩陣 F 將一副圖像中的點映射到另一幅圖像中的線（極線）上。這是通過匹配兩幅圖像上的點來實現的。要計算基礎矩陣至少需要 8 個點（使用 8 點算法）。點越多越好，可以使用 RANSAC 算法得到更加穩定的結果。
44.2 代碼
為了得到基礎矩陣我們應該在兩幅圖像中找到盡量多的匹配點。我們可以使用 SIFT 描述符，FLANN 匹配器和比值檢測。

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

img1 = cv2.imread('myleft.jpg',0)  #queryimage # left image
img2 = cv2.imread('myright.jpg',0) #trainimage # right image

sift = cv2.SIFT()

# find the keypoints and descriptors with SIFT
kp1, des1 = sift.detectAndCompute(img1,None)
kp2, des2 = sift.detectAndCompute(img2,None)

# FLANN parameters
FLANN_INDEX_KDTREE = 0
index_params = dict(algorithm = FLANN_INDEX_KDTREE, trees = 5)
search_params = dict(checks=50)

flann = cv2.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 = cv2.findFundamentalMat(pts1,pts2,cv2.FM_LMEDS)

# We select only inlier points
pts1 = pts1[mask.ravel()==1]
pts2 = pts2[mask.ravel()==1]

下一步我們要找到極線。我們會得到一個包含很多線的數組。所以我們要定義一個新的函數將這些線繪制到圖像中。

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 = cv2.cvtColor(img1,cv2.COLOR_GRAY2BGR)
    img2 = cv2.cvtColor(img2,cv2.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 = 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

現在我們兩幅圖像中計算並繪制極線。

# Find epilines corresponding to points in right image (second image) and
# drawing its lines on left image
lines1 = cv2.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 = cv2.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()

下面是我得到的結果：
　　　　
從上圖可以看出所有的極線都匯聚以圖像外的一點，這個點就是極點。為了得到更好的結果，我們應該使用分辨率比較高的圖像和 non-planar點。

45 立體圖像中的深度地圖

目標
　　• 本節我們要學習為立體圖像制作深度地圖

45.1 基礎
　　在上一節中我們學習了對極約束的基本概念和相關術語。如果同一場景有兩幅圖像的話我們在直覺上就可以獲得圖像的深度信息。下面是的這幅圖和其中的數學公式證明我們的直覺是對的。（圖像來源 image courtesy）

　　　　　　
The above diagram contains equivalent triangles. Writing their
equivalent equations will yield us following result:
　　　　　　
x 和 x' 分別是圖像中的點到 3D 空間中的點和到攝像機中心的距離。B 是這兩個攝像機之間的距離，f 是攝像機的焦距。上邊的等式告訴我們點的深度與x 和 x' 的差成反比。所以根據這個等式我們就可以得到圖像中所有點的深度圖。
這樣就可以找到兩幅圖像中的匹配點了。前面我們已經知道了對極約束可以使這個操作更快更准。一旦找到了匹配，就可以計算出 disparity 了。讓我們看看在 OpenCV 中怎樣做吧。

45.2 代碼
　　下面的代碼顯示了構建深度圖的簡單過程。

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

imgL = cv2.imread('tsukuba_l.png',0)
imgR = cv2.imread('tsukuba_r.png',0)

stereo = cv2.createStereoBM(numDisparities=16, blockSize=15)
disparity = stereo.compute(imgL,imgR)
plt.imshow(disparity,'gray')
plt.show()

下圖左側為原始圖像，右側為深度圖像。如圖所示，結果中有很大的噪音。

　　　　　　
通過調整 numDisparities 和 blockSize 的值，我們會得到更好的結果。