images/block_info.c

   1 /*#include <stdio.h>
   2 #include <magick/api.h>*/
   3
   4 #include <assert.h>
   5 #include <stdlib.h>
   6 #include <string.h>
   7 #include <math.h>
   8 #include "img.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   uns i;
  27   for (i = 0; i < 3; i++)
  28     if (srgb[i] > 0.04045)
  29       a[i] = pow((srgb[i] + 0.055) * (1 / 1.055), 2.4);
  30     else
  31       a[i] = srgb[i] * (1 / 12.92);
  32   xyz[0] = 0.412424 * a[0] + 0.357579 * a[1] + 0.180464 * a[2];
  33   xyz[1] = 0.212656 * a[0] + 0.715158 * a[1] + 0.072186 * a[2];
  34   xyz[2] = 0.019332 * a[0] + 0.119193 * a[1] + 0.950444 * a[2];
  35 }
  36
  37 /* XYZ to CIE-Luv */
  38 static void
  39 xyz_to_luv_slow(double xyz[3], double luv[3])
  40 {
  41   double sum = xyz[0] + 15 * xyz[1] + 3 * xyz[2];
  42   if (sum < 0.000001)
  43     luv[0] = luv[1] = luv[2] = 0;
  44   else
  45     {
  46       double var_u = 4 * xyz[0] / sum;
  47       double var_v = 9 * xyz[1] / sum;
  48       if (xyz[1] > 0.008856)
  49         luv[0] = 116 * pow(xyz[1], 1 / 3.) - 16;
  50       else
  51         luv[0] = (116 * 7.787) * xyz[1];
  52       luv[1] = luv[0] * (13 * (var_u - 4 * REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
  53       luv[2] = luv[0] * (13 * (var_v - 9 * REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
  54       /* intervals [0..100], [-134..220], [-140..122] */
  55     }
  56 }
  57
  58 struct BlockInfo*
  59 computeBlockInfo(PixelPacket *pixels, uns width, uns height, uns *count)
  60 {
  61   assert(width >= 4 && height >= 4);
  62
  63   uns w = width >> 2;
  64   uns h = height >> 2;
  65   fprintf(stderr, "Computing signature for image %dx%d... %dx%d blocks", width, height, w, h);
  66   uns blocks_count = w * h;
  67   struct BlockInfo *blocks = malloc(blocks_count * sizeof(struct BlockInfo)), *block = blocks; /* FIXME: use mempool */
  68
  69   /* Every 4x4 block (FIXME: deal with smaller blocks near the edges) */
  70   PixelPacket *p = pixels;
  71   for (uns block_y = 0; block_y < h; block_y++, p += (width & 3) + 3*width){
  72     for (uns block_x = 0; block_x < w; block_x++, p += 4 - 4*width, block++){
  73         int t[16], s[16], *tp = t;
  74
  75         /* Convert pixels to Luv color space and compute average coefficients
  76          * FIXME:
  77          * - could be MUCH faster with precomputed tables and integer arithmetic...
  78          *   I will propably use interpolation in 3-dim array */
  79         uns l_sum = 0;
  80         uns u_sum = 0;
  81         uns v_sum = 0;
  82         for (uns y = 0; y < 4; y++, p += width - 4){
  83           for (uns x = 0; x < 4; x++, p++)
  84             {
  85               double rgb[3], luv[3], xyz[3];
  86               rgb[0] = (p->red >> (QuantumDepth - 8)) / 255.;
  87               rgb[1] = (p->green >> (QuantumDepth - 8)) / 255.;
  88               rgb[2] = (p->blue >> (QuantumDepth - 8)) / 255.;
  89               srgb_to_xyz_slow(rgb, xyz);
  90               xyz_to_luv_slow(xyz, luv);
  91               l_sum += *tp++ = luv[0];
  92               u_sum += luv[1] + 150;
  93               v_sum += luv[2] + 150;
  94               /*fprintf(stderr, "'%u, %u'; ", (p - pixels)%width ,  (p - pixels)/width);*/
  95             }
  96             /*fprintf(stderr, "\n---\n");*/
  97         }
  98         block->l = l_sum;
  99         block->u = u_sum;
 100         block->v = v_sum;
 101
 102         /* Apply Daubechies wavelet transformation
 103          * FIXME:
 104          * - MMX/SSE instructions or tables could be faster
 105          * - maybe it would be better to compute Luv and wavelet separately because of processor cache or MMX/SSE
 106          * - eliminate slow square roots
 107          * - what about Haar transformation? */
 108
 109 #define DAUB_0  31651 /* (1 + sqrt 3) / (4 * sqrt 2) */
 110 #define DAUB_1  54822 /* (3 + sqrt 3) / (4 * sqrt 2) */
 111 #define DAUB_2  14689 /* (3 - sqrt 3) / (4 * sqrt 2) */
 112 #define DAUB_3  -8481 /* (1 - sqrt 3) / (4 * sqrt 2) */
 113
 114         /* ... to the rows */
 115         uns i;
 116         for (i = 0; i < 16; i += 4)
 117           {
 118             s[i + 0] = (DAUB_0 * t[i + 2] + DAUB_1 * t[i + 3] + DAUB_2 * t[i + 0] + DAUB_3 * t[i + 1]) / 0x10000;
 119             s[i + 1] = (DAUB_0 * t[i + 0] + DAUB_1 * t[i + 1] + DAUB_2 * t[i + 2] + DAUB_3 * t[i + 3]) / 0x10000;
 120             s[i + 2] = (DAUB_3 * t[i + 2] - DAUB_2 * t[i + 3] + DAUB_1 * t[i + 0] - DAUB_0 * t[i + 1]) / 0x10000;
 121             s[i + 3] = (DAUB_3 * t[i + 0] - DAUB_2 * t[i + 1] + DAUB_1 * t[i + 2] - DAUB_0 * t[i + 3]) / 0x10000;
 122           }
 123
 124         /* ... and to the columns... skip LL band */
 125         for (i = 0; i < 2; i++)
 126           {
 127             t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
 128             t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
 129           }
 130         for (; i < 4; i++)
 131           {
 132             t[i + 0] = (DAUB_0 * s[i + 8] + DAUB_1 * s[i +12] + DAUB_2 * s[i + 0] + DAUB_3 * s[i + 4]) / 0x1000;
 133             t[i + 4] = (DAUB_0 * s[i + 0] + DAUB_1 * s[i + 4] + DAUB_2 * s[i + 8] + DAUB_3 * s[i +12]) / 0x1000;
 134             t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
 135             t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
 136           }
 137
 138         /* Extract energies in LH, HL and HH bands */
 139         block->lh = sqrt(t[8] * t[8] + t[9] * t[9] + t[12] * t[12] + t[13] * t[13]);
 140         block->hl = sqrt(t[2] * t[2] + t[3] * t[3] + t[6] * t[6] + t[7] * t[7]);
 141         block->hh = sqrt(t[10] * t[10] + t[11] * t[11] + t[14] * t[14] + t[15] * t[15]);
 142       }
 143   }
 144   *count=blocks_count;
 145   return blocks;
 146 }