//-----------------------------------------------------------------
// clustering.cpp:
//    クラスタリング用関数
//               Last Update: <2004/12/22 22:08:50 A.Murakami>
//-----------------------------------------------------------------
#include "clustering.h"
//-----------------------------------------------------------------
// 画像中の白画素のデータ取得
//-----------------------------------------------------------------
int obtain_points_data(LPBYTE inb,ClutData& cltd)
{
    int x,y;
    cltd.num_pnts=0; // 取得ポイント数
    for(y=0;y<iHeight;y++) for(x=0;x<iWidth;x++){
		if(inb[x+y*iWidth] == iWHITE){
			if(cltd.num_pnts<MAX_PNTS-2){
				cltd.num_pnts ++;
				cltd.points[cltd.num_pnts-1].x=x;
				cltd.points[cltd.num_pnts-1].y=y;
			}
		}
    }
    if(cltd.num_pnts<5){
        MessageBox(NULL,"点列が少な過ぎます. 画像を作り直して下さい",0,0);
        return FALSE;
    }
	if(cltd.num_pnts>=MAX_PNTS){
        MessageBox(NULL,"点列が多すぎます. 画像を作り直して下さい",0,0);
        return FALSE;
    }
	return TRUE;
}
//-----------------------------------------------------------------
// 取得した点データをランダムに入れ替える.
//-----------------------------------------------------------------
void shuffle_points_data(ClutData& cltd)
{
    int i,num1,num2;
    int sw_max;
    // 乱数の初期化
	init_random();
    // 点の総数の２倍の回数の入れ替え
    sw_max = cltd.num_pnts * 2;
    for(i=0;i<sw_max;i++){
        num1 = random_int(cltd.num_pnts);
        do { 
            num2 = random_int(cltd.num_pnts);
        } while (num2==num1);
        // No.num1 と No.num2 の点を入れ替え
		swap_point(cltd.points[num1],cltd.points[num2]);
    }
}
//-----------------------------------------------------------------
// [y][x]ブロックのベクトルとNo.n 代表ベクトル間の距離計算
//-----------------------------------------------------------------
double calc_vdistance(int y,int x,int n,VQData& vqd)
{
    int i; double sum=0.0;
    for(i=0;i<VQ_DIM;i++)
        sum += (vqd.vector[y][x][i]-vqd.main_vector[n][i]) *
			(vqd.vector[y][x][i]-vqd.main_vector[n][i]);
    return (sqrt(sum));
}
//-----------------------------------------------------------------
// 画像を 4x4 画素のブロックに分割してベクトル化
//-----------------------------------------------------------------
void obtain_vector_data(LPBYTE inb,VQData& vqd)
{
    int x,y,xs,ys,i,j;
    for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
		xs=VQ_BLKSZ * x; ys=VQ_BLKSZ * y;
		for(i=0;i<VQ_BLKSZ;i++) for(j=0;j<VQ_BLKSZ;j++){
			if(ip_safe(xs+j,ys+i,iWidth,iHeight))
				vqd.vector[y][x][i*VQ_BLKSZ+j] =
					inb[XY1D(xs+j,ys+i,iWidth)];
			else
				vqd.vector[y][x][i*VQ_BLKSZ+j] = 0;
        }
    }
}
//-----------------------------------------------------------------
// [y][x]ブロックのベクトルとNo.n 代表ベクトル間の距離計算
//-----------------------------------------------------------------
double calc_c_vdistance(int y,int x,int n,CVQData& vqd)
{
    int i; double sum=0.0;
    for(i=0;i<CVQ_DIM;i++)
        sum += (vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[i]-vqd.main_vector[n][i]) *
			(vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[i]-vqd.main_vector[n][i]);
    return (sqrt(sum));
}
//-----------------------------------------------------------------
// 画像を 4x4 画素のブロックに分割してベクトル化
//-----------------------------------------------------------------
void obtain_c_vector_data(LPBYTE inb,CVQData& vqd)
{
    int x,y,xs,ys,i,j;
    for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
		xs=VQ_BLKSZ * x; ys=VQ_BLKSZ * y;
		for(i=0;i<VQ_BLKSZ;i++) for(j=0;j<VQ_BLKSZ;j++){
			if(ip_safe(xs+j,ys+i,iWidth,iHeight)){
				vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[XY1D(j*3,i,VQ_BLKSZ*3)+0] =
					inb[XY1D((xs+j)*3,ys+i,iLength)+2];
				vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[XY1D(j*3,i,VQ_BLKSZ*3)+1] =
					inb[XY1D((xs+j)*3,ys+i,iLength)+1];
				vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[XY1D(j*3,i,VQ_BLKSZ*3)+2] =
					inb[XY1D((xs+j)*3,ys+i,iLength)+0];
			} else {
				vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[XY1D(j*3,i,VQ_BLKSZ*3)+0] = 0;
				vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[XY1D(j*3,i,VQ_BLKSZ*3)+1] = 0;
				vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[XY1D(j*3,i,VQ_BLKSZ*3)+2] = 0;
			}
        }
    }
}
//-----------------------------------------------------------------
// NN(Nearest Neighbour) 法によるクラスタリング
//-----------------------------------------------------------------
void clustering_NNmethod(LPBYTE inb,int clstr_r,int dopt)
{
    int i,j,min_num;
    double min_distance,distance;
	ClutData cltd;
	//--------------------------------------------------
	// 画像中の点を取得
	//--------------------------------------------------
	if(!obtain_points_data(inb,cltd)) return;
	//--------------------------------------------------
	// No.0 の点を No.0 のクラスタ代表点に設定
	//--------------------------------------------------
	cltd.num_clstr=0; // クラスタの総数 - 1
	cltd.cclstr[cltd.num_clstr].x = cltd.points[0].x;
	cltd.cclstr[cltd.num_clstr].y = cltd.points[0].y;
	cltd.p_clstr[0]=0;
	cltd.num_clstr ++;
	//--------------------------------------------------
	// No.1 から No.num_pnts-1 までの点を分類する
	//--------------------------------------------------
	for(i=1;i<cltd.num_pnts;i++){
		// 最も近いクラスタの計算
		//   No: min_num, 距離: min_distance
		min_distance=10000.0;
		for(j=0;j<cltd.num_clstr;j++){
			distance=calc_distance(cltd.points[i],cltd.cclstr[j]);
			if(distance<min_distance){
				min_distance=distance;
				min_num=j;
			}
		}
		// No.i の点をクラスタに加える
		if(min_distance<clstr_r){
			cltd.p_clstr[i]=min_num;
			// クラスタ重心の移動
			cltd.cclstr[min_num].x+=cltd.points[i].x;
			cltd.cclstr[min_num].x/=2;
			cltd.cclstr[min_num].y+=cltd.points[i].y;
			cltd.cclstr[min_num].y/=2;
		}
		// No.i の点は新しいクラスタに加える
		else {
			cltd.cclstr[cltd.num_clstr].x=cltd.points[i].x;
			cltd.cclstr[cltd.num_clstr].y=cltd.points[i].y;
			cltd.p_clstr[i]=cltd.num_clstr;
			cltd.num_clstr ++;
		}
	}
	//--------------------------------------------------
	// クラスタ位置の表示
	//--------------------------------------------------
	// 入力点と各点までの線分を描く
	if(dopt){
		for(i=0;i<cltd.num_pnts;i++){
			draw_circle(inb,cltd.points[i],5);
			if(dopt==CLUT_DALL)
				draw_line(inb,cltd.points[i],cltd.cclstr[cltd.p_clstr[i]]);
		}
	}
	// クラスタ中心の表示
	for(i=0;i<cltd.num_clstr;i++){
		draw_circle(inb,cltd.cclstr[i],clstr_r);
	}
}
//-----------------------------------------------------------------
// 階層的クラスタリング: 樹形図的併合方法によるクラスタリング
//-----------------------------------------------------------------
//   1. パターン間の距離を計算
//   2. 最も近接する2つのクラスタを併合し新しいクラスタをつくる
//   3. 新しいクラスタと残りのクラスタとの距離を計算
//   4. 最小のクラスタ間距離が閾値より大きくなれば終了(2,3の反復)
//-----------------------------------------------------------------
void clustering_Tree(LPBYTE inb,int clstr_r,int dopt)
{
	int end_loop = 1000;                // 最大反復回数
	int min_num1,min_num2;              // 併合するクラスタ番号
	int radius_clstr[MAX_PNTS];         // クラスタの意味上の半径
	LPUINT dist_clstr;                  // クラスタの意味上の半径
	int i,j,dx,dy,counter,finish,sum_number;
	double min_distance,distance;
	double sum_x,sum_y;
	ClutData cltd;
	//--------------------------------------------------
	// 画像中の点を取得
	//--------------------------------------------------
	if(!obtain_points_data(inb,cltd)) return;
	dist_clstr = (LPUINT)GlobalAlloc(GPTR,sizeof(UINT)*MAX_PNTS*MAX_PNTS);
	//--------------------------------------------------
	// 初期のクラスタ中心の設定
	//--------------------------------------------------
	cltd.num_clstr=cltd.num_pnts;
    for(i=0;i<cltd.num_clstr;i++){
		cltd.cclstr[i].x=cltd.points[i].x;
		cltd.cclstr[i].y=cltd.points[i].y;
		radius_clstr[i]=0; // クラスタ半径
		cltd.p_clstr[i] =i;
	}
	// クラスタ間の距離を求める
	for(i=0;i<cltd.num_clstr;i++) for(j=i+1;j<cltd.num_clstr;j++){
		distance=calc_distance(cltd.cclstr[i],cltd.cclstr[j]);
		dist_clstr[j+i*cltd.num_clstr] = rint(distance);
	}
	//--------------------------------------------------
	// クラスタリング処理
	//--------------------------------------------------
    finish=0; counter=0;
    while(!finish){
        counter++;
        // クラスタ間の最小距離を求める
		min_distance=10000.0;
		for(i=0;i<cltd.num_clstr;i++) for(j=i+1;j<cltd.num_clstr;j++){
			if(dist_clstr[j+i*cltd.num_clstr]<min_distance){
				min_distance=dist_clstr[j+i*cltd.num_clstr];
				min_num1=i; min_num2=j;
			}
		}
		// 最小距離が許容範囲以下ならクラスタを併合
		if(min_distance < clstr_r){
			// 併合されたクラスタの除去
			for(i=min_num2;i<cltd.num_clstr-1;i++){
				cltd.cclstr[i].x = cltd.cclstr[i+1].x;
				cltd.cclstr[i].y = cltd.cclstr[i+1].y;
			}
			// パターンの所属クラスタの変更
			for(i=0;i<cltd.num_pnts;i++){
				if(cltd.p_clstr[i]>min_num2)
					cltd.p_clstr[i]--;
				else if(cltd.p_clstr[i]==min_num2)
					cltd.p_clstr[i]=min_num1;
			}
			// 新規クラスタ重心の設定
			sum_number = 0;
			sum_x = sum_y = 0;
			for(i=0;i<cltd.num_pnts;i++){
				if(cltd.p_clstr[i]==min_num1){
					sum_number++;
					sum_x += cltd.points[i].x; sum_y += cltd.points[i].y;
				}
			}
			cltd.cclstr[min_num1].x=rint(sum_x/(double)sum_number);
            cltd.cclstr[min_num1].y=rint(sum_y/(double)sum_number);
			// 併合された距離配列の削除
			for(i=0;i<cltd.num_clstr-1;i++){
				for(j=i+1;j<cltd.num_clstr-1;j++){
					dx = (j>=min_num2)?1:0;
					dy = (i>=min_num2)?1:0;
					dist_clstr[j+i*(cltd.num_clstr-1)] =
						dist_clstr[(j+dx)+(i+dy)*cltd.num_clstr];
				}
			}
			// クラスタ数を減らす
			cltd.num_clstr--;
			// 距離の再計算
			for(i=0;i<=min_num1;i++) for(j=i+1;j<cltd.num_clstr;j++){
				if(i==min_num1 || j==min_num1){
					distance=calc_distance(cltd.cclstr[i],
										   cltd.cclstr[j]);
					dist_clstr[j+i*cltd.num_clstr] = rint(distance);
				}
			}
		} else finish = 1;
		if(counter > end_loop) finish = 1;
	}
	//--------------------------------------------------
	// 各クラスタの意味上の半径を求める
	//--------------------------------------------------
	for(i=0;i<cltd.num_clstr;i++) for(j=0;j<cltd.num_pnts;j++){
		if(cltd.p_clstr[j] == i){
			distance=calc_distance(cltd.cclstr[i],cltd.points[j]);
			if(distance > radius_clstr[i])
				radius_clstr[i]=rint(distance);
		}
	}
	//--------------------------------------------------
	// クラスタ位置の表示
	//--------------------------------------------------
	// 入力点と各点データまでの線分を描く
	if(dopt){
		for(i=0;i<cltd.num_pnts;i++){
			draw_circle(inb,cltd.points[i],5);
			if(dopt==CLUT_DALL)
				draw_line(inb,cltd.points[i],
						  cltd.cclstr[cltd.p_clstr[i]]);
		}
	}
	// クラスタ中心の表示
	for(i=0;i<cltd.num_clstr;i++){
		draw_circle(inb,cltd.cclstr[i],radius_clstr[i]);
	}
	// 後片付け
	GlobalFree(dist_clstr);
}
//-----------------------------------------------------------------
// K-平均法によるクラスタリング
//-----------------------------------------------------------------
//   1. K個のクラスタの初期中心位置を与える(ランダム)
//   2. 各点を近接するクラスタに分類
//   3. 各クラスタに属する点の中心位置の計算
//   4. 初期中心と比較し、一致・許容範囲内なら終了(2,3の反復)
//-----------------------------------------------------------------
void clustering_Kmean(LPBYTE inb,int clstr_r,int K_number,int dopt)
{
    int i,j,min_num,sum_number;
	int counter,finish;
	double min_distance,distance,sum_x,sum_y;
	int radius_clstr[MAX_PNTS];       // クラスタの意味上の半径
    POINT init_cntr_clstr[MAX_PNTS];  // クラスタの初期位置
	ClutData cltd;
	//--------------------------------------------------
	if(K_number>MAX_CLSTR) return;
	//--------------------------------------------------
	// 画像中の点を取得
	//--------------------------------------------------
	if(!obtain_points_data(inb,cltd)) return;
	//--------------------------------------------------
	// 初期のクラスタ中心の設定
	//--------------------------------------------------
	cltd.num_clstr=K_number;
	shuffle_points_data(cltd); // ランダムに入れ替え
    for(i=0;i<cltd.num_clstr;i++){
		cltd.cclstr[i].x=cltd.points[i].x;
		cltd.cclstr[i].y=cltd.points[i].y;
	}
	//--------------------------------------------------
	// クラスタリング処理
	//--------------------------------------------------
    finish=0; counter=0;
    while(!finish){
        counter++;
        // クラスタの初期代表点への代入 
        for(i=0;i<cltd.num_clstr;i++){
			init_cntr_clstr[i].x=cltd.cclstr[i].x;
			init_cntr_clstr[i].y=cltd.cclstr[i].y;
            radius_clstr[i]=clstr_r; // defalt 値
        }
        // 各点の所属クラスタ番号の初期化
        for(i=0;i<cltd.num_pnts;i++)
            cltd.p_clstr[i]=-1;
        // 各点の所属クラスタ
        for(i=0;i<cltd.num_pnts;i++){
            // 各クラスタ代表点との距離
            min_distance=10000.0;
            for(j=0;j<cltd.num_clstr;j++){
                distance=calc_distance(cltd.points[i],cltd.cclstr[j]);
                if(distance<min_distance){
                    min_distance=distance;
                    min_num=j;
                }
            }
            // No.i の点をクラスタ No.min_num に加える
            cltd.p_clstr[i]=min_num;
        }
        // 各クラスタの代表点を修正
        for(i=0;i<cltd.num_clstr;i++){
			sum_number=0;
            sum_x=sum_y=0.0;
            for(j=0;j<cltd.num_pnts;j++){
                if(cltd.p_clstr[j] == i){
                    sum_number++;
                    sum_x=sum_x + cltd.points[j].x;
                    sum_y=sum_y + cltd.points[j].y;
                }
            }
            cltd.cclstr[i].x=rint(sum_x/(double)sum_number);
            cltd.cclstr[i].y=rint(sum_y/(double)sum_number);
        }
        // 各クラスタの意味上の半径を求める
        for(i=0;i<cltd.num_clstr;i++){
            for(j=0;j<cltd.num_pnts;j++){
                if(cltd.p_clstr[j] == i){
                    distance=calc_distance(cltd.cclstr[i],cltd.points[j]);
                    if(distance > radius_clstr[i])
                        radius_clstr[i]=rint(distance);
                }
            }
        }
        // クラスタ中心が移動したかどうかのチェック
        finish=1;
        for(i=0;i<cltd.num_clstr;i++){
            distance=calc_distance(init_cntr_clstr[i],cltd.cclstr[i]);
            if(distance > KC_MINERR) finish=0; // 許容誤差
        }
	}
	//--------------------------------------------------
	// クラスタ位置の表示
	//--------------------------------------------------
	// 入力点と各点データまでの線分を描く
	if(dopt){
		for(i=0;i<cltd.num_pnts;i++){
			draw_circle(inb,cltd.points[i],5);
			if(dopt==CLUT_DALL)
				draw_line(inb,cltd.points[i],cltd.cclstr[cltd.p_clstr[i]]);
		}
	}
	// クラスタ中心の表示
	for(i=0;i<cltd.num_clstr;i++){
		draw_circle(inb,cltd.cclstr[i],radius_clstr[i]);
	}
}

