文章目錄

• 單目相機(jī)標(biāo)定(基于Python OpenCV)
• • 1.上期填坑
  • 2.單目相機(jī)標(biāo)定
  • • 2.1 數(shù)據(jù)采集
    • 2.2 角點(diǎn)提取
    • 2.3 參數(shù)求解
    • 2.4 參數(shù)評(píng)估(重投影誤差)
    • 2.5 相機(jī)位姿(棋盤位姿)可視化
    • 2.6 同Matlab標(biāo)定結(jié)果比較

單目相機(jī)標(biāo)定(基于Python OpenCV)

1.上期填坑

在開始本篇博客之前，先填一下上一篇博客【計(jì)算機(jī)視覺】基于ORB角點(diǎn)+RANSAC算法實(shí)現(xiàn)圖像全景拼接的坑(不算填吧，算做一個(gè)記錄，因?yàn)椴]有解決問題，留著看以后有沒有空解決👀)，不想看的可以直接跳到下一節(jié)。

首先解決如何拼接多張圖像(上篇博客只能拼接兩張圖像，多張圖像需要保存兩張圖像的匹配結(jié)果，再重新讀取計(jì)算角點(diǎn))

改進(jìn)的方法的基本思路是，僅計(jì)算當(dāng)前兩張圖像的單應(yīng)變換M_current，通過先前的單應(yīng)矩陣M_before再計(jì)算一個(gè)累積變換，最終通過矩陣乘法得到的M的變換就延續(xù)了當(dāng)前變換的累積變換矩陣:

# 一開始無先前變換，因此設(shè)置為單位陣M_before = np.eye(3,3)result = cv2.imread('datas/1.jpg')# result = CalcUndistort(result, mtx, dist)result,_,_ = utils.auto_reshape(result, 1080)img2 = resultcors2, desc2= extraORBfromImg(orb, img2)for i in range(1,6):print(i)img1 = cv2.imread('datas/'+str(i+1)+'.jpg')# img1 = CalcUndistort(img1, mtx, dist)img1,_,_ = utils.auto_reshape(img1, 1080)cors1, desc1= extraORBfromImg(orb, img1)match_dist, match_idx = ORBMatch(match, desc1, desc2)# 得到匹配點(diǎn)對(duì)的坐標(biāo)match_pts = findMatchCord(match_idx, cors1, cors2)# 可視化匹配點(diǎn)# utils.drawMatches(img1, img2, cors1, cors2, match_idx)# RANSAC迭代去除異常點(diǎn)對(duì)update_match_pts = RANSAC(match_pts)# 最小二乘方法計(jì)算單應(yīng)矩陣M = calc_homography(update_match_pts)# 圖像拼接結(jié)果可視化, 并傳遞累積單應(yīng)變換矩陣M ,返回先前累積拼接結(jié)果resultresult, M = homography_trans(M, M_before, img1, result)M_before = M# 不用再提取一遍:img2 = img1cors2, desc2 = cors1, desc1

值得注意的是，代碼采用的匹配順序是從右至左，為了保證每次只提取最左側(cè)圖像的角點(diǎn)，img1對(duì)于圖像的拍攝時(shí)序應(yīng)該要在img2的左側(cè)。

相應(yīng)的，函數(shù)homography_trans也要做修改:

