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