Mercurial > hg > wm

comparison Meerwald/wavelet.c @ 0:be303a3f5ea8

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

Repositories maintained by Peter Meerwald, pmeerw@pmeerw.net.