DIRS+=images
-PROGS+=$(addprefix $(o)/images/,image-idx image-test)
+PROGS+=$(addprefix $(o)/images/,image-idx image-test decomp)
$(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-idx: LIBS+=-lGraphicsMagick -ljpeg
$(o)/images/image-test: $(o)/images/images.o $(o)/images/image-test.o $(LIBSH)
+
+# By :;DF
+$(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
+
--- /dev/null
+/*#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;
+}
--- /dev/null
+#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;
+}
--- /dev/null
+#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
+
+
+
--- /dev/null
+#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;
+}
+
+
+
+
projects / libucw.git / commitdiff