]> mj.ucw.cz Git - libucw.git/commitdiff
Merge with git+ssh://cvs.ucw.cz/projects/sherlock/GIT/sherlock.git
authorPavel Charvat <pavel.charvat@netcentrum.cz>
Tue, 25 Jul 2006 07:39:51 +0000 (09:39 +0200)
committerPavel Charvat <pavel.charvat@netcentrum.cz>
Tue, 25 Jul 2006 07:39:51 +0000 (09:39 +0200)
32 files changed:
images/Makefile [new file with mode: 0644]
images/block_info.c [new file with mode: 0644]
images/color.c [new file with mode: 0644]
images/color.h [new file with mode: 0644]
images/color.t [new file with mode: 0644]
images/decomp.c [new file with mode: 0644]
images/dup-cmp.c [new file with mode: 0644]
images/dup-cmp.h [new file with mode: 0644]
images/hilbert-test.c [new file with mode: 0644]
images/hilbert-test.t [new file with mode: 0644]
images/hilbert.h [new file with mode: 0644]
images/image-idx.c [new file with mode: 0644]
images/image-obj.c [new file with mode: 0644]
images/image-obj.h [new file with mode: 0644]
images/image-sig.c [new file with mode: 0644]
images/image-sig.h [new file with mode: 0644]
images/image-test.c [new file with mode: 0644]
images/image-tool.c [new file with mode: 0644]
images/image-walk.h [new file with mode: 0644]
images/image.c [new file with mode: 0644]
images/images.h [new file with mode: 0644]
images/img.h [new file with mode: 0644]
images/io-libjpeg.c [new file with mode: 0644]
images/io-libmagick.c [new file with mode: 0644]
images/io-libpng.c [new file with mode: 0644]
images/io-libungif.c [new file with mode: 0644]
images/io-main.c [new file with mode: 0644]
images/k_means.c [new file with mode: 0644]
images/kd-tree.c [new file with mode: 0644]
images/kd-tree.h [new file with mode: 0644]
images/scale-gen.h [new file with mode: 0644]
images/scale.c [new file with mode: 0644]

