[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 的值，咱们会获得更好的结果。