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