images/k_means.c

   1 #include <stdlib.h>
   2 #include <string.h>
   3 #include <assert.h>
   4 #include <alloca.h>
   5 #include "lib/config.h"
   6 #include "img.h"
   7
   8 #define stack_alloc alloca
   9 #define xdiff(x, y) ( ((x)>(y)) ? (x)-(y) : (y)-(x) )
  10
  11 //void *memset(void *s, int c, size_t n);
  12
  13 uns
  14 bi_dist(struct BlockInfo *bi, struct ClassInfo *ci){
  15   return xdiff(bi->l, ci->l) +
  16          xdiff(bi->u, ci->u) +
  17          xdiff(bi->v, ci->v) +
  18          xdiff(bi->lh, ci->lh) +
  19          xdiff(bi->hl, ci->hl) +
  20          xdiff(bi->hh, ci->hh);
  21 }
  22
  23
  24 /*"simulated annealing norm" - hodnocenní úspì¹nosti simulovaného ¾íhání*/
  25 uns
  26 sa_norm(uns bi_len, struct BlockInfo *bi, u8 cls_num, struct ClassInfo *ci){
  27   uns ret=0;
  28   struct BlockInfo *pbi;
  29
  30   for(pbi=bi; pbi < bi+bi_len; pbi++)
  31     ret+=bi_dist(pbi, &ci[pbi->cls_num]);
  32   return ret;
  33 }
  34
  35 void
  36 showBlockInfoAsImage(struct BlockInfo *bi, uns bi_len, uns width, uns height, Image **__image,
  37                      ExceptionInfo *exception){
  38 /*  ImageInfo *image_info=CloneImageInfo((ImageInfo *) NULL);*/
  39 /*  Image *image=AllocateImage(image_info);*/
  40   unsigned char *p, *q;
  41   uns x, y;
  42
  43   width=(width>>2);
  44   height=(height>>2);
  45   assert(bi_len==width*height);
  46   p=(unsigned char*) malloc(3*(width<<2)*(height<<2)*sizeof(unsigned char));
  47   for (q=p, y=0; y < (height<<2); y++){
  48     for (x=0; x < (width<<2); x++){
  49       uns index=(y>>2)*width + (x>>2);
  50       assert(index<bi_len);
  51       *q++=(255/3) * (bi[index].cls_num>>2);
  52       *q++=(255) * ((bi[index].cls_num>>1)&0x1);
  53       *q++=(255) * ((bi[index].cls_num>>0)&0x1);
  54     }
  55   }
  56   *__image=ConstituteImage(width<<2, height<<2, "RGB", CharPixel, p, exception);
  57 }
  58
  59 void
  60 init_cls(uns cls_num, uns bi_len, u8 *cls){
  61   u8 *p;
  62   for(p=cls; p<cls+bi_len; p++)
  63     *p = (u8) (((double)cls_num)*rand()/(RAND_MAX+1.0));
  64 }
  65
  66 void
  67 eval_centers(uns bi_len, struct BlockInfo *bi, u8 cls_num, struct ClassInfo *ci){
  68   struct BlockInfo *pbi;
  69   struct ClassInfo *pci;
  70
  71   memset((void*) ci, 0, cls_num*sizeof(struct ClassInfo));
  72   for(pbi=bi; pbi<bi+bi_len; pbi++){
  73     struct ClassInfo *pcls=&ci[pbi->cls_num];
  74     pcls->count++;
  75     pcls->l+=pbi->l;
  76     pcls->u+=pbi->u;
  77     pcls->v+=pbi->v;
  78     pcls->lh+=pbi->lh;
  79     pcls->hl+=pbi->hl;
  80     pcls->hh+=pbi->hh;
  81   }
  82
  83   for(pci=ci; pci<ci+cls_num; pci++){
  84     uns count=pci->count;
  85     if(! count) continue;
  86     pci->l  /=count;
  87     pci->u  /=count;
  88     pci->v  /=count;
  89     pci->lh /=count;
  90     pci->hl /=count;
  91     pci->hh /=count;
  92   }
  93 }
  94
  95 #ifdef DEBUG_KMEANS
  96 void
  97 write_BlockInfo(uns len, struct BlockInfo *bi){
  98   struct BlockInfo *pbi;
  99   for(pbi=bi; pbi<bi+len; pbi++){
 100     fprintf(stderr, "(%u, %u, %u, %u, %u, %u; '%u')\n", pbi->l, pbi->u, pbi->v, pbi->lh, pbi->hl, pbi->hh, pbi->cls_num);
 101   }
 102 }
 103
 104 void
 105 write_ClassInfo(uns len, struct ClassInfo *ci){
 106   struct ClassInfo *pci;
 107   for(pci=ci; pci<ci+len; pci++){
 108     fprintf(stderr, "(%u, %u, %u, %u, %u, %u; %u)\n", pci->l, pci->u, pci->v, pci->lh, pci->hl, pci->hh, pci->count);
 109   }
 110 }
 111 #endif
 112
 113 /*bi does not mean an array, but pointer to block we search the nearest Class for*/
 114 u8
 115 nearestClass(uns cls_num, struct ClassInfo *ci, struct BlockInfo *bi){
 116   u8 ret;
 117   struct ClassInfo *pci;
 118   uns min_dist=bi_dist(bi, ci);
 119   ret=0;
 120
 121   for(pci=ci+1; pci<ci+cls_num; pci++){
 122     uns dist=bi_dist(bi, pci);
 123     if(dist<min_dist){
 124       min_dist=dist;
 125       ret = (u8) (pci-ci);
 126     }
 127   }
 128   return ret;
 129 }
 130
 131 uns
 132 __k_means(uns cls_num, struct DecomposeImageInfo* dii, struct BlockInfo* bi){
 133   struct ClassInfo *ci=(struct ClassInfo*) stack_alloc(cls_num*sizeof(struct ClassInfo));
 134   uns ret=~0U;
 135   uns cycle=0;
 136   eval_centers(dii->bi_len, bi, cls_num, ci);
 137   for(cycle=0; cycle<dii->max_cycles; cycle++){
 138     struct BlockInfo* pbi;
 139     for(pbi=bi; pbi<bi+dii->bi_len; pbi++){
 140       pbi->cls_num=nearestClass(cls_num, ci, pbi);
 141     }
 142     eval_centers(dii->bi_len, bi, cls_num, ci);
 143     if((ret=sa_norm(dii->bi_len, bi, cls_num, ci))<dii->threshold)
 144       break;
 145   }
 146   return ret;
 147 }
 148
 149 void
 150 BlockInfoTou8(struct BlockInfo *bi, u8 *cls, uns len){
 151   struct BlockInfo* pbi;
 152   u8* pcls;
 153
 154   for(pbi=bi, pcls=cls; pbi<bi+len; pbi++)
 155     *pcls++ = pbi->cls_num;
 156 }
 157
 158 void
 159 u8ToBlockInfo(u8 *cls, struct BlockInfo *bi, uns len){
 160   struct BlockInfo* pbi;
 161   u8* pcls;
 162
 163   for(pbi=bi, pcls=cls; pbi<bi+len; pbi++)
 164     pbi->cls_num = *pcls++;
 165 }
 166
 167 uns
 168 k_means(uns cls_num, struct BlockInfo *bi, struct DecomposeImageInfo* dii){
 169   u8 *act_cls=(u8*) stack_alloc(dii->bi_len*sizeof(u8));
 170   u8 *best_cls=(u8*) stack_alloc(dii->bi_len*sizeof(u8));
 171   uns best_diff=~0U, act_diff;
 172   uns i;
 173
 174   for(i=0; i<dii->init_decomp_num; i++){
 175     if(i)
 176       init_cls(cls_num, dii->bi_len, act_cls);
 177     else{
 178     /*usually, the decompozition, that imply from here is not the best, but ensures, that the
 179       return values from this fuction forms nonincreasing sequence for increasing cls_num*/
 180       BlockInfoTou8(bi, act_cls, dii->bi_len);
 181       bi[(uns) (dii->bi_len*rand()/(RAND_MAX+1.0))].cls_num=cls_num-1;
 182     }
 183
 184
 185     u8ToBlockInfo(act_cls, bi, dii->bi_len);
 186     act_diff=__k_means(cls_num, dii, bi);
 187     if(act_diff<best_diff){
 188       best_diff=act_diff;
 189       BlockInfoTou8(bi, best_cls, dii->bi_len);
 190     }
 191     /*fprintf(stderr, "...'%u'\n", act_diff);*/
 192   }
 193   u8ToBlockInfo(best_cls, bi, dii->bi_len);
 194   return best_diff;
 195 }
 196
 197 /*
 198   return: final number of classes in decomposition
 199 */
 200 uns
 201 decomposeImage(struct DecomposeImageInfo* dii, struct BlockInfo *bi){
 202   uns cls_num=1;
 203   uns err=0, old_err=0;
 204
 205   {
 206     struct ClassInfo ci;
 207     struct BlockInfo *pbi;
 208     for(pbi=bi; pbi<bi+dii->bi_len; pbi++)
 209       pbi->cls_num=0;
 210     eval_centers(dii->bi_len, bi, 1, &ci);
 211     old_err=sa_norm(dii->bi_len, bi, 1, &ci);
 212   }
 213   dii->threshold=old_err>>6;
 214   dii->diff_threshold=old_err>>8;
 215   fprintf(stderr, "\n%u; --\n", old_err);
 216   if(old_err<dii->threshold)
 217     return old_err;
 218   cls_num=1;
 219   while(cls_num<dii->max_cls_num){
 220     cls_num++;
 221     err=k_means(cls_num, bi, dii);
 222     fprintf(stderr, "\n(%u); %d\n", err, err-old_err);
 223     if(err<dii->threshold) break;
 224     if(err<old_err && err+dii->diff_threshold>old_err) break;
 225     old_err=err;
 226   }
 227   return err;
 228 }
 229
 230
 231
 232