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