wavelet.c

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

	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		*/
/*	                           		           		*/
/************************************************************************/
extern 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 */
/*	                           		           		*/
/************************************************************************/
extern 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.		*/
/************************************************************************/
extern 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, siz