# 可視化圖像對(duì)映射效果 def homography_trans(M, M_before, img1, img2):M = M_before @ M# out_img 第一張圖像映射到第二張x_min, x_max, y_min, y_max, M2 = calc_border(M, img1.shape)# 透視變換+平移變換(使得圖像在正中央)M = M2 @ Mw, h = int(round(x_max)-round(x_min)), int(round(y_max)-round(y_min))out_img = cv2.warpPerspective(img1, M, (w, h))# print(out_img.shape)# cv2.imshow('ww',out_img)# cv2.waitKey(0)# 調(diào)整兩張圖像位姿一致：# x方向out_img_blank_x = np.zeros((out_img.shape[0], abs(round(x_min)), 3)).astype(np.uint8)img2_blank_x = np.zeros((img2.shape[0], abs(round(x_min)), 3)).astype(np.uint8)if(x_min>0):# print(1)out_img = cv2.hconcat((out_img_blank_x, out_img))if(x_min<0):# print(2)img2 = cv2.hconcat((img2_blank_x, img2))# y方向out_img_blank_y = np.zeros((abs(round(y_min)), out_img.shape[1], 3)).astype(np.uint8)img2_blank_y = np.zeros((abs(round(y_min)), img2.shape[1], 3)).astype(np.uint8)if(y_min>0):# print(3)out_img = cv2.vconcat((out_img, out_img_blank_y))if(y_min<0):# print(4)img2 = cv2.vconcat((img2_blank_y, img2))# 調(diào)整兩張圖像尺度一致(邊緣填充)：if(img2.shape[0]<out_img.shape[0]):blank_y = np.zeros((out_img.shape[0]-img2.shape[0], img2.shape[1], 3)).astype(np.uint8)img2 = cv2.vconcat((img2, blank_y)) else:blank_y = np.zeros((img2.shape[0]-out_img.shape[0], out_img.shape[1], 3)).astype(np.uint8)out_img = cv2.vconcat((out_img, blank_y)) if(img2.shape[1]<out_img.shape[1]):blank_x = np.zeros((img2.shape[0], out_img.shape[1]-img2.shape[1], 3)).astype(np.uint8)img2 = cv2.hconcat((img2, blank_x)) else:blank_x = np.zeros((out_img.shape[0], img2.shape[1]-out_img.shape[1], 3)).astype(np.uint8)out_img = cv2.hconcat((out_img, blank_x)) # cv2.imwrite('out_img.jpg',out_img)# 疊加result = addMatches(out_img, img2)# 圖像背景白mask = 255*np.ones(result.shape).astype(np.uint8)gray_res = cv2.cvtColor(result, cv2.COLOR_BGR2GRAY)mask[gray_res==0]=0cv2.imwrite('mask.jpg',mask)# result[result==0]=255cv2.imwrite('result.jpg',result)return result, M

重點(diǎn)是改了這里(實(shí)現(xiàn)單應(yīng)變換的傳遞):

... ... M = M_before @ M ... ... M = M2 @ M

(Ps:如果從左至右拼接，好像有點(diǎn)問題…)

然后另一個(gè)問題是:往后的拼接圖像計(jì)算出的單應(yīng)矩陣就會(huì)越扭曲：整一個(gè)視覺效果就會(huì)像這樣

五張圖像的拼接效果:

事實(shí)上，當(dāng)前拼的結(jié)果還不完整，但一整個(gè)圖像的像素已經(jīng)非常大了，由于單應(yīng)變換導(dǎo)致的圖像扭曲占據(jù)了較大的圖像空間，因此到后面系統(tǒng)直接報(bào)了內(nèi)存分配不足的錯(cuò)誤：

numpy.core._exceptions.MemoryError: Unable to allocate 7.71 GiB for an array with shape (17705, 19486, 3) and data type float64

還記得的上一篇博客提出的解決方案是，利用相機(jī)矩陣的畸變矩陣對(duì)圖像進(jìn)行一個(gè)徑向畸變扭曲，目前本人已經(jīng)初略了解了相機(jī)矩陣以及畸變參數(shù)，對(duì)圖像做一個(gè)畸變處理(需要自己相機(jī)的內(nèi)參)并不是很大問題，但是發(fā)現(xiàn)仍然不能很好的解決問題:

# 為圖像添加徑向畸變, 避免透視作用下可視化效果不美觀 def CalcUndistort(img, mtx, dist):return cv2.undistort(img, mtx, dist*12, dst=None, newCameraMatrix=None)

少量圖片下，問題不大，但是拼接圖像一多，仍然會(huì)有明顯的透視效果:

原以為是畸變得不夠，增加畸變系數(shù)后，拼接效果反而變差，(可能是因?yàn)閳D像畸變破壞了單應(yīng)變換的基本假設(shè))：

所以，OpenCV, 還得是你

OpenCV為圖像拼接專門封裝了一個(gè)stitching類，最終能夠?qū)崿F(xiàn)可以看的效果的拼接大致需要經(jīng)過以下流程：

(上次博客中提到的論文中的rectangling矩形化方法就是解決整個(gè)圖像拼接中扭曲圖像這一小步)

(記得課上蔡老師說過，這學(xué)期學(xué)了計(jì)算機(jī)視覺這門課程，就可以去網(wǎng)上接項(xiàng)目賺錢，現(xiàn)在看起來，和小時(shí)候我看個(gè)宇宙紀(jì)錄片，就想成為物理學(xué)家有點(diǎn)類似。當(dāng)然，如果是單純只學(xué)了課上的內(nèi)容)

