相机参数矩阵计算（内外矩阵）以及光束法平差

单应性矩阵由配对的内点计算出来

u,v是原始图片左边，对应得到变换后的图片坐标x,y,其中

变换矩阵可以拆成4部分，表示线性变换，比如scaling，shearing和ratotion。用于平移，产生透视变换。

重写之前的变换公式可以得到：

OpenCV 代码实现是getPerspectiveTransform()和findHomography()

/* Calculates coefficients of perspectivetransformation

 *which maps (xi,yi) to (ui,vi), (i=1,2,3,4):

 *

 *     c00*xi + c01*yi + c02

 * ui= ---------------------

 *     c20*xi + c21*yi + c22

 *

 *      c10*xi + c11*yi + c12

 * vi= ---------------------

 *     c20*xi + c21*yi + c22

 *

 *Coefficients are calculated by solving linear system:

 * /x0 y0  1 0  0  0 -x0*u0 -y0*u0   /c00 /u0

 * |x1 y1  1 0  0  0 -x1*u1 -y1*u1 | |c01| |u1|

 * |x2 y2  1 0  0  0 -x2*u2 -y2*u2 | |c02| |u2|

 * |x3 y3  1 0  0  0 -x3*u3 -y3*u3 |.|c10|=|u3|,

 *|  0 0  0 x0 y0  1 -x0*v0 -y0*v0 | |c11| |v0|

 *|  0 0  0 x1 y1  1 -x1*v1 -y1*v1 | |c12| |v1|

 *|  0 0  0 x2 y2  1 -x2*v2 -y2*v2 | |c20| |v2|

 *  0 0  0 x3 y3  1 -x3*v3 -y3*v3 /  c21/ v3/

 *

 *where:

 *   cij- matrix coefficients, c22 = 1

 */

cv::Mat cv::getPerspectiveTransform( constPoint2f src[], const Point2f dst[] )

{

}

更完善的函数是findHomography()

pairwise_matches_ contains the matching relation.  BestOf2NearestMatcher.BestOf2NearestMatcher::match() which is derived from FeaturesMatcher which hasone operator(), the () below will invoke BestOf2NearestMatcher::match() whichwill calculate H matrix. 单应性矩阵计算从匹配的内点

matches_info.H = findHomography(src_points,dst_points, matches_info.inliers_mask, RANSAC);

都用于计算单应矩阵，即解一个线性方程组。由于单应矩阵有8个未知数（3*3，其中第9个数为1），所以至少需要4个点（每个点-x,y,提供2个约束方程）。

发现getPerspectiveTransform用的是SVD分解，getPerspectiveTransform只会拿前4个点去计算，getPerspectiveTransform()比较简单粗暴。

findHomography则会拿一堆点(>=4)去计算(其是不断从一堆点中重复拿出4个点去计算出一个结果，再采用一些优化算法RANSAC/LMEDS去筛选出最优解)。

如果参数被设置为0，那么这个函数使用所有的点和一个简单的最小二乘算法来计算最初的单应性估计，但是，如果不是所有的点对都完全符合透视变换，那么这个初始的估计会很差，在这种情况下，你可以使用两个robust算法中的一个。RANSAC和LMeDS , 使用坐标点对生成了很多不同的随机组合子集（每四对一组），使用这些子集和一个简单的最小二乘法来估计变换矩阵，然后计算出单应性的质量，最好的子集被用来产生初始单应性的估计和掩码。

RANSAC方法几乎可以处理任何异常，但是需要一个阈值， LMeDS 方法不需要任何阈值，但是只有在inliers大于50%时才能计算正确，最后，如果没有outliers和噪音非常小，则可以使用默认的方法。

相机参数估计(R,T,K 矩阵)

背景知识部分请参考相机针孔成像模型(三大坐标系变换，外参数矩阵(R*T),内参数矩阵K)

和设备相关基础知识

光束法平差：相机参数计算的基础

//摄像测量学的基本原理Bookmark

http://www.cnblogs.com/libing64/archive/2012/02/12/2878731.html

//3Dreconstruction

http://archimede.bibl.ulaval.ca/archimede/fichiers/25229/ch06.html     

// BundleAdjustment 光束法平差详解中文版本

http://m.blog.csdn.net/blog/junshen1314/48860951

// BundleAdjustment 光束法平差详解英文版本

http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/ZISSERMAN/bundle/bundle.html

//光束法平差， 多核计算优化

http://www.doc88.com/p-7854008874656.html

要通过用摄像机、照相机拍摄的图像得到目标实体的空间信息，首先需要建立描述目标实体与图像之间映射关系的成像模型。为了使图像与实际景象具有相似性，特别是与人眼的成像机理相似，绝大部分成像设备都采用如图 1所示针孔模型，或称为中心透视投影模型来描述成像关系的。图 1中，目标点 P 的空间坐标为（X，Y，Z），像机光心为Oc，像平面为 Π，P 对应像点的图像坐标为p(x，y)。中心透视投影成像关系是：所有的成像光线都通过光心，目标点、像点与光心共线。从而得到摄像测量的最基本数学关系式，即共线方程。其中，（Cx，Cy）为主点坐标，是像机光轴与成像平面的交点；（Fx，Fy）为等效焦距，是实际焦距分别与单个像元宽度和高度的比值；R为旋转矩阵，T 为平移向量，用于描述像机的位置和姿态；λ为比例常数。

