錄影機标定(Camera calibration)筆記

一 作用

建立3D到2D的映射關系，一旦标定後，對于一個錄影機内部參數K(光心焦距變形參數等，簡化的情況是隻有f錯切=0,變比=1，光心位置簡單假設為圖像中心)，參數已知，那麼根據2D投影，就可以估計出R t；

空間3D點所在的線就确定了，根據多視圖（多視圖可以是運動圖像）可以重建3D。

如果場景已知，則可以把場景中的虛拟物體投影到2D圖像平面（DLT，隻要知道M即可）。或者根據世界坐标與錄影機坐标的相對關系，R,T，直接在Wc位置渲染3D圖形，這是AR的應用。

因為是離線的，是以可以用手工定位點的方法。

二 方法

1 Direct linear transformation (DLT) method 

2 A classical approach is "Roger Y. Tsai Algorithm".It is a 2-stage algorithm, calculating the pose (3D Orientation, and x-axis and y-axis translation) in first stage. In second

stage it computes the focal length, distortion coefficients and the z-axis translation. 

3 Zhengyou Zhang's "a flexible new technique for camera calibration" based on a planar chess board. It is based on constrains on homography.

4個空間：世界，錄影機(鏡頭所在位置)，圖像（原點到鏡頭中心位置距離f），圖像仿射空間（最終成像空間），這些空間原點在不同應用場合可以假設重合。

具體說來，前兩個空間是R t關系，6個變量搞定；第二到第三空間用相似三角形投影到平面上，如果沒有變形，光心在圖像中心上，到這裡為止其實就夠了，K就一個參數。

K是錄影機坐标系到圖像平面的轉換函數，叫内部參數(intrinsic parameter)。标定就是求上三角矩陣K的5個參數的過程。

用齊次（Homogeneous）坐标表示這個關系: （實際二維點是[u,v]'=[U/W , V/W]'）

       [U]    [a b -u0][-fxc/zc]   [-fa -fb -u0][xc/zc]        [xc/zc]

ub=[V]  =[0 c -v0][-fyc/zc ]= [0   -fc -v0][yc/zc]= K * [yc/zc]

       [W]   [0 0   1][  1    ]       [0   0     1]  [ 1   ]           [ 1   ]

世界坐标系到圖像坐标系的轉換有平移矢量t和旋轉R矩陣組成，R,t是外部參數(extrinsic parameters)，表示錄影機姿态，包含6個自由度,R t 各3個。Xc=[xc yc zc]'= R(Xw-t) ， ' 表示轉置。

把兩個合起來 

zc*ub=zc*[U V W]'= K*R(Xw-t)

把三維點Xw表示成齊次坐标，Xwb=[Xw,1]',則

       [U] 

