Python--利用張正友法進行相機标定

介紹

 ”張正友标定”是指張正友教授1998年提出的單平面棋盤格的錄影機标定方法。該方法介于攝影标定法和自标定法之間，克服了傳統标定法需要的高精度标定物的缺點，又解決了自标定法魯棒性差的難題。傳統标定法需要三個兩兩垂直的平面，實驗器材在操作過程中難以實作，而僅需使用一個列印出來的棋盤格就可以，并從不同方向拍攝幾組照片即可。相對于自标定而言，提高了精度，便于操作。是以張氏标定法被廣泛應用于計算機視覺方面。

一、基本原理

 相機标定的目的之一是為了建立物體從三維世界到二維成像平面上各坐标點的對應關系，是以首先我們需要定義這樣幾個坐标系：

（1）世界坐标系

世界坐标系是系統的絕對坐标系，是一個三維坐标系，為了描述目标物在真實世界裡的位置而被引入。

（2）相機坐标系

在相機上建立的坐标系，為了從相機的角度描述物體位置而定義，是以光軸與圖像平面的交點為圖像坐标系的原點所構成的直角坐标系。是一個作為溝通世界坐标系和圖像坐标系的中間坐标系。

（3）像平面坐标系

二維坐标系，在一個圖像平面中，以平面的中心像主點P為原點和坐标軸x,y組成了圖檔坐标系。是為了描述成像過程中物體從相機坐标系到圖像坐标系的投影透射關系

（4）像素坐标系

為了描述物體成像後的像點在數字圖像上（相片）的坐标而引入，是真正從相機内讀取到的資訊所在的坐标系。

圖檔來源于

https://blog.csdn.net/a083614/article/details/78579163

 那麼相機标定的一般過程為：

1、從世界坐标系轉為相機坐标系，這一步是三維點到三維點的轉換，包括相機外參R，t等參數；

2、從相機坐标系轉為成像平面坐标系（像素坐标系），這一步是三維點到二維點的轉換，包括相機内參K等參數；

接下來介紹一下實驗所用到的棋盤

在黑白相間的棋盤格上，二維圖像點很容易通過角點檢測找到

實驗中使用的棋盤是129的 棋盤格，每個方格邊長是50mm，即含有118個角點

拍攝工具為：小米8

二、相機标定步驟

1、列印一張棋盤格，把它貼在一個平面上，作為标定物。

2、通過調整标定物或錄影機的方向，為标定物拍攝一些不同方向（這裡拍攝15-20張）的照片。

3、從照片中提取棋盤格角點。（Harris角點）

4、估算理想無畸變的情況下，五個内參和六個外參。

5、應用最小二乘法估算實際存在徑向畸變下的畸變系數。

6、極大似然法，優化估計，提升估計精度。

三、數學原理

（一）外參的計算

設三維世界坐标的點坐标為

二維相機平面的像素點為

齊次坐标為：

是以标定用的棋盤格平面到圖像平面的單應性關系為：

其中：

S：世界坐标系到像平面坐标系的尺度因子

A：相機的内參矩陣

(u0,v0)： 像主點坐标

α, β: 焦距與像素橫縱比的融合

γ: 徑向畸變參數

 因為張正友法标定隻需要一張棋盤格，是以我們不妨假設棋盤格位于Z=0，是以平移向量隻有r1，r2 。。定義旋轉矩陣R的第i列為 ri, 則有：

于是空間到圖像的映射可改為：

 在這裡描述的是空間中平面三維點和相機平面二維點之間的關系。因為相機平面中點的坐标可以通過圖像處理的方式（如Harris角點）擷取，而空間平面中三維點可以通過事先做好的棋盤擷取。是以也就是說每張圖檔都可以計算出一個H矩陣

 其中H 是描述Homographic矩陣。H是一個齊次矩陣，是以有8個未知數，至少需要8個方程，每對對應點能提供兩個方程，是以至少需要四個對應點，就可以算出世界平面到圖像平面的單應性矩陣H

