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