images/image-sig.c

   1 #define LOCAL_DEBUG
   2
   3 #include "sherlock/sherlock.h"
   4 #include "lib/math.h"
   5 #include "lib/fastbuf.h"
   6 #include "images/images.h"
   7
   8 /*
   9  * Color spaces
  10  *
  11  * http://www.tecgraf.puc-rio.br/~mgattass/color/ColorIndex.html
  12  *
  13  */
  14
  15 #define REF_WHITE_X 0.96422
  16 #define REF_WHITE_Y 1.
  17 #define REF_WHITE_Z 0.82521
  18
  19 /* sRGB to XYZ */
  20 static void
  21 srgb_to_xyz_slow(double srgb[3], double xyz[3])
  22 {
  23   double a[3];
  24   for (uns i = 0; i < 3; i++)
  25     if (srgb[i] > 0.04045)
  26       a[i] = pow((srgb[i] + 0.055) * (1 / 1.055), 2.4);
  27     else
  28       a[i] = srgb[i] * (1 / 12.92);
  29   xyz[0] = 0.412424 * a[0] + 0.357579 * a[1] + 0.180464 * a[2];
  30   xyz[1] = 0.212656 * a[0] + 0.715158 * a[1] + 0.072186 * a[2];
  31   xyz[2] = 0.019332 * a[0] + 0.119193 * a[1] + 0.950444 * a[2];
  32 }
  33
  34 /* XYZ to CIE-Luv */
  35 static void
  36 xyz_to_luv_slow(double xyz[3], double luv[3])
  37 {
  38   double sum = xyz[0] + 15 * xyz[1] + 3 * xyz[2];
  39   if (sum < 0.000001)
  40     luv[0] = luv[1] = luv[2] = 0;
  41   else
  42     {
  43       double var_u = 4 * xyz[0] / sum;
  44       double var_v = 9 * xyz[1] / sum;
  45       if (xyz[1] > 0.008856)
  46         luv[0] = 116 * pow(xyz[1], 1 / 3.) - 16;
  47       else
  48         luv[0] = (116 * 7.787) * xyz[1];
  49       luv[1] = luv[0] * (13 * (var_u - 4 * REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
  50       luv[2] = luv[0] * (13 * (var_v - 9 * REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
  51       /* intervals [0..100], [-134..220], [-140..122] */
  52     }
  53 }
  54
  55 struct block {
  56   uns l, u, v;          /* average Luv coefficients */
  57   uns lh, hl, hh;       /* energies in Daubechies wavelet bands */
  58 };
  59
  60 static void
  61 compute_image_area_signature(PixelPacket *pixels, uns width, uns height, struct image_signature *sig)
  62 {
  63   ASSERT(width >= 4 && height >= 4);
  64
  65   uns w = width >> 2;
  66   uns h = height >> 2;
  67   DBG("Computing signature for image %dx%d... %dx%d blocks", width, height, w, h);
  68   uns blocks_count = w * h;
  69   struct block *blocks = xmalloc(blocks_count * sizeof(struct block)), *block = blocks; /* FIXME: use mempool */
  70
  71   /* Every 4x4 block (FIXME: deal with smaller blocks near the edges) */
  72   PixelPacket *p = pixels;
  73   for (uns block_y = 0; block_y < h; block_y++, p += width & 3 + width * 3)
  74     for (uns block_x = 0; block_x < w; block_x++, p -= 4 * width - 4, block++)
  75       {
  76         int t[16], s[16], *tp = t;
  77
  78         /* Convert pixels to Luv color space and compute average coefficients
  79          * FIXME:
  80          * - could be MUCH faster with precomputed tables and integer arithmetic...
  81          *   I will propably use interpolation in 3-dim array */
  82         uns l_sum = 0;
  83         uns u_sum = 0;
  84         uns v_sum = 0;
  85         for (uns y = 0; y < 4; y++, p += width - 4)
  86           for (uns x = 0; x < 4; x++, p++)
  87             {
  88               double rgb[3], luv[3], xyz[3];
  89               rgb[0] = (p->red >> (QuantumDepth - 8)) / 255.;
  90               rgb[1] = (p->green >> (QuantumDepth - 8)) / 255.;
  91               rgb[2] = (p->blue >> (QuantumDepth - 8)) / 255.;
  92               srgb_to_xyz_slow(rgb, xyz);
  93               xyz_to_luv_slow(xyz, luv);
  94               l_sum += *tp++ = luv[0];
  95               u_sum += luv[1] + 150;
  96               v_sum += luv[2] + 150;
  97             }
  98
  99         block->l = l_sum;
 100         block->u = u_sum;
 101         block->v = v_sum;
 102
 103         /* Apply Daubechies wavelet transformation
 104          * FIXME:
 105          * - MMX/SSE instructions or tables could be faster
 106          * - maybe it would be better to compute Luv and wavelet separately because of processor cache or MMX/SSE
 107          * - eliminate slow square roots
 108          * - what about Haar transformation? */
 109
 110 #define DAUB_0  31651 /* (1 + sqrt 3) / (4 * sqrt 2) */
 111 #define DAUB_1  54822 /* (3 + sqrt 3) / (4 * sqrt 2) */
 112 #define DAUB_2  14689 /* (3 - sqrt 3) / (4 * sqrt 2) */
 113 #define DAUB_3  -8481 /* (1 - sqrt 3) / (4 * sqrt 2) */
 114
 115         /* ... to the rows */
 116         uns i;
 117         for (i = 0; i < 16; i += 4)
 118           {
 119             s[i + 0] = (DAUB_0 * t[i + 2] + DAUB_1 * t[i + 3] + DAUB_2 * t[i + 0] + DAUB_3 * t[i + 1]) / 0x10000;
 120             s[i + 1] = (DAUB_0 * t[i + 0] + DAUB_1 * t[i + 1] + DAUB_2 * t[i + 2] + DAUB_3 * t[i + 3]) / 0x10000;
 121             s[i + 2] = (DAUB_3 * t[i + 2] - DAUB_2 * t[i + 3] + DAUB_1 * t[i + 0] - DAUB_0 * t[i + 1]) / 0x10000;
 122             s[i + 3] = (DAUB_3 * t[i + 0] - DAUB_2 * t[i + 1] + DAUB_1 * t[i + 2] - DAUB_0 * t[i + 3]) / 0x10000;
 123           }
 124
 125         /* ... and to the columns... skip LL band */
 126         for (i = 0; i < 2; i++)
 127           {
 128             t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
 129             t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
 130           }
 131         for (; i < 4; i++)
 132           {
 133             t[i + 0] = (DAUB_0 * s[i + 8] + DAUB_1 * s[i +12] + DAUB_2 * s[i + 0] + DAUB_3 * s[i + 4]) / 0x1000;
 134             t[i + 4] = (DAUB_0 * s[i + 0] + DAUB_1 * s[i + 4] + DAUB_2 * s[i + 8] + DAUB_3 * s[i +12]) / 0x1000;
 135             t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
 136             t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
 137           }
 138
 139         /* Extract energies in LH, HL and HH bands */
 140         block->lh = sqrt(t[8] * t[8] + t[9] * t[9] + t[12] * t[12] + t[13] * t[13]);
 141         block->hl = sqrt(t[2] * t[2] + t[3] * t[3] + t[6] * t[6] + t[7] * t[7]);
 142         block->hh = sqrt(t[10] * t[10] + t[11] * t[11] + t[14] * t[14] + t[15] * t[15]);
 143       }
 144
 145   /* FIXME: simple average is for testing pusposes only */
 146   uns l_sum = 0;
 147   uns u_sum = 0;
 148   uns v_sum = 0;
 149   uns lh_sum = 0;
 150   uns hl_sum = 0;
 151   uns hh_sum = 0;
 152   for (uns i = 0; i < blocks_count; i++)
 153     {
 154       l_sum += blocks[i].l;
 155       u_sum += blocks[i].u;
 156       v_sum += blocks[i].v;
 157       lh_sum += blocks[i].lh;
 158       hl_sum += blocks[i].hl;
 159       hh_sum += blocks[i].hh;
 160     }
 161
 162   sig->vec.f[0] = l_sum / blocks_count;
 163   sig->vec.f[1] = u_sum / blocks_count;
 164   sig->vec.f[2] = v_sum / blocks_count;
 165   sig->vec.f[3] = lh_sum / blocks_count;
 166   sig->vec.f[4] = hl_sum / blocks_count;
 167   sig->vec.f[5] = hh_sum / blocks_count;
 168
 169   sig->len = 0;
 170
 171   xfree(blocks);
 172
 173   DBG("Resulting signature is (%s)", stk_print_image_vector(&sig->vec));
 174 }
 175
 176 static ExceptionInfo exception;
 177 static QuantizeInfo quantize_info;
 178 static ImageInfo *image_info;
 179
 180 void
 181 compute_image_signature_prepare(void)
 182 {
 183   InitializeMagick(NULL);
 184   GetExceptionInfo(&exception);
 185   image_info = CloneImageInfo(NULL);
 186   image_info->subrange = 1;
 187   GetQuantizeInfo(&quantize_info);
 188   quantize_info.colorspace = RGBColorspace;
 189 }
 190
 191 void
 192 compute_image_signature_finish(void)
 193 {
 194   DestroyImageInfo(image_info);
 195   DestroyExceptionInfo(&exception);
 196   DestroyMagick();
 197 }
 198
 199 int
 200 compute_image_signature(void *data, uns len, struct image_signature *sig)
 201 {
 202   Image *image = BlobToImage(image_info, data, len, &exception); /* Damn slow... takes most of CPU time :-/ */
 203   if (!image)
 204     die("Invalid image format");
 205   if (image->columns < 4 || image->rows < 4)
 206     {
 207       DBG("Image too small (%dx%d)", (int)image->columns, (int)image->rows);
 208       DestroyImage(image);
 209       return -1;
 210     }
 211   QuantizeImage(&quantize_info, image); /* Also slow... and propably not necessary... */
 212   PixelPacket *pixels = (PixelPacket *) AcquireImagePixels(image, 0, 0, image->columns, image->rows, &exception);
 213   compute_image_area_signature(pixels, image->columns, image->rows, sig);
 214   DestroyImage(image);
 215   return 0;
 216 }
 217