射影変換のパラメータを求める

4点なら

射影変換の式として以下の式を用いる場合、変換元と変換先の対応点が4組であればその座標値から一意にパラメータを求めることができます。

のとき、

//射影変換パラメータを求める。引数は、変換元と変換先の座標値
void CalcProjParam
(
    int    naOrig[4][2],
    int    naTran[4][2],
    double    *dA,
    double    *dB,
    double    *dC,
    double    *dD,
    double    *dE,
    double    *dF,
    double    *dG,
    double    *dH
)
{
    double **dATA, **dATA_I;
    int i;

    dATA = new double*[8];
    dATA_I = new double*[8];

    for ( i = 0; i < 8; i++ )
    {
        dATA[i] = new double[8];
        dATA_I[i] = new double[8];
    }

    dATA[0][0] =                 naOrig[0][0];
    dATA[0][1] =                 naOrig[0][1];
    dATA[0][2] =                            1;
    dATA[0][3] =                            0;
    dATA[0][4] =                            0;
    dATA[0][5] =                            0;
    dATA[0][6] = -1*naTran[0][0]*naOrig[0][0];
    dATA[0][7] = -1*naTran[0][0]*naOrig[0][1];
    dATA[1][0] =                 naOrig[1][0];
    dATA[1][1] =                 naOrig[1][1];
    dATA[1][2] =                            1;
    dATA[1][3] =                            0;
    dATA[1][4] =                            0;
    dATA[1][5] =                            0;
    dATA[1][6] = -1*naTran[1][0]*naOrig[1][0];
    dATA[1][7] = -1*naTran[1][0]*naOrig[1][1];
    dATA[2][0] =                 naOrig[2][0];
    dATA[2][1] =                 naOrig[2][1];
    dATA[2][2] =                            1;
    dATA[2][3] =                            0;
    dATA[2][4] =                            0;
    dATA[2][5] =                            0;
    dATA[2][6] = -1*naTran[2][0]*naOrig[2][0];
    dATA[2][7] = -1*naTran[2][0]*naOrig[2][1];
    dATA[3][0] =                 naOrig[3][0];
    dATA[3][1] =                 naOrig[3][1];
    dATA[3][2] =                            1;
    dATA[3][3] =                            0;
    dATA[3][4] =                            0;
    dATA[3][5] =                            0;
    dATA[3][6] = -1*naTran[3][0]*naOrig[3][0];
    dATA[3][7] = -1*naTran[3][0]*naOrig[3][1];
    dATA[4][0] =                            0;
    dATA[4][1] =                            0;
    dATA[4][2] =                            0;
    dATA[4][3] =                 naOrig[0][0];
    dATA[4][4] =                 naOrig[0][1];
    dATA[4][5] =                            1;
    dATA[4][6] = -1*naTran[0][1]*naOrig[0][0];
    dATA[4][7] = -1*naTran[0][1]*naOrig[0][1];
    dATA[5][0] =                            0;
    dATA[5][1] =                            0;
    dATA[5][2] =                            0;
    dATA[5][3] =                 naOrig[1][0];
    dATA[5][4] =                 naOrig[1][1];
    dATA[5][5] =                            1;
    dATA[5][6] = -1*naTran[1][1]*naOrig[1][0];
    dATA[5][7] = -1*naTran[1][1]*naOrig[1][1];
    dATA[6][0] =                            0;
    dATA[6][1] =                            0;
    dATA[6][2] =                            0;
    dATA[6][3] =                 naOrig[2][0];
    dATA[6][4] =                 naOrig[2][1];
    dATA[6][5] =                            1;
    dATA[6][6] = -1*naTran[2][1]*naOrig[2][0];
    dATA[6][7] = -1*naTran[2][1]*naOrig[2][1];
    dATA[7][0] =                            0;
    dATA[7][1] =                            0;
    dATA[7][2] =                            0;
    dATA[7][3] =                 naOrig[3][0];
    dATA[7][4] =                 naOrig[3][1];
    dATA[7][5] =                            1;
    dATA[7][6] = -1*naTran[3][1]*naOrig[3][0];
    dATA[7][7] = -1*naTran[3][1]*naOrig[3][1];

    // 逆行列
    matinv( 8, dATA, dATA_I );

    *dA = dATA_I[0][0]*naTran[0][0] + dATA_I[0][1]*naTran[1][0] + 
          dATA_I[0][2]*naTran[2][0] + dATA_I[0][3]*naTran[3][0] + 
          dATA_I[0][4]*naTran[0][1] + dATA_I[0][5]*naTran[1][1] + 
          dATA_I[0][6]*naTran[2][1] + dATA_I[0][7]*naTran[3][1];
    *dB = dATA_I[1][0]*naTran[0][0] + dATA_I[1][1]*naTran[1][0] + 
          dATA_I[1][2]*naTran[2][0] + dATA_I[1][3]*naTran[3][0] + 
          dATA_I[1][4]*naTran[0][1] + dATA_I[1][5]*naTran[1][1] + 
          dATA_I[1][6]*naTran[2][1] + dATA_I[1][7]*naTran[3][1];
    *dC = dATA_I[2][0]*naTran[0][0] + dATA_I[2][1]*naTran[1][0] + 
          dATA_I[2][2]*naTran[2][0] + dATA_I[2][3]*naTran[3][0] + 
          dATA_I[2][4]*naTran[0][1] + dATA_I[2][5]*naTran[1][1] + 
          dATA_I[2][6]*naTran[2][1] + dATA_I[2][7]*naTran[3][1];
    *dD = dATA_I[3][0]*naTran[0][0] + dATA_I[3][1]*naTran[1][0] + 
          dATA_I[3][2]*naTran[2][0] + dATA_I[3][3]*naTran[3][0] + 
          dATA_I[3][4]*naTran[0][1] + dATA_I[3][5]*naTran[1][1] + 
          dATA_I[3][6]*naTran[2][1] + dATA_I[3][7]*naTran[3][1];
    *dE = dATA_I[4][0]*naTran[0][0] + dATA_I[4][1]*naTran[1][0] + 
          dATA_I[4][2]*naTran[2][0] + dATA_I[4][3]*naTran[3][0] + 
          dATA_I[4][4]*naTran[0][1] + dATA_I[4][5]*naTran[1][1] + 
          dATA_I[4][6]*naTran[2][1] + dATA_I[4][7]*naTran[3][1];
    *dF = dATA_I[5][0]*naTran[0][0] + dATA_I[5][1]*naTran[1][0] + 
          dATA_I[5][2]*naTran[2][0] + dATA_I[5][3]*naTran[3][0] + 
          dATA_I[5][4]*naTran[0][1] + dATA_I[5][5]*naTran[1][1] + 
          dATA_I[5][6]*naTran[2][1] + dATA_I[5][7]*naTran[3][1];
    *dG = dATA_I[6][0]*naTran[0][0] + dATA_I[6][1]*naTran[1][0] + 
          dATA_I[6][2]*naTran[2][0] + dATA_I[6][3]*naTran[3][0] + 
          dATA_I[6][4]*naTran[0][1] + dATA_I[6][5]*naTran[1][1] + 
          dATA_I[6][6]*naTran[2][1] + dATA_I[6][7]*naTran[3][1];
    *dH = dATA_I[7][0]*naTran[0][0] + dATA_I[7][1]*naTran[1][0] + 
          dATA_I[7][2]*naTran[2][0] + dATA_I[7][3]*naTran[3][0] + 
          dATA_I[7][4]*naTran[0][1] + dATA_I[7][5]*naTran[1][1] + 
          dATA_I[7][6]*naTran[2][1] + dATA_I[7][7]*naTran[3][1];

    for ( i = 0; i < 8; i++ )
    {
        delete dATA[i];
        delete dATA_I[i];
    }
    delete dATA;
    delete dATA_I;
}
//射影変換を行う
void Projection
(
    int    nOrigX,
    int    nOrigY,
    double    dA,
    double    dB,
    double    dC,
    double    dD,
    double    dE,
    double    dF,
    double    dG,
    double    dH,
    double    *dTranX,
    double    *dTranY
)
{
    *dTranX = ((double)nOrigX*dA + (double)nOrigY*dB + dC) / 
        ((double)nOrigX*dG + (double)nOrigY*dH + 1);
    *dTranY = ((double)nOrigX*dD + (double)nOrigY*dE + dF) / 
        ((double)nOrigX*dG + (double)nOrigY*dH + 1);
}

