// Snake.v スネークモデルを用いた輪郭検出

#define P_NUM 60      // 円周のプロット数
#define BOX_SIZE 8    // ワーキングセットの大きさ
#define NO_USE -1     // エッジ情報外の領域にアクセスした際の画素値

// 節点のつながりを滑らかにするためのエネルギーを計算します
float CalcAngle(int pre_x, int pre_y,
		 int x,     int y,
		 int next_x,int next_y)
{
    int Vec1_x,Vec1_y;   // preとの間に引かれるベクトルの成分
    float Vec1_size;     // 同大きさ
    
    int Vec2_x,Vec2_y;   // nextとの間に引かれるベクトルの成分
    float Vec2_size;     // 同大きさ
    
    float I_product;     // 2つのベクトルのなす内積
    float CosTheta;       // 2つのベクトルのなす角度のコサイン値
    
    // 隣接2節点との間に引かれるベクトルの成分から、そのベクトルの
    // 大きさを計算します
    Vec1_x = x-pre_x;
    Vec1_y = y-pre_y;
    Vec1_size = sqrt((float)(Vec1_x*Vec1_x + Vec1_y*Vec1_y));
    // 隣接点と重なっている場合は0を返します
    if(Vec1_size == 0.){
	return 0.;
    }
    Vec2_x = x-next_x;
    Vec2_y = y-next_y;
    Vec2_size = sqrt((float)(Vec2_x*Vec2_x + Vec2_y*Vec2_y));
    // 隣接点と重なっている場合は0を返します
    if(Vec2_size == 0.) return 0.;

    // 2つのベクトルの成分から、その内積を計算します
    I_product = (float)(Vec1_x*Vec2_x + Vec1_y*Vec2_y);
    
    // 2ベクトルの内積と大きさから、コサイン値を求めます
    CosTheta = I_product/(Vec1_size*Vec2_size);
    
    // 角度から求められるエネルギーとして値を返します
    if(CosTheta > 0) return CosTheta;
    else             return CosTheta/2;
}

// 節点間の距離を等しくしようとするエネルギーを計算します
float CalcDist(int pre_x, int pre_y,
	       int x,     int y,
	       int next_x,int next_y)
{
    int Vec1_x,Vec1_y;   // preとの間に引かれるベクトルの成分
    float Vec1_size;     // 同大きさ
    
    int Vec2_x,Vec2_y;   // nextとの間に引かれるベクトルの成分
    float Vec2_size;     // 同大きさ
    
    float Diff;          // 2つのベクトルの大きさの差
    
    // 隣接2節点との間に引かれるベクトルの成分から、そのベクトルの
    // 大きさを計算します
    Vec1_x = x-pre_x;
    Vec1_y = y-pre_y;
    Vec1_size = sqrt((float)(Vec1_x*Vec1_x + Vec1_y*Vec1_y));
    
    Vec2_x = x-next_x;
    Vec2_y = y-next_y;
    Vec2_size = sqrt((float)(Vec2_x*Vec2_x + Vec2_y*Vec2_y));
    
    // 2つのベクトルの大きさの差を正規化したものをエネルギーとして返します
    return fabs(Vec1_size - Vec2_size)/(Vec1_size + Vec2_size);
}

// 節点を円の外側に移動させようとするエネルギー、すなわち
// 円を徐々に膨張させるエネルギーを計算します
float CalcMove(int pre_x,int pre_y,
	       int new_x, int new_y,	       
	       int old_x, int old_y)
{
    int Vec1_x,Vec1_y;   // newとの間に引かれるベクトルの成分
    float Vec1_size;     // 同大きさ
    
    int Vec2_x,Vec2_y;   // oldとの間に引かれるベクトルの成分
    float Vec2_size;     // 同大きさ
    
    float O_product;     // 2つのベクトルのなす外積
    float SinTheta;       // 2つのベクトルのなす角度のサイン値
    
    // 隣接2節点との間に引かれるベクトルの成分から、そのベクトルの
    // 大きさを計算します
    Vec1_x = new_x-pre_x;
    Vec1_y = new_y-pre_y;
    Vec1_size = sqrt((float)(Vec1_x*Vec1_x + Vec1_y*Vec1_y));
    // 隣接点と重なっている場合は0を返します
    if(Vec1_size == 0.) return 0.;
    
    Vec2_x = old_x-pre_x;
    Vec2_y = old_y-pre_y;
    Vec2_size = sqrt((float)(Vec2_x*Vec2_x + Vec2_y*Vec2_y));
    // 隣接点と重なっている場合は0を返します
    if(Vec2_size == 0.) return 0.;
    
    // 2つのベクトルの成分から、その外積を計算します
    O_product = (float)(Vec1_x*Vec2_y - Vec1_y*Vec2_x);
    
    // 2ベクトルの外積と大きさから、サイン値を求めます
    SinTheta = O_product/(Vec1_size*Vec2_size);
    
    // 求められた値をエネルギーとして
    return -SinTheta;
}

