Homography


ホモグラフィ変換用の行列(転写関数)を算出する方法に関して、覚書として以下を残しておく。
下記画像は、乱数関数によりランダムに生成された出力領域に、ホモグラフィ転写した結果。

各画像をクリックすると、実際の変換画像が表示される。

Homography
Homography001 Homography002 Homography003 Homography004 Homography005 Homography006 Homography007 Homography008 Homography009 Homography010 Homography011 Homography012 Homography013 Homography014 Homography015 Homography016 Other areas



//////////////////////////////////////////////////////////////////////
// Homography
までは、
乱数より生成された座標点で、出力領域になり得る座標群であるかを確認する部分。
以降が、変換用の3x3の行列を生成する手順。(異なる4点の対応条件が必要)

原画上の四隅の点が、SA, SB, SC, SD の時、転写先の対応点は、DA, DB, DC, DD。
    SA => DA
    SB => DB
    SC => DC
    SD => DD
HomographyTrans() は、逆行列を用いた画像転写処理ルーチン


typedef struct POINT { float x, y; } POINT;

int Homography( PIXEL **s, int slx, int sly, PIXEL **d, PIXEL **c, int lx, int ly ){
    srand( (unsigned int)time(NULL) );
    for(int e=0; e<16; ) {////   (POINT){ float x, y; }
        POINT A = (POINT){ rand() % MAX_SIZE, rand() % MAX_SIZE };
        POINT B = (POINT){ rand() % MAX_SIZE, rand() % MAX_SIZE };
        POINT C = (POINT){ rand() % MAX_SIZE, rand() % MAX_SIZE };
        POINT D = (POINT){ rand() % MAX_SIZE, rand() % MAX_SIZE };
        if( (A.x-C.x)*(B.y-C.y) == (B.x-C.x)*(A.y-C.y) ) continue;
        if( (B.x-D.x)*(C.y-D.y) == (C.x-D.x)*(B.y-D.y) ) continue;
        if( (C.x-A.x)*(D.y-A.y) == (D.x-A.x)*(C.y-A.y) ) continue;
        if( (D.x-B.x)*(A.y-B.y) == (A.x-B.x)*(D.y-B.y) ) continue;
        //////////////////////////////////////////////////////////////////////
        ////       (PTSEL){   x    y  id   norm                     flag}
        PTSEL rawp[4] = { { A.x, A.y, 'A', sqrt( A.x*A.x + A.y*A.y ), 0 },
                          { B.x, B.y, 'B', sqrt( B.x*B.x + B.y*B.y ), 0 },
                          { C.x, C.y, 'C', sqrt( C.x*C.x + C.y*C.y ), 0 },
                          { D.x, D.y, 'D', sqrt( D.x*D.x + D.y*D.y ), 0 },
                        };  if( sort_by_order(rawp) ) continue;;
        e++;
        //////////////////////////////////////////////////////////////////////
        // Homography
        POINT  SA = (POINT){ slx,   0 }, DA = (POINT){ rawp[0].x, rawp[0].y };
        POINT  SB = (POINT){ slx, sly }, DB = (POINT){ rawp[1].x, rawp[1].y };
        POINT  SC = (POINT){   0, sly }, DC = (POINT){ rawp[2].x, rawp[2].y };
        POINT  SD = (POINT){   0,   0 }, DD = (POINT){ rawp[3].x, rawp[3].y };
        double tranp[8]  = { DA.x, DA.y, DB.x, DB.y, DC.x, DC.y, DD.x, DD.y };
        double abcdefgh1[3*3]={0,0,0,0,0,0,0,0,1}; //Coefficient of 3x3 matrix
        double hg[8][8] = {
                 { SA.x, SA.y,  1,     0,    0,  0,  -SA.x*DA.x, -SA.y*DA.x },
                 {    0,    0,  0,  SA.x, SA.y,  1,  -SA.x*DA.y, -SA.y*DA.y },
                 { SB.x, SB.y,  1,     0,    0,  0,  -SB.x*DB.x, -SB.y*DB.x },
                 {    0,    0,  0,  SB.x, SB.y,  1,  -SB.x*DB.y, -SB.y*DB.y },
                 { SC.x, SC.y,  1,     0,    0,  0,  -SC.x*DC.x, -SC.y*DC.x },
                 {    0,    0,  0,  SC.x, SC.y,  1,  -SC.x*DC.y, -SC.y*DC.y },
                 { SD.x, SD.y,  1,     0,    0,  0,  -SD.x*DD.x, -SD.y*DD.x },
                 {    0,    0,  0,  SD.x, SD.y,  1,  -SD.x*DD.y, -SD.y*DD.y },
               }, invhg[8][8], hkl[3][3], invhkl[3][3];
        //////////////////////////////////////////////////////////////////////
        int x, y, z,  b=255,g=255,r=255; // Clear canvas for drawing images
        for(y=0;y<ly;y++)for(x=0;x<lx;x++)c[y][x]=(PIXEL){b,g,r,0};
        //////////////////////////////////////////////////////////////////////
        inverseMatrix( 8, invhg, hg );   // Calculate the inverse matrix
        // { {a,b,c}, {d,e,f}, {g,h,1} } abcdefgh1[3*3]={0,0,0, 0,0,0, 0,0,1};
        for(y=0;y<8;y++)for(x=0;x<8;x++) abcdefgh1[y] += invhg[y][x]*tranp[x];
        for(z=0,y=0;y<3;y++)for(x=0;x<3;x++) hkl[y][x] = abcdefgh1[z++];
        //////////////////////////////////////////////////////////////////////
        inverseMatrix( 3, invhkl, hkl ); // Calculate the inverse matrix
        HomographyTrans( c, lx, ly, s, slx, sly, 3, invhkl );
        /////////////////////////////////////////////////////
        char fname[88]; snprintf(fname,sizeof(fname)-1,"Homography%3.3d.bmp",e);
        saveArray2bmp32( fname, c, lx, ly );  // Save the image as a BMP file
    }// Consider the adjugate matrix because it's processed by a computer.
    return 0;
}
 ※ inverseMatrix(...) コンピュータによる逆行列の算出は、余因子行列から逆行列を求める
   旧式のnotePCでも、再帰呼び出しで8x8の逆行列の算出は 30ms程度なので処理時間は許容の範囲
   ただし少し工夫をした改良版にすれば、処理時間が 数ms以下になる。
   処理方法の改善により、13x13の逆行列の算出が30s程度なので、実用的には改良版の利用が前提



