resample.cpp

基于小波的图像配准

#include "resample.h"

/* a local function, just called by resample function */
long c_arrayRead(unsigned char *image_in,
			     long img_width,long img_height,
                 long xBegin,long yBegin,long xLen,long yLen,
			     unsigned char *stenArray);
long d_arrayRead(double *image_in, 
				 long img_width,long img_height,
		         long xBegin,long yBegin,long xLen,long yLen,
			     double *stenArray);

/* This is the resample function for affine transform */
/*
 * 1. image_in and image_out are source and target image data
 * 2. img_width and img_height are image size;
 * 3. trans_x and trans_y are transform parameters
 * 4. angle is rotated angle
 * 5. interface starting by "c" means this function is for char* dataArray;
 *    while statrting by "d" means it is for double* dataArray.
 */

void c_resample(unsigned char *image_in, unsigned char *image_out,
			    long img_width, long img_height,
                double trans_tx, double trans_ty,double angle)
{
	long xBe,yBe,xLen,yLen;                       /* template data size */
	unsigned char *stenArray,*st;                 /* template data */
	int tx,ty;                                    /* index */
	double max_u,min_u,max_v,min_v;               /* target image domain */
	double tmp_u,tmp_v;                           /* tmp var */
	double x,y;                                   /* x,y var */
	double x_fra,y_fra;                           /* x,y fractal part */
	long x_int,y_int;                             /* x,y integer part */
	long f_x,f_y;                                 /* filter x,y */
	double tab_gap;                               /* the tab gap */
	int filterType = CUBIC_TABLE_FILTER;          /* filter type */
	double *cubicTable,*ct;                       /* filter template */
	size_t mem_size;                              /* temp variable */
	int cw,ch,tn;                                 /* temp varaible */
	char pixel;                                   /* pixel value */
	double pixel_d =0.;                           /* temp var */
	long pixels_num = 0;                          /* resampled pixels number */

//=================================================================================
//      the works related to filter
	mem_size = CUBIC_WIDTH * CUBIC_HEIGHT * TAB_NUM * TAB_NUM * sizeof(double);
	cubicTable = (double *)malloc(mem_size);
	tab_gap = 1. / TAB_NUM;
    cubicFunctionTable(CUBIC_WIDTH,CUBIC_HEIGHT,TAB_NUM,cubicTable);

    xLen = CUBIC_WIDTH;
    yLen = CUBIC_HEIGHT;	      
    stenArray = (unsigned char *)malloc(xLen*yLen);
//================================================================================
//     Main Loop
    mem_size = CUBIC_WIDTH * CUBIC_HEIGHT;
	for(ty = 0; ty < img_height;ty++ ) {
		for(tx = 0;tx <img_width;tx++ ) {
                //+++++++++++++++ 
                // The inverse function
		x = cos( angle ) * (tx - trans_tx - img_width/2 ) 
                   +sin( angle ) * (ty - trans_ty - img_height/2) + img_width/2;
		y =-sin( angle ) * (tx - trans_tx - img_width/2 ) 
                   +cos( angle ) * (ty - trans_ty - img_height/2) + img_height/2;
                //+++++++++++++++

		x_int = (int) x;
		y_int = (int) y;
		x_fra = x - x_int;
		y_fra = y - y_int;
		f_x = (int)(x_fra/tab_gap);
		f_y = (int)(y_fra/tab_gap);   
		ct = cubicTable + (f_y*TAB_NUM+f_x)*mem_size;
	    if( (x_int >0) && (x_int < (img_width-2)) && (y_int > 0)
			&& (y_int < (img_height-2)) ) {
			c_arrayRead(image_in,img_width,img_height,x_int-1,
				      y_int-1,xLen,yLen,stenArray);
			st = stenArray;
			pixel_d = 0.;
			for(ch = 0;ch <yLen;ch++) {
				for(cw = 0;cw <xLen;cw++) { 

					pixel_d += (*ct)*(*st);
					ct++;
					st++;
				}
			}
			if(pixel_d >255.) { pixel_d = 255.; }
			if(pixel_d <0. )  { pixel_d = 0.; }
			pixel = (char)pixel_d;
			pixels_num++;
            *(image_out+ty*img_width+tx) = pixel;
			//imgWrite(tfPtr,tx,ty,1,1,&pixel);
			
		}
	  } /* for(tx... */
	} /* for(ty.. */
       free(cubicTable);
       free(stenArray);

}