// 入力されたエッジ情報、円周情報から注目円周節点の
// 新たな位置を決定します
void SelectNewPoint(
		    int l_image[3][3],int v,
		    int pre_x,int pre_y,
		    int x[1],int y[1],
		    int next_x,int next_y)
{
    int   i,j;
    int   old_x,old_y;    
    float Efield,Eangle,Edist,Emove;  // 輪郭抽出に用いるエネルギーです
    float Eall,Eall_min = 1000000.;  // エネルギーの総和とその最小値です
    float MU_FIELD      = 40.;  // 正則化理論に基づき新たな点の 
    float MU_ANGLE      = 5.;  // 位置を決める際に利用する各  
    float MU_DISTANCE   = 5.;  // エネルギー関数の重みとなる  
    float MU_MOVE       = 1.;  // 係数を定義します            


    old_x = x[0];
    old_y = y[0];
    // 近傍8箇所と元の位置のうちでこの時点で最も総エネルギーが低い
    // 位置を探します
    for(j = old_y - 1; j < old_y + 2; j++){
	for(i = old_x - 1; i < old_x + 2; i++){
	    // エッジ画像領域内にある点についてのみ計算します
	    if(l_image[i - old_x + 1][j - old_y + 1] != NO_USE){
		// 位置エネルギーを計算します
		Efield = (float)l_image[i-old_x+1][j-old_y+1]*MU_FIELD/255;
		// 隣接2節点となす角度を直線に近付ける、すなわち、円周の
		// その部分を滑らかにするためのエネルギーを計算します
		Eangle = CalcAngle(pre_x,pre_y,i,j,next_x,next_y)*MU_ANGLE;
		// 注目節点を隣接2節点の中点に近付けようとするエネルギー
		// を計算します
		Edist = CalcDist(pre_x,pre_y,i,j,next_x,next_y)*MU_DISTANCE;
		// 円周を徐々に大きくするため、つまり、膨張させるための
		// エネルギーを計算します

		Emove = -CalcMove(pre_x,pre_y,i,j,old_x,old_y)*MU_MOVE;
		Emove += CalcMove(next_x,next_y,i,j,old_x,old_y)*MU_MOVE;
		
		// 全てのエネルギーの総和を計算します
		Eall = Efield + Eangle + Edist + Emove;
	
		// 過去のEallと比較し、小さければ、その値を最小値として保存します
		if(Eall < Eall_min){
		    Eall_min = Eall;
		    x[0] = i;
		    y[0] = j;
		}
	    }
	}
    }
}

// スネークモデルを用いた輪郭検出プログラム
// 画像処理モジュール本体部分