追加

変形した画像から、正方形の画像を再変換する例

画像(Homography011.bmp)に対して、以下の様な転写座標を設定して変換を行った結果
ただし、座標点は左下を原点とした時の位置を示している。

   SA(1055,1372)  =>  DA( 200,1200)
   SB(1109,1165)  =>  DB(1200,1200)
   SC( 333, 826)  =>  DC(1200, 200)
   SD( 175, 925)  =>  DD( 200, 200)

再変換された画像
Reconverted_thumbnail

Reference example coding
int inverseMatrix( int k, double inverse[k][k], double matrix[k][k] ){ int x, y, below = k-1; long double detA, adjugate[k][k]; for(y=0;y<k;y++)for(x=0;x<k;x++) inverse[y][x] = (x==y) ? 1 : 0; if( k<=1 || (detA=determinant( k, matrix ))==0 ) return -1; for(y=0;y<k;y++)for(x=0;x<k;x++){ double cofactor[below][below]; for(int m=0, n=0, u=0;u<k;u++)for(int v=0;v<k;v++) { if( u==y || v==x ) continue; cofactor[m][n++] = matrix[u][v]; // Extraction of cofactors if( below <= n ){ m++, n=0; } } long double dtw = determinant( below, cofactor ); adjugate[x][y] = (y+x)&1 ? -dtw : dtw;//[y][x]->[x][y] Symmetric matrix } for(y=0;y<k;y++)for(x=0;x<k;x++) inverse[y][x]=adjugate[y][x]/detA; return 0; } ※ IEEE double の指数部(11ビット)の表現力不足が原因で、高次の行列式の計算は失敗する ※ 行列式の記憶領域double から long double (多分:64ビットから80ビット)に変更

ランダム生成された 13x13行列,64x64行列,92x92行列,128x128行列の、
逆行列を算出した例   (ホモグラフィ変換とは無関係な別の処理用)

・ 4-core, 8-Thread 2.30GHz Note PC 13x13 改良版  約30秒
・16-core,32-Thread 3.49GHz Desktop PC 13x13 改良版  約 6秒
・ 4-core, 8-Thread 2.30GHz Note PC 64x64 最適化版  約1.5秒
・16-core,32-Thread 3.49GHz Desktop PC 64x64 最適化版  約300ms
 92x92行列 現行(double型変数)の上限, stack 8Mにアップ 
・16-core,32-Thread 3.49GHz Desktop PC 92x92 最適化版  約990ms
 128x1282行列 double -> long double, stack 12Mにアップ 
・16-core,32-Thread 3.49GHz Desktop PC 128x128 最適化版  約11秒
下10桁版


スタック利用が原因で動作しないので、使わない様に変更し高速化させたが、
512x512の逆行列が約20分では使えない。(スマートな行列式算出法が必要)
750x750の逆行列が約3時間30分では問題外

Matrix Size: 512
START_try: 512x512 14:25:36.416269
E N D_try: 512x512 14:47:13.905295

START_try: 512x512 14:47:14.321943
E N D_try: 512x512 15:08:36.838259

START_try: 512x512 15:08:37.255890
E N D_try: 512x512 15:30:12.540603

START_try: 512x512 15:30:12.957733
E N D_try: 512x512 15:51:41.640397

START_try: 512x512 15:51:42.063563
E N D_try: 512x512 16:13:20.993235

 再計算結果:64x64,128x128,256x256,512x512


※ long double は環境依存、スタックサイズの変更は editbin.exe を使用
  GPU等の補助計算機は不使用。CPU内の演算器,マルチスレッド機構のみを利用
  256x256行列の逆行列は求め無いとは思うが M256_16C32T_L5_Optimization
  又 256x256逆行列の算出時間(主に行列式算出)は、最適化版でも10分は必要
  逆行列の要素が小数点以下10桁の精度なら、積の精度は小数点以下6桁程度
  実用性で言うなら、Octave かな。