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