一 作用
建立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
http://www.cs.berkeley.edu/~ug/slide/pipeline/assignments/as5/rotation.html
2 OpenCV函數說明
http://opencv.willowgarage.com/documentation/camera_calibration_and_3d_reconstruction.html
3 Camera Calibration Toolbox for Matlab
http://www.vision.caltech.edu/bouguetj/calib_doc/
--------------------------------------------------
如果只有平面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 一書。
http://blog.sciencenet.cn/home.php?mod=space&uid=465130&do=blog&id=365366