對(duì)于計(jì)算機(jī)視覺而言，實(shí)踐是必要的，如果沒有親手實(shí)現(xiàn)，單純的被動(dòng)輸入或許永遠(yuǎn)也不會(huì)知道某種看似可行的方法會(huì)存在的缺陷。

關(guān)于具體實(shí)現(xiàn)或者相關(guān)函數(shù)的API，可以觀摩這位大佬的博客OpenCV總結(jié)3——圖像拼接Stitching，我上面那張圖也是從他那偷的，respect！

2.單目相機(jī)標(biāo)定

進(jìn)入正題

在圖像測(cè)量過程以及機(jī)器視覺應(yīng)用中，為確定空間物體表面某點(diǎn)的三維幾何位置與其在圖像中對(duì)應(yīng)點(diǎn)之間的相互關(guān)系，必須建立相機(jī)成像的幾何模型，這些幾何模型參數(shù)就是相機(jī)參數(shù)。在大多數(shù)條件下這些參數(shù)必須通過實(shí)驗(yàn)與計(jì)算才能得到，這個(gè)求解參數(shù)（內(nèi)參、外參、畸變參數(shù)）的過程就稱之為相機(jī)標(biāo)定（或攝像機(jī)標(biāo)定）。無論是在圖像測(cè)量或者機(jī)器視覺應(yīng)用中，相機(jī)參數(shù)的標(biāo)定都是非常關(guān)鍵的環(huán)節(jié)，其標(biāo)定結(jié)果的精度及算法的穩(wěn)定性直接影響相機(jī)工作產(chǎn)生結(jié)果的準(zhǔn)確性。因此，做好相機(jī)標(biāo)定是做好后續(xù)工作的前提，提高標(biāo)定精度是科研工作的重點(diǎn)所在。

OpenCV封裝了單目相機(jī)標(biāo)定所需的各種方法。本篇博客的實(shí)現(xiàn)將基于這些封裝好的函數(shù)。

本次單目相機(jī)的標(biāo)定方法基于張正友博士于1998年提出的棋盤格標(biāo)定方法,

論文鏈接：A Flexible New Technique for Camera Calibration，收錄于2000年P(guān)AMI期刊。

具體的數(shù)學(xué)原理以及公式的推導(dǎo)可以參考網(wǎng)絡(luò)上各類詳細(xì)的博客以及教程，本篇博客不再贅述。

2.1 數(shù)據(jù)采集

標(biāo)定板規(guī)格: GP340 12x9 25mm

相機(jī): HUAWEI MATE 30 PRO（LIO AL00）

焦距: 27mm

分辨率: 3648x2736

拍攝數(shù)量: 20

2.2 角點(diǎn)提取

(獲取角點(diǎn)兩個(gè)坐標(biāo)系下的坐標(biāo)，作為求解方程中的已知量)

首先，基于cv2.findChessboardCorners尋找棋盤格角點(diǎn)：

# 定位角點(diǎn) def find_chessboard_cor(img):# 轉(zhuǎn)為灰度圖gray_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)# OpenCV內(nèi)置函數(shù)提取棋盤格角點(diǎn), (11,8)為棋盤格尺寸-1(12x9)is_success, corner = cv2.findChessboardCorners(gray_img, (11,8), None)# 計(jì)算亞像素時(shí)停止迭代的標(biāo)準(zhǔn)# 后者表示迭代次數(shù)達(dá)到了最大次數(shù)時(shí)停止，前者表示角點(diǎn)位置變化的最小值已經(jīng)達(dá)到最小時(shí)停止迭代criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 100, 0.001)# 亞像素角點(diǎn)檢測(cè)，提高角點(diǎn)定位精度, (7, 7)為考慮角點(diǎn)周圍區(qū)域的大小corner = cv2.cornerSubPix(gray_img, corner, (7, 7), (-1, -1), criteria) return is_success, corner

值得注意的是，這里我們必須使用cv2.cornerSubPix再進(jìn)行一次迭代定位亞像素角點(diǎn)，使得角點(diǎn)的定位更為精確。

基于cv2.drawChessboardCorners可視化棋盤格角點(diǎn)：

# 可視化角點(diǎn) def draw_chessboard_cor(img, cor, is_success):cv2.drawChessboardCorners(img, (11,8), cor, is_success)cv2.imshow('cor', img)cv2.waitKey(50)

OpenCV提取棋盤角點(diǎn)可視化:

別忘了我們提取了棋盤格角點(diǎn)坐標(biāo)的目的：這一步是為了獲取角點(diǎn)在圖像坐標(biāo)系中的坐標(biāo)。為了求解相機(jī)內(nèi)外參矩陣，緊接著我們還需要獲取棋盤格在原始世界坐標(biāo)系中的坐標(biāo)。

