|
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>
|
|
|
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 }
|
