ちょっと気合いを入れて、と。
 画像処理の講義の内容が難しくなり(というか、数学の知識と画像処理の知識とをうまく結びつけられていないのであるが)、すでに2回分の講義内容をできないままになっている。
 さらに、来週までにレポート、とのこと。

 というわけで、今週は『画像処理強化週間(仮称)』といきます。

 できた人に聞けば(ヒントでも良いし、アドバイスだけでも良いから)いいんだろうけど、つまらない性格のため、よう聞かないんだよね。それに、答えを見てレポートが書けても、その場しのぎにしかならないし、今までも必死で講義について行っていただけで、自分のものになっていない、ということ等から、きちっと理解してレポートを作ろう!と思ったんだよね(だからといって、できるとは限らないが)。

 ということで、S氏とK氏の協力を仰ぎ、なんとかプログラムを完成したいな、と思っているのであった。



●rotate-nn.c
これがもとのプログラム。
●rotate-nn2.c
a = cos (th); b = sin (th); c = 127.0;
d = -sin (th); e = cos (th); f = 127.0;
30,31行目のc,fに127.0変更。しかし、回転する原点の位置を変えただけ。
※ここで、薄々c,fを変えるだけではないことに気づく(回転の中心がやっぱり原点だから)。
●rotate-nn3.c
for (y = -127; y < y_size1/2; y ++) {
for (x = -127; x < x_size1/2; x ++) {
34,35行目を変更。だめ。
●rotate-nn4.c
x1 = x-(x_size1/2) - c; y1 = y+(y_size1/2) - f;
38行目を変更。やっぱりだめ。
●rotate-nn5.c
a = cos (th); b = sin (th); c = x_size1/2;
d = -sin (th); e = cos (th); f = -y_size1/2;
30,31行目を変更。だめ。
●rotate-nn6.c
a = cos (th); b = sin (th); c = -x_size1/2;
d = -sin (th); e = cos (th); f = y_size1/2;
30,31行目を変更。だめ。
※ここで、c,f,x1,y1等を変更してもだめなことに気づく。

 アフィン変換(平行移動、回転、拡大縮小等の処理に対応)について調べる。ここで、c,fは平行移動する、ということなんだと理解する。
 変換前の画像の座標系を(x、y)、変換後の座標系を(x*,y*)とすると、座標変換は
 x*=g(x,y)
 y*=h(x,y)
で表される。特に、
 x*=ax+by+c
 y*=dx+ey+f
 ここで、回転の場合は、
 a=cosθ
 b=sinθ
 c=0
 d=-sinθ
 e=cosθ
 f=0
 回転におけるθは座標原点を中心とする回転角を表す。(参考文献2pp114,pp115)

 アフィン変換等の幾何学的変換により出力画像を求めるには、画素の位置を格子状配列に配列し直す再配列(resampling)が必要(画素の座標が変換されたときに変換後の画素の間に隙間が空くという問題が生じるため)。これをさけるために、出力画像の各画素に対応する入力画像の座標値を逆変換により求め、そこから画素を持ってくるようにする。
 このとき、逆変換により計算された入力画像の座標は必ずしも整数値とは限らない(格子状配列の画素の座標は整数値をとる)ため、求めたい点の濃度値をその周辺の画素の濃度値から内挿(interpolation)して近似的に求める必要がある。
 [1]最近隣内挿法(nearest neighbor interpolation)
 内挿したい点に最も近い格子点の画素値を求める点の画素値とする。処理アルゴリズムが簡単であり、オリジナルの画素値のデータが破壊されないという利点がある。ただし、位置誤差が最大1/2画素生じ、変換後の画像の境界線等に細かいギザギザ(ジャギー)が生じることがある。

 逆行列をする理由はなんとなくわかったけど、結局まだ画像の中心で回転させる方法はわからず。

 ここで、S氏より「できた!!!!」との声が!!!しかし、なぜできたのかはわからない模様。そこで、プログラムrotate-nn7.cをじっくりと観察。
 
●rotate-nn7.c
a = cos (th); b = sin (th); c = 127.0;/*変更した部分*/
d = -sin (th); e = cos (th); f = 127.0;/*変更した部分*/
det = a * e - b * d;

for (y = 0; y < y_size1; y ++) {
for (x = 0; x < x_size1; x ++) {
double x1, y1, x2, y2;
int ix, iy;
x1 = x - c; y1 = y - f;/*ここで、元画像を左上に移動(x方向に-127、y方向に-127)*/
x2 = ( e * x1 - b * y1) / det;
y2 = (-d * x1 + a * y1) / det;
ix = (int)(x2 + 0.5)+127; iy = (int)(y2 + 0.5)+127;/*変更した部分。ここで回転された画像を基準にして画像を戻す(x方向に127、y方向に127)*/

 上記の注釈のように変更すると、なんと!うまくいきました!!!!!
a0014712_1833434.jpg

 なんとまぁ、灯台もと暗し、て感じ?しかし、これは本来目指していた解法ではないため、再び本来の解法を目指す。

 ここで、参考文献4のpp151を読むと・・・
 ・・・例えば、(x0,y0)を中心にした回転は、・・・画像をまず(-x0,-y0)だけ移動させて(x0,y0)点を原点に持ってきたあとにθ°回転させて、さらに(x0,y0)だけ移動させればよいのです。・・・・
と書いているではないか!!!!なんだ、rotate-nn7.cのままじゃん。
 しかし、あきらめずに考える。
●rotate-nn8.c
for (y = -y_size1/2; y < y_size1/2; y ++) {
for (x = -x_size1/2; x < x_size1/2; x ++) {
y軸の始点が-y辺/2、終点がy辺/2・・・としてみた。しかし、画像を左上に移動しただけではだめ。
●rotate-nn9.c
for (y = -y_size1/2; y < y_size1/2; y ++) {
for (x = -x_size1/2; x < x_size1/2; x ++) {
double x1, y1, x2, y2;
int ix, iy;
x1 = x - c; y1 = y - f;
x2 = ( e * x1 - b * y1) / det;
y2 = (-d * x1 + a * y1) / det;
ix = (int)(x2 + 0.5); iy = (int)(y2 + 0.5);
if (-x_size1/2 <= ix && ix < x_size1/2 &&
-y_size1/2 <= iy && iy < y_size1/2) {
image2[y+y_size1/2][x+x_size1/2] = image1[iy+y_size1/2][ix+x_size1/2];
} else {
image2[y+y_size1/2][x+x_size1/2] = 0;
 これで、画像を左上に移動して回転し、最後に戻す。これでOK!!!やった!!!
a0014712_639635.jpg
 

 ということで、本題の共一次内挿法にいくことにする。最近隣内挿法では、ジャギーが目立つことがあるので、目立たない方法ということで、共一次内挿法を使ってプログラムを作ることになった。

 [2]共一次内挿法(bi-lenear interpolation)
 周囲の4画素の濃度値をf(u,v)からの距離の比率により線形補完して濃度値を求める。i,jをu,vの値を超えない最大の整数とし、
p=u-i、q=v-j
とすれば、次式により濃度値が求められる。(略)
 この方法では変換後の画像に平滑化の効果があるが、画像がぼやけた感じになることがある。

 講義である程度のプログラムが提示された。
double x,p,q,f0,f1,f2,f3,g0,g1;

p=x2-(int)x2;
q=y2-(int)y2;

f0=image1[iy][ix];
f1=image2[iy][ix+1];
f2=image1[iy+1][ix];
f3=image2[iy+1][ix+1];
g0=
g1=
image2[y][x]=q*g1+(1.0-q)*g0+0.5;


 ここで考えなければならないことは、上記の式をプログラムのどこに入れるのか、また、どこを削除するのか。そして、追加する部分。

 プログラムをじっくりと見ていくと、
x2 = ( e * x1 - b * y1) / det;
y2 = (-d * x1 + a * y1) / det;
 ここで、回転行列。そして、
ix = (int)(x2 + 0.5)+127; iy = (int)(y2 + 0.5)+127;
これがいわゆる最近隣内挿法の式(回転行列によって出された座標を整数値(小数点以下四捨五入)にしている)。つまり、この式の変わりに、上記の式を入力することになると考えられる。
 与えられた式と、講義ノート、テキスト(5・23式)等を参考にして、
内挿したい点
f(u,v)=(1-q){f[i,j]・(1-p)+f[i+1,j]・p}+q{f[i,j+1]・(1-p)+f[i+1,j+1]・p}
    =(1-q){f0・(1-p)+f1・p}+q{f2(1-p)+f3・p}
 つまり、
g0=p*f1+(1-p)*f0;
g1=p*f3+(1-p)*f2;
となる。早速、これらの式を代入してみる。
 
●rotate-bl.c
> gcc -o rotate-bl rotate-bl.c -lm
rotate-bl.c: In function `main':
rotate-bl.c:15: conflicting types for `x'
rotate-bl.c:13: previous declaration of `x'
rotate-bl.c:56: array subscript is not an integer
rotate-bl.c:61: array subscript is not an integer
rotate-bl.c:63: array subscript is not an integer
 エラー多発!!!
 原因はintとdoubleの両方にxがあったためであった。
 そのほか、以下のように変更。
f0=image1[iy][ix];
f1=image1[iy][ix+1];/*image2は出力画像なのですべて1にした*/
f2=image1[iy+1][ix];
f3=image1[iy+1][ix+1];

g0=p*f1+(1-p)*f0;
g1=p*f3+(1-p)*f2;

image2[y][x]=q*g1+(1.0-q)*g0+0.5;/*整数値の座標を出す*/
/*image2は出力なので、以下の部分を省略*/
/* ix = (int)(x2 + 0.5); iy = (int)(y2 + 0.5); 最近隣内挿法*/
/* if (0 <= ix && ix < x_size1 &&
0 <= iy && iy < y_size1) {
image2[y][x] = image1[iy][ix];
} else {
image2[y][x] = 0;
}*/
 すると、一応画像が表示されたものの、灰色。
a0014712_5541057.jpg

<現在、朝の6時42分・・・> 
●rotate-nnbl.c
/* ix = (int)(x2 + 0.5); iy = (int)(y2 + 0.5);*/
if (-x_size1/2 <= ix && ix < x_size1/2 &&
-y_size1/2 <= iy && iy < y_size1/2) {
/*ここから共一次内挿法のメインプログラム*/
p=x2-(int)x2;/*小数点以下を出す*/
q=y2-(int)y2;/*小数点以下を出す*/

f0=image1[iy][ix];
f1=image1[iy][ix+1];
f2=image1[iy+1][ix];
f3=image1[iy+1][ix+1];

g0=p*f1+(1-p)*f0;
g1=p*f3+(1-p)*f2;

image2[y+y_size1/2][x+x_size1/2]=q*g1+(1.0-q)*g0+0.5;/*整数値の座標を出す*/

/*image2[y+y_size1/2][x+x_size1/2] = image1[iy+y_size1/2][ix+x_size1/2];*/
} else {
image2[y+y_size1/2][x+x_size1/2] = 0;
 エラーはでなかったが、真っ黒。がく。
●rotate-nnbl2.c
image1[iy][ix]=(int)(q*g1+(1.0-q)*g0+0.5);/*整数値の座標を出す*/
if (-x_size1/2 <= x && x < x_size1/2 &&
-y_size1/2 <= y && y < y_size1/2) {
image2[y+y_size1/2][x+x_size1/2] = image1[iy+y_size1/2][ix+x_size1/2];
} else {
image2[y+y_size1/2][x+x_size1/2] = 0;
 ここで、Segmentation fault (core dumped)がでるようになった。
●rotate-nnbl3.c
 試行錯誤中・・・
●rotate-nnbl4.c
f0=image1[iy][ix];
f1=image1[iy][ix+1];
f2=image1[iy+1][ix];
f3=image1[iy+1][ix+1];

g0=p*f1+(1-p)*f0;
g1=p*f3+(1-p)*f2;

if (-x_size1/2 <= x && x < x_size1/2 &&
-y_size1/2 <= y && y < y_size1/2) {
image2[y+y_size1/2][x+x_size1/2]=(int)(q*g1+(1.0-q)*g0+0.5);/*整数値の座標を出す*/
} else {
image2[y+y_size1/2][x+x_size1/2] = 0;
 これでエラーはなくなり、画像も一応出力された。しかし、画像の端が白く、半分欠けている。
a0014712_725821.jpg

●rotate-nnbl5.c
f0=image1[iy+y_size1/2][ix+x_size1/2];
f1=image1[iy+y_size1/2][ix+x_size1/2+1];
f2=image1[iy+y_size1/2+1][ix+x_size1/2];
f3=image1[iy+y_size1/2+1][ix+x_size1/2+1];

g0=p*f1+(1-p)*f0;
g1=p*f3+(1-p)*f2;

if (-x_size1/2 <= x && x < x_size1/2 &&
-y_size1/2 <= y && y < y_size1/2) {
image2[y+y_size1/2][x+x_size1/2]=(int)(q*g1+(1.0-q)*g0+0.5);/*整数値の座標を出す*/
 右下へ画像が移動しすぎていると考え、f0-f3までの配列にx,y辺/2を入れた。
 すると、エラーはなくなり、位置もOK.しかし、画像がぼけるはずなのに、さらにジャギーが目立つように・・・
a0014712_744840.jpg

 小手先のごまかしではうまくいかない。参考文献4のpp156のプログラムを参考に、ひとつひとつ確認していくことにした。
●rotate-nnblx.c
/* 回転変換のプログラム rotate-nnblx.c */
/* 共一次内挿法(Bi-lenear) */

#include
#include
#include
#include"mypgm.h"

main( )
{
double a, b, c, d, e, f, th;
double det;
int x, y;
double pi;
double p,q,f0,f1,f2,f3,g0,g1;

pi = 4.0 * atan(1.0);

printf("*** 回転変換のプログラム ***\n");
printf("*** 共一次内挿法 ***\n");

printf("\n★ 元画像を読み込みます.\n");
load_image_data( );

printf ("回転角(度)= "); scanf ("%lf", &th);
printf ("\n★ %10.2lf度回転します。\n", th);

th = pi / 180.0 * th;

/* アフィン変換行列 */
a = cos (th); b = sin (th); c = 0.0;
d = -sin (th); e = cos (th); f = 0.0;
det = a * e - b * d;

for (y = -y_size1/2; y < y_size1/2; y ++) {
for (x = -x_size1/2; x < x_size1/2; x ++) {
double x1, y1, x2, y2;
int ix, iy;
x1 = x - c; y1 = y - f;
x2 = ( e * x1 - b * y1) / det;
y2 = (-d * x1 + a * y1) / det;

if (y2>0) iy=y2; else iy=y2-1;・・・・・・・・・・・・・・追加
if (x2>0) ix=x2; else ix=x2-1;・・・・・・・・・・・・・・追加
/* ix = (int)(x2 + 0.5); iy = (int)(y2 + 0.5);*/・・・・・・以下、ix,iyに置き換えて引数を整理していく。
p=x2-ix;/*小数点以下を出す*/
q=y2-iy;/*小数点以下を出す*/
if (-x_size1/2 <= ix && ix < x_size1/2 &&
-y_size1/2 <= iy && iy < y_size1/2) {

f0=image1[iy+y_size1/2][ix+x_size1/2];
f1=image1[iy+y_size1/2][ix+x_size1/2+1];
f2=image1[iy+y_size1/2+1][ix+x_size1/2];
f3=image1[iy+y_size1/2+1][ix+x_size1/2+1];

g0=p*f1+(1-p)*f0;
g1=p*f3+(1-p)*f2;

image2[y+y_size1/2][x+x_size1/2]=(int)(q*g1+(1.0-q)*g0+0.5);/*整数値の座標を出す*/
} else {
image2[y+y_size1/2][x+x_size1/2] = 0;
}
}
}

x_size2 = x_size1; y_size2 = y_size1;

printf("\n★ 回転した画像を保存します.\n");
save_image_data( ); /* 画像を保存する */
return 0;
}


 12月11日、午前8時22分、ついに完成・・・・(たぶん)。
a0014712_8405559.jpg

 とりあえず、プログラム的には完成したようである。ただ、まだプログラムを完全に理解できたわけではない。参考になるようなプログラムがなければ完成はしなかっただろう。さぁ、あとはレポートだ。レポートを作りながら、もう一度プログラムを見直そうっと。(・x・)

◎参考文献・Web
1.数学ナビゲータ
2.画像処理工学;末松良一・山田宏尚共著;コロナ社;2000;pp114-pp117
3.行列と一次変換・数列;細野真宏著;中経出版;1993;pp55-pp57
4.C言語で学ぶ実践画像処理;八木伸行等;オーム社;1992;pp150-pp159
[PR]
by viewtleaf | 2004-12-10 16:49 | 画像処理
<< 徹夜 情報応用演習 サーバについて >>



毎日起こった出来事を記入し、「振り返り」に活用したいと思います。写真は我が子の作品。
by viewtleaf
S M T W T F S
1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30