From 59f66141d8aaa587c8274e005575f41b21e5159c Mon Sep 17 00:00:00 2001 From: Daniel Fiala Date: Sun, 16 Apr 2006 19:56:40 +0200 Subject: [PATCH] Adding files for image decompostition. --- images/Makefile | 11 ++- images/block_info.c | 146 ++++++++++++++++++++++++++++ images/decomp.c | 73 ++++++++++++++ images/img.h | 44 +++++++++ images/k_means.c | 232 ++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 505 insertions(+), 1 deletion(-) create mode 100644 images/block_info.c create mode 100644 images/decomp.c create mode 100644 images/img.h create mode 100644 images/k_means.c diff --git a/images/Makefile b/images/Makefile index de77b745..cd3af5ff 100644 --- a/images/Makefile +++ b/images/Makefile @@ -2,7 +2,7 @@ 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 @@ -10,3 +10,12 @@ $(o)/images/image-idx: $(o)/images/images.o $(o)/images/image-idx.o $(o)/images/ $(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 + diff --git a/images/block_info.c b/images/block_info.c new file mode 100644 index 00000000..a8de2bb3 --- /dev/null +++ b/images/block_info.c @@ -0,0 +1,146 @@ +/*#include +#include */ + +#include +#include +#include +#include +#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/decomp.c b/images/decomp.c new file mode 100644 index 00000000..1659efac --- /dev/null +++ b/images/decomp.c @@ -0,0 +1,73 @@ +#include +#include +#include +#include +#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/img.h b/images/img.h new file mode 100644 index 00000000..e28775e8 --- /dev/null +++ b/images/img.h @@ -0,0 +1,44 @@ +#ifndef __IMG_H__ +# define __IMG_H__ + +# include +# include +# 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/k_means.c b/images/k_means.c new file mode 100644 index 00000000..68a840e1 --- /dev/null +++ b/images/k_means.c @@ -0,0 +1,232 @@ +#include +#include +#include +#include +#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>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; pcls_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; pcicount; + 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; pbil, 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; pcil, 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; pcibi_len, bi, cls_num, ci); + for(cycle=0; cyclemax_cycles; cycle++){ + struct BlockInfo* pbi; + for(pbi=bi; pbibi_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))threshold) + break; + } + return ret; +} + +void +BlockInfoTou8(struct BlockInfo *bi, u8 *cls, uns len){ + struct BlockInfo* pbi; + u8* pcls; + + for(pbi=bi, pcls=cls; pbicls_num; +} + +void +u8ToBlockInfo(u8 *cls, struct BlockInfo *bi, uns len){ + struct BlockInfo* pbi; + u8* pcls; + + for(pbi=bi, pcls=cls; pbicls_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; iinit_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_diffbi_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; pbibi_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_errthreshold) + return old_err; + cls_num=1; + while(cls_nummax_cls_num){ + cls_num++; + err=k_means(cls_num, bi, dii); + fprintf(stderr, "\n(%u); %d\n", err, err-old_err); + if(errthreshold) break; + if(errdiff_threshold>old_err) break; + old_err=err; + } + return err; +} + + + + -- 2.39.2