值得一提的是，無論是固定相機(jī)移動(dòng)棋盤格進(jìn)行拍攝，還是固定棋盤格移動(dòng)相機(jī)進(jìn)行拍攝，在張正友標(biāo)定法的假設(shè)中，棋盤格永遠(yuǎn)處在世界坐標(biāo)系Z=0的平面上，棋盤格最左上角的角點(diǎn)為世界坐標(biāo)系原點(diǎn)，往右為x軸正方向，往下為y軸正方向：

由于同一張標(biāo)定板的規(guī)格不會(huì)改變，因此，世界坐標(biāo)系通過人為定義即可:

# 注:相機(jī)參數(shù)的計(jì)算只要求角點(diǎn)之間的世界坐標(biāo)比例一致,因此可以單位化world_coord = np.zeros((w * h, 3), np.float32)world_coord[:, :2] = np.mgrid[0:w, 0:h].T.reshape(-1, 2)

世界坐標(biāo)(單位化):

[[ 0. 0. 0.]
[ 1. 0. 0.]
[ 2. 0. 0.]
[ 3. 0. 0.]
[ 4. 0. 0.]
[ 5. 0. 0.]
[ 6. 0. 0.]
[ 7. 0. 0.]

… …

]]

shape = (88, 3) # 每個(gè)棋盤格上提取11*8=88個(gè)角點(diǎn)

一個(gè)特別值得注意且容易迷惑的地方是，由于相機(jī)參數(shù)的求解只依賴于角點(diǎn)與角點(diǎn)之間世界坐標(biāo)的比例關(guān)系以及圖像的大小(分辨率)，而與具體的尺度無關(guān)，因此，角點(diǎn)世界坐標(biāo)系的設(shè)置無需嚴(yán)格按照實(shí)際的尺度(25mm)，只需按照角點(diǎn)之間距離的比例(單位化)(相鄰角點(diǎn)的x，y坐標(biāo)差均為1:1)

如果不信可以試試將角點(diǎn)的世界坐標(biāo)賦予不同但比例一致的值，最終的求解結(jié)果都是完全一致的(len可以表示具體的尺度):

# world_coord[:,:2] = np.mgrid[0:w*len:len,0:h*len:len].T.reshape(-1,2)

所以說，后續(xù)求得的相機(jī)內(nèi)外參數(shù)也并不代表實(shí)際的尺度，它的單位是像素，如果要轉(zhuǎn)化為實(shí)際的尺度，還需要知道相機(jī)每個(gè)像素代表現(xiàn)實(shí)中的實(shí)際尺度大小。

2.3 參數(shù)求解

通過cv2.calibrateCamera求解相機(jī)內(nèi)外參數(shù)

# 求解攝像機(jī)的內(nèi)在參數(shù)和外在參數(shù)# ret 非0表示標(biāo)定成功 mtx 內(nèi)參數(shù)矩陣，dist 畸變系數(shù)，rvecs 旋轉(zhuǎn)向量，tvecs 平移向量# 注:求解的結(jié)果的單位為像素,若想化為度量單位還需乘上每個(gè)像素代表的實(shí)際尺寸(如:毫米/像素)ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(world, cam, (img.shape[1], img.shape[0]), None, None)

當(dāng)然，只基于一張圖像的標(biāo)定結(jié)果并不可靠，因此，需要自定義一個(gè)函數(shù)進(jìn)行批量操作：

