MVG读书笔记——单应矩阵估计这件小事（三）

坐标变换下的不变性

实际应用中，不同程序对于坐标原点的规定可能不同，坐标系的尺度也可能不同。于是很自然会有一个问题，在这些情况下我们通过两幅图像的对应点求解出的单应矩阵还是同一个么？即单应矩阵有没有不变性？

用更加正式点语言表述如下：

假设一幅图像中的点x在另外一套坐标系中表示方法为$\widetilde x$，它的在另一幅图像中的对应点x’在另一套坐标系中表示方法为$\widetilde x'$。其中

$$\widetilde x_i=Tx_i,\widetilde x_i'=T'x_i'$$
考虑 $x'=Hx$则有 $\widetilde x'=T'HT^{-1}\widetilde x$

于是可以得到$\widetilde H = T'HT^{-1}$

于是，从新的坐标$\widetilde x\leftrightarrow\widetilde x'$我们可以得到$\widetilde H$。进而得到$H=T'^{-1}\widetilde HT$

那么问题来了，对于任意的$T,T'$我们是否总能求出同样的H呢？

DLT算法的不变性

考虑一系列对应点集合$x_i\leftrightarrow x'_i$，使用DLT算法求解出它们的单应矩阵为$H$。再考虑对$\widetilde x\leftrightarrow\widetilde x'$，其中$\widetilde x_i = Tx_i,\widetilde x'_i=T'x_i'$。令$\widetilde H=T'HT^{-1}$。上一节的的问题就变成了：

DLT算法是否能从$\widetilde x\leftrightarrow\widetilde x'$解出上述的H？

如上节一样，定义$\widetilde\epsilon_i=\widetilde x'_i\times \widetilde H\widetilde x_i$,可以计算出