令 H 為 H = [h1 h2 h3]

通過上述等式的矩陣運算，根據正交和歸一化的限制可以得到如下等式：

外部參數可通過Homography求解，由 H = [h1 h2 h3] = λA[r1 r2 t]，可推出

（二）内參的計算

由r1和r2正交，且r1和r2的模相等，可以得到：

正交：

模相等：

B 中的未知量可表示為6D 向量 b，由于B矩陣是個對稱矩陣，是以可以寫成一個6維向量形式：

設H中的第i列為 hi

根據b的定義，可以推導出如下公式

可以推導出

如果有n組觀察圖像，則V 是 2n x 6 的矩陣

根據推到的結果可知如果有n組觀察圖像，則V 是 2n x 6 的矩陣

根據最小二乘定義，V b = 0 的解是 VTV 最小特征值對應的特征向量。

是以, 可以直接估算出 b，可以通過b求解内參

因為B中的未知量為6個，是以當觀測平面 n ≥ 3 時，可以得到b的唯一解

當 n = 2時, 一般可令畸變參數γ = 0

當 n = 1時, 僅能估算出α 與 β, 此時一般可假定像主點坐标 u0 與 v0 為0

B = ATA-1

B 是通過b構造的對稱矩陣

内部參數可通過如下公式計算(cholesky分解)：

(三)最大似然估計

 上述的推導結果是理想情況下的解，但在實際操作過程中，由于可能存在噪聲等一系列的影響，為了增加結果的可靠性，可以使用最大似然估計(Maximum likelihood estimation)來對上述的結果進行優化。

 假設同一個相機從n個不同的角度采集了n副包含棋盤格的圖像進行定标，每個圖像裡有棋盤格角點m個。Mij表示第i幅圖像上第j個像點對應的标定闆上的三維點，則：

其中Ri和ti是第i副圖對應的旋轉矩陣和平移向量，K是内參數矩陣

則角點mij的機率密度函數為：

構造似然函數：

為了能夠讓L取得最大值，需要最小化下面的值

這裡使用的是多參數非線性系統優化問題的Levenberg-Marquardt算法進行疊代求最優解。

（四）畸變

 拍照時通常要在相機的鏡頭前添加透鏡，在相機成像的過程中，透鏡會對光線的傳播産生影響，進而影響相機的成像效果，産生了畸變

 畸變主要分為兩種：徑向畸變、切向畸變（畸變是相機本身的特性，與内參無關，标定一次後即可）

徑向畸變來自于透鏡形狀。

切向畸變來自于整個錄影機的組裝過程。

1、徑向畸變：

 透鏡自身的形狀對才光線的傳播産生影響，實際錄影機的透鏡總是在成像儀的邊緣産生顯著的畸變，這種現象來源于“筒形”或“魚眼”的影響。

 在針孔相機模型中，一條直線在成像平面上的像仍然是直線。但是在實際拍攝的過程中，由于透鏡的存在，一條直線往往被投影成了曲線，當拍攝的角度越靠近圖像的邊緣，這種現象越明顯。因為透鏡是中心對稱的，是以這種不規則的畸變通常是徑向對稱的。主要有兩大類：桶形畸變和枕形畸變。如圖

對于徑向畸變可以通過下面的泰勒級數展開式進行校正。這裡不做闡述。

2、切向畸變：

由于透鏡制造上的缺陷使得透鏡本身與圖像平面不平行而産生的。

3、徑向畸變估計

參考部落格

https://www.cnblogs.com/wangguchangqing/p/8335131.html#autoid-2-3-0

 在張氏定标法中，隻關注了影響最大的徑向畸變。

假設，(μ,ν)是理想的無畸變的像素坐标；(μ^,ν ^)是畸變後的像素坐标；(μ0,ν0)是相機的像主點；(x,y)和(x ^,y ^)理想的無畸變的歸一化的圖像坐标和畸變後的歸一化圖像坐标，數學表達為：