4点以上の参照点から最小二乗法を用いる場合

この場合、観測方程式は以下のようになります。

行列表現

パラメータpの最小二乗解を求める式は以下のようになります。

ここで、


void Calc_ProjectionParam
(
    vector<Dpoint> &vOrig,
    vector<Dpoint> &vTrans,
    double aParam[8]
)
{
    UINT i, j;
    double dA[8][8];
    double dV[8];

    for ( i = 0; i < 8; i++ )
    {
        for ( j = 0; j < 8; j++ )
        {
            dA[i][j] = 0.0;
        }
        dV[i] = 0.0;
        aParam[i] = 0.0;
    }

    for ( i = 0; i < vOrig.size(); i++ )
    {
        // AT*A
        dA[0][0] += vOrig[i].x * vOrig[i].x;
        dA[3][3] += vOrig[i].x * vOrig[i].x;
        dA[0][1] += vOrig[i].x * vOrig[i].y;
        dA[1][0] += vOrig[i].x * vOrig[i].y;
        dA[3][4] += vOrig[i].x * vOrig[i].y;
        dA[4][3] += vOrig[i].x * vOrig[i].y;
        dA[1][1] += vOrig[i].y * vOrig[i].y;
        dA[4][4] += vOrig[i].y * vOrig[i].y;
        dA[0][2] += vOrig[i].x;
        dA[2][0] += vOrig[i].x;
        dA[3][5] += vOrig[i].x;
        dA[5][3] += vOrig[i].x;
        dA[1][2] += vOrig[i].y;
        dA[2][1] += vOrig[i].y;
        dA[4][5] += vOrig[i].y;
        dA[5][4] += vOrig[i].y;
        dA[2][2] += 1;
        dA[5][5] += 1;
        dA[0][6] += -1 * vTrans[i].x * vOrig[i].x * vOrig[i].x;
        dA[6][0] += -1 * vTrans[i].x * vOrig[i].x * vOrig[i].x;
        dA[0][7] += -1 * vTrans[i].x * vOrig[i].x * vOrig[i].y;
        dA[1][6] += -1 * vTrans[i].x * vOrig[i].x * vOrig[i].y;
        dA[6][1] += -1 * vTrans[i].x * vOrig[i].x * vOrig[i].y;
        dA[7][0] += -1 * vTrans[i].x * vOrig[i].x * vOrig[i].y;
        dA[1][7] += -1 * vTrans[i].x * vOrig[i].y * vOrig[i].y;
        dA[7][1] += -1 * vTrans[i].x * vOrig[i].y * vOrig[i].y;
        dA[2][6] += -1 * vTrans[i].x * vOrig[i].x;
        dA[6][2] += -1 * vTrans[i].x * vOrig[i].x;
        dA[2][7] += -1 * vTrans[i].x * vOrig[i].y;
        dA[7][2] += -1 * vTrans[i].x * vOrig[i].y;
        dA[3][6] += -1 * vTrans[i].y * vOrig[i].x * vOrig[i].x;
        dA[6][3] += -1 * vTrans[i].y * vOrig[i].x * vOrig[i].x;
        dA[3][7] += -1 * vTrans[i].y * vOrig[i].x * vOrig[i].y;
        dA[4][6] += -1 * vTrans[i].y * vOrig[i].x * vOrig[i].y;
        dA[6][4] += -1 * vTrans[i].y * vOrig[i].x * vOrig[i].y;
        dA[7][3] += -1 * vTrans[i].y * vOrig[i].x * vOrig[i].y;
        dA[4][7] += -1 * vTrans[i].y * vOrig[i].y * vOrig[i].y;
        dA[7][4] += -1 * vTrans[i].y * vOrig[i].y * vOrig[i].y;
        dA[5][6] += -1 * vTrans[i].y * vOrig[i].x;
        dA[6][5] += -1 * vTrans[i].y * vOrig[i].x;
        dA[5][7] += -1 * vTrans[i].y * vOrig[i].y;
        dA[7][5] += -1 * vTrans[i].y * vOrig[i].y;
        dA[6][6] += vTrans[i].x * vTrans[i].x * vOrig[i].x * vOrig[i].x + 
                    vTrans[i].y * vTrans[i].y * vOrig[i].x * vOrig[i].x;
        dA[6][7] += vTrans[i].x * vTrans[i].x * vOrig[i].x * vOrig[i].y + 
                    vTrans[i].y * vTrans[i].y * vOrig[i].x * vOrig[i].y;
        dA[7][6] += vTrans[i].x * vTrans[i].x * vOrig[i].x * vOrig[i].y + 
                    vTrans[i].y * vTrans[i].y * vOrig[i].x * vOrig[i].y;
        dA[7][7] += vTrans[i].x * vTrans[i].x * vOrig[i].y * vOrig[i].y + 
                    vTrans[i].y * vTrans[i].y * vOrig[i].y * vOrig[i].y;

        // AT*V
        dV[0] += vTrans[i].x * vOrig[i].x;
        dV[1] += vTrans[i].x * vOrig[i].y;
        dV[2] += vTrans[i].x;
        dV[3] += vTrans[i].y * vOrig[i].x;
        dV[4] += vTrans[i].y * vOrig[i].y;
        dV[5] += vTrans[i].y;
        dV[6] += -1 * ( vTrans[i].x * vTrans[i].x * vOrig[i].x + 
                        vTrans[i].y * vTrans[i].y * vOrig[i].x );
        dV[7] += -1 * ( vTrans[i].x * vTrans[i].x * vOrig[i].y + 
                        vTrans[i].y * vTrans[i].y * vOrig[i].y );
    }

    // AT*A-1 * AT*V
    matinv( 8, dA );

    for ( i = 0; i < 8; i++ )
    {
        for ( j = 0; j < 8; j++ )
        {
            aParam[i] += dA[i][j] * dV[j];
        }
    }
}

//射影変換を行う
void Projection
(
    Dpoint &sOrig,
    Dpoint &sTrans,
    double dParam[8]
)
{
    sTrans.x = (sOrig.x*dParam[0] + sOrig.y*dParam[1] + dParam[2]) / 
               (sOrig.x*dParam[6] + sOrig.y*dParam[7] + 1);
    sTrans.y = (sOrig.x*dParam[3] + sOrig.y*dParam[4] + dParam[5]) / 
               (sOrig.x*dParam[6] + sOrig.y*dParam[7] + 1);
}

実行結果(例)

正方形を茶色の点を参照点として変換した例

参考文献

アーカイブ