+++ /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