//-----------------------------------------------------------------
// ベクトル量子化[K平均法によるクラスタリング]: 濃淡画像
//-----------------------------------------------------------------
double vector_quantize(LPBYTE inb,int K_number)
// int K_number: K値(生成するクラスタの総数),default=100
// 返り値: 圧縮率[圧縮後のサイズ/元サイズ]
{
    int i,j,x,y,xs,ys;
    int counter,min_num,finish;
	int original_data_size,coded_data_size;
	double min_distance,distance,compression_rate;
    double sum[VQ_DIM],sum_number,difference;
	BYTE init_main_vector[MAX_CLSTR][VQ_DIM]; // 代表ベクトル
	VQData vqd;
	//-----------------------------------------------------------------
	if(K_number>MAX_CLSTR) return 0;
	//-----------------------------------------------------------------
	// ブロック数
	//-----------------------------------------------------------------
	vqd.x_block = iWidth  / VQ_BLKSZ;
    vqd.y_block = iHeight / VQ_BLKSZ;
	if(iWidth  % VQ_BLKSZ != 0) vqd.x_block ++;
	if(iHeight % VQ_BLKSZ != 0) vqd.y_block ++;
	//-----------------------------------------------------------------
	// ベクトルデータの取得
	//-----------------------------------------------------------------
	obtain_vector_data(inb,vqd);
	//-----------------------------------------------------------------
    // 仮代表ベクトルの設定
	//-----------------------------------------------------------------
	init_random(); // 乱数の初期化
    for(i=0;i<K_number;i++){
        x=random_int(vqd.x_block);
        y=random_int(vqd.y_block);
        for(j=0;j<VQ_DIM;j++){
            vqd.main_vector[i][j] = vqd.vector[y][x][j];
        }
    }
	//-----------------------------------------------------------------
    // K-平均法によるベクトルのクラスタリング
	//-----------------------------------------------------------------
    finish=0; counter=0;
    while(!finish){
        counter++;
		//--------------------------------------------------
        // クラスタの初期代表点への代入
		//--------------------------------------------------
        for(i=0;i<K_number;i++) for(j=0;j<VQ_DIM;j++)
			init_main_vector[i][j] = vqd.main_vector[i][j];
		//--------------------------------------------------
        // 各ブロックの所属クラスタ番号の初期化
		//--------------------------------------------------
        for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++)
			vqd.b_clstr[y][x] = -1;
		//--------------------------------------------------
        // 各点の所属クラスタを決定
		//--------------------------------------------------
        for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
			// 各クラスタ代表点との距離
			min_distance=iWHITE * iWHITE * VQ_DIM;
			for(i=0;i<K_number;i++){
				distance=calc_vdistance(y,x,i,vqd);
				if(distance<min_distance){
					min_distance = distance;
					min_num = i;
				}
			}
			// [y][x]ブロックをクラスタ No.min_num へ
			vqd.b_clstr[y][x] = min_num;
        }
		//--------------------------------------------------
		// 各クラスタの代表点を修正する
		//--------------------------------------------------
        for(i=0;i<K_number;i++){
            for(j=0;j<VQ_DIM;j++) sum[j]=0.0;
            sum_number=0.0;
            for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
				if(vqd.b_clstr[y][x]==i){
					sum_number = sum_number + 1.0;
					for(j=0;j<VQ_DIM;j++)
						sum[j] += vqd.vector[y][x][j];
                }
            }
            for(j=0;j<VQ_DIM;j++){
                vqd.main_vector[i][j]=(BYTE)(sum[j] / sum_number);
            }
        }
		//--------------------------------------------------
        // クラスタ中心が移動したかどうかのチェック
		//--------------------------------------------------
        finish=1;
        for(i=0;i<K_number;i++){
            difference=0.0;   
            for(j=0;j<VQ_DIM;j++){
                difference +=
					abs(init_main_vector[i][j]-vqd.main_vector[i][j]);
            }
            if(difference>VQ_MINERR) finish=0;
        }
        if(finish == 1){ // 処理終了
            // 各ブロックを代表ベクトルに置換
            for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
				xs=VQ_BLKSZ * x;   
				ys=VQ_BLKSZ * y;
				for(i=0;i<VQ_BLKSZ;i++) for(j=0;j<VQ_BLKSZ;j++){
					if(ip_safe(xs+j,ys+i,iWidth,iHeight))
						inb[XY1D(xs+j,ys+i,iWidth)] =
							vqd.main_vector[vqd.b_clstr[y][x]][i*VQ_BLKSZ+j];
					else
						inb[XY1D(xs+j,ys+i,iWidth)] = 0;
                }
            }
        }
    }
    // データ量・圧縮率関係の計算と表示
    original_data_size = iSize * 8;
    i=K_number-1; j=0; // j:代表ベクトルの指定に必要なビット数
    do{ i>>=1; j++; } while(i > 0);
	// 圧縮後のサイズ
    coded_data_size = VQ_BLKSZ * VQ_BLKSZ * K_number * 8 +
		j * vqd.x_block * vqd.y_block; // 代表ベクトル数
	// 圧縮率
    compression_rate=(double)coded_data_size /
		(double)original_data_size * 100.0;
	
	return compression_rate;
}

