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 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   uns cols = image->cols;
 154   uns rows = image->rows;
 155
 156   /* FIXME: deal with smaller images */
 157   if (cols < 4 || cols < 4)
 158     {
 159       image_thread_err_format(thread, IMAGE_ERR_INVALID_DIMENSIONS, "Image too small... %ux%u", cols, rows);
 160       return 0;
 161     }
 162
 163   uns w = cols >> 2;
 164   uns h = rows >> 2;
 165   DBG("Computing signature for image %dx%d... %dx%d blocks", cols, rows, w, h);
 166   uns blocks_count = w * h;
 167   struct block *blocks = xmalloc(blocks_count * sizeof(struct block)), *block = blocks; /* FIXME: use mempool */
 168
 169   /* Every block of 4x4 pixels (FIXME: deal with smaller blocks near the edges) */
 170   byte *p = image->pixels;
 171   for (uns block_y = 0; block_y < h; block_y++, p += 3 * ((cols & 3) + cols * 3))
 172     for (uns block_x = 0; block_x < w; block_x++, p -= 3 * (4 * cols - 4), block++)
 173       {
 174         block->x = block_x;
 175         block->y = block_y;
 176
 177         int t[16], s[16], *tp = t;
 178
 179         /* Convert pixels to Luv color space and compute average coefficients
 180          * FIXME:
 181          * - could be MUCH faster with precomputed tables and integer arithmetic...
 182          *   I will propably use interpolation in 3-dim array */
 183         uns l_sum = 0;
 184         uns u_sum = 0;
 185         uns v_sum = 0;
 186         for (uns y = 0; y < 4; y++, p += 3 * (cols - 4))
 187           for (uns x = 0; x < 4; x++, p += 3)
 188             {
 189               byte luv[3];
 190               srgb_to_luv_pixel(luv, p);
 191               l_sum += *tp++ = luv[0];
 192               u_sum += luv[1];
 193               v_sum += luv[2];
 194             }
 195
 196         block->l = (l_sum >> 4);
 197         block->u = (u_sum >> 4);
 198         block->v = (v_sum >> 4);
 199
 200         /* Apply Daubechies wavelet transformation
 201          * FIXME:
 202          * - MMX/SSE instructions or tables could be faster
 203          * - maybe it would be better to compute Luv and wavelet separately because of processor cache or MMX/SSE
 204          * - eliminate slow square roots
 205          * - what about Haar transformation? */
 206
 207 #define DAUB_0  31651 /* (1 + sqrt 3) / (4 * sqrt 2) */
 208 #define DAUB_1  54822 /* (3 + sqrt 3) / (4 * sqrt 2) */
 209 #define DAUB_2  14689 /* (3 - sqrt 3) / (4 * sqrt 2) */
 210 #define DAUB_3  -8481 /* (1 - sqrt 3) / (4 * sqrt 2) */
 211
 212         /* ... to the rows */
 213         uns i;
 214         for (i = 0; i < 16; i += 4)
 215           {
 216             s[i + 0] = (DAUB_0 * t[i + 2] + DAUB_1 * t[i + 3] + DAUB_2 * t[i + 0] + DAUB_3 * t[i + 1]) / 0x10000;
 217             s[i + 1] = (DAUB_0 * t[i + 0] + DAUB_1 * t[i + 1] + DAUB_2 * t[i + 2] + DAUB_3 * t[i + 3]) / 0x10000;
 218             s[i + 2] = (DAUB_3 * t[i + 2] - DAUB_2 * t[i + 3] + DAUB_1 * t[i + 0] - DAUB_0 * t[i + 1]) / 0x10000;
 219             s[i + 3] = (DAUB_3 * t[i + 0] - DAUB_2 * t[i + 1] + DAUB_1 * t[i + 2] - DAUB_0 * t[i + 3]) / 0x10000;
 220           }
 221
 222         /* ... and to the columns... skip LL band */
 223         for (i = 0; i < 2; i++)
 224           {
 225             t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
 226             t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
 227           }
 228         for (; i < 4; i++)
 229           {
 230             t[i + 0] = (DAUB_0 * s[i + 8] + DAUB_1 * s[i +12] + DAUB_2 * s[i + 0] + DAUB_3 * s[i + 4]) / 0x1000;
 231             t[i + 4] = (DAUB_0 * s[i + 0] + DAUB_1 * s[i + 4] + DAUB_2 * s[i + 8] + DAUB_3 * s[i +12]) / 0x1000;
 232             t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
 233             t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
 234           }
 235
 236         /* Extract energies in LH, HL and HH bands */
 237         block->lh = CLAMP((int)(sqrt(t[8] * t[8] + t[9] * t[9] + t[12] * t[12] + t[13] * t[13]) / 16), 0, 255);
 238         block->hl = CLAMP((int)(sqrt(t[2] * t[2] + t[3] * t[3] + t[6] * t[6] + t[7] * t[7]) / 16), 0, 255);
 239         block->hh = CLAMP((int)(sqrt(t[10] * t[10] + t[11] * t[11] + t[14] * t[14] + t[15] * t[15]) / 16), 0, 255);
 240       }
 241
 242   /* FIXME: simple average is for testing pusposes only */
 243   uns l_sum = 0;
 244   uns u_sum = 0;
 245   uns v_sum = 0;
 246   uns lh_sum = 0;
 247   uns hl_sum = 0;
 248   uns hh_sum = 0;
 249   for (uns i = 0; i < blocks_count; i++)
 250     {
 251       l_sum += blocks[i].l;
 252       u_sum += blocks[i].u;
 253       v_sum += blocks[i].v;
 254       lh_sum += blocks[i].lh;
 255       hl_sum += blocks[i].hl;
 256       hh_sum += blocks[i].hh;
 257     }
 258
 259   sig->vec.f[0] = l_sum / blocks_count;
 260   sig->vec.f[1] = u_sum / blocks_count;
 261   sig->vec.f[2] = v_sum / blocks_count;
 262   sig->vec.f[3] = lh_sum / blocks_count;
 263   sig->vec.f[4] = hl_sum / blocks_count;
 264   sig->vec.f[5] = hh_sum / blocks_count;
 265
 266   /* Quantize blocks to image regions */
 267   struct region regions[IMAGE_REG_MAX];
 268   sig->len = compute_k_means(blocks, blocks_count, regions, MIN(blocks_count, IMAGE_REG_MAX));
 269
 270   /* For each region */
 271   u64 w_total = 0;
 272   uns w_border = (MIN(w, h) + 3) / 4;
 273   uns w_mul = 127 * 256 / w_border;
 274   for (uns i = 0; i < sig->len; i++)
 275     {
 276       struct region *r = regions + i;
 277       DBG("Processing region %u: count=%u", i, r->count);
 278       ASSERT(r->count);
 279
 280       /* Copy texture properties */
 281       sig->reg[i].f[0] = r->l;
 282       sig->reg[i].f[1] = r->u;
 283       sig->reg[i].f[2] = r->v;
 284       sig->reg[i].f[3] = r->lh;
 285       sig->reg[i].f[4] = r->hl;
 286       sig->reg[i].f[5] = r->hh;
 287
 288       /* Compute coordinates centroid and region weight */
 289       u64 x_avg = 0, y_avg = 0, w_sum = 0;
 290       for (struct block *b = r->blocks; b; b = b->next)
 291         {
 292           x_avg += b->x;
 293           y_avg += b->y;
 294           uns d = b->x;
 295           d = MIN(d, b->y);
 296           d = MIN(d, w - b->x - 1);
 297           d = MIN(d, h - b->y - 1);
 298           if (d >= w_border)
 299             w_sum += 128;
 300           else
 301             w_sum += 128 + (d - w_border) * w_mul / 256;
 302         }
 303       w_total += w_sum;
 304       r->w_sum = w_sum;
 305       x_avg /= r->count;
 306       y_avg /= r->count;
 307       DBG("  centroid=(%u %u)", (uns)x_avg, (uns)y_avg);
 308
 309       /* Compute normalized inertia */
 310       u64 sum1 = 0, sum2 = 0, sum3 = 0;
 311       for (struct block *b = r->blocks; b; b = b->next)
 312         {
 313           uns inc2 = dist(x_avg, b->x) + dist(y_avg, b->y);
 314           uns inc1 = sqrt(inc2);
 315           sum1 += inc1;
 316           sum2 += inc2;
 317           sum3 += inc1 * inc2;
 318         }
 319       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);
 320       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);
 321       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);
 322
 323     }
 324
 325   /* Compute average differences */
 326   u64 df = 0, dh = 0;
 327   uns cnt = 0;
 328   for (uns i = 0; i < sig->len; i++)
 329     for (uns j = i + 1; j < sig->len; j++)
 330       {
 331         uns d = 0;
 332         for (uns k = 0; k < IMAGE_REG_F; k++)
 333           d += dist(sig->reg[i].f[k], sig->reg[j].f[k]);
 334         df += sqrt(d);
 335         d = 0;
 336         for (uns k = 0; k < IMAGE_REG_H; k++)
 337           d += dist(sig->reg[i].h[k], sig->reg[j].h[k]);
 338         dh += sqrt(d);
 339         cnt++;
 340       }
 341   sig->df = df / cnt;
 342   sig->dh = dh / cnt;
 343   DBG("Average regions difs: df=%u dh=%u", sig->df, sig->dh);
 344
 345   /* Compute normalized weights */
 346   uns wa = 128, wb = 128;
 347   for (uns i = sig->len; --i > 0; )
 348     {
 349       struct region *r = regions + i;
 350       wa -= sig->reg[i].wa = CLAMP(r->count * 128 / blocks_count, 1, (int)(wa - i));
 351       wb -= sig->reg[i].wb = CLAMP(r->w_sum * 128 / w_total, 1, (int)(wa - i));
 352     }
 353   sig->reg[0].wa = wa;
 354   sig->reg[0].wb = wb;
 355
 356   /* Dump regions features */
 357 #ifdef LOCAL_DEBUG
 358   for (uns i = 0; i < sig->len; i++)
 359     {
 360       byte buf[IMAGE_REGION_DUMP_MAX];
 361       image_region_dump(buf, sig->reg + i);
 362       DBG("region %u: features=%s", i, buf);
 363     }
 364 #endif
 365
 366   xfree(blocks);
 367
 368   return 1;
 369 }
 370