diff --git a/images/Makefile b/images/Makefile
new file mode 100644 (file)
index 0000000..412d959
--- /dev/null
@@ -0,0 +1,64 @@
+# Testing dir... code will be moved somewhere else... maybe to trash :-)
+
+DIRS+=images
+
+PROGS+=$(addprefix $(o)/images/,image-tool)
+
+LIBIMAGES_MODS=image scale io-main
+LIBIMAGES=$(o)/images/libimages.$(LS)
+LIBIMAGES_LIBS=
+
+ifdef CONFIG_LIBJPEG
+LIBIMAGES_MODS+=io-libjpeg
+LIBIMAGES_LIBS+=-ljpeg
+endif
+
+ifdef CONFIG_LIBPNG
+LIBIMAGES_MODS+=io-libpng
+LIBIMAGES_LIBS+=-lpng
+endif
+
+ifdef CONFIG_LIBUNGIF
+LIBIMAGES_MODS+=io-libungif
+LIBIMAGES_LIBS+=-lungif
+endif
+
+ifdef CONFIG_LIBMAGICK
+LIBIMAGES_MODS+=io-libmagick
+MAGICK_LIBS:=$(shell GraphicsMagick-config --libs)
+LIBIMAGES_LIBS+=$(MAGICK_LIBS)
+$(o)/images/io-libmagick.o: CFLAGS+=-I/usr/include/GraphicsMagick
+endif
+
+$(o)/images/libimages.a: $(addsuffix .o,$(addprefix $(o)/images/,$(LIBIMAGES_MODS)))
+$(o)/images/libimages.so: $(addsuffix .oo,$(addprefix $(o)/images/,$(LIBIMAGES_MODS)))
+
+$(o)/images/image-tool: $(o)/images/image-tool.o $(LIBIMAGES) $(LIBUCW)
+$(o)/images/image-tool: LIBS+=$(LIBIMAGES_LIBS)
+
+TESTS+=$(o)/images/hilbert-test.test
+$(o)/images/hilbert-test: LIBS+=-lm $(LIBSH)
+$(o)/images/hilbert-test.test: $(o)/images/hilbert-test
+
+#$(o)/images/image-test: $(o)/images/image-test.o $(LIBIMAGES) $(LIBUCW)
+#$(o)/images/image-test: LIBS+=$(LIBIMAGES_LIBS)
+
+#$(o)/images/image-sig.o $(o)/images/image-sig.oo: CFLAGS+=-I/usr/include/GraphicsMagick
+#$(o)/images/image-idx.o $(o)/images/image-idx.oo: CFLAGS+=-I/usr/include/GraphicsMagick
+#$(o)/images/image-obj.o $(o)/images/image-obj.oo: CFLAGS+=-I/usr/include/GraphicsMagick
+#$(o)/images/image-idx: $(o)/images/image-idx.o $(o)/images/image-obj.o $(o)/images/dup-cmp.o $(o)/indexer/iconfig.o $(o)/images/image-sig.o $(o)/images/kd-tree.o $(o)/images/color.o $(o)/images/image-io.o $(LIBSH) $(LIBLANG) $(LIBCHARSET)
+#$(o)/images/image-idx: LIBS+=-lGraphicsMagick -ljpeg -lpng
+
+#$(o)/images/color-t: LIBS+=-lm
+#$(o)/images/color.test: $(o)/images/color-t
+
+# By :;DF
+#PROGS+=$(addprefix $(o)/images/,decomp)
+#
+#$(o)/images/block_info.o $(o)/images/block_info.oo: CFLAGS+=-I/usr/include/GraphicsMagick
+#$(o)/images/k_means.o $(o)/images/k_means.oo: CFLAGS+=-I/usr/include/GraphicsMagick
+#$(o)/images/decomp.o $(o)/images/decomp.oo: CFLAGS+=-I/usr/include/GraphicsMagick
+
+#$(o)/images/decomp: $(o)/images/decomp.o $(o)/images/block_info.o $(o)/images/k_means.o $(LIBSH) $(LIBLANG) $(LIBCHARSET)
+#$(o)/images/decomp: LIBS+=-lGraphicsMagick -ljpeg -lm
+
diff --git a/images/block_info.c b/images/block_info.c
new file mode 100644 (file)
index 0000000..a8de2bb
--- /dev/null
@@ -0,0 +1,146 @@
+/*#include <stdio.h>
+#include <magick/api.h>*/
+
+#include <assert.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "img.h"
+
+/*
+ * Color spaces
+ * 
+ * http://www.tecgraf.puc-rio.br/~mgattass/color/ColorIndex.html
+ * 
+ */
+
+#define REF_WHITE_X 0.96422
+#define REF_WHITE_Y 1.
+#define REF_WHITE_Z 0.82521
+
+/* sRGB to XYZ */
+static void
+srgb_to_xyz_slow(double srgb[3], double xyz[3])
+{
+  double a[3];
+  uns i;
+  for (i = 0; i < 3; i++)
+    if (srgb[i] > 0.04045)
+      a[i] = pow((srgb[i] + 0.055) * (1 / 1.055), 2.4);
+    else
+      a[i] = srgb[i] * (1 / 12.92);
+  xyz[0] = 0.412424 * a[0] + 0.357579 * a[1] + 0.180464 * a[2];
+  xyz[1] = 0.212656 * a[0] + 0.715158 * a[1] + 0.072186 * a[2];
+  xyz[2] = 0.019332 * a[0] + 0.119193 * a[1] + 0.950444 * a[2];
+}
+
+/* XYZ to CIE-Luv */
+static void
+xyz_to_luv_slow(double xyz[3], double luv[3])
+{
+  double sum = xyz[0] + 15 * xyz[1] + 3 * xyz[2];
+  if (sum < 0.000001)
+    luv[0] = luv[1] = luv[2] = 0;
+  else
+    {
+      double var_u = 4 * xyz[0] / sum;
+      double var_v = 9 * xyz[1] / sum;
+      if (xyz[1] > 0.008856)
+        luv[0] = 116 * pow(xyz[1], 1 / 3.) - 16;
+      else
+        luv[0] = (116 * 7.787) * xyz[1];
+      luv[1] = luv[0] * (13 * (var_u - 4 * REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
+      luv[2] = luv[0] * (13 * (var_v - 9 * REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
+      /* intervals [0..100], [-134..220], [-140..122] */
+    }
+}
+
+struct BlockInfo*
+computeBlockInfo(PixelPacket *pixels, uns width, uns height, uns *count)
+{
+  assert(width >= 4 && height >= 4);
+
+  uns w = width >> 2;
+  uns h = height >> 2;
+  fprintf(stderr, "Computing signature for image %dx%d... %dx%d blocks", width, height, w, h);
+  uns blocks_count = w * h;
+  struct BlockInfo *blocks = malloc(blocks_count * sizeof(struct BlockInfo)), *block = blocks; /* FIXME: use mempool */
+  
+  /* Every 4x4 block (FIXME: deal with smaller blocks near the edges) */
+  PixelPacket *p = pixels;
+  for (uns block_y = 0; block_y < h; block_y++, p += (width & 3) + 3*width){
+    for (uns block_x = 0; block_x < w; block_x++, p += 4 - 4*width, block++){
+        int t[16], s[16], *tp = t;
+        
+        /* Convert pixels to Luv color space and compute average coefficients 
+         * FIXME:
+         * - could be MUCH faster with precomputed tables and integer arithmetic... 
+         *   I will propably use interpolation in 3-dim array */
+        uns l_sum = 0;
+        uns u_sum = 0;
+        uns v_sum = 0;
+        for (uns y = 0; y < 4; y++, p += width - 4){
+          for (uns x = 0; x < 4; x++, p++)
+            {
+              double rgb[3], luv[3], xyz[3];
+              rgb[0] = (p->red >> (QuantumDepth - 8)) / 255.;
+              rgb[1] = (p->green >> (QuantumDepth - 8)) / 255.;
+              rgb[2] = (p->blue >> (QuantumDepth - 8)) / 255.;
+              srgb_to_xyz_slow(rgb, xyz);
+              xyz_to_luv_slow(xyz, luv);
+              l_sum += *tp++ = luv[0];
+              u_sum += luv[1] + 150;
+              v_sum += luv[2] + 150;
+              /*fprintf(stderr, "'%u, %u'; ", (p - pixels)%width ,  (p - pixels)/width);*/
+            }
+            /*fprintf(stderr, "\n---\n");*/
+        }
+        block->l = l_sum;
+        block->u = u_sum;
+        block->v = v_sum;
+        
+        /* Apply Daubechies wavelet transformation 
+         * FIXME:
+         * - MMX/SSE instructions or tables could be faster 
+         * - maybe it would be better to compute Luv and wavelet separately because of processor cache or MMX/SSE 
+         * - eliminate slow square roots 
+         * - what about Haar transformation? */
+
+#define DAUB_0  31651 /* (1 + sqrt 3) / (4 * sqrt 2) */
+#define DAUB_1  54822 /* (3 + sqrt 3) / (4 * sqrt 2) */
+#define DAUB_2  14689 /* (3 - sqrt 3) / (4 * sqrt 2) */
+#define DAUB_3  -8481 /* (1 - sqrt 3) / (4 * sqrt 2) */
+
+        /* ... to the rows */
+        uns i;
+        for (i = 0; i < 16; i += 4)
+          {
+            s[i + 0] = (DAUB_0 * t[i + 2] + DAUB_1 * t[i + 3] + DAUB_2 * t[i + 0] + DAUB_3 * t[i + 1]) / 0x10000;
+            s[i + 1] = (DAUB_0 * t[i + 0] + DAUB_1 * t[i + 1] + DAUB_2 * t[i + 2] + DAUB_3 * t[i + 3]) / 0x10000;
+            s[i + 2] = (DAUB_3 * t[i + 2] - DAUB_2 * t[i + 3] + DAUB_1 * t[i + 0] - DAUB_0 * t[i + 1]) / 0x10000;
+            s[i + 3] = (DAUB_3 * t[i + 0] - DAUB_2 * t[i + 1] + DAUB_1 * t[i + 2] - DAUB_0 * t[i + 3]) / 0x10000;
+          }
+
+        /* ... and to the columns... skip LL band */
+        for (i = 0; i < 2; i++)
+          {
+            t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
+            t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
+          }
+        for (; i < 4; i++)
+          {
+            t[i + 0] = (DAUB_0 * s[i + 8] + DAUB_1 * s[i +12] + DAUB_2 * s[i + 0] + DAUB_3 * s[i + 4]) / 0x1000;
+            t[i + 4] = (DAUB_0 * s[i + 0] + DAUB_1 * s[i + 4] + DAUB_2 * s[i + 8] + DAUB_3 * s[i +12]) / 0x1000;
+            t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
+            t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
+          }
+
+        /* Extract energies in LH, HL and HH bands */
+        block->lh = sqrt(t[8] * t[8] + t[9] * t[9] + t[12] * t[12] + t[13] * t[13]);
+        block->hl = sqrt(t[2] * t[2] + t[3] * t[3] + t[6] * t[6] + t[7] * t[7]);
+        block->hh = sqrt(t[10] * t[10] + t[11] * t[11] + t[14] * t[14] + t[15] * t[15]);
+      }
+  }
+  *count=blocks_count;
+  return blocks;
+}
diff --git a/images/color.c b/images/color.c
new file mode 100644 (file)
index 0000000..e684660
--- /dev/null
@@ -0,0 +1,353 @@
+/*
+ *     Image Library -- Color Spaces
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU Lesser General Public License.
+ */
+
+#undef LOCAL_DEBUG
+
+#include "sherlock/sherlock.h"
+#include "lib/math.h"
+#include "images/color.h"
+
+
+/********************* EXACT CONVERSION ROUTINES **********************/
+
+/* sRGB to XYZ */
+void
+srgb_to_xyz_slow(double xyz[3], double srgb[3])
+{
+  double a[3];
+  for (uns i = 0; i < 3; i++)
+    if (srgb[i] > 0.04045)
+      a[i] = pow((srgb[i] + 0.055) * (1 / 1.055), 2.4);
+    else
+      a[i] = srgb[i] * (1 / 12.92);
+  xyz[0] = SRGB_XYZ_XR * a[0] + SRGB_XYZ_XG * a[1] + SRGB_XYZ_XB * a[2];
+  xyz[1] = SRGB_XYZ_YR * a[0] + SRGB_XYZ_YG * a[1] + SRGB_XYZ_YB * a[2];
+  xyz[2] = SRGB_XYZ_ZR * a[0] + SRGB_XYZ_ZG * a[1] + SRGB_XYZ_ZB * a[2];
+}
+
+/* XYZ to CIE-Luv */
+void
+xyz_to_luv_slow(double luv[3], double xyz[3])
+{
+  double sum = xyz[0] + 15 * xyz[1] + 3 * xyz[2];
+  if (sum < 0.000001)
+    luv[0] = luv[1] = luv[2] = 0;
+  else
+    {
+      double var_u = 4 * xyz[0] / sum;
+      double var_v = 9 * xyz[1] / sum;
+      if (xyz[1] > 0.008856)
+        luv[0] = 116 * pow(xyz[1], 1 / 3.) - 16;
+      else
+        luv[0] = (116 * 7.787) * xyz[1];
+     luv[1] = luv[0] * (13 * (var_u - 4 * REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
+     luv[2] = luv[0] * (13 * (var_v - 9 * REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z)));
+     /* intervals [0..100], [-134..220], [-140..122] */
+   }
+}
+
+
+/***************** OPTIMIZED SRGB -> LUV CONVERSION *********************/
+
+u16 srgb_to_luv_tab1[256];
+u16 srgb_to_luv_tab2[9 << SRGB_TO_LUV_TAB2_SIZE];
+u32 srgb_to_luv_tab3[20 << SRGB_TO_LUV_TAB3_SIZE];
+
+void
+srgb_to_luv_init(void)
+{
+  DBG("Initializing sRGB -> Luv table");
+  for (uns i = 0; i < 256; i++)
+    {
+      double t = i / 255.;
+      if (t > 0.04045)
+        t = pow((t + 0.055) * (1 / 1.055), 2.4);
+      else
+        t = t * (1 / 12.92);
+      srgb_to_luv_tab1[i] = CLAMP(t * 0xfff + 0.5, 0, 0xfff);
+    }
+  for (uns i = 0; i < (9 << SRGB_TO_LUV_TAB2_SIZE); i++)
+    {
+      double t = i / (double)((9 << SRGB_TO_LUV_TAB2_SIZE) - 1);
+      if (t > 0.008856)
+        t = 1.16 * pow(t, 1 / 3.) - 0.16;
+      else
+        t = (1.16 * 7.787) * t;
+      srgb_to_luv_tab2[i] =
+       CLAMP(t * ((1 << SRGB_TO_LUV_TAB2_SCALE) - 1) + 0.5,
+          0, (1 << SRGB_TO_LUV_TAB2_SCALE) - 1);
+    }
+  for (uns i = 0; i < (20 << SRGB_TO_LUV_TAB3_SIZE); i++)
+    {
+      srgb_to_luv_tab3[i] = i ? (13 << (SRGB_TO_LUV_TAB3_SCALE + SRGB_TO_LUV_TAB3_SIZE)) / i : 0;
+    }
+}
+
+void
+srgb_to_luv_pixels(byte *dest, byte *src, uns count)
+{
+  while (count--)
+    {
+      srgb_to_luv_pixel(dest, src);
+      dest += 3;
+      src += 3;
+    }
+}
+
+
+/************************ GRID INTERPOLATION ALGORITHM ************************/
+
+struct color_grid_node *srgb_to_luv_grid;
+struct color_interpolation_node *color_interpolation_table;
+
+/* Returns volume of a given tetrahedron multiplied by 6 */
+static inline uns
+tetrahedron_volume(uns *v1, uns *v2, uns *v3, uns *v4)
+{
+  int a[3], b[3], c[3];
+  for (uns i = 0; i < 3; i++)
+    {
+      a[i] = v2[i] - v1[i];
+      b[i] = v3[i] - v1[i];
+      c[i] = v4[i] - v1[i];
+    }
+  int result =
+    a[0] * (b[1] * c[2] - b[2] * c[1]) -
+    a[1] * (b[0] * c[2] - b[2] * c[0]) +
+    a[2] * (b[0] * c[1] - b[1] * c[0]);
+  return (result > 0) ? result : -result;
+}
+
+static void
+interpolate_tetrahedron(struct color_interpolation_node *n, uns *p, const uns *c)
+{
+  uns v[4][3];
+  for (uns i = 0; i < 4; i++)
+    {
+      v[i][0] = (c[i] & 0001) ? (1 << COLOR_CONV_OFS) : 0;
+      v[i][1] = (c[i] & 0010) ? (1 << COLOR_CONV_OFS) : 0;
+      v[i][2] = (c[i] & 0100) ? (1 << COLOR_CONV_OFS) : 0;
+      n->ofs[i] =
+       ((c[i] & 0001) ? 1 : 0) +
+       ((c[i] & 0010) ? (1 << COLOR_CONV_SIZE) : 0) +
+       ((c[i] & 0100) ? (1 << (COLOR_CONV_SIZE * 2)) : 0);
+    }
+  uns vol = tetrahedron_volume(v[0], v[1], v[2], v[3]);
+  n->mul[0] = ((tetrahedron_volume(p, v[1], v[2], v[3]) << 8) + (vol >> 1)) / vol;
+  n->mul[1] = ((tetrahedron_volume(v[0], p, v[2], v[3]) << 8) + (vol >> 1)) / vol;
+  n->mul[2] = ((tetrahedron_volume(v[0], v[1], p, v[3]) << 8) + (vol >> 1)) / vol;
+  n->mul[3] = ((tetrahedron_volume(v[0], v[1], v[2], p) << 8) + (vol >> 1)) / vol;
+  uns j;
+  for (j = 0; j < 4; j++)
+    if (n->mul[j])
+      break;
+  for (uns i = 0; i < 4; i++)
+    if (n->mul[i] == 0)
+      n->ofs[i] = n->ofs[j];
+}
+
+static void
+interpolation_table_init(void)
+{
+  DBG("Initializing color interpolation table");
+  struct color_interpolation_node *n = color_interpolation_table =
+    xmalloc(sizeof(struct color_interpolation_node) << (COLOR_CONV_OFS * 3));
+  uns p[3];
+  for (p[2] = 0; p[2] < (1 << COLOR_CONV_OFS); p[2]++)
+    for (p[1] = 0; p[1] < (1 << COLOR_CONV_OFS); p[1]++)
+      for (p[0] = 0; p[0] < (1 << COLOR_CONV_OFS); p[0]++)
+        {
+         uns index;
+          static const uns tetrahedrons[5][4] = {
+            {0000, 0001, 0010, 0100},
+            {0110, 0111, 0100, 0010},
+            {0101, 0100, 0111, 0001},
+            {0011, 0010, 0001, 0111},
+            {0111, 0001, 0010, 0100}};
+         if (p[0] + p[1] + p[2] <= (1 << COLOR_CONV_OFS))
+           index = 0;
+         else if ((1 << COLOR_CONV_OFS) + p[0] <= p[1] + p[2])
+           index = 1;
+         else if ((1 << COLOR_CONV_OFS) + p[1] <= p[0] + p[2])
+           index = 2;
+         else if ((1 << COLOR_CONV_OFS) + p[2] <= p[0] + p[1])
+           index = 3;
+         else
+           index = 4;
+         interpolate_tetrahedron(n, p, tetrahedrons[index]);
+         n++;
+       }
+}
+
+typedef void color_conv_func(double dest[3], double src[3]);
+
+static void
+conv_grid_init(struct color_grid_node **grid, color_conv_func func)
+{
+  if (*grid)
+    return;
+  struct color_grid_node *g = *grid = xmalloc((sizeof(struct color_grid_node)) << (COLOR_CONV_SIZE * 3));
+  double src[3], dest[3];
+  for (uns k = 0; k < (1 << COLOR_CONV_SIZE); k++)
+    {
+      src[2] = k * (255 / (double)((1 << COLOR_CONV_SIZE) - 1));
+      for (uns j = 0; j < (1 << COLOR_CONV_SIZE); j++)
+        {
+          src[1] = j * (255/ (double)((1 << COLOR_CONV_SIZE) - 1));
+          for (uns i = 0; i < (1 << COLOR_CONV_SIZE); i++)
+            {
+              src[0] = i * (255 / (double)((1 << COLOR_CONV_SIZE) - 1));
+             func(dest, src);
+             g->val[0] = CLAMP(dest[0] + 0.5, 0, 255);
+             g->val[1] = CLAMP(dest[1] + 0.5, 0, 255);
+             g->val[2] = CLAMP(dest[2] + 0.5, 0, 255);
+             g++;
+           }
+       }
+    }
+}
+
+static void
+srgb_to_luv_func(double dest[3], double src[3])
+{
+  double srgb[3], xyz[3], luv[3];
+  srgb[0] = src[0] / 255.;
+  srgb[1] = src[1] / 255.;
+  srgb[2] = src[2] / 255.;
+  srgb_to_xyz_slow(xyz, srgb);
+  xyz_to_luv_slow(luv, xyz);
+  dest[0] = luv[0] * 2.55;
+  dest[1] = luv[1] * (2.55 / 4) + 128;
+  dest[2] = luv[2] * (2.55 / 4) + 128;
+}
+
+void
+color_conv_init(void)
+{
+  interpolation_table_init();
+  conv_grid_init(&srgb_to_luv_grid, srgb_to_luv_func);
+}
+
+void
+color_conv_pixels(byte *dest, byte *src, uns count, struct color_grid_node *grid)
+{
+  while (count--)
+    {
+      color_conv_pixel(dest, src, grid);
+      dest += 3;
+      src += 3;
+    }
+}
+
+
+/**************************** TESTS *******************************/
+
+#ifdef TEST
+#include <string.h>
+
+static double
+conv_error(u32 color, struct color_grid_node *grid, color_conv_func func)
+{
+  byte src[3], dest[3];
+  src[0] = color & 255;
+  src[1] = (color >> 8) & 255;
+  src[2] = (color >> 16) & 255;
+  color_conv_pixel(dest, src, grid);
+  double src2[3], dest2[3];
+  for (uns i = 0; i < 3; i++)
+    src2[i] = src[i];
+  func(dest2, src2);
+  double err = 0;
+  for (uns i = 0; i < 3; i++)
+    err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]);
+  return err;
+}
+
+typedef void test_fn(byte *dest, byte *src);
+
+static double
+func_error(u32 color, test_fn test, color_conv_func func)
+{
+  byte src[3], dest[3];
+  src[0] = color & 255;
+  src[1] = (color >> 8) & 255;
+  src[2] = (color >> 16) & 255;
+  test(dest, src);
+  double src2[3], dest2[3];
+  for (uns i = 0; i < 3; i++)
+    src2[i] = src[i];
+  func(dest2, src2);
+  double err = 0;
+  for (uns i = 0; i < 3; i++)
+    err += (dest[i] - dest2[i]) * (dest[i] - dest2[i]);
+  return err;
+}
+
+static void
+test_grid(byte *name, struct color_grid_node *grid, color_conv_func func)
+{
+  double max_err = 0, sum_err = 0;
+  uns count = 100000;
+  for (uns i = 0; i < count; i++)
+    {
+      double err = conv_error(random_max(0x1000000), grid, func);
+      max_err = MAX(err, max_err);
+      sum_err += err;
+    }
+  DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count);
+  if (max_err > 12)
+    die("Too large error in %s conversion", name);
+}
+
+static void
+test_func(byte *name, test_fn test, color_conv_func func)
+{
+  double max_err = 0, sum_err = 0;
+  uns count = 100000;
+  for (uns i = 0; i < count; i++)
+    {
+      double err = func_error(random_max(0x1000000), test, func);
+      max_err = MAX(err, max_err);
+      sum_err += err;
+    }
+  DBG("%s: error max=%f avg=%f", name, max_err, sum_err / count);
+  if (max_err > 12)
+    die("Too large error in %s conversion", name);
+}
+
+int
+main(void)
+{
+  srgb_to_luv_init();
+  test_func("func sRGB -> Luv", srgb_to_luv_pixel, srgb_to_luv_func);
+  color_conv_init();
+  test_grid("grid sRGB -> Luv", srgb_to_luv_grid, srgb_to_luv_func);
+#ifdef LOCAL_DEBUG
+#define CNT 1000000
+#define TESTS 10
+  byte *a = xmalloc(3 * CNT), *b = xmalloc(3 * CNT);
+  for (uns i = 0; i < 3 * CNT; i++)
+    a[i] = random_max(256);
+  init_timer();
+  for (uns i = 0; i < TESTS; i++)
+    memcpy(b, a, CNT * 3);
+  DBG("memcpy time=%d", (uns)get_timer());
+  init_timer();
+  for (uns i = 0; i < TESTS; i++)
+    srgb_to_luv_pixels(b, a, CNT);
+  DBG("direct time=%d", (uns)get_timer());
+  init_timer();
+  for (uns i = 0; i < TESTS; i++)
+    color_conv_pixels(b, a, CNT, srgb_to_luv_grid);
+  DBG("grid time=%d", (uns)get_timer());
+#endif
+  return 0;
+}
+#endif
+
diff --git a/images/color.h b/images/color.h
new file mode 100644 (file)
index 0000000..ec3cd4a
--- /dev/null
@@ -0,0 +1,143 @@
+/*
+ *     Image Library -- Color Spaces
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU Lesser General Public License.
+ *
+ *
+ *     References:
+ *     - http://www.tecgraf.puc-rio.br/~mgattass/color/ColorIndex.html
+ *
+ *     FIXME:
+ *     - fix theoretical problems with rounding errors in srgb_to_luv_pixel()
+ *     - SIMD should help to speed up conversion of large arrays
+ *     - maybe try to generate a long switch in color_conv_pixel()
+ *       with optimized entries instead of access to interpolation table
+ *     - most of multiplications in srgb_to_luv_pixels can be replaced
+ *       with tables lookup... tests shows almost the same speed for random
+ *       input and cca 40% gain when input colors fit in CPU chache
+ */
+
+#ifndef _IMAGES_COLOR_H
+#define _IMAGES_COLOR_H
+
+/* Exact slow conversion routines */
+void srgb_to_xyz_slow(double dest[3], double src[3]);
+void xyz_to_luv_slow(double dest[3], double src[3]);
+
+/* Reference white */
+#define REF_WHITE_X 0.96422
+#define REF_WHITE_Y 1.
+#define REF_WHITE_Z 0.82521
+
+/* sRGB -> XYZ matrix */
+#define SRGB_XYZ_XR 0.412424
+#define SRGB_XYZ_XG 0.357579
+#define SRGB_XYZ_XB 0.180464
+#define SRGB_XYZ_YR 0.212656
+#define SRGB_XYZ_YG 0.715158
+#define SRGB_XYZ_YB 0.072186
+#define SRGB_XYZ_ZR 0.019332
+#define SRGB_XYZ_ZG 0.119193
+#define SRGB_XYZ_ZB 0.950444
+
+
+/*********************** OPTIMIZED CONVERSION ROUTINES **********************/
+
+/* sRGB -> Luv parameters */
+#define SRGB_TO_LUV_TAB2_SIZE 9
+#define SRGB_TO_LUV_TAB2_SCALE 11
+#define SRGB_TO_LUV_TAB3_SIZE 8
+#define SRGB_TO_LUV_TAB3_SCALE (39 - SRGB_TO_LUV_TAB2_SCALE - SRGB_TO_LUV_TAB3_SIZE)
+
+extern u16 srgb_to_luv_tab1[256];
+extern u16 srgb_to_luv_tab2[9 << SRGB_TO_LUV_TAB2_SIZE];
+extern u32 srgb_to_luv_tab3[20 << SRGB_TO_LUV_TAB3_SIZE];
+
+void srgb_to_luv_init(void);
+void srgb_to_luv_pixels(byte *dest, byte *src, uns count);
+
+/* L covers the interval [0..255]; u and v are centered to 128 and scaled by 1/4 in respect of L */
+static inline void
+srgb_to_luv_pixel(byte *dest, byte *src)
+{
+  uns r = srgb_to_luv_tab1[src[0]];
+  uns g = srgb_to_luv_tab1[src[1]];
+  uns b = srgb_to_luv_tab1[src[2]];
+  uns x =
+    (uns)(4 * SRGB_XYZ_XR * 0xffff) * r +
+    (uns)(4 * SRGB_XYZ_XG * 0xffff) * g +
+    (uns)(4 * SRGB_XYZ_XB * 0xffff) * b;
+  uns y =
+    (uns)(9 * SRGB_XYZ_YR * 0xffff) * r +
+    (uns)(9 * SRGB_XYZ_YG * 0xffff) * g +
+    (uns)(9 * SRGB_XYZ_YB * 0xffff) * b;
+  uns l = srgb_to_luv_tab2[y >> (28 - SRGB_TO_LUV_TAB2_SIZE)];
+    dest[0] = l >> (SRGB_TO_LUV_TAB2_SCALE - 8);
+  uns sum =
+    (uns)((SRGB_XYZ_XR + 15 * SRGB_XYZ_YR + 3 * SRGB_XYZ_ZR) * 0x7fff) * r +
+    (uns)((SRGB_XYZ_XG + 15 * SRGB_XYZ_YG + 3 * SRGB_XYZ_ZG) * 0x7fff) * g +
+    (uns)((SRGB_XYZ_XB + 15 * SRGB_XYZ_YB + 3 * SRGB_XYZ_ZB) * 0x7fff) * b;
+  uns s = srgb_to_luv_tab3[sum >> (27 - SRGB_TO_LUV_TAB3_SIZE)];
+  int xs = ((u64)x * s) >> 32;
+  int ys = ((u64)y * s) >> 32;
+  int xw = ((4 * 13) << (SRGB_TO_LUV_TAB3_SCALE - 4)) *
+    REF_WHITE_X / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z);
+  int yw = ((9 * 13) << (SRGB_TO_LUV_TAB3_SCALE - 4)) *
+    REF_WHITE_Y / (REF_WHITE_X + 15 * REF_WHITE_Y + 3 * REF_WHITE_Z);
+  int u = (int)(l) * (xs - xw);
+  int v = (int)(l) * (ys - yw);
+  dest[1] = 128 + (u >> (SRGB_TO_LUV_TAB3_SCALE + SRGB_TO_LUV_TAB2_SCALE - 10));
+  dest[2] = 128 + (v >> (SRGB_TO_LUV_TAB3_SCALE + SRGB_TO_LUV_TAB2_SCALE - 10));
+}
+
+
+/****************** GENERAL INTERPOLATION IN 3D GRID ********************/
+
+#define COLOR_CONV_SIZE        5  /* 128K conversion grid size */
+#define COLOR_CONV_OFS 3  /* 8K interpolation table size */
+
+struct color_grid_node {
+  byte val[4];
+};
+
+struct color_interpolation_node {
+  u16 ofs[4];
+  u16 mul[4];
+};
+
+extern struct color_grid_node *srgb_to_luv_grid;
+extern struct color_interpolation_node *color_interpolation_table;
+
+void color_conv_init(void);
+void color_conv_pixels(byte *dest, byte *src, uns count, struct color_grid_node *grid);
+
+#define COLOR_CONV_SCALE_CONST (((((1 << COLOR_CONV_SIZE) - 1) << 16) + (1 << (16 - COLOR_CONV_OFS))) / 255)
+
+static inline void
+color_conv_pixel(byte *dest, byte *src, struct color_grid_node *grid)
+{
+  uns s0 = src[0] * COLOR_CONV_SCALE_CONST;
+  uns s1 = src[1] * COLOR_CONV_SCALE_CONST;
+  uns s2 = src[2] * COLOR_CONV_SCALE_CONST;
+  struct color_grid_node *g0, *g1, *g2, *g3, *g = grid +
+    ((s0 >> 16) + ((s1 >> 16) << COLOR_CONV_SIZE) + ((s2 >> 16) << (2 * COLOR_CONV_SIZE)));
+  struct color_interpolation_node *n = color_interpolation_table +
+    (((s0 & (0x10000 - (0x10000 >> COLOR_CONV_OFS))) >> (16 - COLOR_CONV_OFS)) +
+    ((s1 & (0x10000 - (0x10000 >> COLOR_CONV_OFS))) >> (16 - 2 * COLOR_CONV_OFS)) +
+    ((s2 & (0x10000 - (0x10000 >> COLOR_CONV_OFS))) >> (16 - 3 * COLOR_CONV_OFS)));
+  g0 = g + n->ofs[0];
+  g1 = g + n->ofs[1];
+  g2 = g + n->ofs[2];
+  g3 = g + n->ofs[3];
+  dest[0] = (g0->val[0] * n->mul[0] + g1->val[0] * n->mul[1] +
+             g2->val[0] * n->mul[2] + g3->val[0] * n->mul[3] + 128) >> 8;
+  dest[1] = (g0->val[1] * n->mul[0] + g1->val[1] * n->mul[1] +
+             g2->val[1] * n->mul[2] + g3->val[1] * n->mul[3] + 128) >> 8;
+  dest[2] = (g0->val[2] * n->mul[0] + g1->val[2] * n->mul[1] +
+             g2->val[2] * n->mul[2] + g3->val[2] * n->mul[3] + 128) >> 8;
+}
+
+#endif
diff --git a/images/color.t b/images/color.t
new file mode 100644 (file)
index 0000000..4259608
--- /dev/null
@@ -0,0 +1,4 @@
+# Tests for color conversion module
+
+Run:   obj/images/color-t
+
diff --git a/images/decomp.c b/images/decomp.c
new file mode 100644 (file)
index 0000000..1659efa
--- /dev/null
@@ -0,0 +1,73 @@
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include <alloca.h>
+#include "lib/config.h"
+#include "img.h"
+
+int
+main(int argc, char *argv[]){
+  if(argc<2)
+    return 0;
+  
+  u8 retval=0;
+  struct DecomposeImageInfo dii;
+  ExceptionInfo exception;
+  Image *image=NULL;
+  ImageInfo *image_info=NULL;
+  struct BlockInfo *bi=NULL;
+
+  Image *out_image=NULL;
+
+  InitializeMagick(NULL);
+  GetExceptionInfo(&exception);
+  
+  image_info=CloneImageInfo((ImageInfo *) NULL);
+  (void) strcpy(image_info->filename, argv[1]);
+  image=ReadImage(image_info,&exception);
+  if(exception.severity != UndefinedException)
+    CatchException(&exception);
+
+  if(!image){
+    fprintf(stderr, "Invalid image format");
+    goto exit;
+  }
+  if(image->columns < 4 || image->rows < 4){
+    fprintf(stderr, "Image too small (%dx%d)", (int)image->columns, (int)image->rows);
+    retval = -1;
+    goto exit;
+  }
+
+  QuantizeInfo quantize_info;
+  GetQuantizeInfo(&quantize_info);
+  quantize_info.colorspace = RGBColorspace;
+  QuantizeImage(&quantize_info, image);
+
+  PixelPacket *pixels = (PixelPacket *) AcquireImagePixels(image, 0, 0, image->columns, image->rows, &exception);
+  if (exception.severity != UndefinedException) CatchException(&exception);
+  bi=computeBlockInfo(pixels, image->columns, image->rows, &dii.bi_len);
+  
+  dii.max_cls_num=16;
+  //dii.threshold=100;
+  //dii.diff_threshold=1000;
+  dii.max_cycles=7;
+  dii.init_decomp_num=5;
+  
+  decomposeImage(&dii, bi);
+
+  showBlockInfoAsImage(bi, dii.bi_len, image->columns, image->rows, &out_image, &exception);
+  if (exception.severity != UndefinedException) CatchException(&exception);
+  
+  image_info=CloneImageInfo((ImageInfo *) NULL);
+  strcpy(out_image->filename, "/proc/self/fd/1"); /*out_image->file=stdout did'n work for me*/
+  out_image->compression=JPEGCompression;
+  if(WriteImage(image_info, out_image)==0)
+    CatchException(&out_image->exception);
+exit:
+  DestroyImage(out_image);
+  DestroyImage(image);
+  DestroyImageInfo(image_info);
+  DestroyExceptionInfo(&exception);
+  DestroyMagick();
+  return retval;
+}
diff --git a/images/dup-cmp.c b/images/dup-cmp.c
new file mode 100644 (file)
index 0000000..fbb2de8
--- /dev/null
@@ -0,0 +1,366 @@
+/*
+ *      Image Library -- Duplicates Comparison
+ *
+ *      (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *      This software may be freely distributed and used according to the terms
+ *      of the GNU Lesser General Public License.
+ *
+ *
+ *      FIXME:
+ *      - many possible optimization
+ *      - compare normalized pictures (brightness, ...)
+ *      - better image scale... now it can completely miss some rows/cols of pixels
+ *      - maybe better/slower last step
+ *      - different thresholds for various transformations
+ *      - do not test all transformations for symetric pictures
+ *      - allocated memory could be easily decreased to about 1/3 
+ *        for aspect ratio threshold near one
+ *      - ... secret ideas :-)
+ */
+
+#undef LOCAL_DEBUG
+
+#include "sherlock/sherlock.h"
+#include "lib/mempool.h"
+#include "images/images.h"
+#include "images/dup-cmp.h"
+
+static uns image_dup_scale_min_size = 16;
+static uns image_dup_ratio_threshold = 140;
+static uns image_dup_error_threshold = 50;
+
+static inline byte *
+image_dup_block(struct image_dup *dup, uns col, uns row)
+{
+  ASSERT(col <= dup->cols && row <= dup->rows);
+  return dup->buf + (dup->line << row) + (3 << (row + col));
+}
+
+static inline void
+pixels_average(byte *dest, byte *src1, byte *src2)
+{
+  dest[0] = ((uns)src1[0] + (uns)src2[0]) >> 1;
+  dest[1] = ((uns)src1[1] + (uns)src2[1]) >> 1;
+  dest[2] = ((uns)src1[2] + (uns)src2[2]) >> 1;
+}
+
+uns
+image_dup_estimate_size(uns width, uns height)
+{
+  uns cols, rows;
+  for (cols = 0; (uns)(2 << cols) < width; cols++);
+  for (rows = 0; (uns)(2 << rows) < height; rows++);
+  return sizeof(struct image_dup) + (12 << (cols + rows));
+}
+
+void
+image_dup_init(struct image_dup *dup, struct image_data *image, struct mempool *pool)
+{
+  ASSERT(image->width && image->height);
+  
+  dup->image = image;
+  dup->width = image->width;
+  dup->height = image->height;
+  for (dup->cols = 0; (uns)(2 << dup->cols) < image->width; dup->cols++);
+  for (dup->rows = 0; (uns)(2 << dup->rows) < image->height; dup->rows++);
+  dup->buf = mp_alloc(pool, dup->buf_size = (12 << (dup->cols + dup->rows)));
+  dup->line = 6 << dup->cols;
+  dup->flags = 0;
+  if (image->width >= image_dup_scale_min_size && image->height >= image_dup_scale_min_size)
+    dup->flags |= IMAGE_DUP_FLAG_SCALE;
+  
+  /* Scale original image to right bottom block */
+  {
+    byte *d = image_dup_block(dup, dup->cols, dup->rows);
+    uns width = 1 << dup->cols;
+    uns height = 1 << dup->rows;
+    uns line_size = 3 * image->width;
+    uns src_y = 0;
+    for (uns y = 0; y < height; y++)
+      {
+       byte *line = image->pixels + line_size * (src_y >> dup->rows);
+        uns src_x = 0;
+        for (uns x = 0; x < width; x++)
+          {
+           byte *s = line + 3 * (src_x >> dup->cols);
+           d[0] = s[0];
+           d[1] = s[1];
+           d[2] = s[2];
+           d += 3;
+           src_x += image->width;
+         }
+       src_y += image->height;
+      }
+  }
+
+  /* Complete bottom row */
+  for (uns i = dup->cols; i--; )
+    {
+      byte *d = image_dup_block(dup, i, dup->rows);
+      byte *s = image_dup_block(dup, i + 1, dup->rows);
+      for (uns y = 0; y < (uns)(1 << dup->rows); y++)
+       for (uns x = 0; x < (uns)(1 << i); x++)
+         {
+           pixels_average(d, s, s + 3);
+           d += 3;
+           s += 6;
+         }
+    }
+  /* Complete remaining blocks */
+  for (uns i = 0; i <= dup->cols; i++)
+    {
+      uns line_size = (3 << i);
+      for (uns j = dup->rows; j--; )
+        {
+          byte *d = image_dup_block(dup, i, j);
+          byte *s = image_dup_block(dup, i, j + 1);
+          for (uns y = 0; y < (uns)(1 << j); y++)
+            {
+              for (uns x = 0; x < (uns)(1 << i); x++)
+                {
+                 pixels_average(d, s, s + line_size);
+                 d += 3;
+                 s += 3;
+               }
+             s += line_size;
+           }
+        }
+    }
+}
+
+static inline uns
+err (int a, int b)
+{
+  a -= b;
+  return a * a;
+}
+
+static inline uns
+err_sum(byte *pos1, byte *end1, byte *pos2)
+{
+  uns e = 0;
+  while (pos1 != end1)
+    e += err(*pos1++, *pos2++);
+  return e;
+}
+
+static inline uns
+err_sum_transformed(byte *pos1, byte *end1, byte *pos2, uns width, int add1, int add2)
+{
+  DBG("err_sum_transformed(): %p %p %p %d %d %d", pos1, end1, pos2, width, add1, add2);
+  uns e = 0;
+  while (pos1 != end1)
+    {
+      for (uns i = 0; i < width; i++, pos2 += add1)
+      {
+       e += err(pos1[0], pos2[0]);
+       e += err(pos1[1], pos2[1]);
+       e += err(pos1[2], pos2[2]);
+       pos1 += 3;
+      }
+      pos2 += add2;
+    }
+  return e;
+}
+
+static inline int
+aspect_ratio_test(uns width1, uns height1, uns width2, uns height2)
+{
+  uns r1 = width1 * height2;
+  uns r2 = height1 * width2;
+  return
+    r1 <= ((r2 * image_dup_ratio_threshold) >> 5) && 
+    r2 <= ((r1 * image_dup_ratio_threshold) >> 5);
+}
+
+static inline int
+average_compare(struct image_dup *dup1, struct image_dup *dup2)
+{
+  byte *block1 = image_dup_block(dup1, 0, 0);
+  byte *block2 = image_dup_block(dup2, 0, 0);
+  uns e =
+    err(block1[0], block2[0]) +
+    err(block1[1], block2[1]) +
+    err(block1[2], block2[2]);
+  return e <= image_dup_error_threshold;
+}
+
+static int
+blocks_compare(struct image_dup *dup1, struct image_dup *dup2, uns col, uns row, uns trans)
+{
+  DBG("blocks_compare(): col=%d row=%d trans=%d", col, row, trans);
+  byte *block1 = image_dup_block(dup1, col, row);
+  byte *block2 = (trans < 4) ? image_dup_block(dup2, col, row) : image_dup_block(dup2, row, col);
+  int add1, add2;
+  switch (trans)
+    {
+      case 0: ;
+       uns err = (err_sum(block1, block1 + (3 << (col + row)), block2) >> (col + row));
+       DBG("average error=%d", err);
+       return err <= image_dup_error_threshold;
+      case 1:
+       add1 = -3;
+       add2 = 6 << col;
+       block2 += (3 << col) - 3;
+       break;
+      case 2:
+       add1 = 1;
+       add2 = -(6 << col);
+       block2 += (3 << (col + row)) - (3 << col);
+       break;
+      case 3:
+       add1 = -3;
+       add2 = 0;
+       block2 += (3 << (col + row)) - 3;
+       break;
+      case 4:
+       add1 = (3 << col);
+       add2 = -(3 << (col + row)) + 3;
+       break;
+      case 5:
+       add1 = -(3 << col);
+       add2 = (3 << (col + row)) + 3;
+       block2 += (3 << (col + row)) - (3 << col);
+       break;
+      case 6:
+       add1 = (3 << col);
+       add2 = -(3 << (col + row)) - 3;
+       block2 += (3 << col) - 3;
+       break;
+      case 7:
+       add1 = -(3 << col);
+       add2 = (3 << (col + row)) - 3;
+       block2 += (3 << (col + row)) - 3;
+       break;
+      default:
+       ASSERT(0);
+    }
+  uns err = (err_sum_transformed(block1, block1 + (3 << (col + row)), block2, (1 << col), add1, add2) >> (col + row));
+  DBG("average error=%d", err);
+  return err <= image_dup_error_threshold;
+}
+
+static int
+same_size_compare(struct image_dup *dup1, struct image_dup *dup2, uns trans)
+{
+  byte *block1 = dup1->image->pixels;
+  byte *block2 = dup2->image->pixels;
+  DBG("same_size_compare(): trans=%d",  trans);
+  int add1, add2;
+  switch (trans)
+    {
+      case 0: ;
+        uns err = (err_sum(block1, block1 + 3 * dup1->width * dup1->height, block2) / (dup1->width * dup1->height));
+       DBG("average error=%d", err);
+       return err <= image_dup_error_threshold;
+      case 1:
+       add1 = -3;
+       add2 = 6 * dup1->width;
+       block2 += 3 * (dup1->width - 1);
+       break;
+      case 2:
+       add1 = 1;
+       add2 = -6 * dup1->width;
+       block2 += 3 * dup1->width * (dup1->height - 1);
+       break;
+      case 3:
+       add1 = -3;
+       add2 = 0;
+       block2 += 3 * (dup1->width * dup1->height - 1);
+       break;
+      case 4:
+       add1 = 3 * dup1->width;
+       add2 = -3 * (dup1->width * dup1->height - 1);
+       break;
+      case 5:
+       add1 = -3 * dup1->width;
+       add2 = 3 * (dup1->width * dup1->height + 1);
+       block2 += 3 * dup1->width * (dup1->height - 1);
+       break;
+      case 6:
+       add1 = 3 * dup1->width;
+       add2 = -3 * (dup1->width * dup1->height + 1);
+       block2 += 3 * (dup1->width - 1);
+       break;
+      case 7:
+       add1 = -3 * dup1->width;
+       add2 = 3 * (dup1->width * dup1->height - 1);
+       block2 += 3 * (dup1->width * dup1->height - 1);
+       break;
+      default:
+       ASSERT(0);
+    }
+  uns err = (err_sum_transformed(block1, block1 + 3 * dup1->width * dup1->height, block2, dup1->width, add1, add2) / (dup1->width * dup1->height));
+  DBG("average error=%d", err);
+  return err <= image_dup_error_threshold;
+}
+
+int
+image_dup_compare(struct image_dup *dup1, struct image_dup *dup2, uns trans)
+{
+  if (!average_compare(dup1, dup2))
+    return 0;
+  if ((dup1->flags & dup2->flags) & IMAGE_DUP_FLAG_SCALE)
+    {
+      DBG("Scale support");
+      if (!aspect_ratio_test(dup1->width, dup1->height, dup2->width, dup2->height))
+       trans &= 0xf0;
+      if (!aspect_ratio_test(dup1->width, dup1->height, dup2->height, dup2->width))
+       trans &= 0x0f;
+    }
+  else
+    {
+      DBG("No scale support");
+      if (!(dup1->width == dup2->width && dup1->height == dup2->height))
+       trans &= 0xf0;
+      if (!(dup1->width == dup2->height && dup1->height == dup2->width))
+       trans &= 0x0f;
+    }
+  if (!trans)
+    return 0;
+  if (trans & 0x0f)
+    {
+      uns cols = MIN(dup1->cols, dup2->cols);
+      uns rows = MIN(dup1->rows, dup2->rows);
+      for (uns t = 0; t < 4; t++)
+       if (trans & (1 << t))
+         {
+           DBG("Testing trans %d", t);
+            for (uns i = MAX(cols, rows); i--; )
+              {
+               uns col = MAX(0, (int)(cols - i));
+               uns row = MAX(0, (int)(rows - i));
+               if (!blocks_compare(dup1, dup2, col, row, t))
+                 break;
+               if (!i &&
+                   (dup1->width != dup2->width || dup1->height != dup2->height ||
+                   same_size_compare(dup1, dup2, t)))
+                 return 1;
+             }
+         }
+    }
+  if (trans & 0xf0)
+    {
+      uns cols = MIN(dup1->cols, dup2->rows);
+      uns rows = MIN(dup1->rows, dup2->cols);
+      for (uns t = 4; t < 8; t++)
+       if (trans & (1 << t))
+         {
+           DBG("Testing trans %d", t);
+            for (uns i = MAX(cols, rows); i--; )
+              {
+               uns col = MAX(0, (int)(cols - i));
+               uns row = MAX(0, (int)(rows - i));
+               if (!blocks_compare(dup1, dup2, col, row, t))
+                 break;
+               if (!i &&
+                   (dup1->width != dup2->height || dup1->height != dup2->width ||
+                   same_size_compare(dup1, dup2, t)) )
+                 return 1;
+             }
+         }
+    }
+  return 0;
+}
diff --git a/images/dup-cmp.h b/images/dup-cmp.h
new file mode 100644 (file)
index 0000000..03d477d
--- /dev/null
@@ -0,0 +1,27 @@
+#ifndef _IMAGES_DUP_CMP_H
+#define _IMAGES_DUP_CMP_H
+
+struct image_data;
+
+struct image_dup {
+  struct image *image;
+  byte *buf;
+  uns buf_size;
+  uns flags;
+  uns cols;
+  uns rows;
+  uns line;
+  uns width;
+  uns height;
+};
+
+#define IMAGE_DUP_FLAG_SCALE   0x1
+
+#define IMAGE_DUP_TRANS_ID     0x01
+#define IMAGE_DUP_TRANS_ALL    0xff
+
+void image_dup_init(struct image_dup *dup, struct image *image, struct mempool *pool);
+int image_dup_compare(struct image_dup *dup1, struct image_dup *dup2, uns trans);
+uns image_dup_estimate_size(uns width, uns height);
+
+#endif
diff --git a/images/hilbert-test.c b/images/hilbert-test.c
new file mode 100644 (file)
index 0000000..9824f09
--- /dev/null
@@ -0,0 +1,124 @@
+/* Tests for multidimensional Hilbert curves */
+
+#define LOCAL_DEBUG
+
+#include "lib/lib.h"
+#include "lib/mempool.h"
+#include "lib/math.h"
+#include <stdlib.h>
+#include <stdio.h>
+
+static struct mempool *pool;
+
+static uns dim;
+static uns order;
+
+static inline void
+rand_vec(uns *vec)
+{
+  for (uns i = 0; i < dim; i++)
+    vec[i] = (uns)rand() >> (32 - order);
+}
+
+static byte *
+print_vec(uns *vec)
+{
+  byte *s = mp_alloc(pool, dim * 16), *res = s;
+  *s++ = '(';
+  for (uns i = 0; i < dim; i++)
+    {
+      if (i)
+       *s++ = ' ';
+      s += sprintf(s, "%x", vec[i]);
+    }
+  *s++ = ')';
+  *s = 0;
+  return res;
+}
+
+static inline int
+cmp_vec(uns *vec1, uns *vec2)
+{
+  for (uns i = dim; i--; )
+    if (vec1[i] < vec2[i])
+      return -1;
+    else if (vec1[i] > vec2[i])
+      return 1;
+  return 0;
+}
+
+#if 0
+static long double
+param_dist(uns *vec1, uns *vec2)
+{
+  long double d1 = 0, d2 = 0;
+  for (uns i = 0; i < dim; i++)
+    {
+      d1 = (d1 + vec1[i]) / ((u64)1 << order);
+      d2 = (d2 + vec2[i]) / ((u64)1 << order);
+    }
+  return fabsl(d1 - d2);
+}
+
+static long double
+vec_dist(uns *vec1, uns *vec2)
+{
+  long double d = 0;
+  for (uns i = 0; i < dim; i++)
+    {
+      long double x = fabsl(vec1[i] - vec2[i]) / ((u64)1 << order);
+      d += x * x;
+    }
+  return sqrtl(d);
+}
+#endif
+
+#define HILBERT_PREFIX(x) test1_##x
+#define HILBERT_DIM dim
+#define HILBERT_ORDER order
+#define HILBERT_WANT_DECODE
+#define HILBERT_WANT_ENCODE
+#include "images/hilbert.h"
+
+static void
+test1(void)
+{
+  uns a[32], b[32], c[32];
+  for (dim = 2; dim <= 8; dim++)
+    for (order = 8; order <= 32; order++)
+      for (uns i = 0; i < 1000; i++)
+        {
+         rand_vec(a);
+          test1_encode(b, a);
+          test1_decode(c, b);
+         if (cmp_vec(a, c))
+           die("Error... dim=%d order=%d testnum=%d ... %s -> %s -> %s", 
+               dim, order, i, print_vec(a), print_vec(b), print_vec(c));
+        }
+}
+
+#if 0
+#include "images/hilbert-origin.h"
+static void
+test_origin(void)
+{
+  Hcode code;
+  Point pt, pt2;
+  pt.hcode[0] = 0x12345678;
+  pt.hcode[1] = 0x654321;
+  pt.hcode[2] = 0x11122233;
+  code = H_encode(pt);
+  pt2 = H_decode(code);
+  DBG("origin: [%08x, %08x, %08x] --> [%08x, %08x %08x] --> [%08x, %08x %08x]", 
+    pt.hcode[0], pt.hcode[1], pt.hcode[2], code.hcode[0], code.hcode[1], code.hcode[2], pt2.hcode[0], pt2.hcode[1], pt2.hcode[2]);
+}
+#endif
+
+int
+main(int argc UNUSED, char **argv UNUSED)
+{
+  pool = mp_new(1 << 16);
+  test1();
+  //test_origin();
+  return 0;
+}
diff --git a/images/hilbert-test.t b/images/hilbert-test.t
new file mode 100644 (file)
index 0000000..f8794e3
--- /dev/null
@@ -0,0 +1,3 @@
+# Tests for multidimensional Hilbert curves
+
+Run:   obj/images/hilbert-test
diff --git a/images/hilbert.h b/images/hilbert.h
new file mode 100644 (file)
index 0000000..4971826
--- /dev/null
@@ -0,0 +1,340 @@
+/*
+ *     Image Library -- multidimensional Hilbert curves
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU Lesser General Public License.
+ *
+ *
+ *     References:
+ *     - http://www.dcs.bbk.ac.uk/~jkl/mapping.c
+ *       (c) 2002 J.K.Lawder
+ *     - J.K. Lawder. Calculation of Mappings between One and n-dimensional Values
+ *       Using the Hilbert Space-Filling Curve. Technical Report JL1/00, Birkbeck
+ *       College, University of London, 2000.
+ *
+ *     FIXME:
+ *     - the algorithm fails for some combinations of HILBERT_DIM and HILBERT_ORDER,
+ *       but it should be safe for HILBERT_DIM = 2..8, HILBERT_ORDER = 8..32
+ *     - clean and optimize the code
+ */
+
+#ifndef HILBERT_PREFIX
+#  error Undefined HILBERT_PREFIX
+#endif
+
+#define P(x) HILBERT_PREFIX(x)
+
+/*
+ * HILBERT_DIM is the number of dimensions in space through which the
+ * Hilbert Curve passes.
+ * Don't use this implementation with values for HILBERT_DIM of > 31!
+ * Also, make sure you use a 32 bit compiler!
+ */
+#ifndef HILBERT_DIM
+#  define HILBERT_DIM 2
+#endif
+
+#ifndef HILBERT_TYPE
+#  define HILBERT_TYPE u32
+#endif
+
+#ifndef HILBERT_ORDER
+#  define HILBERT_ORDER (8 * sizeof(HILBERT_TYPE))
+#endif
+
+typedef HILBERT_TYPE P(t);
+
+/*
+ * retained for historical reasons: the number of bits in an attribute value:
+ * effectively the order of a curve
+ */
+#define NUMBITS HILBERT_ORDER
+
+/*
+ * the number of bits in a word used to store an hcode (or in an element of
+ * an array that's used)
+ */
+#define WORDBITS HILBERT_ORDER
+
+#ifdef HILBERT_WANT_ENCODE
+/*
+ * given the coordinates of a point, it finds the sequence number of the point
+ * on the Hilbert Curve
+ */
+static void
+P(encode) (P(t) *dest, P(t) *src)
+{
+  P(t) mask = (P(t))1 << WORDBITS - 1, element, temp1, temp2,
+    A, W = 0, S, tS, T, tT, J, P = 0, xJ;
+  uns i = NUMBITS * HILBERT_DIM - HILBERT_DIM, j;
+
+  for (j = 0; j < HILBERT_DIM; j++)
+    dest[j] = 0;
+  for (j = A = 0; j < HILBERT_DIM; j++)
+    if (src[j] & mask)
+      A |= (1 << HILBERT_DIM - 1 - j);
+
+  S = tS = A;
+
+  P |= S & (1 << HILBERT_DIM - 1);
+  for (j = 1; j < HILBERT_DIM; j++)
+    if( S & (1 << HILBERT_DIM - 1 - j) ^ (P >> 1) & (1 << HILBERT_DIM - 1 - j))
+      P |= (1 << HILBERT_DIM - 1 - j);
+
+  /* add in HILBERT_DIM bits to hcode */
+  element = i / WORDBITS;
+  if (i % WORDBITS > WORDBITS - HILBERT_DIM)
+    {
+      dest[element] |= P << i % WORDBITS;
+      dest[element + 1] |= P >> WORDBITS - i % WORDBITS;
+    }
+  else
+    dest[element] |= P << i - element * WORDBITS;
+
+  J = HILBERT_DIM;
+  for (j = 1; j < HILBERT_DIM; j++)
+    if ((P >> j & 1) == (P & 1))
+      continue;
+    else
+      break;
+  if (j != HILBERT_DIM)
+    J -= j;
+  xJ = J - 1;
+
+  if (P < 3)
+    T = 0;
+  else
+    if (P % 2)
+      T = (P - 1) ^ (P - 1) / 2;
+    else
+      T = (P - 2) ^ (P - 2) / 2;
+  tT = T;
+
+  for (i -= HILBERT_DIM, mask >>= 1; (int)i >= 0; i -= HILBERT_DIM, mask >>= 1)
+    {
+      for (j = A = 0; j < HILBERT_DIM; j++)
+        if (src[j] & mask)
+          A |= (1 << HILBERT_DIM - 1 - j);
+
+      W ^= tT;
+      tS = A ^ W;
+      if (xJ % HILBERT_DIM != 0)
+        {
+          temp1 = tS << xJ % HILBERT_DIM;
+          temp2 = tS >> HILBERT_DIM - xJ % HILBERT_DIM;
+         S = temp1 | temp2;
+         S &= ((P(t))1 << HILBERT_DIM) - 1;
+       }
+      else
+       S = tS;
+
+      P = S & (1 << HILBERT_DIM - 1);
+      for (j = 1; j < HILBERT_DIM; j++)
+       if( S & (1 << HILBERT_DIM - 1 - j) ^ (P >> 1) & (1 << HILBERT_DIM - 1 - j))
+         P |= (1 << HILBERT_DIM - 1 - j);
+
+      /* add in HILBERT_DIM bits to hcode */
+      element = i / WORDBITS;
+      if (i % WORDBITS > WORDBITS - HILBERT_DIM)
+       {
+         dest[element] |= P << i % WORDBITS;
+         dest[element + 1] |= P >> WORDBITS - i % WORDBITS;
+       }
+      else
+       dest[element] |= P << i - element * WORDBITS;
+
+      if (i > 0)
+       {
+         if (P < 3)
+           T = 0;
+         else
+           if (P % 2)
+             T = (P - 1) ^ (P - 1) / 2;
+           else
+             T = (P - 2) ^ (P - 2) / 2;
+
+         if (xJ % HILBERT_DIM != 0)
+           {
+             temp1 = T >> xJ % HILBERT_DIM;
+             temp2 = T << HILBERT_DIM - xJ % HILBERT_DIM;
+             tT = temp1 | temp2;
+             tT &= ((P(t))1 << HILBERT_DIM) - 1;
+           }
+         else
+           tT = T;
+
+         J = HILBERT_DIM;
+         for (j = 1; j < HILBERT_DIM; j++)
+           if ((P >> j & 1) == (P & 1))
+             continue;
+           else
+             break;
+         if (j != HILBERT_DIM)
+           J -= j;
+
+         xJ += J - 1;
+          /* J %= HILBERT_DIM; */
+       }
+    }
+  for (j = 0; j < HILBERT_DIM; j++)
+    dest[j] &= ~(P(t))0 >> (8 * sizeof(P(t)) - WORDBITS);
+}
+#endif
+
+#ifdef HILBERT_WANT_DECODE
+/*
+ * given the sequence number of a point, it finds the coordinates of the point
+ * on the Hilbert Curve
+ */
+static void
+P(decode) (P(t) *dest, P(t) *src)
+{
+  P(t) mask = (P(t))1 << WORDBITS - 1, element, temp1, temp2,
+    A, W = 0, S, tS, T, tT, J, P = 0, xJ;
+  uns i = NUMBITS * HILBERT_DIM - HILBERT_DIM, j;
+
+  for (j = 0; j < HILBERT_DIM; j++)
+    dest[j] = 0;
+
+  /*--- P ---*/
+  element = i / WORDBITS;
+  P = src[element];
+  if (i % WORDBITS > WORDBITS - HILBERT_DIM)
+    {
+      temp1 = src[element + 1];
+      P >>= i % WORDBITS;
+      temp1 <<= WORDBITS - i % WORDBITS;
+      P |= temp1;
+    }
+  else
+    P >>= i % WORDBITS;        /* P is a HILBERT_DIM bit hcode */
+
+  /* the & masks out spurious highbit values */
+  if (HILBERT_DIM < WORDBITS)
+    P &= (1 << HILBERT_DIM) -1;
+
+  /*--- xJ ---*/
+  J = HILBERT_DIM;
+  for (j = 1; j < HILBERT_DIM; j++)
+    if ((P >> j & 1) == (P & 1))
+      continue;
+    else
+      break;
+  if (j != HILBERT_DIM)
+    J -= j;
+  xJ = J - 1;
+
+  /*--- S, tS, A ---*/
+  A = S = tS = P ^ P / 2;
+
+
+  /*--- T ---*/
+  if (P < 3)
+    T = 0;
+  else
+    if (P % 2)
+      T = (P - 1) ^ (P - 1) / 2;
+    else
+      T = (P - 2) ^ (P - 2) / 2;
+
+  /*--- tT ---*/
+  tT = T;
+
+  /*--- distrib bits to coords ---*/
+  for (j = HILBERT_DIM - 1; P > 0; P >>=1, j--)
+    if (P & 1)
+      dest[j] |= mask;
+
+
+  for (i -= HILBERT_DIM, mask >>= 1; (int)i >= 0; i -= HILBERT_DIM, mask >>= 1)
+    {
+      /*--- P ---*/
+      element = i / WORDBITS;
+      P = src[element];
+      if (i % WORDBITS > WORDBITS - HILBERT_DIM)
+       {
+         temp1 = src[element + 1];
+         P >>= i % WORDBITS;
+         temp1 <<= WORDBITS - i % WORDBITS;
+         P |= temp1;
+       }
+      else
+       P >>= i % WORDBITS;     /* P is a HILBERT_DIM bit hcode */
+
+      /* the & masks out spurious highbit values */
+      if (HILBERT_DIM < WORDBITS)
+        P &= (1 << HILBERT_DIM) -1;
+
+      /*--- S ---*/
+      S = P ^ P / 2;
+
+      /*--- tS ---*/
+      if (xJ % HILBERT_DIM != 0)
+       {
+         temp1 = S >> xJ % HILBERT_DIM;
+         temp2 = S << HILBERT_DIM - xJ % HILBERT_DIM;
+         tS = temp1 | temp2;
+         tS &= ((P(t))1 << HILBERT_DIM) - 1;
+       }
+      else
+       tS = S;
+
+      /*--- W ---*/
+      W ^= tT;
+
+      /*--- A ---*/
+      A = W ^ tS;
+
+      /*--- distrib bits to coords ---*/
+      for (j = HILBERT_DIM - 1; A > 0; A >>=1, j--)
+       if (A & 1)
+         dest[j] |= mask;
+
+      if (i > 0)
+       {
+         /*--- T ---*/
+         if (P < 3)
+           T = 0;
+         else
+           if (P % 2)
+             T = (P - 1) ^ (P - 1) / 2;
+           else
+             T = (P - 2) ^ (P - 2) / 2;
+
+         /*--- tT ---*/
+         if (xJ % HILBERT_DIM != 0)
+           {
+             temp1 = T >> xJ % HILBERT_DIM;
+             temp2 = T << HILBERT_DIM - xJ % HILBERT_DIM;
+             tT = temp1 | temp2;
+             tT &= ((P(t))1 << HILBERT_DIM) - 1;
+           }
+         else
+           tT = T;
+
+         /*--- xJ ---*/
+         J = HILBERT_DIM;
+         for (j = 1; j < HILBERT_DIM; j++)
+           if ((P >> j & 1) == (P & 1))
+             continue;
+           else
+             break;
+         if (j != HILBERT_DIM)
+           J -= j;
+         xJ += J - 1;
+       }
+    }
+}
+#endif
+
+#undef P
+#undef HILBERT_PREFIX
+#undef HILBERT_DIM
+#undef HILBERT_TYPE
+#undef HILBERT_ORDER
+#undef HILBERT_WANT_DECODE
+#undef HILBERT_WANT_ENCODE
+#undef NUMBITS
+#undef WORDBITS
diff --git a/images/image-idx.c b/images/image-idx.c
new file mode 100644 (file)
index 0000000..9e5ece0
--- /dev/null
@@ -0,0 +1,944 @@
+// FIXME: this file is full of experiments... will be completely different in final version
+
+#undef LOCAL_DEBUG
+
+#include "sherlock/sherlock.h"
+#include "lib/mempool.h"
+#include "lib/conf.h"
+#include "lib/getopt.h"
+#include "lib/fastbuf.h"
+#include "lib/chartype.h"
+#include "sherlock/object.h"
+#include "lib/url.h"
+#include "lib/unicode.h"
+#include "sherlock/lizard-fb.h"
+#include "sherlock/tagged-text.h"
+#include "charset/charconv.h"
+#include "charset/unicat.h"
+#include "charset/fb-charconv.h"
+#include "indexer/indexer.h"
+#include "indexer/lexicon.h"
+#include "indexer/params.h"
+#include "utils/dumpconfig.h"
+#include "lang/lang.h"
+#include "lib/base224.h"
+#include "lib/bbuf.h"
+#include "lib/clists.h"
+
+#include "images/images.h"
+#include "images/image-obj.h"
+#include "images/image-sig.h"
+#include "images/dup-cmp.h"
+#include "images/kd-tree.h"
+#include "images/color.h"
+
+#include <stdlib.h>
+#include <fcntl.h>
+#include <string.h>
+
+static struct fastbuf *fb_cards;
+static struct fastbuf *fb_card_attrs;
+static struct buck2obj_buf *buck2obj;
+
+/* This should happen in gatherer or scanner */
+static void
+generate_signatures(uns limit)
+{
+  fb_cards = index_bopen("cards", O_RDONLY);
+  fb_card_attrs = index_bopen("card-attrs", O_RDONLY);
+  struct fastbuf *fb_signatures = index_bopen("image-sig", O_CREAT | O_WRONLY | O_TRUNC);
+  struct card_attr ca;
+  struct image_signature sig;
+  struct mempool *pool = mp_new(1 << 16);
+  struct buck2obj_buf *bob = buck2obj_alloc();
+  uns count = 0;
+
+  if (limit == ~0U)
+    log(L_INFO, "Generating image signatures");
+  else
+    log(L_INFO, "Generating at most %d image signatures", limit);
+  bputl(fb_signatures, 0);
+  imo_decompress_thumbnails_init();
+
+  for (oid_t oid = 0; bread(fb_card_attrs, &ca, sizeof(ca)); oid++)
+    if ((uns)((ca.type_flags >> 4) - 8) < 4)
+      {
+        bsetpos(fb_cards, (sh_off_t)ca.card << CARD_POS_SHIFT);
+        uns buck_len = bgetl(fb_cards) - (LIZARD_COMPRESS_HEADER - 1);
+        uns buck_type = bgetc(fb_cards) + BUCKET_TYPE_PLAIN;
+        mp_flush(pool);
+        struct odes *obj = obj_read_bucket(bob, pool, buck_type, buck_len, fb_cards, NULL);
+        struct oattr *attr;
+        if (!obj)
+          die("Failed to read card");
+        if (attr = obj_find_attr(obj, 'N'))
+          {
+#ifdef LOCAL_DEBUG
+           byte *url = obj_find_aval(obj_find_attr(obj, 'U' + OBJ_ATTR_SON)->son, 'U');
+           DBG("Reading oid=%d url=%s", oid, url);
+#endif
+           struct image_obj imo;
+           imo_init(&imo, pool, obj);
+           if (imo_decompress_thumbnail(&imo))
+             {
+               if (compute_image_signature(&imo.thumb, &sig))
+                 {
+                   bwrite(fb_signatures, &oid, sizeof(oid));
+                   bwrite(fb_signatures, &sig.vec, sizeof(struct image_vector));
+                   bputc(fb_signatures, sig.len);
+                   if (sig.len)
+                     bwrite(fb_signatures, sig.reg, sig.len * sizeof(struct image_region));
+                   count++;
+                   if (count % 10000 == 0)
+                     log(L_DEBUG, "... passed %d images", count);
+                   if (count >= limit)
+                     break;
+                 }
+               else
+                 DBG("Cannot create signature");
+             }
+           else
+             DBG("Cannot decompress thumbnail");
+         }
+      }
+  brewind(fb_signatures);
+  bputl(fb_signatures, count);
+  DBG("%d signatures written", count);
+
+  imo_decompress_thumbnails_done();
+  buck2obj_free(bob);
+  mp_delete(pool);
+  bclose(fb_cards);
+  bclose(fb_card_attrs);
+  bclose(fb_signatures);
+}
+
+/*********************************************************************************/
+
+struct vectors_node {
+  oid_t oid;
+  u32 temp;
+  struct image_vector vec;
+};
+
+static uns vectors_count;
+static struct vectors_node *vectors;
+
+static void
+vectors_read(void)
+{
+  log(L_DEBUG, "Reading signature vectors");
+  struct fastbuf *fb = index_bopen("image-sig", O_RDONLY);
+  vectors_count = bgetl(fb);
+  if (vectors_count)
+    {
+      vectors = xmalloc(vectors_count * sizeof(struct vectors_node));
+      for (uns i = 0; i < vectors_count; i++)
+        {
+         bread(fb, &vectors[i].oid, sizeof(oid_t));
+         bread(fb, &vectors[i].vec, sizeof(struct image_vector));
+         bskip(fb, bgetc(fb) * sizeof(struct image_region));
+       }
+    }
+  bclose(fb);
+}
+
+static void
+vectors_cleanup(void)
+{
+  log(L_DEBUG, "Freeing signature vectors");
+  if (vectors_count)
+    xfree(vectors);
+}
+
+/*********************************************************************************/
+
+static u64 random_clusters_max_size = 500000;
+static uns random_clusters_max_count = 1000;
+
+#define RANDOM_CLUSTERS_SIZE   0x7fffffff
+#define RANDOM_CLUSTERS_LAST   0x80000000
+
+static struct random_clusters_node {
+  struct vectors_node *node;
+  s32 dot_prod;
+} *random_clusters_temp;
+static uns random_clusters_count;
+
+#define ASORT_PREFIX(x) random_clusters_##x
+#define ASORT_KEY_TYPE s32
+#define ASORT_ELT(i) start[i].dot_prod
+#define ASORT_SWAP(i,j) do { struct random_clusters_node _s = start[i]; start[i] = start[j]; start[j] = _s; } while(0)
+#define ASORT_EXTRA_ARGS , struct random_clusters_node *start
+#include "lib/arraysort.h"
+
+static void
+random_clusters_init(void)
+{
+  if (!vectors_count)
+    return;
+  log(L_INFO, "Initializing random clusters generator");
+  random_clusters_temp = xmalloc(vectors_count * sizeof(struct random_clusters_node));
+  for (uns i = 0; i < vectors_count; i++)
+    random_clusters_temp[i].node = vectors + i;
+}
+
+static void
+random_clusters_build(void)
+{
+  random_clusters_count = 0;
+  if (!vectors_count)
+    return;
+
+  log(L_INFO, "Generating random clusters for duplicates comparision");
+
+  for (uns i = 0; i < vectors_count; i++)
+    vectors[i].temp &= RANDOM_CLUSTERS_SIZE;
+
+  /* Initialize recursion */
+  struct stk {
+    uns count;
+    struct random_clusters_node *start;
+  } stk_top[64], *stk = stk_top + 1;
+  stk->start = random_clusters_temp;
+  stk->count = vectors_count;
+
+  /* Main loop */
+  while (stk != stk_top)
+    {
+      /* Split conditions */
+      uns split;
+      if (stk->count < 2)
+       split = 0;
+      else if (stk->count > random_clusters_max_count)
+       split = 1;
+      else
+        {
+          s64 size = random_clusters_max_size;
+          for (uns i = 0; i < stk->count && size >= 0; i++)
+           size -= stk->start[i].node->temp;
+         split = size < 0;
+       }
+
+      /* BSP leaf node */
+      if (!split)
+        {
+         stk->start[stk->count - 1].node->temp |= RANDOM_CLUSTERS_LAST;
+         random_clusters_count++;
+         stk--;
+       }
+
+      /* BSP internal node */
+      else
+        {
+         /* Generate random normal vector of the splitting plane */
+         int normal[IMAGE_VEC_K];
+         for (uns i = 0; i < IMAGE_VEC_K; i++)
+           normal[i] = random_max(0x20001) - 0x10000;
+
+         /* Compute dot produts */
+         for (uns i = 0; i < stk->count; i++)
+           {
+             stk->start[i].dot_prod = 0;
+             for (uns j = 0; j < IMAGE_VEC_K; j++)
+               stk->start[i].dot_prod += normal[j] * stk->start[i].node->vec.f[j];
+           }
+
+         /* Sort... could be faster, because we only need the median */
+         random_clusters_sort(stk->count, stk->start);
+
+         /* Split in the middle */
+         stk[1].count = stk[0].count >> 1;
+         stk[0].count -= stk[1].count;
+         stk[1].start = stk[0].start;
+         stk[0].start += stk[1].count;
+         stk++;
+       }
+    }
+  log(L_INFO, "Generated %u clusters", random_clusters_count);
+}
+
+static void
+random_clusters_cleanup(void)
+{
+  if (vectors_count)
+    xfree(random_clusters_temp);
+}
+
+/*********************************************************************************/
+
+// FIXME: use vectors_read()... duplicate code
+
+struct signature_record {
+  oid_t oid;
+  struct image_vector vec;
+};
+
+#define ASORT_PREFIX(x) build_search_tree_##x
+#define ASORT_KEY_TYPE struct signature_record *
+#define ASORT_ELT(i) rec[i]
+#define ASORT_LT(x,y) x->vec.f[dim] < y->vec.f[dim]
+#define ASORT_EXTRA_ARGS , uns dim, struct signature_record **rec
+#include "lib/arraysort.h"
+
+#if 0
+#define DBG_KD(x...) DBG(x)
+#else
+#define DBG_KD(x...) do{}while(0)
+#endif
+
+static struct image_tree tree;
+static struct signature_record *records;
+static struct signature_record **precords;
+
+static void
+build_kd_tree(void)
+{
+  log(L_INFO, "Building KD-tree");
+
+  struct fastbuf *fb_signatures = index_bopen("image-sig", O_RDONLY);
+  tree.count = bgetl(fb_signatures);
+  ASSERT(tree.count < 0x80000000);
+  if (!tree.count)
+    {
+      /* FIXME */
+      bclose(fb_signatures);
+      die("There are no signatures");
+    }
+  else
+    {
+      DBG("Reading %d signatures", tree.count);
+      records = xmalloc(tree.count * sizeof(struct signature_record));
+      precords = xmalloc(tree.count * sizeof(void *));
+      for (uns i = 0; i < tree.count; i++)
+        {
+         bread(fb_signatures, &records[i].oid, sizeof(oid_t));
+         bread(fb_signatures, &records[i].vec, sizeof(struct image_vector));
+         uns len = bgetc(fb_signatures);
+         bskip(fb_signatures, len * sizeof(struct image_region));
+         precords[i] = records + i;
+         if (likely(i))
+           for (uns j = 0; j < IMAGE_VEC_K; j++)
+             {
+               tree.bbox.vec[0].f[j] = MIN(tree.bbox.vec[0].f[j], records[i].vec.f[j]);
+               tree.bbox.vec[1].f[j] = MAX(tree.bbox.vec[1].f[j], records[i].vec.f[j]);
+             }
+         else
+            tree.bbox.vec[0] = tree.bbox.vec[1] = records[0].vec;
+       }
+      bclose(fb_signatures);
+
+      for (tree.depth = 1; (uns)(2 << tree.depth) < tree.count; tree.depth++);
+      DBG("depth=%d nodes=%d bbox=[(%s), (%s)]", tree.depth, 1 << tree.depth,
+         stk_print_image_vector(tree.bbox.vec + 0), stk_print_image_vector(tree.bbox.vec + 1));
+      uns leaves_index = 1 << (tree.depth - 1);
+      tree.nodes = xmalloc_zero((1 << tree.depth) * sizeof(struct image_node));
+      tree.leaves = xmalloc_zero(tree.count * sizeof(struct image_leaf));
+
+      /* Initialize recursion */
+      struct stk {
+       struct image_bbox bbox;
+       uns index, count;
+       struct signature_record **start;
+      } stk_top[32], *stk = stk_top + 1;
+      stk->index = 1;
+      stk->start = precords;
+      stk->count = tree.count;
+      stk->bbox.vec[0] = tree.bbox.vec[0];
+      for (uns i = 0; i < IMAGE_VEC_K; i++)
+       stk->bbox.vec[1].f[i] = tree.bbox.vec[1].f[i] - tree.bbox.vec[0].f[i];
+      uns entry_index = 0;
+
+      /* Main loop */
+      while (stk != stk_top)
+        {
+         DBG_KD("Main loop... depth=%d index=%d count=%d, start=%d, min=%s dif=%s",
+             stk - stk_top, stk->index, stk->count, stk->start - precords,
+             stk_print_image_vector(stk->bbox.vec + 0), stk_print_image_vector(stk->bbox.vec + 1));
+         ASSERT(stk->count);
+
+         /* Create leaf node */
+         if (stk->index >= leaves_index || stk->count < 2)
+           {
+             tree.nodes[stk->index].val = IMAGE_NODE_LEAF | entry_index;
+             for (; stk->count--; stk->start++)
+               {
+                 struct image_leaf *leaf = &tree.leaves[entry_index++];
+                 struct signature_record *record = *stk->start;
+                 leaf->oid = record->oid;
+                 leaf->flags = 0;
+                 for (uns i = IMAGE_VEC_K; i--; )
+                   {
+                     uns bits = IMAGE_LEAF_BITS(i);
+                     leaf->flags <<= bits;
+                     if (stk->bbox.vec[1].f[i])
+                       {
+                         uns value =
+                           (record->vec.f[i] - stk->bbox.vec[0].f[i]) *
+                           ((1 << bits) - 1) / stk->bbox.vec[1].f[i];
+                         ASSERT(value < (uns)(1 << bits));
+                         leaf->flags |= value;
+                       }
+                   }
+                 if (!stk->count)
+                   leaf->flags |= IMAGE_LEAF_LAST;
+                 DBG_KD("Creating leaf node; oid=%d vec=(%s) flags=0x%08x",
+                     leaf->oid, stk_print_image_vector(&record->vec), leaf->flags);
+               }
+             stk--;
+           }
+
+         /* Create internal node */
+         else
+           {
+             /* Select dimension to splis */
+             uns dim = 0;
+             for (uns i = 1; i < IMAGE_VEC_K; i++)
+               if (stk->bbox.vec[1].f[i] > stk->bbox.vec[1].f[dim])
+                 dim = i;
+
+             /* Sort... FIXME: we only need the median */
+             build_search_tree_sort(stk->count, dim, stk->start);
+
+             /* Split in the middle */
+             uns index = stk->index;
+             stk[1].index = stk[0].index * 2;
+             stk[0].index = stk[1].index + 1;
+             stk[1].count = stk[0].count >> 1;
+             stk[0].count -= stk[1].count;
+             stk[1].start = stk[0].start;
+             stk[0].start += stk[1].count;
+
+             /* Choose split value */
+             uns lval = stk->start[-1]->vec.f[dim];
+             uns rval = stk->start[0]->vec.f[dim];
+             uns pivot = stk->bbox.vec[0].f[dim] + (stk->bbox.vec[1].f[dim] >> 1);
+             if (pivot <= lval)
+               pivot = lval;
+             else if (pivot >= rval)
+               pivot = rval;
+
+             DBG_KD("Created internal node; dim=%d pivot=%d", dim, pivot);
+
+             /* Split the box */
+             stk[1].bbox = stk[0].bbox;
+              stk[1].bbox.vec[1].f[dim] = pivot - stk[0].bbox.vec[0].f[dim];
+             stk[0].bbox.vec[0].f[dim] += stk[1].bbox.vec[1].f[dim];
+             stk[0].bbox.vec[1].f[dim] -= stk[1].bbox.vec[1].f[dim];
+
+             /* Fill the node structure */
+             tree.nodes[index].val = dim + (pivot << 8);
+             stk++;
+           }
+       }
+
+      DBG("Tree constructed, saving...");
+
+      struct fastbuf *fb_tree = index_bopen("image-tree", O_CREAT | O_WRONLY | O_TRUNC);
+      bputl(fb_tree, tree.count);
+      bputl(fb_tree, tree.depth);
+      bwrite(fb_tree, &tree.bbox, sizeof(struct image_bbox));
+      bwrite(fb_tree, tree.nodes + 1, ((1 << tree.depth) - 1) * sizeof(struct image_node));
+      bwrite(fb_tree, tree.leaves, tree.count * sizeof(struct image_leaf));
+      bclose(fb_tree);
+
+      //xfree(tree.leaves);
+      //xfree(tree.nodes);
+      //xfree(precords);
+      //xfree(records);
+    }
+}
+
+/*********************************************************************************/
+
+#if 0
+
+struct pass1_hilbert {
+  u32 index;
+  struct image_vector vec;
+};
+
+struct pass1_node {
+  cnode lru_node;
+  cnode buf_node;
+  uns buf_size;
+  byte *buf;
+  oid_t oid;
+  byte *url;
+  struct image image;
+  struct image_dup dup;
+};
+
+static uns pass1_buf_size = 400 << 20;
+static uns pass1_max_count = 100000;
+static uns pass1_search_dist = 40;
+static uns pass1_search_count = 500;
+
+static struct mempool *pass1_pool;
+static struct pass1_hilbert *pass1_hilbert_list;
+static byte *pass1_buf_start;
+static byte *pass1_buf_pos;
+static uns pass1_buf_free;
+static uns pass1_buf_used;
+static clist pass1_buf_list;
+static clist pass1_lru_list;
+static u64 pass1_lookups;
+static u64 pass1_reads;
+static u64 pass1_pairs;
+static u64 pass1_dups;
+static u64 pass1_shrinks;
+static u64 pass1_alloc_sum;
+
+#define HILBERT_PREFIX(x) pass1_hilbert_##x
+#define HILBERT_TYPE byte
+#define HILBERT_ORDER 8
+#define HILBERT_DIM IMAGE_VEC_K
+#define HILBERT_WANT_ENCODE
+#include "images/hilbert.h"
+
+#define ASORT_PREFIX(x) pass1_hilbert_sort_##x
+#define ASORT_KEY_TYPE struct image_vector *
+#define ASORT_ELT(i) (&pass1_hilbert_list[i].vec)
+#define ASORT_LT(x,y) (memcmp(x, y, sizeof(*x)) < 0)
+#define ASORT_SWAP(i,j) do { struct pass1_hilbert _s;          \
+               _s = pass1_hilbert_list[i];                     \
+               pass1_hilbert_list[i] = pass1_hilbert_list[j];  \
+               pass1_hilbert_list[j] = _s; } while(0)
+#include "lib/arraysort.h"
+
+static void
+pass1_hilbert_sort(void)
+{
+  DBG("Computing positions on the Hilbert curve");
+  pass1_hilbert_list = xmalloc(tree.count * sizeof(struct pass1_hilbert));
+  for (uns i = 0; i < tree.count; i++)
+    {
+      struct pass1_hilbert *h = pass1_hilbert_list + i;
+      h->index = i;
+      byte vec[IMAGE_VEC_K];
+      pass1_hilbert_encode(vec, precords[i]->vec.f);
+      for (uns j = 0; j < IMAGE_VEC_K; j++)
+       h->vec.f[j] = vec[IMAGE_VEC_K - 1 - j];
+    }
+  DBG("Sorting signatures in order of incresing parameters on the Hilbert curve");
+  pass1_hilbert_sort_sort(tree.count);
+#if 0
+  for (uns i = 0; i < tree.count; i++)
+    {
+      if (i)
+        {
+         byte *v1 = precords[pass1_hilbert_list[i - 1].index]->vec.f;
+         byte *v2 = precords[pass1_hilbert_list[i].index]->vec.f;
+#define SQR(x) ((x)*(x))
+          uns dist = 0;
+         for (uns j = 0; j < 6; j++)
+           dist += SQR(v1[j] - v2[j]);
+         DBG("dist %d", dist);
+       }
+      DBG("index %d", pass1_hilbert_list[i].index);
+    }
+#endif
+}
+
+static void
+pass1_hilbert_cleanup(void)
+{
+  xfree(pass1_hilbert_list);
+}
+
+#define HASH_PREFIX(x) pass1_hash_##x
+#define HASH_NODE struct pass1_node
+#define HASH_KEY_ATOMIC oid
+#define HASH_WANT_CLEANUP
+#define HASH_WANT_FIND
+#define HASH_WANT_NEW
+#define HASH_WANT_REMOVE
+#include "lib/hashtable.h"
+
+static inline void
+pass1_buf_init(void)
+{
+  //DBG("pass1_buf_init()");
+  pass1_buf_free = pass1_buf_size;
+  pass1_buf_start = pass1_buf_pos = xmalloc(pass1_buf_size);
+  pass1_buf_used = 0;
+}
+
+static inline void
+pass1_buf_cleanup(void)
+{
+  //DBG("pass1_buf_cleanup()");
+  xfree(pass1_buf_start);
+}
+
+static void
+pass1_node_free(struct pass1_node *node)
+{
+  //DBG("pass1_node_free(%d)", (uns)node->oid);
+  if (node->buf_size)
+    {
+      pass1_buf_used -= node->buf_size;
+      clist_remove(&node->buf_node);
+    }
+  clist_remove(&node->lru_node);
+  pass1_hash_remove(node);
+}
+
+static inline void
+pass1_node_free_lru(void)
+{
+  ASSERT(!clist_empty(&pass1_lru_list));
+  pass1_node_free(SKIP_BACK(struct pass1_node, lru_node, clist_head(&pass1_lru_list)));
+}
+
+static inline void
+pass1_node_after_move(struct pass1_node *node, addr_int_t move)
+{
+  //DBG("pass1_node_after_mode(%d, %d)", (uns)node->oid, (uns)move);
+  /* adjust internal pointers */
+#define MOVE(x) x = (byte *)(x) - move
+  MOVE(node->url);
+  MOVE(node->image.pixels);
+  MOVE(node->dup.buf);
+#undef MOVE
+}
+
+static inline void
+pass1_buf_shrink(void)
+{
+  DBG("pass1_buf_shrink()");
+  pass1_shrinks++;
+  pass1_buf_free = pass1_buf_size;
+  pass1_buf_pos = pass1_buf_start;
+  CLIST_FOR_EACH(void *, p, pass1_buf_list)
+    {
+      struct pass1_node *node = SKIP_BACK(struct pass1_node, buf_node, p);
+      if (node->buf != pass1_buf_pos)
+        {
+          memmove(pass1_buf_pos, node->buf, node->buf_size);
+          pass1_node_after_move(node, node->buf - pass1_buf_pos);
+          node->buf = pass1_buf_pos;
+       }
+      pass1_buf_pos += node->buf_size;
+      pass1_buf_free -= node->buf_size;
+    }
+}
+
+static void *
+pass1_buf_alloc(uns size)
+{
+  //DBG("pass1_buf_alloc(%d)", size);
+
+  /* if there is not enough free space at the end of the buffer */
+  if (size > pass1_buf_free)
+    {
+      /* free some lru nodes */
+      //DBG("freeing lru nodes");
+      while (size > pass1_buf_size - pass1_buf_used || pass1_buf_used > pass1_buf_size / 2)
+        {
+         if (unlikely(clist_empty(&pass1_lru_list))) // FIXME
+           die("Buffer too small");
+          pass1_node_free_lru();
+       }
+
+      pass1_buf_shrink();
+    }
+
+  /* final allocation */
+  void *result = pass1_buf_pos;
+  pass1_buf_pos += size;
+  pass1_buf_free -= size;
+  pass1_buf_used += size;
+  pass1_alloc_sum += size;
+  return result;
+}
+
+static struct pass1_node *
+pass1_node_new(oid_t oid)
+{
+  DBG("pass1_node_new(%d)", (uns)oid);
+  if (pass1_hash_table.hash_count == pass1_max_count)
+    pass1_node_free_lru();
+  struct pass1_node *node = pass1_hash_new(oid);
+  mp_flush(pass1_pool);
+  pass1_reads++;
+
+  /* read object */
+  struct card_attr ca;
+  bsetpos(fb_card_attrs, (sh_off_t)oid * sizeof(ca)); /* FIXME: these seeks can be easily removed */
+  bread(fb_card_attrs, &ca, sizeof(ca));
+
+  bsetpos(fb_cards, (sh_off_t)ca.card << CARD_POS_SHIFT); /* FIXME: maybe a presort should handle these random seeks */
+  uns buck_len = bgetl(fb_cards) - (LIZARD_COMPRESS_HEADER - 1);
+  uns buck_type = bgetc(fb_cards) + BUCKET_TYPE_PLAIN;
+  struct odes *obj = obj_read_bucket(buck2obj, pass1_pool, buck_type, buck_len, fb_cards, NULL);
+  if (unlikely(!obj))
+    die("Failed to read card");
+  byte *url = obj_find_aval(obj_find_attr(obj, 'U' + OBJ_ATTR_SON)->son, 'U');
+  uns url_len = strlen(url);
+
+  /* decompress thumbnail */
+  struct image_obj imo;
+  imo_init(&imo, pass1_pool, obj);
+  if (unlikely(!imo_decompress_thumbnail(&imo)))
+    die("Cannot decompress thumbnail");
+  node->image = imo.thumb;
+
+  /* create duplicates comparision object */
+  image_dup_init(&node->dup, &node->image, pass1_pool);
+
+  /* copy data */
+  //DBG("loaded image %s s=%d d=%d", url, node->image.size, node->dup.buf_size);
+  node->buf_size = node->image.size + node->dup.buf_size + url_len + 1;
+  if (node->buf_size)
+    {
+      byte *buf = node->buf = pass1_buf_alloc(node->buf_size);
+      clist_add_tail(&pass1_buf_list, &node->buf_node);
+#define COPY(ptr, size) ({ void *_p=buf; uns _size=(size); buf+=_size; memcpy(_p,(ptr),_size); _p; })
+      node->url = COPY(url, url_len + 1);
+      node->image.pixels = COPY(node->image.pixels, node->image.size);
+      node->dup.buf = COPY(node->dup.buf, node->dup.buf_size);
+#undef COPY
+    }
+
+  /* add to lru list */
+  return node;
+}
+
+static inline struct pass1_node *
+pass1_node_lock(oid_t oid)
+{
+  DBG("pass1_node_lock(%d)", (uns)oid);
+  pass1_lookups++;
+  struct pass1_node *node = pass1_hash_find(oid);
+  if (node)
+    {
+      clist_remove(&node->lru_node);
+      return node;
+    }
+  else
+    return pass1_node_new(oid);
+}
+
+static inline void
+pass1_node_unlock(struct pass1_node *node)
+{
+  //DBG("pass1_node_unlock(%d)", (uns)node->oid);
+  clist_add_tail(&pass1_lru_list, &node->lru_node);
+}
+
+static void
+pass1_show_stats(void)
+{
+  log(L_INFO, "%d count, %Ld lookups, %Ld reads, %Ld pairs, %Ld dups, %Ld shrinks", tree.count,
+    (long long int)pass1_lookups, (long long int)pass1_reads,
+    (long long int)pass1_pairs, (long long int)pass1_dups, (long long int)pass1_shrinks);
+}
+
+static void
+pass1(void)
+{
+  log(L_INFO, "Looking for duplicates");
+  ASSERT(tree.nodes);
+
+  /* initialization */
+  pass1_lookups = pass1_reads = pass1_pairs = pass1_dups = pass1_shrinks = pass1_alloc_sum = 0;
+  fb_cards = bopen("index/cards", O_RDONLY, 10000); // FIXME
+  fb_card_attrs = bopen("index/card-attrs", O_RDONLY, sizeof(struct card_attr)); // FIXME
+  buck2obj = buck2obj_alloc();
+  imo_decompress_thumbnails_init();
+  clist_init(&pass1_lru_list);
+  clist_init(&pass1_buf_list);
+  pass1_hash_init();
+  pass1_buf_init();
+  pass1_pool = mp_new(1 << 20);
+
+  /* Hilbert sort */
+  pass1_hilbert_sort();
+  pass1_hilbert_cleanup();
+
+  /* main loop */
+  for (uns i = 0; i < tree.count; )
+    {
+      /* lookup next image */
+      oid_t oid = tree.leaves[i].oid;
+      struct pass1_node *node = pass1_node_lock(oid);
+
+      /* compare with all near images */
+      struct image_search search;
+      image_search_init(&search, &tree, &precords[i]->vec, pass1_search_dist);
+      /* FIXME: can be faster than general search in KD-tree */
+      oid_t oid2;
+      uns dist;
+      for (uns j = 0; j < pass1_search_count && image_search_next(&search, &oid2, &dist); j++)
+        {
+         if (oid < oid2)
+            {
+             struct pass1_node *node2 = pass1_node_lock(oid2);
+             DBG("comparing %d and %d", oid, oid2);
+             if (image_dup_compare(&node->dup, &node2->dup, IMAGE_DUP_TRANS_ID))
+               {
+                 pass1_dups++;
+                 log(L_DEBUG, "*** Found duplicates oid1=0x%x oid=0x%x", (uns)node->oid, (uns)node2->oid);
+                 log(L_DEBUG, "  %s", node->url);
+                 log(L_DEBUG, "  %s", node2->url);
+               }
+             pass1_pairs++;
+             pass1_node_unlock(node2);
+           }
+       }
+      image_search_done(&search);
+      pass1_node_unlock(node);
+      i++;
+      if (i % 1000 == 0)
+        log(L_DEBUG, "... passed %d images", i);
+    }
+
+  /* clean up */
+  pass1_hash_cleanup();
+  pass1_buf_cleanup();
+  mp_delete(pass1_pool);
+  bclose(fb_cards);
+  bclose(fb_card_attrs);
+  buck2obj_free(buck2obj);
+  imo_decompress_thumbnails_done();
+
+  /* print statistics */
+  pass1_show_stats();
+}
+
+/*********************************************************************************/
+
+static uns pass2_clusterings_count = 1;
+
+static void
+pass2_estimate_sizes(void)
+{
+  if (!vectors_count)
+    return;
+  log(L_DEBUG, "Reading image sizes");
+
+  /* FIXME: hack, these reads are not necessary, can be done in previous phases */
+  struct fastbuf *fb_cards = index_bopen("cards", O_RDONLY);
+  struct fastbuf *fb_card_attrs = index_bopen("card-attrs", O_RDONLY);
+  struct mempool *pool = mp_new(1 << 16);
+  struct buck2obj_buf *bob = buck2obj_alloc();
+
+  for (uns i = 0; i < vectors_count; i++)
+    {
+      oid_t oid = vectors[i].oid;
+      struct card_attr ca;
+      bsetpos(fb_card_attrs, (sh_off_t)oid * sizeof(ca));
+      bread(fb_card_attrs, &ca, sizeof(ca));
+      bsetpos(fb_cards, (sh_off_t)ca.card << CARD_POS_SHIFT);
+      uns buck_len = bgetl(fb_cards) - (LIZARD_COMPRESS_HEADER - 1);
+      uns buck_type = bgetc(fb_cards) + BUCKET_TYPE_PLAIN;
+      mp_flush(pool);
+      struct odes *obj = obj_read_bucket(bob, pool, buck_type, buck_len, fb_cards, NULL);
+      byte *attr = obj_find_aval(obj, 'G');
+      ASSERT(attr);
+      uns image_width, image_height, image_colors, thumb_width, thumb_height;
+      byte color_space[MAX_ATTR_SIZE];
+      sscanf(attr, "%d%d%s%d%d%d", &image_width, &image_height, color_space, &image_colors, &thumb_width, &thumb_height);
+      vectors[i].temp = image_dup_estimate_size(thumb_width, thumb_height) +
+       sizeof(struct image) + thumb_width * thumb_height * 3;
+    }
+  buck2obj_free(bob);
+  mp_delete(pool);
+  bclose(fb_cards);
+  bclose(fb_card_attrs);
+}
+
+static void
+pass2(void)
+{
+  // FIXME: presorts, much allocated memory when not needed
+  vectors_read();
+  pass2_estimate_sizes();
+  random_clusters_init();
+  for (uns clustering = 0; clustering < pass2_clusterings_count; clustering++)
+    {
+      random_clusters_build();
+      // FIXME
+      // - external sort
+      // - generate and compare pairs in clusters
+    }
+  random_clusters_cleanup();
+  vectors_cleanup();
+}
+#endif
+
+/*********************************************************************************/
+
+static char *shortopts = CF_SHORT_OPTS "";
+static struct option longopts[] =
+{
+  CF_LONG_OPTS
+  { NULL, 0, 0, 0 }
+};
+
+static char *help = "\
+Usage: image-indexer [<options>]\n\
+\n\
+Options:\n" CF_USAGE;
+
+static void NONRET
+usage(byte *msg)
+{
+  if (msg)
+  {
+    fputs(msg, stderr);
+    fputc('\n', stderr);
+  }
+  fputs(help, stderr);
+  exit(1);
+}
+
+
+int
+main(int argc UNUSED, char **argv)
+{
+  int opt;
+
+  log_init(argv[0]);
+  while ((opt = cf_getopt(argc, argv, shortopts, longopts, NULL)) >= 0)
+    switch (opt)
+    {
+      default:
+      usage("Invalid option");
+    }
+  if (optind != argc)
+    usage("Invalid usage");
+
+  srgb_to_luv_init();
+
+#if 0
+  while (1)
+  {
+  struct mempool *pool = mp_new(1024);
+  struct fastbuf *fb = bopen("a.jpg", O_RDONLY, 1024);
+  struct image_io io;
+  log(L_DEBUG, "opening");
+  image_open(&io, fb, pool);
+  io.format = IMAGE_FORMAT_JPEG;
+  log(L_DEBUG, "reading");
+  int i;
+  i = image_read(&io);
+  log(L_DEBUG, "done %d %d %d", i, io.image.width, io.image.height);
+  for (i = 0; i < 1000000000; i++)
+    ;
+  image_close(&io);
+  mp_delete(pool);
+  bclose(fb);
+  }
+#endif  
+
+#if 0  
+  generate_signatures(20000);
+  build_kd_tree();
+  //pass1();
+  pass2();
+#endif  
+
+  return 0;
+}
diff --git a/images/image-obj.c b/images/image-obj.c
new file mode 100644 (file)
index 0000000..5823048
--- /dev/null
@@ -0,0 +1,518 @@
+/*
+ *     Image Library -- Image Cards Manipulations
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU Lesser General Public License.
+ *
+ *     FIXME:
+ *     - improve thumbnail creation in gatherer... faster compression,
+ *       only grayscale/RGB colorspaces and maybe fixed headers (abbreviated datastreams in libjpeg)
+ *     - hook memory allocation managers, get rid of multiple initializations
+ *     - supply background color to transparent PNG images
+ *     - optimize decompression parameters
+ *     - create interface for thumbnail compression (for gatherer) and reading (MUX)
+ *     - benchmatk libraries
+ */
+
+#undef LOCAL_DEBUG
+
+#include "sherlock/sherlock.h"
+#include "lib/base224.h"
+#include "lib/mempool.h"
+#include "sherlock/object.h"
+#include "images/images.h"
+#include "images/image-obj.h"
+
+#include <stdio.h>
+#include <alloca.h>
+
+/* Selection of libraries to use */
+#define USE_LIBPNG
+#define USE_LIBJPEG
+#define USE_MAGICK
+
+#if defined(USE_LIBPNG) && defined(USE_LIBJPEG)
+#undef USE_MAGICK
+#endif
+
+
+/*********************************  LIBPNG Library ****************************************/
+
+#ifdef USE_LIBPNG
+
+#include <png.h>
+#include <setjmp.h>
+
+static struct mempool *libpng_pool;
+static byte *libpng_buf;
+static uns libpng_len;
+
+static png_voidp
+libpng_malloc(png_structp png_ptr UNUSED, png_size_t size)
+{
+  DBG("libpng_malloc(): size=%d", (uns)size);
+  return mp_alloc(libpng_pool, size);
+}
+
+static void
+libpng_free(png_structp png_ptr UNUSED, png_voidp ptr UNUSED)
+{
+  DBG("libpng_free()");
+}
+
+static void NONRET
+libpng_error(png_structp png_ptr, png_const_charp msg UNUSED)
+{
+  DBG("libpng_error(): msg=%s", (byte *)msg);
+  longjmp(png_jmpbuf(png_ptr), 1);
+}
+
+static void
+libpng_warning(png_structp png_ptr UNUSED, png_const_charp msg UNUSED)
+{
+  DBG("libpng_warning(): msg=%s", (byte *)msg);
+}
+
+static void
+libpng_read_data(png_structp png_ptr UNUSED, png_bytep data, png_size_t length)
+{
+  DBG("libpng_read_data(): len=%d", (uns)length);
+  if (unlikely(libpng_len < length))
+    png_error(png_ptr, "Incomplete data");
+  memcpy(data, libpng_buf, length);
+  libpng_buf += length;
+  libpng_len -= length;
+}
+
+static inline void
+libpng_decompress_thumbnails_init(void)
+{
+}
+
+static inline void
+libpng_decompress_thumbnails_done(void)
+{
+}
+
+static int
+libpng_decompress_thumbnail(struct image_obj *imo)
+{
+  /* create libpng read structure */
+  DBG("Creating libpng read structure");
+  libpng_pool = imo->pool;
+  libpng_buf = imo->thumb_data;
+  libpng_len = imo->thumb_size;
+  png_structp png_ptr = png_create_read_struct_2(PNG_LIBPNG_VER_STRING,
+      NULL, libpng_error, libpng_warning,
+      NULL, libpng_malloc, libpng_free);
+  if (unlikely(!png_ptr))
+    return 0;
+  png_infop info_ptr = png_create_info_struct(png_ptr);
+  if (unlikely(!info_ptr))
+    {
+      png_destroy_read_struct(&png_ptr, NULL, NULL);
+      return 0;
+    }
+  png_infop end_ptr = png_create_info_struct(png_ptr);
+  if (unlikely(!end_ptr))
+    {
+      png_destroy_read_struct(&png_ptr, &info_ptr, NULL);
+      return 0;
+    }
+  if (setjmp(png_jmpbuf(png_ptr)))
+    {
+      DBG("Libpng failed to read the image, longjump saved us");
+      png_destroy_read_struct(&png_ptr, &info_ptr, &end_ptr);
+      return 0;
+    }
+  png_set_read_fn(png_ptr, NULL, libpng_read_data);
+
+  /* Read image info */
+  DBG("Reading image info");
+  png_read_info(png_ptr, info_ptr);
+  png_uint_32 width, height;
+  int bit_depth, color_type;
+  png_get_IHDR(png_ptr, info_ptr, &width, &height, &bit_depth, &color_type, NULL, NULL, NULL);
+  ASSERT(width == imo->thumb.width && height == imo->thumb.height);
+
+  /* Apply transformations */
+  imo->thumb.flags = 0;
+  if (bit_depth == 16)
+    png_set_strip_16(png_ptr);
+  switch (color_type)
+    {
+      case PNG_COLOR_TYPE_PALETTE:
+       png_set_palette_to_rgb(png_ptr);
+        png_set_strip_alpha(png_ptr);
+       break;
+      case PNG_COLOR_TYPE_GRAY:
+       imo->thumb.flags |= IMAGE_GRAYSCALE;
+        png_set_gray_to_rgb(png_ptr);
+       break;
+      case PNG_COLOR_TYPE_GRAY_ALPHA:
+       imo->thumb.flags |= IMAGE_GRAYSCALE;
+        png_set_gray_to_rgb(png_ptr);
+        png_set_strip_alpha(png_ptr);
+       break;
+      case PNG_COLOR_TYPE_RGB:
+       break;
+      case PNG_COLOR_TYPE_RGB_ALPHA:
+        png_set_strip_alpha(png_ptr);
+       break;
+      default:
+       ASSERT(0);
+    }
+  png_read_update_info(png_ptr, info_ptr);
+  ASSERT(png_get_channels(png_ptr, info_ptr) == 3);
+
+  /* Read image data */
+  DBG("Reading image data");
+  byte *pixels = imo->thumb.pixels = mp_alloc(imo->pool, imo->thumb.size = width * height * 3);
+  png_bytep rows[height];
+  for (uns i = 0; i < height; i++, pixels += width * 3)
+    rows[i] = (png_bytep)pixels;
+  png_read_image(png_ptr, rows);
+  png_read_end(png_ptr, end_ptr);
+  /* Destroy libpng read structure */
+  png_destroy_read_struct(&png_ptr, &info_ptr, &end_ptr);
+  return 1;
+}
+
+#endif /* USE_LIBPNG */
+
+
+
+/*******************************  LIBJPEG Library *************************************/
+
+#ifdef USE_LIBJPEG
+
+#include <jpeglib.h>
+#include <setjmp.h>
+
+struct libjpeg_err {
+  struct jpeg_error_mgr pub;
+  jmp_buf setjmp_buf;
+};
+
+static void NONRET
+libjpeg_error_exit(j_common_ptr cinfo)
+{
+  DBG("libjpeg_error_exit()");
+  longjmp(((struct libjpeg_err *)(cinfo)->err)->setjmp_buf, 1);
+}
+
+static void
+libjpeg_emit_message(j_common_ptr cinfo UNUSED, int msg_level UNUSED)
+{
+  DBG("libjpeg_emit_message(): level=%d", msg_level);
+  /* if (unlikely(msg_level == -1))
+    longjmp(((struct libjpeg_err *)(cinfo)->err)->setjmp_buf, 1); */
+}
+
+static void
+libjpeg_init_source(j_decompress_ptr cinfo UNUSED)
+{
+  DBG("libjpeg_init_source()");
+}
+
+static boolean NONRET
+libjpeg_fill_input_buffer(j_decompress_ptr cinfo)
+{ 
+  DBG("libjpeg_fill_input_buffer()");
+  longjmp(((struct libjpeg_err *)(cinfo)->err)->setjmp_buf, 1);
+}
+
+static void
+libjpeg_skip_input_data(j_decompress_ptr cinfo, long num_bytes) 
+{
+  DBG("libjpeg_skip_input_data(): len=%d", (int)num_bytes);
+  if (num_bytes > 0)
+    {
+      cinfo->src->next_input_byte += num_bytes;
+      cinfo->src->bytes_in_buffer -= num_bytes;
+    }
+}
+
+static inline void
+libjpeg_decompress_thumbnails_init(void)
+{
+}
+
+static inline void
+libjpeg_decompress_thumbnails_done(void)
+{
+}
+
+static int
+libjpeg_decompress_thumbnail(struct image_obj *imo)
+{
+  /* Create libjpeg read structure */
+  DBG("Creating libjpeg read structure");
+  struct jpeg_decompress_struct cinfo;
+  struct libjpeg_err err;
+  cinfo.err = jpeg_std_error(&err.pub);
+  err.pub.error_exit = libjpeg_error_exit;
+  err.pub.emit_message = libjpeg_emit_message;
+  if (setjmp(err.setjmp_buf))
+    {
+      DBG("Libjpeg failed to read the image, longjump saved us");
+      jpeg_destroy_decompress(&cinfo);
+      return 0;
+    }
+  jpeg_create_decompress(&cinfo);
+
+  /* Initialize source manager */
+  struct jpeg_source_mgr src;
+  cinfo.src = &src;
+  src.next_input_byte = imo->thumb_data;
+  src.bytes_in_buffer = imo->thumb_size;
+  src.init_source = libjpeg_init_source;
+  src.fill_input_buffer = libjpeg_fill_input_buffer;
+  src.skip_input_data = libjpeg_skip_input_data;
+  src.resync_to_restart = jpeg_resync_to_restart;
+  src.term_source = libjpeg_init_source;
+
+  /* Read JPEG header and setup decompression options */
+  DBG("Reading image header");
+  jpeg_read_header(&cinfo, TRUE);
+  imo->thumb.flags = 0;
+  if (cinfo.out_color_space == JCS_GRAYSCALE)
+    imo->thumb.flags |= IMAGE_GRAYSCALE;
+  else
+    cinfo.out_color_space = JCS_RGB;
+
+  /* Decompress the image */
+  DBG("Reading image data");
+  jpeg_start_decompress(&cinfo);
+  ASSERT(imo->thumb.width == cinfo.output_width && imo->thumb.height == cinfo.output_height);
+  ASSERT(sizeof(JSAMPLE) == 1);
+  byte *pixels = imo->thumb.pixels = mp_alloc(imo->pool, imo->thumb.size = cinfo.output_width * cinfo.output_height * 3);
+  if (cinfo.out_color_space == JCS_RGB)
+    { /* Read RGB pixels */
+      uns size = cinfo.output_width * 3;
+      while (cinfo.output_scanline < cinfo.output_height)
+        {
+          jpeg_read_scanlines(&cinfo, (JSAMPLE **)&pixels, 1);
+         pixels += size;
+        }
+    }
+  else
+    { /* Read grayscale pixels */
+      JSAMPLE buf[cinfo.output_width], *buf_end = buf + cinfo.output_width;
+      while (cinfo.output_scanline < cinfo.output_height)
+        {
+          JSAMPLE *p = buf;
+          jpeg_read_scanlines(&cinfo, &p, 1);
+          for (; p != buf_end; p++)
+            {
+              pixels[0] = pixels[1] = pixels[2] = p[0];
+             pixels += 3;
+           }
+        }
+    }
+  jpeg_finish_decompress(&cinfo);
+
+  /* Destroy libjpeg object and leave */
+  jpeg_destroy_decompress(&cinfo);
+  return 1;
+}
+
+#endif /* USE_LIBJPEG */
+
+
+
+/******************************  GraphicsMagick Library ******************************/
+
+#ifdef USE_MAGICK
+
+#include <magick/api.h>
+
+static ExceptionInfo magick_exception;
+static QuantizeInfo magick_quantize;
+static ImageInfo *magick_info;
+
+static void
+magick_decompress_thumbnails_init(void)
+{
+  DBG("Initializing magick thumbnails decompression");
+  InitializeMagick(NULL);
+  GetExceptionInfo(&magick_exception);
+  magick_info = CloneImageInfo(NULL);
+  magick_info->subrange = 1;
+  GetQuantizeInfo(&magick_quantize);
+  magick_quantize.colorspace = RGBColorspace;
+}
+
+static void
+magick_decompress_thumbnails_done(void)
+{
+  DBG("Finalizing magick thumbnails decompression");
+  DestroyImageInfo(magick_info);
+  DestroyExceptionInfo(&magick_exception);
+  DestroyMagick();
+}
+
+static int
+magick_decompress_thumbnail(struct image_obj *imo)
+{
+  DBG("Reading image data");
+  Image *image = BlobToImage(magick_info, imo->thumb_data, imo->thumb_size, &magick_exception);
+  if (unlikely(!image))
+    return 0;
+  ASSERT(image->columns == imo->thumb.width && image->rows == imo->thumb.height);
+  DBG("Quantizing image");
+  QuantizeImage(&magick_quantize, image);
+  DBG("Converting pixels");
+  PixelPacket *pixels = (PixelPacket *)AcquireImagePixels(image, 0, 0, image->columns, image->rows, &magick_exception);
+  ASSERT(pixels);
+  uns size = image->columns * image->rows;
+  byte *p = imo->thumb.pixels = mp_alloc(imo->pool, imo->thumb.size = size * 3);
+  for (uns i = 0; i < size; i++)
+    {
+      p[0] = pixels->red >> (QuantumDepth - 8);
+      p[1] = pixels->green >> (QuantumDepth - 8);
+      p[2] = pixels->blue >> (QuantumDepth - 8);
+      p += 3;
+      pixels++;
+    }
+  DestroyImage(image);
+  return 1;
+}
+
+#endif /* USE_MAGICK */
+
+
+
+/*************************************************************************************/
+
+static int
+extract_image_info(struct image_obj *imo)
+{
+  DBG("Parsing image info attribute");
+  ASSERT(!(imo->flags & IMAGE_OBJ_VALID_INFO));
+  imo->flags |= IMAGE_OBJ_VALID_INFO;
+  byte *info = obj_find_aval(imo->obj, 'G');
+  if (!info)
+    {
+      DBG("Attribute G not found");
+      return 0;
+    }
+  uns colors;
+  byte color_space[MAX_ATTR_SIZE], thumb_format[MAX_ATTR_SIZE];
+  UNUSED uns cnt = sscanf(info, "%d%d%s%d%d%d%s", &imo->width, &imo->height, color_space, &colors, &imo->thumb.width, &imo->thumb.height, thumb_format);
+  ASSERT(cnt == 7);
+  switch (*thumb_format)
+    {
+      case 'j':
+       imo->thumb_format = IMAGE_OBJ_FORMAT_JPEG;
+       break;
+      case 'p':
+       imo->thumb_format = IMAGE_OBJ_FORMAT_PNG;
+       break;
+      default:
+       ASSERT(0);
+    }
+  return 1;
+}
+
+static int
+extract_thumb_data(struct image_obj *imo)
+{
+  DBG("Extracting thumbnail data");
+  ASSERT(!(imo->flags & IMAGE_OBJ_VALID_DATA) && 
+      (imo->flags & IMAGE_OBJ_VALID_INFO));
+  imo->flags |= IMAGE_OBJ_VALID_DATA;
+  struct oattr *attr = obj_find_attr(imo->obj, 'N');
+  if (!attr)
+    {
+      DBG("There is no thumbnail attribute N");
+      return 0;
+    }
+  uns count = 0;
+  for (struct oattr *a = attr; a; a = a->same)
+    count++;
+  byte b224[count * MAX_ATTR_SIZE], *b = b224;
+  for (struct oattr *a = attr; a; a = a->same)
+    for (byte *s = a->val; *s; )
+      *b++ = *s++;
+  ASSERT(b != b224);
+  uns size = b - b224;
+  imo->thumb_data = mp_alloc(imo->pool, size);
+  imo->thumb_size = base224_decode(imo->thumb_data, b224, size);
+  DBG("Thumbnail data size is %d", imo->thumb_size);
+  return 1;
+}
+
+static int
+extract_thumb_image(struct image_obj *imo)
+{
+  DBG("Decompressing thumbnail image");
+  ASSERT(!(imo->flags & IMAGE_OBJ_VALID_IMAGE) &&
+      (imo->flags & IMAGE_OBJ_VALID_INFO) &&
+      (imo->flags & IMAGE_OBJ_VALID_DATA));
+  imo->flags |= IMAGE_OBJ_VALID_IMAGE;
+  switch (imo->thumb_format)
+    {
+      case IMAGE_OBJ_FORMAT_JPEG:
+#if defined(USE_LIBJPEG)
+        return libjpeg_decompress_thumbnail(imo);
+#elif defined(USE_MAGICK)
+        return magick_decompress_thumbnail(imo);
+#else
+        DBG("JPEG not supported");
+        return 0;
+#endif
+      case IMAGE_OBJ_FORMAT_PNG:
+#if defined(USE_LIBPNG)
+        return libpng_decompress_thumbnail(imo);
+#elif defined(USE_MAGICK)
+        return magick_decompress_thumbnail(imo);
+#else
+        DBG("PNG not supported");
+        return 0;
+#endif      
+      default:
+       ASSERT(0);
+    }
+}
+
+void
+imo_decompress_thumbnails_init(void)
+{
+#ifdef USE_LIBPNG
+  libpng_decompress_thumbnails_init();
+#endif
+#ifdef USE_LIBJPEG
+  libjpeg_decompress_thumbnails_init();
+#endif
+#ifdef USE_MAGICK
+  magick_decompress_thumbnails_init();
+#endif
+}
+
+void
+imo_decompress_thumbnails_done(void)
+{
+#ifdef USE_MAGICK
+  magick_decompress_thumbnails_done();
+#endif
+#ifdef USE_LIBJPEG
+  libjpeg_decompress_thumbnails_done();
+#endif
+#ifdef USE_LIBPNG
+  libpng_decompress_thumbnails_done();
+#endif
+}
+
+int
+imo_decompress_thumbnail(struct image_obj *imo)
+{
+  return
+    extract_image_info(imo) &&
+    extract_thumb_data(imo) &&
+    extract_thumb_image(imo);
+}
+
diff --git a/images/image-obj.h b/images/image-obj.h
new file mode 100644 (file)
index 0000000..95976bc
--- /dev/null
@@ -0,0 +1,53 @@
+/*
+ *     Image Library -- Image Cards Manipulations
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU Lesser General Public License.
+ */
+
+#ifndef _IMAGES_IMAGE_OBJ_H
+#define _IMAGES_IMAGE_OBJ_H
+
+#include "images/images.h"
+
+struct mempool;
+struct odes;
+
+enum image_obj_format {
+  IMAGE_OBJ_FORMAT_JPEG,
+  IMAGE_OBJ_FORMAT_PNG
+};
+
+enum image_obj_flag {
+  IMAGE_OBJ_VALID_INFO = 0x1,
+  IMAGE_OBJ_VALID_DATA = 0x2,
+  IMAGE_OBJ_VALID_IMAGE = 0x4
+};
+
+struct image_obj {
+  struct odes *obj;
+  struct mempool *pool;
+  uns flags;
+  uns width;
+  uns height;
+  uns thumb_format;
+  byte *thumb_data;
+  uns thumb_size;
+  struct image thumb;
+};
+
+static inline void
+imo_init(struct image_obj *imo, struct mempool *pool, struct odes *obj)
+{
+  imo->obj = obj;
+  imo->pool = pool;
+  imo->flags = 0;
+}
+
+void imo_decompress_thumbnails_init(void);
+void imo_decompress_thumbnails_done(void);
+int imo_decompress_thumbnail(struct image_obj *imo);
+
+#endif
diff --git a/images/image-sig.c b/images/image-sig.c
new file mode 100644 (file)
index 0000000..cd07b8a
--- /dev/null
@@ -0,0 +1,137 @@
+#undef LOCAL_DEBUG
+
+#include "sherlock/sherlock.h"
+#include "lib/math.h"
+#include "lib/fastbuf.h"
+#include "images/images.h"
+#include "images/image-obj.h"
+#include "images/image-sig.h"
+#include "images/color.h"
+
+#include <alloca.h>
+
+struct block {
+  uns l, u, v;         /* average Luv coefficients */
+  uns lh, hl, hh;      /* energies in Daubechies wavelet bands */
+};
+
+int
+compute_image_signature(struct image_data *image, struct image_signature *sig)
+{
+  uns width = image->width;
+  uns height = image->height;
+
+  if (width < 4 || height < 4)
+    {
+      DBG("Image too small... %dx%d", width, height);
+      return 0;
+    }
+  
+  uns w = width >> 2;
+  uns h = height >> 2;
+  DBG("Computing signature for image %dx%d... %dx%d blocks", width, height, w, h);
+  uns blocks_count = w * h;
+  struct block *blocks = xmalloc(blocks_count * sizeof(struct block)), *block = blocks; /* FIXME: use mempool */
+  
+  /* Every 4x4 block (FIXME: deal with smaller blocks near the edges) */
+  byte *p = image->pixels;
+  for (uns block_y = 0; block_y < h; block_y++, p += 3 * ((width & 3) + width * 3))
+    for (uns block_x = 0; block_x < w; block_x++, p -= 3 * (4 * width - 4), block++)
+      {
+        int t[16], s[16], *tp = t;
+
+       /* Convert pixels to Luv color space and compute average coefficients 
+        * FIXME:
+        * - could be MUCH faster with precomputed tables and integer arithmetic... 
+        *   I will propably use interpolation in 3-dim array */
+       uns l_sum = 0;
+       uns u_sum = 0;
+       uns v_sum = 0;
+       for (uns y = 0; y < 4; y++, p += 3 * (width - 4))
+         for (uns x = 0; x < 4; x++, p += 3)
+           {
+             byte luv[3];
+             srgb_to_luv_pixel(luv, p);
+             l_sum += *tp++ = luv[0];
+             u_sum += luv[1];
+             v_sum += luv[2];
+           }
+
+       block->l = (l_sum >> 4);
+       block->u = (u_sum >> 4);
+       block->v = (v_sum >> 4);
+
+       /* Apply Daubechies wavelet transformation 
+        * FIXME:
+        * - MMX/SSE instructions or tables could be faster 
+        * - maybe it would be better to compute Luv and wavelet separately because of processor cache or MMX/SSE 
+        * - eliminate slow square roots 
+        * - what about Haar transformation? */
+
+#define DAUB_0 31651 /* (1 + sqrt 3) / (4 * sqrt 2) */
+#define DAUB_1 54822 /* (3 + sqrt 3) / (4 * sqrt 2) */
+#define DAUB_2 14689 /* (3 - sqrt 3) / (4 * sqrt 2) */
+#define DAUB_3 -8481 /* (1 - sqrt 3) / (4 * sqrt 2) */
+
+       /* ... to the rows */
+       uns i;
+        for (i = 0; i < 16; i += 4)
+          {
+           s[i + 0] = (DAUB_0 * t[i + 2] + DAUB_1 * t[i + 3] + DAUB_2 * t[i + 0] + DAUB_3 * t[i + 1]) / 0x10000;
+           s[i + 1] = (DAUB_0 * t[i + 0] + DAUB_1 * t[i + 1] + DAUB_2 * t[i + 2] + DAUB_3 * t[i + 3]) / 0x10000;
+           s[i + 2] = (DAUB_3 * t[i + 2] - DAUB_2 * t[i + 3] + DAUB_1 * t[i + 0] - DAUB_0 * t[i + 1]) / 0x10000;
+           s[i + 3] = (DAUB_3 * t[i + 0] - DAUB_2 * t[i + 1] + DAUB_1 * t[i + 2] - DAUB_0 * t[i + 3]) / 0x10000;
+         }
+
+       /* ... and to the columns... skip LL band */
+       for (i = 0; i < 2; i++)
+         {
+           t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
+           t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
+         }
+       for (; i < 4; i++)
+         {
+           t[i + 0] = (DAUB_0 * s[i + 8] + DAUB_1 * s[i +12] + DAUB_2 * s[i + 0] + DAUB_3 * s[i + 4]) / 0x1000;
+           t[i + 4] = (DAUB_0 * s[i + 0] + DAUB_1 * s[i + 4] + DAUB_2 * s[i + 8] + DAUB_3 * s[i +12]) / 0x1000;
+           t[i + 8] = (DAUB_3 * s[i + 8] - DAUB_2 * s[i +12] + DAUB_1 * s[i + 0] - DAUB_0 * s[i + 4]) / 0x1000;
+           t[i +12] = (DAUB_3 * s[i + 0] - DAUB_2 * s[i + 4] + DAUB_1 * s[i + 8] - DAUB_0 * s[i +12]) / 0x1000;
+         }
+
+       /* Extract energies in LH, HL and HH bands */
+       block->lh = CLAMP((int)(sqrt(t[8] * t[8] + t[9] * t[9] + t[12] * t[12] + t[13] * t[13]) / 16), 0, 255);
+       block->hl = CLAMP((int)(sqrt(t[2] * t[2] + t[3] * t[3] + t[6] * t[6] + t[7] * t[7]) / 16), 0, 255);
+       block->hh = CLAMP((int)(sqrt(t[10] * t[10] + t[11] * t[11] + t[14] * t[14] + t[15] * t[15]) / 16), 0, 255);
+      }
+
+  /* FIXME: simple average is for testing pusposes only */
+  uns l_sum = 0;
+  uns u_sum = 0;
+  uns v_sum = 0;
+  uns lh_sum = 0;
+  uns hl_sum = 0;
+  uns hh_sum = 0;
+  for (uns i = 0; i < blocks_count; i++)
+    {
+      l_sum += blocks[i].l;
+      u_sum += blocks[i].u;
+      v_sum += blocks[i].v;
+      lh_sum += blocks[i].lh;
+      hl_sum += blocks[i].hl;
+      hh_sum += blocks[i].hh;
+    }
+
+  sig->vec.f[0] = l_sum / blocks_count;
+  sig->vec.f[1] = u_sum / blocks_count;
+  sig->vec.f[2] = v_sum / blocks_count;
+  sig->vec.f[3] = lh_sum / blocks_count;
+  sig->vec.f[4] = hl_sum / blocks_count;
+  sig->vec.f[5] = hh_sum / blocks_count;
+
+  sig->len = 0;
+
+  xfree(blocks);
+
+  DBG("Resulting signature is (%s)", stk_print_image_vector(&sig->vec));
+  return 1;
+}
+
diff --git a/images/image-sig.h b/images/image-sig.h
new file mode 100644 (file)
index 0000000..3532262
--- /dev/null
@@ -0,0 +1,65 @@
+#ifndef _IMAGES_IMAGE_SIG_H
+#define _IMAGES_IMAGE_SIG_H
+
+#include "images/images.h"
+
+#define IMAGE_VEC_K    6
+#define IMAGE_REG_K    9
+#define IMAGE_REG_MAX  4
+
+typedef byte image_feature_t;  /* 8 or 16 bits precision */
+
+/* K-dimensional feature vector */
+struct image_vector {
+  image_feature_t f[IMAGE_VEC_K];
+};
+
+/* K-dimensional interval */
+struct image_bbox {
+  struct image_vector vec[2];
+};
+
+/* Fetures for image regions */
+struct image_region {
+  image_feature_t f[IMAGE_REG_K];
+};
+
+/* Image signature */
+struct image_signature {
+  struct image_vector vec;     /* Combination of all regions... simplier signature */
+  image_feature_t len;         /* Number of regions */
+  struct image_region reg[IMAGE_REG_MAX];/* Feature vector for every region */
+};
+
+/* Similarity search tree... will be changed */
+struct image_tree {
+  uns count;                   /* Number of images in the tree */
+  uns depth;                   /* Tree depth */
+  struct image_bbox bbox;      /* Bounding box containing all the */
+  struct image_node *nodes;    /* Internal nodes */
+  struct image_leaf *leaves;   /* Leaves */
+};
+
+/* Internal node in the search tree */
+#define IMAGE_NODE_LEAF                0x80000000              /* Node contains pointer to leaves array */
+#define IMAGE_NODE_DIM         0xff                    /* Split dimension */
+struct image_node {
+  u32 val;
+};
+
+/* Leaves in the search tree */
+#define IMAGE_LEAF_LAST                0x80000000              /* Last entry in the list */
+#define IMAGE_LEAF_BITS(i)     (31 / IMAGE_VEC_K)      /* Number of bits for relative position in i-th dimension */
+struct image_leaf {
+  u32 flags;           /* Relative position in bbox and last node flag */ 
+  oid_t oid;
+};
+
+#define stk_print_image_vector(v) ({ struct image_vector *_v = v; \
+    byte *_s = (byte *) alloca(IMAGE_VEC_K * 6), *_p = _s + sprintf(_s, "%d", _v->f[0]); \
+    for (uns _i = 1; _i < IMAGE_VEC_K; _i++) _p += sprintf(_p, " %d", _v->f[_i]); _s; })
+
+int compute_image_signature(struct image *image, struct image_signature *sig);
+
+#endif
+
diff --git a/images/image-test.c b/images/image-test.c
new file mode 100644 (file)
index 0000000..b909ed6
--- /dev/null
@@ -0,0 +1,121 @@
+#undef LOCAL_DEBUG
+
+#include "sherlock/sherlock.h"
+#include "lib/fastbuf.h"
+#include "images/images.h"
+//#include "images/image-sig.h"
+//#include "images/kd-tree.h"
+#include "sherlock/index.h"
+#include "lib/mempool.h"
+#include "sherlock/object.h"
+#include "sherlock/lizard-fb.h"
+#include <fcntl.h>
+#include <stdio.h>
+
+#include "sherlock/sherlock.h"
+#include "lib/fastbuf.h"
+//#include "images/images.h"
+#include "sherlock/index.h"
+
+#include <stdio.h>
+#include <fcntl.h>
+#include <alloca.h>
+
+#define BEST_CNT 30
+
+struct image_tree image_tree;
+
+static void
+image_tree_init(void)
+{
+  DBG("Initializing image search structures");
+  struct fastbuf *fb = bopen("index/image-tree", O_RDONLY, 1 << 16); /* FIXME: filename hack */
+  image_tree.count = bgetl(fb);
+  image_tree.depth = bgetl(fb);
+  ASSERT(image_tree.count < 0x80000000 && image_tree.depth > 0 && image_tree.depth < 30);
+  image_tree.nodes = xmalloc((1 << image_tree.depth) * sizeof(struct image_node));
+  image_tree.leaves = xmalloc(image_tree.count * sizeof(struct image_leaf));
+  bread(fb, &image_tree.bbox, sizeof(struct image_bbox));
+  bread(fb, image_tree.nodes + 1, ((1 << image_tree.depth) - 1) * sizeof(struct image_node));
+  bread(fb, image_tree.leaves, image_tree.count * sizeof(struct image_leaf));
+  DBG("Search tree with depth %d and %d leaves loaded", image_tree.depth, image_tree.count);
+  bclose(fb);
+}
+
+static void
+image_tree_done(void)
+{
+  DBG("Freeing image search structures");
+  xfree(image_tree.nodes);
+  xfree(image_tree.leaves);
+}
+  
+int
+main(int argc, char **argv)
+{
+  struct image_vector query;
+  if (argc != IMAGE_VEC_K + 1)
+    die("Invalid number of arguments");
+
+  for (uns i = 0; i < IMAGE_VEC_K; i++)
+    {
+      uns v;
+      if (sscanf(argv[i + 1], "%d", &v) != 1)
+       die("Invalid numeric format");
+      query.f[i] = v;
+    }
+  
+  
+  struct image_search is;
+  oid_t best[BEST_CNT];
+  uns dist[BEST_CNT];
+  uns cardpos[BEST_CNT];
+  uns best_n = 0;
+  
+  image_tree_init();
+  log(L_INFO, "Executing query (%s)", stk_print_image_vector(&query));
+  image_search_init(&is, &image_tree, &query, IMAGE_SEARCH_DIST_UNLIMITED);
+  for (uns i = 0; i < BEST_CNT; i++)
+    {
+      if (!image_search_next(&is, best + i, dist + i))
+        {
+          log(L_INFO, "No more images");
+          break;
+        }
+      DBG("*** Found %d. best image with oid=%d", i + 1, best[i]);
+      best_n++;
+    }
+  image_search_done(&is);
+  image_tree_done();
+
+  log(L_INFO, "Resolving URLs");
+  struct mempool *pool = mp_new(1 << 16);
+  struct buck2obj_buf *bob = buck2obj_alloc();
+  struct fastbuf *fb = bopen("index/card-attrs", O_RDONLY, 1 << 10);
+  for (uns i = 0; i < best_n; i++)
+    {
+      bsetpos(fb, best[i] * sizeof(struct card_attr));
+      struct card_attr ca;
+      bread(fb, &ca, sizeof(ca));
+      cardpos[i] = ca.card;
+    }
+  bclose(fb);
+  fb = bopen("index/cards", O_RDONLY, 1 << 14);
+  for (uns i = 0; i < best_n; i++)
+    {
+      bsetpos(fb, (sh_off_t)cardpos[i] << CARD_POS_SHIFT);
+      uns buck_len = bgetl(fb) - (LIZARD_COMPRESS_HEADER - 1);
+      uns buck_type = bgetc(fb) + BUCKET_TYPE_PLAIN;
+      mp_flush(pool);
+      struct odes *obj = obj_read_bucket(bob, pool, buck_type, buck_len, fb, NULL);
+
+      printf("%2d. match: dist=%-8d oid=%-8d url=%s\n", i + 1, dist[i], best[i],
+         obj_find_aval(obj_find_attr(obj, 'U' + OBJ_ATTR_SON)->son, 'U'));
+    }
+  bclose(fb);
+  buck2obj_free(bob);
+  mp_delete(pool);
+     
+  return 0;
+}
+
diff --git a/images/image-tool.c b/images/image-tool.c
new file mode 100644 (file)
index 0000000..94fdd44
--- /dev/null
@@ -0,0 +1,162 @@
+/*
+ *     Simple image manupulation utility
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU General Public License.
+ */
+
+#include "lib/lib.h"
+#include "lib/getopt.h"
+#include "lib/fastbuf.h"
+#include "images/images.h"
+#include <stdlib.h>
+#include <fcntl.h>
+
+static void NONRET
+usage(void)
+{
+  fputs("\
+Usage: image-tool [options] infile [outfile]\n\
+\n\
+-q --quiet          no progress messages\n\
+-f --input-format   input image format (jpeg, gif, png)\n\
+-F --output-format  output image format\n\
+-s --size           force output dimensions (100x200)\n\
+-b --fit-box        scale to fit the box (100x200)\n\
+-c --colorspace     force output colorspace (gray, grayalpha, rgb, rgbalpha)\n\
+", stderr);
+  exit(1);
+}
+
+static char *shortopts = "qf:F:s:b:c:" CF_SHORT_OPTS;
+static struct option longopts[] =
+{
+  CF_LONG_OPTS
+  { "quiet",           0, 0, 'q' },
+  { "input-format",    0, 0, 'f' },
+  { "output-format",   0, 0, 'F' },
+  { "size",            0, 0, 's' },
+  { "fit-box",         0, 0, 'b' },
+  { "colorspace",      0, 0, 'c' },
+  { NULL,              0, 0, 0 }
+};
+                                                         
+static uns verbose = 1;
+static byte *input_file_name;
+static enum image_format input_format;
+static byte *output_file_name;
+static enum image_format output_format;
+static uns cols;
+static uns rows;
+static uns fit_box;
+static uns channels_format;
+
+#define MSG(x...) do{ if (verbose) log(L_INFO, ##x); }while(0)
+
+int
+main(int argc, char **argv)
+{
+  log_init(argv[0]);
+  int opt;
+  while ((opt = cf_getopt(argc, argv, shortopts, longopts, NULL)) >= 0)
+    switch (opt)
+      {
+       case 'q':
+         verbose = 0;
+         break;
+       case 'f':
+         if (!(input_format = image_extension_to_format(optarg)))
+           usage();
+         break;
+       case 'F':
+         if (!(output_format = image_extension_to_format(optarg)))
+           usage();
+         break;
+       case 's':
+         {
+           byte *r = strchr(optarg, 'x');
+           if (!r)
+             usage();
+           *r++ = 0;
+           if (!(cols = atoi(optarg)) || !(rows = atoi(r)))
+             usage();
+           fit_box = 0;
+           break;
+         }
+       case 'b':
+         {
+           byte *r = strchr(optarg, 'x');
+           if (!r)
+             usage();
+           *r++ = 0;
+           if (!(cols = atoi(optarg)) || !(rows = atoi(r)))
+             usage();
+           fit_box = 1;
+           break;
+         }
+       case 'c':
+         if (!(channels_format = image_name_to_channels_format(optarg)))
+           usage();
+         break;
+       default:
+         usage();
+      }
+
+  if (argc != optind + 1 && argc != optind + 2)
+    usage();
+  input_file_name = argv[optind++];
+  if (argc > optind)
+    output_file_name = argv[optind];
+  
+#define TRY(x) do{ if (!(x)) die("Error: %s", it.err_msg); }while(0)
+  MSG("Initializing image library");
+  struct image_thread it;
+  struct image_io io;
+  image_thread_init(&it);
+  image_io_init(&it, &io);
+
+  MSG("Reading %s", input_file_name);
+  io.fastbuf = bopen(input_file_name, O_RDONLY, 1 << 18);
+  io.format = input_format ? : image_file_name_to_format(input_file_name);
+  TRY(image_io_read_header(&io));
+  if (!output_file_name)
+    {
+      bclose(io.fastbuf);
+      printf("Format:      %s\n", image_format_to_extension(io.format) ? : (byte *)"?");
+      printf("Dimensions:  %dx%d\n", io.cols, io.rows);
+      printf("Colorspace:  %s\n", image_channels_format_to_name(io.flags & IMAGE_CHANNELS_FORMAT));
+    }
+  else
+    {
+      MSG("%s %dx%d %s", image_format_to_extension(io.format) ? : (byte *)"?", io.cols, io.rows,
+         image_channels_format_to_name(io.flags & IMAGE_CHANNELS_FORMAT));
+      if (cols)
+        if (fit_box)
+         {
+            image_dimensions_fit_to_box(&io.cols, &io.rows, MIN(cols, 0xffff), MIN(rows, 0xffff), 0);
+         }
+        else
+          {
+            io.cols = cols;
+            io.rows = rows;
+          }
+      if (channels_format)
+        io.flags = io.flags & ~IMAGE_PIXEL_FORMAT | channels_format;
+      TRY(image_io_read_data(&io, 0));
+      bclose(io.fastbuf);
+      MSG("Writing %s", output_file_name);
+      io.fastbuf = bopen(output_file_name, O_WRONLY | O_CREAT | O_TRUNC, 1 << 18);
+      io.format = output_format ? : image_file_name_to_format(output_file_name);
+      MSG("%s %dx%d %s", image_format_to_extension(io.format) ? : (byte *)"?", io.cols, io.rows,
+         image_channels_format_to_name(io.flags & IMAGE_CHANNELS_FORMAT));
+      TRY(image_io_write(&io));
+      bclose(io.fastbuf);
+    }
+  
+  image_io_cleanup(&io);
+  image_thread_cleanup(&it);
+  MSG("Done.");
+  return 0;
+}
diff --git a/images/image-walk.h b/images/image-walk.h
new file mode 100644 (file)
index 0000000..a276b08
--- /dev/null
@@ -0,0 +1,157 @@
+/*
+ *     Image Library -- Pixels iteration
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU Lesser General Public License.
+ */
+
+#if !defined(IMAGE_WALK_INLINE) && !defined(IMAGE_WALK_STATIC)
+#  error Missing IMAGE_WALK_INLINE or IMAGE_WALK_STATIC
+#endif
+
+#if !defined(IMAGE_WALK_UNROLL)
+#  define IMAGE_WALK_UNROLL 1
+#elif IMAGE_WALK_UNROLL != 1 && IMAGE_WALK_UNROLL != 2 && IMAGE_WALK_UNROLL != 4
+#  error IMAGE_WALK_UNROLL must be 1, 2 or 4
+#endif
+
+#ifndef IMAGE_WALK_PIXELS
+#  define IMAGE_WALK_PIXELS (img->pixels)
+#endif
+#ifndef IMAGE_WALK_COLS
+#  define IMAGE_WALK_COLS (img->cols)
+#endif
+#ifndef IMAGE_WALK_ROWS
+#  define IMAGE_WALK_ROWS (img->rows)
+#endif
+#ifndef IMAGE_WALK_COL_STEP
+#  define IMAGE_WALK_COL_STEP (img->pixel_size)
+#endif
+#ifndef IMAGE_WALK_ROW_STEP
+#  define IMAGE_WALK_ROW_STEP (img->row_size)
+#endif
+
+#ifdef IMAGE_WALK_DOUBLE
+#  ifndef IMAGE_WALK_SEC_PIXELS
+#    define IMAGE_WALK_SEC_PIXELS (sec_img->pixels)
+#  endif
+#  ifndef IMAGE_WALK_SEC_COLS
+#    define IMAGE_WALK_SEC_COLS (sec_img->cols)
+#  endif
+#  ifndef IMAGE_WALK_SEC_ROWS
+#    define IMAGE_WALK_SEC_ROWS (sec_img->rows)
+#  endif
+#  ifndef IMAGE_WALK_SEC_COL_STEP
+#    define IMAGE_WALK_SEC_COL_STEP (sec_img->pixel_size)
+#  endif
+#  ifndef IMAGE_WALK_SEC_ROW_STEP
+#    define IMAGE_WALK_SEC_ROW_STEP (sec_img->row_size)
+#  endif
+#  define STEP IMAGE_WALK_DO_STEP; pos += col_step; sec_pos += sec_col_step
+#else
+#  define STEP IMAGE_WALK_DO_STEP; pos += col_step
+#endif
+
+#ifndef IMAGE_WALK_DO_START
+#  define IMAGE_WALK_DO_START
+#endif
+
+#ifndef IMAGE_WALK_DO_END
+#  define IMAGE_WALK_DO_END
+#endif
+
+#ifndef IMAGE_WALK_DO_ROW_START
+#  define IMAGE_WALK_DO_ROW_START
+#endif
+
+#ifndef IMAGE_WALK_DO_ROW_END
+#  define IMAGE_WALK_DO_ROW_END
+#endif
+
+#ifndef IMAGE_WALK_DO_STEP
+#  define IMAGE_WALK_DO_STEP
+#endif
+
+#ifndef IMAGE_WALK_INLINE
+static void IMAGE_WALK_STATIC
+    (struct image *img
+#   ifdef IMAGE_WALK_DOUBLE
+    , struct image *sec_img
+#   endif
+#   ifdef IMAGE_WALK_EXTRA_ARGS
+    , IMAGE_WALK_EXTRA_ARGS
+#   endif
+    )
+#endif
+{
+  uns cols = IMAGE_WALK_COLS;
+  uns rows = IMAGE_WALK_ROWS;
+# if IMAGE_WALK_UNROLL > 1
+  uns cols_unroll_block_count = cols / IMAGE_WALK_UNROLL;
+  uns cols_unroll_end_count = cols % IMAGE_WALK_UNROLL;
+# endif
+  byte *pos = IMAGE_WALK_PIXELS, *row_start = pos;
+  int col_step = IMAGE_WALK_COL_STEP;
+  int row_step = IMAGE_WALK_ROW_STEP;
+# ifdef IMAGE_WALK_DOUBLE
+  byte *sec_pos = IMAGE_WALK_SEC_PIXELS, *sec_row_start = sec_pos;
+  int sec_col_step = IMAGE_WALK_SEC_COL_STEP;
+  int sec_row_step = IMAGE_WALK_SEC_ROW_STEP;
+# endif
+  IMAGE_WALK_DO_START;
+  while (rows--)
+    {
+      IMAGE_WALK_DO_ROW_START;
+#     if IMAGE_WALK_UNROLL == 1
+      for (uns i = cols; i--; )
+#     else
+      for (uns i = cols_unroll_block_count; i--; )
+#     endif
+        {
+#         if IMAGE_WALK_UNROLL >= 4
+         STEP;
+         STEP;
+#         endif
+#         if IMAGE_WALK_UNROLL >= 2
+         STEP;
+#         endif
+         STEP;
+       }
+#     if IMAGE_WALK_UNROLL > 1
+      for (uns i = cols_unroll_end_count; i--; )
+        {
+         STEP;
+       }
+#     endif
+      IMAGE_WALK_DO_ROW_END;
+      pos = (row_start += row_step);
+#     ifdef IMAGE_WALK_DOUBLE
+      sec_pos = (sec_row_start += sec_row_step);
+#     endif
+    }
+  IMAGE_WALK_DO_END;
+}
+
+#undef IMAGE_WALK_INLINE
+#undef IMAGE_WALK_STATIC
+#undef IMAGE_WALK_UNROLL
+#undef IMAGE_WALK_DOUBLE
+#undef IMAGE_WALK_EXTRA_ARGS
+#undef IMAGE_WALK_PIXELS
+#undef IMAGE_WALK_COLS
+#undef IMAGE_WALK_ROWS
+#undef IMAGE_WALK_COL_STEP
+#undef IMAGE_WALK_ROW_STEP
+#undef IMAGE_WALK_SEC_PIXELS
+#undef IMAGE_WALK_SEC_COLS
+#undef IMAGE_WALK_SEC_ROWS
+#undef IMAGE_WALK_SEC_COL_STEP
+#undef IMAGE_WALK_SEC_ROW_STEP
+#undef IMAGE_WALK_DO_START
+#undef IMAGE_WALK_DO_END
+#undef IMAGE_WALK_DO_ROW_START
+#undef IMAGE_WALK_DO_ROW_END
+#undef IMAGE_WALK_DO_STEP
+#undef STEP
diff --git a/images/image.c b/images/image.c
new file mode 100644 (file)
index 0000000..1db02b1
--- /dev/null
@@ -0,0 +1,211 @@
+/*
+ *     Image Library -- Basic image manipulation
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU Lesser General Public License.
+ */
+
+#undef LOCAL_DEBUG
+
+#include "lib/lib.h"
+#include "lib/mempool.h"
+#include "images/images.h"
+#include <string.h>
+
+#define MAX_IMAGE_BYTES (1 << 30)
+
+void
+image_thread_init(struct image_thread *it)
+{
+  DBG("image_thread_init()");
+  bzero(it, sizeof(*it));
+  it->pool = mp_new(1024);
+}
+
+void
+image_thread_cleanup(struct image_thread *it)
+{
+  DBG("image_thread_cleanup()");
+  mp_delete(it->pool);
+}
+
+void
+image_thread_err_format(struct image_thread *it, uns code, char *msg, ...)
+{
+  va_list args;
+  va_start(args, msg);
+  it->err_code = code;
+  it->err_msg = mp_vprintf(it->pool, msg, args);
+  va_end(args);
+}
+
+struct image *
+image_new(struct image_thread *it, uns cols, uns rows, uns flags, struct mempool *pool)
+{
+  DBG("image_new(cols=%u rows=%u flags=0x%x pool=%p)", cols, rows, flags, pool);
+  if (cols > IMAGE_MAX_SIZE || rows > IMAGE_MAX_SIZE)
+    {
+      image_thread_err(it, IMAGE_ERR_INVALID_DIMENSIONS, "Image dimension(s) too large");
+      return NULL;
+    }
+  struct image *img;
+  uns pixel_size, row_size, image_size, align;
+  switch (flags & IMAGE_COLOR_SPACE)
+    {
+      case COLOR_SPACE_GRAYSCALE:
+       pixel_size = 1;
+       break;
+      case COLOR_SPACE_RGB:
+       pixel_size = 3;
+       break;
+      default:
+       ASSERT(0);
+    }
+  if (flags & IMAGE_ALPHA)
+    pixel_size++;
+  switch (pixel_size)
+    {
+      case 1:
+      case 2:
+      case 4:
+       flags |= IMAGE_PIXELS_ALIGNED;
+       break;
+      case 3:
+       if (flags & IMAGE_PIXELS_ALIGNED)
+         pixel_size = 4;
+       break;
+      default:
+       ASSERT(0);
+    }
+  if (flags & IMAGE_SSE_ALIGNED)
+    align = IMAGE_SSE_ALIGN_SIZE;
+  else if (flags & IMAGE_PIXELS_ALIGNED)
+    align = pixel_size;
+  else
+    align = 1;
+  row_size = cols * pixel_size;
+  row_size = ALIGN(row_size, align);
+  u64 image_size_64 = (u64)row_size * rows;
+  u64 bytes_64 = image_size_64 + (sizeof(struct image) + IMAGE_SSE_ALIGN_SIZE - 1 + sizeof(uns));
+  if (bytes_64 > MAX_IMAGE_BYTES)
+    {
+      image_thread_err(it, IMAGE_ERR_INVALID_DIMENSIONS, "Image does not fit in memory");
+      return NULL;
+    }
+  if (!(image_size = image_size_64))
+    {
+      image_thread_err(it, IMAGE_ERR_INVALID_DIMENSIONS, "Zero dimension(s)");
+      return NULL;
+    }
+  img = pool ? mp_alloc(pool, (uns)bytes_64) : xmalloc((uns)bytes_64);
+  bzero(img, sizeof(struct image));
+  byte *p = (byte *)img + sizeof(struct image);
+  img->pixels = ALIGN_PTR(p, IMAGE_SSE_ALIGN_SIZE);
+  img->flags = flags;
+  img->pixel_size = pixel_size;
+  img->cols = cols;
+  img->rows = rows;
+  img->row_size = row_size;
+  img->image_size = image_size;
+  DBG("img=%p flags=0x%x pixel_size=%u row_size=%u image_size=%u pixels=%p",
+    img, img->flags, img->pixel_size, img->row_size, img->image_size, img->pixels);
+  return img;
+}
+
+struct image *
+image_clone(struct image_thread *it, struct image *src, uns flags, struct mempool *pool)
+{
+  DBG("image_clone(src=%p flags=0x%x pool=%p)", src, src->flags, pool);
+  struct image *img;
+  flags &= ~IMAGE_CHANNELS_FORMAT;
+  flags |= src->flags & IMAGE_CHANNELS_FORMAT;
+  if (!(img = image_new(it, src->cols, src->rows, flags, pool)))
+    return NULL;
+  if (img->image_size)
+    {
+      if (src->pixel_size != img->pixel_size)
+        {
+         struct image *sec_img = src;
+#         define IMAGE_WALK_INLINE
+#         define IMAGE_WALK_DOUBLE
+#         define IMAGE_WALK_DO_STEP do{ *(u32 *)pos = *(u32 *)sec_pos; }while(0)
+#         include "images/image-walk.h"
+       }
+      else if (src->row_size != img->row_size)
+        {
+          byte *s = src->pixels;
+          byte *d = img->pixels;
+          uns bytes = src->cols * img->pixel_size;
+         for (uns row = src->rows; row--; )
+            {
+             memcpy(d, s, bytes);
+             d += img->row_size;
+             s += src->row_size;
+           }
+        }
+      else
+        memcpy(img->pixels, src->pixels, img->image_size);
+    }
+  return img;
+}
+
+void
+image_destroy(struct image_thread *it UNUSED, struct image *img)
+{
+  DBG("image_destroy(img=%p)", img);
+  xfree((byte *)img);
+}
+
+void
+image_clear(struct image_thread *it UNUSED, struct image *img)
+{
+  DBG("image_clear(img=%p)", img);
+  if (img->image_size)
+    bzero(img->pixels, img->image_size);
+}
+
+byte *
+color_space_to_name(enum color_space cs)
+{
+  return image_channels_format_to_name(cs);
+}
+
+byte *
+image_channels_format_to_name(uns format)
+{
+  switch (format)
+    {
+      case COLOR_SPACE_GRAYSCALE:
+       return "Gray";
+      case COLOR_SPACE_GRAYSCALE | IMAGE_ALPHA:
+       return "GrayAlpha";
+      case COLOR_SPACE_RGB:
+       return "RGB";
+      case COLOR_SPACE_RGB | IMAGE_ALPHA:
+       return "RGBAlpha";
+      default:
+       return NULL;
+    }
+}
+
+uns
+image_name_to_channels_format(byte *name)
+{
+  if (!strcasecmp(name, "gray"))
+    return COLOR_SPACE_GRAYSCALE;
+  if (!strcasecmp(name, "grayscale"))
+    return COLOR_SPACE_GRAYSCALE;
+  if (!strcasecmp(name, "grayalpha"))
+    return COLOR_SPACE_GRAYSCALE | IMAGE_ALPHA;
+  if (!strcasecmp(name, "grayscalealpha"))
+    return COLOR_SPACE_GRAYSCALE | IMAGE_ALPHA;
+  if (!strcasecmp(name, "rgb"))
+    return COLOR_SPACE_RGB;
+  if (!strcasecmp(name, "rgbalpha"))
+    return COLOR_SPACE_RGB + IMAGE_ALPHA;
+  if (!strcasecmp(name, "rgba"))
+    return COLOR_SPACE_RGB + IMAGE_ALPHA;
+  return 0;
+}
diff --git a/images/images.h b/images/images.h
new file mode 100644 (file)
index 0000000..91c43d1
--- /dev/null
@@ -0,0 +1,168 @@
+#ifndef _IMAGES_IMAGES_H
+#define _IMAGES_IMAGES_H
+
+#include "lib/mempool.h"
+
+/* image.c */
+
+/* error handling */
+
+enum image_error {
+  IMAGE_ERR_OK = 0,
+  IMAGE_ERR_UNSPECIFIED,
+  IMAGE_ERR_NOT_IMPLEMENTED,
+  IMAGE_ERR_INVALID_DIMENSIONS,
+  IMAGE_ERR_INVALID_FILE_FORMAT,
+  IMAGE_ERR_INVALID_PIXEL_FORMAT,
+  IMAGE_ERR_READ_FAILED,
+  IMAGE_ERR_WRITE_FAILED,
+  IMAGE_ERR_MAX
+};
+
+struct image_thread {
+  byte *err_msg;
+  enum image_error err_code;
+  struct mempool *pool;
+};
+
+void image_thread_init(struct image_thread *thread);
+void image_thread_cleanup(struct image_thread *thread);
+
+static inline void
+image_thread_flush(struct image_thread *thread)
+{
+  thread->err_code = 0;
+  thread->err_msg = NULL;
+  mp_flush(thread->pool);
+}
+
+static inline void
+image_thread_err(struct image_thread *thread, uns code, char *msg)
+{
+  thread->err_code = code;
+  thread->err_msg = (byte *)msg;
+}
+
+void image_thread_err_format(struct image_thread *thread, uns code, char *msg, ...);
+
+/* basic image manupulation */
+
+#define IMAGE_MAX_SIZE         0xffffU /* maximum number of cols/rows, must be <(1<<16) */
+#define IMAGE_SSE_ALIGN_SIZE   (MAX(16, sizeof(uns)))
+
+enum color_space {
+  COLOR_SPACE_UNKNOWN,
+  COLOR_SPACE_GRAYSCALE,
+  COLOR_SPACE_RGB,
+  COLOR_SPACE_MAX
+};
+
+enum image_flag {
+  IMAGE_COLOR_SPACE = 0x7,     /* mask for enum color_space */
+  IMAGE_ALPHA = 0x8,           /* alpha channel */
+  IMAGE_PIXELS_ALIGNED = 0x10, /* align pixel size to the nearest power of two  */
+  IMAGE_SSE_ALIGNED = 0x20,    /* align scanlines to multiples of 16 bytes (both start and size) */
+  IMAGE_CHANNELS_FORMAT = IMAGE_COLOR_SPACE | IMAGE_ALPHA,
+  IMAGE_PIXEL_FORMAT = IMAGE_CHANNELS_FORMAT | IMAGE_PIXELS_ALIGNED,
+  IMAGE_ALIGNED = IMAGE_PIXELS_ALIGNED | IMAGE_SSE_ALIGNED,
+};
+
+struct image {
+  byte *pixels;                        /* left top pixel, there are at least sizeof(uns)
+                                  unsed bytes after the buffer (possible optimizations) */
+  u32 cols;                    /* number of columns */
+  u32 rows;                    /* number of rows */
+  u32 pixel_size;              /* size of pixel (1, 2, 3 or 4) */
+  u32 row_size;                        /* scanline size in bytes */
+  u32 image_size;              /* size of pixels buffer (rows * rows_size) */
+  u32 flags;                   /* enum image_flag */
+};
+
+struct image *image_new(struct image_thread *it, uns cols, uns rows, uns flags, struct mempool *pool);
+struct image *image_clone(struct image_thread *it, struct image *src, uns flags, struct mempool *pool);
+void image_destroy(struct image_thread *it, struct image *img); /* only with NULL mempool */
+void image_clear(struct image_thread *it, struct image *img);
+
+byte *color_space_to_name(enum color_space cs);
+byte *image_channels_format_to_name(uns format);
+uns image_name_to_channels_format(byte *name);
+
+/* scale.c */
+
+int image_scale(struct image_thread *thread, struct image *dest, struct image *src);
+void image_dimensions_fit_to_box(u32 *cols, u32 *rows, u32 max_cols, u32 max_rows, uns upsample);
+
+/* image-io.c */
+
+enum image_format {
+  IMAGE_FORMAT_UNDEFINED,
+  IMAGE_FORMAT_JPEG,
+  IMAGE_FORMAT_PNG,
+  IMAGE_FORMAT_GIF,
+  IMAGE_FORMAT_MAX
+};
+
+struct image_io {
+  struct image *image;
+  struct fastbuf *fastbuf;
+  enum image_format format;
+  struct mempool *pool;
+  u32 cols;
+  u32 rows;
+  u32 flags;
+  /* internals */
+  struct image_thread *thread;
+  struct mempool *internal_pool;
+  int image_destroy;
+  void *read_data;
+  void (*read_cancel)(struct image_io *io);
+  union {
+    struct {
+    } jpeg;
+    struct {
+    } png;
+    struct {
+    } gif;
+  };
+};
+
+void image_io_init(struct image_thread *it, struct image_io *io);
+void image_io_cleanup(struct image_io *io);
+void image_io_reset(struct image_io *io);
+
+int image_io_read_header(struct image_io *io);
+struct image *image_io_read_data(struct image_io *io, int ref);
+struct image *image_io_read(struct image_io *io, int ref);
+
+int image_io_write(struct image_io *io);
+
+byte *image_format_to_extension(enum image_format format);
+enum image_format image_extension_to_format(byte *extension);
+enum image_format image_file_name_to_format(byte *file_name);
+
+/* internals */
+
+#ifdef CONFIG_LIBJPEG
+int libjpeg_read_header(struct image_io *io);
+int libjpeg_read_data(struct image_io *io);
+int libjpeg_write(struct image_io *io);
+#endif
+
+#ifdef CONFIG_LIBPNG
+int libpng_read_header(struct image_io *io);
+int libpng_read_data(struct image_io *io);
+int libpng_write(struct image_io *io);
+#endif
+
+#ifdef CONFIG_LIBUNGIF
+int libungif_read_header(struct image_io *io);
+int libungif_read_data(struct image_io *io);
+#endif
+
+#ifdef CONFIG_LIBMAGICK
+int libmagick_read_header(struct image_io *io);
+int libmagick_read_data(struct image_io *io);
+int libmagick_write(struct image_io *io);
+#endif
+
+#endif
diff --git a/images/img.h b/images/img.h
new file mode 100644 (file)
index 0000000..e28775e
--- /dev/null
@@ -0,0 +1,44 @@
+#ifndef __IMG_H__
+#  define __IMG_H__
+
+#  include <stdio.h>
+#  include <magick/api.h>
+#  include "lib/config.h"
+
+struct ClassInfo{
+  uns l, u, v;                /* average Luv coefficients */
+  uns lh, hl, hh;             /* energies in Daubechies wavelet bands */
+  uns count;                  /*number of blocks in this class*/
+};
+
+struct
+BlockInfo{
+  uns l, u, v;                /* average Luv coefficients */
+  uns lh, hl, hh;             /* energies in Daubechies wavelet bands */
+  u8 cls_num;                 /* number of class for this block*/
+};
+
+struct
+PerturbHisInfo{
+  unsigned block_num : 24;    /*24 bits for number of picture's block should be enough*/
+  unsigned cls_num : 8;
+};
+
+struct
+DecomposeImageInfo{
+  uns max_cls_num;              /*self explaining*/
+  uns threshold;                /*stopping condition*/
+  uns diff_threshold;           /*stopping condition*/
+  uns max_cycles;               /*max number of loops of k_means clustering algorithm*/
+  uns init_decomp_num;          /*number of init decompositios */
+  uns bi_len;                   /*number of image blocks*/
+};
+
+struct BlockInfo*
+computeBlockInfo(PixelPacket*, uns, uns,  uns*);
+uns
+decomposeImage(struct DecomposeImageInfo* dii, struct BlockInfo *bi);
+#endif
+
+
+
diff --git a/images/io-libjpeg.c b/images/io-libjpeg.c
new file mode 100644 (file)
index 0000000..bcfa135
--- /dev/null
@@ -0,0 +1,454 @@
+/*
+ *     Image Library -- libjpeg
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU Lesser General Public License.
+ */
+
+#undef LOCAL_DEBUG
+
+#include "lib/lib.h"
+#include "lib/mempool.h"
+#include "lib/fastbuf.h"
+#include "images/images.h"
+
+#include <stdio.h>
+#include <sys/types.h>
+#include <jpeglib.h>
+#include <setjmp.h>
+
+struct libjpeg_err {
+  struct jpeg_error_mgr pub;
+  jmp_buf setjmp_buf;
+  struct image_io *io;
+};
+
+struct libjpeg_read_internals {
+  struct jpeg_decompress_struct cinfo;
+  struct jpeg_source_mgr src;
+  struct libjpeg_err err;
+  struct fastbuf *fastbuf;
+  byte *fastbuf_pos;
+};
+
+struct libjpeg_write_internals {
+  struct jpeg_compress_struct cinfo;
+  struct jpeg_destination_mgr dest;
+  struct libjpeg_err err;
+  struct fastbuf *fastbuf;
+  byte *fastbuf_pos;
+};
+
+static void NONRET
+libjpeg_read_error_exit(j_common_ptr cinfo)
+{
+  DBG("libjpeg_error_exit()");
+  struct libjpeg_err *e = (struct libjpeg_err *)cinfo->err;
+  byte buf[JMSG_LENGTH_MAX];
+  e->pub.format_message(cinfo, buf);
+  image_thread_err(e->io->thread, IMAGE_ERR_READ_FAILED, buf);
+  longjmp(e->setjmp_buf, 1);
+}
+
+static void NONRET
+libjpeg_write_error_exit(j_common_ptr cinfo)
+{
+  DBG("libjpeg_error_exit()");
+  struct libjpeg_err *e = (struct libjpeg_err *)cinfo->err;
+  byte buf[JMSG_LENGTH_MAX];
+  e->pub.format_message(cinfo, buf);
+  image_thread_err(e->io->thread, IMAGE_ERR_WRITE_FAILED,  buf);
+  longjmp(e->setjmp_buf, 1);
+}
+
+static void
+libjpeg_emit_message(j_common_ptr cinfo UNUSED, int msg_level UNUSED)
+{
+#ifdef LOCAL_DEBUG  
+  byte buf[JMSG_LENGTH_MAX];
+  cinfo->err->format_message(cinfo, buf);
+  DBG("libjpeg_emit_message(): [%d] %s", msg_level, buf);
+#endif  
+  if (unlikely(msg_level == -1))
+    longjmp(((struct libjpeg_err *)(cinfo)->err)->setjmp_buf, 1);
+}
+
+static inline uns
+libjpeg_fastbuf_read_prepare(struct libjpeg_read_internals *i)
+{
+  byte *start;
+  uns len = bdirect_read_prepare(i->fastbuf, &start);
+  i->fastbuf_pos = start + len;
+  i->src.next_input_byte = start;
+  i->src.bytes_in_buffer = len;
+  return len;
+}
+
+static inline void
+libjpeg_fastbuf_read_commit(struct libjpeg_read_internals *i)
+{
+  bdirect_read_commit(i->fastbuf, i->fastbuf_pos);
+}
+
+static void
+libjpeg_init_source(j_decompress_ptr cinfo)
+{
+  DBG("libjpeg_init_source()");
+  libjpeg_fastbuf_read_prepare((struct libjpeg_read_internals *)cinfo);
+}
+
+static void
+libjpeg_term_source(j_decompress_ptr cinfo UNUSED)
+{
+  DBG("libjpeg_term_source()");
+  //libjpeg_fastbuf_read_commit((struct libjpeg_read_internals *)cinfo);
+}
+
+static boolean
+libjpeg_fill_input_buffer(j_decompress_ptr cinfo)
+{
+  DBG("libjpeg_fill_input_buffer()");
+  struct libjpeg_read_internals *i = (struct libjpeg_read_internals *)cinfo;
+  libjpeg_fastbuf_read_commit(i);
+  return !!libjpeg_fastbuf_read_prepare(i);
+}
+
+static void
+libjpeg_skip_input_data(j_decompress_ptr cinfo, long num_bytes)
+{
+  DBG("libjpeg_skip_input_data(num_bytes=%d)", (int)num_bytes);
+  if (num_bytes > 0)
+    {
+      struct libjpeg_read_internals *i = (struct libjpeg_read_internals *)cinfo;
+      if ((unsigned long)num_bytes <= i->src.bytes_in_buffer)
+        {
+          i->src.next_input_byte += num_bytes;
+          i->src.bytes_in_buffer -= num_bytes;
+       }
+      else
+        {
+         num_bytes -= i->src.bytes_in_buffer;
+         libjpeg_fastbuf_read_commit(i);
+         bskip(i->fastbuf, num_bytes);
+         libjpeg_fastbuf_read_prepare(i);
+       }
+    }
+}
+
+static inline void
+libjpeg_fastbuf_write_prepare(struct libjpeg_write_internals *i)
+{
+  byte *start;
+  uns len = bdirect_write_prepare(i->fastbuf, &start);
+  i->fastbuf_pos = start + len;
+  i->dest.next_output_byte = start;
+  i->dest.free_in_buffer = len;
+  if (!len)
+    {
+      image_thread_err(i->err.io->thread, IMAGE_ERR_WRITE_FAILED, "Unexpected end of stream");
+      longjmp(i->err.setjmp_buf, 1);
+    }
+}
+
+static void
+libjpeg_init_destination(j_compress_ptr cinfo)
+{
+  DBG("libjpeg_init_destination()");
+  libjpeg_fastbuf_write_prepare((struct libjpeg_write_internals *)cinfo);
+}
+
+static void
+libjpeg_term_destination(j_compress_ptr cinfo)
+{
+  DBG("libjpeg_term_destination()");
+  struct libjpeg_write_internals *i = (struct libjpeg_write_internals *)cinfo;
+  bdirect_write_commit(i->fastbuf, (byte *)i->dest.next_output_byte);
+}
+
+static boolean
+libjpeg_empty_output_buffer(j_compress_ptr cinfo)
+{
+  DBG("libjpeg_empty_output_buffer()");
+  struct libjpeg_write_internals *i = (struct libjpeg_write_internals *)cinfo;
+  bdirect_write_commit(i->fastbuf, i->fastbuf_pos);
+  libjpeg_fastbuf_write_prepare(i);
+  return TRUE;
+}
+
+static void
+libjpeg_read_cancel(struct image_io *io)
+{
+  DBG("libjpeg_read_cancel()");
+  struct libjpeg_read_internals *i = io->read_data;
+  jpeg_destroy_decompress(&i->cinfo);
+}
+
+int
+libjpeg_read_header(struct image_io *io)
+{
+  DBG("libjpeg_read_header()");
+  struct libjpeg_read_internals *i = io->read_data = mp_alloc(io->internal_pool, sizeof(*i));
+  i->fastbuf = io->fastbuf;
+
+  /* Create libjpeg read structure */
+  DBG("Creating libjpeg read structure");
+  i->cinfo.err = jpeg_std_error(&i->err.pub);
+  i->err.pub.error_exit = libjpeg_read_error_exit;
+  i->err.pub.emit_message = libjpeg_emit_message;
+  i->err.io = io;
+  if (setjmp(i->err.setjmp_buf))
+    {
+      DBG("Libjpeg failed to read the image, longjump saved us");
+      jpeg_destroy_decompress(&i->cinfo);
+      return 0;
+    }
+  jpeg_create_decompress(&i->cinfo);
+
+  /* Initialize source manager */
+  i->cinfo.src = &i->src;
+  i->src.init_source = libjpeg_init_source;
+  i->src.fill_input_buffer = libjpeg_fill_input_buffer;
+  i->src.skip_input_data = libjpeg_skip_input_data;
+  i->src.resync_to_restart = jpeg_resync_to_restart;
+  i->src.term_source = libjpeg_term_source;
+
+  /* Read JPEG header and setup decompression options */
+  DBG("Reading image header");
+  jpeg_read_header(&i->cinfo, TRUE);
+  if (!(io->flags & IMAGE_COLOR_SPACE))
+    switch (i->cinfo.jpeg_color_space)
+      {
+        case JCS_GRAYSCALE:
+         io->flags |= COLOR_SPACE_GRAYSCALE;
+         break;
+        default:
+         io->flags |= COLOR_SPACE_RGB;
+         break;
+      }
+  if (!io->cols)
+    io->cols = i->cinfo.image_width;
+  if (!io->rows)
+    io->rows = i->cinfo.image_height;
+
+  io->read_cancel = libjpeg_read_cancel;
+  return 1;
+}
+
+int
+libjpeg_read_data(struct image_io *io)
+{
+  DBG("libjpeg_read_data()");
+
+  struct libjpeg_read_internals *i = io->read_data;
+  /* Select color space */ 
+  switch (io->flags & IMAGE_COLOR_SPACE)
+    {
+      case COLOR_SPACE_GRAYSCALE:
+       i->cinfo.out_color_space = JCS_GRAYSCALE;
+       break;
+      case COLOR_SPACE_RGB:
+       i->cinfo.out_color_space = JCS_RGB;
+       break;
+      default:
+       jpeg_destroy_decompress(&i->cinfo);
+       image_thread_err(io->thread, IMAGE_ERR_INVALID_PIXEL_FORMAT, "Unsupported color space.");
+       return 0;
+    }
+
+  /* Allocate the image... FIXME: use libjpeg feature to speed up downscale */
+  volatile int need_scale = io->cols != i->cinfo.image_width || io->rows != i->cinfo.image_height;
+  struct image * volatile img = need_scale ?
+    image_new(io->thread, i->cinfo.image_width, i->cinfo.image_height, io->flags & IMAGE_PIXEL_FORMAT, NULL) :
+    image_new(io->thread, i->cinfo.image_width, i->cinfo.image_height, io->flags, io->pool);
+  if (!img)
+    {
+      image_thread_err(io->thread, IMAGE_ERR_INVALID_PIXEL_FORMAT, "Unsupported color space.");
+      return 0;
+    }
+  
+  /* Setup fallback */
+  if (setjmp(i->err.setjmp_buf))
+    {
+      DBG("Libjpeg failed to read the image, longjump saved us");
+      jpeg_destroy_decompress(&i->cinfo);
+      if (need_scale || !io->pool)
+       image_destroy(io->thread, img);
+      return 0;
+    }
+
+  /* Decompress the image */
+  jpeg_start_decompress(&i->cinfo);
+  switch (img->pixel_size)
+    {
+      /* grayscale or RGB */
+      case 1:
+      case 3:
+       {
+          byte *pixels = img->pixels;
+         for (uns r = img->rows; r--; )
+            {
+              jpeg_read_scanlines(&i->cinfo, (JSAMPLE **)&pixels, 1);
+              pixels += img->row_size;
+            }
+       }
+       break;
+      /* garscale with alpha */
+      case 2:
+       {
+         byte buf[img->cols], *src;
+#         define IMAGE_WALK_INLINE
+#         define IMAGE_WALK_UNROLL 4
+#         define IMAGE_WALK_COL_STEP 2
+#         define IMAGE_WALK_DO_ROW_START do{ src = buf; jpeg_read_scanlines(&i->cinfo, (JSAMPLE **)&src, 1); }while(0)
+#         define IMAGE_WALK_DO_STEP do{ pos[0] = *src++; pos[1] = 255; }while(0)
+#         include "images/image-walk.h"
+       }
+       break;
+      /* RGBA or aligned RGB */
+      case 4:
+       {
+         byte buf[img->cols * 3], *src;
+#         define IMAGE_WALK_INLINE
+#         define IMAGE_WALK_UNROLL 4
+#         define IMAGE_WALK_COL_STEP 4
+#         define IMAGE_WALK_DO_ROW_START do{ src = buf; jpeg_read_scanlines(&i->cinfo, (JSAMPLE **)&src, 1); }while(0)
+#         define IMAGE_WALK_DO_STEP do{ *(u32 *)pos = *(u32 *)src; pos[3] = 255; src += 3; }while(0)
+#         include "images/image-walk.h"
+       }
+       break;
+      default:
+       ASSERT(0);
+    }
+  ASSERT(i->cinfo.output_scanline == i->cinfo.output_height);
+  
+  /* Destroy libjpeg object */
+  jpeg_finish_decompress(&i->cinfo);
+  jpeg_destroy_decompress(&i->cinfo);
+
+  /* Scale result if necessary */
+  if (need_scale)
+    {
+      DBG("Scaling image");
+      struct image *dest = image_new(io->thread, io->cols, io->rows, io->flags, io->pool);
+      if (!dest)
+        {
+         image_destroy(io->thread, img);
+         return 0;
+       }
+      if (!(image_scale(io->thread, dest, img)))
+        {
+         image_destroy(io->thread, img);
+         if (!io->pool)
+           image_destroy(io->thread, dest);
+         return 0;
+       }
+      image_destroy(io->thread, img);
+      io->image = dest;
+    }
+  else
+    io->image = img;
+  io->image_destroy = !io->pool;
+
+  DBG("Image readed");
+  return 1;
+}
+
+int
+libjpeg_write(struct image_io *io)
+{
+  DBG("libjpeg_write()");
+  struct libjpeg_write_internals i;
+  i.fastbuf = io->fastbuf;
+
+  /* Create libjpeg write structure */
+  DBG("Creating libjpeg write structure");
+  i.cinfo.err = jpeg_std_error(&i.err.pub);
+  i.err.pub.error_exit = libjpeg_write_error_exit;
+  i.err.pub.emit_message = libjpeg_emit_message;
+  i.err.io = io;
+  if (setjmp(i.err.setjmp_buf))
+    {
+      DBG("Libjpeg failed to write the image, longjump saved us");
+      jpeg_destroy_compress(&i.cinfo);
+      return 0;
+    }
+  jpeg_create_compress(&i.cinfo);
+  
+  /* Initialize destination manager */
+  i.cinfo.dest = &i.dest;
+  i.dest.init_destination = libjpeg_init_destination;
+  i.dest.term_destination = libjpeg_term_destination;
+  i.dest.empty_output_buffer = libjpeg_empty_output_buffer;
+
+  /* Set output parameters */
+  struct image *img = io->image;
+  i.cinfo.image_width = img->cols;
+  i.cinfo.image_height = img->rows;
+  switch (io->flags & IMAGE_COLOR_SPACE)
+    {
+      case COLOR_SPACE_GRAYSCALE:
+       i.cinfo.input_components = 1;
+       i.cinfo.in_color_space = JCS_GRAYSCALE;
+       break;
+      case COLOR_SPACE_RGB:
+       i.cinfo.input_components = 3;
+       i.cinfo.in_color_space = JCS_RGB;
+       break;
+      default:
+       jpeg_destroy_compress(&i.cinfo);
+       image_thread_err(io->thread, IMAGE_ERR_INVALID_PIXEL_FORMAT, "Unsupported pixel format.");
+       return 0;
+    }
+  jpeg_set_defaults(&i.cinfo);
+
+  /* Compress the image */
+  jpeg_start_compress(&i.cinfo, TRUE);
+  switch (img->pixel_size)
+    {
+      /* grayscale or RGB */
+      case 1:
+      case 3:
+       {
+          byte *pixels = img->pixels;
+         for (uns r = img->rows; r--; )
+           {
+              jpeg_write_scanlines(&i.cinfo, (JSAMPLE **)&pixels, 1);
+              pixels += img->row_size;
+            }
+       }
+       break;
+      /* grayscale with alpha (ignore alpha) */
+      case 2:
+       {
+         byte buf[img->cols], *dest = buf;
+#         define IMAGE_WALK_INLINE
+#         define IMAGE_WALK_UNROLL 4
+#         define IMAGE_WALK_COL_STEP 2
+#         define IMAGE_WALK_DO_ROW_END do{ dest = buf; jpeg_write_scanlines(&i.cinfo, (JSAMPLE **)&dest, 1); }while(0)
+#         define IMAGE_WALK_DO_STEP do{ *dest++ = pos[0]; }while(0)
+#         include "images/image-walk.h"
+       }
+       break;
+      /* RGBA (ignore alpha) or aligned RGB */
+      case 4:
+       {
+         byte buf[img->cols * 3], *dest = buf;
+#         define IMAGE_WALK_INLINE
+#         define IMAGE_WALK_UNROLL 4
+#         define IMAGE_WALK_COL_STEP 4
+#         define IMAGE_WALK_DO_ROW_END do{ dest = buf; jpeg_write_scanlines(&i.cinfo, (JSAMPLE **)&dest, 1); }while(0)
+#         define IMAGE_WALK_DO_STEP do{ *dest++ = pos[0]; *dest++ = pos[1]; *dest++ = pos[2]; }while(0)
+#         include "images/image-walk.h"
+       }
+       break;
+      default:
+       ASSERT(0);
+    }
+  ASSERT(i.cinfo.next_scanline == i.cinfo.image_height);
+  jpeg_finish_compress(&i.cinfo);
+  jpeg_destroy_compress(&i.cinfo);
+  return 1;
+}
diff --git a/images/io-libmagick.c b/images/io-libmagick.c
new file mode 100644 (file)
index 0000000..1bbdfc0
--- /dev/null
@@ -0,0 +1,391 @@
+/*
+ *     Image Library -- GraphicsMagick (slow fallback library)
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU Lesser General Public License.
+ */
+
+#define LOCAL_DEBUG
+
+#include "lib/lib.h"
+#include "lib/mempool.h"
+#include "lib/fastbuf.h"
+#include "images/images.h"
+#include <sys/types.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <magick/api.h>
+
+#define MAX_FILE_SIZE (1 << 30)
+#define QUANTUM_SCALE (QuantumDepth - 8)
+#define QUANTUM_TO_BYTE(x) ((uns)(x) >> QUANTUM_SCALE)
+#define BYTE_TO_QUANTUM(x) ((uns)(x) << QUANTUM_SCALE)
+#define OPACITY_MAX ((1 << QuantumDepth) - 1)
+
+struct magick_read_data {
+  ExceptionInfo exception;
+  ImageInfo *info;
+  Image *image;
+};
+
+static inline void
+libmagick_destroy_read_data(struct magick_read_data *rd)
+{
+  if (rd->image)
+    DestroyImage(rd->image);
+  DestroyImageInfo(rd->info);
+  DestroyExceptionInfo(&rd->exception);
+  DestroyMagick();
+}
+
+static void
+libmagick_read_cancel(struct image_io *io)
+{
+  DBG("libmagick_read_cancel()");
+
+  struct magick_read_data *rd = io->read_data;
+
+  DestroyImage(rd->image);
+  libmagick_destroy_read_data(rd);
+}
+
+int
+libmagick_read_header(struct image_io *io)
+{
+  DBG("libmagick_read_header()");
+
+  /* Read entire stream */
+  sh_off_t file_size = bfilesize(io->fastbuf) - btell(io->fastbuf);
+  if (unlikely(file_size > MAX_FILE_SIZE))
+    {
+      image_thread_err(io->thread, IMAGE_ERR_READ_FAILED, "Too long stream.");
+      return 0;
+    }
+  uns buf_size = file_size;
+  byte *buf = xmalloc(buf_size);
+  bread(io->fastbuf, buf, buf_size);
+
+  /* Allocate read structure */
+  struct magick_read_data *rd = io->read_data = mp_alloc(io->internal_pool, sizeof(*rd));
+
+  /* Initialize GraphicsMagick */
+  InitializeMagick(NULL);
+  GetExceptionInfo(&rd->exception);
+  rd->info = CloneImageInfo(NULL);
+  rd->info->subrange = 1;
+
+  /* Read the image */
+  rd->image = BlobToImage(rd->info, buf, buf_size, &rd->exception);
+  xfree(buf);
+  if (unlikely(!rd->image))
+    {
+      image_thread_err(io->thread, IMAGE_ERR_READ_FAILED, "GraphicsMagick failed to read the image.");
+      goto err;
+    }
+  if (unlikely(rd->image->columns > IMAGE_MAX_SIZE || rd->image->rows > IMAGE_MAX_SIZE))
+    {
+      image_thread_err(io->thread, IMAGE_ERR_INVALID_DIMENSIONS, "Image too large.");
+      goto err;
+    }
+
+  /* Fill image parameters */
+  if (!io->cols)
+    io->cols = rd->image->columns;
+  if (!io->rows)
+    io->rows = rd->image->rows;
+  if (!(io->flags & IMAGE_CHANNELS_FORMAT))
+    {
+      switch (rd->image->colorspace)
+        {
+         case GRAYColorspace:
+           io->flags |= COLOR_SPACE_GRAYSCALE | IMAGE_ALPHA;
+           break;
+         default:
+           io->flags |= COLOR_SPACE_RGB | IMAGE_ALPHA;
+           break;
+       }
+    }
+
+  io->read_cancel = libmagick_read_cancel;
+  return 1;
+
+err:
+  libmagick_destroy_read_data(rd);
+  return 0;
+}
+
+static inline byte
+libmagick_pixel_to_gray(PixelPacket *pixel)
+{
+  return ((uns)pixel->red * 19660 + (uns)pixel->green * 38666 + (uns)pixel->blue * 7210) >> (16 + QUANTUM_SCALE);
+}
+
+int
+libmagick_read_data(struct image_io *io)
+{
+  DBG("libmagick_read_data()");
+
+  struct magick_read_data *rd = io->read_data;
+
+  /* Quantize image */
+  switch (rd->image->colorspace)
+    {
+      case RGBColorspace:
+      case GRAYColorspace:
+        break;
+      default: ;
+        QuantizeInfo quantize;
+        GetQuantizeInfo(&quantize);
+        quantize.colorspace = RGBColorspace;
+        QuantizeImage(&quantize, rd->image);
+       break;
+    }
+
+  /* Allocate image for conversion */
+  int need_scale = io->cols != rd->image->columns || io->rows != rd->image->rows;
+  int need_destroy = need_scale || !io->pool;
+  struct image *img = need_scale ?
+    image_new(io->thread, rd->image->columns, rd->image->rows, io->flags & IMAGE_CHANNELS_FORMAT, NULL) :
+    image_new(io->thread, io->cols, io->rows, io->flags, io->pool);
+  if (unlikely(!img))
+    goto err;
+
+  /* Acquire pixels */
+  PixelPacket *src = (PixelPacket *)AcquireImagePixels(rd->image, 0, 0, rd->image->columns, rd->image->rows, &rd->exception);
+  if (unlikely(!src))
+    {
+      image_thread_err(io->thread, IMAGE_ERR_READ_FAILED, "Cannot acquire image pixels.");
+      goto err;
+    }
+
+  /* Convert pixels */
+  switch (img->pixel_size)
+    {
+      case 1:
+#       define IMAGE_WALK_INLINE
+#       define IMAGE_WALK_UNROLL 4
+#       define IMAGE_WALK_COL_STEP 1
+#       define IMAGE_WALK_DO_STEP do{ \
+         pos[0] = libmagick_pixel_to_gray(src); \
+         src++; }while(0)
+#       include "images/image-walk.h"
+       break;
+
+      case 2:
+#       define IMAGE_WALK_INLINE
+#       define IMAGE_WALK_UNROLL 4
+#       define IMAGE_WALK_COL_STEP 2
+#       define IMAGE_WALK_DO_STEP do{ \
+         pos[0] = libmagick_pixel_to_gray(src); \
+         pos[1] = QUANTUM_TO_BYTE(src->opacity); \
+         src++; }while(0)
+#       include "images/image-walk.h"
+       break;
+
+      case 3:
+#       define IMAGE_WALK_INLINE
+#       define IMAGE_WALK_UNROLL 4
+#       define IMAGE_WALK_COL_STEP 3
+#       define IMAGE_WALK_DO_STEP do{ \
+         pos[0] = QUANTUM_TO_BYTE(src->red); \
+         pos[1] = QUANTUM_TO_BYTE(src->green); \
+         pos[2] = QUANTUM_TO_BYTE(src->blue); \
+         src++; }while(0)
+#       include "images/image-walk.h"
+       break;
+
+      case 4:
+#       define IMAGE_WALK_INLINE
+#       define IMAGE_WALK_UNROLL 4
+#       define IMAGE_WALK_COL_STEP 4
+#       define IMAGE_WALK_DO_STEP do{ \
+         pos[0] = QUANTUM_TO_BYTE(src->red); \
+         pos[1] = QUANTUM_TO_BYTE(src->green); \
+         pos[2] = QUANTUM_TO_BYTE(src->blue); \
+         pos[3] = QUANTUM_TO_BYTE(src->opacity); \
+         src++; }while(0)
+#       include "images/image-walk.h"
+       break;
+
+      default:
+       ASSERT(0);
+    }
+
+  /* Free GraphicsMagick structures */
+  libmagick_destroy_read_data(rd);
+
+  /* Scale image */
+  if (need_scale)
+    {
+      struct image *img2 = image_new(io->thread, io->cols, io->rows, io->flags, io->pool);
+      if (unlikely(!img2))
+        goto err2;
+      int result = image_scale(io->thread, img2, img);
+      image_destroy(io->thread, img);
+      img = img2;
+      need_destroy = !io->pool;
+      if (unlikely(!result))
+       goto err2;
+    }
+
+  /* Success */
+  io->image = img;
+  io->image_destroy = need_destroy;
+  return 1;
+
+  /* Free structures */
+err:
+  libmagick_destroy_read_data(rd);
+err2:
+  if (need_destroy)
+    image_destroy(io->thread, img);
+  return 0;
+}
+
+int
+libmagick_write(struct image_io *io)
+{
+  DBG("libmagick_write()");
+
+  /* Initialize GraphicsMagick */
+  int result = 0;
+  ExceptionInfo exception;
+  ImageInfo *info;
+  InitializeMagick(NULL);
+  GetExceptionInfo(&exception);
+  info = CloneImageInfo(NULL);
+
+  /* Setup image parameters and allocate the image*/
+  switch (io->flags & IMAGE_COLOR_SPACE)
+    {
+      case COLOR_SPACE_GRAYSCALE:
+       info->colorspace = GRAYColorspace;
+       break;
+      case COLOR_SPACE_RGB:
+        info->colorspace = RGBColorspace;
+        break;
+      default:
+        ASSERT(0);
+    }
+  switch (io->format)
+    {
+      case IMAGE_FORMAT_JPEG:
+       strcpy(info->magick, "JPEG");
+       break;
+      case IMAGE_FORMAT_PNG:
+       strcpy(info->magick, "PNG");
+       break;
+      case IMAGE_FORMAT_GIF:
+        strcpy(info->magick, "GIF");
+       break;
+      default:
+        ASSERT(0);
+    }
+  Image *image = AllocateImage(info);
+  if (unlikely(!image))
+    {
+      image_thread_err(io->thread, IMAGE_ERR_WRITE_FAILED, "GraphicsMagick failed to allocate the image.");
+      goto err;
+    }
+  image->columns = io->cols;
+  image->rows = io->rows;
+
+  /* Get pixels */
+  PixelPacket *pixels = SetImagePixels(image, 0, 0, io->cols, io->rows), *dest = pixels;
+  if (unlikely(!pixels))
+    {
+      image_thread_err(io->thread, IMAGE_ERR_WRITE_FAILED, "Cannot get GraphicsMagick pixels.");
+      goto err2;
+    }
+
+  /* Convert pixels */
+  struct image *img = io->image;
+  switch (img->pixel_size)
+    {
+      case 1:
+#       define IMAGE_WALK_INLINE
+#       define IMAGE_WALK_UNROLL 4
+#       define IMAGE_WALK_COL_STEP 1
+#       define IMAGE_WALK_DO_STEP do{ \
+         dest->red = BYTE_TO_QUANTUM(pos[0]); \
+         dest->green = BYTE_TO_QUANTUM(pos[0]); \
+         dest->blue = BYTE_TO_QUANTUM(pos[0]); \
+         dest->opacity = OPACITY_MAX; \
+         dest++; }while(0)
+#       include "images/image-walk.h"
+       break;
+      case 2:
+#       define IMAGE_WALK_INLINE
+#       define IMAGE_WALK_UNROLL 4
+#       define IMAGE_WALK_COL_STEP 2
+#       define IMAGE_WALK_DO_STEP do{ \
+         dest->red = BYTE_TO_QUANTUM(pos[0]); \
+         dest->green = BYTE_TO_QUANTUM(pos[0]); \
+         dest->blue = BYTE_TO_QUANTUM(pos[0]); \
+         dest->opacity = BYTE_TO_QUANTUM(pos[1]); \
+         dest++; }while(0)
+#       include "images/image-walk.h"
+       break;
+      case 3:
+#       define IMAGE_WALK_INLINE
+#       define IMAGE_WALK_UNROLL 4
+#       define IMAGE_WALK_COL_STEP 3
+#       define IMAGE_WALK_DO_STEP do{ \
+         dest->red = BYTE_TO_QUANTUM(pos[0]); \
+         dest->green = BYTE_TO_QUANTUM(pos[1]); \
+         dest->blue = BYTE_TO_QUANTUM(pos[2]); \
+         dest->opacity = OPACITY_MAX; \
+         dest++; }while(0)
+#       include "images/image-walk.h"
+       break;
+      case 4:
+#       define IMAGE_WALK_INLINE
+#       define IMAGE_WALK_UNROLL 4
+#       define IMAGE_WALK_COL_STEP 4
+#       define IMAGE_WALK_DO_STEP do{ \
+         dest->red = BYTE_TO_QUANTUM(pos[0]); \
+         dest->green = BYTE_TO_QUANTUM(pos[1]); \
+         dest->blue = BYTE_TO_QUANTUM(pos[2]); \
+         dest->opacity = BYTE_TO_QUANTUM(pos[3]); \
+         dest++; }while(0)
+#       include "images/image-walk.h"
+       break;
+    }
+
+  /* Store pixels */
+  if (unlikely(!SyncImagePixels(image)))
+    {
+      image_thread_err(io->thread, IMAGE_ERR_WRITE_FAILED, "Cannot sync GraphicsMagick pixels.");
+      goto err2;
+    }
+
+  /* Write image */
+  size_t buf_len = 0;
+  void *buf = ImageToBlob(info, image, &buf_len, &exception);
+  if (unlikely(!buf))
+    {
+      image_thread_err(io->thread, IMAGE_ERR_WRITE_FAILED, "GraphicsMagick failed to compress the image.");
+      goto err2;
+    }
+  if (unlikely(buf_len > MAX_FILE_SIZE))
+    {
+      image_thread_err(io->thread, IMAGE_ERR_WRITE_FAILED, "Image too large.");
+      goto err2;
+    }
+
+  /* Write to stream */
+  bwrite(io->fastbuf, buf, buf_len);
+
+  /* Success */
+  result = 1;
+
+err2:
+  DestroyImage(image);
+err:
+  DestroyImageInfo(info);
+  DestroyExceptionInfo(&exception);
+  DestroyMagick();
+  return result;
+}
diff --git a/images/io-libpng.c b/images/io-libpng.c
new file mode 100644 (file)
index 0000000..6ed0516
--- /dev/null
@@ -0,0 +1,268 @@
+/*
+ *     Image Library -- libpng
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU Lesser General Public License.
+ */
+
+#undef LOCAL_DEBUG
+
+#include "lib/lib.h"
+#include "lib/mempool.h"
+#include "lib/fastbuf.h"
+#include "images/images.h"
+#include <png.h>
+#include <setjmp.h>
+
+struct libpng_internals {
+  png_structp png_ptr;
+  png_infop info_ptr;
+  png_infop end_ptr;
+  png_uint_32 cols;
+  png_uint_32 rows;
+  int bit_depth;
+  int color_type;
+};
+
+static png_voidp
+libpng_malloc(png_structp png_ptr, png_size_t size)
+{
+  DBG("libpng_malloc(size=%u)", (uns)size);
+  return mp_alloc(png_get_mem_ptr(png_ptr), size);
+}
+
+static void
+libpng_free(png_structp png_ptr UNUSED, png_voidp ptr UNUSED)
+{
+  DBG("libpng_free()");
+}
+
+static void NONRET
+libpng_error(png_structp png_ptr, png_const_charp msg)
+{
+  DBG("libpng_error()");
+  image_thread_err(png_get_error_ptr(png_ptr), IMAGE_ERR_READ_FAILED, (byte *)msg);
+  longjmp(png_jmpbuf(png_ptr), 1);
+}
+
+static void
+libpng_warning(png_structp png_ptr UNUSED, png_const_charp msg UNUSED)
+{
+  DBG("libpng_warning(): %s", (byte *)msg);
+}
+
+static void
+libpng_read(png_structp png_ptr, png_bytep data, png_size_t length)
+{
+  DBG("libpng_read(): len=%d", (uns)length);
+  if (unlikely(bread(png_get_io_ptr(png_ptr), data, length) < length))
+    png_error(png_ptr, "Incomplete data");
+}
+
+static void
+libpng_read_cancel(struct image_io *io)
+{
+  DBG("libpng_read_cancel()");
+  struct libpng_internals *i = io->read_data;
+  png_destroy_read_struct(&i->png_ptr, &i->info_ptr, &i->end_ptr);
+}
+
+int
+libpng_read_header(struct image_io *io)
+{
+  DBG("libpng_read_header()");
+  struct libpng_internals *i = io->read_data = mp_alloc(io->internal_pool, sizeof(*i));
+  i->png_ptr = png_create_read_struct_2(PNG_LIBPNG_VER_STRING,
+      io->thread, libpng_error, libpng_warning,
+      io->internal_pool, libpng_malloc, libpng_free);
+  if (unlikely(!i->png_ptr))
+    goto err_create;
+  i->info_ptr = png_create_info_struct(i->png_ptr);
+  if (unlikely(!i->info_ptr))
+    {
+      png_destroy_read_struct(&i->png_ptr, NULL, NULL);
+      goto err_create;
+    }
+  i->end_ptr = png_create_info_struct(i->png_ptr);
+  if (unlikely(!i->end_ptr))
+    {
+      png_destroy_read_struct(&i->png_ptr, &i->info_ptr, NULL);
+      goto err_create;
+    }
+  if (setjmp(png_jmpbuf(i->png_ptr)))
+    {
+      DBG("Libpng failed to read the image, longjump saved us");
+      png_destroy_read_struct(&i->png_ptr, &i->info_ptr, &i->end_ptr);
+      return 0;
+    }
+  png_set_read_fn(i->png_ptr, io->fastbuf, libpng_read);
+  png_set_user_limits(i->png_ptr, 0xffff, 0xffff);
+
+  DBG("Reading image info");
+  png_read_info(i->png_ptr, i->info_ptr);
+  png_get_IHDR(i->png_ptr, i->info_ptr, &i->cols, &i->rows, &i->bit_depth, &i->color_type, NULL, NULL, NULL);
+
+  if (!io->cols)
+    io->cols = i->cols;
+  if (!io->rows)
+    io->rows = i->rows;
+  if (!(io->flags & IMAGE_CHANNELS_FORMAT))
+    switch (i->color_type)
+      {
+       case PNG_COLOR_TYPE_GRAY:
+         io->flags |= COLOR_SPACE_GRAYSCALE;
+         break;
+       case PNG_COLOR_TYPE_GRAY_ALPHA:
+         io->flags |= COLOR_SPACE_GRAYSCALE | IMAGE_ALPHA;
+         break;
+       case PNG_COLOR_TYPE_RGB:
+         io->flags |= COLOR_SPACE_RGB;
+         break;
+       case PNG_COLOR_TYPE_RGB_ALPHA:
+       case PNG_COLOR_TYPE_PALETTE:
+         io->flags |= COLOR_SPACE_RGB | IMAGE_ALPHA;
+         break;
+       default:
+         png_destroy_read_struct(&i->png_ptr, &i->info_ptr, &i->end_ptr);
+         image_thread_err(io->thread, IMAGE_ERR_READ_FAILED, "Unknown color type");
+         break;
+      }
+
+  io->read_cancel = libpng_read_cancel;
+  return 1;
+
+err_create:
+  image_thread_err(io->thread, IMAGE_ERR_READ_FAILED, "Cannot create libpng read structure.");
+  return 0;
+}
+
+int
+libpng_read_data(struct image_io *io)
+{
+  DBG("libpng_read_data()");
+
+  struct libpng_internals *i = io->read_data;
+
+  /* Test supported pixel formats */
+  switch (io->flags & IMAGE_COLOR_SPACE)
+    {
+      case COLOR_SPACE_GRAYSCALE:
+      case COLOR_SPACE_RGB:
+       break;
+      default:
+        png_destroy_read_struct(&i->png_ptr, &i->info_ptr, &i->end_ptr);
+       image_thread_err(io->thread, IMAGE_ERR_INVALID_PIXEL_FORMAT, "Unsupported color space.");
+        return 0;
+    }
+
+  volatile int need_scale = io->cols != i->cols || io->rows != i->rows;
+  struct image * volatile img = need_scale ? 
+    image_new(io->thread, i->cols, i->rows, io->flags & IMAGE_PIXEL_FORMAT, NULL) :
+    image_new(io->thread, i->cols, i->rows, io->flags, io->pool);
+  if (!img)
+    {
+      png_destroy_read_struct(&i->png_ptr, &i->info_ptr, &i->end_ptr);
+      return 0;
+    }
+
+  if (setjmp(png_jmpbuf(i->png_ptr)))
+    {
+      DBG("Libpng failed to read the image, longjump saved us");
+      png_destroy_read_struct(&i->png_ptr, &i->info_ptr, &i->end_ptr);
+      if (need_scale || !io->pool)
+       image_destroy(io->thread, img);
+      return 0;
+    }
+
+  /* Apply transformations */
+  if (i->bit_depth == 16)
+    png_set_strip_16(i->png_ptr);
+  switch (i->color_type)
+    {
+      case PNG_COLOR_TYPE_PALETTE:
+       if ((io->flags & IMAGE_COLOR_SPACE) == COLOR_SPACE_GRAYSCALE)
+         {
+           png_set_palette_to_rgb(i->png_ptr);
+           png_set_rgb_to_gray_fixed(i->png_ptr, 1, 21267, 71514);
+         }
+       else
+         png_set_palette_to_rgb(i->png_ptr);
+       if ((io->flags & IMAGE_ALPHA) || (io->flags & IMAGE_PIXEL_FORMAT) == (COLOR_SPACE_RGB | IMAGE_PIXELS_ALIGNED))
+          png_set_add_alpha(i->png_ptr, 255, PNG_FILLER_AFTER);
+       else
+         png_set_strip_alpha(i->png_ptr);
+       break;
+      case PNG_COLOR_TYPE_GRAY:
+       if ((io->flags & IMAGE_COLOR_SPACE) == COLOR_SPACE_RGB)
+          png_set_gray_to_rgb(i->png_ptr);
+       if (io->flags & IMAGE_ALPHA)
+         png_set_add_alpha(i->png_ptr, 255, PNG_FILLER_AFTER);
+       break;
+      case PNG_COLOR_TYPE_GRAY_ALPHA:
+       if ((io->flags & IMAGE_COLOR_SPACE) == COLOR_SPACE_RGB)
+          png_set_gray_to_rgb(i->png_ptr);
+       if (!(io->flags & IMAGE_ALPHA))
+          png_set_strip_alpha(i->png_ptr);
+       break;
+      case PNG_COLOR_TYPE_RGB:
+       if ((io->flags & IMAGE_COLOR_SPACE) == COLOR_SPACE_GRAYSCALE)
+         png_set_rgb_to_gray_fixed(i->png_ptr, 1, 21267, 71514);
+       if ((io->flags & IMAGE_ALPHA) || (io->flags & IMAGE_PIXEL_FORMAT) == (COLOR_SPACE_RGB | IMAGE_PIXELS_ALIGNED))
+         png_set_add_alpha(i->png_ptr, 255, PNG_FILLER_AFTER);
+       break;
+      case PNG_COLOR_TYPE_RGB_ALPHA:
+       if ((io->flags & IMAGE_COLOR_SPACE) == COLOR_SPACE_GRAYSCALE)
+         png_set_rgb_to_gray_fixed(i->png_ptr, 1, 21267, 71514);
+       if (!(io->flags & IMAGE_ALPHA) && (io->flags & IMAGE_PIXEL_FORMAT) != (COLOR_SPACE_RGB | IMAGE_PIXELS_ALIGNED))
+          png_set_strip_alpha(i->png_ptr);
+       break;
+      default:
+       ASSERT(0);
+    }
+  png_read_update_info(i->png_ptr, i->info_ptr);
+
+  /* Read image data */
+  DBG("Reading image data");
+  byte *pixels = img->pixels;
+  png_bytep rows[img->rows];
+  for (uns r = 0; r < img->rows; r++, pixels += img->row_size)
+    rows[r] = (png_bytep)pixels;
+  png_read_image(i->png_ptr, rows);
+  png_read_end(i->png_ptr, i->end_ptr);
+
+  /* Destroy libpng read structure */
+  png_destroy_read_struct(&i->png_ptr, &i->info_ptr, &i->end_ptr);
+
+  /* Scale and store the resulting image */
+  if (need_scale)
+    {
+      struct image *dest = image_new(io->thread, io->cols, io->rows, io->flags, io->pool);
+      if (!dest)
+        {
+         image_destroy(io->thread, img);
+         return 0;
+       }
+      if (!image_scale(io->thread, dest, img))
+        {
+         image_destroy(io->thread, img);
+         if (!io->pool)
+           image_destroy(io->thread, dest);
+         return 0;
+       }
+      io->image = dest;
+    }
+  else
+    io->image = img;
+  io->image_destroy = !io->pool;
+  
+  return 1;
+}
+
+int
+libpng_write(struct image_io *io)
+{
+  image_thread_err(io->thread, IMAGE_ERR_NOT_IMPLEMENTED, "Libpng writing not implemented.");
+  return 0;
+}
diff --git a/images/io-libungif.c b/images/io-libungif.c
new file mode 100644 (file)
index 0000000..0cc2622
--- /dev/null
@@ -0,0 +1,26 @@
+/*
+ *     Image Library -- libungif
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU Lesser General Public License.
+ */
+
+#define LOCAL_DEBUG
+
+#include "lib/lib.h"
+#include "images/images.h"
+
+int
+libungif_read_header(struct image_io *io)
+{
+  image_thread_err(io->thread, IMAGE_ERR_NOT_IMPLEMENTED, "Libungif read not implemented.");
+  return 0;
+}
+
+int
+libungif_read_data(struct image_io *io UNUSED)
+{
+  ASSERT(0);
+}
diff --git a/images/io-main.c b/images/io-main.c
new file mode 100644 (file)
index 0000000..a930450
--- /dev/null
@@ -0,0 +1,234 @@
+/*
+ *     Image Library -- Image compression/decompression interface
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU Lesser General Public License.
+ */
+
+#undef LOCAL_DEBUG
+
+#include "lib/lib.h"
+#include "images/images.h"
+#include <string.h>
+
+void
+image_io_init(struct image_thread *it, struct image_io *io)
+{
+  DBG("image_io_init()");
+  bzero(io, sizeof(*io));
+  io->thread = it;
+  io->internal_pool = mp_new(1024);
+}
+
+static inline void
+image_io_read_cancel(struct image_io *io)
+{
+  if (io->read_cancel)
+    {
+      io->read_cancel(io);
+      io->read_cancel = NULL;
+    }
+}
+
+static inline void
+image_io_image_destroy(struct image_io *io)
+{
+  if (io->image_destroy)
+    {
+      image_destroy(io->thread, io->image);
+      io->image_destroy = 0;
+      io->image = NULL;
+    }
+}
+
+void
+image_io_cleanup(struct image_io *io)
+{
+  DBG("image_io_cleanup()");
+  image_io_read_cancel(io);
+  image_io_image_destroy(io);
+  mp_delete(io->internal_pool);
+}
+
+void
+image_io_reset(struct image_io *io)
+{
+  DBG("image_io_reset()");
+  image_io_read_cancel(io);
+  image_io_image_destroy(io);
+  struct mempool *pool = io->internal_pool;
+  mp_flush(pool);
+  bzero(io, sizeof(*io));
+  io->internal_pool = pool;
+}
+
+int
+image_io_read_header(struct image_io *io)
+{
+  DBG("image_io_read_header()");
+  image_io_read_cancel(io);
+  image_io_image_destroy(io);
+  switch (io->format) {
+    case IMAGE_FORMAT_JPEG:
+#if defined(CONFIG_LIBJPEG)
+      return libjpeg_read_header(io);
+#elif defined(CONFIG_LIBMAGICK)
+      return libmagick_read_header(io);
+#endif
+      break;
+
+    case IMAGE_FORMAT_PNG:
+#if defined(CONFIG_LIBPNG)
+      return libpng_read_header(io);
+#elif defined(CONFIG_LIBMAGICK)
+      return libmagick_read_header(io);
+#endif
+      break;
+
+    case IMAGE_FORMAT_GIF:
+#if defined(CONFIG_LIBUNGIG)
+      return libungif_read_header(io);
+#elif defined(CONFIG_LIBMAGICK)
+      return libmagick_read_header(io);
+#endif
+      break;
+
+    case IMAGE_FORMAT_UNDEFINED:
+      // FIXME: auto-detect
+      break;
+
+    default:
+      ASSERT(0);
+  }
+  image_thread_err(io->thread, IMAGE_ERR_INVALID_FILE_FORMAT, "Image format not supported.");
+  return 0;
+}
+
+struct image *
+image_io_read_data(struct image_io *io, int ref)
+{
+  DBG("image_io_read_data()");
+  ASSERT(io->read_cancel);
+  io->read_cancel = NULL;
+  int result;
+  switch (io->format) {
+    case IMAGE_FORMAT_JPEG:
+#if defined(CONFIG_LIBJPEG)
+      result = libjpeg_read_data(io);
+#elif defined(CONFIG_LIBMAGICK)
+      result = libmagick_read_data(io);
+#else
+      ASSERT(0);
+#endif
+      break;
+
+    case IMAGE_FORMAT_PNG:
+#if defined(CONFIG_LIBPNG)
+      result = libpng_read_data(io);
+#elif defined(CONFIG_LIBMAGICK)
+      result = libmagick_read_data(io);
+#else
+      ASSERT(0);
+#endif
+      break;
+
+    case IMAGE_FORMAT_GIF:
+#if defined(CONFIG_LIBMAGICK)
+      result = libmagick_read_data(io);
+#else
+      ASSERT(0);
+#endif
+      break;
+
+    default:
+      ASSERT(0);
+  }
+  if (result)
+    {
+      if (ref)
+       io->image_destroy = 0;
+      return io->image;
+    }
+  else
+    return NULL;
+}
+
+struct image *
+image_io_read(struct image_io *io, int ref)
+{
+  if (!image_io_read_header(io))
+    return NULL;
+  return image_io_read_data(io, ref);
+}
+
+int
+image_io_write(struct image_io *io)
+{
+  DBG("image_io_write()");
+  image_io_read_cancel(io);
+  switch (io->format) {
+    case IMAGE_FORMAT_JPEG:
+#if defined(CONFIG_LIBJPEG)
+      return libjpeg_write(io);
+#elif defined(CONFIG_LIBMAGICK)
+      return libmagick_write(io);
+#endif
+      break;
+
+    case IMAGE_FORMAT_PNG:
+#if defined(CONFIG_LIBMAGICK)
+      return libmagick_write(io);
+#endif
+      break;
+
+    case IMAGE_FORMAT_GIF:
+#if defined(CONFIG_LIBMAGICK)
+      return libmagick_write(io);
+#endif
+      break;
+
+    default:
+      break;
+  }
+  image_thread_err(io->thread, IMAGE_ERR_INVALID_FILE_FORMAT, "Image format not supported.");
+  return 0;
+}
+
+byte *
+image_format_to_extension(enum image_format format)
+{
+  switch (format)
+    {
+      case IMAGE_FORMAT_JPEG:
+       return "jpg";
+      case IMAGE_FORMAT_PNG:
+       return "png";
+      case IMAGE_FORMAT_GIF:
+       return "gif";
+      default:
+       return NULL;
+    }
+}
+
+enum image_format
+image_extension_to_format(byte *extension)
+{
+  if (!strcasecmp(extension, "jpg"))
+    return IMAGE_FORMAT_JPEG;
+  if (!strcasecmp(extension, "jpeg"))
+    return IMAGE_FORMAT_JPEG;
+  if (!strcasecmp(extension, "png"))
+    return IMAGE_FORMAT_PNG;
+  if (!strcasecmp(extension, "gif"))
+    return IMAGE_FORMAT_GIF;
+  return IMAGE_FORMAT_UNDEFINED;
+}
+
+enum image_format
+image_file_name_to_format(byte *file_name)
+{
+  byte *extension = strrchr(file_name, '.');
+  return extension ? image_extension_to_format(extension + 1) : IMAGE_FORMAT_UNDEFINED;
+}
diff --git a/images/k_means.c b/images/k_means.c
new file mode 100644 (file)
index 0000000..68a840e
--- /dev/null
@@ -0,0 +1,232 @@
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include <alloca.h>
+#include "lib/config.h"
+#include "img.h"
+
+#define stack_alloc alloca
+#define xdiff(x, y) ( ((x)>(y)) ? (x)-(y) : (y)-(x) )
+
+//void *memset(void *s, int c, size_t n);
+
+uns
+bi_dist(struct BlockInfo *bi, struct ClassInfo *ci){
+  return xdiff(bi->l, ci->l) +
+         xdiff(bi->u, ci->u) +
+         xdiff(bi->v, ci->v) +
+         xdiff(bi->lh, ci->lh) +
+         xdiff(bi->hl, ci->hl) +
+         xdiff(bi->hh, ci->hh);
+}
+
+
+/*"simulated annealing norm" - hodnocenní Ãºspì¹nosti simulovaného Â¾Ã­hání*/
+uns
+sa_norm(uns bi_len, struct BlockInfo *bi, u8 cls_num, struct ClassInfo *ci){
+  uns ret=0;
+  struct BlockInfo *pbi;
+  
+  for(pbi=bi; pbi < bi+bi_len; pbi++)
+    ret+=bi_dist(pbi, &ci[pbi->cls_num]);
+  return ret;
+}
+
+void
+showBlockInfoAsImage(struct BlockInfo *bi, uns bi_len, uns width, uns height, Image **__image,
+                     ExceptionInfo *exception){
+/*  ImageInfo *image_info=CloneImageInfo((ImageInfo *) NULL);*/
+/*  Image *image=AllocateImage(image_info);*/
+  unsigned char *p, *q;
+  uns x, y;
+  
+  width=(width>>2);
+  height=(height>>2);
+  assert(bi_len==width*height);
+  p=(unsigned char*) malloc(3*(width<<2)*(height<<2)*sizeof(unsigned char));
+  for (q=p, y=0; y < (height<<2); y++){
+    for (x=0; x < (width<<2); x++){
+      uns index=(y>>2)*width + (x>>2);
+      assert(index<bi_len);
+      *q++=(255/3) * (bi[index].cls_num>>2);
+      *q++=(255) * ((bi[index].cls_num>>1)&0x1);
+      *q++=(255) * ((bi[index].cls_num>>0)&0x1);
+    }
+  }
+  *__image=ConstituteImage(width<<2, height<<2, "RGB", CharPixel, p, exception);
+}
+
+void
+init_cls(uns cls_num, uns bi_len, u8 *cls){
+  u8 *p;
+  for(p=cls; p<cls+bi_len; p++)
+    *p = (u8) (((double)cls_num)*rand()/(RAND_MAX+1.0));
+}
+
+void
+eval_centers(uns bi_len, struct BlockInfo *bi, u8 cls_num, struct ClassInfo *ci){
+  struct BlockInfo *pbi;
+  struct ClassInfo *pci;
+
+  memset((void*) ci, 0, cls_num*sizeof(struct ClassInfo));
+  for(pbi=bi; pbi<bi+bi_len; pbi++){
+    struct ClassInfo *pcls=&ci[pbi->cls_num];
+    pcls->count++;
+    pcls->l+=pbi->l;
+    pcls->u+=pbi->u;
+    pcls->v+=pbi->v;
+    pcls->lh+=pbi->lh;
+    pcls->hl+=pbi->hl;
+    pcls->hh+=pbi->hh;
+  }
+
+  for(pci=ci; pci<ci+cls_num; pci++){
+    uns count=pci->count;
+    if(! count) continue;
+    pci->l  /=count;
+    pci->u  /=count;
+    pci->v  /=count;
+    pci->lh /=count;
+    pci->hl /=count;
+    pci->hh /=count;
+  }
+}
+
+#ifdef DEBUG_KMEANS
+void
+write_BlockInfo(uns len, struct BlockInfo *bi){
+  struct BlockInfo *pbi;
+  for(pbi=bi; pbi<bi+len; pbi++){
+    fprintf(stderr, "(%u, %u, %u, %u, %u, %u; '%u')\n", pbi->l, pbi->u, pbi->v, pbi->lh, pbi->hl, pbi->hh, pbi->cls_num);
+  }
+}
+
+void
+write_ClassInfo(uns len, struct ClassInfo *ci){
+  struct ClassInfo *pci;
+  for(pci=ci; pci<ci+len; pci++){
+    fprintf(stderr, "(%u, %u, %u, %u, %u, %u; %u)\n", pci->l, pci->u, pci->v, pci->lh, pci->hl, pci->hh, pci->count);
+  }
+}
+#endif
+
+/*bi does not mean an array, but pointer to block we search the nearest Class for*/
+u8
+nearestClass(uns cls_num, struct ClassInfo *ci, struct BlockInfo *bi){
+  u8 ret;
+  struct ClassInfo *pci;
+  uns min_dist=bi_dist(bi, ci);
+  ret=0;
+
+  for(pci=ci+1; pci<ci+cls_num; pci++){
+    uns dist=bi_dist(bi, pci);
+    if(dist<min_dist){
+      min_dist=dist;
+      ret = (u8) (pci-ci);
+    }
+  }
+  return ret;
+}
+
+uns
+__k_means(uns cls_num, struct DecomposeImageInfo* dii, struct BlockInfo* bi){
+  struct ClassInfo *ci=(struct ClassInfo*) stack_alloc(cls_num*sizeof(struct ClassInfo));
+  uns ret=~0U;
+  uns cycle=0;
+  eval_centers(dii->bi_len, bi, cls_num, ci);
+  for(cycle=0; cycle<dii->max_cycles; cycle++){
+    struct BlockInfo* pbi;
+    for(pbi=bi; pbi<bi+dii->bi_len; pbi++){
+      pbi->cls_num=nearestClass(cls_num, ci, pbi);
+    }
+    eval_centers(dii->bi_len, bi, cls_num, ci);
+    if((ret=sa_norm(dii->bi_len, bi, cls_num, ci))<dii->threshold)
+      break;
+  }
+  return ret;
+}
+
+void
+BlockInfoTou8(struct BlockInfo *bi, u8 *cls, uns len){
+  struct BlockInfo* pbi;
+  u8* pcls;
+
+  for(pbi=bi, pcls=cls; pbi<bi+len; pbi++)
+    *pcls++ = pbi->cls_num;
+}
+
+void
+u8ToBlockInfo(u8 *cls, struct BlockInfo *bi, uns len){
+  struct BlockInfo* pbi;
+  u8* pcls;
+
+  for(pbi=bi, pcls=cls; pbi<bi+len; pbi++)
+    pbi->cls_num = *pcls++;
+}
+
+uns
+k_means(uns cls_num, struct BlockInfo *bi, struct DecomposeImageInfo* dii){
+  u8 *act_cls=(u8*) stack_alloc(dii->bi_len*sizeof(u8));
+  u8 *best_cls=(u8*) stack_alloc(dii->bi_len*sizeof(u8));
+  uns best_diff=~0U, act_diff;
+  uns i;
+
+  for(i=0; i<dii->init_decomp_num; i++){
+    if(i)
+      init_cls(cls_num, dii->bi_len, act_cls);
+    else{
+    /*usually, the decompozition, that imply from here is not the best, but ensures, that the
+      return values from this fuction forms nonincreasing sequence for increasing cls_num*/
+      BlockInfoTou8(bi, act_cls, dii->bi_len);
+      bi[(uns) (dii->bi_len*rand()/(RAND_MAX+1.0))].cls_num=cls_num-1;
+    }
+      
+    
+    u8ToBlockInfo(act_cls, bi, dii->bi_len);
+    act_diff=__k_means(cls_num, dii, bi);
+    if(act_diff<best_diff){
+      best_diff=act_diff;
+      BlockInfoTou8(bi, best_cls, dii->bi_len);
+    }
+    /*fprintf(stderr, "...'%u'\n", act_diff);*/
+  }
+  u8ToBlockInfo(best_cls, bi, dii->bi_len);
+  return best_diff;
+}
+
+/*
+  return: final number of classes in decomposition
+*/
+uns
+decomposeImage(struct DecomposeImageInfo* dii, struct BlockInfo *bi){
+  uns cls_num=1;
+  uns err=0, old_err=0;
+
+  {
+    struct ClassInfo ci;
+    struct BlockInfo *pbi;
+    for(pbi=bi; pbi<bi+dii->bi_len; pbi++)
+      pbi->cls_num=0;
+    eval_centers(dii->bi_len, bi, 1, &ci);
+    old_err=sa_norm(dii->bi_len, bi, 1, &ci);
+  }
+  dii->threshold=old_err>>6;
+  dii->diff_threshold=old_err>>8;  
+  fprintf(stderr, "\n%u; --\n", old_err);
+  if(old_err<dii->threshold)
+    return old_err;
+  cls_num=1;
+  while(cls_num<dii->max_cls_num){
+    cls_num++;
+    err=k_means(cls_num, bi, dii);
+    fprintf(stderr, "\n(%u); %d\n", err, err-old_err);
+    if(err<dii->threshold) break;
+    if(err<old_err && err+dii->diff_threshold>old_err) break;
+    old_err=err;
+  }
+  return err;
+}
+
+
+
+
diff --git a/images/kd-tree.c b/images/kd-tree.c
new file mode 100644 (file)
index 0000000..f482854
--- /dev/null
@@ -0,0 +1,168 @@
+#undef LOCAL_DEBUG
+
+#include "sherlock/sherlock.h"
+#include "lib/heap.h"
+#include "images/images.h"
+#include "images/image-sig.h"
+#include "images/kd-tree.h"
+
+#include <alloca.h>
+
+#define SQR(x) ((x)*(x))
+#define IMAGE_SEARCH_CMP(x,y) (is->buf[x].dist < is->buf[y].dist)
+
+void
+image_search_init(struct image_search *is, struct image_tree *tree, struct image_vector *query, uns max_dist)
+{
+  // FIXME: empty tree
+  is->tree = tree;
+  is->nodes = tree->nodes;
+  is->leaves = tree->leaves;
+  is->query = *query;
+  is->max_dist = max_dist;
+  is->size = 0x1000;
+  is->buf = xmalloc((is->size + 1) * sizeof(struct image_search_item));
+  is->heap = xmalloc((is->size + 1) * sizeof(u32));
+  is->visited = is->count = 1;
+  is->heap[1] = 1;
+  struct image_search_item *item = is->buf + 1;
+  item->index = 1;
+  item->bbox = tree->bbox;
+  item->dist = 0;
+  for (uns i = 0; i < IMAGE_VEC_K; i++)
+    {
+      if (query->f[i] < item->bbox.vec[0].f[i])
+       item->dist += SQR(item->bbox.vec[0].f[i] - query->f[i]);
+      else if (query->f[i] > item->bbox.vec[1].f[i])
+       item->dist += SQR(query->f[i] - item->bbox.vec[0].f[i]);
+      else
+        {
+         item->dist = 0;
+         break;
+       }
+    }
+}
+
+void
+image_search_done(struct image_search *is)
+{
+  xfree(is->buf);
+  xfree(is->heap);
+}
+
+static void
+image_search_grow_slow(struct image_search *is)
+{
+  is->size *= 2;
+  is->buf = xrealloc(is->buf, (is->size + 1) * sizeof(struct image_search_item));
+  is->heap = xrealloc(is->heap, (is->size + 1) * sizeof(u32));
+}
+
+static inline struct image_search_item *
+image_search_grow(struct image_search *is)
+{
+  if (is->count == is->visited)
+  {
+    if (is->count == is->size)
+      image_search_grow_slow(is);
+    is->visited++;
+    is->heap[is->visited] = is->visited;
+  }
+  return is->buf + is->heap[++is->count];
+}
+
+static inline uns
+image_search_leaf_dist(struct image_search *is, struct image_bbox *bbox, struct image_leaf *leaf)
+{
+  uns dist = 0;
+  uns flags = leaf->flags; 
+  for (uns i = 0; i < IMAGE_VEC_K; i++)
+    {
+      uns bits = IMAGE_LEAF_BITS(i);
+      uns mask = (1 << bits) - 1;
+      uns value = flags & mask;
+      flags >>= bits;
+      int dif = bbox->vec[0].f[i] + (bbox->vec[1].f[i] - bbox->vec[0].f[i]) * value / ((1 << bits) - 1) - is->query.f[i];
+      dist += dif * dif;
+    }
+  return dist;
+}
+
+int
+image_search_next(struct image_search *is, oid_t *oid, uns *dist)
+{
+  while (likely(is->count))
+    {
+      struct image_search_item *item = is->buf + is->heap[1];
+      DBG("Main loop... dist=%d count=%d visited=%d size=%d index=0x%08x bbox=[(%s),(%s)]", 
+         item->dist, is->count, is->visited, is->size, item->index, 
+         stk_print_image_vector(&item->bbox.vec[0]), stk_print_image_vector(&item->bbox.vec[1]));
+      if (unlikely(item->dist > is->max_dist))
+        {
+         DBG("Maximum distance reached");
+         return 0;
+       }
+      
+      /* Expand leaf */
+      if (item->index & IMAGE_SEARCH_ITEM_TYPE)
+        {
+         *oid = item->index & ~IMAGE_SEARCH_ITEM_TYPE;
+         *dist = item->dist;
+         DBG("Found item %d at distance %d", *oid, *dist);
+         HEAP_DELMIN(u32, is->heap, is->count, IMAGE_SEARCH_CMP, HEAP_SWAP);
+         return 1;
+       }
+      
+      /* Expand node with leaves */
+      else if (is->nodes[item->index].val & IMAGE_NODE_LEAF)
+        {
+         DBG("Expanding node to list of leaves");
+         struct image_leaf *leaf = is->leaves + (is->nodes[item->index].val & ~IMAGE_NODE_LEAF);
+         item->dist = image_search_leaf_dist(is, &item->bbox, leaf);
+         item->index = IMAGE_SEARCH_ITEM_TYPE | leaf->oid;
+         HEAP_INCREASE(u32, is->heap, is->count, IMAGE_SEARCH_CMP, HEAP_SWAP, 1);
+         while (!((leaf++)->flags & IMAGE_LEAF_LAST))
+           {
+             struct image_search_item *nitem = image_search_grow(is);
+             nitem->dist = image_search_leaf_dist(is, &item->bbox, leaf);
+             nitem->index = IMAGE_SEARCH_ITEM_TYPE | leaf->oid;
+             HEAP_INSERT(u32, is->heap, is->count, IMAGE_SEARCH_CMP, HEAP_SWAP);
+           }
+       }
+
+      /* Expand internal node */
+      else
+        {
+         DBG("Expanding internal node");
+         struct image_search_item *nitem = image_search_grow(is);
+         uns dim = is->nodes[item->index].val & IMAGE_NODE_DIM;
+         uns pivot = is->nodes[item->index].val >> 8;
+         item->index *= 2;
+         nitem->bbox = item->bbox;
+         nitem->dist = item->dist;
+         uns query = is->query.f[dim];
+         int dif = query - pivot;
+         if (dif > 0)
+           {
+             nitem->index = item->index++;
+             item->bbox.vec[0].f[dim] = pivot;
+             nitem->bbox.vec[1].f[dim] = pivot;
+             if (query > item->bbox.vec[1].f[dim])
+               nitem->dist -= SQR(query - item->bbox.vec[1].f[dim]);
+           }
+         else
+           {
+             nitem->index = item->index + 1;
+             item->bbox.vec[1].f[dim] = pivot;
+             nitem->bbox.vec[0].f[dim] = pivot;
+             if (query < item->bbox.vec[0].f[dim])
+               nitem->dist -= SQR(item->bbox.vec[0].f[dim] - query);
+           }
+         nitem->dist += SQR(dif);
+         HEAP_INSERT(u32, is->heap, is->count, IMAGE_SEARCH_CMP, HEAP_SWAP);
+       }
+    }
+  DBG("Heap is empty");
+  return 0;
+}
+
diff --git a/images/kd-tree.h b/images/kd-tree.h
new file mode 100644 (file)
index 0000000..1cd1095
--- /dev/null
@@ -0,0 +1,28 @@
+#ifndef _IMAGES_KD_TREE_H
+#define _IMAGES_KD_TREE_H
+
+#define IMAGE_SEARCH_DIST_UNLIMITED    (~0U)
+
+/* FIXME: support full length of oid_t, currently must be <2^31 */
+#define IMAGE_SEARCH_ITEM_TYPE         0x80000000U
+struct image_search_item {
+  u32 dist;
+  u32 index;
+  struct image_bbox bbox;
+};
+
+struct image_search {
+  struct image_tree *tree;
+  struct image_node *nodes;
+  struct image_leaf *leaves;
+  struct image_vector query;
+  struct image_search_item *buf;
+  u32 *heap;
+  uns count, visited, size, max_dist;
+};
+
+void image_search_init(struct image_search *is, struct image_tree *tree, struct image_vector *query, uns max_dist);
+void image_search_done(struct image_search *is);
+int image_search_next(struct image_search *is, oid_t *oid, uns *dist);
+
+#endif
diff --git a/images/scale-gen.h b/images/scale-gen.h
new file mode 100644 (file)
index 0000000..ff51a51
--- /dev/null
@@ -0,0 +1,186 @@
+/*
+ *     Image Library -- Image scaling algorithms
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU Lesser General Public License.
+ */
+
+#ifndef IMAGE_SCALE_CHANNELS
+#  define IMAGE_SCALE_CHANNELS IMAGE_SCALE_PIXEL_SIZE
+#endif
+
+static void
+IMAGE_SCALE_PREFIX(downsample)(struct image *dest, struct image *src)
+{
+  /* FIXME slow */
+  byte *rsrc = src->pixels, *psrc;
+  byte *rdest = dest->pixels, *pdest;
+  uns x_inc = (dest->cols << 16) / src->cols, x_pos, x_inc_frac = 0xffffff / x_inc;
+  uns y_inc = (dest->rows << 16) / src->rows, y_pos = 0, y_inc_frac = 0xffffff / y_inc;
+  uns final_mul = ((u64)x_inc * y_inc) >> 16;
+  uns buf_size = dest->cols * IMAGE_SCALE_CHANNELS;
+  u32 buf[buf_size], *pbuf;
+  buf_size *= sizeof(u32);
+  bzero(buf, buf_size);
+  for (uns rows_counter = src->rows; rows_counter--; )
+    {
+      pbuf = buf;
+      psrc = rsrc;
+      rsrc += src->row_size;
+      x_pos = 0;
+      y_pos += y_inc;
+      if (y_pos <= 0x10000)
+        {
+          for (uns cols_counter = src->cols; cols_counter--; )
+            {
+             x_pos += x_inc;
+             if (x_pos <= 0x10000)
+               {
+                 pbuf[0] += psrc[0];
+#                if IMAGE_SCALE_CHANNELS >= 2
+                 pbuf[1] += psrc[1];
+#                 endif
+#                if IMAGE_SCALE_CHANNELS >= 3
+                 pbuf[2] += psrc[2];
+#                 endif
+#                if IMAGE_SCALE_CHANNELS >= 4
+                 pbuf[3] += psrc[3];
+#                 endif
+               }
+             else
+               {
+                 x_pos -= 0x10000;
+                 uns mul2 = x_pos * x_inc_frac;
+                 uns mul1 = 0xffffff - mul2;
+                 pbuf[0] += (psrc[0] * mul1) >> 24;
+                 pbuf[0 + IMAGE_SCALE_CHANNELS] += (psrc[0] * mul2) >> 24;
+#                if IMAGE_SCALE_CHANNELS >= 2
+                 pbuf[1] += (psrc[1] * mul1) >> 24;
+                 pbuf[1 + IMAGE_SCALE_CHANNELS] += (psrc[1] * mul2) >> 24;
+#                 endif
+#                if IMAGE_SCALE_CHANNELS >= 3
+                 pbuf[2] += (psrc[2] * mul1) >> 24;
+                 pbuf[2 + IMAGE_SCALE_CHANNELS] += (psrc[2] * mul2) >> 24;
+#                 endif
+#                if IMAGE_SCALE_CHANNELS >= 4
+                 pbuf[3] += (psrc[3] * mul1) >> 24;
+                 pbuf[3 + IMAGE_SCALE_CHANNELS] += (psrc[3] * mul2) >> 24;
+#                 endif
+                 pbuf += IMAGE_SCALE_CHANNELS;
+               }
+             psrc += IMAGE_SCALE_PIXEL_SIZE;
+           }
+       }
+      else
+        {
+         y_pos -= 0x10000;
+          pdest = rdest;
+          rdest += dest->row_size;
+         uns mul2 = y_pos * y_inc_frac;
+         uns mul1 = 0xffffff - mul2;
+         uns a0 = 0;
+#        if IMAGE_SCALE_CHANNELS >= 2
+         uns a1 = 0;
+#        endif
+#        if IMAGE_SCALE_CHANNELS >= 3
+         uns a2 = 0;
+#        endif
+#        if IMAGE_SCALE_CHANNELS >= 4
+         uns a3 = 0;
+#        endif
+          for (uns cols_counter = src->cols; cols_counter--; )
+            {
+             x_pos += x_inc;
+             if (x_pos <= 0x10000)
+               {
+                 pbuf[0] += ((psrc[0] * mul1) >> 24);
+                 a0 += (psrc[0] * mul2) >> 24;
+#                if IMAGE_SCALE_CHANNELS >= 2
+                 pbuf[1] += ((psrc[1] * mul1) >> 24);
+                 a1 += (psrc[1] * mul2) >> 24;
+#                 endif
+#                if IMAGE_SCALE_CHANNELS >= 3
+                 pbuf[2] += ((psrc[2] * mul1) >> 24);
+                 a2 += (psrc[2] * mul2) >> 24;
+#                 endif
+#                if IMAGE_SCALE_CHANNELS >= 4
+                 pbuf[3] += ((psrc[3] * mul1) >> 24);
+                 a3 += (psrc[3] * mul2) >> 24;
+#                 endif
+               }
+             else
+               {
+                 x_pos -= 0x10000;
+                 uns mul4 = x_pos * x_inc_frac;
+                 uns mul3 = 0xffffff - mul4;
+                 uns mul13 = ((u64)mul1 * mul3) >> 24;
+                 uns mul23 = ((u64)mul2 * mul3) >> 24;
+                 uns mul14 = ((u64)mul1 * mul4) >> 24;
+                 uns mul24 = ((u64)mul2 * mul4) >> 24;
+                 pdest[0] = ((((psrc[0] * mul13) >> 24) + pbuf[0]) * final_mul) >> 16;
+                 pbuf[0] = ((psrc[0] * mul23) >> 24) + a0;
+                 pbuf[0 + IMAGE_SCALE_CHANNELS] += ((psrc[0 + IMAGE_SCALE_PIXEL_SIZE] * mul14) >> 24);
+                 a0 = ((psrc[0 + IMAGE_SCALE_PIXEL_SIZE] * mul24) >> 24);
+#                if IMAGE_SCALE_CHANNELS >= 2
+                 pdest[1] = ((((psrc[1] * mul13) >> 24) + pbuf[1]) * final_mul) >> 16;
+                 pbuf[1] = ((psrc[1] * mul23) >> 24) + a1;
+                 pbuf[1 + IMAGE_SCALE_CHANNELS] += ((psrc[1 + IMAGE_SCALE_PIXEL_SIZE] * mul14) >> 24);
+                 a1 = ((psrc[1 + IMAGE_SCALE_PIXEL_SIZE] * mul24) >> 24);
+#                 endif
+#                if IMAGE_SCALE_CHANNELS >= 3
+                 pdest[2] = ((((psrc[2] * mul13) >> 24) + pbuf[2]) * final_mul) >> 16;
+                 pbuf[2] = ((psrc[2] * mul23) >> 24) + a2;
+                 pbuf[2 + IMAGE_SCALE_CHANNELS] += ((psrc[2 + IMAGE_SCALE_PIXEL_SIZE] * mul14) >> 24);
+                 a2 = ((psrc[2 + IMAGE_SCALE_PIXEL_SIZE] * mul24) >> 24);
+#                 endif
+#                if IMAGE_SCALE_CHANNELS >= 4
+                 pdest[3] = ((((psrc[3] * mul13) >> 24) + pbuf[3]) * final_mul) >> 16;
+                 pbuf[3] = ((psrc[3] * mul23) >> 24) + a3;
+                 pbuf[3 + IMAGE_SCALE_CHANNELS] += ((psrc[3 + IMAGE_SCALE_PIXEL_SIZE] * mul14) >> 24);
+                 a3 = ((psrc[3 + IMAGE_SCALE_PIXEL_SIZE] * mul24) >> 24);
+#                 endif
+                 pbuf += IMAGE_SCALE_CHANNELS;
+                 pdest += IMAGE_SCALE_PIXEL_SIZE;
+               }
+             psrc += IMAGE_SCALE_PIXEL_SIZE;
+           }
+         pdest[0] = (pbuf[0] * final_mul) >> 16;
+         pbuf[0] = a0;
+#         if IMAGE_SCALE_CHANNELS >= 2
+         pdest[1] = (pbuf[1] * final_mul) >> 16;
+         pbuf[1] = a1;
+#        endif
+#         if IMAGE_SCALE_CHANNELS >= 3
+         pdest[2] = (pbuf[2] * final_mul) >> 16;
+         pbuf[2] = a2;
+#        endif
+#         if IMAGE_SCALE_CHANNELS >= 4
+         pdest[3] = (pbuf[3] * final_mul) >> 16;
+         pbuf[3] = a3;
+#        endif
+       }
+    }
+  pdest = rdest;
+  pbuf = buf;
+  for (uns cols_counter = dest->cols; cols_counter--; )
+    {
+      pdest[0] = (pbuf[0] * final_mul) >> 16;
+#     if IMAGE_SCALE_CHANNELS >= 2
+      pdest[1] = (pbuf[1] * final_mul) >> 16;
+#     endif
+#     if IMAGE_SCALE_CHANNELS >= 3
+      pdest[2] = (pbuf[2] * final_mul) >> 16;
+#     endif
+#     if IMAGE_SCALE_CHANNELS >= 4
+      pdest[3] = (pbuf[3] * final_mul) >> 16;
+#     endif
+      pbuf += IMAGE_SCALE_CHANNELS;
+      pdest += IMAGE_SCALE_PIXEL_SIZE;
+    }
+}
+
+#undef IMAGE_SCALE_PREFIX
+#undef IMAGE_SCALE_PIXEL_SIZE
+#undef IMAGE_SCALE_CHANNELS
diff --git a/images/scale.c b/images/scale.c
new file mode 100644 (file)
index 0000000..4eceac6
--- /dev/null
@@ -0,0 +1,104 @@
+/*
+ *     Image Library -- Image scaling algorithms
+ *
+ *     (c) 2006 Pavel Charvat <pchar@ucw.cz>
+ *
+ *     This software may be freely distributed and used according to the terms
+ *     of the GNU Lesser General Public License.
+ */
+
+#undef LOCAL_DEBUG
+
+#include "lib/lib.h"
+#include "images/images.h"
+#include <string.h>
+
+#define IMAGE_SCALE_PREFIX(x) image_scale_1_##x
+#define IMAGE_SCALE_PIXEL_SIZE 1
+#include "images/scale-gen.h"
+
+#define IMAGE_SCALE_PREFIX(x) image_scale_2_##x
+#define IMAGE_SCALE_PIXEL_SIZE 2
+#include "images/scale-gen.h"
+
+#define IMAGE_SCALE_PREFIX(x) image_scale_3_##x
+#define IMAGE_SCALE_PIXEL_SIZE 3
+#include "images/scale-gen.h"
+
+#define IMAGE_SCALE_PREFIX(x) image_scale_4_##x
+#define IMAGE_SCALE_PIXEL_SIZE 4
+#include "images/scale-gen.h"
+
+int
+image_scale(struct image_thread *it, struct image *dest, struct image *src)
+{
+  if (src->cols < dest->cols || src->rows < dest->rows)
+    {
+      image_thread_err(it, IMAGE_ERR_INVALID_DIMENSIONS, "Upsampling not supported.");
+      return 0;
+    }
+  if ((src->flags & IMAGE_PIXEL_FORMAT) != (dest->flags & IMAGE_PIXEL_FORMAT))
+    {
+      image_thread_err(it, IMAGE_ERR_INVALID_PIXEL_FORMAT, "Different pixel format not supported.");
+      return 0;
+    }
+  switch (src->pixel_size)
+    {
+      /* Gray */
+      case 1:
+       image_scale_1_downsample(dest, src);
+       return 1;
+      /* GrayA */
+      case 2:
+       image_scale_2_downsample(dest, src);
+       return 1;
+      /* RGB */
+      case 3:
+       image_scale_3_downsample(dest, src);
+       return 1;
+      /* RGBA or aligned RGB */
+      case 4:
+       image_scale_4_downsample(dest, src);
+       return 1;
+      default:
+       ASSERT(0);
+    }
+}
+
+void
+image_dimensions_fit_to_box(u32 *cols, u32 *rows, u32 max_cols, u32 max_rows, uns upsample)
+{
+  ASSERT(*cols && *rows && *cols <= IMAGE_MAX_SIZE && *rows <= IMAGE_MAX_SIZE);
+  ASSERT(max_cols && max_rows && max_cols <= IMAGE_MAX_SIZE && max_rows <= IMAGE_MAX_SIZE);
+  if (*cols <= max_cols && *rows <= max_rows)
+    {
+      if (!upsample)
+       return;
+      if (max_cols / *cols > max_rows / *rows)
+        {
+         *cols = *cols * max_rows / *rows;
+         *cols = MIN(*cols, max_cols);
+         *rows = max_rows;
+       }
+      else
+        {
+         *rows = *rows * max_cols / *cols;
+         *rows = MIN(*rows, max_rows);
+         *cols = max_cols;
+       }
+    }
+  else if (*cols <= max_cols)
+    goto down_cols;
+  else if (*rows <= max_rows || max_rows * *cols > max_cols * *rows)
+    goto down_rows;
+down_cols:
+  *cols = *cols * max_rows / *rows;
+  *cols = MAX(*cols, 1);
+  *rows = max_rows;
+  return;
+down_rows:  
+  *rows = *rows * max_cols / *cols;
+  *rows = MAX(*rows, 1);
+  *cols = max_cols;
+  return;
+}