ガウス消去法の射影変換(座標変換の計算について)
Posted: 2014年2月26日(水) 14:32
現在、ガウス消去法による射影変換を勉強しています。
パラメータの計算までは求められているのですが、変換前の座標と変換パラメータを用いて変換後の座標を計算するところで詰まっています。
以下、現在書いているコードです。効率の良し悪しはとりあえず放っておいてます。
最後の変換の計算部分についてよくないところ、間違っているところを指摘頂ければと思います。
よろしくお願い致します。
パラメータの計算までは求められているのですが、変換前の座標と変換パラメータを用いて変換後の座標を計算するところで詰まっています。
以下、現在書いているコードです。効率の良し悪しはとりあえず放っておいてます。
最後の変換の計算部分についてよくないところ、間違っているところを指摘頂ければと思います。
よろしくお願い致します。
/* Gaussの単純消去法(部分ピボット選択) */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define X_SIZE 5000 /*画像の横サイズを定義*/
#define Y_SIZE 5000 /*画像の縦サイズを定義*/
#define MAXN 32
unsigned char screen[Y_SIZE][X_SIZE][3],img_work[Y_SIZE][X_SIZE][3],projection_work[Y_SIZE][X_SIZE][3];
int main(void)
{
/*--------------------【画像取得プログラム】用変数型--------------------*/
FILE *fp; /*ファイルポインタの宣言*/
unsigned char header[54]; /*unsigned:符号無しの整数型、char:文字データ、header[54]:0~53の1次配列*/
char input_file[128],gray_file[128],projection_file[128]; /*ファイル名の配列を定義*/
/*--------------------【射影変換プログラム】用変数型--------------------*/
int n=8, m, i, j, k, l, p, imax, x1, x2, x3, x0, y1, y2, y3, y0, u1, u2, u3, u0, v1, v2, v3, v0, U, V, z;
double b, t, t1, a[MAXN][MAXN + 1];
double x[4],y[4],u[4],v[4];
/*--------------------ここから【画像取得プログラム】--------------------*/
/*----------画像ファイルを読み込む----------*/
printf("処理を行う画像ファイル名:(拡張子:bmp)\n---");
scanf("%s",input_file); /*入力ファイル名を指定*/
printf("処理を行う画像ファイル:「 %s 」.\n\n",input_file);
fp=fopen(input_file,"rb"); /* 処理する画像をオープンする
Windowsビットマップ形式 X_SIZE*Y_SIZEピクセル,24ビットカラー、rb:バイナリデータの読み込みに使用 */
if((fp=fopen(input_file,"rb"))==NULL){ /*ファイルが開けない場合のエラー表示*/
printf("画像ファイルが正常に開けませんでした.\n");
exit(1);
}
fread(header,1,54,fp); /* ヘッダ(54バイト)を飛ばす */
fread(screen,1,X_SIZE*Y_SIZE*3,fp); /* 残りはデータ(最下行から順に入る) */
fclose(fp);
/*----------カラー画像をグレースケールに変換する----------*/
int X,Y,col;
for(X=0;X<X_SIZE;X++){
for(Y=0;Y<Y_SIZE;Y++){
for(col=0;col<=2;col++){
img_work[Y][X][col]=0.1145*screen[Y][X][0]
+0.5866*screen[Y][X][1]+0.2989*screen[Y][X][2];
}
}
}
/*----------グレースケール画像を保存する----------*/
printf("グレースケールにする画像ファイル名を入力してください.(拡張子:bmp)\n---");
scanf("%s",gray_file); /*グレースケールとするファイル名を指定*/
fp=fopen(gray_file,"wb"); /*処理した画像を書き込むファイルをオープンする*/
if((fp=fopen(gray_file,"wb"))==NULL){ /*ファイルが開けない場合のエラー表示*/
printf("画像ファイルが正常に開けませんでした.\n");
exit(1);
}
fwrite(header,1,54,fp); /* ヘッダ */
fwrite(img_work,1,X_SIZE*Y_SIZE*3,fp); /* データ */
fclose(fp);
printf("グレースケールの画像ファイル「 %s 」を保存しました.\n\n",gray_file);
/*--------------------ここまで【画像取得プログラム】--------------------*/
/*--------------------ここから【射影変換プログラム】--------------------*/
/* 行列の大きさをセット */
printf("ガウスの消去法により射影変換のパラメータを求めます。\n");
printf("↓ 8元連立方程式の結果 ↓\n\n");
//一般画像なら・・・ →〈1〉
//テンプレ画像なら・・・ →〈2〉
/*-----------↓ 〈1〉 ↓-----------*/
/*
printf("変換前後の座標をそれぞれ入力してください。(左上から時計回りの順)\n");
for(j = 0; j<4 ; j++){
printf("x%d=", j+1);
scanf("%lf", &x[j]);
printf("y%d=", j+1);
scanf("%lf", &y[j]);
}
for(m = 0; m<4 ; m++){
printf("u%d=", m+1);
scanf("%lf", &u[m]);
printf("v%d=", m+1);
scanf("%lf", &v[m]);
}
for(j = 0; j<4 ; j++){
printf("x%d=%lf\n",j,x[j]);
printf("y%d=%lf\n",j,y[j]);
}
for(j = 0; j<4 ; j++){
printf("u%d=%lf\n",j,u[j]);
printf("v%d=%lf\n",j,v[j]);
}
*/
/*-----------↑ 〈1〉 ↑-----------*/
/*-----------↓ 〈2〉 ↓-----------*/
x[0]=229;
y[0]=84;
x[1]=333;
y[1]=84;
x[2]=438;
y[2]=0;
x[3]=0;
y[3]=0;
u[0]=0;
v[0]=700;
u[1]=203;
v[1]=700;
u[2]=203;
v[2]=0;
u[3]=0;
v[3]=0;
a[0][0]=x[0];
a[0][1]=y[0];
a[0][2]=1;
a[0][3]=0;
a[0][4]=0;
a[0][5]=0;
a[0][6]=-x[0]*u[0];
a[0][7]=-y[0]*u[0];
a[0][8]=u[0];
a[1][0]=x[1];
a[1][1]=y[1];
a[1][2]=1;
a[1][3]=0;
a[1][4]=0;
a[1][5]=0;
a[1][6]=-x[1]*u[1];
a[1][7]=-y[1]*u[1];
a[1][8]=u[1];
a[2][0]=x[2];
a[2][1]=y[2];
a[2][2]=1;
a[2][3]=0;
a[2][4]=0;
a[2][5]=0;
a[2][6]=-x[2]*u[2];
a[2][7]=-y[2]*u[2];
a[2][8]=u[2];
a[3][0]=x[3];
a[3][1]=y[3];
a[3][2]=1;
a[3][3]=0;
a[3][4]=0;
a[3][5]=0;
a[3][6]=-x[3]*u[3];
a[3][7]=-y[3]*u[3];
a[3][8]=u[3];
a[4][0]=0;
a[4][1]=0;
a[4][2]=0;
a[4][3]=x[0];
a[4][4]=y[0];
a[4][5]=1;
a[4][6]=-x[0]*v[0];
a[4][7]=-y[0]*v[0];
a[4][8]=v[0];
a[5][0]=0;
a[5][1]=0;
a[5][2]=0;
a[5][3]=x[1];
a[5][4]=y[1];
a[5][5]=1;
a[5][6]=-x[1]*v[1];
a[5][7]=-y[1]*v[1];
a[5][8]=v[1];
a[6][0]=0;
a[6][1]=0;
a[6][2]=0;
a[6][3]=x[2];
a[6][4]=y[2];
a[6][5]=1;
a[6][6]=-x[2]*v[2];
a[6][7]=-y[2]*v[2];
a[6][8]=v[2];
a[7][0]=0;
a[7][1]=0;
a[7][2]=0;
a[7][3]=x[3];
a[7][4]=y[3];
a[7][5]=1;
a[7][6]=-x[3]*v[3];
a[7][7]=-y[3]*v[3];
a[7][8]=v[3];
/*-----------↑ 〈2〉 ↑-----------*/
for (k = 0; k < n; k++) {
/* ピポットを走査する */
t = 0; imax = k;
for (l = k; l < n; l++) {
if (fabs(a[l][k]) > t) {
imax = l;
t = a[l][k];
}
}
/* 見つかったピボット行をk行と入れ替える */
for (l = 0; l < n + 1; l++) {
t = a[k][l];
a[k][l] = a[imax][l];
a[imax][l] = t; // t ---> t1
}
/* 後は同じ */
b = a[k][k];
for (j = 0; j < n + 1; j++) a[k][j] /= b; // a /= b [a = a/b]
/* ピポット行よりも下の行について */
for (i = k + 1; i < n; i++) {
/* ピポットと同じ列を0に */
b = -a[i][k];
for (j = 0; j < n + 1; j++) a[i][j] += b * a[k][j];
}
}
/* 次に後置代入 */
for (k = n - 2; k >= 0; k--)
for (j = k + 1; j < n; j++)
a[k][n] -= a[k][j] * a[j][n];
/* 答を表示 */
printf("Gauss単純消去法(部分ピボット選択)による解\n");
for (i = 0; i < n; i++){
printf("--- x[%d]=%lf\n", i, a[i][n]);
}
/*--------------------画像の変換処理--------------------*/
// (ここまでは一応正しい出力がされていることを確認済みです。)
// (ここからの計算がうまくいきません。)
for(X=0;X<X_SIZE;X++){
for(Y=0;Y<Y_SIZE;Y++){
for(col=0;col<=2;col++){
U=((int)a[0][8]*X+(int)a[1][8]*Y+(int)a[2][8])/((int)a[6][8]*X+(int)a[7][8]*Y+1);
V=((int)a[3][8]*X+(int)a[4][8]*Y+(int)a[5][8])/((int)a[6][8]*X+(int)a[7][8]*Y+1);
projection_work[Y][X][col]=img_work[V][U][col];
}
}
}
printf("射影変換した画像ファイル名を入力してください.(拡張子:bmp)\n---");
scanf("%s",projection_file); /*グレースケールとするファイル名を指定*/
fp=fopen(projection_file,"wb"); /*処理した画像を書き込むファイルをオープンする*/
if((fp=fopen(projection_file,"wb"))==NULL){ /*ファイルが開けない場合のエラー表示*/
printf("画像ファイルが正常に開けませんでした.\n");
exit(1);
}
fwrite(header,1,54,fp); /* ヘッダ */
fwrite(projection_work,1,X_SIZE*Y_SIZE*3,fp); /* データ */
fclose(fp);
printf("射影変換した画像ファイル「 %s 」を保存しました.\n\n",projection_file);
/*--------------------ここまで【射影変換プログラム】--------------------*/
return 0;
}