OpenCV-Python] OpenCV 中攝像機標定和 3D 重構 部分 VII
https://www.cnblogs.com/Undo-self-blog/p/8448500.html
42 攝像機標定
目標
• 學習攝像機畸變以及攝像機的內部引數和外部引數
• 學習找到這些引數，對畸變影象進行修復
42.1 基礎
今天的低價單孔攝像機（照相機）會給影象帶來很多畸變。畸變主要有兩種：徑向畸變和切想畸變。如下圖所示，用紅色直線將棋盤的兩個邊標註出來，但是你會發現棋盤的邊界並不和紅線重合。所有我們認為應該是直線的也都凸出來了。你可以通過訪問Distortion (optics)獲得更多相關細節。
這種畸變可以通過下面的方程組進行糾正：
於此相似，另外一個畸變是切向畸變，這是由於透鏡與成像平面不可能絕對平行造成的。這種畸變會造成影象中的某些點看上去的位置會比我們認為的位置要近一些。它可以通過下列方程組進行校正：
![x_{corrected} = x + [ 2p_1xy + p_2(r^2+2x^2)] \\
y_{corrected} = y + [ p_1(r^2+ 2y^2)+ 2p_2xy]](http://opencv-python-tutroals.readthedocs.io/en/latest/_images/math/7817dd9b0769c056c23aeefc392cf06907a990a9.png)
簡單來說，如果我們想對畸變的影象進行校正就必須找到五個造成畸變的係數：
除此之外，我們還需要再找到一些資訊，比如攝像機的內部和外部引數。
內部引數是攝像機特異的。它包括的資訊有焦距（f x ,f y ），光學中心（c x ,c y ）
等。這也被稱為攝像機矩陣。它完全取決於攝像機自身，只需要計算一次，以後就可以已知使用了。可以用下面的 3x3 的矩陣表示：
![camera \; matrix = \left [ \begin{matrix} f_x & 0 & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{matrix} \right ]](http://opencv-python-tutroals.readthedocs.io/en/latest/_images/math/5e9d2db8b3561f68fa0688992fa45c13766a382c.png)
外部引數與旋轉和變換向量相對應，它可以將 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 的值，我們會得到更好的結果。