images/sig-init.c

   1 /*
   2  *      Image Library -- Computation of image signatures
   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 static double image_sig_inertia_scale[3] = { 3, 1, 0.3 };
  25
  26 struct block {
  27   u32 area;             /* block area in pixels (usually 16) */
  28   u32 v[IMAGE_VEC_F];
  29   u32 x, y;             /* block position */
  30   struct block *next;
  31 };
  32
  33 struct region {
  34   struct block *blocks;
  35   u32 count;
  36   u32 a[IMAGE_VEC_F];
  37   u32 b[IMAGE_VEC_F];
  38   u32 c[IMAGE_VEC_F];
  39   u64 e;
  40   u64 w_sum;
  41 };
  42
  43 static inline uns
  44 dist(uns a, uns b)
  45 {
  46   int d = a - b;
  47   return d * d;
  48 }
  49
  50 #ifdef LOCAL_DEBUG
  51 static void
  52 dump_segmentation(struct region *regions, uns regions_count, uns cols, uns rows)
  53 {
  54   uns size = (cols + 1) * rows;
  55   byte buf[size];
  56   bzero(buf, size);
  57   for (uns i = 0; i < regions_count; i++)
  58     {
  59       byte c = (i < 10) ? '0' + i : 'A' - 10 + i;
  60       for (struct block *b = regions[i].blocks; b; b = b->next)
  61         buf[b->x + b->y * (cols + 1)] = c;
  62     }
  63   for (uns i = 0; i < rows; i++)
  64     log(L_DEBUG, "%s", &buf[i * (cols + 1)]);
  65 }
  66 #endif
  67
  68 /* Pre-quantization - recursively split groups of blocks with large error */
  69
  70 static inline void
  71 prequant_init_region(struct region *region)
  72 {
  73   bzero(region, sizeof(*region));
  74 }
  75
  76 static inline void
  77 prequant_add_block(struct region *region, struct block *block)
  78 {
  79   block->next = region->blocks;
  80   region->blocks = block;
  81   region->count++;
  82   for (uns i = 0; i < IMAGE_VEC_F; i++)
  83     {
  84       region->b[i] += block->v[i];
  85       region->c[i] += isqr(block->v[i]);
  86     }
  87 }
  88
  89 static void
  90 prequant_finish_region(struct region *region)
  91 {
  92   if (region->count < 2)
  93     memcpy(region->c, region->a, sizeof(region->c));
  94   else
  95     {
  96       u64 a = 0;
  97       region->e = 0;
  98       for (uns i = 0; i < IMAGE_VEC_F; i++)
  99         {
 100           region->e += region->c[i];
 101           a += (u64)region->b[i] * region->b[i];
 102         }
 103       region->e -= a / region->count;
 104     }
 105 }
 106
 107 static inline uns
 108 prequant_heap_cmp(struct region *a, struct region *b)
 109 {
 110   return a->e > b->e;
 111 }
 112
 113 #define ASORT_PREFIX(x) prequant_##x
 114 #define ASORT_KEY_TYPE u32
 115 #define ASORT_ELT(i) val[i]
 116 #define ASORT_EXTRA_ARGS , u32 *val
 117 #include "lib/arraysort.h"
 118
 119 static uns
 120 prequant(struct block *blocks, uns blocks_count, struct region *regions)
 121 {
 122   DBG("Starting pre-quantization");
 123
 124   uns regions_count, heap_count, axis, cov;
 125   struct block *blocks_end = blocks + blocks_count, *block, *block2;
 126   struct region *heap[IMAGE_REG_MAX + 1], *region, *region2;
 127
 128   /* Initialize single region with all blocks */
 129   regions_count = heap_count = 1;
 130   heap[1] = regions;
 131   prequant_init_region(regions);
 132   for (block = blocks; block != blocks_end; block++)
 133     prequant_add_block(regions, block);
 134   prequant_finish_region(regions);
 135
 136   /* Main cycle */
 137   while (regions_count < IMAGE_REG_MAX &&
 138       regions_count <= DARY_LEN(image_sig_prequant_thresholds) && heap_count)
 139     {
 140       region = heap[1];
 141       DBG("Step... regions_count=%u heap_count=%u region->count=%u, region->e=%u",
 142           regions_count, heap_count, region->count, (uns)region->e);
 143       if (region->count < 2 ||
 144           region->e < image_sig_prequant_thresholds[regions_count - 1] * region->count)
 145         {
 146           HEAP_DELMIN(struct region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP);
 147           continue;
 148         }
 149
 150       /* Select axis to split - the one with maximum covariance */
 151       axis = 0;
 152       cov = region->count * region->c[0] - region->b[0];
 153       for (uns i = 1; i < 6; i++)
 154         {
 155           uns j = region->count * region->c[i] - region->b[i];
 156           if (j > cov)
 157             {
 158               axis = i;
 159               cov = j;
 160             }
 161         }
 162       DBG("Splitting axis %u with covariance %u", axis, cov / region->count);
 163
 164       /* Sort values on the split axis */
 165       u32 val[region->count];
 166       block = region->blocks;
 167       for (uns i = 0; i < region->count; i++, block = block->next)
 168         val[i] = block->v[axis];
 169       prequant_sort(region->count, val);
 170
 171       /* Select split value - to minimize error */
 172       uns b1 = 0, c1 = 0, b2 = region->b[axis], c2 = region->c[axis];
 173       uns i = 0, j = region->count, best_v = 0;
 174       u64 best_err = 0xffffffffffffffff;
 175       while (i < region->count)
 176         {
 177           u64 err = (u64)i * c1 - (u64)b1 * b1 + (u64)j * c2 - (u64)b2 * b2;
 178           if (err < best_err)
 179             {
 180               best_err = err;
 181               best_v = val[i];
 182             }
 183           uns sqr = isqr(val[i]);
 184           b1 += val[i];
 185           b2 -= val[i];
 186           c1 += sqr;
 187           c2 -= sqr;
 188           i++;
 189           j--;
 190         }
 191       uns split_val = best_v;
 192       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]);
 193
 194       /* Split region */
 195       block = region->blocks;
 196       region2 = regions + regions_count++;
 197       prequant_init_region(region);
 198       prequant_init_region(region2);
 199       while (block)
 200         {
 201           block2 = block->next;
 202           if (block->v[axis] < split_val)
 203             prequant_add_block(region, block);
 204           else
 205             prequant_add_block(region2, block);
 206           block = block2;
 207         }
 208       prequant_finish_region(region);
 209       prequant_finish_region(region2);
 210       HEAP_INCREASE(struct region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP, 1);
 211       heap[++heap_count] = region2;
 212       HEAP_INSERT(struct region *, heap, heap_count, prequant_heap_cmp, HEAP_SWAP);
 213     }
 214
 215   DBG("Pre-quantized to %u regions", regions_count);
 216
 217   return regions_count;
 218 }
 219
 220 /* Post-quantization - run a few K-mean iterations to improve pre-quantized regions */
 221
 222 static uns
 223 postquant(struct block *blocks, uns blocks_count, struct region *regions, uns regions_count)
 224 {
 225   DBG("Starting post-quantization");
 226
 227   struct block *blocks_end = blocks + blocks_count, *block;
 228   struct region *regions_end = regions + regions_count, *region;
 229   uns error = 0, last_error;
 230
 231   /* Initialize regions and initial segmentation error */
 232   for (region = regions; region != regions_end; )
 233     {
 234       uns inv = 0xffffffffU / region->count;
 235       for (uns i = 0; i < IMAGE_VEC_F; i++)
 236         {
 237           region->a[i] = ((u64)region->b[i] * inv) >> 32;
 238           error += region->c[i] - region->a[i] * region->b[i];
 239         }
 240       region++;
 241     }
 242
 243   /* Convergation cycle */
 244   for (uns step = 0; step < image_sig_postquant_max_steps; step++)
 245     {
 246       DBG("Step...");
 247
 248       /* Clear regions */
 249       for (region = regions; region != regions_end; region++)
 250         {
 251           region->blocks = NULL;
 252           region->count = 0;
 253           bzero(region->b, sizeof(region->b));
 254           bzero(region->c, sizeof(region->c));
 255         }
 256
 257       /* Assign each block to its nearest pivot and accumulate region variables */
 258       for (block = blocks; block != blocks_end; block++)
 259         {
 260           struct region *best_region = NULL;
 261           uns best_dist = ~0U;
 262           for (region = regions; region != regions_end; region++)
 263             {
 264               uns dist =
 265                 isqr(block->v[0] - region->a[0]) +
 266                 isqr(block->v[1] - region->a[1]) +
 267                 isqr(block->v[2] - region->a[2]) +
 268                 isqr(block->v[3] - region->a[3]) +
 269                 isqr(block->v[4] - region->a[4]) +
 270                 isqr(block->v[5] - region->a[5]);
 271               if (dist <= best_dist)
 272                 {
 273                   best_dist = dist;
 274                   best_region = region;
 275                 }
 276             }
 277           region = best_region;
 278           region->count++;
 279           block->next = region->blocks;
 280           region->blocks = block;
 281           for (uns i = 0; i < IMAGE_VEC_F; i++)
 282             {
 283               region->b[i] += block->v[i];
 284               region->c[i] += isqr(block->v[i]);
 285             }
 286         }
 287
 288       /* Finish regions, delete empty ones (should appear rarely), compute segmentation error */
 289       last_error = error;
 290       error = 0;
 291       for (region = regions; region != regions_end; )
 292         if (region->count)
 293           {
 294             uns inv = 0xffffffffU / region->count;
 295             for (uns i = 0; i < IMAGE_VEC_F; i++)
 296               {
 297                 region->a[i] = ((u64)region->b[i] * inv) >> 32;
 298                 error += region->c[i] - region->a[i] * region->b[i];
 299               }
 300             region++;
 301           }
 302         else
 303           {
 304             regions_end--;
 305             *region = *regions_end;
 306           }
 307
 308       DBG("last_error=%u error=%u", last_error, error);
 309
 310       /* Convergation criteria */
 311       if (step >= image_sig_postquant_min_steps)
 312         {
 313           if (error > last_error)
 314             break;
 315           u64 dif = last_error - error;
 316           if (dif * image_sig_postquant_threshold < last_error * 100)
 317             break;
 318         }
 319     }
 320
 321   DBG("Post-quantized to %u regions with average square error %u", regions_end - regions, error / blocks_count);
 322
 323   return regions_end - regions;
 324 }
 325
 326 static inline uns
 327 segmentation(struct block *blocks, uns blocks_count, struct region *regions, uns width UNUSED, uns height UNUSED)
 328 {
 329   uns regions_count;
 330   regions_count = prequant(blocks, blocks_count, regions);
 331 #ifdef LOCAL_DEBUG
 332   dump_segmentation(regions, regions_count, width, height);
 333 #endif
 334   regions_count = postquant(blocks, blocks_count, regions, regions_count);
 335 #ifdef LOCAL_DEBUG
 336   dump_segmentation(regions, regions_count, width, height);
 337 #endif
 338   return regions_count;
 339 }
 340
 341 int
 342 compute_image_signature(struct image_thread *thread UNUSED, struct image_signature *sig, struct image *image)
 343 {
 344   bzero(sig, sizeof(*sig));
 345   ASSERT((image->flags & IMAGE_PIXEL_FORMAT) == COLOR_SPACE_RGB);
 346   uns cols = image->cols;
 347   uns rows = image->rows;
 348   uns row_size = image->row_size;
 349
 350   uns w = (cols + 3) >> 2;
 351   uns h = (rows + 3) >> 2;
 352
 353   DBG("Computing signature for image of %ux%u pixels (%ux%u blocks)", cols, rows, w, h);
 354
 355   uns blocks_count = w * h;
 356   struct block *blocks = xmalloc(blocks_count * sizeof(struct block)), *block = blocks;
 357
 358   /* Every block of 4x4 pixels */
 359   byte *row_start = image->pixels;
 360   for (uns block_y = 0; block_y < h; block_y++, row_start += row_size * 4)
 361     {
 362       byte *p = row_start;
 363       for (uns block_x = 0; block_x < w; block_x++, p += 12, block++)
 364         {
 365           int t[16], s[16], *tp = t;
 366           block->x = block_x;
 367           block->y = block_y;
 368
 369           /* Convert pixels to Luv color space and compute average coefficients */
 370           uns l_sum = 0;
 371           uns u_sum = 0;
 372           uns v_sum = 0;
 373           byte *p2 = p;
 374           if ((!(cols & 3) || block_x < w - 1) && (!(rows & 3) || block_y < h - 1))
 375             {
 376               for (uns y = 0; y < 4; y++, p2 += row_size - 12)
 377                 for (uns x = 0; x < 4; x++, p2 += 3)
 378                   {
 379                     byte luv[3];
 380                     srgb_to_luv_pixel(luv, p2);
 381                     l_sum += *tp++ = luv[0];
 382                     u_sum += luv[1];
 383                     v_sum += luv[2];
 384                   }
 385               block->area = 16;
 386               block->v[0] = (l_sum >> 4);
 387               block->v[1] = (u_sum >> 4);
 388               block->v[2] = (v_sum >> 4);
 389             }
 390           /* Incomplete square near the edge */
 391           else
 392             {
 393               uns x, y;
 394               uns square_cols = (block_x < w - 1 || !(cols & 3)) ? 4 : cols & 3;
 395               uns square_rows = (block_y < h - 1 || !(rows & 3)) ? 4 : rows & 3;
 396               for (y = 0; y < square_rows; y++, p2 += row_size)
 397                 {
 398                   byte *p3 = p2;
 399                   for (x = 0; x < square_cols; x++, p3 += 3)
 400                     {
 401                       byte luv[3];
 402                       srgb_to_luv_pixel(luv, p3);
 403                       l_sum += *tp++ = luv[0];
 404                       u_sum += luv[1];
 405                       v_sum += luv[2];
 406                     }
 407                   for (; x < 4; x++)
 408                     {
 409                       *tp = tp[-square_cols];
 410                       tp++;
 411                     }
 412                 }
 413               for (; y < 4; y++)
 414                 for (x = 0; x < 4; x++)
 415                   {
 416                     *tp = tp[-square_rows * 4];
 417                     tp++;
 418                   }
 419               block->area = square_cols * square_rows;
 420               uns div = 0x10000 / block->area;
 421               block->v[0] = (l_sum * div) >> 16;
 422               block->v[1] = (u_sum * div) >> 16;
 423               block->v[2] = (v_sum * div) >> 16;
 424             }
 425
 426           /* Apply Daubechies wavelet transformation */
 427
 428 #         define DAUB_0 31651   /* (1 + sqrt 3) / (4 * sqrt 2) * 0x10000 */
 429 #         define DAUB_1 54822   /* (3 + sqrt 3) / (4 * sqrt 2) * 0x10000 */
 430 #         define DAUB_2 14689   /* (3 - sqrt 3) / (4 * sqrt 2) * 0x10000 */
 431 #         define DAUB_3 -8481   /* (1 - sqrt 3) / (4 * sqrt 2) * 0x10000 */
 432
 433           /* ... to the rows */
 434           uns i;
 435           for (i = 0; i < 16; i += 4)
 436             {
 437               s[i + 0] = (DAUB_0 * t[i + 2] + DAUB_1 * t[i + 3] + DAUB_2 * t[i + 0] + DAUB_3 * t[i + 1]) / 0x10000;
 438               s[i + 1] = (DAUB_0 * t[i + 0] + DAUB_1 * t[i + 1] + DAUB_2 * t[i + 2] + DAUB_3 * t[i + 3]) / 0x10000;
 439               s[i + 2] = (DAUB_3 * t[i + 2] - DAUB_2 * t[i + 3] + DAUB_1 * t[i + 0] - DAUB_0 * t[i + 1]) / 0x10000;
 440               s[i + 3] = (DAUB_3 * t[i + 0] - DAUB_2 * t[i + 1] + DAUB_1 * t[i + 2] - DAUB_0 * t[i + 3]) / 0x10000;
 441             }
 442
 443           /* ... and to the columns... skip LL band */
 444           for (i = 0; i < 2; i++)
 445             {
 446               t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x2000;
 447               t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x2000;
 448             }
 449           for (; i < 4; i++)
 450             {
 451               t[i + 0] = (DAUB_0 * s[i + 8] + DAUB_1 * s[i +12] + DAUB_2 * s[i + 0] + DAUB_3 * s[i + 4]) / 0x2000;
 452               t[i + 4] = (DAUB_0 * s[i + 0] + DAUB_1 * s[i + 4] + DAUB_2 * s[i + 8] + DAUB_3 * s[i +12]) / 0x2000;
 453               t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x2000;
 454               t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x2000;
 455             }
 456
 457           /* Extract energies in LH, HL and HH bands */
 458           block->v[3] = fast_sqrt_u16(isqr(t[8]) + isqr(t[9]) + isqr(t[12]) + isqr(t[13]));
 459           block->v[4] = fast_sqrt_u16(isqr(t[2]) + isqr(t[3]) + isqr(t[6]) + isqr(t[7]));
 460           block->v[5] = fast_sqrt_u16(isqr(t[10]) + isqr(t[11]) + isqr(t[14]) + isqr(t[15]));
 461         }
 462     }
 463
 464   /* FIXME: simple average is for testing pusposes only */
 465   uns l_sum = 0;
 466   uns u_sum = 0;
 467   uns v_sum = 0;
 468   uns lh_sum = 0;
 469   uns hl_sum = 0;
 470   uns hh_sum = 0;
 471   for (uns i = 0; i < blocks_count; i++)
 472     {
 473       l_sum += blocks[i].v[0];
 474       u_sum += blocks[i].v[1];
 475       v_sum += blocks[i].v[2];
 476       lh_sum += blocks[i].v[3];
 477       hl_sum += blocks[i].v[4];
 478       hh_sum += blocks[i].v[5];
 479     }
 480
 481   sig->vec.f[0] = l_sum / blocks_count;
 482   sig->vec.f[1] = u_sum / blocks_count;
 483   sig->vec.f[2] = v_sum / blocks_count;
 484   sig->vec.f[3] = lh_sum / blocks_count;
 485   sig->vec.f[4] = hl_sum / blocks_count;
 486   sig->vec.f[5] = hh_sum / blocks_count;
 487
 488   if (cols < image_sig_min_width || rows < image_sig_min_height)
 489     {
 490       xfree(blocks);
 491       return 1;
 492     }
 493
 494   /* Quantize blocks to image regions */
 495   struct region regions[IMAGE_REG_MAX];
 496   sig->len = segmentation(blocks, blocks_count, regions, w, h);
 497
 498   /* For each region */
 499   u64 w_total = 0;
 500   uns w_border = (MIN(w, h) + 3) / 4;
 501   uns w_mul = 127 * 256 / w_border;
 502   for (uns i = 0; i < sig->len; i++)
 503     {
 504       struct region *r = regions + i;
 505       DBG("Processing region %u: count=%u", i, r->count);
 506       ASSERT(r->count);
 507
 508       /* Copy texture properties */
 509       sig->reg[i].f[0] = r->a[0];
 510       sig->reg[i].f[1] = r->a[1];
 511       sig->reg[i].f[2] = r->a[2];
 512       sig->reg[i].f[3] = r->a[3];
 513       sig->reg[i].f[4] = r->a[4];
 514       sig->reg[i].f[5] = r->a[5];
 515
 516       /* Compute coordinates centroid and region weight */
 517       u64 x_avg = 0, y_avg = 0, w_sum = 0;
 518       for (struct block *b = r->blocks; b; b = b->next)
 519         {
 520           x_avg += b->x;
 521           y_avg += b->y;
 522           uns d = b->x;
 523           d = MIN(d, b->y);
 524           d = MIN(d, w - b->x - 1);
 525           d = MIN(d, h - b->y - 1);
 526           if (d >= w_border)
 527             w_sum += 128;
 528           else
 529             w_sum += 128 + (d - w_border) * w_mul / 256;
 530         }
 531       w_total += w_sum;
 532       r->w_sum = w_sum;
 533       x_avg /= r->count;
 534       y_avg /= r->count;
 535       DBG("  centroid=(%u %u)", (uns)x_avg, (uns)y_avg);
 536
 537       /* Compute normalized inertia */
 538       u64 sum1 = 0, sum2 = 0, sum3 = 0;
 539       for (struct block *b = r->blocks; b; b = b->next)
 540         {
 541           uns inc2 = dist(x_avg, b->x) + dist(y_avg, b->y);
 542           uns inc1 = sqrt(inc2);
 543           sum1 += inc1;
 544           sum2 += inc2;
 545           sum3 += inc1 * inc2;
 546         }
 547       sig->reg[i].h[0] = CLAMP(image_sig_inertia_scale[0] * sum1 * ((3 * M_PI * M_PI) / 2) * pow(r->count, -1.5), 0, 65535);
 548       sig->reg[i].h[1] = CLAMP(image_sig_inertia_scale[1] * sum2 * ((4 * M_PI * M_PI * M_PI) / 2) / ((u64)r->count * r->count), 0, 65535);
 549       sig->reg[i].h[2] = CLAMP(image_sig_inertia_scale[2] * sum3 * ((5 * M_PI * M_PI * M_PI * M_PI) / 2) * pow(r->count, -2.5), 0, 65535);
 550
 551     }
 552
 553   /* Compute average differences */
 554   u64 df = 0, dh = 0;
 555
 556   if (sig->len < 2)
 557     {
 558       sig->df = 1;
 559       sig->dh = 1;
 560     }
 561   else
 562     {
 563       uns cnt = 0;
 564       for (uns i = 0; i < sig->len; i++)
 565         for (uns j = i + 1; j < sig->len; j++)
 566           {
 567             uns d = 0;
 568             for (uns k = 0; k < IMAGE_REG_F; k++)
 569               d += dist(sig->reg[i].f[k], sig->reg[j].f[k]);
 570             df += sqrt(d);
 571             d = 0;
 572             for (uns k = 0; k < IMAGE_REG_H; k++)
 573               d += dist(sig->reg[i].h[k], sig->reg[j].h[k]);
 574             dh += sqrt(d);
 575             cnt++;
 576           }
 577       sig->df = CLAMP(df / cnt, 1, 255);
 578       sig->dh = CLAMP(dh / cnt, 1, 65535);
 579     }
 580   DBG("Average regions difs: df=%u dh=%u", sig->df, sig->dh);
 581
 582   /* Compute normalized weights */
 583   uns wa = 128, wb = 128;
 584   for (uns i = sig->len; --i > 0; )
 585     {
 586       struct region *r = regions + i;
 587       wa -= sig->reg[i].wa = CLAMP(r->count * 128 / blocks_count, 1, (int)(wa - i));
 588       wb -= sig->reg[i].wb = CLAMP(r->w_sum * 128 / w_total, 1, (int)(wa - i));
 589     }
 590   sig->reg[0].wa = wa;
 591   sig->reg[0].wb = wb;
 592
 593   /* Dump regions features */
 594 #ifdef LOCAL_DEBUG
 595   for (uns i = 0; i < sig->len; i++)
 596     {
 597       byte buf[IMAGE_REGION_DUMP_MAX];
 598       image_region_dump(buf, sig->reg + i);
 599       DBG("region %u: features=%s", i, buf);
 600     }
 601 #endif
 602
 603   xfree(blocks);
 604
 605   return 1;
 606 }
 607