wavelet.c

包含了多个matlab编程在图像中加入水印的处理代码

/* 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, *f_data;
	int fstart,fend; /* Source interval */
	int in_pos,in_dir,in_div,in_mod;

	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 compute_levels(Image_tree tree,double *entropies,enum Information_Cost cost,double epsilon);
static 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 			*/
/*	                           		           		*/
/************************************************************************/
extern 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;