# 求解內(nèi)外參數(shù) def CamCalibrate(w, h, num, root):# 圖像縮放比例(如果你的圖像進(jìn)行了縮放，與實(shí)際拍攝的分辨率不一致，最終求得的參數(shù)需要乘上這個(gè)比例進(jìn)行校正)ratio = 1 # 3648 / 1920 world, cam = [], []# 多張圖像進(jìn)行標(biāo)定，減小誤差:for i in range(num):img = cv2.imread(root + str(i+1)+'.jpg')# img,_,_ = utils.auto_reshape(img, 1920)# 定位角點(diǎn)is_success, cam_coord = find_chessboard_cor(img)print('第'+str(i+1)+'張角點(diǎn)提取完畢, 角點(diǎn)數(shù) =',cam_coord.shape[0])# 可視化角點(diǎn)# draw_chessboard_cor(img, cam_coord, is_success)# 角點(diǎn)的世界坐標(biāo):# 注:相機(jī)參數(shù)的計(jì)算只要求角點(diǎn)之間的世界坐標(biāo)比例一致,因此可以單位化world_coord = np.zeros((w * h, 3), np.float32)world_coord[:, :2] = np.mgrid[0:w, 0:h].T.reshape(-1, 2)world_coord[:,1] = -world_coord[:,1]# world_coord[:,:2] = np.mgrid[0:w*len:len,0:h*len:len].T.reshape(-1,2)# 將世界坐標(biāo)與像素坐標(biāo)加入待求解系數(shù)矩陣world.append(world_coord)cam.append(cam_coord)# 求解攝像機(jī)的內(nèi)在參數(shù)和外在參數(shù)# ret 非0表示標(biāo)定成功 mtx 內(nèi)參數(shù)矩陣，dist 畸變系數(shù)，rvecs 旋轉(zhuǎn)向量，tvecs 平移向量# 注:求解的結(jié)果的單位為像素,若想化為度量單位還需乘上每個(gè)像素代表的實(shí)際尺寸(如:毫米/像素)ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(world, cam, (img.shape[1], img.shape[0]), None, None)rvecs = np.array(rvecs).reshape(-1,3)tvecs = np.array(tvecs).reshape(-1,3)# 單位:像素(1像素=??mm)print("標(biāo)定結(jié)果 ret:", ret)print("內(nèi)參矩陣 mtx:\n", mtx) # 內(nèi)參數(shù)矩陣print("畸變系數(shù) dist:\n", dist) # 畸變系數(shù) distortion cofficients = (k_1,k_2,p_1,p_2,k_3)print("旋轉(zhuǎn)向量(外參) rvecs:\n", rvecs) # 旋轉(zhuǎn)向量 # 外參數(shù)(歐拉角)print("平移向量(外參) tvecs:\n", tvecs) # 平移向量 # 外參數(shù)np.save('./param/mtx.npy',mtx)np.save('./param/dist.npy',dist)np.save('./param/rvecs.npy',rvecs)np.save('./param/tvecs.npy',tvecs)np.save('./param/world.npy',np.array(world))np.save('./param/cam.npy',np.array(cam))return ret, mtx, dist, rvecs, tvecs

最終結(jié)果:

標(biāo)定結(jié)果 ret: 1.4146770040984205
內(nèi)參矩陣 mtx:
[[2.88760072e+03 0.00000000e+00 1.82151647e+03]
[0.00000000e+00 2.88741538e+03 1.37233666e+03]
[0.00000000e+00 0.00000000e+00 1.00000000e+00]]
畸變系數(shù) dist:
[[ 0.01488259 0.03388964 -0.00044845 -0.00178608 0.03024054]]
旋轉(zhuǎn)向量(外參) rvecs:
[[ 2.77705978 -0.46989843 0.1686412 ]
[-2.8145649 -0.91512141 -0.43737601]
… …
[ 2.67447733 0.66658815 -0.11254641]
[ 2.818184 -0.02291704 0.51413859]]
平移向量(外參) tvecs:
[[-5.06957615 -2.43006408 16.33928058]
[ 0.29829902 -4.17287918 13.05665933]
… …
[-3.8811355 -4.39040007 19.95916291]
[-4.53034049 -4.16521408 14.12011733]]

當(dāng)然，單目相機(jī)標(biāo)定的流程到這里就可以結(jié)束了

但是，評(píng)估求得的相機(jī)矩陣是否可靠也同樣重要。一個(gè)方法是計(jì)算重投影誤差，將角點(diǎn)的世界坐標(biāo)通過已經(jīng)求得的外參矩陣和相機(jī)內(nèi)參矩陣重新轉(zhuǎn)化到像素坐標(biāo)下，與原始提取的圖像角點(diǎn)坐標(biāo)計(jì)算RMSE，得到的結(jié)果作為重投影誤差，誤差越小，結(jié)果越可靠。

2.4 參數(shù)評(píng)估(重投影誤差)

我們可以使用OpenCV內(nèi)置的cv2.projectPoints()直接計(jì)算重投影誤差

當(dāng)然，也可以自己實(shí)現(xiàn)：

