From: Pavel Charvat Date: Mon, 7 Aug 2006 20:54:07 +0000 (+0200) Subject: Dan's image decomposition moved to a new directory X-Git-Tag: holmes-import~609 X-Git-Url: http://mj.ucw.cz/gitweb/?a=commitdiff_plain;h=57a04b2f4e5c827a2e28f1e9599443c6308dfdcc;p=libucw.git Dan's image decomposition moved to a new directory --- diff --git a/images/block_info.c b/images/block_info.c deleted file mode 100644 index a8de2bb3..00000000 --- a/images/block_info.c +++ /dev/null @@ -1,146 +0,0 @@ -/*#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/danf/block_info.c b/images/danf/block_info.c new file mode 100644 index 00000000..a8de2bb3 --- /dev/null +++ b/images/danf/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/danf/decomp.c b/images/danf/decomp.c new file mode 100644 index 00000000..1659efac --- /dev/null +++ b/images/danf/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/danf/img.h b/images/danf/img.h new file mode 100644 index 00000000..e28775e8 --- /dev/null +++ b/images/danf/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/danf/k_means.c b/images/danf/k_means.c new file mode 100644 index 00000000..68a840e1 --- /dev/null +++ b/images/danf/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; +} + + + + diff --git a/images/decomp.c b/images/decomp.c deleted file mode 100644 index 1659efac..00000000 --- a/images/decomp.c +++ /dev/null @@ -1,73 +0,0 @@ -#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 deleted file mode 100644 index e28775e8..00000000 --- a/images/img.h +++ /dev/null @@ -1,44 +0,0 @@ -#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 deleted file mode 100644 index 68a840e1..00000000 --- a/images/k_means.c +++ /dev/null @@ -1,232 +0,0 @@ -#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; -} - - - -