摄像测量的任务主要是从二维图像获取目标的三维信息，其中最基本的是目标的位置、姿态和结构参数，而更多的速度、加速度、角速度、变形参数等则在位置、姿态、结构参数的基础上进行综合分析得到。测量中最常用的确定目标三维信息的方式是利用两台或多台像机同时成像，进行交会测量。其基本原理如图2所示，两台像机成像所得到的像点 p、p′分别与对应光心 Oc、O′c确定了一条光线，这两条光线必然都通过目标点 P，因而它们的交点就是目标点的空间位置。

基本求解算法具体解算时，典型的代表算法是光束法平差（bundle adjustment）。光束法平差的基本原理是对被测物体从不同角度、位置拍摄多幅图像，则物体上所有特征点的物点、像点、光心都应落在相应的光线上，即满足成像的共线关系。平差的过程就是通过优化使这些共线关系的总体偏差为最小。 如下图。

光束法平差的先决条件是，找到不同相机拍摄到图片的对应的点（相同物理点的图片投影），这里就涉及到特征提取和配对，以及配对角点的选取。譬如。

A 3D point  is projected atpoint and  in correspondingleft and right image planes of a stereo pair. Suppose and are corresponding projections matrices used in theprojection, two rays can be reconstructed from the corresponding projectedpoints and to the initial 3D point . In theory, the two lines should intersect perfectly inspace but due to noise and other errors it is not the case. The followingsections describe different approaches that have been proposed for 3Dreconstruction.

Bundle adjustment is the problemof refining a visual reconstruction to produce jointly optimal 3D structuresand viewing parameters (camera pose and/or calibration) estimates. This methodsearches for the optimal 3D point and the optimal projection matricessimultaneously, by minimizing the distance between observed points in the imageplanes and projected 3D scene points from the projection matrices as shown inequation (6.16).

(6.16)

where is the observed image point, the unknown projection matrix, the 3D reconstructed point and the Euclidian distance. To resolve this problem, manyalternative optimization algorithms such as Levenberg-Marquardt’s Algorithm orDog Leg Algorithm can be used [53].

4

假设我们有一个3D空间中的一个点Xj，他可以被多个位于不同角度的摄像机看到，设第i个摄像机看到的第j个点坐标为xij ,由3D点像二维相机平面转换的矩阵为Pi .我们现在考虑这样的问题，给出一系列的的坐标xij，找到相机转换矩阵Pi使得PiXj=xij

and

OpenCV stitching模块的相机参数和单应性

http://www.360doc.com/content/14/0410/14/10724725_367760906.shtml

http://www.360doc.com/content/14/0410/14/10724725_367760675.shtml

http://blog.csdn.net/chenfeiyang2009/article/details/7212091

摄像机的透视投影模型（即针孔成像模型）：
       设OXYZ为世界坐标系，uv为以像素为单位的图像坐标系。如果物点P在世界坐标系下的坐标为（X,Y,Z）,对应的图像点p在图像坐标系的坐标为（u，v），可以有：
                 

在该式中，fu、fv、u0、v0只与摄像机内部参数有关，故称矩阵M1为内参数矩阵。
     其中fu ＝ f/dX ,fv = f/dY ,分别称为u轴和v轴上的归一化焦距；f是相机的焦距，dX和dY分别表示传感器u轴和v轴上单位像素的尺寸大小。
    u0和v0则表示的是光学中心，即摄像机光轴与图像平面的交点，通常位于图像中心处，故其值常取分辨率的一半。
现以NiKon D700相机为例进行求解其内参数矩阵：
就算大家身边没有这款相机也无所谓，可以在网上百度一下，很方便的就知道其一些参数——        焦距 f = 35mm  最高分辨率：4256×2832    传感器尺寸：36.0×23.9 mm
根据以上定义可以有：
u0 = 4256/2 =2128   v0 = 2832/2 = 1416  dx = 36.0/4256   dy =23.9/2832
fu = f/dx = 4137.8   fv = f/dy = 4147.3
该模型中其他一些参数的含义解释：
  为投影深度，其几何意义是目标点P在摄像机坐标系中坐标的Z分量，一般情况下，该值令为1.
R，t则构成了摄像机的外参数矩阵，分别称为旋转矩阵和平移矩阵。R是一个3X3的矩阵，t则是一个3X1的矩阵。R可以使用四元素法表示也可以使用欧拉角表示，各有各的好处。
摄像机的仿真，可以简单的认为是对内外参数的仿真。只要理解了以上模型中各个量的含义，就可以很容易的模拟一个摄像机。

焦距计算（K矩阵）

K矩阵有[焦距, ppx(图片的宽)，ppy（图片的高），aspect（宽高比）]等计算出来的。 这里焦距需要计算。

焦距计算从单应性矩阵里面计算出来，一个H矩阵计算出两个焦距。单应性矩阵由配对的内点计算出来。

HomographyBasedEstimator::estimate() =>

estimateFocal(features, pairwise_matches,focals);

R矩阵计算

根据匹配的内点对两个K矩阵和内点对的H矩阵计算出R矩阵

HomographyBasedEstimator::estimate() =>

span_tree.walkBreadthFirst(span_tree_centers[0],CalcRotation(num_images, pairwise_matches, cameras)); =>

struct CalcRotation

{

   CalcRotation(int _num_images, const std::vector<MatchesInfo>&_pairwise_matches, std::vector<CameraParams> &_cameras)

        :num_images(_num_images), pairwise_matches(&_pairwise_matches[0]),cameras(&_cameras[0]) {}

    void operator()(const GraphEdge &edge)

    {

        intpair_idx = edge.from * num_images + edge.to;

……………………………………………

                            //yueqiang:calculate the R matrix based on K and H

        Mat R =K_from.inv() * pairwise_matches[pair_idx].H.inv() * K_to;

       cameras[edge.to].R = cameras[edge.from].R * R;

}