# 計(jì)算重投影坐標(biāo) def CalcProjectPoints(world, rvecs, tvecs, mtx, dist):# 旋轉(zhuǎn)向量轉(zhuǎn)旋轉(zhuǎn)矩陣M = Rodriguez(rvecs)# c = Rw + t (世界坐標(biāo)系轉(zhuǎn)相機(jī)坐標(biāo)系)R_t = (M @ world.T).T + tvecs# (相機(jī)坐標(biāo)系到圖像坐標(biāo)系)# print(R_t)plain_pts = (mtx @ R_t.T)plain_pts = (plain_pts / plain_pts[2,:]).T[:,:2]# 去畸變c_xy = np.array([mtx[0,2], mtx[1,2]])f_xy = np.array([mtx[0,0], mtx[1,1]])k1, k2, p1, p2, k3 = dist[0]x_y = (plain_pts - c_xy) / f_xyr = np.sum(x_y * x_y, 1)x_distorted = x_y[:,0] * (1 + k1*r + k2*r*r + k3*r*r*r) + 2*p1*x_y[:,0]*x_y[:,1] + p2*(r + 2*x_y[:,0]*x_y[:,0])y_distorted = x_y[:,1] * (1 + k1*r + k2*r*r + k3*r*r*r) + 2*p2*x_y[:,0]*x_y[:,1] + p1*(r + 2*x_y[:,1]*x_y[:,1])u_distorted = f_xy[0]*x_distorted + c_xy[0]v_distorted = f_xy[1]*y_distorted + c_xy[1]plain_pts = np.array([u_distorted,v_distorted]).Treturn plain_pts

這里面有很多細(xì)節(jié)需要了解，就比如OpenCV的旋轉(zhuǎn)向量轉(zhuǎn)旋轉(zhuǎn)矩陣基于Rodriguez變換(我之前一直以為是歐拉角，算出來的結(jié)果讓我迷了一段時(shí)間):

# 旋轉(zhuǎn)向量轉(zhuǎn)旋轉(zhuǎn)矩陣 def Rodriguez(rvecs):# 旋轉(zhuǎn)向量模長(zhǎng)θ = (rvecs[0] * rvecs[0] + rvecs[1] * rvecs[1] + rvecs[2] * rvecs[2])**(1/2)# 旋轉(zhuǎn)向量的單位向量r = rvecs / θ# 旋轉(zhuǎn)向量單位向量的反對(duì)稱矩陣anti_r = np.array([[0, -r[2], r[1]],[r[2], 0, -r[0]],[-r[1], r[0], 0]])# 旋轉(zhuǎn)向量轉(zhuǎn)旋轉(zhuǎn)矩陣(Rodriguez公式) # np.outer(r, r) = r @ r.T 向量外積M = np.eye(3) * np.cos(θ) + (1 - np.cos(θ)) * np.outer(r, r) + np.sin(θ) * anti_rreturn M

還有就是如何重投影時(shí)添加畸變，這些網(wǎng)上都可以查閱得到，這里不再贅述。

OpenCV同樣內(nèi)置cv2.norm()求解RMSE，當(dāng)然也可以自己實(shí)現(xiàn)：

# RMSEerror = original_pts - reprojec_ptserror = np.sum(error*error)**(1/2) / reproj.shape[0]

我們也可以將重投影結(jié)果在圖像中可視化出來：

# 重投影可視化 def drawReprojCor(img, original_pts, reprojec_pts, idx):r,g = (0,0,255),(0,255,0)for i in range(original_pts.shape[0]):# 原始角點(diǎn)x0, y0 = int(round(original_pts[i,0])), int(round(original_pts[i,1]))cv2.circle(img, (x0, y0), 10, g, 2, lineType=cv2.LINE_AA)# 重投影角點(diǎn)x1, y1 = int(round(reprojec_pts[i,0])), int(round(reprojec_pts[i,1]))cv2.circle(img, (x1, y1), 10, r, 2, lineType=cv2.LINE_AA)cv2.imwrite('./reproject_result/'+ str(idx+1)+'.jpg', img)

紅圈表示重投影點(diǎn)，綠圈表示原始提取的角點(diǎn)，來可視化兩張誤差比較大的：

對(duì)應(yīng)序列第6張: RMSE=0.2501498893076572

對(duì)應(yīng)序列第9張:RMSE=0.32640338317469253

同樣，封裝成批處理函數(shù)：

