images/image-sig.c

   1 #undef LOCAL_DEBUG
   2
   3 #include "sherlock/sherlock.h"
   4 #include "lib/math.h"
   5 #include "lib/fastbuf.h"
   6 #include "images/images.h"
   7 #include "images/image-obj.h"
   8 #include "images/image-sig.h"
   9 #include "images/color.h"
  10
  11 #include <alloca.h>
  12
  13 struct block {
  14   uns l, u, v;          /* average Luv coefficients */
  15   uns lh, hl, hh;       /* energies in Daubechies wavelet bands */
  16 };
  17
  18 int
  19 compute_image_signature(struct image_data *image, struct image_signature *sig)
  20 {
  21   uns width = image->width;
  22   uns height = image->height;
  23
  24   if (width < 4 || height < 4)
  25     {
  26       DBG("Image too small... %dx%d", width, height);
  27       return 0;
  28     }
  29
  30   uns w = width >> 2;
  31   uns h = height >> 2;
  32   DBG("Computing signature for image %dx%d... %dx%d blocks", width, height, w, h);
  33   uns blocks_count = w * h;
  34   struct block *blocks = xmalloc(blocks_count * sizeof(struct block)), *block = blocks; /* FIXME: use mempool */
  35
  36   /* Every 4x4 block (FIXME: deal with smaller blocks near the edges) */
  37   byte *p = image->pixels;
  38   for (uns block_y = 0; block_y < h; block_y++, p += 3 * ((width & 3) + width * 3))
  39     for (uns block_x = 0; block_x < w; block_x++, p -= 3 * (4 * width - 4), block++)
  40       {
  41         int t[16], s[16], *tp = t;
  42
  43         /* Convert pixels to Luv color space and compute average coefficients
  44          * FIXME:
  45          * - could be MUCH faster with precomputed tables and integer arithmetic...
  46          *   I will propably use interpolation in 3-dim array */
  47         uns l_sum = 0;
  48         uns u_sum = 0;
  49         uns v_sum = 0;
  50         for (uns y = 0; y < 4; y++, p += 3 * (width - 4))
  51           for (uns x = 0; x < 4; x++, p += 3)
  52             {
  53               byte luv[3];
  54               srgb_to_luv_pixel(luv, p);
  55               l_sum += *tp++ = luv[0];
  56               u_sum += luv[1];
  57               v_sum += luv[2];
  58             }
  59
  60         block->l = (l_sum >> 4);
  61         block->u = (u_sum >> 4);
  62         block->v = (v_sum >> 4);
  63
  64         /* Apply Daubechies wavelet transformation
  65          * FIXME:
  66          * - MMX/SSE instructions or tables could be faster
  67          * - maybe it would be better to compute Luv and wavelet separately because of processor cache or MMX/SSE
  68          * - eliminate slow square roots
  69          * - what about Haar transformation? */
  70
  71 #define DAUB_0  31651 /* (1 + sqrt 3) / (4 * sqrt 2) */
  72 #define DAUB_1  54822 /* (3 + sqrt 3) / (4 * sqrt 2) */
  73 #define DAUB_2  14689 /* (3 - sqrt 3) / (4 * sqrt 2) */
  74 #define DAUB_3  -8481 /* (1 - sqrt 3) / (4 * sqrt 2) */
  75
  76         /* ... to the rows */
  77         uns i;
  78         for (i = 0; i < 16; i += 4)
  79           {
  80             s[i + 0] = (DAUB_0 * t[i + 2] + DAUB_1 * t[i + 3] + DAUB_2 * t[i + 0] + DAUB_3 * t[i + 1]) / 0x10000;
  81             s[i + 1] = (DAUB_0 * t[i + 0] + DAUB_1 * t[i + 1] + DAUB_2 * t[i + 2] + DAUB_3 * t[i + 3]) / 0x10000;
  82             s[i + 2] = (DAUB_3 * t[i + 2] - DAUB_2 * t[i + 3] + DAUB_1 * t[i + 0] - DAUB_0 * t[i + 1]) / 0x10000;
  83             s[i + 3] = (DAUB_3 * t[i + 0] - DAUB_2 * t[i + 1] + DAUB_1 * t[i + 2] - DAUB_0 * t[i + 3]) / 0x10000;
  84           }
  85
  86         /* ... and to the columns... skip LL band */
  87         for (i = 0; i < 2; i++)
  88           {
  89             t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
  90             t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
  91           }
  92         for (; i < 4; i++)
  93           {
  94             t[i + 0] = (DAUB_0 * s[i + 8] + DAUB_1 * s[i +12] + DAUB_2 * s[i + 0] + DAUB_3 * s[i + 4]) / 0x1000;
  95             t[i + 4] = (DAUB_0 * s[i + 0] + DAUB_1 * s[i + 4] + DAUB_2 * s[i + 8] + DAUB_3 * s[i +12]) / 0x1000;
  96             t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
  97             t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
  98           }
  99
 100         /* Extract energies in LH, HL and HH bands */
 101         block->lh = CLAMP((int)(sqrt(t[8] * t[8] + t[9] * t[9] + t[12] * t[12] + t[13] * t[13]) / 16), 0, 255);
 102         block->hl = CLAMP((int)(sqrt(t[2] * t[2] + t[3] * t[3] + t[6] * t[6] + t[7] * t[7]) / 16), 0, 255);
 103         block->hh = CLAMP((int)(sqrt(t[10] * t[10] + t[11] * t[11] + t[14] * t[14] + t[15] * t[15]) / 16), 0, 255);
 104       }
 105
 106   /* FIXME: simple average is for testing pusposes only */
 107   uns l_sum = 0;
 108   uns u_sum = 0;
 109   uns v_sum = 0;
 110   uns lh_sum = 0;
 111   uns hl_sum = 0;
 112   uns hh_sum = 0;
 113   for (uns i = 0; i < blocks_count; i++)
 114     {
 115       l_sum += blocks[i].l;
 116       u_sum += blocks[i].u;
 117       v_sum += blocks[i].v;
 118       lh_sum += blocks[i].lh;
 119       hl_sum += blocks[i].hl;
 120       hh_sum += blocks[i].hh;
 121     }
 122
 123   sig->vec.f[0] = l_sum / blocks_count;
 124   sig->vec.f[1] = u_sum / blocks_count;
 125   sig->vec.f[2] = v_sum / blocks_count;
 126   sig->vec.f[3] = lh_sum / blocks_count;
 127   sig->vec.f[4] = hl_sum / blocks_count;
 128   sig->vec.f[5] = hh_sum / blocks_count;
 129
 130   sig->len = 0;
 131
 132   xfree(blocks);
 133
 134   DBG("Resulting signature is (%s)", stk_print_image_vector(&sig->vec));
 135   return 1;
 136 }
 137