基本矩阵的估计和图像矫正

0引言F矩阵描述的是从不同拍摄点拍摄同一场景的图像时两图中点的对应关系，它还提供有关位置和摄像机参数的信息。计算F矩阵最常用的是八点算法。要得到至少八个对应的图像点对。在这种算法下，我们可以很容易的使用线性方法来估计矩阵。但是八点算法对额外噪声十分敏感。为了克服这一缺点。哈特利（《计算机视觉中的多视图几何》的作者）引入一个坐标使图像标准化归一化。但是这种线性方法有一个缺点，那就是排除了秩等于2的约束。为了强加约束，这个F矩阵是利用奇异值分解重建获得的。这就造成了另一个问题：这个重建的过程会产生额外的噪声。为了减少噪声，哈特利又应用非线性优化设置七个变化参量来满足秩等于2的约束。关于那七个变量参数有很多相关的方法，我提出用一个标量的双四元数来描述矩阵。这些参数代表对一个图像缩放和旋转。也就是说，我们对原图像进行了旋转。通过使用这些旋转，我们创造了一个新的图像矫正方法。对于立体匹配来说，图像矫正是不可或缺的部分。它把任务由二维搜索变为一维搜索，减轻了任务量。流行的矫正方法大致可以分为两类：平面法和极性法。平面矫正法是对图像进行透视变换，这是容易理解且易于实现。然而，这种方法不能处理那些极点位于图像内部的图像对。极点矫正法能够处理这种情况，但他们需要根据极点的位置进行一些繁琐的预处理。F矩阵最麻烦的地方在于它的秩等于2，没有逆矩阵。为了计算F矩阵，人们提出了八点算法。这是一种易于实现的线性算法，但是对噪声太过敏感。人们又提出一种基于图像坐标归一化和非线性最优化的八点算法。还有很多其他成本函数最小化的方案。一些人采用另一种参数化的七自由度的方法来满足秩等于2的约束。对于线性求解，七点匹配是可以用于F矩阵的。图像矫正技术与F矩阵关系密切，所以密集点的匹配是不可或缺的。通过矫正，两图片中所有的极线会水平对齐，将二维搜索变为一维搜索，使得任务大大减轻。本质上来讲，图像矫正是一种旋转。另一方面，一些矫正方法已经适用于非标定图像对。非标定图像的矫正比标定图像更有用，因为准确的估计的相机参数是困难的，哪怕只有一个像素的错误就不利于准确匹配。1F矩阵的估计1.1对极几何(两视几何）对极几何描述的是场景中的三维点和一对图像点之间的关系。给定一个场景中三维点P被映射到I和I0两图像中，这三个点就构成一对极几何（图1）。点C和C0表示分别图像I和I0相机中心。三点P，C和C0组成平面，这就是所谓的极平面。每个图像平面与极平面相交，称为极线；一点投影在一个图像会有一个对应点位于其他图像平面的极线上（极线约束）。这个约束是指当我们有两个图像I和I0，和一点P在I上，如果我们要寻找在I0上对应点P0，仅沿极线搜索就可以了。假设I平面上的点TvuP,对图像和I0平面上的点TvuP''',都是是P的投影点，增加了齐次坐标的二维图像坐标，这样笛卡尔坐标下P和P旋转为射影坐标Tvum1,,和Tvum1,,'''。采用这种画法，极线约束表现为以下形式：0'FmmT（1）图1。对极几何Fig.1.Epipolargeometry在两幅图像之间，F矩阵将点m映射为对应的对极线，将对极点映射为0。不能提供对应点间的一一对应。F是一个3x3的矩阵。F矩阵的有秩等于2的约束，这是由于所有极线必须通过图像的极点。1.2八点算法的重新审视对于给定的图像对，常用八点算法来线性的估计F矩阵。对应的点ivu,i和''i,ivu，由线性系统得到方程：（2）其中（3）（4）其中f由F矩阵获得：（5）同时有约束条件1f233212211fff需要获得至少八个对，我们才能求解出（2）的估计矩阵。在刚才的过程中我们没有考虑“矩阵秩等于2”的约束。为了满足这一约束，选取满足约束和'FF的Frobenius范数最小的'F作为F矩阵。具体的方法如下：首先，对没有秩等于2约束条件下的矩阵进行奇异值分解分解（SVD）如下：（6）如果有在上述过程中无噪音，第三个奇异值，3应该是0。因此，以下矩阵'F可以取代作为约束矩阵的秩等于2：（7）这个重构矩阵'F可以作为F矩阵。这个过程实现了估计矩阵的秩为2的约束，但它引入附加噪声。为了解决这个问题，大多数工程在下一步采用非线性优化方法。由于对应点必须位于极线的约束，我们可以设计下面的几何成本函数：（8）mFmd'',表示极线上点P0和对应点P在图像的坐标系统下的距离的平方。1.3秩等于2的参数化方法我们所提出的方法的出发点是公式的形式（7）。这种形式表示的矩阵由两正交矩阵SO（3）和秩等于2的对角矩阵。如果不失一般性，我们可以指定对角矩阵为（1，s，0），其中12s。一般情况下，当一个正交矩阵的行列式是1，矩阵是一个旋转矩阵，可以用一个单位四元数描述。四元数3210,,,qqqqq其中123222120qqqq，旋转矩阵可以表示如下：（9）独立变量的数目是三；例如q1，q2和q3，如果公式（7）中的正交矩阵U和V是旋转矩阵，我们可以采用双四元（DQ）来参数化。不幸的是，不是所有的正交矩阵都是旋转矩阵。正交矩阵的行列式可以是1个或—1，行列式是—1的正交矩阵并不是旋转矩阵。然而，我们总是能把公式（7）中的两个正交矩阵旋变为两旋转矩阵，方法是通过改变对角矩阵中的另一个参数的标志s。有四种情况要考虑，根据组合两者的决定因素：应用上述程序，方程（7）成为两旋转矩阵和对角矩阵的乘积。旋转矩阵用式（9）那样的有三个参数的单位四元数描述，而对角矩阵只需要单参数。因此，估计矩阵一个七自由度的矩阵是：（10）1s（11）采用上述参数，几何成本函数可以改写为如下形式：（12）然后，上述的成本函数可以通过参数使其最小，同时保持秩等于2约束：（13）图2一个原始图像变换；图像I使用TR旋转，图像'I使用TR'旋转。Fig.2.Transformationofanoriginalimage;theimageIwithhomogeneouscoordinatesisrotatedbyTRandtheimage'IistransformedbyTR'.图3极线。它们被表示为变换后的图像和平面之间的交叉线，包括线x=y=0；平面包括线x=y=0被称为极平面。在这个图中，我们展示了两个极平面和两条极线。Fig.3.Epipolarlines.Theyarerepresentedasintersectionlinesbetweenthetransformedimagesandplanesthatincludethelinex=y=0;theplanesincludingthelinex=y=0areinterpretedasepipolarplanes.Inthisfigure,weshowtwoepipolarplanesandtwoelipolarlinesbydotedlines.2图片矫正由两个旋转矩阵和对角的参数化产生一种新的使所有极线平行的矫正方法。（14）令TzyxTnnnmRn,,，TzyxTnnnmRn'''''',,则式（14）改写如下：（15）上述方程表示如下比例关系：（16）更具体（17）符号xy代表二维矢量组成的前两个元素，即x和y作为的三位矢量的元素。方程（16）的几何解释如下：首先，图像I使用TR旋转而来，'I图像的旋转是由TR转化来的。我们称由TR和TR映射而来的图像平面作为变换后的图像。值得注意的是，三维空间的原点不是相机的中心，而是投影坐标系统的（0,0,0）。然后，变换后的图像进行垂直投影到xy平面。这个正交投影是使得式（17）中z分量消失。方程（16）表明，任何一对对应点在投影平面上位于同一穿过原点的直线上。换句话，由TR旋转后的I平面和一个过x=y=0的平面有一条交线，这条交线就是极线。由TR''变换后的平面'I和同一含x=y=0的平面相交于一线，这就是相对应的极线。那个包括线x=y=0平面被称为极平面。让我们考虑一个任意平面平行线x=y=0：z轴。我们称这种平面为参考平面。所有参考平面和极线之间的交线（其中包括线x=y=0），应平行于彼此。将旋转后图像平面上的一个点沿着线（极平面与平面z=常数之间的交线）投影到参考平面，即图4中箭头所示的投影。然后再投影点在参考平面上满足式（16）。因此，任何从旋转后平面再投影到参考平面的投影点满足方程（14）并且这个在参考平面上的再投影的图像和极线平行。通过改变在再投影过程中的z值，所有再投影点在参考平面上的构成矫正图像。矫正过程是可逆的：通过改变z的值，由一个矫正后的点可以找回它对应的唯一的原始点。图4。矫正图像：图像再投影到参考平面的矫正图像。Fig.4.Therectifiedimage:theimagere-projectedontothereferenceplanebecomestherectifiedimage.Inthere-projectionprocess,apixelonthetransformedimageplaneismappedalongtheepipolarplaneandperpendicularlytothez-axis.3实验3.1通过模拟数据评价矩阵为了评价DQ方法，同时为了比较其他现有的方法，我们利用一些人工数据集。我们选择了以下三种方法：极点参数法（EP），两极点参数法（BEP）和巴托丽参数法（BP）。3.1.1极点参数法（EP）（18）为了确保一个矩阵是奇异的，一列被描述为其他两列的线性组合。方程（18）有八参数；a~f最大绝对值被设置为1。3.1.2两极点参数法（BEP）（19）图5。相机配置。配置1（平行）和配置2（同一直线）。Fig.5.Cameraconfigurations.Configuration1(parallel)andConfiguration2(straight).图6。算法仿真结果。前两张是配置1，后两张是配置2。Fig.6.Simulationresultsofthealgorithms.Upper:Configuration1,Lower:Configuration2.在这个参数化，两极点位于T1,,T1,,''，a~d中绝对值最大的置1。3.1.3巴托丽参数法（BP）（20）这个参数是基于式（7）。巴托丽提出使用一个标量的六个旋转矩阵表示欧拉角。例如，R（x）代表一个欧拉角x绕x轴旋转矩阵。这四个参数化方法（EP，EP，BP和DQ），我们使用LevenburgMaquardt算法来最小化成本函数（8）。我们根据式（7）利用线性解的秩等于2约束求解非线性最小化的初始解。表1平均计算时间（毫秒）计算一个矩阵的结构1，N=30=2。Table1Themeancomputationaltimes(msec)forcalculatinganF-MatrixunderConfiguration1withN=30and=2.0.为了模拟，我们准备了一些人工数据集。我们随机产生N个三维空间的点并由图5中所示的两种配置产生一对480x640像素的图像。首先，我们用线性解来估计初始解与虚拟图像上的再投影点。然后，我们用四参数化方法得到最终的F矩阵。误差被定义为一个点的极线和其他图像上的对应点之间的距离。（21）通过改变噪声的水平和数量，我们可以比较这些方法的准确性。在本文中，我们重复100次，并利用平均残余误差。图6显示结果。上边是下配置1，下边的是在配置2。“SVD”指方程（7）的初始解。关于噪声的水平，我们将点数固定为N=30，并在图像中增加g高斯噪声来协调投影点。虽然误差是随着g上升增加，但是这四种方法（EP，EP，BP和DQ）在两种配置下有几乎相同的精度。有关用于计算F矩阵点数量，我们使用g=2。我们发现这四种参数化方法有很近似的精度。这些结果表明，我们的参数化方法是有竞争力的与其他现有方法。为了明确四种方法的性能，我们进行了T检验。根据结果证实，显著性差异不超过1%。表1显示了计算时间使用Athlon64X2运行