# 計(jì)算重投影誤差: def CalcReprojError(num, world, cam, mtx, dist, rvecs, tvecs):errors = []for i in range(num):# 計(jì)算重投影坐標(biāo)(dist=0不考慮畸變)# reproj, _ = cv2.projectPoints(world[i], rvecs[i], tvecs[i], mtx, dist)reproj = CalcProjectPoints(world[i], rvecs[i], tvecs[i], mtx, dist)# 原始坐標(biāo)original_pts = cam[i].reshape(cam[i].shape[0], 2)# 重投影坐標(biāo)reprojec_pts = reproj.reshape(reproj.shape[0], 2)# RMSEerror = original_pts - reprojec_ptserror = np.sum(error*error)**(1/2) / reproj.shape[0]# 等價(jià):# error = cv2.norm(cam[i],reproj, cv2.NORM_L2) / reproj.shape[0]errors.append(error)# 重投影可視化img = cv2.imread('./0/' + str(i+1)+'.jpg')# img,_,_ = utils.auto_reshape(img, 1920)drawReprojCor(img, original_pts, reprojec_pts, i)

如果想讓可視化更直觀些，可以用條形圖統(tǒng)計(jì)每一張的誤差：

# 誤差條形圖# 可視化每張圖像的誤差(單位:像素)plt.bar(range(num), errors, width=0.8, label='reproject error', color='#87cefa')# 誤差平均值(單位:像素)mean_error = sum(errors) / numplt.plot([-1,num], [mean_error,mean_error], color='r', linestyle='--',label='overall RMSE:%.3f'%(mean_error))plt.xticks(range(num), range(1,num+1))plt.ylabel('RMSE Error in Pixels')plt.xlabel('Images')plt.legend()

那這樣標(biāo)定的完整流程應(yīng)該就可以結(jié)束了吧

好吧，求解的參數(shù)中還有外參數(shù)沒用，為什么不來可視化看看相機(jī)位姿呢👀

2.5 相機(jī)位姿(棋盤位姿)可視化

可視化相機(jī)位姿的核心思路就是將相機(jī)坐標(biāo)系下的相機(jī)坐標(biāo)通過外參矩陣(旋轉(zhuǎn)+平移)轉(zhuǎn)換到世界坐標(biāo)系：

# 可視化標(biāo)定過程中的相機(jī)位姿 def show_cam_pose(rvecs, tvecs):# 相機(jī)坐標(biāo)系下基向量vec = np.array([[1,0,0],[0,1,0],[0,0,1]]) # 相機(jī)位姿模型cam = (2/3)*np.array([[ 1, 1, 2],[-1, 1, 2],[-1,-1, 2],[ 1,-1, 2],[ 1, 1, 2],])fig = plt.figure()ax = fig.add_subplot(111, projection='3d')# 繪制每一個(gè)角度拍攝的相機(jī)for i in range(rvecs.shape[0]): # 旋轉(zhuǎn)向量轉(zhuǎn)旋轉(zhuǎn)矩陣M = Rodriguez(rvecs[i,:])# 相機(jī)原點(diǎn) w = R^(-1)(0 - t)x0,y0,z0 = M.T @ (-tvecs[i,:])c = ['r','g','b']# 隨機(jī)顏色hex = '0123456789abcdef'rand_c = '#'+''.join([hex[random.randint(0,15)] for _ in range(6)])# 繪制相機(jī)坐標(biāo)系for j in range(3):# 相機(jī)位姿(相機(jī)坐標(biāo)系轉(zhuǎn)世界坐標(biāo)系)# w = R^(-1)(c - t)x1,y1,z1 = M.T @ (vec[j,:] - tvecs[i,:])# 相機(jī)坐標(biāo)系ax.plot([x0,x1],[y0,y1],[z0,z1],color=c[j])C = (M.T @ (cam - tvecs[i,:]).T).T # 繪制相機(jī)位姿for k in range(4):ax.plot([C[k,0],C[k+1,0]],[C[k,1],C[k+1,1]],[C[k,2],C[k+1,2]], color=rand_c)ax.plot([x0,C[k+1,0]],[y0,C[k+1,1]],[z0,C[k+1,2]], color=rand_c)# 相機(jī)編號(hào)ax.text(x0,y0,z0,i+1)# 繪制棋盤格for i in range(9):ax.plot([0, 11],[-i, -i],[0,0],color="black")for i in range(12):ax.plot([i, i],[-8, 0],[0,0],color="black")# 繪制世界(棋盤格)坐標(biāo)系for i in range(3):ax.plot([0, 3*vec[i,0]],[0, -3*vec[i,1]],[0, 2*vec[i,2]],color=c[i],linewidth=3)plt.xlim(-3,14)plt.ylim(-13,4)

