7 動画像処理 intx, y;/* 制御変数 * / intvx,vy,step;/* 速度べクトル , べクトルを求める点の間隔 * / printf("frame No. 1 を読み込みます . Vn"); / * imagel - - > image2 copy—image1_t0—image2( ) ; / * frame No. 2 ー > imagel ☆ / load—image—data( ) ; printf("frame NO. 2 を読み込みます . yn"); printf("= = 動きべクトルの検出 = main( ) draw-vector( x, y, x-vctr[y][x] , y-vctr[y][x]); if ( ( y % step ) = = 0 & & ( x % step ) = for ( x = 0 ; x く x—size2; x + + ) for ( y = 0 ; y く y_size2; y 十 + ) image2[y][x] = 0 ; for (x = 0 ; x く x—sizeX x 十 + ) for ( y = 0 ; y < y—sizeX y + + ) だ計測結果を画像にする * / x vctr[y][x] = vx; y-vctr [y][x] = vy ; calc_motion_vector( x, y, &vx, &vy ) ; if ( ( x % step ) = = 0 ) { for(x=O; x<x_sizel; x + + ) { if ( ( y % step ) = printf("line %dYn",y); for()= 0 ; y<y—sizel; y + + ) { printf(" 画像 1 から画像 2 への速度べクトルを求めています . yn"); だ処理開始 y-vctr[y][x] = 0 ; x-vctr[y][x] = 0 ; for()= 0 ; x<x—sizel; x + + ) { for(y=0; y<y_sizel; y 十十 ) だ速度場の画像の初期化 * / } while ( step く 2 Ⅱ step > 20 ) ; scanf("%d" &step); printf(" 速度べクトルを計算する画素間の距離 ( 2 ~ 20 ) : " ) ; do { 163

6 立体認識 解答例 / * 画像を簡易 3 次元表示するプログラム draw3D. c * / #include<stdio. h> #include<stdlib. h> #include<math. h> #include"mypgm. h" 139 void draw a_straight_line( int xl, int yl, int x2, int y2, int brightness ) だ image2[y][x] 上に ( xl , yl ) , ( x2 , y2 ) 間の直線を描画する . だ線の明るさは brightness ( 0 ~ 255 ) で指定する . double distance, step, t; だ t : 内分点発生用制御変数 * / intx,y; だ描く点の座標 * / distance = sqrt( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) ) ; / * 2 点間の距離 * / step= 1.0 / ( 1.5 * distance); / * 発生させる内分点の総数 * / t= 0.0 ; while ()< 1.0 ) { x = (int)( い x2 + ( 1.0 - t 戸 xl ) ; y = (int)( t * y2 + ( 1.0 - t 戸 yl ) ; if( 0 く = x & & x く x_size2 & & 0 <= y & & y < y_size2 ) image2 [y][x] = (unsigned char)brightness; t = t 十 step; void draw_3D_graph( ) int y, x; int plot_step; / * ループ変数 * / だ線の描画間隔 [ 画素 ] 2 ~ 10 程度 float plot-ratio; / * 階調方向の倍率 0.0 ~ 1.0 ( = 255 画素戸 / int plot—height; だ階調方向の高さ = 256.0 * plot_ratio * / / * グラフの想定横幅 ( 画素 ) 響 x_size2 = x_sizel 十 y_sizel / 2 ; だグラフの想定縦幅 ( 画素 ) * / y—size2 = y—sizel 十 256 ; if ( x—size2 > MAX IMAGESIZE Ⅱ y_size2 > MAX IMAGESIZE ) { printf ( " 画像が想定最大サイズを超えています . vn"); exit(l); } else { / * 描画パラメータ ( plot_step, plot_ratio ) の指定 * / printf(" 線の間隔 ( 画素 , 2 ~ 10 程度 ) : " ) ; scanf("%d" &plot_step);

6 立体認識 では後ろの曲線を消すときにも用いるため , 描く直線の明るさを引数 brightness で指定し , 0 を指定すると黒い線を描くようにしている . 関数 draw-3D-graph() が前述の方法でグラフを描く関数である . は変数 plot-step によって , 描く階調変化曲線の間隔を指定し , 変数 plot—ratio で , 階調方向の倍率を指定している . ただし plot-ratio は 0.0 ~ 1.0 の間の実数で指定し , 1.0 にすると階調方向の最大値を 256 画素にする . 図 6.10 ( a ) に示す , 描画ソフトウェアを用いて作成したテスト画像を原画 像とし , 本プログラムを実行した結果得られた画像を同図 ( b ) , ( c ) に示し , 各図で用いたパラメータの値を画像の下に示す . 同図から , 非常に簡単な方 法ではあるが , 画像の階調変化が直観的に分かり易く表現されていることが わかる . 141 (a) 原画像 ( 256X256 画素 ) (c) step = 3 , ratio = 0.3 図 6. 10 画像の 3 次元グラフ表示の例 (b) step = 5 , ratio = 0.5

160 7 動画像処理 例題 7-4 : プロックマッチングによるオプティカルフローの決定 図 7.9 に示したプロックマッチングに基づいて , 連続する 2 枚のフレー ム間のオプティカルフローを決定するプログラムを作成せよ . このとき , 速度べクトルを求める点を , 第 1 フレームの横方向 , 縦方向とも一定画 素間隔で設定せよ . 解答例 : だ速度べクトルを求めるプログラム motion. c 響 #include く stdio. h> #include く stdlib. h> #include<math. h> #include"mypgm. h" void copy—image 1—t0—image2( ) /*imagel を image2 へコピーする * / int x, y; だループ変数 * / X Size2 = X sizel; y_size2 = y—sizel; for(y=O; y<y_sizel; y 十十 ) for()= 0 ; x<x_sizel; x 十 + ) image2[y][x] = imagel[y][x] ; void draw_vector( int xl, int yl, int vx, int vy ) / * image2[y][x] に ( xl , yl ) の速度べクトル ( vx , vy ) を描画する . * / だ t : 内分点発生用制御変数 * / double distance, step, t; int x2, y2; だ終 * / だ描く点の座標 * / int x, y; X2 = xl 十 vx; y2 =yl 十 vy; distance = sqrt( ( x2- xl 戸 ( x2- xl ) + ( y2- yl 戸 ( y2- yl ) ) ; step = 1.0 / ( 1.5 * distance ) ; t = 0.0 ; while ( t く 1.0 ) { x = (int)( に x2 + ( 1.0 - t 戸 xl ) ; y = (int)( い y2 + ( 1.0 - t 戸 yl ) ; if( x > = 0 & & x < x_size2 & & y > = 0 & & y く y—size2 ) image2 [y][x] = (unsigned char)MAX BRIGHTNESS; t = t 十 step;

5 パターン認識 double calc_distance( int xl, int yl, int x2, int y2 ) / * (x1,y1),(x2,y2) 間の距離を返す関数 return sqrt( ( x2- xl 当 x2- xl ) + ( y2- yl 戸 ( y2- yl ) ) ; void draw a_straight—line( int xl, int yl, int x2, int y2 ) 層 image2[y][x] 上に ( xl , yl ) , ( x2 , y2 ) 間の直線を発生させる . * / / * t: 内分点発生用制御変数ツ double distance, step, t; int x, y; / * 描く点の座標 * / distance = calc—distance( xl, yl, x2, y2 ) ; だ 2 点間の距離 * / step = 1.0 / ( 1.5 * distance); / * 発生させる内分点の総数 * / t = 0.0 ; while ()< 1.0 ) { x = (int)( い x2 + ( 1.0 - t 戸 xl ) ; y = (int)( い y2 + ( 1.0 - t 戸 yl ) ; if( x > = 0 & & x < x_size2 & & y > = 0 & & y く y_size2 ) image2[y][x] = DARK PIXEL; だ暗い点 * / t = t 十 step; int random int( int n ) / * 0 以上 n - 1 以下の整数の乱数を返す * / points[num—pnts - 1 Ⅱ 0 ] = x; num_pnts 十十 ; if ( num—pnts く MAX PNTS - 2 ) { if ( imagel[y][x] = = MAX BRIGHTNESS ) { for(x = 0 ; x く x—sizel; x + + ) { for ( y = 0 ; y < y—sizel; y + + ) { num—pnts = 0 ; int x, y; だループ変数 * / 層 imagel[y][x] 中の白画素のデータを points[][] に代入 * / void obtain—points data( ) return (int)( (double)rand( ) / ( (double)RAND MAX + 1.0 戸 n ) ; 109

66 3 2 値画像処理 例題 3-6 : Houqh 変換 あらかじめいくつかの直線を含む 2 値原画像 ( 白 : 255 , 黒 : 0 ) を作成 した後 , その画像を読み込んで H 。 ugh 変換を施し , の平面上の分布を 画像として保存するプログラムを作成せよ . ただし , の平面を表す 2 次元配列の要素数は日方向は 360 , 〃方向は 720 とし , 各要素の値は 255 を上限値にせよ . 解答例 : だ Hough 変換のプログラム Hough. c * / #include<stdio. h> #inclu de<stdlib. h> #include<math. h> #include"mypgm. h" #define PI 3.141592653589 / * 円周率 #define MAX THETA 360 だ〃軸のサイズ ( 0.5 [ deg ト 1 画素 ) * / 720 / * p 軸のサイズ #define MAX RHO void draw a_straight_line( int xl, int yl, int x2, int y2 ) / * image2[y][x] 上の 2 点 ( xl , yl ) , ( x2 , y2 ) 間の直線上の配列要素☆ / / * の値をインクリメントする . * / だ t : 内分点発生用制御変数☆ / double distance, step, t; intx,y; だ描く点の座標 * / distance = sqrt( (x2-x1)*(x2-x1) + ( y2- yl 戸 ( y2- yl ) ) ; だ 2 点間の距離 * / step= 1.0 / ( 1.5 * distance); / * 発生させる内分点の総数 * / t = 0.0 ; while ( t く 1.0 ) { x = (int)( い x2 + ( 1.0 - t 戸 xl ) ; y = (int)( い y2 + ( 1.0 - t 戸 yl ) ; if( x > = 0 & & x < x_size2 & & y > = 0 & & y く y_size2 t = t 十 step; image2[y][x] + + ; & & image2[y][x] く MAX BRIGHTNESS ) { / ☆ p を計算するルーチン (draw_a_curve で使われている ) * / int —calc rho( double px, double py, double r_mg, int theta )

画像表示 / 出力装置によっては表現可能な階調数が限られており , そのよう な装置に画像を表示・出力する際は , 原画像を限られた階調数を用いてできる だけ忠実に表現するための手法が必要になる . このような擬似階調表示法とし て , 本節では , 組織的デイザ法 1 ) と誤差拡散法について述べる . 組織的デイザ法は , 例えば図 1.7 に示すようなデイザ行列を用いて , 階調画 像を 2 値画像 ( 各画素が黒画素あるいは白画素の画像 ) で表す手法である . 1 画像の入力・補正・出力 X Size2 = X S1ze1; y—size2 = y—sizel; gray—step = (double)GRAYLEVEL / FINAL—LEVEL; for (y = 0 ; y く y_size2; y 十 + ) for ( x = 0 ; x < x_size2; x 十十 ) image2 [y][x] = (unsigned char) ( trans—table[ imagel[y][x] ] * gray_step ) ; main( ) load—image—data( ) ; / * 原画像を imagel に読み込む * / smooth-histgram( ) ; ヒストグラムを平滑化して image2 へ響 save—image—data( ) ; / * image2 を保存する * / return 0 ; 1.3 画像の擬似階調表現法 1.3.1 組織的デイザ法 15 13 15 14 13 12 11 0 8 2 12 4 14 3 11 1 15 7 13 (a) Bayer 型 10 6 9 5 11 12 7 3 4 0 8 6 2 10 9 14 5 1 6 5 4 7 0 3 8 9 1 10 2 11 (b) 網点型 図 1.7 デイザ行列の例 (c) 渦巻型 組織的デイザ法の手順を次に示す .

140 6 立体認識 printf ( " 階調方向の倍率 ( 0 以上 1 以下 ) : " ) ; scanf("%f' &plot—ratio); plot—height = (int)( 256.0 * plot—ratio ) ; y—size2 = y—sizel 十 plot—height; だ描画ェリア ( image2[y][x]) の初期化☆ / for ( y = 0 ; y く y—size2; y 十十 ) for ( x = 0 ; x く x—size2; x + + ) image2[y][x] = 0 ; だグラフの表示 * / printf ( " 画像生成中です . しばらくお待ち下さい . *n"); for(y=O; y<y—sizel; y + + ) { if ( y % plot—step = = 0 ) { /*y が plot-step で割り切れるときだけ描く * / だ始めに曲線部分を黒で描いて隠面消去する * / for ( x = 0 ; x く x—sizel; x + + ) { draw—a—straight—line( x + Y / 2 , plot—height + Y, x + y / 2 , plot—height + y - (int)(imagel[y][x] * plot—ratio), 0 ) ; だ曲線部分を白で描く * / for ()= 0 ; x<x—sizel - 1 ; x + + ) { draw a—straight—line( x + Y / 2 , plot—height + Y - (int)(imagel[y][x] ま plot—ratio), x + y / 2 + 1 , plot—height + Y - (int)(imagel[y][x + 1] * plot—ratio)' MAX BRIGHTNESS ) ; main( ) load—image—data( ) ; draw_3D—graph( ) ; save—image—data( ) ; return 0 ; 実行結果 だ原画像の読み込み * / だ 3 次元グラフ画像の生成響 / * 画像の出力 * / 関数 draw-a-straight—line( ) は画像上に直線を描く関数であるが ,

36 2 画像のフィルタリング dft = dft / 2 ; void FFTI( double *re_part, double *im_part, int num of_data, int flag ) だデータの FFT(flag=1), 逆 FFT(flag=-1) を行う関数 だ double *re-part, *im-part; 元データ ( 実数部 , 虚数部 ) の先頭アドレス * / データ総数 , FFT ・逆 FFT を決めるフラグ / * int num—of—data, flag; int i, 」 , k, power; double unit_angle, step_angle; だ制御変数 , べき乗 * / だ角度用の変数 * / int dft, half, num—of_dft; だ DFT 点数 , その半数 , DFT 実行回数 * / だ DFT の入出力信号線の番号 * / int num out, num inl, num in2; double re-buf, im-buf, angle; だ実数部 , 虚数部 , 角度用作業変数 * / / * 次は計算結果を一時的に保存するための作業用配列 static double re—part—new[FFT MAX], im—part—new[FFT MAX] ; だ逆 FFT(flg=-I) のとき各データを num_of_data で割り * / だ複素共役をとる if(flag = = -1){ for ( i = 0 ; i く num of_data; i + + ) { * ( re_part 十 i ) = * ( re_part 十 i ) / num of_data; * ( im—part 十 i ) = * ( im—part 十 i ) / num of—data; だ num—of—data が 2 の何乗 ( power ) であるか調べる * / power = calc_power of_two( num of_data ) ; だ初期データのための元データの順番の入れ替えま / make initial_data( re_part, im_part, num of_data, power ) ; だ2 点 DFT,4 点 DFT,8 点 DFT, …の順に実行する unit—angle = 2.0 * PI / num of—data; /*=DFT の点数の初期値 , 2 倍してゆく * / dft= 2 ; for ( i = 0 ; i < power; i + + ) { /*DFT 点の DFT を行う / * i=0 - > 2 点 , i=l - > 4 点 , i=2 - > 8 点 , . … * / /*DFT の実行回数ツ num of dft = num_of data / dft; step—angle = unit—angle ま num of—dft; half = dft / 2 ; for (j = 0 ; j < num—of—dft; j + + ) { angle = 0.0 ;

2 画像のフィルタリング for ( k = 0 ; k く dft; k + + ) { num in1 = num out; num in2 = num out; num—out = j * dft 十 k; だ ( 係数が 1 の方の信号を num_inl としている ) * / だ num_inl と num—in2 から num_out を出力 * / num in1 } else num in2 ー if( k く half){ だ実数部 ( re ー ) ・虚数部 ( im つに分けて計算 * / = num out - halfi ー num_inl + half; 37 ( re_part 十 num in2 ) ; re buf = ま im_buf = * ( im_part 十 num in2 ) ; re—part_new[num out] = * ( re—part 十 num inl ) + re buf * cos(angle) + im buf * sin(angle); im—part_new[num out] = * ( im part 十 num inl ) + im buf * cos(angle) - re—buf * sin(angle); だ角度を更新 angle = angle 十 step—angle; / * 計算後のデータを元の配列に戻す * / for (j = 0 ; j く num of—data; j + + ) { * ( re—part 十 j ) = re—part—new[j] ; * ( im—part 十 j ) = im—part—new[j]; dft=dft*2;/*DFT の点数を 2 倍に更新 * / だ逆 FFT ( flg = -1 ) のとき各データの複素共役をとる響 if ( flag = for (j = 0 ; j く num of—data;j + + ) { * ( im—part 十 j ) = * ( im—part 十 j ) ; 補足説明 : プログラム中で本バッケージ "FFTI. h" をインクルードした後 , 元データ の実数部 , 虚数部をそれぞれ double data[FFT_MA 刈 , jdata[FFT MAX]; と 宣言すれば FFT を利用できる . ただし配列の最大値である FFT MAX は本