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