//-----------------------------------------------------------------
// ベクトル量子化[K平均法によるクラスタリング]: カラー画像
//-----------------------------------------------------------------
double c_vector_quantize(LPBYTE inb,int K_number)
// int K_number: K値(生成するクラスタの総数),default=100
// 返り値: 圧縮率[圧縮後のサイズ/元サイズ]
{
#if 0 // 各RGB成分ごと
	double comp_rate;
	LPBYTE iR=GetGray(COLOR_R);
	LPBYTE iG=GetGray(COLOR_G);
	LPBYTE iB=GetGray(COLOR_B);
	comp_rate = 0;
	comp_rate += vector_quantize(iR,K_number);
	comp_rate += vector_quantize(iG,K_number);
	comp_rate += vector_quantize(iB,K_number);
	comp_rate /= 3.0;
	RGBToColor(inb,iR,iG,iB);
	// 後片付け
	GlobalFree(iR); GlobalFree(iG); GlobalFree(iB);
	return comp_rate;
#else // RGB成分をベクトルとして扱う
    int i,j,x,y,xs,ys;
    int counter,min_num,finish;
	int original_data_size,coded_data_size;
	double min_distance,distance,compression_rate;
    double sum[CVQ_DIM],sum_number,difference;
	BYTE init_main_vector[MAX_CLSTR][CVQ_DIM]; // 代表ベクトル
	CVQData vqd;
	//-----------------------------------------------------------------
	if(K_number>MAX_CLSTR) return 0;
	//-----------------------------------------------------------------
	// ブロック数
	//-----------------------------------------------------------------
	vqd.x_block = iWidth  / VQ_BLKSZ;
    vqd.y_block = iHeight / VQ_BLKSZ;
	if(iWidth  % VQ_BLKSZ != 0) vqd.x_block ++;
	if(iHeight % VQ_BLKSZ != 0) vqd.y_block ++;
	//-----------------------------------------------------------------
	// ベクトルデータの取得
	//-----------------------------------------------------------------
	vqd.cv = (C_vector*)GlobalAlloc(GPTR,sizeof(C_vector)*VQ_MAXBLK*VQ_MAXBLK);
	obtain_c_vector_data(inb,vqd);
	//-----------------------------------------------------------------
    // 仮代表ベクトルの設定
	//-----------------------------------------------------------------
	init_random(); // 乱数の初期化
    for(i=0;i<K_number;i++){
        x=random_int(vqd.x_block);
        y=random_int(vqd.y_block);
        for(j=0;j<CVQ_DIM;j++){
            vqd.main_vector[i][j] = vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[j];
        }
    }
	//-----------------------------------------------------------------
    // K-平均法によるベクトルのクラスタリング
	//-----------------------------------------------------------------
    finish=0; counter=0;
    while(!finish){
        counter++;
		//--------------------------------------------------
        // クラスタの初期代表点への代入
		//--------------------------------------------------
        for(i=0;i<K_number;i++) for(j=0;j<CVQ_DIM;j++)
			init_main_vector[i][j] = vqd.main_vector[i][j];
		//--------------------------------------------------
        // 各ブロックの所属クラスタ番号の初期化
		//--------------------------------------------------
        for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++)
			vqd.b_clstr[y][x] = -1;
		//--------------------------------------------------
        // 各点の所属クラスタを決定
		//--------------------------------------------------
        for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
			// 各クラスタ代表点との距離
			min_distance=iWHITE * iWHITE * CVQ_DIM;
			for(i=0;i<K_number;i++){
				distance=calc_c_vdistance(y,x,i,vqd);
				if(distance<min_distance){
					min_distance = distance;
					min_num = i;
				}
			}
			// [y][x]ブロックをクラスタ No.min_num へ
			vqd.b_clstr[y][x] = min_num;
        }
		//--------------------------------------------------
		// 各クラスタの代表点を修正する
		//--------------------------------------------------
        for(i=0;i<K_number;i++){
            for(j=0;j<CVQ_DIM;j++) sum[j]=0.0;
            sum_number=0.0;
            for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
				if(vqd.b_clstr[y][x]==i){
					sum_number = sum_number + 1.0;
					for(j=0;j<CVQ_DIM;j++)
						sum[j] += vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[j];
                }
            }
            for(j=0;j<CVQ_DIM;j++){
                vqd.main_vector[i][j]=(BYTE)(sum[j] / sum_number);
            }
        }
		//--------------------------------------------------
        // クラスタ中心が移動したかどうかのチェック
		//--------------------------------------------------
        finish=1;
        for(i=0;i<K_number;i++){
            difference=0.0;   
            for(j=0;j<CVQ_DIM;j++){
                difference +=
					abs(init_main_vector[i][j]-vqd.main_vector[i][j]);
            }
            if(difference>VQ_MINERR) finish=0;
        }
        if(finish == 1){ // 処理終了
            // 各ブロックを代表ベクトルに置換
            for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
				xs=VQ_BLKSZ * x;   
				ys=VQ_BLKSZ * y;
				for(i=0;i<VQ_BLKSZ;i++) for(j=0;j<VQ_BLKSZ;j++){
					if(ip_safe(xs+j,ys+i,iWidth,iHeight)){
						inb[XY1D((xs+j)*3,ys+i,iLength)+0] =
							vqd.main_vector[vqd.b_clstr[y][x]][XY1D(j*3,i,VQ_BLKSZ*3)+2];
						inb[XY1D((xs+j)*3,ys+i,iLength)+1] =
							vqd.main_vector[vqd.b_clstr[y][x]][XY1D(j*3,i,VQ_BLKSZ*3)+1];
						inb[XY1D((xs+j)*3,ys+i,iLength)+2] =
							vqd.main_vector[vqd.b_clstr[y][x]][XY1D(j*3,i,VQ_BLKSZ*3)+0];
					} else {
						inb[XY1D((xs+j)*3,ys+i,iLength)+0] = 0;
						inb[XY1D((xs+j)*3,ys+i,iLength)+1] = 0;
						inb[XY1D((xs+j)*3,ys+i,iLength)+2] = 0;
					}
                }
            }
        }
    }
    // データ量・圧縮率関係の計算と表示
    original_data_size = iSize * 8 * 3;
    i=K_number-1; j=0; // j:代表ベクトルの指定に必要なビット数
    do{ i>>=1; j++; } while(i > 0);
	// 圧縮後のサイズ
    coded_data_size = VQ_BLKSZ * VQ_BLKSZ * K_number * 8 * 3+
		j * vqd.x_block * vqd.y_block; // 代表ベクトル数
	// 圧縮率
    compression_rate=(double)coded_data_size /
		(double)original_data_size * 100.0;
	// 後片付け
	GlobalFree(vqd.cv);
	
	return compression_rate;
#endif
}