其中，k1,k2表示徑向畸變的系數。

矩陣形式：

記做：Dk=d

則可得：

計算得到畸變系數k。

使用最大似然的思想優化得到的結果，即像上一步一樣，LM法計算下列函數值最小的參數值：

得到畸變參數k1,k2後，可以先将圖像進行去畸變處理，然後用去畸變後的圖像坐标估計相機的内參數。

四、實作

（一） matlab實作

1、應用程式中找到Camera Calibration

2、添加标定闆拍攝圖檔

3、輸入棋盤格每格的尺寸大小

4、點選Calibration，開始标定。

5、得到标定結果

6、結果

每一幅圖像都有一個外參矩陣，它通過将棋盤格上坐标系上的角點坐标（x,y,0,1）（計算過程中設定參考坐标系在棋盤格上左上方角點，Z軸垂直棋盤格平面，是以z為0,1表示齊次坐标）左乘外參矩陣再左乘内參矩陣可以得到該外參矩陣對應圖像中的相應角點坐标

（二） python實作

# -*- coding: cp936 -*-
import cv2
import numpy as np
import glob

# 設定尋找亞像素角點的參數，采用的停止準則是最大循環次數30和最大誤差容限0.001
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)

# 擷取标定闆角點的位置
objp = np.zeros((5*9,3), np.float32)
objp[:,:2] = np.mgrid[0:9,0:5].T.reshape(-1,2)  # 将世界坐标系建在标定闆上，所有點的Z坐标全部為0，是以隻需要指派x和y

obj_points = []  # 存儲3D點
img_points = []  # 存儲2D點

images = glob.glob(r"C:\Users\WeiLinLin\PycharmProjects\untitled\data\chess\*.jpg")
i=0;
for fname in images:
    img = cv2.imread(fname)
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    size = gray.shape[::-1]
    ret, corners = cv2.findChessboardCorners(gray, (11, 8), None)
    #print(corners)

    if ret:

        obj_points.append(objp)

        corners2 = cv2.cornerSubPix(gray, corners, (5, 5), (-1, -1), criteria)  # 在原角點的基礎上尋找亞像素角點
        #print(corners2)
        if [corners2]:
            img_points.append(corners2)
        else:
            img_points.append(corners)

        cv2.drawChessboardCorners(img, (11, 8), corners, ret)  # 記住，OpenCV的繪制函數一般無傳回值
        i+=1;
        cv2.imwrite('conimg'+str(i)+'.jpg', img)
        cv2.waitKey(4000)

print(len(img_points))
cv2.destroyAllWindows()

# 标定
ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(obj_points, img_points, size, None, None)

print("ret:", ret)
print("mtx:\n", mtx) # 内參數矩陣
print("dist:\n", dist)  # 畸變系數   distortion cofficients = (k_1,k_2,p_1,p_2,k_3)
print("rvecs:\n", rvecs)  # 旋轉向量  # 外參數
print("tvecs:\n", tvecs ) # 平移向量  # 外參數

print("-----------------------------------------------------")

img = cv2.imread(images[2])
h, w = img.shape[:2]
newcameramtx, roi = cv2.getOptimalNewCameraMatrix(mtx,dist,(w,h),1,(w,h))#顯示更大範圍的圖檔（正常重映射之後會删掉一部分圖像）
print (newcameramtx)
print("------------------使用undistort函數-------------------")
dst = cv2.undistort(img,mtx,dist,None,newcameramtx)
x,y,w,h = roi
dst1 = dst[y:y+h,x:x+w]
cv2.imwrite('calibresult3.jpg', dst1)
print ("dst的大小為:", dst1.shape)
           

實驗結果：

mtx為内參矩陣，dist為畸變參數：

rvecs為旋轉向量：

畸變矯正：