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