Mercurial > hg > wm
annotate Meerwald/wavelet.c @ 22:d8551fb39a5e default tip
Added tag v0.6 for changeset 1c4ccd635a68
author | Peter Meerwald-Stadler <pmeerw@pmeerw.net> |
---|---|
date | Sat, 28 Jan 2023 23:57:51 +0100 |
parents | bd669312f068 |
children |
rev | line source |
---|---|
0 | 1 #include <stdio.h> |
2 #include <stdlib.h> | |
3 #include <string.h> | |
4 #include <math.h> | |
5 #include "wavelet.h" | |
6 #include <ctype.h> | |
8 | 7 #include <values.h> |
0 | 8 |
9 static int read_char(FILE *fp); | |
10 static int read_int(FILE *fp); | |
11 | |
12 IntImage new_intimage(int width, int height) | |
13 { | |
14 IntImage image; | |
15 | |
16 image = (IntImage) calloc(1,sizeof(struct IntImage_struct)); | |
17 if (image==NULL) goto error; | |
18 image->data = (IntPixel*) calloc(width*height,sizeof(IntPixel)); | |
19 if (image->data==NULL) goto error; | |
20 image->width = width; | |
21 image->height = height; | |
22 image->size = width*height; | |
23 image->bpp = 0; | |
24 image->min_val = (IntPixel) 0; | |
25 image->max_val = (IntPixel) 0; | |
26 | |
27 return image; | |
28 | |
29 error: | |
30 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
31 return NULL; | |
32 } | |
33 | |
34 Image new_image(int width, int height) | |
35 { | |
36 Image image; | |
37 | |
38 image = (Image) calloc(1,sizeof(struct Image_struct)); | |
39 if (image==NULL) goto error; | |
40 image->data = (Pixel*) calloc(width*height,sizeof(Pixel)); | |
41 if (image->data==NULL) goto error; | |
42 image->width = width; | |
43 image->height = height; | |
44 image->size = width*height; | |
45 image->bpp = 0; | |
46 image->min_val = (Pixel) 0; | |
47 image->max_val = (Pixel) 0; | |
48 | |
49 return image; | |
50 | |
51 error: | |
52 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
53 return NULL; | |
54 } | |
55 | |
56 | |
57 void free_intimage(IntImage img) | |
58 { | |
59 if (img) { | |
60 if (img->data) free(img->data); | |
61 free(img); | |
62 } | |
63 } | |
64 | |
65 void free_image(Image img) | |
66 { | |
67 if (img) { | |
68 if (img->data) free(img->data); | |
69 free(img); | |
70 } | |
71 } | |
72 | |
73 /************************************************************************ | |
74 * Functionname: load_intimage | |
75 * -------------------------------------------------------------------- | |
76 * PARAMETER: | |
77 * file: filename of image | |
78 * max_val: scaling of grey values to [0..max_val] | |
79 * | |
80 * RETURN: | |
81 * Image that shoud be loaded, if not possible return NULL | |
82 * -------------------------------------------------------------------- | |
83 * DESCRIPTION: | |
84 * This function loads an IntImage (PGM, ASCII or binary | |
85 * encoded (P5 or P3 format) ) from a file. | |
86 ************************************************************************/ | |
87 | |
88 IntImage load_intimage(char *file, int max_val) | |
89 { | |
90 IntImage img; | |
91 FILE *fp; | |
92 IntPixel *data; | |
20
bd669312f068
suppress warnings, fix link errors
Peter Meerwald-Stadler <pmeerw@pmeerw.net>
parents:
8
diff
changeset
|
93 int width, height, i, ch1, ch2; |
0 | 94 |
95 fp=fopen(file,"rb"); | |
96 if (fp==NULL) goto error; | |
97 | |
98 ch1=getc(fp); | |
99 ch2=getc(fp); | |
100 if (ch1!='P' || (ch2!='5' && ch2!='2')) goto error1; | |
101 | |
102 width=read_int(fp); | |
103 height=read_int(fp); | |
104 if ((width==0) || (height==0) ) goto error1; | |
20
bd669312f068
suppress warnings, fix link errors
Peter Meerwald-Stadler <pmeerw@pmeerw.net>
parents:
8
diff
changeset
|
105 read_int(fp); |
0 | 106 |
107 img=new_intimage(width,height); | |
108 | |
109 img->bpp=8; | |
110 | |
111 data=img->data; | |
112 for (i=0; i<img->size; i++) | |
113 { if (ch2=='5') | |
114 *data=getc(fp); | |
115 else | |
116 *data=read_int(fp); | |
117 data++; | |
118 } | |
119 fclose(fp); | |
120 return img; | |
121 | |
122 error1: | |
123 err_SimpleMessage(err_GetErrorMessage(Error_WrongFileFormat)); | |
124 return NULL; | |
125 error: | |
126 err_SimpleMessage(err_GetErrorMessage(Error_CantOpenFile)); | |
127 return NULL; | |
128 } | |
129 | |
130 /************************************************************************ | |
131 * Functionname: load_image | |
132 * -------------------------------------------------------------------- | |
133 * PARAMETER: | |
134 * file: filename of image | |
135 * max_val: scaling of grey values to [0..max_val] | |
136 * | |
137 * RETURN: | |
138 * Image that shoud be loaded, if not possible return NULL | |
139 * -------------------------------------------------------------------- | |
140 * DESCRIPTION: | |
141 * This function loads an IntImage with load_intimage | |
142 * and then converts to Image. | |
143 ************************************************************************/ | |
8 | 144 Image load_image(char *file, int max_val) |
0 | 145 { |
146 Image img; | |
147 IntImage intimg; | |
148 | |
149 intimg = load_intimage(file, max_val); | |
150 if (!intimg) return NULL; | |
151 img = intimage_to_image(intimg); | |
152 if (!intimg) return NULL; | |
153 | |
154 return img; | |
155 } | |
156 | |
157 /************************************************************************/ | |
158 /* Functionname: save_image_P5 */ | |
159 /* -------------------------------------------------------------------- */ | |
160 /* Parameter: */ | |
161 /* img: Image that shoud be saved */ | |
162 /* file: filename of image */ | |
163 /* -------------------------------------------------------------------- */ | |
164 /* Description: save an image as PGM (P5 binary decoded) file */ | |
165 /* */ | |
166 /************************************************************************/ | |
167 | |
8 | 168 int save_image_P5(char *file, Image img) |
0 | 169 { FILE *fp; |
170 Pixel *data; | |
171 long i; | |
172 int p; | |
173 | |
174 fp=fopen(file,"wb"); | |
175 if (fp==NULL) | |
176 goto error; | |
177 fprintf(fp,"P5\n"); | |
178 fprintf(fp,"%d %d\n%d ",img->width,img->height,255); | |
179 data=img->data; | |
180 for (i=0;i<img->size;i++) { | |
181 p=floor(*data+0.5); | |
182 if (p<0) p=0; | |
183 if (p>255) p=255; | |
184 /* putc(*data,fp); */ | |
185 putc(p,fp); | |
186 data++; | |
187 } | |
188 fclose(fp); | |
189 return 1; | |
190 | |
191 error: | |
192 err_SimpleMessage(err_GetErrorMessage(Error_CantOpenFile)); | |
193 return 0; | |
194 } | |
195 | |
196 void clear_image(Image img) | |
197 { | |
198 int i; | |
199 | |
200 PreCondition(img!=NULL,"image==NULL"); | |
201 | |
202 for (i=0;i<img->size;i++) | |
203 (img->data)[i]=(Pixel) 0; | |
204 } | |
205 | |
8 | 206 void copy_into_image(Image img1,Image img2,int x,int y) |
0 | 207 /* copy img2 into img1 at position (x,y)*/ |
208 { | |
209 int start,i,j,aim; | |
210 Pixel *temp; | |
211 | |
212 temp=img2->data; | |
213 start=img1->width*y+x; | |
214 for (i=0;i<img2->height;i++) { | |
215 for (j=0;j<img2->width;j++) { | |
216 aim=start+j+img1->width*i; | |
217 img1->data[aim]=*temp; | |
218 temp++; | |
219 } | |
220 } | |
221 } | |
222 | |
8 | 223 void copy_into_intimage(IntImage img1,IntImage img2,int x,int y) |
0 | 224 {/* copy img2 into img1 at position (x,y)*/ |
225 | |
226 int start,i,j,aim; | |
227 IntPixel *temp; | |
228 | |
229 temp=img2->data; | |
230 start=img1->width*y+x; | |
231 for (i=0;i<img2->height;i++) | |
232 { | |
233 for (j=0;j<img2->width;j++) | |
234 { | |
235 aim=start+j+img1->width*i; | |
236 img1->data[aim]=*temp; | |
237 temp++; | |
238 } | |
239 } | |
240 } | |
241 | |
242 void copy_part_of_image_into_image(Image dest_img, int dest_x, int dest_y, | |
243 Image src_img, int src_x, int src_y, | |
244 int width, int height) | |
245 { | |
246 Pixel *sp,*dp; | |
247 int y,siz; | |
248 | |
249 sp=get_pixel_adr(src_img,src_x,src_y); | |
250 dp=get_pixel_adr(dest_img,dest_x,dest_y); | |
251 | |
252 siz=width*sizeof(Pixel); | |
253 | |
254 for (y=0;y<height;y++) | |
255 { | |
256 memcpy(dp,sp,siz); | |
257 sp+=src_img->width; | |
258 dp+=dest_img->width; | |
259 } | |
260 } | |
261 | |
8 | 262 void copy_part_of_image(Image img1,Image img2,int x,int y) |
0 | 263 /* copy part of img2 begining at position (x,y) into img1 */ |
264 { int i,j,width,height,start,step; | |
265 Pixel *data; | |
266 | |
267 width=img1->width; | |
268 height=img1->height; | |
269 start=img2->width*y+x; | |
270 data=img1->data; | |
271 for (i=0;i<height;i++) { | |
272 step=i*img2->width; | |
273 for (j=0;j<width;j++){ | |
274 *data=img2->data[start+j+step]; | |
275 data++; | |
276 } | |
277 } | |
278 } | |
279 | |
280 | |
8 | 281 void scale_image(Image img, int maximum) |
0 | 282 /* scale image to [0..maximum]*/ |
283 { int i; | |
8 | 284 Pixel max = MINDOUBLE, min = MAXDOUBLE, multi; |
0 | 285 |
286 for (i=0;i<img->size;i++) { | |
287 if (img->data[i]<min) min=img->data[i]; | |
288 else if (img->data[i]>max) max=img->data[i]; | |
289 } | |
290 | |
291 multi=(Pixel)maximum/(max-min); | |
292 for (i=0;i<img->size;i++) img->data[i]=multi*(img->data[i]-min); | |
293 img->min_val=0; | |
294 img->max_val=maximum; | |
295 } | |
296 | |
297 int string_to_pixel(char *str, Pixel *p) | |
298 { | |
299 float ppp; | |
300 if (sscanf(str,"%f",&ppp)==1) | |
301 { | |
302 *p=(Pixel) ppp; | |
303 return 1; | |
304 } | |
305 else | |
306 { | |
307 *p=0.0; | |
308 return 0; | |
309 } | |
310 } | |
311 | |
312 Image intimage_to_image(IntImage i) | |
313 { Image img; | |
314 int j; | |
315 | |
316 img=new_image(i->width,i->height); | |
317 if (img==NULL) goto error; | |
318 for (j=0;j<i->size;j++) | |
319 img->data[j]=(Pixel)i->data[j]; | |
320 img->bpp=i->bpp; | |
321 return img; | |
322 | |
323 error: | |
324 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
325 return NULL; | |
326 } | |
327 IntImage image_to_intimage(Image i) | |
328 { IntImage img; | |
329 int j,multi=1,max,min; | |
330 | |
331 img=(IntImage)calloc(1,sizeof(struct IntImage_struct)); | |
332 if (img==NULL) goto error; | |
333 img->data=(IntPixel *)calloc(i->size,sizeof(IntPixel)); | |
334 if (img->data==NULL) goto error; | |
335 img->width=i->width; | |
336 img->height=i->height; | |
337 img->size=i->size; | |
338 img->bpp=i->bpp; | |
339 | |
340 max=i->max_val; | |
341 min=i->min_val; | |
342 if ((max-min)!=0) | |
343 multi=255.0/(max-min); | |
344 | |
345 for (j=0;j<img->size;j++) | |
346 img->data[j]=(int)((i->data[j]-min)*multi+0.5); | |
347 return img; | |
348 | |
349 error: | |
350 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
351 return NULL; | |
352 | |
353 } | |
354 | |
355 static int read_char(FILE *fp) | |
356 /*read a character from file, but skip any comments*/ | |
357 { int ch; | |
358 | |
359 ch=getc(fp); | |
360 if (ch=='#'){ | |
361 do { | |
362 ch=getc(fp); | |
363 } while (ch!='\n' && ch!=EOF); | |
364 } | |
365 return ch; | |
366 } | |
367 | |
368 | |
369 static int read_int(FILE *fp) | |
370 /*read an ascii integer from file, and skip leading tabstops,new lines ...*/ | |
371 { int r,ch; | |
372 | |
373 do { | |
374 ch=read_char(fp); | |
375 } while (ch==' ' || ch=='\n' || ch=='\t'); | |
376 | |
377 if (ch<'0' || ch>'9') | |
378 goto error; | |
379 | |
380 r=ch-'0'; | |
381 while ( (ch=read_char(fp)) >='0' && (ch <= '9') ) { | |
382 r*=10; | |
383 r+=ch-'0'; | |
384 } | |
385 return r; | |
386 error: | |
387 return 0; | |
388 } | |
389 | |
390 Image_tree new_image_tree() | |
391 { | |
392 Image_tree t; | |
393 t=(Image_tree) calloc(1,sizeof(struct Image_tree_struct)); | |
394 t->entropy=0.0; | |
395 t->image=NULL; | |
396 t->coarse=t->horizontal=t->vertical=t->diagonal=t->doubletree=NULL; | |
397 t->level=0; | |
398 t->codec_data=NULL; | |
399 t->significance_map=NULL; | |
400 return t; | |
401 } | |
402 | |
403 void free_image_tree(Image_tree t) | |
404 { | |
405 if (t->coarse) free_image_tree(t->coarse); | |
406 if (t->horizontal) free_image_tree(t->horizontal); | |
407 if (t->vertical) free_image_tree(t->vertical); | |
408 if (t->diagonal) free_image_tree(t->diagonal); | |
409 if (t->doubletree) free_image_tree(t->doubletree); | |
410 if (t->image) free_image(t->image); | |
411 if (t->significance_map) free_intimage(t->significance_map); | |
412 if (t->codec_data) free(t->codec_data); | |
413 t->image=NULL; | |
414 free(t); | |
415 } | |
416 | |
417 /*********************************************************************** | |
418 * Functionname: get_image_infos | |
419 * -------------------------------------------------------------------- | |
420 * Parameter: | |
421 * Image image: input image | |
422 * Pixel *min,*max,*avg,*var: return minimum, maximum, | |
423 * average and variance of current image | |
424 * -------------------------------------------------------------------- | |
425 * Description: | |
426 * get statistical information of Image | |
427 ************************************************************************/ | |
428 | |
429 void get_image_infos(Image image, Image_info info) | |
430 { | |
431 int x,y; | |
432 Pixel p,sp,sp2; | |
433 | |
434 sp=sp2=(Pixel)0.0; | |
435 | |
436 p=get_pixel(image,0,0); | |
437 | |
438 info->min=info->max=p; | |
439 | |
440 for (y=0;y<image->height;y++) | |
441 for (x=0;x<image->width;x++) | |
442 { | |
443 p=get_pixel(image,x,y); | |
444 info->max=MAX(info->max,p); | |
445 info->min=MIN(info->min,p); | |
446 sp+=p; | |
447 sp2+=p*p; | |
448 } | |
449 sp=sp/image->width/image->height; | |
450 sp2=sp2/image->width/image->height; | |
451 | |
452 info->mean=sp; | |
453 info->var=sp2-sp*sp; | |
454 info->rms=sqrt(sp2); | |
455 } | |
456 | |
457 /*********************************************************************** | |
458 * Functionname: get_difference_image | |
459 * -------------------------------------------------------------------- | |
460 * Parameter: | |
461 * Image image1, image 2: input images | |
462 * | |
463 * Return: | |
464 * Image : difference of image1 and image 2, | |
465 * NULL if error occurs | |
466 ************************************************************************/ | |
467 | |
468 Image get_difference_image(Image image1, Image image2) | |
469 { | |
470 Image diff; | |
471 int i,max,w,h; | |
472 Pixel *d,*i1,*i2; | |
473 | |
474 if ((!image1) || (!image2)) return NULL; | |
475 | |
476 w=image1->width; | |
477 h=image1->height; | |
478 | |
479 if (image2->width != w || image2->height != h) return NULL; | |
480 | |
481 diff=new_image(w,h); | |
482 max=w*h; | |
483 | |
484 d=diff->data; | |
485 i1=image1->data; | |
486 i2=image2->data; | |
487 | |
488 for (i=0;i<max;i++) | |
489 d[i]=i2[i]-i1[i]; | |
490 | |
491 return diff; | |
492 } | |
493 | |
494 | |
495 /************************************************************************/ | |
496 /* Functionname: get_intimage_infos */ | |
497 /* -------------------------------------------------------------------- */ | |
498 /* Parameter: */ | |
499 /* IntImage image: input image */ | |
500 /* IntPixel *min,*max, return minimum, maximum */ | |
501 /* Pixel *avg,*var: average and variance of current image */ | |
502 /* average and variance of current image */ | |
503 /* -------------------------------------------------------------------- */ | |
504 /* Description: */ | |
505 /* get statistical information of Image */ | |
506 /************************************************************************/ | |
507 | |
508 void get_intimage_infos(IntImage image, IntPixel *min, IntPixel *max, Pixel *avg, Pixel *var) | |
509 { | |
510 int x,y; | |
511 Pixel p,sp,sp2; | |
512 | |
513 sp=sp2=(Pixel)0.0; | |
514 | |
515 p= (Pixel) get_intpixel(image,0,0); | |
516 *min=*max=p; | |
517 | |
518 for (y=0;y<image->height;y++) | |
519 for (x=0;x<image->width;x++) | |
520 { | |
521 p= (Pixel) get_intpixel(image,x,y); | |
522 *max=MAX(*max, (IntPixel) p); | |
523 *min=MIN(*min, (IntPixel) p); | |
524 sp+=p; | |
525 sp2+=p*p; | |
526 } | |
527 sp=sp/image->width/image->height; | |
528 sp2=sp2/image->width/image->height; | |
529 | |
530 *avg=sp; | |
531 *var=sp2-sp*sp; | |
532 } | |
533 | |
534 /************************************************************************/ | |
535 /* Functionname: init_zigzag */ | |
536 /* -------------------------------------------------------------------- */ | |
537 /* Parameter: */ | |
538 /* Zigzag_data_struct: */ | |
539 /* output: will be initialized, x/y hold coordinates of */ | |
540 /* the first pixel */ | |
541 /* int width,height: */ | |
542 /* input: width/height of image: */ | |
543 /* -------------------------------------------------------------------- */ | |
544 /* Description: */ | |
545 /* initializes Zigzag_data structure for use with next_zigzag */ | |
546 /************************************************************************/ | |
547 | |
548 void init_zigzag(Zigzag_data zz, int width, int height) | |
549 { | |
550 zz->x=0; | |
551 zz->y=0; | |
552 zz->dir=zigzag_up; | |
553 zz->w=width; | |
554 zz->h=height; | |
555 } | |
556 | |
557 /************************************************************************/ | |
558 /* Functionname: next_zigzag */ | |
559 /* -------------------------------------------------------------------- */ | |
560 /* Parameter: */ | |
561 /* Zigzag_data_struct: */ | |
562 /* int x,y: */ | |
563 /* input: current position of zigzag-scan */ | |
564 /* output: next position of zigzag-scan */ | |
565 /* int w,h: width and height of image */ | |
566 /* enum zigzag_direction *dir: i/o: */ | |
567 /* direction moving thru the image */ | |
568 /* -------------------------------------------------------------------- */ | |
569 /* Description: */ | |
570 /* calculates the next point (x',y') of the zigzag-scan */ | |
571 /* through the image with size (w,h) */ | |
572 /************************************************************************/ | |
573 | |
574 | |
575 void next_zigzag(Zigzag_data zz) | |
576 { | |
577 switch(zz->dir) | |
578 { | |
579 case zigzag_up: | |
580 if (zz->y==0) | |
581 { | |
582 if (zz->x==zz->w-1) | |
583 { | |
584 (zz->y)++; zz->dir=zigzag_down; | |
585 } | |
586 else | |
587 { | |
588 (zz->x)++; zz->dir=zigzag_down; | |
589 } | |
590 } | |
591 else | |
592 { | |
593 if (zz->x==zz->w-1) | |
594 { | |
595 (zz->y)++; zz->dir=zigzag_down; | |
596 } | |
597 else | |
598 { | |
599 (zz->x)++; (zz->y)--; | |
600 } | |
601 } | |
602 break; | |
603 | |
604 case zigzag_down: | |
605 | |
606 if (zz->x==0) | |
607 { | |
608 if (zz->y==zz->h-1) | |
609 { | |
610 (zz->x)++; zz->dir=zigzag_up; | |
611 } | |
612 else | |
613 { | |
614 (zz->y)++; zz->dir=zigzag_up; | |
615 } | |
616 } | |
617 else | |
618 { | |
619 if (zz->y==zz->h-1) | |
620 { | |
621 (zz->x)++; zz->dir=zigzag_up; | |
622 } | |
623 else | |
624 { | |
625 (zz->x)--;(zz->y)++; | |
626 } | |
627 } | |
628 break; | |
629 } | |
630 } | |
631 | |
632 Image get_absolute_image_scaled(Image img) | |
633 { | |
634 Image out; | |
635 | |
636 int x,y; | |
637 | |
638 struct Image_info_struct info; | |
639 Pixel scale,p; | |
640 | |
641 out=new_image(img->width,img->height); | |
642 | |
643 get_image_infos(img, &info); | |
644 | |
645 scale=255/MAX(fabs(info.min),fabs(info.max)); | |
646 | |
647 for(y=0;y<img->height;y++) | |
648 for(x=0;x<img->width;x++) | |
649 { | |
650 p=get_pixel(img,x,y)*scale; | |
651 set_pixel(out,x,y,p); | |
652 } | |
653 return out; | |
654 } | |
655 | |
656 #define FLOOR_HALF(x) ((x)&1 ? ((x)-1)/2 : (x)/2) | |
657 #define CEILING_HALF(x) ((x)&1 ? ((x)+1)/2 : (x)/2) | |
658 | |
659 #define MOD(a,b) ( (a)<0 ? ((b)-((-(a))%(b))) : (a)%(b) ) | |
660 | |
661 Filter new_filter(int size) | |
662 { | |
663 Filter f; | |
664 | |
665 Entering; | |
666 f=(Filter) calloc(1,sizeof(struct FilterStruct)); | |
667 f->data=(Pixel *)calloc(size,sizeof(Pixel)); | |
668 f->len=size; | |
669 f->hipass=0; | |
670 | |
671 Leaving; | |
672 return f; | |
673 } | |
674 | |
675 Pixel get_filter_center(Filter f) | |
676 { | |
677 int i; | |
678 Pixel p, sum, norm; | |
679 | |
680 if (f==NULL) return 0; | |
681 | |
682 sum=norm=0; | |
683 | |
684 for (i=0;i<f->len;i++) | |
685 { | |
686 p=f->data[i]; | |
687 p=p*p; | |
688 norm += p; | |
689 sum += (i+f->start)*p; | |
690 } | |
691 p=sum/norm; | |
692 | |
693 return p; | |
694 } | |
695 int filter_cutoff(Image in, int in_start, int in_len, int in_step, | |
696 Image out, int out_start, int out_len, int out_step, | |
697 Filter f) | |
698 { | |
699 int i,i2,j; | |
700 Pixel *out_pix, *in_pix, *f_data; | |
701 int fstart,fend; /* Source interval */ | |
702 | |
703 Entering; | |
704 | |
705 PreCondition(out_len == in_len/2,"out_len != in_len/2 !!!"); | |
706 | |
707 /* convolution: out[i]=sum_{j=start}^{end} (in[2*i-j]*f[j]) | |
708 | |
709 boundaries: image in [in_start ... in_start + in_len-1] | |
710 image out [out_start ... out_start + out_len-1] | |
711 filter f [0..f->len-1] = [f->start .. f->end] | |
712 cutoff at: | |
713 */ | |
714 | |
715 for (i=0;i<out_len;i++) | |
716 { | |
717 i2=2*i; | |
718 | |
719 fstart=i2-(in_len-1); | |
720 fstart=MAX(fstart,f->start); | |
721 fend=MIN(i2,f->end); | |
722 | |
723 #ifdef TRACE | |
724 sprintf(dbgstr,"i=%d fstart=%d fend=%d\n",i,fstart,fend); | |
725 Trace(dbgstr); | |
726 #endif | |
727 | |
728 out_pix=out->data+out_start+i*out_step; | |
729 | |
730 in_pix=in->data+in_start+(i2-fstart)*in_step; | |
731 | |
732 f_data=f->data-f->start+fstart; | |
733 | |
734 for (j=fstart;j<=fend;j++,in_pix-=in_step,f_data++) | |
735 { | |
736 *out_pix += (*f_data) * (*in_pix); | |
737 #ifdef TRACE | |
738 | |
739 sprintf(dbgstr," j=%2d in: %4.2f filter: %4.2f [%4d] [%4d] opt : %4.2f %4.2f\n", | |
740 j, | |
741 in->data[in_start+in_step*(i2-j)], | |
742 f->data[j-f->start], | |
743 in_start+in_step*(i2-j), | |
744 j-f->start, | |
745 *in_pix, *f_data); | |
746 Trace(dbgstr); | |
747 #endif | |
748 } | |
749 } | |
750 | |
751 Leaving; | |
752 return 1; | |
753 } | |
754 | |
755 | |
756 int filter_inv_cutoff(Image in, int in_start, int in_len, int in_step, | |
757 Image out, int out_start, int out_len, int out_step, | |
758 Filter f) | |
759 { | |
760 int i,j; | |
761 Pixel *out_pix, *in_pix, *f_data; | |
762 int fstart,fend; /* Source interval */ | |
763 Entering; | |
764 PreCondition(out_len == in_len*2,"out_len != in_len*2 !!!"); | |
765 | |
766 /* convolution: out[i]=sum_{j=start}^{end} (f[2*j-i]*in[j]) | |
767 | |
768 boundaries: image in [in_start ... in_start + in_len-1] | |
769 image out [out_start ... out_start + out_len-1] | |
770 filter f [0..f->len-1] = [f->start .. f->end] | |
771 cutoff at: | |
772 */ | |
773 | |
774 for (i=0;i<out_len;i++) | |
775 { | |
776 fstart=CEILING_HALF(f->start+i); | |
777 fend=FLOOR_HALF(f->end+i); | |
778 fstart=MAX(fstart,0); | |
779 fend=MIN(fend,in_len-1); | |
780 | |
781 #ifdef TRACE | |
782 sprintf(dbgstr,"i=%d fstart=%d fend=%d\n",i,fstart,fend); | |
783 Trace(dbgstr); | |
784 #endif | |
785 out_pix=out->data+out_start+i*out_step; | |
786 | |
787 in_pix=in->data+in_start+fstart*in_step; | |
788 | |
789 f_data=f->data-f->start+2*fstart-i; | |
790 | |
791 for (j=fstart;j<=fend;j++,in_pix+=in_step,f_data+=2) | |
792 { | |
793 *out_pix += (*f_data) * (*in_pix); | |
794 #ifdef TRACE | |
795 sprintf(dbgstr," j=%2d in: %4.2f filter: %4.2f [%4d] [%4d] opt : %4.2f %4.2f\n", | |
796 j, | |
797 in->data[in_start+j*in_step], | |
798 f->data[2*j-i-f->start], | |
799 in_start+j*in_step, | |
800 2*j-i-f->start, | |
801 *in_pix, *f_data); | |
802 Trace(dbgstr); | |
803 #endif | |
804 } | |
805 } | |
806 | |
807 Leaving; | |
808 return 1; | |
809 } | |
810 | |
811 int filter_periodical(Image in, int in_start, int in_len, int in_step, | |
812 Image out, int out_start, int out_len, int out_step, | |
813 Filter f) | |
814 { | |
815 int i,i2,j; | |
816 Pixel *out_pix, *in_pix, *f_data; | |
817 int fstart,fend; | |
818 int istart; | |
819 int ipix_len; | |
820 | |
821 Entering; | |
822 PreCondition(out_len == in_len/2,"out_len != in_len/2 !!!"); | |
823 | |
824 /* convolution: out[i]=sum_{j=start}^{end} (in[2*i-j]*f[j]) | |
825 | |
826 boundaries: image in [in_start ... in_start + in_len-1] | |
827 image out [out_start ... out_start + out_len-1] | |
828 filter f [0..f->len-1] = [f->start .. f->end] | |
829 */ | |
830 | |
831 ipix_len=in_len*in_step; | |
832 | |
833 for (i=0;i<out_len;i++) | |
834 { | |
835 i2=2*i; | |
836 | |
837 fstart=f->start; | |
838 fend=f->end; | |
839 | |
840 istart=(i2-fstart); | |
841 istart=MOD(istart,in_len); | |
842 | |
843 #ifdef TRACE | |
844 sprintf(dbgstr,"i=%d istart=%d\n",i,istart); | |
845 Trace(dbgstr); | |
846 #endif | |
847 | |
848 out_pix=out->data+out_start+i*out_step; | |
849 | |
850 in_pix=in->data+in_start+istart*in_step; | |
851 | |
852 f_data=f->data; | |
853 | |
854 for (j=fstart;j<=fend;j++,f_data++) | |
855 { | |
856 *out_pix += (*f_data) * (*in_pix); | |
857 #ifdef TRACE | |
858 | |
859 sprintf(dbgstr," j=%2d in: %4.2f filter: %4.2f [%4d] [%4d] opt : %4.2f %4.2f\n", | |
860 j, | |
861 in->data[in_start+in_step*((i2-j+in_len)%in_len)], | |
862 f->data[j-f->start], | |
863 in_start+in_step*((i2-j+in_len)%in_len), | |
864 j-f->start, | |
865 *in_pix, *f_data); | |
866 Trace(dbgstr); | |
867 #endif | |
868 in_pix-=in_step; | |
869 istart--; | |
870 if (istart<0) | |
871 { | |
872 istart+=in_len; | |
873 in_pix+=ipix_len; | |
874 } | |
875 } | |
876 } | |
877 | |
878 Leaving; | |
879 return 1; | |
880 } | |
881 | |
882 int filter_inv_periodical(Image in, int in_start, int in_len, int in_step, | |
883 Image out, int out_start, int out_len, int out_step, | |
884 Filter f) | |
885 { | |
886 int i,j; | |
887 Pixel *out_pix, *in_pix, *f_data; | |
888 int fstart,fend; /* Source interval */ | |
889 int istart; | |
890 int ipix_len; | |
891 Entering; | |
892 PreCondition(out_len == in_len*2,"out_len != in_len*2 !!!"); | |
893 | |
894 /* convolution: out[i]=sum_{j=start}^{end} (f[2*j-i]*in[j]) | |
895 | |
896 boundaries: image in [in_start ... in_start + in_len-1] | |
897 image out [out_start ... out_start + out_len-1] | |
898 filter f [0..f->len-1] = [f->start .. f->end] | |
899 */ | |
900 | |
901 ipix_len=in_len*in_step; | |
902 | |
903 for (i=0;i<out_len;i++) | |
904 { | |
905 fstart=CEILING_HALF(f->start+i); | |
906 fend=FLOOR_HALF(f->end+i); | |
907 | |
908 istart=MOD(fstart,in_len); | |
909 #ifdef TRACE | |
910 sprintf(dbgstr,"i=%d fstart=%d fend=%d istart=%d\n",i,fstart,fend,istart); | |
911 Trace(dbgstr); | |
912 #endif | |
913 out_pix=out->data+out_start+i*out_step; | |
914 | |
915 in_pix=in->data+in_start+istart*in_step; | |
916 | |
917 f_data=f->data-f->start+2*fstart-i; | |
918 | |
919 for (j=fstart;j<=fend;j++,f_data+=2) | |
920 { | |
921 *out_pix += (*f_data) * (*in_pix); | |
922 #ifdef TRACE | |
923 sprintf(dbgstr," j=%2d in: %4.2f filter: %4.2f [%4d] [%4d] opt : %4.2f %4.2f\n", | |
924 j, | |
925 in->data[in_start+(j % in_len)*in_step], | |
926 f->data[2*j-i-f->start], | |
927 in_start+(j%in_len)*in_step, | |
928 2*j-i-f->start, | |
929 *in_pix, *f_data); | |
930 Trace(dbgstr); | |
931 #endif | |
932 in_pix+=in_step; | |
933 istart++; | |
934 if (istart>=in_len) | |
935 { | |
936 istart-=in_len; | |
937 in_pix-=ipix_len; | |
938 } | |
939 } | |
940 } | |
941 | |
942 Leaving; | |
943 return 1; | |
944 } | |
945 | |
946 int filter_mirror(Image in, int in_start, int in_len, int in_step, | |
947 Image out, int out_start, int out_len, int out_step, | |
948 Filter f) | |
949 { | |
950 int i,i2,j; | |
951 Pixel *out_pix, *in_pix, *f_data; | |
952 int fstart,fend; | |
953 int in_pos; | |
954 | |
955 Entering; | |
956 PreCondition(out_len == in_len/2,"out_len != in_len/2 !!!"); | |
957 | |
958 /* convolution: out[i]=sum_{j=start}^{end} (in[2*i-j]*f[j]) | |
959 | |
960 boundaries: image in [in_start ... in_start + in_len-1] | |
961 image out [out_start ... out_start + out_len-1] | |
962 filter f [0..f->len-1] = [f->start .. f->end] | |
963 */ | |
964 | |
965 in_pix=in->data+in_start; | |
966 | |
967 for (i=0;i<out_len;i++) | |
968 { | |
969 i2=2*i; | |
970 | |
971 fstart=f->start; | |
972 fend=f->end; | |
973 | |
974 out_pix=out->data+out_start+i*out_step; | |
975 | |
976 f_data=f->data; | |
977 | |
978 for (j=fstart;j<=fend;j++) | |
979 { | |
980 in_pos=(i2-j); | |
981 if (in_pos<0) | |
982 { | |
983 in_pos=-in_pos; | |
984 if (in_pos>=in_len) continue; | |
985 } | |
986 if (in_pos>=in_len) | |
987 { | |
988 in_pos=2*in_len-2-in_pos; | |
989 if (in_pos<0) continue; | |
990 } | |
991 *out_pix += (f_data[j-fstart]) * (in_pix[in_pos*in_step]); | |
992 } | |
993 } | |
994 | |
995 Leaving; | |
996 return 1; | |
997 } | |
998 | |
999 int filter_inv_mirror(Image in, int in_start, int in_len, int in_step, | |
1000 Image out, int out_start, int out_len, int out_step, | |
1001 Filter f) | |
1002 { | |
1003 int i,j; | |
8 | 1004 Pixel *out_pix, *in_pix; |
0 | 1005 int fstart,fend; /* Source interval */ |
8 | 1006 int in_pos; |
0 | 1007 |
1008 Entering; | |
1009 PreCondition(out_len == in_len*2,"out_len != in_len*2 !!!"); | |
1010 | |
1011 /* convolution: out[i]=sum_{j=start}^{end} (f[2*j-i]*in[j]) | |
1012 | |
1013 boundaries: image in [in_start ... in_start + in_len-1] | |
1014 image out [out_start ... out_start + out_len-1] | |
1015 filter f [0..f->len-1] = [f->start .. f->end] | |
1016 */ | |
1017 | |
1018 /*fprintf(stderr,"inv started\n");*/ | |
1019 for (i=0;i<out_len;i++) | |
1020 { | |
1021 fstart=CEILING_HALF(f->start+i); | |
1022 fend=FLOOR_HALF(f->end+i); | |
1023 | |
1024 out_pix=out->data+out_start+i*out_step; | |
1025 | |
1026 in_pix=in->data+in_start; | |
1027 | |
1028 /* | |
1029 printf("in: %4d - %4d flt: %4d - %4d [%s]\n",fstart,fend,2*fstart-i,2*fend-i, | |
1030 (2*fstart-i<f->start || 2*fend-i>f->end) ? "error":"ok"); | |
1031 */ | |
1032 /*fprintf(stderr,"inv[%d]\n",i);*/ | |
1033 for (j=fstart;j<=fend;j++) | |
1034 { | |
1035 in_pos=j; | |
1036 if (in_pos<0) | |
1037 { | |
1038 if (f->hipass) | |
1039 in_pos=-in_pos-1; | |
1040 else | |
1041 in_pos=-in_pos; | |
1042 if (in_pos>=in_len) continue; | |
1043 } | |
1044 if (in_pos>=in_len) | |
1045 { | |
1046 if (f->hipass) | |
1047 in_pos=2*in_len-2-in_pos; | |
1048 else | |
1049 in_pos=2*in_len-1-in_pos; | |
1050 if (in_pos<0) continue; | |
1051 } | |
1052 /*fprintf(stderr,"out+= %7.2f * %7.2f = %7.2f\n",f->data[2*j-i-f->start],in_pix[in_pos*in_step],f->data[2*j-i-f->start]*in_pix[in_pos*in_step]);*/ | |
1053 *out_pix += f->data[2*j-i-f->start] * (in_pix[in_pos*in_step]); | |
1054 } | |
1055 } | |
1056 | |
1057 Leaving; | |
1058 return 1; | |
1059 } | |
1060 | |
1061 #define MAX_LINE 256 | |
1062 | |
1063 #define skip_blank(str) { while(isspace(*(str))) (str)++; } | |
1064 | |
1065 static int get_next_line(FILE *f, char *c) | |
1066 { | |
1067 char *str,string[200]; | |
1068 int len; | |
1069 do | |
1070 { | |
1071 str=string; | |
1072 if (!fgets(str,200,f)) | |
1073 { | |
1074 Trace("get_next_line: eof\n"); | |
1075 goto error; | |
1076 } | |
1077 len=strlen(str); | |
1078 while (len>=1 && isspace(str[len-1])) str[--len]=0; | |
1079 while (isspace(*str)) str++; | |
1080 } | |
1081 while (strlen(str)==0 || *str=='#'); | |
1082 strcpy(c,str); | |
1083 return 1; | |
1084 error: | |
1085 return 0; | |
1086 } | |
1087 | |
1088 static int next_line_str(FILE *f, char *tag, char *out) | |
1089 { | |
1090 char str[MAX_LINE],*t_str; | |
1091 | |
1092 if (!get_next_line(f,str)) goto error; | |
1093 t_str=strtok(str," "); | |
1094 if (!t_str || strcmp(t_str,tag)) goto error; | |
1095 t_str=strtok(NULL,"\n"); | |
1096 if (!t_str) goto error; | |
1097 skip_blank(t_str); | |
1098 | |
1099 strcpy(out,t_str); | |
1100 return 1; | |
1101 error: | |
1102 return 0; | |
1103 } | |
1104 | |
1105 static int next_line_str_alloc(FILE *f, char *tag, char **out) | |
1106 { | |
1107 char str[MAX_LINE]; | |
1108 if (!next_line_str(f,tag,str)) goto error; | |
1109 | |
1110 *out=malloc(strlen(str)+1); | |
1111 strcpy(*out,str); | |
1112 | |
1113 return 1; | |
1114 error: | |
1115 return 0; | |
1116 } | |
1117 | |
1118 static int next_line_int(FILE *f, char *tag, int *out) | |
1119 { | |
1120 char str[MAX_LINE]; | |
1121 if (next_line_str(f,tag,str) && sscanf(str,"%d",out)==1) | |
1122 return 1; | |
1123 else | |
1124 return 0; | |
1125 } | |
1126 | |
1127 | |
1128 static Filter read_filter(FILE *f) | |
1129 { | |
1130 char str[MAX_LINE]; | |
1131 Filter filter; | |
1132 int i; | |
1133 | |
1134 Entering; | |
1135 | |
1136 filter=calloc(1,sizeof(struct FilterStruct)); | |
1137 | |
1138 if (!next_line_str(f,"Type",str)) goto error1; | |
1139 | |
1140 if (!strcmp(str,"nosymm")) | |
1141 { | |
1142 filter->type=FTNoSymm; | |
1143 } | |
1144 else if (!strcmp(str,"symm")) | |
1145 { | |
1146 filter->type=FTSymm; | |
1147 } | |
1148 else if (!strcmp(str,"antisymm")) | |
1149 { | |
1150 filter->type=FTAntiSymm; | |
1151 } | |
1152 else | |
1153 goto error1; | |
1154 | |
1155 if (!next_line_int(f,"Length",&(filter->len))) goto error1; | |
1156 if (!next_line_int(f,"Start",&(filter->start))) goto error1; | |
1157 if (!next_line_int(f,"End",&(filter->end))) goto error1; | |
1158 | |
1159 if ((filter->end-filter->start+1!=filter->len)) | |
1160 { | |
1161 Trace("error: len != end-start+1\n"); | |
1162 goto error1; | |
1163 } | |
1164 | |
1165 filter->data=calloc(filter->len,sizeof(Pixel)); | |
1166 | |
1167 for (i=0;i<filter->len;i++) | |
1168 { | |
1169 if (!get_next_line(f,str)) goto error2; | |
1170 if (!string_to_pixel(str,filter->data+i)) | |
1171 { | |
1172 Trace("error: invalid filter-value\n"); | |
1173 goto error2; | |
1174 } | |
1175 } | |
1176 if (!get_next_line(f,str)) goto error2; | |
1177 if (*str!='}') | |
1178 { | |
1179 Trace("error: '}' not found\n"); | |
1180 goto error2; | |
1181 } | |
1182 | |
1183 Leaving; | |
1184 return filter; | |
1185 | |
1186 error2: | |
1187 free(filter->data); | |
1188 | |
1189 error1: | |
1190 free(filter); | |
1191 | |
1192 LeavingErr; | |
1193 return NULL; | |
1194 | |
1195 } | |
1196 | |
1197 static FilterGH read_filter_gh(FILE *f) | |
1198 { | |
1199 char str[MAX_LINE]; | |
1200 FilterGH fgh; | |
1201 Filter filter; | |
1202 int i,max; | |
1203 | |
1204 Entering; | |
1205 | |
1206 fgh=calloc(1,sizeof(struct FilterGHStruct)); | |
1207 | |
1208 if (!next_line_str_alloc(f,"Name",&(fgh->name))) | |
1209 { | |
1210 Trace("error: 'Name' tag not found\n"); | |
1211 goto error1; | |
1212 } | |
1213 | |
1214 if (!next_line_str(f,"Type",str)) | |
1215 { | |
1216 Trace("error: 'Type' tag not found\n"); | |
1217 goto error1; | |
1218 } | |
1219 | |
1220 if (!strcmp(str,"orthogonal")) | |
1221 { | |
1222 fgh->type=FTOrtho; | |
1223 max=2; | |
1224 } | |
1225 else if (!strcmp(str,"biorthogonal")) | |
1226 { | |
1227 fgh->type=FTBiOrtho; | |
1228 max=4; | |
1229 } | |
1230 else if (!strcmp(str,"other")) | |
1231 { | |
1232 fgh->type=FTOther; | |
1233 max=4; | |
1234 } | |
1235 else | |
1236 { | |
1237 Trace("error: expecting 'orthogonal', 'biorthogonal' or 'other' type-tag\n"); | |
1238 goto error1; | |
1239 } | |
1240 | |
1241 for (i=0;i<max;i++) | |
1242 { | |
1243 if (!get_next_line(f,str)) goto error2; | |
1244 if (*str!='{') | |
1245 { | |
1246 Trace("error: '{' not found\n"); | |
1247 goto error2; | |
1248 } | |
1249 if (!(filter=read_filter(f))) | |
1250 { | |
1251 Trace("error: read_filter failed\n"); | |
1252 goto error2; | |
1253 } | |
1254 filter->hipass = !(i&1); | |
1255 switch(i) | |
1256 { | |
1257 case 0: fgh->g=filter; break; | |
1258 case 1: fgh->h=filter; break; | |
1259 case 2: fgh->gi=filter; break; | |
1260 case 3: fgh->hi=filter; break; | |
1261 } | |
1262 } | |
1263 if (!get_next_line(f,str)) goto error2; | |
1264 if (*str!='}') | |
1265 { | |
1266 Trace("error: '}' not found\n"); | |
1267 goto error2; | |
1268 } | |
1269 | |
1270 Leaving; | |
1271 return fgh; | |
1272 | |
1273 error2: | |
1274 if (fgh->g) free(fgh->g); | |
1275 if (fgh->h) free(fgh->h); | |
1276 if (fgh->gi) free(fgh->gi); | |
1277 if (fgh->hi) free(fgh->hi); | |
1278 | |
1279 error1: | |
1280 free(fgh); | |
1281 | |
1282 LeavingErr; | |
1283 return NULL; | |
1284 } | |
1285 | |
1286 | |
1287 AllFilters load_filters(char *name) | |
1288 { | |
1289 FILE *f; | |
1290 char str[MAX_LINE]; | |
1291 AllFilters a; | |
1292 FilterGH fgh; | |
1293 | |
1294 Entering; | |
1295 | |
1296 PreCondition(name!=NULL,"name=NULL!"); | |
1297 | |
1298 f=fopen(name,"rt"); | |
1299 if (!f) | |
1300 { | |
1301 Trace("error: fopen failed\n"); | |
1302 goto error1; | |
1303 } | |
1304 | |
1305 a=calloc(1,sizeof(struct AllFilterStruct)); | |
1306 a->count=0; | |
1307 | |
1308 while (get_next_line(f,str)) | |
1309 { | |
1310 if (*str=='{') | |
1311 { | |
1312 fgh=read_filter_gh(f); | |
1313 if (!fgh) | |
1314 { | |
1315 Trace("error: read_filter returned NULL\n"); | |
1316 goto error2; | |
1317 } | |
1318 if (a->count) | |
1319 a->filter=realloc(a->filter,(a->count+1)*sizeof(FilterGH)); | |
1320 else | |
1321 a->filter=malloc(sizeof(FilterGH)); | |
1322 (a->filter)[a->count]=fgh; | |
1323 a->count++; | |
1324 } | |
1325 else | |
1326 { | |
1327 Trace("error: '{' not found\n"); | |
1328 goto error2; | |
1329 } | |
1330 } | |
1331 | |
1332 fclose(f); | |
1333 | |
1334 Leaving; | |
1335 return a; | |
1336 | |
1337 error2: | |
1338 fclose(f); | |
1339 error1: | |
1340 | |
1341 LeavingErr; | |
1342 return 0; | |
1343 } | |
1344 | |
1345 #define doubletree_min 32 | |
1346 #define best_basis_min 8 | |
1347 | |
1348 static int convolute_lines(Image output,Image input,Filter flt,enum FilterMethod method); | |
1349 static int convolute_rows(Image output,Image input,Filter flt,enum FilterMethod method); | |
1350 static int decomposition(Image t_img,Image coarse,Image horizontal,Image vertical, | |
1351 Image diagonal,Filter g,Filter h,enum FilterMethod method); | |
1352 static int compute_best(Image_tree tree,int level,int max_level,FilterGH *flt, | |
1353 enum FilterMethod method,enum Information_Cost cost,double epsilon); | |
1354 static double compute_entropy(Image img,enum Information_Cost cost,double epsilon); | |
8 | 1355 static void compute_levels(Image_tree tree,double *entropies,enum Information_Cost cost,double epsilon); |
1356 static void free_levels(Image_tree tree,int best); | |
0 | 1357 |
1358 static Pixel sumationq(Image img); | |
1359 static Pixel normq(Image_tree tree); | |
1360 static Pixel sumation_down(Image_tree tree, Pixel normq); | |
1361 static Pixel compute_non_additive(Image_tree tree,int size,enum Information_Cost cost,double | |
1362 epsilon,int down); | |
1363 | |
1364 /************************************************************************/ | |
1365 /* Functionname: wavelettransform */ | |
1366 /* -------------------------------------------------------------------- */ | |
1367 /* Parameter: */ | |
1368 /* original: Image that should be transformed */ | |
1369 /* level: transform down to level */ | |
1370 /* flt: transform with filters flt[0..level] */ | |
1371 /* method: method of filtering */ | |
1372 /* -------------------------------------------------------------------- */ | |
1373 /* Description: Carries out the wavelettransform */ | |
1374 /* */ | |
1375 /************************************************************************/ | |
8 | 1376 Image_tree wavelettransform(Image original,int level,FilterGH *flt,enum FilterMethod method) |
0 | 1377 { int i,width,height,min,max_level,e; |
1378 Image coarsei,horizontali,verticali,diagonali,tempi; | |
1379 Image_tree ret_tree,temp_tree; | |
1380 | |
1381 width=original->width; | |
1382 height=original->height; | |
1383 | |
1384 tempi=new_image(width,height); | |
1385 if(!tempi) goto error; | |
1386 | |
1387 copy_into_image(tempi,original,0,0); | |
1388 | |
1389 ret_tree=new_image_tree(); | |
1390 if(!ret_tree) goto error; | |
1391 | |
1392 temp_tree=ret_tree; | |
1393 ret_tree->level=0; | |
1394 | |
1395 min=original->width; | |
1396 if (original->height<min) min=original->height; | |
1397 max_level=log(min)/log(2)-2; | |
1398 if (max_level<level) level=max_level; | |
1399 | |
1400 if (level<1) /* do not transform */ | |
1401 { | |
1402 ret_tree->image=tempi; | |
1403 return ret_tree; | |
1404 } | |
1405 | |
1406 /* decomposition */ | |
1407 | |
1408 for (i=0;i<level;i++) | |
1409 { | |
1410 | |
1411 width=(width+1)/2; | |
1412 height=(height+1)/2; | |
1413 | |
1414 coarsei=new_image(width,height); | |
1415 horizontali=new_image(width,height); | |
1416 verticali=new_image(width,height); | |
1417 diagonali=new_image(width,height); | |
1418 if(!coarsei||!horizontali||!verticali||!diagonali) goto error; | |
1419 | |
1420 e=decomposition(tempi,coarsei,horizontali,verticali,diagonali, | |
1421 flt[i]->g,flt[i]->h,method); | |
1422 if (!e) return NULL; | |
1423 | |
1424 temp_tree->coarse=new_image_tree(); | |
1425 temp_tree->horizontal=new_image_tree(); | |
1426 temp_tree->vertical=new_image_tree(); | |
1427 temp_tree->diagonal=new_image_tree(); | |
1428 | |
1429 temp_tree->coarse->level=i+1; | |
1430 temp_tree->horizontal->level=i+1; | |
1431 temp_tree->vertical->level=i+1; | |
1432 temp_tree->diagonal->level=i+1; | |
1433 | |
1434 temp_tree->horizontal->image=horizontali; | |
1435 temp_tree->vertical->image=verticali; | |
1436 temp_tree->diagonal->image=diagonali; | |
1437 free_image(tempi); | |
1438 | |
1439 if (i!=(level-1)) | |
1440 { | |
1441 tempi=new_image(width,height); | |
1442 copy_into_image(tempi,coarsei,0,0); | |
1443 free_image(coarsei); | |
1444 /*if i=level coarsei is inserted into the image tree | |
1445 so we should not free coarsei on level-1*/ | |
1446 } | |
1447 | |
1448 temp_tree=temp_tree->coarse; | |
1449 | |
1450 } | |
1451 | |
1452 temp_tree->image=coarsei; | |
1453 return ret_tree; | |
1454 | |
1455 error: | |
1456 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
1457 return NULL; | |
1458 } | |
1459 | |
1460 static Image_tree wavelettransform__wp(Image original, int current_level, int level, FilterGH *flt, enum FilterMethod method) | |
1461 { | |
1462 int i, width, height, min, max_level, e; | |
1463 Image coarse_image,horizontal_image,vertical_image,diagonal_image,temp_image; | |
1464 Image_tree return_tree, temp_tree; | |
1465 | |
1466 width = original->width; | |
1467 height = original->height; | |
1468 | |
1469 temp_image = new_image(width, height); | |
1470 if (!temp_image) goto error; | |
1471 | |
1472 copy_into_image(temp_image, original, 0, 0); | |
1473 | |
1474 temp_tree = return_tree = new_image_tree(); | |
1475 if (!return_tree) goto error; | |
1476 | |
1477 temp_tree->level = current_level; | |
1478 | |
1479 min = original->width; | |
1480 if (original->height < min) min = original->height; | |
1481 max_level = log(min) / log(2) - 2; | |
1482 if (max_level < level) level = max_level; | |
1483 | |
1484 if (current_level >= level) { | |
1485 return_tree->image = temp_image; | |
1486 return return_tree; | |
1487 } | |
1488 | |
1489 for (i = current_level; i < level; i++) { | |
1490 width = (width + 1) / 2; | |
1491 height = (height + 1) / 2; | |
1492 | |
1493 coarse_image = new_image(width, height); | |
1494 horizontal_image = new_image(width, height); | |
1495 vertical_image = new_image(width, height); | |
1496 diagonal_image = new_image(width, height); | |
1497 | |
1498 if (!coarse_image || !horizontal_image || | |
1499 !vertical_image || !diagonal_image) goto error; | |
1500 | |
1501 e = decomposition(temp_image, coarse_image, horizontal_image, | |
1502 vertical_image, diagonal_image, | |
1503 flt[i]->g, flt[i]->h, method); | |
1504 if (!e) return NULL; | |
1505 | |
1506 temp_tree->coarse = new_image_tree(); | |
1507 temp_tree->coarse->level = i+1; | |
1508 temp_tree->horizontal = wavelettransform__wp(horizontal_image, i+1, level, flt, method); | |
1509 temp_tree->vertical = wavelettransform__wp(vertical_image, i+1, level, flt, method); | |
1510 temp_tree->diagonal = wavelettransform__wp(diagonal_image, i+1, level, flt, method); | |
1511 | |
1512 free_image(horizontal_image); | |
1513 free_image(vertical_image); | |
1514 free_image(diagonal_image); | |
1515 free_image(temp_image); | |
1516 | |
1517 if (i != (level - 1)) { | |
1518 temp_image = new_image(width, height); | |
1519 copy_into_image(temp_image, coarse_image, 0, 0); | |
1520 free_image(coarse_image); | |
1521 } | |
1522 | |
1523 temp_tree = temp_tree->coarse; | |
1524 } | |
1525 | |
1526 temp_tree->image = coarse_image; | |
1527 return return_tree; | |
1528 | |
1529 error: | |
1530 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
1531 return NULL; | |
1532 } | |
1533 | |
8 | 1534 Image_tree wavelettransform_wp(Image original, int level, FilterGH *flt, enum FilterMethod method) { |
1535 return wavelettransform__wp(original, 0, level, flt, method); | |
0 | 1536 } |
1537 | |
1538 | |
1539 /************************************************************************/ | |
1540 /* Functionname: best_basis */ | |
1541 /* -------------------------------------------------------------------- */ | |
1542 /* Parameter: */ | |
1543 /* original: Image to be transformed */ | |
1544 /* level: best basis selection down to this level */ | |
1545 /* flt: transform with filters flt[0..level] */ | |
1546 /* method: transform with filter method */ | |
1547 /* cost: carry best basis selection out with this costfunc */ | |
1548 /* epsilon: limit for threshold method */ | |
1549 /* -------------------------------------------------------------------- */ | |
1550 /* Description: carries best basis and near best basis selection */ | |
1551 /* out */ | |
1552 /************************************************************************/ | |
8 | 1553 Image_tree best_basis(Image original,int level,FilterGH *flt, |
0 | 1554 enum FilterMethod method,enum Information_Cost cost,double epsilon) |
1555 | |
1556 { Image_tree tree; | |
1557 Image img; | |
20
bd669312f068
suppress warnings, fix link errors
Peter Meerwald-Stadler <pmeerw@pmeerw.net>
parents:
8
diff
changeset
|
1558 int min,max_level; |
0 | 1559 |
1560 tree=new_image_tree(); | |
1561 if(!tree) goto error; | |
1562 | |
1563 img=new_image(original->width,original->height); | |
1564 if(!img) goto error; | |
1565 | |
1566 copy_into_image(img,original,0,0); | |
1567 | |
1568 tree->image=img; | |
1569 | |
1570 min=original->width; | |
1571 if (original->height<min) min=original->height; | |
1572 max_level=log10((float) min/best_basis_min)/log10(2); | |
1573 if (max_level>level) max_level=level; | |
1574 | |
20
bd669312f068
suppress warnings, fix link errors
Peter Meerwald-Stadler <pmeerw@pmeerw.net>
parents:
8
diff
changeset
|
1575 compute_best(tree,0,max_level,flt,method,cost,epsilon); |
0 | 1576 |
1577 if (!tree->image) free(img); | |
1578 | |
1579 return tree; | |
1580 | |
1581 error: | |
1582 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
1583 return NULL; | |
1584 | |
1585 } | |
1586 /************************************************************************/ | |
1587 /* Functionname: best_level_selection */ | |
1588 /* -------------------------------------------------------------------- */ | |
1589 /* Parameter: */ | |
1590 /* original: Image that should be transformed */ | |
1591 /* maxlevel: transform down to level */ | |
1592 /* flt: transform with filters flt[0..level] */ | |
1593 /* method: transform with filter method */ | |
1594 /* cost: carry best basis selection out with this costfunc */ | |
1595 /* epsilon: limit for threshold method */ | |
1596 /* -------------------------------------------------------------------- */ | |
1597 /* Description: Carries out the best level selection */ | |
1598 /* */ | |
1599 /************************************************************************/ | |
8 | 1600 Image_tree best_level(Image original,int maxlevel,int *bestlevel,FilterGH *flt,enum FilterMethod method, |
0 | 1601 enum Information_Cost cost,double epsilon) |
1602 { Image_tree tree; | |
1603 Image img; | |
1604 double *entropies,min; | |
1605 int best=0,i,e; | |
1606 | |
1607 img=new_image(original->width,original->height); | |
1608 copy_into_image(img,original,0,0); | |
1609 | |
1610 tree=new_image_tree(); | |
1611 tree->image=img; | |
1612 | |
1613 entropies=(double *)calloc(maxlevel+1,sizeof(double)); | |
1614 if(!entropies) goto error; | |
1615 | |
1616 /* decompose down to maxlevel */ | |
1617 e=decompose_all(tree,maxlevel,flt,method,cost,epsilon); | |
1618 if (!e) return NULL; | |
1619 | |
1620 /* compute costs of each level and store it in entropies array*/ | |
1621 compute_levels(tree,entropies,cost,epsilon); | |
1622 | |
1623 min=entropies[0]; | |
1624 for (i=1;i<=maxlevel;i++) | |
1625 { | |
1626 if (entropies[i]<min) | |
1627 { | |
1628 min=entropies[i]; | |
1629 best=i; | |
1630 } | |
1631 } | |
1632 | |
1633 /* free all other levels */ | |
1634 free_levels(tree,best); | |
1635 | |
1636 *bestlevel=best; | |
1637 | |
1638 return tree; | |
1639 | |
1640 error: | |
1641 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
1642 return NULL; | |
1643 } | |
1644 | |
1645 /************************************************************************/ | |
1646 /* Functionname: compute_best */ | |
1647 /* -------------------------------------------------------------------- */ | |
1648 /* Parameter: */ | |
1649 /* img: Image that should be transformed */ | |
1650 /* level: transform level */ | |
1651 /* max_level: transform to maximum level */ | |
1652 /* flt: transform with filters flt[0..level] */ | |
1653 /* method: transform with filter method */ | |
1654 /* cost: carry best basis selection out with this costfunc */ | |
1655 /* epsilon: limit for threshold method */ | |
1656 /* -------------------------------------------------------------------- */ | |
1657 /* Description: Carries out the best basis selection (recursivly) */ | |
1658 /* (is needed by the waveletpacket procedure) */ | |
1659 /************************************************************************/ | |
1660 static int compute_best(Image_tree tree,int level,int max_level,FilterGH *flt, | |
1661 enum FilterMethod method,enum Information_Cost cost,double epsilon) | |
1662 | |
1663 { Image coarse,horizontal,vertical,diagonal; | |
1664 int e,width,height; | |
1665 double sum; | |
1666 | |
1667 tree->level=level; | |
1668 | |
1669 /* non additive cost function*/ | |
1670 if (cost>=shanon) | |
1671 { | |
1672 tree->entropy=compute_non_additive(tree,tree->image->size,cost,epsilon,0); | |
1673 } | |
1674 /*additive cost function*/ | |
1675 else tree->entropy=compute_entropy(tree->image,cost,epsilon); | |
1676 | |
1677 if (level<max_level) { | |
1678 width=(tree->image->width+1)/2; | |
1679 height=(tree->image->height+1)/2; | |
1680 | |
1681 tree->coarse=new_image_tree(); | |
1682 tree->horizontal=new_image_tree(); | |
1683 tree->vertical=new_image_tree(); | |
1684 tree->diagonal=new_image_tree(); | |
1685 | |
1686 coarse=new_image(width,height); | |
1687 horizontal=new_image(width,height); | |
1688 vertical=new_image(width,height); | |
1689 diagonal=new_image(width,height); | |
1690 if(!coarse||!vertical||!horizontal||!diagonal) goto error; | |
1691 | |
1692 e=decomposition(tree->image,coarse,horizontal,vertical,diagonal,flt[0]->g,flt[0]->h,method); | |
1693 if (!e) return 0; | |
1694 | |
1695 tree->coarse->image=coarse; | |
1696 tree->horizontal->image=horizontal; | |
1697 tree->vertical->image=vertical; | |
1698 tree->diagonal->image=diagonal; | |
1699 | |
1700 e=compute_best(tree->coarse,level+1,max_level,flt+1,method,cost,epsilon); | |
1701 e=compute_best(tree->horizontal,level+1,max_level,flt+1,method,cost,epsilon); | |
1702 e=compute_best(tree->vertical,level+1,max_level,flt+1,method,cost,epsilon); | |
1703 e=compute_best(tree->diagonal,level+1,max_level,flt+1,method,cost,epsilon); | |
1704 if (!e) return 0; | |
1705 | |
1706 /*going back in recursion*/ | |
1707 | |
1708 if (cost>=shanon) { | |
1709 sum=compute_non_additive(tree,tree->image->size,cost,epsilon,1); | |
1710 } | |
1711 else sum=(tree->coarse->entropy)+(tree->horizontal->entropy) | |
1712 +(tree->vertical->entropy)+(tree->diagonal->entropy); | |
1713 | |
1714 if (tree->entropy>sum) | |
1715 { | |
1716 tree->entropy=sum; | |
1717 free_image(tree->image); /* take down tree */ | |
1718 tree->image=NULL; | |
1719 | |
1720 } | |
1721 else | |
1722 { /* delete the tree downwards */ | |
1723 free_image_tree(tree->coarse); | |
1724 free_image_tree(tree->horizontal); | |
1725 free_image_tree(tree->vertical); | |
1726 free_image_tree(tree->diagonal); | |
1727 | |
1728 tree->coarse=tree->vertical=tree->horizontal=tree->diagonal=NULL; | |
1729 } | |
1730 } | |
1731 | |
1732 return 1; | |
1733 | |
1734 error: | |
1735 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
1736 return 0; | |
1737 | |
1738 } | |
1739 | |
1740 /************************************************************************/ | |
1741 /* Functionname: decompose_all */ | |
1742 /* -------------------------------------------------------------------- */ | |
1743 /* Parameter: */ | |
1744 /* tree: Image tree to be decomposed */ | |
1745 /* maxlevel: decompose down to level */ | |
1746 /* flt: transform with filters flt[0..maxlevel] */ | |
1747 /* method: transform with filter method */ | |
1748 /* cost: cost function for entropy computing */ | |
1749 /* epsilon: limit for threshold method */ | |
1750 /* -------------------------------------------------------------------- */ | |
1751 /* Description: whole decompositing down to maxlevel */ | |
1752 /* The original image must be in tree->image */ | |
1753 /************************************************************************/ | |
8 | 1754 int decompose_all(Image_tree tree,int maxlevel,FilterGH *flt,enum FilterMethod method, |
0 | 1755 enum Information_Cost cost,double epsilon) |
1756 { | |
1757 Image original,coarse,horizontal,vertical,diagonal; | |
1758 int e,width,height,level; | |
1759 | |
1760 if (tree->level<maxlevel) | |
1761 { | |
1762 tree->coarse=new_image_tree(); | |
1763 tree->horizontal=new_image_tree(); | |
1764 tree->vertical=new_image_tree(); | |
1765 tree->diagonal=new_image_tree(); | |
1766 | |
1767 original=tree->image; | |
1768 width=(original->width+1)/2; | |
1769 height=(original->height+1)/2; | |
1770 level=tree->level; | |
1771 | |
1772 coarse=new_image(width,height); | |
1773 horizontal=new_image(width,height); | |
1774 vertical=new_image(width,height); | |
1775 diagonal=new_image(width,height); | |
1776 if(!coarse||!vertical||!horizontal||!diagonal) goto error; | |
1777 | |
1778 | |
1779 e=decomposition(tree->image,coarse,horizontal,vertical,diagonal,flt[0]->g,flt[0]->h,method); | |
1780 if (!e) return 0; | |
1781 | |
1782 tree->coarse->image=coarse; | |
1783 tree->horizontal->image=horizontal; | |
1784 tree->vertical->image=vertical; | |
1785 tree->diagonal->image=diagonal; | |
1786 | |
1787 tree->coarse->entropy=compute_entropy(coarse,cost,epsilon); | |
1788 tree->horizontal->entropy=compute_entropy(horizontal,cost,epsilon); | |
1789 tree->vertical->entropy=compute_entropy(vertical,cost,epsilon); | |
1790 tree->diagonal->entropy=compute_entropy(diagonal,cost,epsilon); | |
1791 | |
1792 tree->coarse->level=tree->horizontal->level= | |
1793 tree->vertical->level=tree->diagonal->level=level+1; | |
1794 | |
1795 e=decompose_all(tree->coarse,maxlevel,flt+1,method,cost,epsilon); | |
1796 e=decompose_all(tree->horizontal,maxlevel,flt+1,method,cost,epsilon); | |
1797 e=decompose_all(tree->vertical,maxlevel,flt+1,method,cost,epsilon); | |
1798 e=decompose_all(tree->diagonal,maxlevel,flt+1,method,cost,epsilon); | |
1799 if (!e) return 0; | |
1800 | |
1801 } | |
1802 | |
1803 return 1; | |
1804 | |
1805 error: | |
1806 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
1807 return 0; | |
1808 } | |
1809 | |
1810 /************************************************************************/ | |
1811 /* Functionname: compute_levels */ | |
1812 /* -------------------------------------------------------------------- */ | |
1813 /* Parameter: */ | |
1814 /* tree: Image tree where the entropy should be computed */ | |
1815 /* entropies : array for entropy */ | |
1816 /* cost: carry best basis selection out with this costfunc */ | |
1817 /* epsilon: limit for threshold method */ | |
1818 /* -------------------------------------------------------------------- */ | |
1819 /* Description: Compute the entropies of all decomposition levels */ | |
1820 /************************************************************************/ | |
8 | 1821 static void compute_levels(Image_tree tree,double *entropies,enum Information_Cost cost,double epsilon) |
0 | 1822 { |
1823 if (tree->image){ | |
1824 entropies[tree->level]+=compute_entropy(tree->image,cost,epsilon); | |
1825 } | |
1826 if (tree->coarse) compute_levels(tree->coarse,entropies,cost,epsilon); | |
1827 if (tree->horizontal) compute_levels(tree->horizontal,entropies,cost,epsilon); | |
1828 if (tree->vertical) compute_levels(tree->vertical,entropies,cost,epsilon); | |
1829 if (tree->diagonal) compute_levels(tree->diagonal,entropies,cost,epsilon); | |
1830 | |
1831 } | |
1832 | |
1833 /************************************************************************/ | |
1834 /* Functionname: free_levels */ | |
1835 /* -------------------------------------------------------------------- */ | |
1836 /* Parameter: */ | |
1837 /* tree: Image tree which should be cleaned */ | |
1838 /* best: best level */ | |
1839 /* -------------------------------------------------------------------- */ | |
1840 /* Description: clean the image tree except the best level */ | |
1841 /************************************************************************/ | |
8 | 1842 static void free_levels(Image_tree tree,int best) |
0 | 1843 { |
1844 if (tree->level<best) | |
1845 { | |
1846 free_image(tree->image); | |
1847 tree->image=NULL; | |
1848 free_levels(tree->coarse,best); | |
1849 free_levels(tree->horizontal,best); | |
1850 free_levels(tree->vertical,best); | |
1851 free_levels(tree->diagonal,best); | |
1852 } | |
1853 else | |
1854 { | |
1855 if (tree->coarse) | |
1856 { | |
1857 free_image_tree(tree->coarse); | |
1858 free_image_tree(tree->horizontal); | |
1859 free_image_tree(tree->vertical); | |
1860 free_image_tree(tree->diagonal); | |
1861 tree->coarse=tree->horizontal=tree->vertical=tree->diagonal=NULL; | |
1862 } | |
1863 } | |
1864 } | |
1865 | |
1866 /************************************************************************/ | |
1867 /* Functionname: decompose_to_level */ | |
1868 /* -------------------------------------------------------------------- */ | |
1869 /* Parameter: */ | |
1870 /* original: original image */ | |
1871 /* level: decompose to level */ | |
1872 /* flt: decompos with filters[0..level] */ | |
1873 /* method: transform with filter method */ | |
1874 /* -------------------------------------------------------------------- */ | |
1875 /* Description: Decomposes an image to an certain level and stores */ | |
1876 /* only this level in the returned quadtree */ | |
1877 /************************************************************************/ | |
8 | 1878 Image_tree decompose_to_level(Image original,int level,FilterGH *flt,enum FilterMethod method) |
0 | 1879 { Image_tree tree; |
1880 int e; | |
1881 | |
1882 tree=new_image_tree(); | |
1883 tree->image=original; | |
1884 | |
1885 e=decompose_all(tree,level,flt,method,entropy,1); | |
1886 if (!e) return NULL; | |
1887 | |
1888 free_levels(tree,level); | |
1889 | |
1890 return tree; | |
1891 | |
1892 } | |
1893 | |
1894 /************************************************************************/ | |
1895 /* Functionname: decomposition */ | |
1896 /* -------------------------------------------------------------------- */ | |
1897 /* Parameter: */ | |
1898 /* t_img: Image which should be decomposed */ | |
1899 /* coarse,horizontal,vertical,diagonal: */ | |
1900 /* decomposed images */ | |
1901 /* method: transform with filter method */ | |
1902 /* g,h: the transformation is carried out with these filters*/ | |
1903 /* -------------------------------------------------------------------- */ | |
1904 /* Description: This carries out one wavelettransformation */ | |
1905 /* using waveletfilters. */ | |
1906 /************************************************************************/ | |
1907 | |
1908 static int decomposition(Image t_img,Image coarse,Image horizontal,Image vertical, | |
1909 Image diagonal,Filter g,Filter h,enum FilterMethod method) | |
1910 { Image temp1; | |
1911 | |
1912 /*coarse*/ | |
1913 temp1=new_image(coarse->width,t_img->height); | |
1914 if(!temp1) goto error; | |
1915 convolute_lines(temp1,t_img,h,method); | |
1916 convolute_rows(coarse,temp1,h,method); | |
1917 | |
1918 /*horizontal*/ | |
1919 convolute_rows(horizontal,temp1,g,method); | |
1920 free_image(temp1); | |
1921 | |
1922 /*vertical*/ | |
1923 temp1=new_image(vertical->width,t_img->height); | |
1924 if(!temp1) goto error; | |
1925 convolute_lines(temp1,t_img,g,method); | |
1926 convolute_rows(vertical,temp1,h,method); | |
1927 | |
1928 /*diagonal*/ | |
1929 convolute_rows(diagonal,temp1,g,method); | |
1930 free_image(temp1); | |
1931 | |
1932 return 1; | |
1933 | |
1934 error: | |
1935 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
1936 return 0; | |
1937 | |
1938 } | |
1939 | |
1940 /************************************************************************/ | |
1941 /* Functionname: inv_decomposition */ | |
1942 /* -------------------------------------------------------------------- */ | |
1943 /* Parameter: */ | |
1944 /* sum: reconstructed image */ | |
1945 /* coarse,horizontal,vertical,diagonal: images to carry out*/ | |
1946 /* the inverse transformation */ | |
1947 /* flt_gh: transform with filters g and h */ | |
1948 /* method: transform with filter method */ | |
1949 /* -------------------------------------------------------------------- */ | |
1950 /* Description: Carries out the wavelettransform */ | |
1951 /* */ | |
1952 /************************************************************************/ | |
1953 static int inv_decomposition(Image sum,Image coarse,Image horizontal,Image vertical, | |
1954 Image diagonal,FilterGH flt_gh,enum FilterMethod method) | |
1955 { Image temp1; | |
1956 Filter g,h; | |
1957 | |
1958 if (flt_gh->type==FTOrtho) { | |
1959 g=flt_gh->g; | |
1960 h=flt_gh->h; | |
1961 } | |
1962 else { | |
1963 g=flt_gh->gi; | |
1964 h=flt_gh->hi; | |
1965 } | |
1966 | |
1967 /*coarse*/ | |
1968 temp1=new_image(coarse->width,sum->height); | |
1969 if(!temp1) goto error; | |
1970 convolute_rows(temp1,coarse,h,method); | |
1971 | |
1972 /*horizontal*/ | |
1973 convolute_rows(temp1,horizontal,g,method); | |
1974 convolute_lines(sum,temp1,h,method); | |
1975 free_image(temp1); | |
1976 | |
1977 /*vertical*/ | |
1978 temp1=new_image(vertical->width,sum->height); | |
1979 if(!temp1) goto error; | |
1980 convolute_rows(temp1,vertical,h,method); | |
1981 | |
1982 /*diagonal*/ | |
1983 convolute_rows(temp1,diagonal,g,method); | |
1984 convolute_lines(sum,temp1,g,method); | |
1985 | |
1986 free_image(temp1); | |
1987 | |
1988 return 1; | |
1989 | |
1990 error: | |
1991 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
1992 return 0; | |
1993 } | |
1994 | |
1995 /************************************************************************/ | |
1996 /* Functionname: build_image */ | |
1997 /* -------------------------------------------------------------------- */ | |
1998 /* Parameter: */ | |
1999 /* quadtree: quadtree with decomposition information */ | |
2000 /* width,height: image width and height */ | |
2001 /* RETURN: returns the build up image */ | |
2002 /* -------------------------------------------------------------------- */ | |
2003 /* Description: builds up an image out of an Image_tree */ | |
2004 /* */ | |
2005 /************************************************************************/ | |
8 | 2006 Image build_image(Image_tree quadtree,int width,int height) |
0 | 2007 { Image ret_img,coarse,horizontal,vertical,diagonal; |
2008 | |
2009 | |
2010 ret_img=new_image(width,height); | |
2011 if(!ret_img) goto error; | |
2012 | |
2013 width=(width+1)/2; | |
2014 height=(height+1)/2; | |
2015 | |
2016 if (!(quadtree->image)) { | |
2017 coarse=build_image(quadtree->coarse,width,height); | |
2018 horizontal=build_image(quadtree->horizontal,width,height); | |
2019 vertical=build_image(quadtree->vertical,width,height); | |
2020 diagonal=build_image(quadtree->diagonal,width,height); | |
2021 if (!coarse||!horizontal||!vertical||!diagonal) return NULL; | |
2022 | |
2023 copy_into_image(ret_img,coarse,0,0); | |
2024 copy_into_image(ret_img,horizontal,width,0); | |
2025 copy_into_image(ret_img,vertical,0,height); | |
2026 copy_into_image(ret_img,diagonal,width,height); | |
2027 | |
2028 if (!quadtree->coarse->image) free_image(coarse); | |
2029 if (!quadtree->horizontal->image) free_image(horizontal); | |
2030 if (!quadtree->vertical->image) free_image(vertical); | |
2031 if (!quadtree->diagonal->image) free_image(diagonal); | |
2032 | |
2033 return ret_img; | |
2034 } | |
2035 else return quadtree->image; | |
2036 | |
2037 error: | |
2038 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
2039 return NULL; | |
2040 } | |
2041 | |
2042 /************************************************************************/ | |
2043 /* Functionname: inv_transform */ | |
2044 /* -------------------------------------------------------------------- */ | |
2045 /* Parameter: */ | |
2046 /* tree: tree with decomposition information */ | |
2047 /* flt_gh: transform with filters g and h */ | |
2048 /* method: transform with filter method */ | |
2049 /* -------------------------------------------------------------------- */ | |
2050 /* Description: Inverts the wavelettransform,best_basis,best_level */ | |
2051 /* */ | |
2052 /************************************************************************/ | |
8 | 2053 Image inv_transform(Image_tree tree,FilterGH *flt, |
0 | 2054 enum FilterMethod method) |
2055 | |
2056 { int er,width,height; | |
2057 Image ret_img,coarse,vertical,horizontal,diagonal; | |
2058 | |
2059 if (!tree->image) { | |
2060 | |
2061 coarse=inv_transform(tree->coarse,flt,method); | |
2062 horizontal=inv_transform(tree->horizontal,flt,method); | |
2063 vertical=inv_transform(tree->vertical,flt,method); | |
2064 diagonal=inv_transform(tree->diagonal,flt,method); | |
2065 if (!coarse||!horizontal||!vertical||!diagonal) return NULL; | |
2066 | |
2067 width=coarse->width+horizontal->width; | |
2068 height=coarse->height+vertical->height; | |
2069 | |
2070 ret_img=new_image(width,height); | |
2071 if(!ret_img) goto error; | |
2072 | |
2073 | |
2074 if (tree->flag==0) /*if flag is set it is a doubletree tiling*/ | |
2075 { | |
2076 // er=inv_decomposition(ret_img,coarse,horizontal,vertical,diagonal,flt[1],method); | |
2077 er=inv_decomposition(ret_img,coarse,horizontal,vertical,diagonal,flt[tree->level],method); | |
2078 if (!er) return NULL; | |
2079 } | |
2080 else | |
2081 { | |
2082 copy_into_image(ret_img,coarse,0,0); | |
2083 copy_into_image(ret_img,horizontal,coarse->width,0); | |
2084 copy_into_image(ret_img,vertical,0,coarse->height); | |
2085 copy_into_image(ret_img,diagonal,coarse->width,coarse->height); | |
2086 } | |
2087 | |
2088 if (!tree->coarse->image) free_image(coarse); | |
2089 if (!tree->horizontal->image) free_image(horizontal); | |
2090 if (!tree->vertical->image) free_image(vertical); | |
2091 if (!tree->diagonal->image) free_image(diagonal); | |
2092 | |
2093 return ret_img; | |
2094 } | |
2095 | |
2096 else return tree->image; | |
2097 | |
2098 error: | |
2099 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
2100 return NULL; | |
2101 } | |
2102 | |
2103 /************************************************************************/ | |
2104 /* Functionname: find_deepest_level */ | |
2105 /* -------------------------------------------------------------------- */ | |
2106 /* Parameter: */ | |
2107 /* -------------------------------------------------------------------- */ | |
2108 /* Description: Finds the deepest possible level where width and */ | |
2109 /* height can divided by two exactly. */ | |
2110 /************************************************************************/ | |
8 | 2111 int find_deepest_level(int width,int height) |
0 | 2112 { |
2113 int level=0,w=width,h=height; | |
2114 | |
2115 while ( !((w%2)||(h%2))) | |
2116 { | |
2117 w=w/2; | |
2118 h=h/2; | |
2119 level++; | |
2120 } | |
2121 | |
2122 return level-1; | |
2123 | |
2124 } | |
2125 | |
2126 /************************************************************************/ | |
2127 /* Functionname: convolute_lines */ | |
2128 /* -------------------------------------------------------------------- */ | |
2129 /* Parameter: */ | |
2130 /* output: output image of wavelettransformation */ | |
2131 /* input: input image for decomposition */ | |
2132 /* flt: transform with filter flt */ | |
2133 /* method: transform with filter method */ | |
2134 /* -------------------------------------------------------------------- */ | |
2135 /* Description: Carries out the wavelettransform for all lines of */ | |
2136 /* the input image */ | |
2137 /************************************************************************/ | |
2138 static int convolute_lines(Image output,Image input,Filter flt,enum FilterMethod method) | |
2139 /*Convolute the lines with filter*/ | |
2140 { int i; | |
2141 | |
2142 for (i=0;i<input->height;i++) { | |
2143 switch(method) { | |
2144 case cutoff: | |
2145 filter_cutoff(input,input->width*i,input->width,1, | |
2146 output,output->width*i,output->width,1,flt); | |
2147 break; | |
2148 case inv_cutoff: | |
2149 filter_inv_cutoff(input,input->width*i,input->width,1, | |
2150 output,output->width*i,output->width,1,flt); | |
2151 break; | |
2152 case periodical: | |
2153 filter_periodical(input,input->width*i,input->width,1, | |
2154 output,output->width*i,output->width,1,flt); | |
2155 break; | |
2156 case inv_periodical: | |
2157 filter_inv_periodical(input,input->width*i,input->width,1, | |
2158 output,output->width*i,output->width,1,flt); | |
2159 break; | |
2160 case mirror: | |
2161 filter_mirror(input,input->width*i,input->width,1, | |
2162 output,output->width*i,output->width,1,flt); | |
2163 break; | |
2164 case inv_mirror: | |
2165 filter_inv_mirror(input,input->width*i,input->width,1, | |
2166 output,output->width*i,output->width,1,flt); | |
2167 break; | |
2168 | |
2169 | |
2170 } | |
2171 } | |
2172 | |
2173 return 1; | |
2174 } | |
2175 | |
2176 /************************************************************************/ | |
2177 /* Functionname: convolute_rows */ | |
2178 /* -------------------------------------------------------------------- */ | |
2179 /* Parameter: */ | |
2180 /* output: output image of wavelettransformation */ | |
2181 /* input: input image for decomposition */ | |
2182 /* flt: transform with filter flt */ | |
2183 /* method: transform with filter method */ | |
2184 /* -------------------------------------------------------------------- */ | |
2185 /* Description: Carries out the wavelettransform for all rows of */ | |
2186 /* the input image */ | |
2187 /************************************************************************/ | |
2188 static int convolute_rows(Image output,Image input,Filter flt,enum FilterMethod method) | |
2189 /*Convolute the rows with filter*/ | |
2190 { int i; | |
2191 | |
2192 for (i=0;i<input->width;i++) | |
2193 { | |
2194 switch (method) | |
2195 { | |
2196 case cutoff: | |
2197 filter_cutoff(input,i,input->height,input->width, | |
2198 output,i,output->height,output->width,flt); | |
2199 break; | |
2200 case inv_cutoff: | |
2201 filter_inv_cutoff(input,i,input->height,input->width, | |
2202 output,i,output->height,output->width,flt); | |
2203 break; | |
2204 case periodical: | |
2205 filter_periodical(input,i,input->height,input->width, | |
2206 output,i,output->height,output->width,flt); | |
2207 break; | |
2208 case inv_periodical: | |
2209 filter_inv_periodical(input,i,input->height,input->width, | |
2210 output,i,output->height,output->width,flt); | |
2211 break; | |
2212 case mirror: | |
2213 filter_mirror(input,i,input->height,input->width, | |
2214 output,i,output->height,output->width,flt); | |
2215 break; | |
2216 case inv_mirror: | |
2217 filter_inv_mirror(input,i,input->height,input->width, | |
2218 output,i,output->height,output->width,flt); | |
2219 break; | |
2220 | |
2221 } | |
2222 } | |
2223 return 1; | |
2224 } | |
2225 | |
2226 /************************************************************************/ | |
2227 /* Functionname: sumationq */ | |
2228 /* -------------------------------------------------------------------- */ | |
2229 /* Parameter: */ | |
2230 /* img: image to compute */ | |
2231 /* -------------------------------------------------------------------- */ | |
2232 /* Description: compute the sum of quadrats of all elements of */ | |
2233 /* the input image */ | |
2234 /************************************************************************/ | |
2235 static Pixel sumationq(Image img) | |
2236 { Pixel sum=0; | |
2237 int i; | |
2238 | |
2239 for (i=0;i<img->size;i++) { | |
2240 sum+=(*img->data+i)*(*img->data+i); | |
2241 } | |
2242 return sum; | |
2243 } | |
2244 | |
2245 /************************************************************************/ | |
2246 /* Functionname: normq */ | |
2247 /* -------------------------------------------------------------------- */ | |
2248 /* Parameter: */ | |
2249 /* tree: tree to compute */ | |
2250 /* -------------------------------------------------------------------- */ | |
2251 /* Description: computes the quadratic norm over all images in */ | |
2252 /* the input tree */ | |
2253 /************************************************************************/ | |
2254 static Pixel normq(Image_tree tree) | |
2255 { Pixel sum=0; | |
2256 | |
2257 if (tree->image) | |
2258 { | |
2259 sum=sumationq(tree->image); | |
2260 } | |
2261 else | |
2262 { | |
2263 if (tree->coarse) sum+=normq(tree->coarse); | |
2264 if (tree->horizontal) sum+=normq(tree->horizontal); | |
2265 if (tree->vertical) sum+=normq(tree->vertical); | |
2266 if (tree->diagonal) sum+=normq(tree->diagonal); | |
2267 } | |
2268 | |
2269 return sum; | |
2270 } | |
2271 | |
2272 /************************************************************************/ | |
2273 /* Functionname: sumation_down */ | |
2274 /* -------------------------------------------------------------------- */ | |
2275 /* Parameter: */ | |
2276 /* tree: tree to compute */ | |
2277 /* normq: norm of the images in the tree */ | |
2278 /* -------------------------------------------------------------------- */ | |
2279 /* Description: computes the Entropy over all (string aded) images */ | |
2280 /* in the input tree */ | |
2281 /************************************************************************/ | |
2282 static Pixel sumation_down(Image_tree tree, Pixel normq) | |
2283 { Pixel sum=0,p; | |
2284 int i; | |
2285 Image img; | |
2286 Pixel *data; | |
2287 | |
2288 if (tree->image) | |
2289 { | |
2290 img=tree->image; | |
2291 data=img->data; | |
2292 for (i=0;i<img->size;i++,data++) | |
2293 { | |
2294 if (*data!=0) | |
2295 { | |
2296 p=(*data)*(*data)/normq; | |
2297 sum+=p*log(1/p); | |
2298 } | |
2299 } | |
2300 } | |
2301 else | |
2302 { | |
2303 if (tree->coarse) sum+=sumation_down(tree->coarse,normq); | |
2304 if (tree->horizontal) sum+=sumation_down(tree->horizontal,normq); | |
2305 if (tree->vertical) sum+=sumation_down(tree->vertical,normq); | |
2306 if (tree->diagonal) sum+=sumation_down(tree->diagonal,normq); | |
2307 } | |
2308 | |
2309 return sum; | |
2310 } | |
2311 | |
2312 /************************************************************************/ | |
2313 /* Functionname: comp */ | |
2314 /* -------------------------------------------------------------------- */ | |
2315 /* Description: used for quicksort for decreasing order */ | |
2316 /************************************************************************/ | |
2317 int comp(const Pixel *x,const Pixel *y) | |
2318 { | |
2319 if (*x<*y) return 1; | |
2320 else if (*x==*y) return 0; | |
2321 else return -1; | |
2322 } | |
2323 | |
2324 /************************************************************************/ | |
2325 /* Functionname: recarea */ | |
2326 /* tree: Image tree to compute */ | |
2327 /* list: target list */ | |
2328 /* list_size: actual size of the list */ | |
2329 /* -------------------------------------------------------------------- */ | |
2330 /* Description: copies all elements within the tree into an list */ | |
2331 /************************************************************************/ | |
2332 static void recarea(Image_tree tree,Pixel *list,int *list_size) | |
2333 { Image img; | |
2334 | |
2335 if (tree->image) | |
2336 { | |
2337 img=tree->image; | |
2338 memcpy(list+(*list_size),img->data,img->size*sizeof(Pixel)); | |
2339 *list_size+=img->size; | |
2340 } | |
2341 else | |
2342 { | |
2343 if (tree->coarse) recarea(tree->coarse,list,list_size); | |
2344 if (tree->horizontal) recarea(tree->horizontal,list,list_size); | |
2345 if (tree->vertical) recarea(tree->vertical,list,list_size); | |
2346 if (tree->diagonal) recarea(tree->diagonal,list,list_size); | |
2347 } | |
2348 | |
2349 } | |
2350 | |
2351 static void abs_list(Pixel *list,int list_size) | |
2352 { | |
2353 int i; | |
2354 | |
2355 for (i=0;i<list_size;i++) list[i]=fabs(list[i]); | |
2356 } | |
2357 | |
2358 /************************************************************************/ | |
2359 /* Functionname: sum_list */ | |
2360 /* -------------------------------------------------------------------- */ | |
2361 /* Description: computes the sum of all poweres list elements */ | |
2362 /************************************************************************/ | |
2363 static Pixel sum_list(Pixel *list,int p,int size) | |
2364 { Pixel sum=0; | |
2365 int i; | |
2366 | |
2367 for (i=0;i<size;i++) | |
2368 { | |
2369 if (p!=1) sum+=pow(list[i],p); | |
2370 else sum+=list[i]; | |
2371 } | |
2372 return sum; | |
2373 } | |
2374 | |
2375 static Pixel weak_lp(Image_tree tree,int size,int p,double epsilon) | |
2376 { Pixel wlp,*list,max=0; | |
2377 int *list_size,k; | |
2378 | |
2379 list_size=(int *)malloc(sizeof(int)); | |
2380 if (!list_size) goto error; | |
2381 | |
2382 *list_size=0; | |
2383 | |
2384 list=(Pixel *)calloc(size,sizeof(Pixel)); | |
2385 if (!list) goto error; | |
2386 | |
2387 recarea(tree,list,list_size); | |
2388 abs_list(list,*list_size); | |
2389 | |
8 | 2390 qsort(list,*list_size, sizeof(Pixel), (int (*)(const void*, const void*)) comp); |
0 | 2391 |
2392 for (k=0;k<size;k++) | |
2393 { | |
2394 if (k!=0) wlp=pow(k,1/p)*list[k]; | |
2395 else wlp=0; | |
2396 if (wlp>max) max=wlp; | |
2397 } | |
2398 | |
2399 free(list); | |
2400 | |
2401 return max; | |
2402 | |
2403 error: | |
2404 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
2405 return 0; | |
2406 } | |
2407 | |
2408 static Pixel comp_number(Image_tree tree,int size,int p,double f) | |
8 | 2409 { Pixel sum=0,*list,min=MAXDOUBLE,norm,npf,normf; |
0 | 2410 int *list_size=0,k; |
2411 | |
2412 list_size=(int *)malloc(sizeof(int)); | |
2413 if (!list_size) goto error; | |
2414 *list_size=0; | |
2415 | |
2416 list=(Pixel *)calloc(size,sizeof(Pixel)); | |
2417 if (!list) goto error; | |
2418 recarea(tree,list,list_size); | |
2419 abs_list(list,*list_size); | |
2420 | |
8 | 2421 qsort(list,*list_size, sizeof(Pixel), (int (*)(const void*, const void*)) comp); |
0 | 2422 |
2423 norm=sum_list(list,p,size); | |
2424 normf=norm*f; | |
2425 | |
2426 for (k=0;k<size;k++) | |
2427 { | |
2428 if (list[k]!=0) | |
2429 { | |
2430 sum+=pow(list[k],p); | |
2431 npf=fabs(sum-normf); | |
2432 if (npf<min) min=npf; | |
2433 } | |
2434 } | |
2435 min=min/norm; | |
2436 | |
2437 free(list); | |
2438 | |
2439 return min; | |
2440 | |
2441 error: | |
2442 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
2443 return 0; | |
2444 } | |
2445 | |
2446 static Pixel comp_area(Image_tree tree,int size,int p,double f) | |
2447 { Pixel sum=0,*list,norm,vk=0; | |
2448 int *list_size=0,k; | |
2449 | |
2450 list_size=(int *)malloc(sizeof(int)); | |
2451 if (!list_size) goto error; | |
2452 *list_size=0; | |
2453 | |
2454 list=(Pixel *)calloc(size,sizeof(Pixel)); | |
2455 if (!list) goto error; | |
2456 | |
2457 recarea(tree,list,list_size); | |
2458 abs_list(list,*list_size); | |
2459 | |
8 | 2460 qsort(list,*list_size, sizeof(Pixel), (int (*)(const void*, const void*)) comp); |
0 | 2461 |
2462 norm=sum_list(list,p,size); | |
2463 | |
2464 for (k=0;k<size;k++) | |
2465 { | |
2466 if (list[k]!=0) | |
2467 { | |
2468 vk+=pow(list[k],p); | |
2469 sum+=vk; | |
2470 | |
2471 } | |
2472 } | |
2473 | |
2474 free(list); | |
2475 | |
2476 return (size-sum/norm); | |
2477 | |
2478 error: | |
2479 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
2480 return 0; | |
2481 } | |
2482 | |
2483 static Pixel compute_sdiscrepancy(Image_tree tree,int size) | |
2484 { Pixel *list,min,max,factor,maximum=0,x; | |
2485 int *list_size=0,k; | |
2486 | |
2487 list_size=(int *)malloc(sizeof(int)); | |
2488 if (!list_size) goto error; | |
2489 *list_size=0; | |
2490 | |
2491 list=(Pixel *)calloc(size,sizeof(Pixel)); | |
2492 if (!list) goto error; | |
2493 | |
2494 recarea(tree,list,list_size); | |
2495 | |
8 | 2496 qsort(list,*list_size, sizeof(Pixel), (int (*)(const void*, const void*)) comp); |
0 | 2497 |
2498 min=list[0]; | |
2499 max=list[size-1]; | |
2500 factor=1/(max-min); | |
2501 | |
2502 /*scaling to [0,1]*/ | |
2503 for (k=0;k<size;k++) | |
2504 { | |
2505 list[k]=factor*(list[k]-min); | |
2506 } | |
2507 | |
2508 for (k=0;k<size;k++) | |
2509 { | |
2510 x=fabs(list[k]-(2*k-1)/(2*size)); | |
2511 if (x>maximum) maximum=x; | |
2512 } | |
2513 | |
2514 free(list); | |
2515 | |
2516 return (1/(2*size)+maximum); | |
2517 | |
2518 error: | |
2519 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
2520 return 0; | |
2521 } | |
2522 | |
2523 static Pixel compute_discrepancy(Image_tree tree,int size) | |
2524 { Pixel *list,min,max,factor,maximum=0,minimum=0,x; | |
2525 int *list_size=0,k; | |
2526 | |
2527 list_size=(int *)malloc(sizeof(int)); | |
2528 if (!list_size) goto error; | |
2529 *list_size=0; | |
2530 | |
2531 list=(Pixel *)calloc(size,sizeof(Pixel)); | |
2532 if (!list) goto error; | |
2533 | |
2534 recarea(tree,list,list_size); | |
2535 | |
8 | 2536 qsort(list,*list_size, sizeof(Pixel), (int (*)(const void*, const void*)) comp); |
0 | 2537 |
2538 min=list[0]; | |
2539 max=list[size-1]; | |
2540 factor=1/(max-min); | |
2541 | |
2542 /*scaling to [0,1]*/ | |
2543 for (k=0;k<size;k++) | |
2544 { | |
2545 list[k]=factor*(list[k]-min); | |
2546 } | |
2547 | |
2548 for (k=0;k<size;k++) | |
2549 { | |
2550 x=((Pixel)k/size-list[k]); | |
2551 if (x>maximum) maximum=x; | |
2552 else if (x<minimum) minimum=x; | |
2553 | |
2554 } | |
2555 | |
2556 free(list); | |
2557 | |
2558 return (1/size+maximum-minimum); | |
2559 | |
2560 error: | |
2561 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
2562 return 0; | |
2563 } | |
2564 | |
2565 static Pixel compute_concentration(Image_tree tree,int size) | |
2566 { Pixel *list,min,max,factor,lkm=0,length,sum=0,value,norm; | |
2567 int *list_size=0,k,next=0,last=0,i=0; | |
2568 | |
2569 list_size=(int *)malloc(sizeof(int)); | |
2570 if (!list_size) goto error; | |
2571 *list_size=0; | |
2572 | |
2573 list=(Pixel *)calloc(size+1,sizeof(Pixel)); | |
2574 if (!list) goto error; | |
2575 | |
2576 recarea(tree,list,list_size); | |
2577 | |
8 | 2578 qsort(list,*list_size, sizeof(Pixel), (int (*)(const void*, const void*)) comp); |
0 | 2579 |
2580 min=list[0]; | |
2581 max=list[size-1]; | |
2582 length=(max-min)/100; | |
2583 | |
2584 factor=1/(max-min); | |
2585 for (k=0;k<size;k++) | |
2586 { | |
2587 list[k]=factor*(list[k]-min); | |
2588 } | |
2589 | |
2590 norm=size*sum_list(list,1,size); | |
2591 length=0.01; | |
2592 value=length; | |
2593 | |
2594 list[size]=max+value; | |
2595 | |
2596 for (k=0;k<100;k++) | |
2597 { | |
2598 while ((list[i]<value)&(i<size)) | |
2599 { | |
2600 sum+=list[i]; | |
2601 next++; | |
2602 i++; | |
2603 } | |
2604 lkm+=(next-last)*sum/norm; | |
2605 value+=length; | |
2606 last=next; | |
2607 sum=0; | |
2608 } | |
2609 | |
2610 return -lkm; | |
2611 | |
2612 error: | |
2613 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
2614 return 0; | |
2615 } | |
2616 /************************************************************************/ | |
2617 /* Functionname: compute_entropy */ | |
2618 /* -------------------------------------------------------------------- */ | |
2619 /* Parameter: */ | |
2620 /* img: Image from which the entropy should be computed */ | |
2621 /* cost: choosen costfunction */ | |
2622 /* epsilon: limit for threshold method */ | |
2623 /* -------------------------------------------------------------------- */ | |
2624 /* Description: computes entropy of an image */ | |
2625 /************************************************************************/ | |
2626 static double compute_entropy(Image img,enum Information_Cost cost,double epsilon) | |
2627 { double sum=0,x=0; | |
2628 int i; | |
2629 Pixel *data; | |
2630 | |
2631 data=img->data; | |
2632 | |
2633 switch(cost) { | |
2634 | |
2635 case threshold: | |
2636 for(i=0;i<img->size;i++) | |
2637 if (fabs(img->data[i])>epsilon) sum++; | |
2638 break; | |
2639 | |
2640 case log_energy: | |
2641 for(i=0;i<img->size;i++,data++) { | |
2642 x=(*data) * (*data); | |
2643 if (x!=0) sum+=(x*log(1/x)); | |
2644 } | |
2645 break; | |
2646 | |
2647 case norml: | |
2648 for(i=0;i<img->size;i++,data++) { | |
2649 x=fabs(*data); | |
2650 sum+=x; | |
2651 } | |
2652 break; | |
2653 | |
2654 case norml2: | |
2655 for(i=0;i<img->size;i++,data++) { | |
2656 x=(*data) * (*data); | |
2657 sum+=x; | |
2658 } | |
2659 sum=pow(sum,0.5); | |
2660 break; | |
2661 | |
2662 case entropy: | |
2663 for(i=0;i<img->size;i++,data++) { | |
2664 x=(*data)*(*data); | |
2665 if (x!=0) sum-=(x*log(x)); | |
2666 } | |
2667 break; | |
2668 | |
2669 case gauss_markov: | |
2670 for(i=0;i<img->size;i++) { | |
2671 x=(img->data[i])*(img->data[i]); | |
2672 if (x!=0) sum+=log(x*x); | |
2673 } | |
2674 break; | |
2675 | |
2676 } | |
2677 | |
2678 return sum; | |
2679 } | |
2680 | |
2681 /************************************************************************/ | |
2682 /* Functionname: compute_non_additive */ | |
2683 /* -------------------------------------------------------------------- */ | |
2684 /* Parameter: */ | |
2685 /* tree: Image tree from which the entropy should be */ | |
2686 /* computed */ | |
2687 /* size : size of the image */ | |
2688 /* cost: choosen costfunction */ | |
2689 /* epsilon: limit for threshold method */ | |
2690 /* down: decides if only the first image should be computed*/ | |
2691 /* -------------------------------------------------------------------- */ | |
2692 /* Description: computes entropy of an image */ | |
2693 /************************************************************************/ | |
2694 static Pixel compute_non_additive(Image_tree tree,int size,enum Information_Cost cost,double epsilon,int down) | |
2695 { Pixel sum=0,normx; | |
8 | 2696 Image img=NULL; |
0 | 2697 |
2698 if (down) | |
2699 { | |
2700 img=tree->image; | |
2701 tree->image=NULL; | |
2702 } | |
2703 switch (cost) | |
2704 { | |
2705 case shanon: | |
2706 normx=normq(tree); | |
2707 sum=-sumation_down(tree,normx); | |
2708 | |
2709 break; | |
2710 case weak_l: | |
2711 sum=weak_lp(tree,size,1,epsilon); | |
2712 break; | |
2713 case weak_lq: | |
2714 sum=weak_lp(tree,size,2,epsilon); | |
2715 break; | |
2716 case compression_number: | |
2717 sum=comp_number(tree,size,1,epsilon); | |
2718 break; | |
2719 case compression_numberq: | |
2720 sum=comp_number(tree,size,2,epsilon); | |
2721 break; | |
2722 case compression_area: | |
2723 sum=comp_area(tree,size,1,epsilon); | |
2724 break; | |
2725 case compression_areaq: | |
2726 sum=comp_area(tree,size,2,epsilon); | |
2727 break; | |
2728 case discrepancy: | |
2729 sum=compute_discrepancy(tree,size); | |
2730 break; | |
2731 case sdiscrepancy: | |
2732 sum=compute_sdiscrepancy(tree,size); | |
2733 break; | |
2734 case concentration: | |
2735 sum=compute_concentration(tree,size); | |
2736 break; | |
2737 | |
2738 | |
2739 } | |
2740 | |
2741 if (down) tree->image=img; | |
2742 | |
2743 return sum; | |
2744 } | |
2745 | |
8 | 2746 int rec_double(Image_tree dtree,int level,FilterGH *flt,enum FilterMethod method,enum Information_Cost cost,double epsilon) |
0 | 2747 |
2748 { int min,width,height; | |
2749 double sum=0; | |
2750 Image c,h,v,d; | |
2751 | |
2752 dtree->level=0; | |
2753 | |
2754 if (cost>=shanon) | |
2755 { | |
2756 dtree->entropy=compute_non_additive(dtree,dtree->image->size,cost,epsilon,0); | |
2757 } | |
2758 else dtree->entropy=compute_entropy(dtree->image,cost,epsilon); | |
2759 | |
2760 dtree->doubletree=best_basis(dtree->image,level,flt,method,cost,epsilon); | |
2761 | |
2762 min=dtree->image->width; | |
2763 if (dtree->image->height<min) min=dtree->image->height; | |
2764 | |
2765 if (doubletree_min<min) | |
2766 { | |
2767 width=(dtree->image->width+1)/2; | |
2768 height=(dtree->image->height+1)/2; | |
2769 | |
2770 dtree->coarse=new_image_tree(); | |
2771 dtree->horizontal=new_image_tree(); | |
2772 dtree->vertical=new_image_tree(); | |
2773 dtree->diagonal=new_image_tree(); | |
2774 | |
2775 c=new_image(width,height); | |
2776 h=new_image(width,height); | |
2777 v=new_image(width,height); | |
2778 d=new_image(width,height); | |
2779 if(!c||!h||!v||!d) goto error; | |
2780 | |
2781 | |
2782 copy_part_of_image(c,dtree->image,0,0); | |
2783 copy_part_of_image(h,dtree->image,width,0); | |
2784 copy_part_of_image(v,dtree->image,0,height); | |
2785 copy_part_of_image(d,dtree->image,width,height); | |
2786 | |
2787 dtree->coarse->image=c; | |
2788 dtree->horizontal->image=h; | |
2789 dtree->vertical->image=v; | |
2790 dtree->diagonal->image=d; | |
2791 | |
2792 rec_double(dtree->coarse,level,flt,method,cost,epsilon); | |
2793 rec_double(dtree->horizontal,level,flt,method,cost,epsilon); | |
2794 rec_double(dtree->vertical,level,flt,method,cost,epsilon); | |
2795 rec_double(dtree->diagonal,level,flt,method,cost,epsilon); | |
2796 | |
2797 /* going back in recursion*/ | |
2798 | |
2799 sum=dtree->coarse->entropy+dtree->horizontal->entropy+ | |
2800 dtree->vertical->entropy+dtree->diagonal->entropy; | |
2801 | |
2802 if (sum>dtree->entropy) | |
2803 { | |
2804 /*take image*/ | |
2805 | |
2806 free_image_tree(dtree->coarse); | |
2807 free_image_tree(dtree->horizontal); | |
2808 free_image_tree(dtree->vertical); | |
2809 free_image_tree(dtree->diagonal); | |
2810 dtree->coarse=dtree->horizontal=dtree->vertical=dtree->diagonal=NULL; | |
2811 } | |
2812 else | |
2813 { /*take tiling*/ | |
2814 dtree->entropy=sum; | |
2815 free_image(dtree->image); | |
2816 dtree->image=NULL; | |
2817 } | |
2818 | |
2819 if (dtree->entropy>dtree->doubletree->entropy) | |
2820 { | |
2821 /*take best basis tree*/ | |
2822 | |
2823 dtree->entropy=dtree->doubletree->entropy; | |
2824 | |
2825 if(dtree->coarse) free_image_tree(dtree->coarse); | |
2826 if(dtree->horizontal) free_image_tree(dtree->horizontal); | |
2827 if(dtree->vertical) free_image_tree(dtree->vertical); | |
2828 if(dtree->diagonal) free_image_tree(dtree->diagonal); | |
2829 | |
2830 dtree->coarse=dtree->doubletree->coarse; | |
2831 dtree->horizontal=dtree->doubletree->horizontal; | |
2832 dtree->vertical=dtree->doubletree->vertical; | |
2833 dtree->diagonal=dtree->doubletree->diagonal; | |
2834 | |
2835 free_image(dtree->image); | |
2836 dtree->image=NULL; | |
2837 free(dtree->doubletree); | |
2838 dtree->doubletree=NULL; | |
2839 | |
2840 } | |
2841 else | |
2842 { | |
2843 dtree->flag=1; | |
2844 if(dtree->doubletree) free_image_tree(dtree->doubletree); | |
2845 dtree->doubletree=NULL; | |
2846 } | |
2847 } | |
2848 | |
2849 return 1; | |
2850 | |
2851 error: | |
2852 err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory)); | |
2853 return 0; | |
2854 } | |
2855 | |
8 | 2856 static void save_structur(Image_tree tree,FILE *fp,int pos) |
0 | 2857 { |
2858 int shift,next_pos,max; | |
2859 | |
2860 if (tree->flag) | |
2861 { | |
2862 fprintf(fp,"%d ",pos); | |
2863 | |
2864 shift=pos-(pow(4,tree->level-1)-1)*4/3-1; | |
2865 max=(int) ((pow(4,tree->level)-1)*4/3); | |
2866 next_pos=max+4*shift+1; | |
2867 if (tree->coarse) save_structur(tree->coarse,fp,next_pos); | |
2868 if (tree->horizontal) save_structur(tree->horizontal,fp,next_pos+1); | |
2869 if (tree->vertical) save_structur(tree->vertical,fp,next_pos+2); | |
2870 if (tree->diagonal) save_structur(tree->diagonal,fp,next_pos+3); | |
2871 } | |
2872 } | |
2873 | |
2874 static int is_in_list(int *list,int len, int x) | |
2875 { | |
2876 int i,found=0; | |
2877 | |
2878 for (i=0;i<len;i++) | |
2879 { | |
2880 if (list[i]==x) | |
2881 { | |
2882 found=1; | |
2883 i=len; | |
2884 } | |
2885 } | |
2886 | |
2887 return found; | |
2888 } | |
2889 | |
8 | 2890 static void write_flags(Image_tree tree,int *list,int len,int pos) |
0 | 2891 { |
2892 int shift,next_pos,max; | |
2893 | |
2894 if (is_in_list(list,len,pos)) | |
2895 { | |
2896 tree->flag=1; | |
2897 | |
2898 shift=pos-(pow(4,tree->level-1)-1)*4/3-1; | |
2899 max=(int) ((pow(4,tree->level)-1)*4/3); | |
2900 next_pos=max+4*shift+1; | |
2901 | |
2902 write_flags(tree->coarse,list,len,next_pos); | |
2903 write_flags(tree->horizontal,list,len,next_pos+1); | |
2904 write_flags(tree->vertical,list,len,next_pos+2); | |
2905 write_flags(tree->diagonal,list,len,next_pos+3); | |
2906 } | |
2907 } | |
2908 | |
2909 /************************************************************************/ | |
2910 /* Functionname: err_simple_message */ | |
2911 /* -------------------------------------------------------------------- */ | |
2912 /* Parameter: */ | |
2913 /* char *: string that contains information about an */ | |
2914 /* error the user should know. */ | |
2915 /* -------------------------------------------------------------------- */ | |
2916 /* Description: */ | |
2917 /* Prints error messages for the user. */ | |
2918 /************************************************************************/ | |
2919 | |
2920 void err_SimpleMessage(char *message) | |
2921 { | |
2922 fprintf(stderr,"%s\n",message); | |
2923 } | |
2924 | |
2925 /************************************************************************/ | |
2926 /* Functionname: err_get_message */ | |
2927 /* -------------------------------------------------------------------- */ | |
2928 /* Return value: Errormessage for this specific error. */ | |
2929 /* Parameter: */ | |
2930 /* Error err: Error whose errormessage should be returned */ | |
2931 /* -------------------------------------------------------------------- */ | |
2932 /* Description: */ | |
2933 /************************************************************************/ | |
2934 char * err_GetErrorMessage(Error err) | |
2935 { | |
2936 | |
2937 switch (err) | |
2938 { | |
2939 case Error_NotImplemented: | |
2940 return "Sorry, this is not implemented yet. "; | |
2941 break; | |
2942 | |
2943 case Error_AssertionFailed: | |
2944 return "Sorry, an internal assertion was violated.\n" | |
2945 "This action can not be completed. :-("; | |
2946 break; | |
2947 | |
2948 case Error_NotEnoughMemory: | |
2949 return "Sorry, there is not enough memory"; | |
2950 break; | |
2951 | |
2952 case Error_Limitation: | |
2953 return "Some limitation of the program exceeded"; | |
2954 break; | |
2955 | |
2956 /* - FILES - */ | |
2957 | |
2958 case Error_CantOpenFile: | |
2959 return "Could not open file"; | |
2960 break; | |
2961 | |
2962 case Error_CantCreateFile: | |
2963 return "Could not create file"; | |
2964 break; | |
2965 | |
2966 case Error_CantCloseFile: | |
2967 return "Could not close file"; | |
2968 break; | |
2969 | |
2970 case Error_InternalError: | |
2971 return "Sorry, an internal error occured.\n" | |
2972 "This action can not be completed. :-("; | |
2973 break; | |
2974 | |
2975 default: | |
2976 return "Sorry, but an unknown error ocurred.\n" | |
2977 "This action can not be completed. :-("; | |
2978 break; | |
2979 | |
2980 | |
2981 } | |
2982 } |