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 "lib/conf.h"
  16 #include "lib/heap.h"
  17 #include "images/math.h"
  18 #include "images/images.h"
  19 #include "images/color.h"
  20 #include "images/signature.h"
  21
  22 #include <alloca.h>
  23
  24 static double image_sig_inertia_scale[3] = { 3, 1, 0.3 };
  25
  26 struct block {
  27   u32 area;             /* block area in pixels (usually 16) */
  28   u32 v[IMAGE_VEC_F];
  29   u32 x, y;             /* block position */
  30   struct block *next;
  31 };
  32
  33 int
  34 compute_image_signature(struct image_thread *thread UNUSED, struct image_signature *sig, struct image *image)
  35 {
  36   bzero(sig, sizeof(*sig));
  37   ASSERT((image->flags & IMAGE_PIXEL_FORMAT) == COLOR_SPACE_RGB);
  38   uns cols = image->cols;
  39   uns rows = image->rows;
  40   uns row_size = image->row_size;
  41
  42   uns w = (cols + 3) >> 2;
  43   uns h = (rows + 3) >> 2;
  44
  45   DBG("Computing signature for image of %ux%u pixels (%ux%u blocks)", cols, rows, w, h);
  46
  47   uns blocks_count = w * h;
  48   struct image_sig_block *blocks = xmalloc(blocks_count * sizeof(struct image_sig_block)), *block = blocks;
  49
  50   /* Every block of 4x4 pixels */
  51   byte *row_start = image->pixels;
  52   for (uns block_y = 0; block_y < h; block_y++, row_start += row_size * 4)
  53     {
  54       byte *p = row_start;
  55       for (uns block_x = 0; block_x < w; block_x++, p += 12, block++)
  56         {
  57           int t[16], s[16], *tp = t;
  58           block->x = block_x;
  59           block->y = block_y;
  60
  61           /* Convert pixels to Luv color space and compute average coefficients */
  62           uns l_sum = 0;
  63           uns u_sum = 0;
  64           uns v_sum = 0;
  65           byte *p2 = p;
  66           if ((!(cols & 3) || block_x < w - 1) && (!(rows & 3) || block_y < h - 1))
  67             {
  68               for (uns y = 0; y < 4; y++, p2 += row_size - 12)
  69                 for (uns x = 0; x < 4; x++, p2 += 3)
  70                   {
  71                     byte luv[3];
  72                     srgb_to_luv_pixel(luv, p2);
  73                     l_sum += *tp++ = luv[0];
  74                     u_sum += luv[1];
  75                     v_sum += luv[2];
  76                   }
  77               block->area = 16;
  78               block->v[0] = (l_sum >> 4);
  79               block->v[1] = (u_sum >> 4);
  80               block->v[2] = (v_sum >> 4);
  81             }
  82           /* Incomplete square near the edge */
  83           else
  84             {
  85               uns x, y;
  86               uns square_cols = (block_x < w - 1 || !(cols & 3)) ? 4 : cols & 3;
  87               uns square_rows = (block_y < h - 1 || !(rows & 3)) ? 4 : rows & 3;
  88               for (y = 0; y < square_rows; y++, p2 += row_size)
  89                 {
  90                   byte *p3 = p2;
  91                   for (x = 0; x < square_cols; x++, p3 += 3)
  92                     {
  93                       byte luv[3];
  94                       srgb_to_luv_pixel(luv, p3);
  95                       l_sum += *tp++ = luv[0];
  96                       u_sum += luv[1];
  97                       v_sum += luv[2];
  98                     }
  99                   for (; x < 4; x++)
 100                     {
 101                       *tp = tp[-square_cols];
 102                       tp++;
 103                     }
 104                 }
 105               for (; y < 4; y++)
 106                 for (x = 0; x < 4; x++)
 107                   {
 108                     *tp = tp[-square_rows * 4];
 109                     tp++;
 110                   }
 111               block->area = square_cols * square_rows;
 112               uns div = 0x10000 / block->area;
 113               block->v[0] = (l_sum * div) >> 16;
 114               block->v[1] = (u_sum * div) >> 16;
 115               block->v[2] = (v_sum * div) >> 16;
 116             }
 117
 118           /* Apply Daubechies wavelet transformation */
 119
 120 #         define DAUB_0 31651   /* (1 + sqrt 3) / (4 * sqrt 2) * 0x10000 */
 121 #         define DAUB_1 54822   /* (3 + sqrt 3) / (4 * sqrt 2) * 0x10000 */
 122 #         define DAUB_2 14689   /* (3 - sqrt 3) / (4 * sqrt 2) * 0x10000 */
 123 #         define DAUB_3 -8481   /* (1 - sqrt 3) / (4 * sqrt 2) * 0x10000 */
 124
 125           /* ... to the rows */
 126           uns i;
 127           for (i = 0; i < 16; i += 4)
 128             {
 129               s[i + 0] = (DAUB_0 * t[i + 2] + DAUB_1 * t[i + 3] + DAUB_2 * t[i + 0] + DAUB_3 * t[i + 1]) / 0x10000;
 130               s[i + 1] = (DAUB_0 * t[i + 0] + DAUB_1 * t[i + 1] + DAUB_2 * t[i + 2] + DAUB_3 * t[i + 3]) / 0x10000;
 131               s[i + 2] = (DAUB_3 * t[i + 2] - DAUB_2 * t[i + 3] + DAUB_1 * t[i + 0] - DAUB_0 * t[i + 1]) / 0x10000;
 132               s[i + 3] = (DAUB_3 * t[i + 0] - DAUB_2 * t[i + 1] + DAUB_1 * t[i + 2] - DAUB_0 * t[i + 3]) / 0x10000;
 133             }
 134
 135           /* ... and to the columns... skip LL band */
 136           for (i = 0; i < 2; i++)
 137             {
 138               t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x2000;
 139               t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x2000;
 140             }
 141           for (; i < 4; i++)
 142             {
 143               t[i + 0] = (DAUB_0 * s[i + 8] + DAUB_1 * s[i +12] + DAUB_2 * s[i + 0] + DAUB_3 * s[i + 4]) / 0x2000;
 144               t[i + 4] = (DAUB_0 * s[i + 0] + DAUB_1 * s[i + 4] + DAUB_2 * s[i + 8] + DAUB_3 * s[i +12]) / 0x2000;
 145               t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x2000;
 146               t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x2000;
 147             }
 148
 149           /* Extract energies in LH, HL and HH bands */
 150           block->v[3] = fast_sqrt_u16(isqr(t[8]) + isqr(t[9]) + isqr(t[12]) + isqr(t[13]));
 151           block->v[4] = fast_sqrt_u16(isqr(t[2]) + isqr(t[3]) + isqr(t[6]) + isqr(t[7]));
 152           block->v[5] = fast_sqrt_u16(isqr(t[10]) + isqr(t[11]) + isqr(t[14]) + isqr(t[15]));
 153         }
 154     }
 155
 156   /* FIXME: simple average is for testing pusposes only */
 157   uns l_sum = 0;
 158   uns u_sum = 0;
 159   uns v_sum = 0;
 160   uns lh_sum = 0;
 161   uns hl_sum = 0;
 162   uns hh_sum = 0;
 163   for (uns i = 0; i < blocks_count; i++)
 164     {
 165       l_sum += blocks[i].v[0];
 166       u_sum += blocks[i].v[1];
 167       v_sum += blocks[i].v[2];
 168       lh_sum += blocks[i].v[3];
 169       hl_sum += blocks[i].v[4];
 170       hh_sum += blocks[i].v[5];
 171     }
 172
 173   sig->vec.f[0] = l_sum / blocks_count;
 174   sig->vec.f[1] = u_sum / blocks_count;
 175   sig->vec.f[2] = v_sum / blocks_count;
 176   sig->vec.f[3] = lh_sum / blocks_count;
 177   sig->vec.f[4] = hl_sum / blocks_count;
 178   sig->vec.f[5] = hh_sum / blocks_count;
 179
 180   if (cols < image_sig_min_width || rows < image_sig_min_height)
 181     {
 182       xfree(blocks);
 183       return 1;
 184     }
 185
 186   /* Quantize blocks to image regions */
 187   struct image_sig_region regions[IMAGE_REG_MAX];
 188   sig->len = image_sig_segmentation(blocks, blocks_count, regions);
 189
 190   /* For each region */
 191   u64 w_total = 0;
 192   uns w_border = (MIN(w, h) + 3) / 4;
 193   uns w_mul = 127 * 256 / w_border;
 194   for (uns i = 0; i < sig->len; i++)
 195     {
 196       struct image_sig_region *r = regions + i;
 197       DBG("Processing region %u: count=%u", i, r->count);
 198       ASSERT(r->count);
 199
 200       /* Copy texture properties */
 201       sig->reg[i].f[0] = r->a[0];
 202       sig->reg[i].f[1] = r->a[1];
 203       sig->reg[i].f[2] = r->a[2];
 204       sig->reg[i].f[3] = r->a[3];
 205       sig->reg[i].f[4] = r->a[4];
 206       sig->reg[i].f[5] = r->a[5];
 207
 208       /* Compute coordinates centroid and region weight */
 209       u64 x_avg = 0, y_avg = 0, w_sum = 0;
 210       for (struct image_sig_block *b = r->blocks; b; b = b->next)
 211         {
 212           x_avg += b->x;
 213           y_avg += b->y;
 214           uns d = b->x;
 215           d = MIN(d, b->y);
 216           d = MIN(d, w - b->x - 1);
 217           d = MIN(d, h - b->y - 1);
 218           if (d >= w_border)
 219             w_sum += 128;
 220           else
 221             w_sum += 128 + (d - w_border) * w_mul / 256;
 222         }
 223       w_total += w_sum;
 224       r->w_sum = w_sum;
 225       x_avg /= r->count;
 226       y_avg /= r->count;
 227       DBG("  centroid=(%u %u)", (uns)x_avg, (uns)y_avg);
 228
 229       /* Compute normalized inertia */
 230       u64 sum1 = 0, sum2 = 0, sum3 = 0;
 231       for (struct image_sig_block *b = r->blocks; b; b = b->next)
 232         {
 233           uns inc2 = isqr(x_avg - b->x) + isqr(y_avg - b->y);
 234           uns inc1 = sqrt(inc2);
 235           sum1 += inc1;
 236           sum2 += inc2;
 237           sum3 += inc1 * inc2;
 238         }
 239       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);
 240       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);
 241       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);
 242
 243     }
 244
 245   /* Compute average differences */
 246   u64 df = 0, dh = 0;
 247
 248   if (sig->len < 2)
 249     {
 250       sig->df = 1;
 251       sig->dh = 1;
 252     }
 253   else
 254     {
 255       uns cnt = 0;
 256       for (uns i = 0; i < sig->len; i++)
 257         for (uns j = i + 1; j < sig->len; j++)
 258           {
 259             uns d = 0;
 260             for (uns k = 0; k < IMAGE_REG_F; k++)
 261               d += isqr(sig->reg[i].f[k] - sig->reg[j].f[k]);
 262             df += sqrt(d);
 263             d = 0;
 264             for (uns k = 0; k < IMAGE_REG_H; k++)
 265               d += isqr(sig->reg[i].h[k] - sig->reg[j].h[k]);
 266             dh += sqrt(d);
 267             cnt++;
 268           }
 269       sig->df = CLAMP(df / cnt, 1, 255);
 270       sig->dh = CLAMP(dh / cnt, 1, 65535);
 271     }
 272   DBG("Average regions difs: df=%u dh=%u", sig->df, sig->dh);
 273
 274   /* Compute normalized weights */
 275   uns wa = 128, wb = 128;
 276   for (uns i = sig->len; --i > 0; )
 277     {
 278       struct image_sig_region *r = regions + i;
 279       wa -= sig->reg[i].wa = CLAMP(r->count * 128 / blocks_count, 1, (int)(wa - i));
 280       wb -= sig->reg[i].wb = CLAMP(r->w_sum * 128 / w_total, 1, (int)(wa - i));
 281     }
 282   sig->reg[0].wa = wa;
 283   sig->reg[0].wb = wb;
 284
 285   /* Dump regions features */
 286 #ifdef LOCAL_DEBUG
 287   for (uns i = 0; i < sig->len; i++)
 288     {
 289       byte buf[IMAGE_REGION_DUMP_MAX];
 290       image_region_dump(buf, sig->reg + i);
 291       DBG("region %u: features=%s", i, buf);
 292     }
 293 #endif
 294
 295   xfree(blocks);
 296
 297   return 1;
 298 }
 299