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