wm: Meerwald-dir/wavelet.c comparison

comparison Meerwald-dir/wavelet.c @ 24:9f20bce6184e v0.7

move directories, support netpbm 11

author	Peter Meerwald-Stadler <pmeerw@pmeerw.net>
date	Fri, 20 Dec 2024 13:08:59 +0100
parents	Meerwald/wavelet.c@bd669312f068
children

comparison

equal deleted inserted replaced

-:71dd4b96221b
+:9f20bce6184e
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "wavelet.h"
+#include <ctype.h>
+#include <values.h>
+static int read_char(FILE *fp);
+static int read_int(FILE *fp);
+IntImage new_intimage(int width, int height)
+{
+	IntImage image;
+	image = (IntImage) calloc(1,sizeof(struct IntImage_struct));
+	if (image==NULL) goto error;
+	image->data = (IntPixel*) calloc(width*height,sizeof(IntPixel));
+	if (image->data==NULL) goto error;
+	image->width = width;
+	image->height = height;
+	image->size = width*height;
+	image->bpp = 0;
+	image->min_val = (IntPixel) 0;
+	image->max_val = (IntPixel) 0;
+	return image;
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	return NULL;
+}
+Image new_image(int width, int height)
+{
+	Image image;
+	image = (Image) calloc(1,sizeof(struct Image_struct));
+	if (image==NULL) goto error;
+	image->data = (Pixel*) calloc(width*height,sizeof(Pixel));
+	if (image->data==NULL) goto error;
+	image->width = width;
+	image->height = height;
+	image->size = width*height;
+	image->bpp = 0;
+	image->min_val = (Pixel) 0;
+	image->max_val = (Pixel) 0;
+	return image;
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	return NULL;
+}
+void free_intimage(IntImage img)
+{
+	if (img) {
+if (img->data) free(img->data);
+free(img);
+	}
+}
+void free_image(Image img)
+{
+	if (img) {
+if (img->data) free(img->data);
+free(img);
+	}
+}
+/************************************************************************
+*	Functionname:		load_intimage
+* --------------------------------------------------------------------
+*	PARAMETER:
+*	file:		filename of image
+*	max_val:	scaling of grey values to [0..max_val]
+*
+*	RETURN:
+*	Image that shoud be loaded, if not possible return NULL
+* --------------------------------------------------------------------
+* 	DESCRIPTION:
+*	This function loads an IntImage (PGM, ASCII or binary
+*	encoded (P5 or P3 format) ) from a file.
+************************************************************************/
+IntImage load_intimage(char *file, int max_val)
+{
+	IntImage	img;
+	FILE		*fp;
+	IntPixel	*data;
+	int		width, height, i, ch1, ch2;
+	fp=fopen(file,"rb");
+	if (fp==NULL) goto error;
+	ch1=getc(fp);
+	ch2=getc(fp);
+	if (ch1!='P' || (ch2!='5' && ch2!='2')) goto error1;
+	width=read_int(fp);
+	height=read_int(fp);
+	if ((width==0) || (height==0) ) goto error1;
+	read_int(fp);
+	img=new_intimage(width,height);
+	img->bpp=8;
+	data=img->data;
+	for (i=0; i<img->size; i++)
+	{ if (ch2=='5')
+	    *data=getc(fp);
+	  else
+	    *data=read_int(fp);
+	  data++;
+	}
+	fclose(fp);
+	return img;
+error1:
+err_SimpleMessage(err_GetErrorMessage(Error_WrongFileFormat));
+return NULL;
+error:
+err_SimpleMessage(err_GetErrorMessage(Error_CantOpenFile));
+return NULL;
+}
+/************************************************************************
+*	Functionname:		load_image
+* --------------------------------------------------------------------
+*	PARAMETER:
+*	file:		filename of image
+*	max_val:	scaling of grey values to [0..max_val]
+*
+*	RETURN:
+*	Image that shoud be loaded, if not possible return NULL
+* --------------------------------------------------------------------
+* 	DESCRIPTION:
+*	This function loads an IntImage with load_intimage
+*	and then converts to Image.
+************************************************************************/
+Image load_image(char *file, int max_val)
+{
+	Image img;
+	IntImage intimg;
+	intimg = load_intimage(file, max_val);
+	if (!intimg) return NULL;
+	img = intimage_to_image(intimg);
+	if (!intimg) return NULL;
+	return img;
+}
+/************************************************************************/
+/*	Functionname:		save_image_P5                           */
+/* -------------------------------------------------------------------- */
+/*	Parameter:							*/
+/*            img: Image that shoud be saved                            */
+/*            file: filename of image                                   */
+/* -------------------------------------------------------------------- */
+/* 	Description: save an image as PGM (P5 binary decoded) file      */
+/*	                                              			*/
+/************************************************************************/
+int save_image_P5(char *file, Image img)
+{       FILE *fp;
+Pixel *data;
+long i;
+	int p;
+fp=fopen(file,"wb");
+if (fp==NULL)
+goto error;
+fprintf(fp,"P5\n");
+fprintf(fp,"%d %d\n%d ",img->width,img->height,255);
+data=img->data;
+for (i=0;i<img->size;i++) {
+	  p=floor(*data+0.5);
+	  if (p<0) p=0;
+	  if (p>255) p=255;
+/*          putc(*data,fp);	*/
+	  putc(p,fp);
+data++;
+}
+fclose(fp);
+return 1;
+error:
+err_SimpleMessage(err_GetErrorMessage(Error_CantOpenFile));
+return 0;
+}
+void clear_image(Image img)
+{
+	int i;
+	PreCondition(img!=NULL,"image==NULL");
+	for (i=0;i<img->size;i++)
+		(img->data)[i]=(Pixel) 0;
+}
+void copy_into_image(Image img1,Image img2,int x,int y)
+/* copy img2 into img1 at position (x,y)*/
+{
+int start,i,j,aim;
+Pixel *temp;
+temp=img2->data;
+start=img1->width*y+x;
+for (i=0;i<img2->height;i++) {
+for (j=0;j<img2->width;j++) {
+	    aim=start+j+img1->width*i;
+img1->data[aim]=*temp;
+temp++;
+}
+}
+}
+void copy_into_intimage(IntImage img1,IntImage img2,int x,int y)
+{/* copy img2 into img1 at position (x,y)*/
+int start,i,j,aim;
+IntPixel *temp;
+temp=img2->data;
+start=img1->width*y+x;
+for (i=0;i<img2->height;i++)
+{
+for (j=0;j<img2->width;j++)
+{
+	    aim=start+j+img1->width*i;
+img1->data[aim]=*temp;
+temp++;
+}
+}
+}
+void copy_part_of_image_into_image(Image dest_img, int dest_x, int dest_y,
+				   Image src_img, int src_x, int src_y,
+				   int width, int height)
+{
+	Pixel *sp,*dp;
+	int y,siz;
+	sp=get_pixel_adr(src_img,src_x,src_y);
+	dp=get_pixel_adr(dest_img,dest_x,dest_y);
+	siz=width*sizeof(Pixel);
+	for (y=0;y<height;y++)
+	{
+		memcpy(dp,sp,siz);
+		sp+=src_img->width;
+		dp+=dest_img->width;
+	}
+}
+void copy_part_of_image(Image img1,Image img2,int x,int y)
+/* copy part of img2 begining at position (x,y) into img1 */
+{	int i,j,width,height,start,step;
+	Pixel *data;
+	width=img1->width;
+	height=img1->height;
+	start=img2->width*y+x;
+	data=img1->data;
+	for (i=0;i<height;i++) {
+	  step=i*img2->width;
+	  for (j=0;j<width;j++){
+	    *data=img2->data[start+j+step];
+	    data++;
+	  }
+	}
+}
+void scale_image(Image img, int maximum)
+/* scale image to [0..maximum]*/
+{ 	int i;
+	Pixel max = MINDOUBLE, min = MAXDOUBLE, multi;
+	for (i=0;i<img->size;i++) {
+	  if (img->data[i]<min) min=img->data[i];
+	  else if (img->data[i]>max) max=img->data[i];
+	}
+	multi=(Pixel)maximum/(max-min);
+	for (i=0;i<img->size;i++) img->data[i]=multi*(img->data[i]-min);
+	img->min_val=0;
+	img->max_val=maximum;
+}
+int string_to_pixel(char *str, Pixel *p)
+{
+	float ppp;
+	if (sscanf(str,"%f",&ppp)==1)
+	{
+		*p=(Pixel) ppp;
+		return 1;
+	}
+	else
+	{
+		*p=0.0;
+		return 0;
+	}
+}
+Image intimage_to_image(IntImage i)
+{	Image img;
+	int j;
+	img=new_image(i->width,i->height);
+	if (img==NULL) goto error;
+	for (j=0;j<i->size;j++)
+	  img->data[j]=(Pixel)i->data[j];
+	img->bpp=i->bpp;
+return img;
+	error:
+err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  return NULL;
+}
+IntImage image_to_intimage(Image i)
+{       IntImage img;
+int j,multi=1,max,min;
+img=(IntImage)calloc(1,sizeof(struct IntImage_struct));
+if (img==NULL) goto error;
+img->data=(IntPixel *)calloc(i->size,sizeof(IntPixel));
+if (img->data==NULL) goto error;
+img->width=i->width;
+img->height=i->height;
+img->size=i->size;
+img->bpp=i->bpp;
+max=i->max_val;
+min=i->min_val;
+if ((max-min)!=0)
+	  multi=255.0/(max-min);
+for (j=0;j<img->size;j++)
+img->data[j]=(int)((i->data[j]-min)*multi+0.5);
+return img;
+error:
+err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  return NULL;
+}
+static int read_char(FILE *fp)
+/*read a character from file, but skip any comments*/
+{       int ch;
+	ch=getc(fp);
+	if (ch=='#'){
+	   do {
+	     ch=getc(fp);
+	   } while (ch!='\n' && ch!=EOF);
+	}
+	return ch;
+}
+static int read_int(FILE *fp)
+/*read an ascii integer from file, and skip leading tabstops,new lines ...*/
+{  	int r,ch;
+	do {
+	  ch=read_char(fp);
+	} while (ch==' ' || ch=='\n' || ch=='\t');
+	if (ch<'0' || ch>'9')
+	 goto error;
+	r=ch-'0';
+	while ( (ch=read_char(fp)) >='0' && (ch <= '9') ) {
+	  r*=10;
+	  r+=ch-'0';
+	}
+	return r;
+	error:
+	  return 0;
+}
+Image_tree new_image_tree()
+{
+	Image_tree t;
+	t=(Image_tree) calloc(1,sizeof(struct Image_tree_struct));
+	t->entropy=0.0;
+	t->image=NULL;
+	t->coarse=t->horizontal=t->vertical=t->diagonal=t->doubletree=NULL;
+	t->level=0;
+	t->codec_data=NULL;
+	t->significance_map=NULL;
+	return t;
+}
+void free_image_tree(Image_tree t)
+{
+	if (t->coarse) free_image_tree(t->coarse);
+	if (t->horizontal) free_image_tree(t->horizontal);
+	if (t->vertical) free_image_tree(t->vertical);
+	if (t->diagonal) free_image_tree(t->diagonal);
+	if (t->doubletree) free_image_tree(t->doubletree);
+	if (t->image) free_image(t->image);
+	if (t->significance_map) free_intimage(t->significance_map);
+	if (t->codec_data) free(t->codec_data);
+	t->image=NULL;
+	free(t);
+}
+/***********************************************************************
+*      Functionname: get_image_infos
+* --------------------------------------------------------------------
+*      Parameter:
+*          Image image: input image
+*          Pixel *min,*max,*avg,*var: return minimum, maximum,
+*                average and variance of current image
+* --------------------------------------------------------------------
+*      Description:
+*          get statistical information of Image
+************************************************************************/
+void get_image_infos(Image image, Image_info info)
+{
+	int x,y;
+	Pixel p,sp,sp2;
+	sp=sp2=(Pixel)0.0;
+	p=get_pixel(image,0,0);
+	info->min=info->max=p;
+	for (y=0;y<image->height;y++)
+		for (x=0;x<image->width;x++)
+		{
+			p=get_pixel(image,x,y);
+			info->max=MAX(info->max,p);
+			info->min=MIN(info->min,p);
+			sp+=p;
+			sp2+=p*p;
+		}
+	sp=sp/image->width/image->height;
+	sp2=sp2/image->width/image->height;
+	info->mean=sp;
+	info->var=sp2-sp*sp;
+	info->rms=sqrt(sp2);
+}
+/***********************************************************************
+*      Functionname: get_difference_image
+* --------------------------------------------------------------------
+*      Parameter:
+*          Image image1, image 2: input images
+*
+*      Return:
+*          Image : difference of image1 and image 2,
+*                  NULL if error occurs
+************************************************************************/
+Image get_difference_image(Image image1, Image image2)
+{
+	Image diff;
+	int i,max,w,h;
+	Pixel *d,*i1,*i2;
+	if ((!image1) || (!image2)) return NULL;
+	w=image1->width;
+	h=image1->height;
+	if (image2->width != w || image2->height != h) return NULL;
+	diff=new_image(w,h);
+	max=w*h;
+	d=diff->data;
+	i1=image1->data;
+	i2=image2->data;
+	for (i=0;i<max;i++)
+		d[i]=i2[i]-i1[i];
+	return diff;
+}
+/************************************************************************/
+/*      Functionname: get_intimage_infos                                */
+/* -------------------------------------------------------------------- */
+/*      Parameter:                                                      */
+/*          IntImage image: input image                                 */
+/*          IntPixel *min,*max, return minimum, maximum			*/
+/*	       Pixel *avg,*var: average and variance of current image   */
+/*                average and variance of current image                 */
+/* -------------------------------------------------------------------- */
+/*      Description:                                                    */
+/*          get statistical information of Image                        */
+/************************************************************************/
+void get_intimage_infos(IntImage image, IntPixel *min, IntPixel *max, Pixel *avg, Pixel *var)
+{
+	int x,y;
+	Pixel p,sp,sp2;
+	sp=sp2=(Pixel)0.0;
+	p= (Pixel) get_intpixel(image,0,0);
+	*min=*max=p;
+	for (y=0;y<image->height;y++)
+		for (x=0;x<image->width;x++)
+		{
+			p= (Pixel) get_intpixel(image,x,y);
+			*max=MAX(*max, (IntPixel) p);
+			*min=MIN(*min, (IntPixel) p);
+			sp+=p;
+			sp2+=p*p;
+		}
+	sp=sp/image->width/image->height;
+	sp2=sp2/image->width/image->height;
+	*avg=sp;
+	*var=sp2-sp*sp;
+}
+/************************************************************************/
+/*      Functionname: init_zigzag                                       */
+/* -------------------------------------------------------------------- */
+/*      Parameter:                                                      */
+/*          Zigzag_data_struct:                                         */
+/*              output: will be initialized, x/y hold coordinates of    */
+/*                      the first pixel                                 */
+/*          int width,height:                                           */
+/*              input: width/height of image:                           */
+/* -------------------------------------------------------------------- */
+/*      Description:                                                    */
+/*          initializes Zigzag_data structure for use with next_zigzag  */
+/************************************************************************/
+void init_zigzag(Zigzag_data zz, int width, int height)
+{
+	zz->x=0;
+	zz->y=0;
+	zz->dir=zigzag_up;
+	zz->w=width;
+	zz->h=height;
+}
+/************************************************************************/
+/*      Functionname: next_zigzag                                       */
+/* -------------------------------------------------------------------- */
+/*      Parameter:                                                      */
+/*          Zigzag_data_struct:                                         */
+/*              int x,y:                                                */
+/*                  input: current position of zigzag-scan              */
+/*                  output: next position of zigzag-scan                */
+/*              int w,h: width and height of image                      */
+/*              enum zigzag_direction *dir: i/o:                        */
+/*                  direction moving thru the image                     */
+/* -------------------------------------------------------------------- */
+/*      Description:                                                    */
+/*          calculates the next point (x',y') of the zigzag-scan        */
+/*          through the image with size (w,h)                           */
+/************************************************************************/
+void next_zigzag(Zigzag_data zz)
+{
+	switch(zz->dir)
+	{
+	case zigzag_up:
+		if (zz->y==0)
+		{
+			if (zz->x==zz->w-1)
+			{
+				(zz->y)++; zz->dir=zigzag_down;
+			}
+			else
+			{
+				(zz->x)++; zz->dir=zigzag_down;
+			}
+		}
+		else
+		{
+			if (zz->x==zz->w-1)
+			{
+				(zz->y)++; zz->dir=zigzag_down;
+			}
+			else
+			{
+				(zz->x)++; (zz->y)--;
+			}
+		}
+		break;
+	case zigzag_down:
+		if (zz->x==0)
+		{
+			if (zz->y==zz->h-1)
+			{
+				(zz->x)++; zz->dir=zigzag_up;
+			}
+			else
+			{
+				(zz->y)++; zz->dir=zigzag_up;
+			}
+		}
+		else
+		{
+			if (zz->y==zz->h-1)
+			{
+				(zz->x)++; zz->dir=zigzag_up;
+			}
+			else
+			{
+				(zz->x)--;(zz->y)++;
+			}
+		}
+		break;
+	}
+}
+Image get_absolute_image_scaled(Image img)
+{
+	Image out;
+	int x,y;
+	struct Image_info_struct info;
+	Pixel scale,p;
+	out=new_image(img->width,img->height);
+	get_image_infos(img, &info);
+	scale=255/MAX(fabs(info.min),fabs(info.max));
+	for(y=0;y<img->height;y++)
+	for(x=0;x<img->width;x++)
+	{
+		p=get_pixel(img,x,y)*scale;
+		set_pixel(out,x,y,p);
+	}
+	return out;
+}
+#define FLOOR_HALF(x) ((x)&1 ? ((x)-1)/2 : (x)/2)
+#define CEILING_HALF(x) ((x)&1 ? ((x)+1)/2 : (x)/2)
+#define MOD(a,b) ( (a)<0 ? ((b)-((-(a))%(b))) : (a)%(b) )
+Filter new_filter(int size)
+{
+	Filter f;
+	Entering;
+	f=(Filter) calloc(1,sizeof(struct FilterStruct));
+	f->data=(Pixel *)calloc(size,sizeof(Pixel));
+	f->len=size;
+	f->hipass=0;
+	Leaving;
+	return f;
+}
+Pixel get_filter_center(Filter f)
+{
+	int i;
+	Pixel p, sum, norm;
+	if (f==NULL) return 0;
+	sum=norm=0;
+	for (i=0;i<f->len;i++)
+	{
+		p=f->data[i];
+		p=p*p;
+		norm += p;
+		sum += (i+f->start)*p;
+	}
+	p=sum/norm;
+	return p;
+}
+int filter_cutoff(Image in, int in_start, int in_len, int in_step,
+		Image out, int out_start, int out_len, int out_step,
+		Filter f)
+{
+	int i,i2,j;
+	Pixel *out_pix, *in_pix, *f_data;
+	int fstart,fend; /* Source interval */
+	Entering;
+	PreCondition(out_len == in_len/2,"out_len != in_len/2 !!!");
+/* convolution: out[i]=sum_{j=start}^{end} (in[2*i-j]*f[j])
+boundaries:	image in [in_start ... in_start + in_len-1]
+		image out [out_start ... out_start + out_len-1]
+		filter f [0..f->len-1] = [f->start .. f->end]
+		cutoff at:
+*/
+	for (i=0;i<out_len;i++)
+	{
+		i2=2*i;
+		fstart=i2-(in_len-1);
+		fstart=MAX(fstart,f->start);
+		fend=MIN(i2,f->end);
+#ifdef TRACE
+		sprintf(dbgstr,"i=%d fstart=%d fend=%d\n",i,fstart,fend);
+		Trace(dbgstr);
+#endif
+		out_pix=out->data+out_start+i*out_step;
+		in_pix=in->data+in_start+(i2-fstart)*in_step;
+		f_data=f->data-f->start+fstart;
+		for (j=fstart;j<=fend;j++,in_pix-=in_step,f_data++)
+		{
+			*out_pix += (*f_data) * (*in_pix);
+#ifdef TRACE
+			sprintf(dbgstr,"     j=%2d in: %4.2f filter: %4.2f [%4d] [%4d]  opt : %4.2f %4.2f\n",
+				j,
+				in->data[in_start+in_step*(i2-j)],
+				f->data[j-f->start],
+				in_start+in_step*(i2-j),
+				j-f->start,
+				*in_pix, *f_data);
+			Trace(dbgstr);
+#endif
+		}
+	}
+	Leaving;
+	return 1;
+}
+int filter_inv_cutoff(Image in, int in_start, int in_len, int in_step,
+		Image out, int out_start, int out_len, int out_step,
+		Filter f)
+{
+	int i,j;
+	Pixel *out_pix, *in_pix, *f_data;
+	int fstart,fend; /* Source interval */
+	Entering;
+	PreCondition(out_len == in_len*2,"out_len != in_len*2 !!!");
+/* convolution: out[i]=sum_{j=start}^{end} (f[2*j-i]*in[j])
+boundaries:	image in [in_start ... in_start + in_len-1]
+		image out [out_start ... out_start + out_len-1]
+		filter f [0..f->len-1] = [f->start .. f->end]
+		cutoff at:
+*/
+	for (i=0;i<out_len;i++)
+	{
+		fstart=CEILING_HALF(f->start+i);
+		fend=FLOOR_HALF(f->end+i);
+		fstart=MAX(fstart,0);
+		fend=MIN(fend,in_len-1);
+#ifdef TRACE
+		sprintf(dbgstr,"i=%d fstart=%d fend=%d\n",i,fstart,fend);
+		Trace(dbgstr);
+#endif
+		out_pix=out->data+out_start+i*out_step;
+		in_pix=in->data+in_start+fstart*in_step;
+		f_data=f->data-f->start+2*fstart-i;
+		for (j=fstart;j<=fend;j++,in_pix+=in_step,f_data+=2)
+		{
+			*out_pix += (*f_data) * (*in_pix);
+#ifdef TRACE
+			sprintf(dbgstr,"     j=%2d in: %4.2f filter: %4.2f [%4d] [%4d]  opt : %4.2f %4.2f\n",
+				j,
+				in->data[in_start+j*in_step],
+				f->data[2*j-i-f->start],
+				in_start+j*in_step,
+				2*j-i-f->start,
+				*in_pix, *f_data);
+			Trace(dbgstr);
+#endif
+		}
+	}
+	Leaving;
+	return 1;
+}
+int filter_periodical(Image in, int in_start, int in_len, int in_step,
+		Image out, int out_start, int out_len, int out_step,
+		Filter f)
+{
+	int i,i2,j;
+	Pixel *out_pix, *in_pix, *f_data;
+	int fstart,fend;
+	int istart;
+	int ipix_len;
+	Entering;
+	PreCondition(out_len == in_len/2,"out_len != in_len/2 !!!");
+/* convolution: out[i]=sum_{j=start}^{end} (in[2*i-j]*f[j])
+boundaries:	image in [in_start ... in_start + in_len-1]
+		image out [out_start ... out_start + out_len-1]
+		filter f [0..f->len-1] = [f->start .. f->end]
+*/
+	ipix_len=in_len*in_step;
+	for (i=0;i<out_len;i++)
+	{
+		i2=2*i;
+		fstart=f->start;
+		fend=f->end;
+		istart=(i2-fstart);
+		istart=MOD(istart,in_len);
+#ifdef TRACE
+		sprintf(dbgstr,"i=%d istart=%d\n",i,istart);
+		Trace(dbgstr);
+#endif
+		out_pix=out->data+out_start+i*out_step;
+		in_pix=in->data+in_start+istart*in_step;
+		f_data=f->data;
+		for (j=fstart;j<=fend;j++,f_data++)
+		{
+			*out_pix += (*f_data) * (*in_pix);
+#ifdef TRACE
+			sprintf(dbgstr,"     j=%2d in: %4.2f filter: %4.2f [%4d] [%4d]  opt : %4.2f %4.2f\n",
+				j,
+				in->data[in_start+in_step*((i2-j+in_len)%in_len)],
+				f->data[j-f->start],
+				in_start+in_step*((i2-j+in_len)%in_len),
+				j-f->start,
+				*in_pix, *f_data);
+			Trace(dbgstr);
+#endif
+			in_pix-=in_step;
+			istart--;
+			if (istart<0)
+			{
+				istart+=in_len;
+				in_pix+=ipix_len;
+			}
+		}
+	}
+	Leaving;
+	return 1;
+}
+int filter_inv_periodical(Image in, int in_start, int in_len, int in_step,
+		Image out, int out_start, int out_len, int out_step,
+		Filter f)
+{
+	int i,j;
+	Pixel *out_pix, *in_pix, *f_data;
+	int fstart,fend; /* Source interval */
+	int istart;
+	int ipix_len;
+	Entering;
+	PreCondition(out_len == in_len*2,"out_len != in_len*2 !!!");
+/* convolution: out[i]=sum_{j=start}^{end} (f[2*j-i]*in[j])
+boundaries:	image in [in_start ... in_start + in_len-1]
+		image out [out_start ... out_start + out_len-1]
+		filter f [0..f->len-1] = [f->start .. f->end]
+*/
+	ipix_len=in_len*in_step;
+	for (i=0;i<out_len;i++)
+	{
+		fstart=CEILING_HALF(f->start+i);
+		fend=FLOOR_HALF(f->end+i);
+		istart=MOD(fstart,in_len);
+#ifdef TRACE
+		sprintf(dbgstr,"i=%d fstart=%d fend=%d istart=%d\n",i,fstart,fend,istart);
+		Trace(dbgstr);
+#endif
+		out_pix=out->data+out_start+i*out_step;
+		in_pix=in->data+in_start+istart*in_step;
+		f_data=f->data-f->start+2*fstart-i;
+		for (j=fstart;j<=fend;j++,f_data+=2)
+		{
+			*out_pix += (*f_data) * (*in_pix);
+#ifdef TRACE
+			sprintf(dbgstr,"     j=%2d in: %4.2f filter: %4.2f [%4d] [%4d]  opt : %4.2f %4.2f\n",
+				j,
+				in->data[in_start+(j % in_len)*in_step],
+				f->data[2*j-i-f->start],
+				in_start+(j%in_len)*in_step,
+				2*j-i-f->start,
+				*in_pix, *f_data);
+			Trace(dbgstr);
+#endif
+			in_pix+=in_step;
+			istart++;
+			if (istart>=in_len)
+			{
+				istart-=in_len;
+				in_pix-=ipix_len;
+			}
+		}
+	}
+	Leaving;
+	return 1;
+}
+int filter_mirror(Image in, int in_start, int in_len, int in_step,
+		Image out, int out_start, int out_len, int out_step,
+		Filter f)
+{
+	int i,i2,j;
+	Pixel *out_pix, *in_pix, *f_data;
+	int fstart,fend;
+	int in_pos;
+	Entering;
+	PreCondition(out_len == in_len/2,"out_len != in_len/2 !!!");
+/* convolution: out[i]=sum_{j=start}^{end} (in[2*i-j]*f[j])
+boundaries:	image in [in_start ... in_start + in_len-1]
+		image out [out_start ... out_start + out_len-1]
+		filter f [0..f->len-1] = [f->start .. f->end]
+*/
+	in_pix=in->data+in_start;
+	for (i=0;i<out_len;i++)
+	{
+		i2=2*i;
+		fstart=f->start;
+		fend=f->end;
+		out_pix=out->data+out_start+i*out_step;
+		f_data=f->data;
+		for (j=fstart;j<=fend;j++)
+		{
+			in_pos=(i2-j);
+			if (in_pos<0)
+			{
+				in_pos=-in_pos;
+				if (in_pos>=in_len) continue;
+			}
+			if (in_pos>=in_len)
+			{
+				in_pos=2*in_len-2-in_pos;
+				if (in_pos<0) continue;
+			}
+			*out_pix += (f_data[j-fstart]) * (in_pix[in_pos*in_step]);
+		}
+	}
+	Leaving;
+	return 1;
+}
+int filter_inv_mirror(Image in, int in_start, int in_len, int in_step,
+		Image out, int out_start, int out_len, int out_step,
+		Filter f)
+{
+	int i,j;
+	Pixel *out_pix, *in_pix;
+	int fstart,fend; /* Source interval */
+	int in_pos;
+	Entering;
+	PreCondition(out_len == in_len*2,"out_len != in_len*2 !!!");
+/* convolution: out[i]=sum_{j=start}^{end} (f[2*j-i]*in[j])
+boundaries:	image in [in_start ... in_start + in_len-1]
+		image out [out_start ... out_start + out_len-1]
+		filter f [0..f->len-1] = [f->start .. f->end]
+*/
+	/*fprintf(stderr,"inv started\n");*/
+	for (i=0;i<out_len;i++)
+	{
+		fstart=CEILING_HALF(f->start+i);
+		fend=FLOOR_HALF(f->end+i);
+		out_pix=out->data+out_start+i*out_step;
+		in_pix=in->data+in_start;
+/*
+		printf("in: %4d - %4d  flt: %4d - %4d   [%s]\n",fstart,fend,2*fstart-i,2*fend-i,
+		(2*fstart-i<f->start || 2*fend-i>f->end) ? "error":"ok");
+*/
+		/*fprintf(stderr,"inv[%d]\n",i);*/
+		for (j=fstart;j<=fend;j++)
+		{
+			in_pos=j;
+			if (in_pos<0)
+			{
+				if (f->hipass)
+					in_pos=-in_pos-1;
+				else
+					in_pos=-in_pos;
+				if (in_pos>=in_len) continue;
+			}
+			if (in_pos>=in_len)
+			{
+				if (f->hipass)
+					in_pos=2*in_len-2-in_pos;
+				else
+					in_pos=2*in_len-1-in_pos;
+				if (in_pos<0) continue;
+			}
+			/*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]);*/
+			*out_pix += f->data[2*j-i-f->start] * (in_pix[in_pos*in_step]);
+		}
+	}
+	Leaving;
+	return 1;
+}
+#define MAX_LINE 256
+#define skip_blank(str) { while(isspace(*(str))) (str)++; }
+static int get_next_line(FILE *f, char *c)
+{
+char *str,string[200];
+int len;
+do
+{
+str=string;
+if (!fgets(str,200,f))
+		{
+			Trace("get_next_line: eof\n");
+			goto error;
+		}
+len=strlen(str);
+while (len>=1 && isspace(str[len-1])) str[--len]=0;
+while (isspace(*str)) str++;
+}
+while (strlen(str)==0 || *str=='#');
+strcpy(c,str);
+	return 1;
+error:
+	return 0;
+}
+static int next_line_str(FILE *f, char *tag, char *out)
+{
+	char str[MAX_LINE],*t_str;
+	if (!get_next_line(f,str)) goto error;
+	t_str=strtok(str," ");
+	if (!t_str || strcmp(t_str,tag)) goto error;
+	t_str=strtok(NULL,"\n");
+	if (!t_str) goto error;
+	skip_blank(t_str);
+	strcpy(out,t_str);
+	return 1;
+error:
+	return 0;
+}
+static int next_line_str_alloc(FILE *f, char *tag, char **out)
+{
+	char str[MAX_LINE];
+	if (!next_line_str(f,tag,str)) goto error;
+	*out=malloc(strlen(str)+1);
+	strcpy(*out,str);
+	return 1;
+error:
+	return 0;
+}
+static int next_line_int(FILE *f, char *tag, int *out)
+{
+	char str[MAX_LINE];
+	if (next_line_str(f,tag,str) && sscanf(str,"%d",out)==1)
+		return 1;
+	else
+		return 0;
+}
+static Filter read_filter(FILE *f)
+{
+	char str[MAX_LINE];
+	Filter filter;
+	int i;
+	Entering;
+	filter=calloc(1,sizeof(struct FilterStruct));
+	if (!next_line_str(f,"Type",str)) goto error1;
+	if (!strcmp(str,"nosymm"))
+	{
+		filter->type=FTNoSymm;
+	}
+	else if (!strcmp(str,"symm"))
+	{
+		filter->type=FTSymm;
+	}
+	else if (!strcmp(str,"antisymm"))
+	{
+		filter->type=FTAntiSymm;
+	}
+	else
+		goto error1;
+	if (!next_line_int(f,"Length",&(filter->len))) goto error1;
+	if (!next_line_int(f,"Start",&(filter->start))) goto error1;
+	if (!next_line_int(f,"End",&(filter->end))) goto error1;
+	if ((filter->end-filter->start+1!=filter->len))
+	{
+		Trace("error: len != end-start+1\n");
+		goto error1;
+	}
+	filter->data=calloc(filter->len,sizeof(Pixel));
+	for (i=0;i<filter->len;i++)
+	{
+		if (!get_next_line(f,str)) goto error2;
+		if (!string_to_pixel(str,filter->data+i))
+		{
+			Trace("error: invalid filter-value\n");
+			goto error2;
+		}
+	}
+	if (!get_next_line(f,str)) goto error2;
+	if (*str!='}')
+	{
+		Trace("error: '}' not found\n");
+		goto error2;
+	}
+	Leaving;
+	return filter;
+error2:
+	free(filter->data);
+error1:
+	free(filter);
+	LeavingErr;
+	return NULL;
+}
+static FilterGH read_filter_gh(FILE *f)
+{
+	char str[MAX_LINE];
+	FilterGH fgh;
+	Filter filter;
+	int i,max;
+	Entering;
+	fgh=calloc(1,sizeof(struct FilterGHStruct));
+	if (!next_line_str_alloc(f,"Name",&(fgh->name)))
+	{
+		Trace("error: 'Name' tag not found\n");
+		goto error1;
+	}
+	if (!next_line_str(f,"Type",str))
+	{
+		Trace("error: 'Type' tag not found\n");
+		goto error1;
+	}
+	if (!strcmp(str,"orthogonal"))
+	{
+		fgh->type=FTOrtho;
+		max=2;
+	}
+	else if (!strcmp(str,"biorthogonal"))
+	{
+		fgh->type=FTBiOrtho;
+		max=4;
+	}
+	else if (!strcmp(str,"other"))
+	{
+		fgh->type=FTOther;
+		max=4;
+	}
+	else
+	{
+		Trace("error: expecting 'orthogonal', 'biorthogonal' or 'other' type-tag\n");
+		goto error1;
+	}
+	for (i=0;i<max;i++)
+	{
+		if (!get_next_line(f,str)) goto error2;
+		if (*str!='{')
+		{
+			Trace("error: '{' not found\n");
+			goto error2;
+		}
+		if (!(filter=read_filter(f)))
+		{
+			Trace("error: read_filter failed\n");
+			goto error2;
+		}
+		filter->hipass = !(i&1);
+		switch(i)
+		{
+		case 0: fgh->g=filter; break;
+		case 1: fgh->h=filter; break;
+		case 2: fgh->gi=filter; break;
+		case 3: fgh->hi=filter; break;
+		}
+	}
+	if (!get_next_line(f,str)) goto error2;
+	if (*str!='}')
+	{
+		Trace("error: '}' not found\n");
+		goto error2;
+	}
+	Leaving;
+	return fgh;
+error2:
+	if (fgh->g) free(fgh->g);
+	if (fgh->h) free(fgh->h);
+	if (fgh->gi) free(fgh->gi);
+	if (fgh->hi) free(fgh->hi);
+error1:
+	free(fgh);
+	LeavingErr;
+	return NULL;
+}
+AllFilters load_filters(char *name)
+{
+	FILE *f;
+	char str[MAX_LINE];
+	AllFilters a;
+	FilterGH fgh;
+	Entering;
+	PreCondition(name!=NULL,"name=NULL!");
+	f=fopen(name,"rt");
+	if (!f)
+	{
+		Trace("error: fopen failed\n");
+		goto error1;
+	}
+	a=calloc(1,sizeof(struct AllFilterStruct));
+	a->count=0;
+	while (get_next_line(f,str))
+	{
+		if (*str=='{')
+		{
+			fgh=read_filter_gh(f);
+			if (!fgh)
+			{
+				Trace("error: read_filter returned NULL\n");
+				goto error2;
+			}
+			if (a->count)
+				a->filter=realloc(a->filter,(a->count+1)*sizeof(FilterGH));
+			else
+				a->filter=malloc(sizeof(FilterGH));
+			(a->filter)[a->count]=fgh;
+			a->count++;
+		}
+		else
+		{
+			Trace("error: '{' not found\n");
+			goto error2;
+		}
+	}
+	fclose(f);
+	Leaving;
+	return a;
+error2:
+	fclose(f);
+error1:
+	LeavingErr;
+	return 0;
+}
+#define doubletree_min 32
+#define best_basis_min 8
+static int convolute_lines(Image output,Image input,Filter flt,enum FilterMethod method);
+static int convolute_rows(Image output,Image input,Filter flt,enum FilterMethod method);
+static int decomposition(Image t_img,Image coarse,Image horizontal,Image vertical,
+Image diagonal,Filter g,Filter h,enum FilterMethod method);
+static int compute_best(Image_tree tree,int level,int max_level,FilterGH *flt,
+enum FilterMethod method,enum Information_Cost cost,double epsilon);
+static double compute_entropy(Image img,enum Information_Cost cost,double epsilon);
+static void compute_levels(Image_tree tree,double *entropies,enum Information_Cost cost,double epsilon);
+static void free_levels(Image_tree tree,int best);
+static Pixel sumationq(Image img);
+static Pixel normq(Image_tree tree);
+static Pixel sumation_down(Image_tree tree, Pixel normq);
+static Pixel compute_non_additive(Image_tree tree,int size,enum Information_Cost cost,double
+epsilon,int down);
+/************************************************************************/
+/*	Functionname: wavelettransform					*/
+/* -------------------------------------------------------------------- */
+/*	Parameter: 							*/
+/*		original: Image that should be transformed  		*/
+/*	        level: transform down to level                 		*/
+/*		flt: transform with filters flt[0..level]		*/
+/*		method: method of filtering				*/
+/* -------------------------------------------------------------------- */
+/* 	Description: Carries out the wavelettransform 			*/
+/*	                           		           		*/
+/************************************************************************/
+Image_tree wavelettransform(Image original,int level,FilterGH *flt,enum FilterMethod method)
+{	int i,width,height,min,max_level,e;
+Image coarsei,horizontali,verticali,diagonali,tempi;
+	Image_tree ret_tree,temp_tree;
+	width=original->width;
+	height=original->height;
+	tempi=new_image(width,height);
+	if(!tempi) goto error;
+copy_into_image(tempi,original,0,0);
+	ret_tree=new_image_tree();
+	if(!ret_tree) goto error;
+	temp_tree=ret_tree;
+	ret_tree->level=0;
+	min=original->width;
+	if (original->height<min) min=original->height;
+	max_level=log(min)/log(2)-2;
+	if (max_level<level) level=max_level;
+	if (level<1)  /* do not transform */
+	{
+		ret_tree->image=tempi;
+		return ret_tree;
+	}
+	/* decomposition */
+	for (i=0;i<level;i++)
+	{
+	        width=(width+1)/2;
+		height=(height+1)/2;
+		coarsei=new_image(width,height);
+		horizontali=new_image(width,height);
+		verticali=new_image(width,height);
+		diagonali=new_image(width,height);
+		if(!coarsei||!horizontali||!verticali||!diagonali) goto error;
+		e=decomposition(tempi,coarsei,horizontali,verticali,diagonali,
+	                       flt[i]->g,flt[i]->h,method);
+		if (!e) return NULL;
+		temp_tree->coarse=new_image_tree();
+		temp_tree->horizontal=new_image_tree();
+		temp_tree->vertical=new_image_tree();
+		temp_tree->diagonal=new_image_tree();
+		temp_tree->coarse->level=i+1;
+		temp_tree->horizontal->level=i+1;
+		temp_tree->vertical->level=i+1;
+		temp_tree->diagonal->level=i+1;
+		temp_tree->horizontal->image=horizontali;
+		temp_tree->vertical->image=verticali;
+		temp_tree->diagonal->image=diagonali;
+	  	free_image(tempi);
+		if (i!=(level-1))
+		{
+	    		tempi=new_image(width,height);
+	    		copy_into_image(tempi,coarsei,0,0);
+	    		free_image(coarsei);
+	    		/*if i=level coarsei is inserted into the image tree
+	    		  so we should not free coarsei on level-1*/
+	 	}
+		temp_tree=temp_tree->coarse;
+	}
+	temp_tree->image=coarsei;
+	return ret_tree;
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  	return NULL;
+}
+static Image_tree wavelettransform__wp(Image original, int current_level, int level, FilterGH *flt, enum FilterMethod method)
+{
+int i, width, height, min, max_level, e;
+Image coarse_image,horizontal_image,vertical_image,diagonal_image,temp_image;
+Image_tree return_tree, temp_tree;
+width = original->width;
+height = original->height;
+temp_image = new_image(width, height);
+if (!temp_image) goto error;
+copy_into_image(temp_image, original, 0, 0);
+temp_tree = return_tree = new_image_tree();
+if (!return_tree) goto error;
+temp_tree->level = current_level;
+min = original->width;
+if (original->height < min) min = original->height;
+max_level = log(min) / log(2) - 2;
+if (max_level < level) level = max_level;
+if (current_level >= level) {
+return_tree->image = temp_image;
+return return_tree;
+}
+for (i = current_level; i < level; i++) {
+width = (width + 1) / 2;
+height = (height + 1) / 2;
+coarse_image = new_image(width, height);
+horizontal_image = new_image(width, height);
+vertical_image = new_image(width, height);
+diagonal_image = new_image(width, height);
+if (!coarse_image || !horizontal_image ||
+!vertical_image || !diagonal_image) goto error;
+e = decomposition(temp_image, coarse_image, horizontal_image,
+vertical_image, diagonal_image,
+	                          flt[i]->g, flt[i]->h, method);
+if (!e) return NULL;
+temp_tree->coarse = new_image_tree();
+temp_tree->coarse->level = i+1;
+temp_tree->horizontal = wavelettransform__wp(horizontal_image, i+1, level, flt, method);
+temp_tree->vertical = wavelettransform__wp(vertical_image, i+1, level, flt, method);
+temp_tree->diagonal = wavelettransform__wp(diagonal_image, i+1, level, flt, method);
+free_image(horizontal_image);
+free_image(vertical_image);
+free_image(diagonal_image);
+free_image(temp_image);
+if (i != (level - 1)) {
+temp_image = new_image(width, height);
+copy_into_image(temp_image, coarse_image, 0, 0);
+free_image(coarse_image);
+}
+temp_tree = temp_tree->coarse;
+}
+temp_tree->image = coarse_image;
+return return_tree;
+error:
+err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+return NULL;
+}
+Image_tree wavelettransform_wp(Image original, int level, FilterGH *flt, enum FilterMethod method) {
+return wavelettransform__wp(original, 0, level, flt, method);
+}
+/************************************************************************/
+/*	Functionname: best_basis					*/
+/* -------------------------------------------------------------------- */
+/*	Parameter:							*/
+/*	    	original: Image to be transformed			*/
+/*		level: best basis selection down to this level		*/
+/*		flt: transform with filters flt[0..level]		*/
+/*		method: transform with filter method			*/
+/*		cost: carry best basis selection out with this costfunc */
+/*		epsilon: limit for threshold method			*/
+/* -------------------------------------------------------------------- */
+/* 	Description: carries best basis and near best basis selection	*/
+/*			out						*/
+/************************************************************************/
+Image_tree best_basis(Image original,int level,FilterGH *flt,
+				enum FilterMethod method,enum Information_Cost cost,double epsilon)
+{       Image_tree tree;
+	Image img;
+	int min,max_level;
+	tree=new_image_tree();
+	if(!tree) goto error;
+	img=new_image(original->width,original->height);
+	if(!img) goto error;
+copy_into_image(img,original,0,0);
+tree->image=img;
+	min=original->width;
+	if (original->height<min) min=original->height;
+	max_level=log10((float) min/best_basis_min)/log10(2);
+	if (max_level>level) max_level=level;
+	compute_best(tree,0,max_level,flt,method,cost,epsilon);
+	if (!tree->image) free(img);
+return tree;
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  	return NULL;
+}
+/************************************************************************/
+/*	Functionname: best_level_selection				*/
+/* -------------------------------------------------------------------- */
+/*	Parameter: 							*/
+/*		original: Image that should be transformed  		*/
+/*	        maxlevel: transform down to level                 	*/
+/*		flt: transform with filters flt[0..level]		*/
+/*		method: transform with filter method			*/
+/*		cost: carry best basis selection out with this costfunc */
+/*		epsilon: limit for threshold method			*/
+/* -------------------------------------------------------------------- */
+/* 	Description: Carries out the best level selection		*/
+/*	                           		           		*/
+/************************************************************************/
+Image_tree best_level(Image original,int maxlevel,int *bestlevel,FilterGH *flt,enum FilterMethod method,
+				enum Information_Cost cost,double epsilon)
+{	Image_tree tree;
+	Image img;
+	double *entropies,min;
+	int best=0,i,e;
+	img=new_image(original->width,original->height);
+	copy_into_image(img,original,0,0);
+	tree=new_image_tree();
+	tree->image=img;
+	entropies=(double *)calloc(maxlevel+1,sizeof(double));
+	if(!entropies) goto error;
+	/* decompose down to maxlevel */
+	e=decompose_all(tree,maxlevel,flt,method,cost,epsilon);
+	if (!e) return NULL;
+	/* compute costs of each level and store it in entropies array*/
+	compute_levels(tree,entropies,cost,epsilon);
+	min=entropies[0];
+	for (i=1;i<=maxlevel;i++)
+	{
+		if (entropies[i]<min)
+		{
+			min=entropies[i];
+			best=i;
+		}
+	}
+	/* free all other levels */
+	free_levels(tree,best);
+	*bestlevel=best;
+	return tree;
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  	return NULL;
+}
+/************************************************************************/
+/*	Functionname: compute_best					*/
+/* -------------------------------------------------------------------- */
+/*	Parameter: 							*/
+/*		img: Image that should be transformed	 		*/
+/*	        level: transform level	                		*/
+/*		max_level: transform to maximum level			*/
+/*		flt: transform with filters flt[0..level]		*/
+/*		method: transform with filter method			*/
+/*		cost: carry best basis selection out with this costfunc */
+/*		epsilon: limit for threshold method			*/
+/* -------------------------------------------------------------------- */
+/* 	Description: Carries out the best basis selection (recursivly)  */
+/*	             (is needed by the waveletpacket procedure)		*/
+/************************************************************************/
+static int compute_best(Image_tree tree,int level,int max_level,FilterGH *flt,
+enum FilterMethod method,enum Information_Cost cost,double epsilon)
+{ 	Image coarse,horizontal,vertical,diagonal;
+int e,width,height;
+	double sum;
+	tree->level=level;
+	/* non additive cost function*/
+	if (cost>=shanon)
+	{
+		tree->entropy=compute_non_additive(tree,tree->image->size,cost,epsilon,0);
+	}
+	/*additive cost function*/
+	else 	tree->entropy=compute_entropy(tree->image,cost,epsilon);
+	if (level<max_level) {
+	width=(tree->image->width+1)/2;
+	height=(tree->image->height+1)/2;
+	  	tree->coarse=new_image_tree();
+	  	tree->horizontal=new_image_tree();
+		tree->vertical=new_image_tree();
+	  	tree->diagonal=new_image_tree();
+	coarse=new_image(width,height);
+	horizontal=new_image(width,height);
+	vertical=new_image(width,height);
+	diagonal=new_image(width,height);
+		if(!coarse||!vertical||!horizontal||!diagonal) goto error;
+	        e=decomposition(tree->image,coarse,horizontal,vertical,diagonal,flt[0]->g,flt[0]->h,method);
+	        if (!e) return 0;
+	        tree->coarse->image=coarse;
+	        tree->horizontal->image=horizontal;
+	        tree->vertical->image=vertical;
+	        tree->diagonal->image=diagonal;
+	e=compute_best(tree->coarse,level+1,max_level,flt+1,method,cost,epsilon);
+	e=compute_best(tree->horizontal,level+1,max_level,flt+1,method,cost,epsilon);
+	e=compute_best(tree->vertical,level+1,max_level,flt+1,method,cost,epsilon);
+	e=compute_best(tree->diagonal,level+1,max_level,flt+1,method,cost,epsilon);
+	if (!e) return 0;
+		/*going back in recursion*/
+		if (cost>=shanon) {
+			 sum=compute_non_additive(tree,tree->image->size,cost,epsilon,1);
+		}
+		else sum=(tree->coarse->entropy)+(tree->horizontal->entropy)
+			+(tree->vertical->entropy)+(tree->diagonal->entropy);
+	  	if (tree->entropy>sum)
+	  	{
+	  		tree->entropy=sum;
+	      		free_image(tree->image);	/* take down tree */
+	      		tree->image=NULL;
+	  	}
+	  	else
+	  	{   				/* delete the tree downwards */
+		free_image_tree(tree->coarse);
+		free_image_tree(tree->horizontal);
+		free_image_tree(tree->vertical);
+		free_image_tree(tree->diagonal);
+		tree->coarse=tree->vertical=tree->horizontal=tree->diagonal=NULL;
+	  	}
+	}
+	return 1;
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	return 0;
+}
+/************************************************************************/
+/*	Functionname:		decompose_all                         */
+/* -------------------------------------------------------------------- */
+/*	Parameter:							*/
+/*            	tree: Image tree to be decomposed			*/
+/*           	maxlevel: decompose down to level	                */
+/*		flt: transform with filters flt[0..maxlevel]		*/
+/*		method: transform with filter method			*/
+/*		cost: cost function for entropy computing		*/
+/*		epsilon: limit for threshold method			*/
+/* -------------------------------------------------------------------- */
+/* 	Description: whole decompositing down to maxlevel          	*/
+/*		The original image must be in tree->image		*/
+/************************************************************************/
+int decompose_all(Image_tree tree,int maxlevel,FilterGH *flt,enum FilterMethod method,
+				enum Information_Cost cost,double epsilon)
+{
+	Image original,coarse,horizontal,vertical,diagonal;
+	int e,width,height,level;
+	if (tree->level<maxlevel)
+	{
+		tree->coarse=new_image_tree();
+		tree->horizontal=new_image_tree();
+		tree->vertical=new_image_tree();
+		tree->diagonal=new_image_tree();
+		original=tree->image;
+		width=(original->width+1)/2;
+		height=(original->height+1)/2;
+		level=tree->level;
+		coarse=new_image(width,height);
+		horizontal=new_image(width,height);
+		vertical=new_image(width,height);
+		diagonal=new_image(width,height);
+		if(!coarse||!vertical||!horizontal||!diagonal) goto error;
+		e=decomposition(tree->image,coarse,horizontal,vertical,diagonal,flt[0]->g,flt[0]->h,method);
+		if (!e) return 0;
+		tree->coarse->image=coarse;
+		tree->horizontal->image=horizontal;
+		tree->vertical->image=vertical;
+		tree->diagonal->image=diagonal;
+		tree->coarse->entropy=compute_entropy(coarse,cost,epsilon);
+		tree->horizontal->entropy=compute_entropy(horizontal,cost,epsilon);
+		tree->vertical->entropy=compute_entropy(vertical,cost,epsilon);
+		tree->diagonal->entropy=compute_entropy(diagonal,cost,epsilon);
+		tree->coarse->level=tree->horizontal->level=
+			tree->vertical->level=tree->diagonal->level=level+1;
+		e=decompose_all(tree->coarse,maxlevel,flt+1,method,cost,epsilon);
+		e=decompose_all(tree->horizontal,maxlevel,flt+1,method,cost,epsilon);
+		e=decompose_all(tree->vertical,maxlevel,flt+1,method,cost,epsilon);
+		e=decompose_all(tree->diagonal,maxlevel,flt+1,method,cost,epsilon);
+		if (!e) return 0;
+	}
+	return 1;
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+		return 0;
+}
+/************************************************************************/
+/*	Functionname:		compute_levels                          */
+/* -------------------------------------------------------------------- */
+/*	Parameter:							*/
+/*            	tree: Image tree where the entropy should be computed	*/
+/*           	entropies : array for entropy          	                */
+/*		cost: carry best basis selection out with this costfunc */
+/*		epsilon: limit for threshold method			*/
+/* -------------------------------------------------------------------- */
+/* 	Description: Compute the entropies of all decomposition	levels	*/
+/************************************************************************/
+static void compute_levels(Image_tree tree,double *entropies,enum Information_Cost cost,double epsilon)
+{
+	if (tree->image){
+		entropies[tree->level]+=compute_entropy(tree->image,cost,epsilon);
+	}
+	if (tree->coarse) compute_levels(tree->coarse,entropies,cost,epsilon);
+	if (tree->horizontal) compute_levels(tree->horizontal,entropies,cost,epsilon);
+	if (tree->vertical) compute_levels(tree->vertical,entropies,cost,epsilon);
+	if (tree->diagonal) compute_levels(tree->diagonal,entropies,cost,epsilon);
+}
+/************************************************************************/
+/*	Functionname:		free_levels	                        */
+/* -------------------------------------------------------------------- */
+/*	Parameter:							*/
+/*            	tree: Image tree which should be cleaned		*/
+/*           	best: best level	         	                */
+/* -------------------------------------------------------------------- */
+/* 	Description: clean the image tree except the best level      	*/
+/************************************************************************/
+static void free_levels(Image_tree tree,int best)
+{
+	if (tree->level<best)
+	{
+		free_image(tree->image);
+		tree->image=NULL;
+		free_levels(tree->coarse,best);
+		free_levels(tree->horizontal,best);
+		free_levels(tree->vertical,best);
+		free_levels(tree->diagonal,best);
+	}
+	else
+	{
+		if (tree->coarse)
+		{
+			free_image_tree(tree->coarse);
+			free_image_tree(tree->horizontal);
+			free_image_tree(tree->vertical);
+			free_image_tree(tree->diagonal);
+			tree->coarse=tree->horizontal=tree->vertical=tree->diagonal=NULL;
+		}
+	}
+}
+/************************************************************************/
+/*	Functionname:		decompose_to_level                       */
+/* -------------------------------------------------------------------- */
+/*	Parameter:							*/
+/*            	original: original image		 		*/
+/*           	level: decompose to level	 	                */
+/*		flt: decompos with filters[0..level]                   	*/
+/*		method: transform with filter method			*/
+/* -------------------------------------------------------------------- */
+/* 	Description: Decomposes an image to an certain level and stores	*/
+/*	     only this level in the returned quadtree      		*/
+/************************************************************************/
+Image_tree decompose_to_level(Image original,int level,FilterGH *flt,enum FilterMethod method)
+{	Image_tree tree;
+	int e;
+	tree=new_image_tree();
+	tree->image=original;
+	e=decompose_all(tree,level,flt,method,entropy,1);
+	if (!e) return NULL;
+	free_levels(tree,level);
+	return tree;
+}
+/************************************************************************/
+/*	Functionname:		decomposition                           */
+/* -------------------------------------------------------------------- */
+/*	Parameter:							*/
+/*            	t_img: Image which should be decomposed  		*/
+/*           	coarse,horizontal,vertical,diagonal: 	                */
+/*			decomposed images                             	*/
+/*		method: transform with filter method			*/
+/*            	g,h: the transformation is carried out with these filters*/
+/* -------------------------------------------------------------------- */
+/* 	Description: This carries out one wavelettransformation		*/
+/*	     using waveletfilters.            				*/
+/************************************************************************/
+static int decomposition(Image t_img,Image coarse,Image horizontal,Image vertical,
+Image diagonal,Filter g,Filter h,enum FilterMethod method)
+{ 	Image temp1;
+	/*coarse*/
+	temp1=new_image(coarse->width,t_img->height);
+	if(!temp1) goto error;
+convolute_lines(temp1,t_img,h,method);
+	convolute_rows(coarse,temp1,h,method);
+	/*horizontal*/
+	convolute_rows(horizontal,temp1,g,method);
+	free_image(temp1);
+/*vertical*/
+	temp1=new_image(vertical->width,t_img->height);
+	if(!temp1) goto error;
+convolute_lines(temp1,t_img,g,method);
+	convolute_rows(vertical,temp1,h,method);
+/*diagonal*/
+	convolute_rows(diagonal,temp1,g,method);
+	free_image(temp1);
+	return 1;
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  	return 0;
+}
+/************************************************************************/
+/*	Functionname: inv_decomposition					*/
+/* -------------------------------------------------------------------- */
+/*	Parameter: 							*/
+/*		sum: reconstructed image		  		*/
+/*		coarse,horizontal,vertical,diagonal: images to carry out*/
+/*			the inverse transformation			*/
+/*		flt_gh: transform with filters g and h			*/
+/* 		method: transform with filter method			*/
+/* -------------------------------------------------------------------- */
+/* 	Description: Carries out the wavelettransform 			*/
+/*	                           		           		*/
+/************************************************************************/
+static int inv_decomposition(Image sum,Image coarse,Image horizontal,Image vertical,
+Image diagonal,FilterGH flt_gh,enum FilterMethod method)
+{       Image temp1;
+	Filter g,h;
+	if (flt_gh->type==FTOrtho) {
+	  g=flt_gh->g;
+	  h=flt_gh->h;
+	}
+	else {
+	  g=flt_gh->gi;
+	  h=flt_gh->hi;
+	}
+	/*coarse*/
+	temp1=new_image(coarse->width,sum->height);
+	if(!temp1) goto error;
+	convolute_rows(temp1,coarse,h,method);
+/*horizontal*/
+	convolute_rows(temp1,horizontal,g,method);
+convolute_lines(sum,temp1,h,method);
+free_image(temp1);
+/*vertical*/
+temp1=new_image(vertical->width,sum->height);
+	if(!temp1) goto error;
+convolute_rows(temp1,vertical,h,method);
+/*diagonal*/
+convolute_rows(temp1,diagonal,g,method);
+convolute_lines(sum,temp1,g,method);
+free_image(temp1);
+	return 1;
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  	return 0;
+}
+/************************************************************************/
+/*	Functionname: build_image					*/
+/* -------------------------------------------------------------------- */
+/*	Parameter: 							*/
+/*		quadtree: quadtree with decomposition information	*/
+/*	        width,height: image width and height           		*/
+/*		RETURN: returns the build up image			*/
+/* -------------------------------------------------------------------- */
+/* 	Description: builds up an image out of an Image_tree		*/
+/*	                           		           		*/
+/************************************************************************/
+Image build_image(Image_tree quadtree,int width,int height)
+{ 	Image ret_img,coarse,horizontal,vertical,diagonal;
+	ret_img=new_image(width,height);
+	if(!ret_img) goto error;
+	width=(width+1)/2;
+	height=(height+1)/2;
+	if (!(quadtree->image)) {
+	coarse=build_image(quadtree->coarse,width,height);
+	horizontal=build_image(quadtree->horizontal,width,height);
+	vertical=build_image(quadtree->vertical,width,height);
+	diagonal=build_image(quadtree->diagonal,width,height);
+	if (!coarse||!horizontal||!vertical||!diagonal) return NULL;
+	copy_into_image(ret_img,coarse,0,0);
+	  	copy_into_image(ret_img,horizontal,width,0);
+	  	copy_into_image(ret_img,vertical,0,height);
+	  	copy_into_image(ret_img,diagonal,width,height);
+	  	if (!quadtree->coarse->image) free_image(coarse);
+	  	if (!quadtree->horizontal->image) free_image(horizontal);
+	  	if (!quadtree->vertical->image) free_image(vertical);
+	  	if (!quadtree->diagonal->image) free_image(diagonal);
+	return ret_img;
+	}
+	else return quadtree->image;
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  	return NULL;
+}
+/************************************************************************/
+/*	Functionname: inv_transform					*/
+/* -------------------------------------------------------------------- */
+/*	Parameter: 							*/
+/*		tree: tree with decomposition information		*/
+/*		flt_gh: transform with filters g and h			*/
+/* 		method: transform with filter method			*/
+/* -------------------------------------------------------------------- */
+/* 	Description: Inverts the wavelettransform,best_basis,best_level */
+/*	                           		           		*/
+/************************************************************************/
+Image inv_transform(Image_tree tree,FilterGH *flt,
+enum FilterMethod method)
+{	int er,width,height;
+	Image ret_img,coarse,vertical,horizontal,diagonal;
+	if (!tree->image) {
+		coarse=inv_transform(tree->coarse,flt,method);
+	horizontal=inv_transform(tree->horizontal,flt,method);
+	        vertical=inv_transform(tree->vertical,flt,method);
+	        diagonal=inv_transform(tree->diagonal,flt,method);
+	if (!coarse||!horizontal||!vertical||!diagonal) return NULL;
+		width=coarse->width+horizontal->width;
+		height=coarse->height+vertical->height;
+		ret_img=new_image(width,height);
+		if(!ret_img) goto error;
+	        if (tree->flag==0)		/*if flag is set it is a doubletree tiling*/
+	        {
+//			er=inv_decomposition(ret_img,coarse,horizontal,vertical,diagonal,flt[1],method);
+			er=inv_decomposition(ret_img,coarse,horizontal,vertical,diagonal,flt[tree->level],method);
+			if (!er) return NULL;
+	        }
+	  	else
+	  	{
+	  		copy_into_image(ret_img,coarse,0,0);
+	  		copy_into_image(ret_img,horizontal,coarse->width,0);
+	  		copy_into_image(ret_img,vertical,0,coarse->height);
+	  		copy_into_image(ret_img,diagonal,coarse->width,coarse->height);
+	  	}
+		if (!tree->coarse->image) free_image(coarse);
+		if (!tree->horizontal->image) free_image(horizontal);
+		if (!tree->vertical->image) free_image(vertical);
+		if (!tree->diagonal->image) free_image(diagonal);
+		return ret_img;
+	}
+	else return tree->image;
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  	return NULL;
+}
+/************************************************************************/
+/*	Functionname: find_deepest_level				*/
+/* -------------------------------------------------------------------- */
+/*	Parameter: 							*/
+/* -------------------------------------------------------------------- */
+/* 	Description: Finds the deepest possible level where width and 	*/
+/*			height can divided by two exactly.		*/
+/************************************************************************/
+int find_deepest_level(int width,int height)
+{
+	int level=0,w=width,h=height;
+	while ( !((w%2)||(h%2)))
+	{
+		w=w/2;
+		h=h/2;
+		level++;
+	}
+	return level-1;
+}
+/************************************************************************/
+/*	Functionname: convolute_lines					*/
+/* -------------------------------------------------------------------- */
+/*	Parameter: 							*/
+/*		output: output image of wavelettransformation		*/
+/*	        input: input image for decomposition            	*/
+/*		flt: transform with filter flt				*/
+/*		method: transform with filter method			*/
+/* -------------------------------------------------------------------- */
+/* 	Description: Carries out the wavelettransform for all lines of  */
+/*	             the input image					*/
+/************************************************************************/
+static int convolute_lines(Image output,Image input,Filter flt,enum FilterMethod method)
+/*Convolute the lines with filter*/
+{       int i;
+for (i=0;i<input->height;i++) {
+	  switch(method) {
+	  case cutoff:
+filter_cutoff(input,input->width*i,input->width,1,
+			output,output->width*i,output->width,1,flt);
+	    break;
+	  case inv_cutoff:
+filter_inv_cutoff(input,input->width*i,input->width,1,
+			output,output->width*i,output->width,1,flt);
+	    break;
+	  case periodical:
+filter_periodical(input,input->width*i,input->width,1,
+			output,output->width*i,output->width,1,flt);
+	    break;
+	  case inv_periodical:
+filter_inv_periodical(input,input->width*i,input->width,1,
+			output,output->width*i,output->width,1,flt);
+	    break;
+	  case mirror:
+filter_mirror(input,input->width*i,input->width,1,
+			output,output->width*i,output->width,1,flt);
+	    break;
+	  case inv_mirror:
+filter_inv_mirror(input,input->width*i,input->width,1,
+			output,output->width*i,output->width,1,flt);
+	    break;
+	  }
+}
+return 1;
+}
+/************************************************************************/
+/*	Functionname: convolute_rows					*/
+/* -------------------------------------------------------------------- */
+/*	Parameter: 							*/
+/*		output: output image of wavelettransformation		*/
+/*	        input: input image for decomposition            	*/
+/*		flt: transform with filter flt				*/
+/*		method: transform with filter method			*/
+/* -------------------------------------------------------------------- */
+/* 	Description: Carries out the wavelettransform for all rows of   */
+/*	             the input image					*/
+/************************************************************************/
+static int convolute_rows(Image output,Image input,Filter flt,enum FilterMethod method)
+/*Convolute the rows with filter*/
+{       int i;
+for (i=0;i<input->width;i++)
+{
+	switch (method)
+	{
+		case cutoff:
+	    		filter_cutoff(input,i,input->height,input->width,
+	           		output,i,output->height,output->width,flt);
+	    		break;
+	  	case inv_cutoff:
+	    		filter_inv_cutoff(input,i,input->height,input->width,
+	           		output,i,output->height,output->width,flt);
+	    		break;
+	  	case periodical:
+	    		filter_periodical(input,i,input->height,input->width,
+	           		output,i,output->height,output->width,flt);
+	    		break;
+		case inv_periodical:
+			filter_inv_periodical(input,i,input->height,input->width,
+	           		output,i,output->height,output->width,flt);
+	    		break;
+	  	case mirror:
+	    		filter_mirror(input,i,input->height,input->width,
+	           		output,i,output->height,output->width,flt);
+	    		break;
+		case inv_mirror:
+			filter_inv_mirror(input,i,input->height,input->width,
+	           		output,i,output->height,output->width,flt);
+	    		break;
+	  	}
+}
+return 1;
+}
+/************************************************************************/
+/*	Functionname: sumationq						*/
+/* -------------------------------------------------------------------- */
+/*	Parameter: 							*/
+/*		img: image to compute					*/
+/* -------------------------------------------------------------------- */
+/* 	Description: compute the sum of quadrats of all elements of  	*/
+/*	             the input image					*/
+/************************************************************************/
+static Pixel sumationq(Image img)
+{	Pixel sum=0;
+	int i;
+	for (i=0;i<img->size;i++) {
+		sum+=(*img->data+i)*(*img->data+i);
+	}
+	return sum;
+}
+/************************************************************************/
+/*	Functionname: normq						*/
+/* -------------------------------------------------------------------- */
+/*	Parameter: 							*/
+/*		tree: tree to compute					*/
+/* -------------------------------------------------------------------- */
+/* 	Description: computes the quadratic norm over all images in 	*/
+/*	             the input tree					*/
+/************************************************************************/
+static Pixel normq(Image_tree tree)
+{	Pixel sum=0;
+	if (tree->image)
+	{
+		sum=sumationq(tree->image);
+	}
+	else
+	{
+		if (tree->coarse) sum+=normq(tree->coarse);
+		if (tree->horizontal) sum+=normq(tree->horizontal);
+		if (tree->vertical) sum+=normq(tree->vertical);
+		if (tree->diagonal) sum+=normq(tree->diagonal);
+	}
+	return sum;
+}
+/************************************************************************/
+/*	Functionname: sumation_down					*/
+/* -------------------------------------------------------------------- */
+/*	Parameter: 							*/
+/*		tree: tree to compute					*/
+/* 		normq: norm of the images in the tree			*/
+/* -------------------------------------------------------------------- */
+/* 	Description: computes the Entropy over all (string aded) images */
+/*	             in the input tree					*/
+/************************************************************************/
+static Pixel sumation_down(Image_tree tree, Pixel normq)
+{	Pixel sum=0,p;
+	int i;
+	Image img;
+	Pixel *data;
+	if (tree->image)
+	{
+		img=tree->image;
+		data=img->data;
+		for (i=0;i<img->size;i++,data++)
+		{
+			if (*data!=0)
+			{
+				p=(*data)*(*data)/normq;
+				sum+=p*log(1/p);
+			}
+		}
+	}
+	else
+	{
+		if (tree->coarse) sum+=sumation_down(tree->coarse,normq);
+		if (tree->horizontal) sum+=sumation_down(tree->horizontal,normq);
+		if (tree->vertical) sum+=sumation_down(tree->vertical,normq);
+		if (tree->diagonal) sum+=sumation_down(tree->diagonal,normq);
+	}
+	return sum;
+}
+/************************************************************************/
+/*	Functionname: comp						*/
+/* -------------------------------------------------------------------- */
+/* 	Description: used for quicksort for decreasing order		*/
+/************************************************************************/
+int comp(const Pixel *x,const Pixel *y)
+{
+	if (*x<*y) return 1;
+	else if (*x==*y) return 0;
+	else return -1;
+}
+/************************************************************************/
+/*	Functionname: recarea						*/
+/*		tree: Image tree to compute				*/
+/*		list: target list					*/
+/*		list_size: actual size of the list			*/
+/* -------------------------------------------------------------------- */
+/* 	Description: copies all elements within the tree into an list	*/
+/************************************************************************/
+static void recarea(Image_tree tree,Pixel *list,int *list_size)
+{	Image img;
+	if (tree->image)
+	{
+		img=tree->image;
+		memcpy(list+(*list_size),img->data,img->size*sizeof(Pixel));
+		*list_size+=img->size;
+	}
+	else
+	{
+		if (tree->coarse) recarea(tree->coarse,list,list_size);
+		if (tree->horizontal) recarea(tree->horizontal,list,list_size);
+		if (tree->vertical) recarea(tree->vertical,list,list_size);
+		if (tree->diagonal) recarea(tree->diagonal,list,list_size);
+	}
+}
+static void abs_list(Pixel *list,int list_size)
+{
+	int i;
+	for (i=0;i<list_size;i++) list[i]=fabs(list[i]);
+}
+/************************************************************************/
+/*	Functionname: sum_list						*/
+/* -------------------------------------------------------------------- */
+/* 	Description: computes the sum of all poweres list elements	*/
+/************************************************************************/
+static Pixel sum_list(Pixel *list,int p,int size)
+{	Pixel sum=0;
+	int i;
+	for (i=0;i<size;i++)
+	{
+		if (p!=1) sum+=pow(list[i],p);
+		else sum+=list[i];
+	}
+	return sum;
+}
+static Pixel weak_lp(Image_tree tree,int size,int p,double epsilon)
+{	Pixel wlp,*list,max=0;
+	int *list_size,k;
+	list_size=(int *)malloc(sizeof(int));
+	if (!list_size) goto error;
+	*list_size=0;
+	list=(Pixel *)calloc(size,sizeof(Pixel));
+	if (!list) goto error;
+	recarea(tree,list,list_size);
+	abs_list(list,*list_size);
+	qsort(list,*list_size, sizeof(Pixel), (int (*)(const void*, const void*)) comp);
+	for (k=0;k<size;k++)
+	{
+		if (k!=0) wlp=pow(k,1/p)*list[k];
+		else wlp=0;
+		if (wlp>max) max=wlp;
+	}
+	free(list);
+	return max;
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  	return 0;
+}
+static Pixel comp_number(Image_tree tree,int size,int p,double f)
+{	Pixel sum=0,*list,min=MAXDOUBLE,norm,npf,normf;
+	int *list_size=0,k;
+	list_size=(int *)malloc(sizeof(int));
+	if (!list_size) goto error;
+	*list_size=0;
+	list=(Pixel *)calloc(size,sizeof(Pixel));
+	if (!list) goto error;
+	recarea(tree,list,list_size);
+	abs_list(list,*list_size);
+	qsort(list,*list_size, sizeof(Pixel), (int (*)(const void*, const void*)) comp);
+	norm=sum_list(list,p,size);
+	normf=norm*f;
+	for (k=0;k<size;k++)
+	{
+		if (list[k]!=0)
+		{
+			sum+=pow(list[k],p);
+			npf=fabs(sum-normf);
+			if (npf<min) min=npf;
+		}
+	}
+	min=min/norm;
+	free(list);
+	return min;
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  	return 0;
+}
+static Pixel comp_area(Image_tree tree,int size,int p,double f)
+{	Pixel sum=0,*list,norm,vk=0;
+	int *list_size=0,k;
+	list_size=(int *)malloc(sizeof(int));
+	if (!list_size) goto error;
+	*list_size=0;
+	list=(Pixel *)calloc(size,sizeof(Pixel));
+	if (!list) goto error;
+	recarea(tree,list,list_size);
+	abs_list(list,*list_size);
+	qsort(list,*list_size, sizeof(Pixel), (int (*)(const void*, const void*)) comp);
+	norm=sum_list(list,p,size);
+	for (k=0;k<size;k++)
+	{
+		if (list[k]!=0)
+		{
+			vk+=pow(list[k],p);
+			sum+=vk;
+		}
+	}
+	free(list);
+	return (size-sum/norm);
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  	return 0;
+}
+static Pixel compute_sdiscrepancy(Image_tree tree,int size)
+{	Pixel *list,min,max,factor,maximum=0,x;
+	int *list_size=0,k;
+	list_size=(int *)malloc(sizeof(int));
+	if (!list_size) goto error;
+	*list_size=0;
+	list=(Pixel *)calloc(size,sizeof(Pixel));
+	if (!list) goto error;
+	recarea(tree,list,list_size);
+	qsort(list,*list_size, sizeof(Pixel), (int (*)(const void*, const void*)) comp);
+	min=list[0];
+	max=list[size-1];
+	factor=1/(max-min);
+				/*scaling to [0,1]*/
+	for (k=0;k<size;k++)
+	{
+		list[k]=factor*(list[k]-min);
+	}
+	for (k=0;k<size;k++)
+	{
+		x=fabs(list[k]-(2*k-1)/(2*size));
+		if (x>maximum) maximum=x;
+	}
+	free(list);
+	return (1/(2*size)+maximum);
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  	return 0;
+}
+static Pixel compute_discrepancy(Image_tree tree,int size)
+{	Pixel *list,min,max,factor,maximum=0,minimum=0,x;
+	int *list_size=0,k;
+	list_size=(int *)malloc(sizeof(int));
+	if (!list_size) goto error;
+	*list_size=0;
+	list=(Pixel *)calloc(size,sizeof(Pixel));
+	if (!list) goto error;
+	recarea(tree,list,list_size);
+	qsort(list,*list_size, sizeof(Pixel), (int (*)(const void*, const void*)) comp);
+	min=list[0];
+	max=list[size-1];
+	factor=1/(max-min);
+				/*scaling to [0,1]*/
+	for (k=0;k<size;k++)
+	{
+		list[k]=factor*(list[k]-min);
+	}
+	for (k=0;k<size;k++)
+	{
+		x=((Pixel)k/size-list[k]);
+		if (x>maximum) maximum=x;
+		else if (x<minimum) minimum=x;
+	}
+	free(list);
+	return (1/size+maximum-minimum);
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  	return 0;
+}
+static Pixel compute_concentration(Image_tree tree,int size)
+{	Pixel *list,min,max,factor,lkm=0,length,sum=0,value,norm;
+	int *list_size=0,k,next=0,last=0,i=0;
+	list_size=(int *)malloc(sizeof(int));
+	if (!list_size) goto error;
+	*list_size=0;
+	list=(Pixel *)calloc(size+1,sizeof(Pixel));
+	if (!list) goto error;
+	recarea(tree,list,list_size);
+	qsort(list,*list_size, sizeof(Pixel), (int (*)(const void*, const void*)) comp);
+	min=list[0];
+	max=list[size-1];
+	length=(max-min)/100;
+	factor=1/(max-min);
+	for (k=0;k<size;k++)
+	{
+		list[k]=factor*(list[k]-min);
+	}
+	norm=size*sum_list(list,1,size);
+	length=0.01;
+	value=length;
+	list[size]=max+value;
+	for (k=0;k<100;k++)
+	{
+		while ((list[i]<value)&(i<size))
+		{
+			sum+=list[i];
+			next++;
+			i++;
+		}
+		lkm+=(next-last)*sum/norm;
+		value+=length;
+		last=next;
+		sum=0;
+	}
+	return -lkm;
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  	return 0;
+}
+/************************************************************************/
+/*	Functionname: compute_entropy					*/
+/* -------------------------------------------------------------------- */
+/*	Parameter: 							*/
+/*		img: Image from which the entropy should be computed	*/
+/*		cost: choosen costfunction				*/
+/*		epsilon: limit for threshold method			*/
+/* -------------------------------------------------------------------- */
+/* 	Description: computes entropy of an image			*/
+/************************************************************************/
+static double compute_entropy(Image img,enum Information_Cost cost,double epsilon)
+{       double sum=0,x=0;
+	int i;
+	Pixel *data;
+	data=img->data;
+	switch(cost) {
+	case threshold:
+		for(i=0;i<img->size;i++)
+	    		if (fabs(img->data[i])>epsilon) sum++;
+	break;
+	case log_energy:
+	for(i=0;i<img->size;i++,data++) {
+		x=(*data) * (*data);
+	  		if (x!=0) sum+=(x*log(1/x));
+	  	}
+	break;
+	case norml:
+	for(i=0;i<img->size;i++,data++) {
+		x=fabs(*data);
+	  		sum+=x;
+	  	}
+	break;
+	case norml2:
+	for(i=0;i<img->size;i++,data++) {
+		x=(*data) * (*data);
+	  		sum+=x;
+	  	}
+	  	sum=pow(sum,0.5);
+	break;
+	case entropy:
+		for(i=0;i<img->size;i++,data++) {
+	    		x=(*data)*(*data);
+	    		if (x!=0) sum-=(x*log(x));
+	  	}
+	break;
+	case gauss_markov:
+		for(i=0;i<img->size;i++) {
+	    		x=(img->data[i])*(img->data[i]);
+		if (x!=0) sum+=log(x*x);
+		}
+	break;
+}
+return sum;
+}
+/************************************************************************/
+/*	Functionname: compute_non_additive				*/
+/* -------------------------------------------------------------------- */
+/*	Parameter: 							*/
+/*		tree: Image tree from which the entropy should be 	*/
+/*			computed					*/
+/*		size :	size of the image				*/
+/*		cost: choosen costfunction				*/
+/*		epsilon: limit for threshold method			*/
+/*		down: decides if only the first image should be computed*/
+/* -------------------------------------------------------------------- */
+/* 	Description: computes entropy of an image			*/
+/************************************************************************/
+static Pixel compute_non_additive(Image_tree tree,int size,enum Information_Cost cost,double epsilon,int down)
+{	Pixel sum=0,normx;
+	Image img=NULL;
+	if (down)
+	{
+		img=tree->image;
+		tree->image=NULL;
+	}
+	switch (cost)
+	{
+		case shanon:
+			normx=normq(tree);
+			sum=-sumation_down(tree,normx);
+			break;
+		case weak_l:
+			sum=weak_lp(tree,size,1,epsilon);
+			break;
+		case weak_lq:
+			sum=weak_lp(tree,size,2,epsilon);
+			break;
+		case compression_number:
+			sum=comp_number(tree,size,1,epsilon);
+			break;
+		case compression_numberq:
+			sum=comp_number(tree,size,2,epsilon);
+			break;
+		case compression_area:
+			sum=comp_area(tree,size,1,epsilon);
+			break;
+		case compression_areaq:
+			sum=comp_area(tree,size,2,epsilon);
+			break;
+		case discrepancy:
+			sum=compute_discrepancy(tree,size);
+			break;
+		case sdiscrepancy:
+			sum=compute_sdiscrepancy(tree,size);
+			break;
+		case concentration:
+			sum=compute_concentration(tree,size);
+			break;
+	}
+	if (down) tree->image=img;
+	return sum;
+}
+int rec_double(Image_tree dtree,int level,FilterGH *flt,enum FilterMethod method,enum Information_Cost cost,double epsilon)
+{	int min,width,height;
+	double sum=0;
+	Image c,h,v,d;
+	dtree->level=0;
+	if (cost>=shanon)
+	{
+		dtree->entropy=compute_non_additive(dtree,dtree->image->size,cost,epsilon,0);
+	}
+	else 	dtree->entropy=compute_entropy(dtree->image,cost,epsilon);
+	dtree->doubletree=best_basis(dtree->image,level,flt,method,cost,epsilon);
+	min=dtree->image->width;
+	if (dtree->image->height<min) min=dtree->image->height;
+	if (doubletree_min<min)
+	{
+		width=(dtree->image->width+1)/2;
+		height=(dtree->image->height+1)/2;
+		dtree->coarse=new_image_tree();
+		dtree->horizontal=new_image_tree();
+		dtree->vertical=new_image_tree();
+		dtree->diagonal=new_image_tree();
+		c=new_image(width,height);
+		h=new_image(width,height);
+		v=new_image(width,height);
+		d=new_image(width,height);
+		if(!c||!h||!v||!d) goto error;
+		copy_part_of_image(c,dtree->image,0,0);
+		copy_part_of_image(h,dtree->image,width,0);
+		copy_part_of_image(v,dtree->image,0,height);
+		copy_part_of_image(d,dtree->image,width,height);
+		dtree->coarse->image=c;
+		dtree->horizontal->image=h;
+		dtree->vertical->image=v;
+		dtree->diagonal->image=d;
+		rec_double(dtree->coarse,level,flt,method,cost,epsilon);
+		rec_double(dtree->horizontal,level,flt,method,cost,epsilon);
+		rec_double(dtree->vertical,level,flt,method,cost,epsilon);
+		rec_double(dtree->diagonal,level,flt,method,cost,epsilon);
+		/* going back in recursion*/
+		sum=dtree->coarse->entropy+dtree->horizontal->entropy+
+			dtree->vertical->entropy+dtree->diagonal->entropy;
+		if (sum>dtree->entropy)
+		{
+			/*take image*/
+			free_image_tree(dtree->coarse);
+			free_image_tree(dtree->horizontal);
+			free_image_tree(dtree->vertical);
+			free_image_tree(dtree->diagonal);
+			dtree->coarse=dtree->horizontal=dtree->vertical=dtree->diagonal=NULL;
+		}
+		else
+		{	/*take tiling*/
+			dtree->entropy=sum;
+			free_image(dtree->image);
+			dtree->image=NULL;
+		}
+		if (dtree->entropy>dtree->doubletree->entropy)
+		{
+			/*take best basis tree*/
+			dtree->entropy=dtree->doubletree->entropy;
+			if(dtree->coarse) free_image_tree(dtree->coarse);
+			if(dtree->horizontal) free_image_tree(dtree->horizontal);
+			if(dtree->vertical) free_image_tree(dtree->vertical);
+			if(dtree->diagonal) free_image_tree(dtree->diagonal);
+			dtree->coarse=dtree->doubletree->coarse;
+			dtree->horizontal=dtree->doubletree->horizontal;
+			dtree->vertical=dtree->doubletree->vertical;
+			dtree->diagonal=dtree->doubletree->diagonal;
+			free_image(dtree->image);
+			dtree->image=NULL;
+			free(dtree->doubletree);
+			dtree->doubletree=NULL;
+		}
+		else
+		{
+			dtree->flag=1;
+			if(dtree->doubletree) free_image_tree(dtree->doubletree);
+			dtree->doubletree=NULL;
+		}
+	}
+	return 1;
+	error:
+	err_SimpleMessage(err_GetErrorMessage(Error_NotEnoughMemory));
+	  	return 0;
+}
+static void save_structur(Image_tree tree,FILE *fp,int pos)
+{
+	int shift,next_pos,max;
+	if (tree->flag)
+	{
+		fprintf(fp,"%d ",pos);
+		shift=pos-(pow(4,tree->level-1)-1)*4/3-1;
+		max=(int) ((pow(4,tree->level)-1)*4/3);
+		next_pos=max+4*shift+1;
+		if (tree->coarse) save_structur(tree->coarse,fp,next_pos);
+		if (tree->horizontal) save_structur(tree->horizontal,fp,next_pos+1);
+		if (tree->vertical) save_structur(tree->vertical,fp,next_pos+2);
+		if (tree->diagonal) save_structur(tree->diagonal,fp,next_pos+3);
+	}
+}
+static int is_in_list(int *list,int len, int x)
+{
+	int i,found=0;
+	for (i=0;i<len;i++)
+	{
+		if (list[i]==x)
+		{
+			found=1;
+			i=len;
+		}
+	}
+	return found;
+}
+static void write_flags(Image_tree tree,int *list,int len,int pos)
+{
+	int shift,next_pos,max;
+	if (is_in_list(list,len,pos))
+	{
+		tree->flag=1;
+		shift=pos-(pow(4,tree->level-1)-1)*4/3-1;
+		max=(int) ((pow(4,tree->level)-1)*4/3);
+		next_pos=max+4*shift+1;
+		write_flags(tree->coarse,list,len,next_pos);
+		write_flags(tree->horizontal,list,len,next_pos+1);
+		write_flags(tree->vertical,list,len,next_pos+2);
+		write_flags(tree->diagonal,list,len,next_pos+3);
+	}
+}
+/************************************************************************/
+/*	Functionname:		err_simple_message			*/
+/* -------------------------------------------------------------------- */
+/*	Parameter:							*/
+/*	    char *: string that contains information about an		*/
+/*		    error the user should know.				*/
+/* -------------------------------------------------------------------- */
+/* 	Description:							*/
+/*	    Prints error messages for the user.				*/
+/************************************************************************/
+void err_SimpleMessage(char *message)
+{
+	fprintf(stderr,"%s\n",message);
+}
+/************************************************************************/
+/*	Functionname: 		err_get_message				*/
+/* -------------------------------------------------------------------- */
+/* 	Return value:   Errormessage for this specific error.		*/
+/*	Parameter:							*/
+/*	    Error err:	Error whose errormessage should be returned	*/
+/* -------------------------------------------------------------------- */
+/*	Description:							*/
+/************************************************************************/
+char * err_GetErrorMessage(Error err)
+{
+switch (err)
+{
+case Error_NotImplemented:
+return "Sorry, this is not implemented yet. ";
+break;
+case Error_AssertionFailed:
+return "Sorry, an internal assertion was violated.\n"
+"This action can not be completed. :-(";
+break;
+case Error_NotEnoughMemory:
+return "Sorry, there is not enough memory";
+break;
+case Error_Limitation:
+return "Some limitation of the program exceeded";
+break;
+	/* - FILES - */
+	    case Error_CantOpenFile:
+		return "Could not open file";
+		break;
+	    case Error_CantCreateFile:
+		return "Could not create file";
+		break;
+	    case Error_CantCloseFile:
+		return "Could not close file";
+		break;
+	    case Error_InternalError:
+return "Sorry, an internal error occured.\n"
+"This action can not be completed. :-(";
+break;
+default:
+return "Sorry, but an unknown error ocurred.\n"
+"This action can not be completed. :-(";
+break;
+	}
+}

Mercurial > hg > wm

comparison Meerwald-dir/wavelet.c @ 24:9f20bce6184e v0.7