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 #define LOCAL_DEBUG
  11
  12 #include "sherlock/sherlock.h"
  13 #include "lib/math.h"
  14 #include "lib/fastbuf.h"
  15 #include "images/math.h"
  16 #include "images/images.h"
  17 #include "images/color.h"
  18 #include "images/signature.h"
  19 #include <alloca.h>
  20
  21 static double image_sig_inertia_scale[3] = { 3, 1, 0.3 };
  22
  23 struct block {
  24   u32 area;             /* block area in pixels (usually 16) */
  25   u32 l, u, v;          /* average Luv coefficients */
  26   u32 lh, hl, hh;       /* energies in Daubechies wavelet bands */
  27   u32 x, y;             /* block position */
  28   struct block *next;
  29 };
  30
  31 struct region {
  32   u32 l, u, v;
  33   u32 lh, hl, hh;
  34   u32 sum_l, sum_u, sum_v;
  35   u32 sum_lh, sum_hl, sum_hh;
  36   u32 count;
  37   u64 w_sum;
  38   struct block *blocks;
  39 };
  40
  41 static inline uns
  42 dist(uns a, uns b)
  43 {
  44   int d = a - b;
  45   return d * d;
  46 }
  47
  48 #ifdef LOCAL_DEBUG
  49 static void
  50 dump_segmentation(struct region *regions, uns regions_count, uns cols, uns rows)
  51 {
  52   uns size = (cols + 1) * rows;
  53   byte buf[size];
  54   bzero(buf, size);
  55   for (uns i = 0; i < regions_count; i++)
  56     {
  57       byte c = (i < 10) ? '0' + i : 'A' - 10 + i;
  58       for (struct block *b = regions[i].blocks; b; b = b->next)
  59         buf[b->x + b->y * (cols + 1)] = c;
  60     }
  61   for (uns i = 0; i < rows; i++)
  62     log(L_DEBUG, "%s", &buf[i * (cols + 1)]);
  63 }
  64 #endif
  65
  66 /* FIXME: SLOW! */
  67 static uns
  68 compute_k_means(struct block *blocks, uns blocks_count, struct region *regions, uns regions_count)
  69 {
  70   ASSERT(regions_count <= blocks_count);
  71   struct block *mean[IMAGE_REG_MAX], *b, *blocks_end = blocks + blocks_count;
  72   struct region *r, *regions_end = regions + regions_count;
  73
  74   /* Select means_count random blocks as initial regions pivots */
  75   if (regions_count <= blocks_count - regions_count)
  76     {
  77       for (b = blocks; b != blocks_end; b++)
  78         b->next = NULL;
  79       for (uns i = 0; i < regions_count; )
  80         {
  81           uns j = random_max(blocks_count);
  82           b = blocks + j;
  83           if (!b->next)
  84             b->next = mean[i++] = b;
  85         }
  86     }
  87   else
  88     {
  89       uns j = blocks_count;
  90       for (uns i = regions_count; i; j--)
  91         if (random_max(j) <= i)
  92           mean[--i] = blocks + j - 1;
  93     }
  94   r = regions;
  95   for (uns i = 0; i < regions_count; i++, r++)
  96     {
  97       b = mean[i];
  98       r->l = b->l;
  99       r->u = b->u;
 100       r->v = b->v;
 101       r->lh = b->lh;
 102       r->hl = b->hl;
 103       r->hh = b->hh;
 104     }
 105
 106   /* Convergation cycle */
 107   for (uns conv_i = 8; ; conv_i--)
 108     {
 109       for (r = regions; r != regions_end; r++)
 110         {
 111           r->sum_l = r->sum_u = r->sum_v = r->sum_lh = r->sum_hl = r->sum_hh = r->count = 0;
 112           r->blocks = NULL;
 113         }
 114
 115       /* Find nearest regions and accumulate averages */
 116       for (b = blocks; b != blocks_end; b++)
 117         {
 118           uns best_d = ~0U;
 119           struct region *best_r = NULL;
 120           for (r = regions; r != regions_end; r++)
 121             {
 122               uns d =
 123                 dist(r->l, b->l) +
 124                 dist(r->u, b->u) +
 125                 dist(r->v, b->v) +
 126                 dist(r->lh, b->lh) +
 127                 dist(r->hl, b->hl) +
 128                 dist(r->hh, b->hh);
 129               if (d < best_d)
 130                 {
 131                   best_d = d;
 132                   best_r = r;
 133                 }
 134             }
 135           best_r->sum_l += b->l;
 136           best_r->sum_u += b->u;
 137           best_r->sum_v += b->v;
 138           best_r->sum_lh += b->lh;
 139           best_r->sum_hl += b->hl;
 140           best_r->sum_hh += b->hh;
 141           best_r->count++;
 142           b->next = best_r->blocks;
 143           best_r->blocks = b;
 144         }
 145
 146       /* Compute new averages */
 147       for (r = regions; r != regions_end; r++)
 148         if (r->count)
 149           {
 150             r->l = r->sum_l / r->count;
 151             r->u = r->sum_u / r->count;
 152             r->v = r->sum_v / r->count;
 153             r->lh = r->sum_lh / r->count;
 154             r->hl = r->sum_hl / r->count;
 155             r->hh = r->sum_hh / r->count;
 156           }
 157
 158       if (!conv_i)
 159         break; // FIXME: convergation criteria
 160     }
 161
 162   /* Remove empty regions */
 163   struct region *r2 = regions;
 164   for (r = regions; r != regions_end; r++)
 165     if (r->count)
 166       *r2++ = *r;
 167   return r2 - regions;
 168 }
 169
 170 int
 171 compute_image_signature(struct image_thread *thread UNUSED, struct image_signature *sig, struct image *image)
 172 {
 173   bzero(sig, sizeof(*sig));
 174   ASSERT((image->flags & IMAGE_PIXEL_FORMAT) == COLOR_SPACE_RGB);
 175   uns cols = image->cols;
 176   uns rows = image->rows;
 177   uns row_size = image->row_size;
 178
 179   uns w = (cols + 3) >> 2;
 180   uns h = (rows + 3) >> 2;
 181
 182   DBG("Computing signature for image of %ux%u pixels (%ux%u blocks)", cols, rows, w, h);
 183
 184   uns blocks_count = w * h;
 185   struct block *blocks = xmalloc(blocks_count * sizeof(struct block)), *block = blocks;
 186
 187   /* Every block of 4x4 pixels */
 188   byte *row_start = image->pixels;
 189   for (uns block_y = 0; block_y < h; block_y++, row_start += row_size * 4)
 190     {
 191       byte *p = row_start;
 192       for (uns block_x = 0; block_x < w; block_x++, p += 12, block++)
 193         {
 194           int t[16], s[16], *tp = t;
 195           block->x = block_x;
 196           block->y = block_y;
 197
 198           /* Convert pixels to Luv color space and compute average coefficients */
 199           uns l_sum = 0;
 200           uns u_sum = 0;
 201           uns v_sum = 0;
 202           byte *p2 = p;
 203           if ((!(cols & 3) || block_x < w - 1) && (!(rows & 3) || block_y < h - 1))
 204             {
 205               for (uns y = 0; y < 4; y++, p2 += row_size - 12)
 206                 for (uns x = 0; x < 4; x++, p2 += 3)
 207                   {
 208                     byte luv[3];
 209                     srgb_to_luv_pixel(luv, p2);
 210                     l_sum += *tp++ = luv[0];
 211                     u_sum += luv[1];
 212                     v_sum += luv[2];
 213                   }
 214               block->area = 16;
 215               block->l = (l_sum >> 4);
 216               block->u = (u_sum >> 4);
 217               block->v = (v_sum >> 4);
 218             }
 219           /* Incomplete square near the edge */
 220           else
 221             {
 222               uns x, y;
 223               uns square_cols = (block_x < w - 1) ? 4 : cols & 3;
 224               uns square_rows = (block_y < h - 1) ? 4 : rows & 3;
 225               for (y = 0; y < square_rows; y++, p2 += row_size)
 226                 {
 227                   byte *p3 = p2;
 228                   for (x = 0; x < square_cols; x++, p3 += 3)
 229                     {
 230                       byte luv[3];
 231                       srgb_to_luv_pixel(luv, p3);
 232                       l_sum += *tp++ = luv[0];
 233                       u_sum += luv[1];
 234                       v_sum += luv[2];
 235                     }
 236                   for (; x < 4; x++)
 237                     {
 238                       *tp = tp[-square_cols];
 239                       tp++;
 240                     }
 241                 }
 242               for (; y < 4; y++)
 243                 for (x = 0; x < 4; x++)
 244                   {
 245                     *tp = tp[-square_rows * 4];
 246                     tp++;
 247                   }
 248               block->area = square_cols * square_rows;
 249               uns div = 0x10000 / block->area;
 250               block->l = (l_sum * div) >> 16;
 251               block->u = (u_sum * div) >> 16;
 252               block->v = (v_sum * div) >> 16;
 253             }
 254
 255           /* Apply Daubechies wavelet transformation */
 256
 257 #         define DAUB_0 31651   /* (1 + sqrt 3) / (4 * sqrt 2) * 0x10000 */
 258 #         define DAUB_1 54822   /* (3 + sqrt 3) / (4 * sqrt 2) * 0x10000 */
 259 #         define DAUB_2 14689   /* (3 - sqrt 3) / (4 * sqrt 2) * 0x10000 */
 260 #         define DAUB_3 -8481   /* (1 - sqrt 3) / (4 * sqrt 2) * 0x10000 */
 261
 262           /* ... to the rows */
 263           uns i;
 264           for (i = 0; i < 16; i += 4)
 265             {
 266               s[i + 0] = (DAUB_0 * t[i + 2] + DAUB_1 * t[i + 3] + DAUB_2 * t[i + 0] + DAUB_3 * t[i + 1]) / 0x10000;
 267               s[i + 1] = (DAUB_0 * t[i + 0] + DAUB_1 * t[i + 1] + DAUB_2 * t[i + 2] + DAUB_3 * t[i + 3]) / 0x10000;
 268               s[i + 2] = (DAUB_3 * t[i + 2] - DAUB_2 * t[i + 3] + DAUB_1 * t[i + 0] - DAUB_0 * t[i + 1]) / 0x10000;
 269               s[i + 3] = (DAUB_3 * t[i + 0] - DAUB_2 * t[i + 1] + DAUB_1 * t[i + 2] - DAUB_0 * t[i + 3]) / 0x10000;
 270             }
 271
 272           /* ... and to the columns... skip LL band */
 273           for (i = 0; i < 2; i++)
 274             {
 275               t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
 276               t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
 277             }
 278           for (; i < 4; i++)
 279             {
 280               t[i + 0] = (DAUB_0 * s[i + 8] + DAUB_1 * s[i +12] + DAUB_2 * s[i + 0] + DAUB_3 * s[i + 4]) / 0x1000;
 281               t[i + 4] = (DAUB_0 * s[i + 0] + DAUB_1 * s[i + 4] + DAUB_2 * s[i + 8] + DAUB_3 * s[i +12]) / 0x1000;
 282               t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
 283               t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
 284             }
 285
 286           /* Extract energies in LH, HL and HH bands */
 287           block->lh = fast_sqrt_u16(isqr(t[8]) + isqr(t[9]) + isqr(t[12]) + isqr(t[13]));
 288           block->hl = fast_sqrt_u16(isqr(t[2]) + isqr(t[3]) + isqr(t[6]) + isqr(t[7]));
 289           block->hh = fast_sqrt_u16(isqr(t[10]) + isqr(t[11]) + isqr(t[14]) + isqr(t[15]));
 290         }
 291     }
 292
 293   /* FIXME: simple average is for testing pusposes only */
 294   uns l_sum = 0;
 295   uns u_sum = 0;
 296   uns v_sum = 0;
 297   uns lh_sum = 0;
 298   uns hl_sum = 0;
 299   uns hh_sum = 0;
 300   for (uns i = 0; i < blocks_count; i++)
 301     {
 302       l_sum += blocks[i].l;
 303       u_sum += blocks[i].u;
 304       v_sum += blocks[i].v;
 305       lh_sum += blocks[i].lh;
 306       hl_sum += blocks[i].hl;
 307       hh_sum += blocks[i].hh;
 308     }
 309
 310   sig->vec.f[0] = l_sum / blocks_count;
 311   sig->vec.f[1] = u_sum / blocks_count;
 312   sig->vec.f[2] = v_sum / blocks_count;
 313   sig->vec.f[3] = lh_sum / blocks_count;
 314   sig->vec.f[4] = hl_sum / blocks_count;
 315   sig->vec.f[5] = hh_sum / blocks_count;
 316
 317   if (cols < image_sig_min_width || rows < image_sig_min_height)
 318     return 1;
 319
 320   /* Quantize blocks to image regions */
 321   struct region regions[IMAGE_REG_MAX];
 322   sig->len = compute_k_means(blocks, blocks_count, regions, MIN(blocks_count, IMAGE_REG_MAX));
 323
 324   /* For each region */
 325   u64 w_total = 0;
 326   uns w_border = (MIN(w, h) + 3) / 4;
 327   uns w_mul = 127 * 256 / w_border;
 328   for (uns i = 0; i < sig->len; i++)
 329     {
 330       struct region *r = regions + i;
 331       DBG("Processing region %u: count=%u", i, r->count);
 332       ASSERT(r->count);
 333
 334       /* Copy texture properties */
 335       sig->reg[i].f[0] = r->l;
 336       sig->reg[i].f[1] = r->u;
 337       sig->reg[i].f[2] = r->v;
 338       sig->reg[i].f[3] = r->lh;
 339       sig->reg[i].f[4] = r->hl;
 340       sig->reg[i].f[5] = r->hh;
 341
 342       /* Compute coordinates centroid and region weight */
 343       u64 x_avg = 0, y_avg = 0, w_sum = 0;
 344       for (struct block *b = r->blocks; b; b = b->next)
 345         {
 346           x_avg += b->x;
 347           y_avg += b->y;
 348           uns d = b->x;
 349           d = MIN(d, b->y);
 350           d = MIN(d, w - b->x - 1);
 351           d = MIN(d, h - b->y - 1);
 352           if (d >= w_border)
 353             w_sum += 128;
 354           else
 355             w_sum += 128 + (d - w_border) * w_mul / 256;
 356         }
 357       w_total += w_sum;
 358       r->w_sum = w_sum;
 359       x_avg /= r->count;
 360       y_avg /= r->count;
 361       DBG("  centroid=(%u %u)", (uns)x_avg, (uns)y_avg);
 362
 363       /* Compute normalized inertia */
 364       u64 sum1 = 0, sum2 = 0, sum3 = 0;
 365       for (struct block *b = r->blocks; b; b = b->next)
 366         {
 367           uns inc2 = dist(x_avg, b->x) + dist(y_avg, b->y);
 368           uns inc1 = sqrt(inc2);
 369           sum1 += inc1;
 370           sum2 += inc2;
 371           sum3 += inc1 * inc2;
 372         }
 373       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);
 374       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);
 375       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);
 376
 377     }
 378
 379   /* Compute average differences */
 380   u64 df = 0, dh = 0;
 381
 382   if (sig->len < 2)
 383     {
 384       sig->df = 1;
 385       sig->dh = 1;
 386     }
 387   else
 388     {
 389       uns cnt = 0;
 390       for (uns i = 0; i < sig->len; i++)
 391         for (uns j = i + 1; j < sig->len; j++)
 392           {
 393             uns d = 0;
 394             for (uns k = 0; k < IMAGE_REG_F; k++)
 395               d += dist(sig->reg[i].f[k], sig->reg[j].f[k]);
 396             df += sqrt(d);
 397             d = 0;
 398             for (uns k = 0; k < IMAGE_REG_H; k++)
 399               d += dist(sig->reg[i].h[k], sig->reg[j].h[k]);
 400             dh += sqrt(d);
 401             cnt++;
 402           }
 403       sig->df = CLAMP(df / cnt, 1, 255);
 404       sig->dh = CLAMP(dh / cnt, 1, 65535);
 405     }
 406   DBG("Average regions difs: df=%u dh=%u", sig->df, sig->dh);
 407
 408   /* Compute normalized weights */
 409   uns wa = 128, wb = 128;
 410   for (uns i = sig->len; --i > 0; )
 411     {
 412       struct region *r = regions + i;
 413       wa -= sig->reg[i].wa = CLAMP(r->count * 128 / blocks_count, 1, (int)(wa - i));
 414       wb -= sig->reg[i].wb = CLAMP(r->w_sum * 128 / w_total, 1, (int)(wa - i));
 415     }
 416   sig->reg[0].wa = wa;
 417   sig->reg[0].wb = wb;
 418
 419   /* Dump regions features */
 420 #ifdef LOCAL_DEBUG
 421   for (uns i = 0; i < sig->len; i++)
 422     {
 423       byte buf[IMAGE_REGION_DUMP_MAX];
 424       image_region_dump(buf, sig->reg + i);
 425       DBG("region %u: features=%s", i, buf);
 426     }
 427   dump_segmentation(regions, sig->len, w, h);
 428 #endif
 429
 430   xfree(blocks);
 431
 432   return 1;
 433 }
 434