ub=[V] = [K*R |-K*R*t][Xw 1]' =M*[Xw 1]' = M* Xwb

       [W]  ([K*R |-K*R*t]也可寫成K[R t']串聯的形式)

3×4矩陣M稱為投影矩陣，這樣三維點P到二維平面的投影就由上式的線性變換表達了。錄影機标定就是求K，K是與R,t一起耦合在M中的。

常用的方法是已知場景3D點和圖像二維點，求M。如Tsai grid方法。   

根據至少6個對應點，建立方程組:

alpha*[u v 1]' = M * [x y z 1]'

寫成2n×12矩陣 l*mij=0的形式(叫Direct Linear Transformation (DLT))：

[P1 | 0 | -u*P1]  [m11]

[0  |P1 | -v*P1]* [m12] = 0.   

[   ...                 ]   ...

[   ...                  ]   [m34]      (矢量Pi表示第i個三維點[Xi,Yi,Zi,1]')

方程解是l'*l的最小eigenvalue對應的eigenvector,|m|2=1，很多書上提到用SVD分解，l=UDV'，V就是特征矢量矩陣。

其實，l'*l=VDU'UDV'=VDDV',V'=inv(V),是以l'*l*V=VD^2，這就是特征矢量的求法。很明顯，l的singular vector是l'*l的特征矢量的平方根。

（PCA,K-L變換 和SVD之間的關系。邊肇琪《模式識别》提到人臉識别降維的例子時，講到如果樣本數小于特征維數，要用SVD提取特征向量。在這裡協方差矩陣為l'*l，是12×12的矩陣，l矩陣的短邊為9，而且計算M(Homography矩陣)時，點對要經過Normalize,完全可以參照人臉的例子。SVD總是和最小二乘，僞逆Ax=b-》A'*Ax=A'b-》x=inv(A'*A)*A'*b，聯系在一起。）

M=[K*R |-K*R*t]=[A|b],A = K*R,用QR分解（或SVD）将A分解為上三角矩陣K和機關正交陣R.

實際上M隻有11個參數（Rt6個，K5個）。如果把攝像頭畸變也考慮進去，mij變成16個值,則需要更多的對應點。

如果3D點在一個平面上,l is singular,那麼退化成9個參數。

已知場景3D點pattern拍一張照片，點分布在2到3個pattern平面上。

以上參考 milan sonka的Image processing一書。

實際上現在大家都是用張正友的方法，将網格點構成的pattern預先列印在一張紙上,對它拍攝不同方向照片（至少兩張）。

網上有matlab calibration包，OpenCV也有對應的實作，可以同時求錄影機的内部參數和外部參數(姿态估計)。

具體用到梯度下降方法，盡管R有3*3=9個數，但隻有3個自由度，用rodrigues formula轉為3維矢量，計算Jacobi更友善。

具體參考：

1 Rodrigues' Formula

<a href="http://www.cs.berkeley.edu/~ug/slide/pipeline/assignments/as5/rotation.html">http://www.cs.berkeley.edu/~ug/slide/pipeline/assignments/as5/rotation.html</a>

2 OpenCV函數說明

<a href="http://opencv.willowgarage.com/documentation/camera_calibration_and_3d_reconstruction.html">http://opencv.willowgarage.com/documentation/camera_calibration_and_3d_reconstruction.html</a>

3 Camera Calibration Toolbox for Matlab

<a href="http://www.vision.caltech.edu/bouguetj/calib_doc/">http://www.vision.caltech.edu/bouguetj/calib_doc/</a>

--------------------------------------------------

如果隻有平面4個點，如何計算構造DLT，求Homography矩陣[mij]？

1 左邊8×9矩陣A，m為9個數。rank 8, i.e. they have a one-dimensional null-space. 

The solution can be determined from the null-space of A.

svd分解 A=UDV',D為singular valuies，mij的解為V中最小的singular value對應列。

2 退化成一個平面到另一個平面的透視變換

u = (L1*X + L2*Y + L3)/(L7*X + L8*Y + 1) ;

v = (L4*X + L5*X + L6)/(L7*X + L8*Y + 1) ;

姿态估計:

基本思想是用3D點align到2D點上，使誤差最小。已知3D場景點，估計錄影機姿态，或者3D場景點關于錄影機坐标系的相對位置。

比如，marker可以構成世界坐标系，4個點坐标已知。它們相對于錄影機空間的坐标（Xa,Ya,Za）求出，xy軸可以求出，用右手法則z軸可以求出。

這樣世界坐标系的3根軸已知，就可以渲染3D場景了。（AR應用）

具體參見Hybrid camera pose estimation combining square fiducials localisation technique and orthogonal iteration algorithm

或 OI算法 都是類似梯度下降的疊代算法，可以利用前一幀的參數作為初值，實時估計R和T。

重建3D場景：

對于多個未經标定的錄影機，也可以根據在每個圖像中的2D點位置重建3D點

x1=P1X

x2=P2X

。。。。

Pi可以由已知場景點實時更新，xi可以根據運動估計更新，這樣可以計算X。

具體方法參見 Multiple View Geometry in Computer Vision 一書。