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

時間 2019-11-10

標籤 opencv python 攝像機標定 3d 部分 vii 欄目 Python 简体版

原文原文鏈接

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

42 攝像機標定

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

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

　　　　
這種畸變能夠經過下面的方程組進行糾正：
　　　　　　 $x_{corrected} = x( 1 + k_1 r^2 + k_2 r^4 + k_3 r^6) \\ y_{corrected} = y( 1 + k_1 r^2 + k_2 r^4 + k_3 r^6)$

於此類似，另一個畸變是切向畸變，這是因爲透鏡與成像平面不可能絕對平行形成的。這種畸變會形成圖像中的某些點看上去的位置會比咱們認爲的位置要近一些。它能夠經過下列方程組進行校訂：
　　　　　　 $x_{corrected} = x + [ 2p_1xy + p_2(r^2+2x^2)] \\ y_{corrected} = y + [ p_1(r^2+ 2y^2)+ 2p_2xy]$
簡單來講，若是咱們想對畸變的圖像進行校訂就必須找到五個形成畸變的係數：
　　　　　　 $Distortion \; coefficients=(k_1 \hspace{10pt} k_2 \hspace{10pt} p_1 \hspace{10pt} p_2 \hspace{10pt} k_3)$
除此以外，咱們還須要再找到一些信息，好比攝像機的內部和外部參數。
內部參數是攝像機特異的。它包括的信息有焦距（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 ]$
外部參數與旋轉和變換向量相對應，它能夠將 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 圖像點被稱爲圖像點。app

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

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()。它會返回攝像機矩陣，畸變係數，旋轉和變換向量等。post

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 這應該屬於「曲線救國」了。首先咱們要找到從畸變圖像到非畸變圖像的映射方程。再使用重映射方程。ui

# 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:
　　　　　　 $disparity = x - x' = \frac{Bf}{Z}$
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()