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); }
４点以上の参照点から最小二乗法を用いる場合
この場合、観測方程式は以下のようになります。
行列表現
パラメータ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); }
実行結果(例)
