【笔记】摄像机标定和3D重构
目录
一、摄像机标定
1.设置
2.标定
3.畸变矫正
4.反向投影差
三、对极几何
四、深度地图
一、摄像机标定
单孔摄像机（照相机）会给图像带来很多畸变，畸变主要有 “径向畸变” 和 “切向畸变”。如下图所示，用红色直线将棋盘的两个边标注出来，棋盘的边界与红线并不重合，我们认为应该是直线的棋盘边界也凸出来了。
所以，我们需要一些图案来进行摄像机标定。OpenCV官方提供了一些棋盘图像，在理论学习时我们可以使用它们，在实践过程时需要根据具体摄像头拍摄图片的内参和外参进行调节，本文使用图像均为OpenCV官方图像。
Ps：这些图像需要在OpenCV官网单独下载，资源包中还会有cpp，java接口的教程等很多内容。本渣渣学习时看跟多教程都说在源文件的“sample/data”目录下，我还以为是安装opencv时会自动下载好的，浪费我很多时间在自己的安装路径的文件夹下找这些图。。。
1.设置
为了找到棋盘的图案，我们要使用函数 cv2.findChessboardCorners()。我们还需要传入图案的类型，比如说 8x8 的格子或 5x5 的格子等。在本例中我们使用的是 7x8 的格子。（通常情况下棋盘都是 8x8 或者 7x7）。它会返回角点，如果得到图像的话返回值类型（ Retval）就会是 True。这些角点会按顺序排列（从左到右，从上到下）。除了使用棋盘之外，我们还可以使用环形格子，但是要使用函数cv2.findCirclesGrid() 来找图案。据说使用环形格子只需要很少的图像就可以了。
在找到这些角点之后我们可以使用函数 cv2.cornerSubPix() 增加准确度。我们使用函数 cv2.drawChessboardCorners() 绘制图案。
2.标定
在得到了这些对象点和图像点之后，我们已经准备好来做摄像机标定了。我们要使用的函数是 cv2.calibrateCamera()。它会返回摄像机矩阵，畸变系数，旋转和变换向量等。
ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(objpoints, imgpoints, gray.shape[::-1],None,None)3.畸变矫正
现在我们找到我们想要的东西了，我们可以找到一幅图像来对他进行校正了。 OpenCV 提供了两种方法。在那之前我们可以使用从函数 cv2.getOptimalNewCameraMatrix() 得到的自由缩放系数对摄像机矩阵进行优化。如果缩放系数 alpha = 0，返回的非畸变图像会带有最少量的不想要的像素。它甚至有可能在图像角点去除一些像素。如果 alpha = 1，所有的像素都会被返回，还有一些黑图像。它还会返回一个 ROI 图像，我们可以用来对结果进行裁剪。
畸变矫正有两种方法，①使用 cv2.undistort() 这是最简单的方法。只需使用这个函数和上边得到的 ROI 对结果进行裁剪。②使用 remapping 这应该属于“曲线救国”了。首先我们要找到从畸变图像到非畸变图像的映射方程。再使用重映射方程。
这两种方法在下方程序中都有提及。
4.反向投影差
我们可以利用反向投影误差对我们找到的参数的准确性进行估计。得到的结果越接近 0 越好。有了内部参数，畸变参数和旋转变换矩阵，我们就可以使用 cv2.projectPoints() 将对象点转换到图像点。然后就可以计算变换得到图像与角点检测算法的绝对差了。然后我们计算所有标定图像的误差平均值。
完整程序如下所示：
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 planeimages = glob.glob('./sources/*.jpg')
# img = cv2.imread('sources\left12.jpg')for fname in images:img = cv2.imread(fname)gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)# Find the chess board cornersret, 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 cornersimg = cv2.drawChessboardCorners(img, (7,6), corners2,ret)# cv2.imshow('img',img)# cv2.waitKey(500)
cv2.destroyAllWindows()ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(objpoints, imgpoints, gray.shape[::-1],None,None)
# print('mtx:',mtx)
# print('dist:',dist)
img_left12 = cv2.imread('sources\left12.jpg')
h, w = img_left12.shape[:2]
newcameramtx, roi=cv2.getOptimalNewCameraMatrix(mtx,dist,(w,h),1,(w,h))# # way 1: use undistort
# undistort = cv2.undistort(img_left12, mtx, dist, None, newcameramtx)
# # crop the image
# x,y,w,h = roi
# undistort = undistort[y:y+h, x:x+w]
# cv2.imshow('undistort',undistort)
# # cv2.imwrite('calibresult.png',undistort)# way 2: use remapping
# undistort
mapx,mapy = cv2.initUndistortRectifyMap(mtx,dist,None,newcameramtx,(w,h),5)
remapping = cv2.remap(img_left12,mapx,mapy,cv2.INTER_LINEAR)
# crop the image
x,y,w,h = roi
remapping = remapping[y:y+h, x:x+w]
cv2.imshow('remapping',remapping)np.savez('B.npz', bmtx=mtx, bdist=dist, brvecs=rvecs , btvecs=tvecs)mean_error = 0
for i in range(len(objpoints)):imgpoints2, _ = cv2.projectPoints(objpoints[i], rvecs[i], tvecs[i], mtx, dist)error = cv2.norm(imgpoints[i],imgpoints2, cv2.NORM_L2)/len(imgpoints2)mean_error += errorprint("total error: ", mean_error/len(objpoints))
cv2.waitKey(0)
cv2.destroyAllWindows()
# cv2.imwrite('calibresult.png',remapping)注：像我一样没有系统学习过python的人需要学习一下将畸变矫正矩阵保存至代码外文件的方法，从而便于其他程序调用。上方程序使用的是np.savez保存为npz文件，此外还可以保存为yaml文件。
二、姿势估计
在上一节的摄像机标定中，我们已经得到了摄像机矩阵，畸变系数等。有了这些信息我们就可以估计图像中图案的姿势，比如目标对象是如何摆放，如何旋转等。对一个平面对象来说，我们可以假设 Z=0，这样问题就转化成摄像机在空间中是如何摆放（然后拍摄）的。所以，如果我们知道对象在空间中的姿势，我们就可以在图像中绘制一些 2D 的线条来产生 3D 的效果。我们来看一下怎么做吧。
我们的问题是，在棋盘的第一个角点绘制 3D 坐标轴（ X， Y， Z 轴）。 X轴为蓝色， Y 轴为绿色， Z 轴为红色。在视觉效果上来看， Z 轴应该是垂直与棋盘平面的。
首先我们要加载前面结果中摄像机矩阵和畸变系数。
import cv2
import numpy as np
import glob# 读取方法一
with np.load('B.npz') as X:mtx, dist, _, _ = [X[i] for i in ('bmtx','bdist','brvecs','btvecs')]
# 读取方法二
# arr = np.load('B.npz')
# mtx = arr['bmtx']
# dist = arr['bdist']
# 与之前文件对比查看是否一致
# print('mtx:',mtx)
# print('dist:',dist)我们来创建一个函数： draw，它的参数有棋盘上的角点（使用cv2.findChessboardCorners() 得到）和要绘制的 3D 坐标轴上的点。
def draw(img, corners, imgpts):corner = tuple(corners[0].ravel())corner = tuple(map(int, corner))# corner = tuple(map(int, corner))img = cv2.line(img, corner, tuple(map(int,imgpts[0].ravel())), (255,0,0), 5)img = cv2.line(img, corner, tuple(map(int,imgpts[1].ravel())), (0,255,0), 5)img = cv2.line(img, corner, tuple(map(int,imgpts[2].ravel())), (0,0,255), 5)return imgPs：对于这一段程序，需要注意原教材中错误程序未将参数元组中每个元素转化为int型，原程序如下，可以上下对比着学习tuple转为int的方法：
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() 从图像上的第一个角点开始绘制连接这些点的直线了。搞定！！！
完整程序如下：
import cv2
import numpy as np
import glob# 读取方法一
with np.load('B.npz') as X:mtx, dist, _, _ = [X[i] for i in ('bmtx','bdist','brvecs','btvecs')]
# 读取方法二
# arr = np.load('B.npz')
# mtx = arr['bmtx']
# dist = arr['bdist']
# 查看是否一致
# print('mtx:',mtx)
# print('dist:',dist)def draw(img, corners, imgpts):corner = tuple(corners[0].ravel())corner = tuple(map(int, corner))# corner = tuple(map(int, corner))img = cv2.line(img, corner, tuple(map(int,imgpts[0].ravel())), (255,0,0), 5)img = cv2.line(img, corner, tuple(map(int,imgpts[1].ravel())), (0,255,0), 5)img = cv2.line(img, corner, tuple(map(int,imgpts[2].ravel())), (0,0,255), 5)return imgcriteria = (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)for fname in glob.glob('./sources/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)# rvecs, tvecs, inliers= cv2.solvePnPRansac(objp, corners2, mtx, dist)# project 3D points to image planeimgpts, jac = cv2.projectPoints(axis, rvecs, tvecs, mtx, dist)# corners2 = int(corners2)img = draw(img,corners2,imgpts)cv2.imshow('img',img)k = cv2.waitKey(0) & 0xffif k == 's':cv2.imwrite(fname[:6]+'.png', img)
cv2.destroyAllWindows()三、对极几何
在我们使用针孔相机时，我们会丢失大量重要的信心，比如说图像的深度，或者说图像上的点和摄像机的距离，因这是一个从 3D 到 2D 的转换。因此一个重要的问题就产生了，使用这样的摄像机我们能否计算除深度信息呢？答案就是使用多个相机。我们的眼睛就是这样工作的，使用两个摄像机（两个眼睛），这被称为立体视觉。
在进入深度图像之前，我们要先掌握一些多视角几何的基本概念。在本节中我们要处理对极几何。下图为使用两台摄像机同时对一个一个场景进行拍摄的示意图。
如果只是用一台摄像机我们不可能知道 3D 空间中的 X 点到图像平面的距离，因为 OX 连线上的每个点投影到图像平面上的点都是相同的。但是如果我们也考虑上右侧图像的话，直线 OX 上的点将投影到右侧图像上的不同位置。所以根据这两幅图像，我们就可以使用三角测量计算出 3D 空间中的点到摄像机的距离（深度）。这就是整个思路。直线 OX 上的不同点投射到右侧图像上形成的线 l′ 被称为与 x 点对应的极线。也就是说，我们可以在右侧图像中沿着这条极线找到 x 点。它可能在这条直线上某个位置（这意味着对两幅图像间匹配特征的二维搜索就转变成了沿着极线的一维搜索。这不仅节省了大量的计算，还允许我们排除许多导致虚假匹配的点）。这被称为对极约束。与此相同，所有的点在其他图像中都有与之对应的极线。平面 XOO' 被称为对极平面。O 和 O' 是摄像机的中心。从上面的示意图可以看出，右侧摄像机的中心O' 投影到左侧图像平面的 e 点，这个点就被称为极点。极点就是摄像机中心连线与图像平面的交点。因此点 e' 是左侧摄像机的极点。有些情况下，我们可能不会在图像中找到极点，它们可能落在了图像之外（这说明这两个摄像机不能拍摄到彼此）。
本节我们的重点就是找到极线和极点。为了找到它们，我们还需要两个元素， 本征矩阵（ E） 和基础矩阵（ F） 。本征矩阵包含了物理空间中两个摄像机相关的旋转和平移信息。
为了得到基础矩阵我们应该在两幅图像中找到尽量多的匹配点。我们可以使用 SIFT 描述符， FLANN 匹配器和比值检测。
完整代码如下所示，需要注意不同版本OpenCV调用SIFT时的细微差距：
import cv2
import numpy as np
from matplotlib import pyplot as pltdef drawlines(img1,img2,lines,pts1,pts2):''' img1 - image on which we draw the epilines for the points in img2lines - corresponding epilines '''r,c = img1.shapeimg1 = 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,img2img1 = cv2.imread('sources/left.jpg',0) #queryimage # left image
img2 = cv2.imread('sources/right.jpg',0) #trainimage # right image# 高版本opencv
sift = cv2.SIFT_create()
# 低版本opencv
# 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]# 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()四、深度地图
刚开始学习，尚未知这一知识点在实践中的作用，所以就不记录了。
【笔记】摄像机标定和3D重构相关推荐
- [OpenCV-Python] OpenCV 中摄像机标定和 3D 重构 部分 VII
部分 VII 摄像机标定和 3D 重构 OpenCV-Python 中文教程(搬运)目录 42 摄像机标定 目标 • 学习摄像机畸变以及摄像机的内部参数和外部参数 • 学习找到这些参数,对畸变图像进行 ...
- 十五天掌握OpenCV——摄像机标定和3D重构！—摄像机标定
魏老师学生--Cecil:学习OpenCV-机器视觉之旅 基础 代码 设置 标定 畸变校正 反向投影误差 代码演示 Aim: 学习摄像机畸变以及摄像机的内部参数和外部参数: 对畸变图像进行修复. 基础 ...
- OpenCV-Python] OpenCV 中摄像机标定和 3D 重构 部分 VII
https://www.cnblogs.com/Undo-self-blog/p/8448500.html 42 摄像机标定 目标 • 学习摄像机畸变以及摄像机的内部参数和外部参数 • 学习找到这些参 ...
- Opencv相机标定与3D重构---使用棋盘格来进行摄像机标定
让我们写一点代码来检测在一幅图像中的棋盘格,并获取他到摄像机的距离. 你可以使用同样的方法来针对任何已知三维几何结构的物体,这个物体可以在一幅图像中被检测到. 创建一个空的控制台项目. 载入一幅图片: ...
- 2019-9-29 opencv摄像机标定与三维重构4-Depth Map from Stereo Images立体图像中的深度图（视差图）
官网参见https://docs.opencv.org/3.4.1/dd/d53/tutorial_py_depthmap.html 上一节中,我们学习了极线约束的概念和相关术语.主要包含:如果我们有 ...
- 双目立体视觉源代码 双目立体视觉匹配程序 双目视觉3d成像（三维重构图像处理） 基于双目视觉的深度计算和三维重建 opencv写的双目视觉摄像机标定和三维重建代码
双目视觉/双目标定源码/图片集标定匹配三维重建坐标计算OpenCV 1.双目立体视觉源代码(包括标定,匹配,三维重建) 2.双目视觉实验图片集(双目立体视觉中使用的标准实验图,适合初学者进 行实验使用 ...
- 机器视觉学习笔记（4）——单目摄像机标定参数说明
机器视觉学习笔记(4)--单目摄像机标定参数说明 标签: 机器视觉 1.针孔摄像机模型 在介绍摄像机标定参数之前,需要先简单说一下针孔摄像机的原理.投影平面到小孔的距离为焦距f,物体到小孔的距离为Z, ...
- Halcon学习笔记：3D_coordinates（3D标定）
Halcon学习笔记:3D_coordinates(3D标定) 欢迎有兴趣的朋友一起学习,代码理解注释有问题的可以告诉我,一起讨论,共同进步. *初始化程序,dev_close_window() *关 ...
- 【Halcon笔记1】基于Halcon软件的【摄像机标定】以及【内部参数】和【外部参数】的求解过程【原理细节详解】
[1]具体请看Learning OpenCv[摄像机标定] [2]马颂德[摄像机定标] [3]<基于Halcon软件的摄像机标定>论文
最新文章
- Ubuntu 进阶命令——长期不定时更新
- DTO数据传输对象详解
- mysql学习笔记-insert扩展
- 升级php7_PhpStorm 2019.3 发布，全面支持 PHP 7.4
- EJB分布式对象实现方式描述
- Spring源码分析之doDispatch分发请求逻辑
- [原创]jQuery Validation范例
- 【数字信号】基于matlab GUI虚拟信号发生器（各种波形）【含Matlab源码 271期】
- 爬虫学习笔记（四）——糗百爬虫遇到的反爬
- Beyond Compare 激活解决办法
- C语言编写飞机大战程序,C语言实现简单飞机大战
- JS实现合并单元格的两种方法
- 初学者学习网页设计用什么软件最好？
- iOS-关于M1芯片可以下载APP使用问题
- 解决文字与图片始终不并排的问题
- 电商直通车主图设计教程
- 【数学建模暑期培训】配送中心选址问题
- (ECCV-2022)GaitEdge：超越普通的端到端步态识别，提高实用性
- SMBMS(超市订单管理系统)
- LJP Little John PalmOS 1.0 Release 最新版 (RC9后的正式版）