值得注意的是，我們需要在自定義的標(biāo)定函數(shù)中將世界坐標(biāo)系的y軸改為負(fù)的，可視化結(jié)果才是正確的，否則可視化的結(jié)果是錯(cuò)誤的(如下圖)，原因其實(shí)在2.2節(jié)就已給出，大家可以想一想這是為什么

同樣的，可視化棋盤格位姿的核心思路就是將世界坐標(biāo)系下的棋盤坐標(biāo)通過外參矩陣(旋轉(zhuǎn)+平移)的逆變換轉(zhuǎn)換到相機(jī)坐標(biāo)系：

# 可視化標(biāo)定過程中的棋盤位姿 def show_chessboard_pose(rvecs, tvecs):# 棋盤坐標(biāo)系下基向量vec = np.array([[1,0,0],[0,1,0],[0,0,1]]) # 相機(jī)位姿模型C = np.array([[ 1, 1, 3],[-1, 1, 3],[-1,-1, 3],[ 1,-1, 3],[ 1, 1, 3],])fig = plt.figure()ax = fig.add_subplot(111, projection='3d')c = ['r','g','b'] # 坐標(biāo)軸顏色# 繪制每一個(gè)角度拍攝的棋盤格for i in range(rvecs.shape[0]): # 旋轉(zhuǎn)向量轉(zhuǎn)旋轉(zhuǎn)矩陣M = Rodriguez(rvecs[i,:])# 棋盤原點(diǎn) c = Rw + tx0,y0,z0 = tvecs[i,:]# 隨機(jī)顏色hex = '0123456789abcdef'rand_c = '#'+''.join([hex[random.randint(0,15)] for _ in range(6)])# 繪制棋盤位姿for k in range(0,9):b = np.array([[0, -k, 0],[11, -k, 0]])b = (M @ b.T + tvecs[i,:].reshape(-1,1)).Tax.plot([b[0,0], b[1,0]],[b[0,1], b[1,1]],[b[0,2], b[1,2]],color=rand_c)for k in range(0,12):b = np.array([[k, -8, 0],[k, 0, 0]])b = (M @ b.T + tvecs[i,:].reshape(-1,1)).Tax.plot([b[0,0], b[1,0]],[b[0,1], b[1,1]],[b[0,2], b[1,2]],color=rand_c)k += 11# 繪制棋盤坐標(biāo)系for j in range(3):# (世界坐標(biāo)系轉(zhuǎn)相機(jī)坐標(biāo)系)# c = Rw + tx1,y1,z1 = M @ vec[j,:] + tvecs[i,:]ax.plot([x0,x1],[y0,y1],[z0,z1],color=c[j])# 棋盤編號(hào)ax.text(x0,y0,z0,i+1)# 繪制世界(相機(jī))坐標(biāo)系for i in range(3):ax.plot([0, vec[i,0]],[0, -vec[i,1]],[0, 2*vec[i,2]],color=c[i],linewidth=2)# 繪制相機(jī)位姿for i in range(4):ax.plot([C[i,0],C[i+1,0]],[C[i,1],C[i+1,1]],[C[i,2],C[i+1,2]], color="black")ax.plot([0,C[i+1,0]],[0,C[i+1,1]],[0,C[i+1,2]], color="black")

2.6 同Matlab標(biāo)定結(jié)果比較

一般可以認(rèn)為，Matlab內(nèi)置算法的標(biāo)定結(jié)果更為精確，因此，我們可以再通過Matlab標(biāo)定得到一組參數(shù)，來驗(yàn)證我們的參數(shù)是否靠譜:

基于Matlab的標(biāo)定：

在標(biāo)定過程中，程序自動(dòng)舍棄了兩張誤差較大的圖像，正好就是序列第6張和第9張。

注意到Matlab的Reprojection Errors和我們自己標(biāo)定的重投影誤差有些不太一樣，這是因?yàn)镸atlab使用MSE，我們的代碼使用RMSE。

Matlab給出的標(biāo)定內(nèi)參矩陣(是轉(zhuǎn)置的形式)，和：

我們的標(biāo)定結(jié)果：

內(nèi)參矩陣 mtx:
[[2.88760072e+03 0.00000000e+00 1.82151647e+03]
[0.00000000e+00 2.88741538e+03 1.37233666e+03]
[0.00000000e+00 0.00000000e+00 1.00000000e+00]]

是有一定的誤差但也不至于太離譜，至少可以認(rèn)為結(jié)果是相對(duì)可靠的。

(完)

總結(jié)

以上是生活随笔為你收集整理的【计算机视觉】OpenCV实现单目相机标定的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。