void d_resample(double *image_in, double *image_out,
			    long img_width, long img_height,
                double trans_tx, double trans_ty,double angle)
{
	long xBe,yBe,xLen,yLen;                       /* template data size */
	double *stenArray,*st;                        /* template data */
	int tx,ty;                                    /* index */
	double max_u,min_u,max_v,min_v;               /* target image domain */
	double tmp_u,tmp_v;                           /* tmp var */
	double x,y;                                   /* x,y var */
	double x_fra,y_fra;                           /* x,y fractal part */
	long x_int,y_int;                             /* x,y integer part */
	long f_x,f_y;                                 /* filter x,y */
	double tab_gap;                               /* the tab gap */
	int filterType = CUBIC_TABLE_FILTER;          /* filter type */
	double *cubicTable,*ct;                       /* filter template */
	size_t mem_size;                              /* temp variable */
	int cw,ch,tn;                                 /* temp varaible */
	double pixel_d =0.;                           /* temp var */
	long pixels_num = 0;                          /* resampled pixels number */

//=================================================================================
//      the works related to filter
	mem_size = CUBIC_WIDTH * CUBIC_HEIGHT * TAB_NUM * TAB_NUM * sizeof(double);
	cubicTable = (double *)malloc(mem_size);
	tab_gap = 1. / TAB_NUM;
    cubicFunctionTable(CUBIC_WIDTH,CUBIC_HEIGHT,TAB_NUM,cubicTable);

    xLen = CUBIC_WIDTH;
    yLen = CUBIC_HEIGHT;	      
    stenArray = (double *)malloc(sizeof(double)*(xLen*yLen));
//================================================================================
//     Main Loop
    mem_size = CUBIC_WIDTH * CUBIC_HEIGHT;
	for(ty = 0; ty < img_height;ty++ ) {
		for(tx = 0;tx <img_width;tx++ ) {
                //+++++++++++++++ 
                // The inverse function
		x = cos( angle ) * (tx - trans_tx - img_width/2 ) 
                   +sin( angle ) * (ty - trans_ty - img_height/2) + img_width/2;
		y =-sin( angle ) * (tx - trans_tx - img_width/2 ) 
                   +cos( angle ) * (ty - trans_ty - img_height/2) + img_height/2;
                //+++++++++++++++

		x_int = (int) x;
		y_int = (int) y;
		x_fra = x - x_int;
		y_fra = y - y_int;
		f_x = (int)(x_fra/tab_gap);
		f_y = (int)(y_fra/tab_gap);   
		ct = cubicTable + (f_y*TAB_NUM+f_x)*mem_size;
	    if( (x_int >0) && (x_int < (img_width-2)) && (y_int > 0)
			&& (y_int < (img_height-2)) ) {
			d_arrayRead(image_in,img_width,img_height,x_int-1,
						y_int-1,xLen,yLen,stenArray);
			st = stenArray;
			pixel_d = 0.;
			for(ch = 0;ch <yLen;ch++) {
				for(cw = 0;cw <xLen;cw++) { 

					pixel_d += (*ct)*(*st);
					ct++;
					st++;
				}
			}
			if(pixel_d >255.) { pixel_d = 255.; }
			if(pixel_d <0. )  { pixel_d = 0.; }
			pixels_num++;
            *(image_out+ty*img_width+tx) = pixel_d;
			//imgWrite(tfPtr,tx,ty,1,1,&pixel);
			
		}
	  } /* for(tx... */
	} /* for(ty.. */
       free(cubicTable);
       free(stenArray);

}

/*
 * local procedure
 */
long c_arrayRead(unsigned char *image_in, 
				 long img_width,long img_height,
		         long xBegin,long yBegin,long xLen,long yLen,
			     unsigned char *stenArray)
{
	long ty,tx;
    image_in +=  yBegin*img_width + xBegin;
	for(ty = 0; ty <yLen ;ty++) {
        for(tx = 0; tx <xLen;tx++) {
			*stenArray = *image_in;
            image_in++;
            stenArray++;            
        }
            image_in += img_width - xLen;
	}
	return 0;
}

long d_arrayRead(double *image_in, 
				 long img_width,long img_height,
		         long xBegin,long yBegin,long xLen,long yLen,
			     double *stenArray)
{
	long ty,tx;
    image_in +=  yBegin*img_width + xBegin;
	for(ty = 0; ty <yLen ;ty++) {
        for(tx = 0; tx <xLen;tx++) {
			*stenArray = *image_in;
            image_in++;
            stenArray++;            
        }
            image_in += img_width - xLen;
	}
	return 0;
}