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