module Snake(edge:input,result:output,circle_x,circle_y,roop_num:parameter)
int edge on box[BOX_SIZE][BOX_SIZE] cache 1;
int result on box[BOX_SIZE][BOX_SIZE];
int circle_x[P_NUM];  // 円周のx座標を格納する配列
int circle_y[P_NUM];  // 円周のy座標を格納する配列
int roop_num[1];      // 円周を収縮させる回数
#vios mutex circle_x off max
#vios mutex circle_y off max
{
    int temp_x[P_NUM];
    int temp_y[P_NUM];
#vios mutex temp_x off max
#vios mutex temp_y off max

    int flag;     // 処理終了を表すフラグ
#vios mutex flag off max
    {
	int i;
	for(i = 0; i < P_NUM; i++){
	    temp_x[i] = 0;
	    temp_y[i] = 0;
	}
    }

    for(flag = 0; flag < roop_num[0]; flag++){
	parallel{
	    int i,j,v;     // ループ変数
	    int my_x,my_y;
	    int x[1],y[1],pre_x,pre_y,next_x,next_y;
	    int l_image[3][3];  // 注目画素と周辺画素を格納する配列
	    int wdt,hgt;        // エッジ情報の画像の大きさ
	    int old_x,old_y;    // 節点の元の座標

	    my_x = hotx(edge);    // ワーキングセットの先頭座標の
	    my_y = hoty(edge);    // 位置を取得します
	    
	    wdt = width(edge);    // エッジ画像の縦、横方向の大きさを
	    hgt = height(edge);   // 取得します

	    // 全ての円周座標に対して、その座標がワーキングセット内に
	    // 含まれているかどうか調べます
	    for(v = 0; v < P_NUM; v++){
		if(my_x <= circle_x[v] && circle_x[v] < my_x+BOX_SIZE){
		    if(my_y <= circle_y[v] && circle_y[v] < my_y+BOX_SIZE){
			// 使用する座標情報をローカル変数に格納します
			x[0] = circle_x[v];
			y[0] = circle_y[v];
			old_x = x[0];
			old_y = y[0];
			if(v == 0){
			    pre_x = circle_x[P_NUM-1];
			    pre_y = circle_y[P_NUM-1];
			    next_x = circle_x[v+1];
			    next_y = circle_y[v+1];
			}
			else if(v == P_NUM-1){
			    pre_x = circle_x[v-1];
			    pre_y = circle_y[v-1];
			    next_x = circle_x[0];
			    next_y = circle_y[0];
			}
			else{
			    pre_x = circle_x[v-1];
			    pre_y = circle_y[v-1];
			    next_x = circle_x[v+1];
			    next_y = circle_y[v+1];
			}
			// 注目座標の画素と、その周囲8近傍の画素を
			// ローカル変数に格納します
			for(j = old_y - 1; j < old_y + 2; j++){
			    // エッジ画像の領域内の場合はその位置の
			    // 画素値を格納します
			    if(0 <= j && j < hgt){
				for(i = old_x - 1; i < old_x + 2; i++){
				    if(0 <= i && i < wdt){
					l_image[i - old_x + 1][j - old_y + 1] = edge[[i]][[j]];
				    }
				    // それ以外の場合はNO_USEを代入します
				    else{
					l_image[i - old_x + 1][j - old_y + 1] = NO_USE;
				    }
				}
			    }
			    else{
			        l_image[i - old_x + 1][j - old_y + 1] = NO_USE;
			    }
			}
			// 正則化理論に基づき、新たな位置を探します
			SelectNewPoint(l_image,v,pre_x,pre_y,x,y,next_x,next_y);
			
			if(x[0] != old_x) temp_x[v] = x[0];
			if(y[0] != old_y) temp_y[v] = y[0];
		    }
		}
	    }
	}//  1st parallel end

	{// 更新された値のみtempに入るため、更新されていればその値をcircleに書き戻します
	    int i;
	    for(i = 0; i < P_NUM; i++){
		if(temp_x[i] != 0){
		    circle_x[i] = temp_x[i];
		    temp_x[i] = 0;
		}
		if(temp_y[i] != 0){
		    circle_y[i] = temp_y[i];
		    temp_y[i] = 0;
		}
	    }
	}
    }
    // 最終的に節点が止まった位置情報をプロットした画像を作成します
    parallel{
	int i;
	int x,y;
	int my_x = hotx(edge);
	int my_y = hoty(edge);
	
	for(y = 0; y < BOX_SIZE; y++){
	    for(x = 0; x < BOX_SIZE; x++){
		if(x < width(result) && y < height(result))
		    result[x][y] = 255;
	
		for(i = 0; i < P_NUM; i++){
		    if(circle_x[i] == my_x+x && circle_y[i] == my_y+y)
			result[x][y] = 0;
		}
	    }
	}
    }
}

スネークメインフローへ戻る

VIOS トップページ