images/sig-seg.c

   1 /*
   2  *      Image Library -- Image segmentation
   3  *
   4  *      (c) 2006 Pavel Charvat <pchar@ucw.cz>
   5  *
   6  *      This software may be freely distributed and used according to the terms
   7  *      of the GNU Lesser General Public License.
   8  */
   9
  10 #undef LOCAL_DEBUG
  11
  12 #include "ucw/lib.h"
  13 #include "ucw/conf.h"
  14 #include "ucw/heap.h"
  15 #include "images/images.h"
  16 #include "images/signature.h"
  17 #include "images/math.h"
  18
  19 #include <string.h>
  20
  21 #ifdef LOCAL_DEBUG
  22 static void
  23 dump_segmentation(struct image_sig_region *regions, uns regions_count)
  24 {
  25   uns cols = 0, rows = 0;
  26   for (uns i = 0; i < regions_count; i++)
  27     for (struct image_sig_block *b = regions[i].blocks; b; b = b->next)
  28       {
  29         cols = MAX(cols, b->x + 1);
  30         rows = MAX(rows, b->y + 1);
  31       }
  32   uns size = (cols + 1) * rows;
  33   byte buf[size];
  34   bzero(buf, size);
  35   for (uns i = 0; i < regions_count; i++)
  36     {
  37       byte c = (i < 10) ? '0' + i : 'A' - 10 + i;
  38       for (struct image_sig_block *b = regions[i].blocks; b; b = b->next)
  39         buf[b->x + b->y * (cols + 1)] = c;
  40     }
  41   for (uns i = 0; i < rows; i++)
  42     log(L_DEBUG, "%s", &buf[i * (cols + 1)]);
  43 }
  44 #endif
  45
  46 /* Pre-quantization - recursively split groups of blocks with large error */
  47
  48 static inline void
  49 prequant_init_region(struct image_sig_region *region)
  50 {
  51   bzero(region, sizeof(*region));
  52 }
  53
  54 static inline void
  55 prequant_add_block(struct image_sig_region *region, struct image_sig_block *block)
  56 {
  57   block->next = region->blocks;
  58   region->blocks = block;
  59   region->count++;
  60   for (uns i = 0; i < IMAGE_VEC_F; i++)
  61     {
  62       region->b[i] += block->v[i];
  63       region->c[i] += isqr(block->v[i]);
  64     }
  65 }
  66
  67 static void
  68 prequant_finish_region(struct image_sig_region *region)
  69 {
  70   if (region->count < 2)
  71     {
  72       region->e = 0;
  73     }
  74   else
  75     {
  76       u64 a = 0;
  77       region->e = 0;
  78       for (uns i = 0; i < IMAGE_VEC_F; i++)
  79         {
  80           region->e += region->c[i];
  81           a += (u64)region->b[i] * region->b[i];
  82         }
  83       region->e -= a / region->count;
  84       DBG("Finished region %u", (uns)region->e / region->count);
  85     }
  86 }
  87
  88 static inline uns
  89 prequant_heap_cmp(struct image_sig_region *a, struct image_sig_region *b)
  90 {
  91   return a->e > b->e;
  92 }
  93
  94 #define ASORT_PREFIX(x) prequant_##x
  95 #define ASORT_KEY_TYPE uns
  96 #include "ucw/sorter/array-simple.h"
  97
  98 static uns
  99 prequant(struct image_sig_block *blocks, uns blocks_count, struct image_sig_region *regions)
 100 {
 101   DBG("Starting pre-quantization");
 102
 103   uns regions_count, heap_count, axis;
 104   struct image_sig_block *blocks_end = blocks + blocks_count, *block, *block2;
 105   struct image_sig_region *heap[IMAGE_REG_MAX + 1], *region, *region2;
 106
 107   /* Initialize single region with all blocks */
 108   regions_count = heap_count = 1;
 109   heap[1] = regions;
 110   prequant_init_region(regions);
 111   for (block = blocks; block != blocks_end; block++)
 112     prequant_add_block(regions, block);
 113   prequant_finish_region(regions);
 114
 115   /* Main cycle */
 116   while (regions_count < IMAGE_REG_MAX &&
 117       regions_count <= DARY_LEN(image_sig_prequant_thresholds) && heap_count)
 118     {
 119       region = heap[1];
 120       DBG("Step... regions_count=%u heap_count=%u region->count=%u, region->e=%u",
 121           regions_count, heap_count, region->count, (uns)region->e);
 122       if (region->count < 2 ||
 123           region->e < image_sig_prequant_thresholds[regions_count - 1] * blocks_count)
 124         {
 125           HEAP_DELMIN(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP);
 126           continue;
 127         }
 128
 129       /* Select axis to split - the one with maximum average quadratic error */
 130       axis = 0;
 131       u64 cov = (u64)region->count * region->c[0] - (u64)region->b[0] * region->b[0];
 132       for (uns i = 1; i < 6; i++)
 133         {
 134           uns j = (u64)region->count * region->c[i] - (u64)region->b[i] * region->b[i];
 135           if (j > cov)
 136             {
 137               axis = i;
 138               cov = j;
 139             }
 140         }
 141       DBG("Splitting axis %u with average quadratic error %u", axis, (uns)(cov / (region->count * region->count)));
 142
 143       /* Sort values on the split axis */
 144       uns val[256], cnt[256], cval;
 145       if (region->count > 64)
 146         {
 147           bzero(cnt, sizeof(cnt));
 148           for (block = region->blocks; block; block = block->next)
 149             cnt[block->v[axis]]++;
 150           cval = 0;
 151           for (uns i = 0; i < 256; i++)
 152             if (cnt[i])
 153               {
 154                 val[cval] = i;
 155                 cnt[cval] = cnt[i];
 156                 cval++;
 157               }
 158         }
 159       else
 160         {
 161           block = region->blocks;
 162           for (uns i = 0; i < region->count; i++, block = block->next)
 163             val[i] = block->v[axis];
 164           prequant_sort(val, region->count);
 165           cval = 1;
 166           cnt[0] = 1;
 167           for (uns i = 1; i < region->count; i++)
 168             if (val[i] == val[cval - 1])
 169               cnt[cval - 1]++;
 170             else
 171               {
 172                 val[cval] = val[i];
 173                 cnt[cval] = 1;
 174                 cval++;
 175               }
 176         }
 177
 178       /* Select split value - to minimize error */
 179       uns b1 = val[0] * cnt[0];
 180       uns c1 = isqr(val[0]) * cnt[0];
 181       uns b2 = region->b[axis] - b1;
 182       uns c2 = region->c[axis] - c1;
 183       uns i = cnt[0], j = region->count - cnt[0];
 184       u64 best_err = c1 - (u64)b1 * b1 / i + c2 - (u64)b2 * b2 / j;
 185       uns split_val = val[0];
 186       for (uns k = 1; k < cval - 1; k++)
 187         {
 188           uns b0 = val[k] * cnt[k];
 189           uns c0 = isqr(val[k]) * cnt[k];
 190           b1 += b0;
 191           b2 -= b0;
 192           c1 += c0;
 193           c2 -= c0;
 194           i += cnt[k];
 195           j -= cnt[k];
 196           u64 err = (u64)c1 - (u64)b1 * b1 / i + (u64)c2 - (u64)b2 * b2 / j;
 197           if (err < best_err)
 198             {
 199               best_err = err;
 200               split_val = val[k];
 201             }
 202         }
 203       DBG("split_val=%u best_err=%llu b[axis]=%u c[axis]=%u", split_val, (long long)best_err, region->b[axis], region->c[axis]);
 204
 205       /* Split region */
 206       block = region->blocks;
 207       region2 = regions + regions_count++;
 208       prequant_init_region(region);
 209       prequant_init_region(region2);
 210       while (block)
 211         {
 212           block2 = block->next;
 213           if (block->v[axis] <= split_val)
 214             prequant_add_block(region, block);
 215           else
 216             prequant_add_block(region2, block);
 217           block = block2;
 218         }
 219       prequant_finish_region(region);
 220       prequant_finish_region(region2);
 221       HEAP_INCREASE(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP, 1);
 222       heap[++heap_count] = region2;
 223       HEAP_INSERT(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP);
 224     }
 225
 226   DBG("Pre-quantized to %u regions", regions_count);
 227
 228   return regions_count;
 229 }
 230
 231 /* Post-quantization - run a few K-mean iterations to improve pre-quantized regions */
 232
 233 static uns
 234 postquant(struct image_sig_block *blocks, uns blocks_count, struct image_sig_region *regions, uns regions_count)
 235 {
 236   DBG("Starting post-quantization");
 237
 238   struct image_sig_block *blocks_end = blocks + blocks_count, *block;
 239   struct image_sig_region *regions_end = regions + regions_count, *region;
 240   uns error = 0, last_error;
 241
 242   /* Initialize regions and initial segmentation error */
 243   for (region = regions; region != regions_end; )
 244     {
 245       uns inv = 0xffffffffU / region->count;
 246       for (uns i = 0; i < IMAGE_VEC_F; i++)
 247         {
 248           region->a[i] = ((u64)region->b[i] * inv) >> 32;
 249           error += region->c[i] - region->a[i] * region->b[i];
 250         }
 251       region++;
 252     }
 253
 254   /* Convergation cycle */
 255   for (uns step = 0; step < image_sig_postquant_max_steps; step++)
 256     {
 257       DBG("Step...");
 258
 259       /* Clear regions */
 260       for (region = regions; region != regions_end; region++)
 261         {
 262           region->blocks = NULL;
 263           region->count = 0;
 264           bzero(region->b, sizeof(region->b));
 265           bzero(region->c, sizeof(region->c));
 266         }
 267
 268       /* Assign each block to its nearest pivot and accumulate region variables */
 269       for (block = blocks; block != blocks_end; block++)
 270         {
 271           struct image_sig_region *best_region = NULL;
 272           uns best_dist = ~0U;
 273           for (region = regions; region != regions_end; region++)
 274             {
 275               uns dist =
 276                 isqr(block->v[0] - region->a[0]) +
 277                 isqr(block->v[1] - region->a[1]) +
 278                 isqr(block->v[2] - region->a[2]) +
 279                 isqr(block->v[3] - region->a[3]) +
 280                 isqr(block->v[4] - region->a[4]) +
 281                 isqr(block->v[5] - region->a[5]);
 282               if (dist <= best_dist)
 283                 {
 284                   best_dist = dist;
 285                   best_region = region;
 286                 }
 287             }
 288           region = best_region;
 289           region->count++;
 290           block->next = region->blocks;
 291           region->blocks = block;
 292           for (uns i = 0; i < IMAGE_VEC_F; i++)
 293             {
 294               region->b[i] += block->v[i];
 295               region->c[i] += isqr(block->v[i]);
 296             }
 297         }
 298
 299       /* Finish regions, delete empty ones (should appear rarely), compute segmentation error */
 300       last_error = error;
 301       error = 0;
 302       for (region = regions; region != regions_end; )
 303         if (region->count)
 304           {
 305             uns inv = 0xffffffffU / region->count;
 306             for (uns i = 0; i < IMAGE_VEC_F; i++)
 307               {
 308                 region->a[i] = ((u64)region->b[i] * inv) >> 32;
 309                 error += region->c[i] - region->a[i] * region->b[i];
 310               }
 311             region++;
 312           }
 313         else
 314           {
 315             regions_end--;
 316             *region = *regions_end;
 317           }
 318
 319       DBG("last_error=%u error=%u", last_error, error);
 320
 321       /* Convergation criteria */
 322       if (step >= image_sig_postquant_min_steps)
 323         {
 324           if (error > last_error)
 325             break;
 326           u64 dif = last_error - error;
 327           if (dif * image_sig_postquant_threshold < (u64)last_error * 100)
 328             break;
 329         }
 330     }
 331
 332   DBG("Post-quantized to %u regions with average square error %u", regions_end - regions, error / blocks_count);
 333
 334   return regions_end - regions;
 335 }
 336
 337 void
 338 image_sig_segmentation(struct image_sig_data *data)
 339 {
 340   data->regions_count = prequant(data->blocks, data->blocks_count, data->regions);
 341 #ifdef LOCAL_DEBUG
 342   dump_segmentation(data->regions, data->regions_count);
 343 #endif
 344   data->regions_count = postquant(data->blocks, data->blocks_count, data->regions, data->regions_count);
 345 #ifdef LOCAL_DEBUG
 346   dump_segmentation(data->regions, data->regions_count);
 347 #endif
 348 }
 349