$$\begin{align}\widetilde\epsilon_i&=\widetilde x'_i\times \widetilde H\widetilde x_i=T'x'_i \times (T'HT^{-1})Tx_i\\&=T'x_i'\times T'Hx_i=T'^*(x_i'\times Hx_i)\\&=T'^*\epsilon_i\end{align}$$

其中$T^*=det(T)T^{-T}$（存在定理$(Mx)\times(My)=M^*(x\times y)$）

对于相似变换$T=\begin{bmatrix}sR&t\\0^T&1\end{bmatrix}$，可以求出$T'^*=s\begin{bmatrix}R&0\\-t^TR&s\end{bmatrix}$.
将$T'^*$带入$\epsilon_i$的前两项有。

$$\widetilde A_i\widetilde h=(\widetilde \epsilon_{i1} \widetilde\epsilon_{i2})^T=sR(\epsilon_{i1} \epsilon_{i2})^T=sRA_ih$$
R为正交矩阵，所以R不改变向量的范数，于是 $||\widetilde A \widetilde h||=s||Ah||$

用代数误差来表示就是

$$d_{alg}(\widetilde x'_i,\widetilde H\widetilde x_i)=sd_{alg}(x'_i,Hx_i)$$

由此，我们发现对于同样的代数误差（相差一个尺度），可以推导出一组对应的H和$\widetilde H$。可能由此我们会认为最小化代数误差解出的H和$\widetilde H$满足$\widetilde H=T'HT^{-1}$。然而事实并非如此。

因为在DLT算法求解时为了使得方程有唯一解我们添加了一项约束$||H||=1$。于是对于相同的代数误差，在||H||=1和$||\widetilde H||=1$下解出的两者并不满足前面所说的等式。

事实上，它们不满足简单的对应关系，就是说

$$minimize \sum_id_{alg}(x'_i,Hx_i)^2 其中||H||=1\\ \Leftrightarrow minimize \sum_id{alg}(\widetilde x_i',\widetilde H\widetilde x_i)^2 其中||H||=1\\ \neq minimize \sum_id{alg}(\widetilde x_i',\widetilde H\widetilde x_i)^2 其中||\widetilde H||=1$$

因此我们说坐标的变换会导致不同的解。只有当||Ah||不随坐标变换而变换时我们可以找出H与$\widetilde H$之间的关系。如何做到这一点，又是另一个问题。接下来我们将探讨几何误差的不变性。

几何误差的不变性

使用最小化几何误差的方法来寻找H，对于相似变换具有不变性。如前面提到的，考虑一系列对应点集合$x_i\leftrightarrow x'_i$，使用DLT算法求解出它们的单应矩阵为$H$。再考虑对$\widetilde x\leftrightarrow\widetilde x'$，其中$\widetilde x_i = Tx_i,\widetilde x'_i=T'x_i'$。令$\widetilde H=T'HT^{-1}$。假设T，T’代表${IP^2}$中的欧氏变换。可以验证

$$d(\widetilde x',\widetilde H\widetilde x)=d(T'x',T'HT^{-1}Tx)=d(T'x',T'Hx)=d(x',Hx)$$

这是显然的，因为欧式变换不改变长度（距离）。因此它们最小化的是同一个东西，求出来的H也是双射的，即最小化几何误差对于欧式变换具有不变性。

对于相似变换，几何误差之间显然相差一个相似系数。计算出的$H$与$\widetilde H$可以如欧式变换的情况中一样简单的一一对应起来。所以几何误差具有相似变换不变性。

标准化变换

正如开始所讲，DLT算法得出的单应取决于坐标框架，它对于相似变换不具有不变性。也就是说某个坐标框架中计算出的H要好于另外一个框架中计算出的H。为了减少甚至消除这一点的影响，我们在执行DLT算法之前需要进行标准化。

对一幅图像中所有点进行标准化，使得它们的中心移到坐标原点，同时使得它们与坐标原点的平均距离等于$\sqrt 2$（这个数值是为了计算的方便）。即所有点在以原点为圆心，$r = \sqrt 2$的圆两侧。

这样的做法使得各种坐标框架上的点被变换到了同一套坐标系下，再进行DLT算法，就使得它有了相似变换不变性。

由此，我们对DLT算法进行改进：
1）对x进行标准化，变换记作T
2）对x’进行标准化，T’
3）DLT算法
4）反标准化。$H=T'^{-1}\widetilde HT$

于是得到了H。

附录:OpenCV里面的求解算法

可以看到，OpenCV使用的单应矩阵的求解算法就是上面所说的方法

---

int runKernel( InputArray _m1, InputArray _m2, OutputArray _model ) const
{
    Mat m1 = _m1.getMat(), m2 = _m2.getMat();
    int i, count = m1.checkVector(2);
    const Point2f* M = m1.ptr<Point2f>();
    const Point2f* m = m2.ptr<Point2f>();

    double LtL[9][9], W[9][1], V[9][9];
    Mat _LtL( 9, 9, CV_64F, &LtL[0][0] );
    Mat matW( 9, 1, CV_64F, W );
    Mat matV( 9, 9, CV_64F, V );
    Mat _H0( 3, 3, CV_64F, V[8] );
    Mat _Htemp( 3, 3, CV_64F, V[7] );
    Point2d cM(0,0), cm(0,0), sM(0,0), sm(0,0);

    for( i = 0; i < count; i++ )
    {
        cm.x += m[i].x; cm.y += m[i].y;
        cM.x += M[i].x; cM.y += M[i].y;
    }

    cm.x /= count;
    cm.y /= count;
    cM.x /= count;
    cM.y /= count;

    for( i = 0; i < count; i++ )
    {
        sm.x += fabs(m[i].x - cm.x);
        sm.y += fabs(m[i].y - cm.y);
        sM.x += fabs(M[i].x - cM.x);
        sM.y += fabs(M[i].y - cM.y);
    }

    if( fabs(sm.x) < DBL_EPSILON || fabs(sm.y) < DBL_EPSILON ||
        fabs(sM.x) < DBL_EPSILON || fabs(sM.y) < DBL_EPSILON )
        return 0;
    sm.x = count/sm.x; sm.y = count/sm.y;
    sM.x = count/sM.x; sM.y = count/sM.y;

    double invHnorm[9] = { 1./sm.x, 0, cm.x, 0, 1./sm.y, cm.y, 0, 0, 1 };
    double Hnorm2[9] = { sM.x, 0, -cM.x*sM.x, 0, sM.y, -cM.y*sM.y, 0, 0, 1 };
    Mat _invHnorm( 3, 3, CV_64FC1, invHnorm );
    Mat _Hnorm2( 3, 3, CV_64FC1, Hnorm2 );

    _LtL.setTo(Scalar::all(0));
    for( i = 0; i < count; i++ )
    {
        double x = (m[i].x - cm.x)*sm.x, y = (m[i].y - cm.y)*sm.y;
        double X = (M[i].x - cM.x)*sM.x, Y = (M[i].y - cM.y)*sM.y;
        double Lx[] = { X, Y, 1, 0, 0, 0, -x*X, -x*Y, -x };
        double Ly[] = { 0, 0, 0, X, Y, 1, -y*X, -y*Y, -y };
        int j, k;
        for( j = 0; j < 9; j++ )
            for( k = j; k < 9; k++ )
                LtL[j][k] += Lx[j]*Lx[k] + Ly[j]*Ly[k];
    }
    completeSymm( _LtL );

    eigen( _LtL, matW, matV );
    _Htemp = _invHnorm*_H0;
    _H0 = _Htemp*_Hnorm2;
    _H0.convertTo(_model, _H0.type(), 1./_H0.at<double>(2,2) );

    return 1;
}