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 + 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 #define ASORT_ELT(i) val[i]
  97 #define ASORT_EXTRA_ARGS , uns *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;
 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] * blocks_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 average quadratic error */
 132       axis = 0;
 133       u64 cov = (u64)region->count * region->c[0] - (u64)region->b[0] * region->b[0];
 134       for (uns i = 1; i < 6; i++)
 135         {
 136           uns j = (u64)region->count * region->c[i] - (u64)region->b[i] * region->b[i];
 137           if (j > cov)
 138             {
 139               axis = i;
 140               cov = j;
 141             }
 142         }
 143       DBG("Splitting axis %u with average quadratic error %u", axis, (uns)(cov / (region->count * region->count)));
 144
 145       /* Sort values on the split axis */
 146       uns val[256], cnt[256], cval;
 147       if (region->count > 64)
 148         {
 149           bzero(cnt, sizeof(cnt));
 150           for (block = region->blocks; block; block = block->next)
 151             cnt[block->v[axis]]++;
 152           cval = 0;
 153           for (uns i = 0; i < 256; i++)
 154             if (cnt[i])
 155               {
 156                 val[cval] = i;
 157                 cnt[cval] = cnt[i];
 158                 cval++;
 159               }
 160         }
 161       else
 162         {
 163           block = region->blocks;
 164           for (uns i = 0; i < region->count; i++, block = block->next)
 165             val[i] = block->v[axis];
 166           prequant_sort(region->count, val);
 167           cval = 1;
 168           cnt[0] = 1;
 169           for (uns i = 1; i < region->count; i++)
 170             if (val[i] == val[cval - 1])
 171               cnt[cval - 1]++;
 172             else
 173               {
 174                 val[cval] = val[i];
 175                 cnt[cval] = 1;
 176                 cval++;
 177               }
 178         }
 179
 180       /* Select split value - to minimize error */
 181       uns b1 = val[0] * cnt[0];
 182       uns c1 = isqr(val[0]) * cnt[0];
 183       uns b2 = region->b[axis] - b1;
 184       uns c2 = region->c[axis] - c1;
 185       uns i = cnt[0], j = region->count - cnt[0];
 186       u64 best_err = c1 - (u64)b1 * b1 / i + c2 - (u64)b2 * b2 / j;
 187       uns split_val = val[0];
 188       for (uns k = 1; k < cval - 1; k++)
 189         {
 190           uns b0 = val[k] * cnt[k];
 191           uns c0 = isqr(val[k]) * cnt[k];
 192           b1 += b0;
 193           b2 -= b0;
 194           c1 += c0;
 195           c2 -= c0;
 196           i += cnt[k];
 197           j -= cnt[k];
 198           u64 err = (u64)c1 - (u64)b1 * b1 / i + (u64)c2 - (u64)b2 * b2 / j;
 199           if (err < best_err)
 200             {
 201               best_err = err;
 202               split_val = val[k];
 203             }
 204         }
 205       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]);
 206
 207       /* Split region */
 208       block = region->blocks;
 209       region2 = regions + regions_count++;
 210       prequant_init_region(region);
 211       prequant_init_region(region2);
 212       while (block)
 213         {
 214           block2 = block->next;
 215           if (block->v[axis] <= split_val)
 216             prequant_add_block(region, block);
 217           else
 218             prequant_add_block(region2, block);
 219           block = block2;
 220         }
 221       prequant_finish_region(region);
 222       prequant_finish_region(region2);
 223       HEAP_INCREASE(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP, 1);
 224       heap[++heap_count] = region2;
 225       HEAP_INSERT(struct image_sig_region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP);
 226     }
 227
 228   DBG("Pre-quantized to %u regions", regions_count);
 229
 230   return regions_count;
 231 }
 232
 233 /* Post-quantization - run a few K-mean iterations to improve pre-quantized regions */
 234
 235 static uns
 236 postquant(struct image_sig_block *blocks, uns blocks_count, struct image_sig_region *regions, uns regions_count)
 237 {
 238   DBG("Starting post-quantization");
 239
 240   struct image_sig_block *blocks_end = blocks + blocks_count, *block;
 241   struct image_sig_region *regions_end = regions + regions_count, *region;
 242   uns error = 0, last_error;
 243
 244   /* Initialize regions and initial segmentation error */
 245   for (region = regions; region != regions_end; )
 246     {
 247       uns inv = 0xffffffffU / region->count;
 248       for (uns i = 0; i < IMAGE_VEC_F; i++)
 249         {
 250           region->a[i] = ((u64)region->b[i] * inv) >> 32;
 251           error += region->c[i] - region->a[i] * region->b[i];
 252         }
 253       region++;
 254     }
 255
 256   /* Convergation cycle */
 257   for (uns step = 0; step < image_sig_postquant_max_steps; step++)
 258     {
 259       DBG("Step...");
 260
 261       /* Clear regions */
 262       for (region = regions; region != regions_end; region++)
 263         {
 264           region->blocks = NULL;
 265           region->count = 0;
 266           bzero(region->b, sizeof(region->b));
 267           bzero(region->c, sizeof(region->c));
 268         }
 269
 270       /* Assign each block to its nearest pivot and accumulate region variables */
 271       for (block = blocks; block != blocks_end; block++)
 272         {
 273           struct image_sig_region *best_region = NULL;
 274           uns best_dist = ~0U;
 275           for (region = regions; region != regions_end; region++)
 276             {
 277               uns dist =
 278                 isqr(block->v[0] - region->a[0]) +
 279                 isqr(block->v[1] - region->a[1]) +
 280                 isqr(block->v[2] - region->a[2]) +
 281                 isqr(block->v[3] - region->a[3]) +
 282                 isqr(block->v[4] - region->a[4]) +
 283                 isqr(block->v[5] - region->a[5]);
 284               if (dist <= best_dist)
 285                 {
 286                   best_dist = dist;
 287                   best_region = region;
 288                 }
 289             }
 290           region = best_region;
 291           region->count++;
 292           block->next = region->blocks;
 293           region->blocks = block;
 294           for (uns i = 0; i < IMAGE_VEC_F; i++)
 295             {
 296               region->b[i] += block->v[i];
 297               region->c[i] += isqr(block->v[i]);
 298             }
 299         }
 300
 301       /* Finish regions, delete empty ones (should appear rarely), compute segmentation error */
 302       last_error = error;
 303       error = 0;
 304       for (region = regions; region != regions_end; )
 305         if (region->count)
 306           {
 307             uns inv = 0xffffffffU / region->count;
 308             for (uns i = 0; i < IMAGE_VEC_F; i++)
 309               {
 310                 region->a[i] = ((u64)region->b[i] * inv) >> 32;
 311                 error += region->c[i] - region->a[i] * region->b[i];
 312               }
 313             region++;
 314           }
 315         else
 316           {
 317             regions_end--;
 318             *region = *regions_end;
 319           }
 320
 321       DBG("last_error=%u error=%u", last_error, error);
 322
 323       /* Convergation criteria */
 324       if (step >= image_sig_postquant_min_steps)
 325         {
 326           if (error > last_error)
 327             break;
 328           u64 dif = last_error - error;
 329           if (dif * image_sig_postquant_threshold < (u64)last_error * 100)
 330             break;
 331         }
 332     }
 333
 334   DBG("Post-quantized to %u regions with average square error %u", regions_end - regions, error / blocks_count);
 335
 336   return regions_end - regions;
 337 }
 338
 339 void
 340 image_sig_segmentation(struct image_sig_data *data)
 341 {
 342   data->regions_count = prequant(data->blocks, data->blocks_count, data->regions);
 343 #ifdef LOCAL_DEBUG
 344   dump_segmentation(data->regions, data->regions_count);
 345 #endif
 346   data->regions_count = postquant(data->blocks, data->blocks_count, data->regions, data->regions_count);
 347 #ifdef LOCAL_DEBUG
 348   dump_segmentation(data->regions, data->regions_count);
 349 #endif
 350 }
 351