001: //-----------------------------------------------------------------
002: // clustering.cpp:
003: //    クラスタリング用関数
004: //               Last Update: <2004/12/22 22:08:50 A.Murakami>
005: //-----------------------------------------------------------------
006: #include "clustering.h"
007: //-----------------------------------------------------------------
008: // 画像中の白画素のデータ取得
009: //-----------------------------------------------------------------
010: int obtain_points_data(LPBYTE inb,ClutData& cltd)
011: {
012:     int x,y;
013:     cltd.num_pnts=0; // 取得ポイント数
014:     for(y=0;y<iHeight;y++) for(x=0;x<iWidth;x++){
015:         if(inb[x+y*iWidth] == iWHITE){
016:             if(cltd.num_pnts<MAX_PNTS-2){
017:                 cltd.num_pnts ++;
018:                 cltd.points[cltd.num_pnts-1].x=x;
019:                 cltd.points[cltd.num_pnts-1].y=y;
020:             }
021:         }
022:     }
023:     if(cltd.num_pnts<5){
024:         MessageBox(NULL,"点列が少な過ぎます. 画像を作り直して下さい",0,0);
025:         return FALSE;
026:     }
027:     if(cltd.num_pnts>=MAX_PNTS){
028:         MessageBox(NULL,"点列が多すぎます. 画像を作り直して下さい",0,0);
029:         return FALSE;
030:     }
031:     return TRUE;
032: }
033: //-----------------------------------------------------------------
034: // 取得した点データをランダムに入れ替える.
035: //-----------------------------------------------------------------
036: void shuffle_points_data(ClutData& cltd)
037: {
038:     int i,num1,num2;
039:     int sw_max;
040:     // 乱数の初期化
041:     init_random();
042:     // 点の総数の２倍の回数の入れ替え
043:     sw_max = cltd.num_pnts * 2;
044:     for(i=0;i<sw_max;i++){
045:         num1 = random_int(cltd.num_pnts);
046:         do { 
047:             num2 = random_int(cltd.num_pnts);
048:         } while (num2==num1);
049:         // No.num1 と No.num2 の点を入れ替え
050:         swap_point(cltd.points[num1],cltd.points[num2]);
051:     }
052: }
053: //-----------------------------------------------------------------
054: // [y][x]ブロックのベクトルとNo.n 代表ベクトル間の距離計算
055: //-----------------------------------------------------------------
056: double calc_vdistance(int y,int x,int n,VQData& vqd)
057: {
058:     int i; double sum=0.0;
059:     for(i=0;i<VQ_DIM;i++)
060:         sum += (vqd.vector[y][x][i]-vqd.main_vector[n][i]) *
061:             (vqd.vector[y][x][i]-vqd.main_vector[n][i]);
062:     return (sqrt(sum));
063: }
064: //-----------------------------------------------------------------
065: // 画像を 4x4 画素のブロックに分割してベクトル化
066: //-----------------------------------------------------------------
067: void obtain_vector_data(LPBYTE inb,VQData& vqd)
068: {
069:     int x,y,xs,ys,i,j;
070:     for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
071:         xs=VQ_BLKSZ * x; ys=VQ_BLKSZ * y;
072:         for(i=0;i<VQ_BLKSZ;i++) for(j=0;j<VQ_BLKSZ;j++){
073:             if(ip_safe(xs+j,ys+i,iWidth,iHeight))
074:                 vqd.vector[y][x][i*VQ_BLKSZ+j] =
075:                     inb[XY1D(xs+j,ys+i,iWidth)];
076:             else
077:                 vqd.vector[y][x][i*VQ_BLKSZ+j] = 0;
078:         }
079:     }
080: }
081: //-----------------------------------------------------------------
082: // [y][x]ブロックのベクトルとNo.n 代表ベクトル間の距離計算
083: //-----------------------------------------------------------------
084: double calc_c_vdistance(int y,int x,int n,CVQData& vqd)
085: {
086:     int i; double sum=0.0;
087:     for(i=0;i<CVQ_DIM;i++)
088:         sum += (vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[i]-vqd.main_vector[n][i]) *
089:             (vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[i]-vqd.main_vector[n][i]);
090:     return (sqrt(sum));
091: }
092: //-----------------------------------------------------------------
093: // 画像を 4x4 画素のブロックに分割してベクトル化
094: //-----------------------------------------------------------------
095: void obtain_c_vector_data(LPBYTE inb,CVQData& vqd)
096: {
097:     int x,y,xs,ys,i,j;
098:     for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
099:         xs=VQ_BLKSZ * x; ys=VQ_BLKSZ * y;
100:         for(i=0;i<VQ_BLKSZ;i++) for(j=0;j<VQ_BLKSZ;j++){
101:             if(ip_safe(xs+j,ys+i,iWidth,iHeight)){
102:                 vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[XY1D(j*3,i,VQ_BLKSZ*3)+0] =
103:                     inb[XY1D((xs+j)*3,ys+i,iLength)+2];
104:                 vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[XY1D(j*3,i,VQ_BLKSZ*3)+1] =
105:                     inb[XY1D((xs+j)*3,ys+i,iLength)+1];
106:                 vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[XY1D(j*3,i,VQ_BLKSZ*3)+2] =
107:                     inb[XY1D((xs+j)*3,ys+i,iLength)+0];
108:             } else {
109:                 vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[XY1D(j*3,i,VQ_BLKSZ*3)+0] = 0;
110:                 vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[XY1D(j*3,i,VQ_BLKSZ*3)+1] = 0;
111:                 vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[XY1D(j*3,i,VQ_BLKSZ*3)+2] = 0;
112:             }
113:         }
114:     }
115: }
116: //-----------------------------------------------------------------
117: // NN(Nearest Neighbour) 法によるクラスタリング
118: //-----------------------------------------------------------------
119: void clustering_NNmethod(LPBYTE inb,int clstr_r,int dopt)
120: {
121:     int i,j,min_num;
122:     double min_distance,distance;
123:     ClutData cltd;
124:     //--------------------------------------------------
125:     // 画像中の点を取得
126:     //--------------------------------------------------
127:     if(!obtain_points_data(inb,cltd)) return;
128:     //--------------------------------------------------
129:     // No.0 の点を No.0 のクラスタ代表点に設定
130:     //--------------------------------------------------
131:     cltd.num_clstr=0; // クラスタの総数 - 1
132:     cltd.cclstr[cltd.num_clstr].x = cltd.points[0].x;
133:     cltd.cclstr[cltd.num_clstr].y = cltd.points[0].y;
134:     cltd.p_clstr[0]=0;
135:     cltd.num_clstr ++;
136:     //--------------------------------------------------
137:     // No.1 から No.num_pnts-1 までの点を分類する
138:     //--------------------------------------------------
139:     for(i=1;i<cltd.num_pnts;i++){
140:         // 最も近いクラスタの計算
141:         //   No: min_num, 距離: min_distance
142:         min_distance=10000.0;
143:         for(j=0;j<cltd.num_clstr;j++){
144:             distance=calc_distance(cltd.points[i],cltd.cclstr[j]);
145:             if(distance<min_distance){
146:                 min_distance=distance;
147:                 min_num=j;
148:             }
149:         }
150:         // No.i の点をクラスタに加える
151:         if(min_distance<clstr_r){
152:             cltd.p_clstr[i]=min_num;
153:             // クラスタ重心の移動
154:             cltd.cclstr[min_num].x+=cltd.points[i].x;
155:             cltd.cclstr[min_num].x/=2;
156:             cltd.cclstr[min_num].y+=cltd.points[i].y;
157:             cltd.cclstr[min_num].y/=2;
158:         }
159:         // No.i の点は新しいクラスタに加える
160:         else {
161:             cltd.cclstr[cltd.num_clstr].x=cltd.points[i].x;
162:             cltd.cclstr[cltd.num_clstr].y=cltd.points[i].y;
163:             cltd.p_clstr[i]=cltd.num_clstr;
164:             cltd.num_clstr ++;
165:         }
166:     }
167:     //--------------------------------------------------
168:     // クラスタ位置の表示
169:     //--------------------------------------------------
170:     // 入力点と各点までの線分を描く
171:     if(dopt){
172:         for(i=0;i<cltd.num_pnts;i++){
173:             draw_circle(inb,cltd.points[i],5);
174:             if(dopt==CLUT_DALL)
175:                 draw_line(inb,cltd.points[i],cltd.cclstr[cltd.p_clstr[i]]);
176:         }
177:     }
178:     // クラスタ中心の表示
179:     for(i=0;i<cltd.num_clstr;i++){
180:         draw_circle(inb,cltd.cclstr[i],clstr_r);
181:     }
182: }
183: //-----------------------------------------------------------------
184: // 階層的クラスタリング: 樹形図的併合方法によるクラスタリング
185: //-----------------------------------------------------------------
186: //   1. パターン間の距離を計算
187: //   2. 最も近接する2つのクラスタを併合し新しいクラスタをつくる
188: //   3. 新しいクラスタと残りのクラスタとの距離を計算
189: //   4. 最小のクラスタ間距離が閾値より大きくなれば終了(2,3の反復)
190: //-----------------------------------------------------------------
191: void clustering_Tree(LPBYTE inb,int clstr_r,int dopt)
192: {
193:     int end_loop = 1000;                // 最大反復回数
194:     int min_num1,min_num2;              // 併合するクラスタ番号
195:     int radius_clstr[MAX_PNTS];         // クラスタの意味上の半径
196:     LPUINT dist_clstr;                  // クラスタの意味上の半径
197:     int i,j,dx,dy,counter,finish,sum_number;
198:     double min_distance,distance;
199:     double sum_x,sum_y;
200:     ClutData cltd;
201:     //--------------------------------------------------
202:     // 画像中の点を取得
203:     //--------------------------------------------------
204:     if(!obtain_points_data(inb,cltd)) return;
205:     dist_clstr = (LPUINT)GlobalAlloc(GPTR,sizeof(UINT)*MAX_PNTS*MAX_PNTS);
206:     //--------------------------------------------------
207:     // 初期のクラスタ中心の設定
208:     //--------------------------------------------------
209:     cltd.num_clstr=cltd.num_pnts;
210:     for(i=0;i<cltd.num_clstr;i++){
211:         cltd.cclstr[i].x=cltd.points[i].x;
212:         cltd.cclstr[i].y=cltd.points[i].y;
213:         radius_clstr[i]=0; // クラスタ半径
214:         cltd.p_clstr[i] =i;
215:     }
216:     // クラスタ間の距離を求める
217:     for(i=0;i<cltd.num_clstr;i++) for(j=i+1;j<cltd.num_clstr;j++){
218:         distance=calc_distance(cltd.cclstr[i],cltd.cclstr[j]);
219:         dist_clstr[j+i*cltd.num_clstr] = rint(distance);
220:     }
221:     //--------------------------------------------------
222:     // クラスタリング処理
223:     //--------------------------------------------------
224:     finish=0; counter=0;
225:     while(!finish){
226:         counter++;
227:         // クラスタ間の最小距離を求める
228:         min_distance=10000.0;
229:         for(i=0;i<cltd.num_clstr;i++) for(j=i+1;j<cltd.num_clstr;j++){
230:             if(dist_clstr[j+i*cltd.num_clstr]<min_distance){
231:                 min_distance=dist_clstr[j+i*cltd.num_clstr];
232:                 min_num1=i; min_num2=j;
233:             }
234:         }
235:         // 最小距離が許容範囲以下ならクラスタを併合
236:         if(min_distance < clstr_r){
237:             // 併合されたクラスタの除去
238:             for(i=min_num2;i<cltd.num_clstr-1;i++){
239:                 cltd.cclstr[i].x = cltd.cclstr[i+1].x;
240:                 cltd.cclstr[i].y = cltd.cclstr[i+1].y;
241:             }
242:             // パターンの所属クラスタの変更
243:             for(i=0;i<cltd.num_pnts;i++){
244:                 if(cltd.p_clstr[i]>min_num2)
245:                     cltd.p_clstr[i]--;
246:                 else if(cltd.p_clstr[i]==min_num2)
247:                     cltd.p_clstr[i]=min_num1;
248:             }
249:             // 新規クラスタ重心の設定
250:             sum_number = 0;
251:             sum_x = sum_y = 0;
252:             for(i=0;i<cltd.num_pnts;i++){
253:                 if(cltd.p_clstr[i]==min_num1){
254:                     sum_number++;
255:                     sum_x += cltd.points[i].x; sum_y += cltd.points[i].y;
256:                 }
257:             }
258:             cltd.cclstr[min_num1].x=rint(sum_x/(double)sum_number);
259:             cltd.cclstr[min_num1].y=rint(sum_y/(double)sum_number);
260:             // 併合された距離配列の削除
261:             for(i=0;i<cltd.num_clstr-1;i++){
262:                 for(j=i+1;j<cltd.num_clstr-1;j++){
263:                     dx = (j>=min_num2)?1:0;
264:                     dy = (i>=min_num2)?1:0;
265:                     dist_clstr[j+i*(cltd.num_clstr-1)] =
266:                         dist_clstr[(j+dx)+(i+dy)*cltd.num_clstr];
267:                 }
268:             }
269:             // クラスタ数を減らす
270:             cltd.num_clstr--;
271:             // 距離の再計算
272:             for(i=0;i<=min_num1;i++) for(j=i+1;j<cltd.num_clstr;j++){
273:                 if(i==min_num1 || j==min_num1){
274:                     distance=calc_distance(cltd.cclstr[i],
275:                                            cltd.cclstr[j]);
276:                     dist_clstr[j+i*cltd.num_clstr] = rint(distance);
277:                 }
278:             }
279:         } else finish = 1;
280:         if(counter > end_loop) finish = 1;
281:     }
282:     //--------------------------------------------------
283:     // 各クラスタの意味上の半径を求める
284:     //--------------------------------------------------
285:     for(i=0;i<cltd.num_clstr;i++) for(j=0;j<cltd.num_pnts;j++){
286:         if(cltd.p_clstr[j] == i){
287:             distance=calc_distance(cltd.cclstr[i],cltd.points[j]);
288:             if(distance > radius_clstr[i])
289:                 radius_clstr[i]=rint(distance);
290:         }
291:     }
292:     //--------------------------------------------------
293:     // クラスタ位置の表示
294:     //--------------------------------------------------
295:     // 入力点と各点データまでの線分を描く
296:     if(dopt){
297:         for(i=0;i<cltd.num_pnts;i++){
298:             draw_circle(inb,cltd.points[i],5);
299:             if(dopt==CLUT_DALL)
300:                 draw_line(inb,cltd.points[i],
301:                           cltd.cclstr[cltd.p_clstr[i]]);
302:         }
303:     }
304:     // クラスタ中心の表示
305:     for(i=0;i<cltd.num_clstr;i++){
306:         draw_circle(inb,cltd.cclstr[i],radius_clstr[i]);
307:     }
308:     // 後片付け
309:     GlobalFree(dist_clstr);
310: }
311: //-----------------------------------------------------------------
312: // K-平均法によるクラスタリング
313: //-----------------------------------------------------------------
314: //   1. K個のクラスタの初期中心位置を与える(ランダム)
315: //   2. 各点を近接するクラスタに分類
316: //   3. 各クラスタに属する点の中心位置の計算
317: //   4. 初期中心と比較し、一致・許容範囲内なら終了(2,3の反復)
318: //-----------------------------------------------------------------
319: void clustering_Kmean(LPBYTE inb,int clstr_r,int K_number,int dopt)
320: {
321:     int i,j,min_num,sum_number;
322:     int counter,finish;
323:     double min_distance,distance,sum_x,sum_y;
324:     int radius_clstr[MAX_PNTS];       // クラスタの意味上の半径
325:     POINT init_cntr_clstr[MAX_PNTS];  // クラスタの初期位置
326:     ClutData cltd;
327:     //--------------------------------------------------
328:     if(K_number>MAX_CLSTR) return;
329:     //--------------------------------------------------
330:     // 画像中の点を取得
331:     //--------------------------------------------------
332:     if(!obtain_points_data(inb,cltd)) return;
333:     //--------------------------------------------------
334:     // 初期のクラスタ中心の設定
335:     //--------------------------------------------------
336:     cltd.num_clstr=K_number;
337:     shuffle_points_data(cltd); // ランダムに入れ替え
338:     for(i=0;i<cltd.num_clstr;i++){
339:         cltd.cclstr[i].x=cltd.points[i].x;
340:         cltd.cclstr[i].y=cltd.points[i].y;
341:     }
342:     //--------------------------------------------------
343:     // クラスタリング処理
344:     //--------------------------------------------------
345:     finish=0; counter=0;
346:     while(!finish){
347:         counter++;
348:         // クラスタの初期代表点への代入 
349:         for(i=0;i<cltd.num_clstr;i++){
350:             init_cntr_clstr[i].x=cltd.cclstr[i].x;
351:             init_cntr_clstr[i].y=cltd.cclstr[i].y;
352:             radius_clstr[i]=clstr_r; // defalt 値
353:         }
354:         // 各点の所属クラスタ番号の初期化
355:         for(i=0;i<cltd.num_pnts;i++)
356:             cltd.p_clstr[i]=-1;
357:         // 各点の所属クラスタ
358:         for(i=0;i<cltd.num_pnts;i++){
359:             // 各クラスタ代表点との距離
360:             min_distance=10000.0;
361:             for(j=0;j<cltd.num_clstr;j++){
362:                 distance=calc_distance(cltd.points[i],cltd.cclstr[j]);
363:                 if(distance<min_distance){
364:                     min_distance=distance;
365:                     min_num=j;
366:                 }
367:             }
368:             // No.i の点をクラスタ No.min_num に加える
369:             cltd.p_clstr[i]=min_num;
370:         }
371:         // 各クラスタの代表点を修正
372:         for(i=0;i<cltd.num_clstr;i++){
373:             sum_number=0;
374:             sum_x=sum_y=0.0;
375:             for(j=0;j<cltd.num_pnts;j++){
376:                 if(cltd.p_clstr[j] == i){
377:                     sum_number++;
378:                     sum_x=sum_x + cltd.points[j].x;
379:                     sum_y=sum_y + cltd.points[j].y;
380:                 }
381:             }
382:             cltd.cclstr[i].x=rint(sum_x/(double)sum_number);
383:             cltd.cclstr[i].y=rint(sum_y/(double)sum_number);
384:         }
385:         // 各クラスタの意味上の半径を求める
386:         for(i=0;i<cltd.num_clstr;i++){
387:             for(j=0;j<cltd.num_pnts;j++){
388:                 if(cltd.p_clstr[j] == i){
389:                     distance=calc_distance(cltd.cclstr[i],cltd.points[j]);
390:                     if(distance > radius_clstr[i])
391:                         radius_clstr[i]=rint(distance);
392:                 }
393:             }
394:         }
395:         // クラスタ中心が移動したかどうかのチェック
396:         finish=1;
397:         for(i=0;i<cltd.num_clstr;i++){
398:             distance=calc_distance(init_cntr_clstr[i],cltd.cclstr[i]);
399:             if(distance > KC_MINERR) finish=0; // 許容誤差
400:         }
401:     }
402:     //--------------------------------------------------
403:     // クラスタ位置の表示
404:     //--------------------------------------------------
405:     // 入力点と各点データまでの線分を描く
406:     if(dopt){
407:         for(i=0;i<cltd.num_pnts;i++){
408:             draw_circle(inb,cltd.points[i],5);
409:             if(dopt==CLUT_DALL)
410:                 draw_line(inb,cltd.points[i],cltd.cclstr[cltd.p_clstr[i]]);
411:         }
412:     }
413:     // クラスタ中心の表示
414:     for(i=0;i<cltd.num_clstr;i++){
415:         draw_circle(inb,cltd.cclstr[i],radius_clstr[i]);
416:     }
417: }
418: 
419: //-----------------------------------------------------------------
420: // ベクトル量子化[K平均法によるクラスタリング]: 濃淡画像
421: //-----------------------------------------------------------------
422: double vector_quantize(LPBYTE inb,int K_number)
423: // int K_number: K値(生成するクラスタの総数),default=100
424: // 返り値: 圧縮率[圧縮後のサイズ/元サイズ]
425: {
426:     int i,j,x,y,xs,ys;
427:     int counter,min_num,finish;
428:     int original_data_size,coded_data_size;
429:     double min_distance,distance,compression_rate;
430:     double sum[VQ_DIM],sum_number,difference;
431:     BYTE init_main_vector[MAX_CLSTR][VQ_DIM]; // 代表ベクトル
432:     VQData vqd;
433:     //-----------------------------------------------------------------
434:     if(K_number>MAX_CLSTR) return 0;
435:     //-----------------------------------------------------------------
436:     // ブロック数
437:     //-----------------------------------------------------------------
438:     vqd.x_block = iWidth  / VQ_BLKSZ;
439:     vqd.y_block = iHeight / VQ_BLKSZ;
440:     if(iWidth  % VQ_BLKSZ != 0) vqd.x_block ++;
441:     if(iHeight % VQ_BLKSZ != 0) vqd.y_block ++;
442:     //-----------------------------------------------------------------
443:     // ベクトルデータの取得
444:     //-----------------------------------------------------------------
445:     obtain_vector_data(inb,vqd);
446:     //-----------------------------------------------------------------
447:     // 仮代表ベクトルの設定
448:     //-----------------------------------------------------------------
449:     init_random(); // 乱数の初期化
450:     for(i=0;i<K_number;i++){
451:         x=random_int(vqd.x_block);
452:         y=random_int(vqd.y_block);
453:         for(j=0;j<VQ_DIM;j++){
454:             vqd.main_vector[i][j] = vqd.vector[y][x][j];
455:         }
456:     }
457:     //-----------------------------------------------------------------
458:     // K-平均法によるベクトルのクラスタリング
459:     //-----------------------------------------------------------------
460:     finish=0; counter=0;
461:     while(!finish){
462:         counter++;
463:         //--------------------------------------------------
464:         // クラスタの初期代表点への代入
465:         //--------------------------------------------------
466:         for(i=0;i<K_number;i++) for(j=0;j<VQ_DIM;j++)
467:             init_main_vector[i][j] = vqd.main_vector[i][j];
468:         //--------------------------------------------------
469:         // 各ブロックの所属クラスタ番号の初期化
470:         //--------------------------------------------------
471:         for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++)
472:             vqd.b_clstr[y][x] = -1;
473:         //--------------------------------------------------
474:         // 各点の所属クラスタを決定
475:         //--------------------------------------------------
476:         for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
477:             // 各クラスタ代表点との距離
478:             min_distance=iWHITE * iWHITE * VQ_DIM;
479:             for(i=0;i<K_number;i++){
480:                 distance=calc_vdistance(y,x,i,vqd);
481:                 if(distance<min_distance){
482:                     min_distance = distance;
483:                     min_num = i;
484:                 }
485:             }
486:             // [y][x]ブロックをクラスタ No.min_num へ
487:             vqd.b_clstr[y][x] = min_num;
488:         }
489:         //--------------------------------------------------
490:         // 各クラスタの代表点を修正する
491:         //--------------------------------------------------
492:         for(i=0;i<K_number;i++){
493:             for(j=0;j<VQ_DIM;j++) sum[j]=0.0;
494:             sum_number=0.0;
495:             for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
496:                 if(vqd.b_clstr[y][x]==i){
497:                     sum_number = sum_number + 1.0;
498:                     for(j=0;j<VQ_DIM;j++)
499:                         sum[j] += vqd.vector[y][x][j];
500:                 }
501:             }
502:             for(j=0;j<VQ_DIM;j++){
503:                 vqd.main_vector[i][j]=(BYTE)(sum[j] / sum_number);
504:             }
505:         }
506:         //--------------------------------------------------
507:         // クラスタ中心が移動したかどうかのチェック
508:         //--------------------------------------------------
509:         finish=1;
510:         for(i=0;i<K_number;i++){
511:             difference=0.0;   
512:             for(j=0;j<VQ_DIM;j++){
513:                 difference +=
514:                     abs(init_main_vector[i][j]-vqd.main_vector[i][j]);
515:             }
516:             if(difference>VQ_MINERR) finish=0;
517:         }
518:         if(finish == 1){ // 処理終了
519:             // 各ブロックを代表ベクトルに置換
520:             for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
521:                 xs=VQ_BLKSZ * x;   
522:                 ys=VQ_BLKSZ * y;
523:                 for(i=0;i<VQ_BLKSZ;i++) for(j=0;j<VQ_BLKSZ;j++){
524:                     if(ip_safe(xs+j,ys+i,iWidth,iHeight))
525:                         inb[XY1D(xs+j,ys+i,iWidth)] =
526:                             vqd.main_vector[vqd.b_clstr[y][x]][i*VQ_BLKSZ+j];
527:                     else
528:                         inb[XY1D(xs+j,ys+i,iWidth)] = 0;
529:                 }
530:             }
531:         }
532:     }
533:     // データ量・圧縮率関係の計算と表示
534:     original_data_size = iSize * 8;
535:     i=K_number-1; j=0; // j:代表ベクトルの指定に必要なビット数
536:     do{ i>>=1; j++; } while(i > 0);
537:     // 圧縮後のサイズ
538:     coded_data_size = VQ_BLKSZ * VQ_BLKSZ * K_number * 8 +
539:         j * vqd.x_block * vqd.y_block; // 代表ベクトル数
540:     // 圧縮率
541:     compression_rate=(double)coded_data_size /
542:         (double)original_data_size * 100.0;
543:     
544:     return compression_rate;
545: }
546: 
547: //-----------------------------------------------------------------
548: // ベクトル量子化[K平均法によるクラスタリング]: カラー画像
549: //-----------------------------------------------------------------
550: double c_vector_quantize(LPBYTE inb,int K_number)
551: // int K_number: K値(生成するクラスタの総数),default=100
552: // 返り値: 圧縮率[圧縮後のサイズ/元サイズ]
553: {
554: #if 0 // 各RGB成分ごと
555:     double comp_rate;
556:     LPBYTE iR=GetGray(COLOR_R);
557:     LPBYTE iG=GetGray(COLOR_G);
558:     LPBYTE iB=GetGray(COLOR_B);
559:     comp_rate = 0;
560:     comp_rate += vector_quantize(iR,K_number);
561:     comp_rate += vector_quantize(iG,K_number);
562:     comp_rate += vector_quantize(iB,K_number);
563:     comp_rate /= 3.0;
564:     RGBToColor(inb,iR,iG,iB);
565:     // 後片付け
566:     GlobalFree(iR); GlobalFree(iG); GlobalFree(iB);
567:     return comp_rate;
568: #else // RGB成分をベクトルとして扱う
569:     int i,j,x,y,xs,ys;
570:     int counter,min_num,finish;
571:     int original_data_size,coded_data_size;
572:     double min_distance,distance,compression_rate;
573:     double sum[CVQ_DIM],sum_number,difference;
574:     BYTE init_main_vector[MAX_CLSTR][CVQ_DIM]; // 代表ベクトル
575:     CVQData vqd;
576:     //-----------------------------------------------------------------
577:     if(K_number>MAX_CLSTR) return 0;
578:     //-----------------------------------------------------------------
579:     // ブロック数
580:     //-----------------------------------------------------------------
581:     vqd.x_block = iWidth  / VQ_BLKSZ;
582:     vqd.y_block = iHeight / VQ_BLKSZ;
583:     if(iWidth  % VQ_BLKSZ != 0) vqd.x_block ++;
584:     if(iHeight % VQ_BLKSZ != 0) vqd.y_block ++;
585:     //-----------------------------------------------------------------
586:     // ベクトルデータの取得
587:     //-----------------------------------------------------------------
588:     vqd.cv = (C_vector*)GlobalAlloc(GPTR,sizeof(C_vector)*VQ_MAXBLK*VQ_MAXBLK);
589:     obtain_c_vector_data(inb,vqd);
590:     //-----------------------------------------------------------------
591:     // 仮代表ベクトルの設定
592:     //-----------------------------------------------------------------
593:     init_random(); // 乱数の初期化
594:     for(i=0;i<K_number;i++){
595:         x=random_int(vqd.x_block);
596:         y=random_int(vqd.y_block);
597:         for(j=0;j<CVQ_DIM;j++){
598:             vqd.main_vector[i][j] = vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[j];
599:         }
600:     }
601:     //-----------------------------------------------------------------
602:     // K-平均法によるベクトルのクラスタリング
603:     //-----------------------------------------------------------------
604:     finish=0; counter=0;
605:     while(!finish){
606:         counter++;
607:         //--------------------------------------------------
608:         // クラスタの初期代表点への代入
609:         //--------------------------------------------------
610:         for(i=0;i<K_number;i++) for(j=0;j<CVQ_DIM;j++)
611:             init_main_vector[i][j] = vqd.main_vector[i][j];
612:         //--------------------------------------------------
613:         // 各ブロックの所属クラスタ番号の初期化
614:         //--------------------------------------------------
615:         for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++)
616:             vqd.b_clstr[y][x] = -1;
617:         //--------------------------------------------------
618:         // 各点の所属クラスタを決定
619:         //--------------------------------------------------
620:         for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
621:             // 各クラスタ代表点との距離
622:             min_distance=iWHITE * iWHITE * CVQ_DIM;
623:             for(i=0;i<K_number;i++){
624:                 distance=calc_c_vdistance(y,x,i,vqd);
625:                 if(distance<min_distance){
626:                     min_distance = distance;
627:                     min_num = i;
628:                 }
629:             }
630:             // [y][x]ブロックをクラスタ No.min_num へ
631:             vqd.b_clstr[y][x] = min_num;
632:         }
633:         //--------------------------------------------------
634:         // 各クラスタの代表点を修正する
635:         //--------------------------------------------------
636:         for(i=0;i<K_number;i++){
637:             for(j=0;j<CVQ_DIM;j++) sum[j]=0.0;
638:             sum_number=0.0;
639:             for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
640:                 if(vqd.b_clstr[y][x]==i){
641:                     sum_number = sum_number + 1.0;
642:                     for(j=0;j<CVQ_DIM;j++)
643:                         sum[j] += vqd.cv[XY1D(x,y,VQ_MAXBLK)].d[j];
644:                 }
645:             }
646:             for(j=0;j<CVQ_DIM;j++){
647:                 vqd.main_vector[i][j]=(BYTE)(sum[j] / sum_number);
648:             }
649:         }
650:         //--------------------------------------------------
651:         // クラスタ中心が移動したかどうかのチェック
652:         //--------------------------------------------------
653:         finish=1;
654:         for(i=0;i<K_number;i++){
655:             difference=0.0;   
656:             for(j=0;j<CVQ_DIM;j++){
657:                 difference +=
658:                     abs(init_main_vector[i][j]-vqd.main_vector[i][j]);
659:             }
660:             if(difference>VQ_MINERR) finish=0;
661:         }
662:         if(finish == 1){ // 処理終了
663:             // 各ブロックを代表ベクトルに置換
664:             for(y=0;y<vqd.y_block;y++) for(x=0;x<vqd.x_block;x++){
665:                 xs=VQ_BLKSZ * x;   
666:                 ys=VQ_BLKSZ * y;
667:                 for(i=0;i<VQ_BLKSZ;i++) for(j=0;j<VQ_BLKSZ;j++){
668:                     if(ip_safe(xs+j,ys+i,iWidth,iHeight)){
669:                         inb[XY1D((xs+j)*3,ys+i,iLength)+0] =
670:                             vqd.main_vector[vqd.b_clstr[y][x]][XY1D(j*3,i,VQ_BLKSZ*3)+2];
671:                         inb[XY1D((xs+j)*3,ys+i,iLength)+1] =
672:                             vqd.main_vector[vqd.b_clstr[y][x]][XY1D(j*3,i,VQ_BLKSZ*3)+1];
673:                         inb[XY1D((xs+j)*3,ys+i,iLength)+2] =
674:                             vqd.main_vector[vqd.b_clstr[y][x]][XY1D(j*3,i,VQ_BLKSZ*3)+0];
675:                     } else {
676:                         inb[XY1D((xs+j)*3,ys+i,iLength)+0] = 0;
677:                         inb[XY1D((xs+j)*3,ys+i,iLength)+1] = 0;
678:                         inb[XY1D((xs+j)*3,ys+i,iLength)+2] = 0;
679:                     }
680:                 }
681:             }
682:         }
683:     }
684:     // データ量・圧縮率関係の計算と表示
685:     original_data_size = iSize * 8 * 3;
686:     i=K_number-1; j=0; // j:代表ベクトルの指定に必要なビット数
687:     do{ i>>=1; j++; } while(i > 0);
688:     // 圧縮後のサイズ
689:     coded_data_size = VQ_BLKSZ * VQ_BLKSZ * K_number * 8 * 3+
690:         j * vqd.x_block * vqd.y_block; // 代表ベクトル数
691:     // 圧縮率
692:     compression_rate=(double)coded_data_size /
693:         (double)original_data_size * 100.0;
694:     // 後片付け
695:     GlobalFree(vqd.cv);
696:     
697:     return compression_rate;
698: #endif
699: }