resample.c

本文档包括了几何矫正的主体内容

/*
 * The program resample images got from remote rense
 * 
 * 1. we use least quadratic multipfy to compute inverse 
 * function of function 'f';
 * 
 * 2. we open image files for input and output
 *
 * 3. Resampling loop use  target space ordinate (i,j) and inverse 
 * function to get input image space ordinate (x,y),then use 
 * convolution template to get grayscale of (x,y) location.
 */
#include <resample.h>
main(argc,argv)
	int argc;
	char **argv;
{
	double *poCoefArray_x,*poCoefArray_y;        /* function coeficient */
       	double *neCoefArray_u,*neCoefArray_v;        /* inverse coef */
	double *uvMaxMin;                            /* u,v max min */
	char sourcefilename[100],targetfilename[100]; /* filename */
	FILE *sfPtr,*tfPtr;                           /* file ptr */
	long img_width,img_height;                    /* source file size */
	long img_twidth,img_theight;                  /* target file size */
	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 */
	unsigned char pixel;                                   /* pixel value */
	unsigned char pixel_max,pixel_min;                     /* max min value of pixels*/
	double pixel_d =0.;                           /* temp var */
	char *tImage,*tL;                            /* target data */
	char *sImage,*sL;                            /* source data */
	long pixels_num = 0;                          /* resampled pixels number */

	DECLARE_TIME(t0);
	DECLARE_TIME(t1);
	DECLARE_TIME(tA);
	DECLARE_TIME(tB);
	memset(timer_ticks,0,VPTIMER_COUNT);
	if(argc < 3) {
		printf("Usage: resample inputfile outputfile [filterType]\n");
		exit(1);
	}
	strcpy(sourcefilename,*(++argv));
	strcpy(targetfilename,*(++argv));
	argc-=3;
	if (argc > 0) {
		filterType = atoi(*(++argv));
		if(( filterType < 0 )||(filterType >= FILTER_NUM)) {
			filterType = CUBIC_TABLE_FILTER;
		}
	}
#if 0
	printf("%s\n%s\n%d\n",sourcefilename,targetfilename,filterType);
#endif
#if 1   
	mem_size = CUBIC_WIDTH * CUBIC_HEIGHT * TAB_NUM * TAB_NUM * sizeof(double);
	cubicTable = (double *)malloc(mem_size);
	tab_gap = 1. / TAB_NUM;
	GET_TIME(t0);
        cubicFunctionTable(CUBIC_WIDTH,CUBIC_HEIGHT,TAB_NUM,cubicTable);
	GET_TIME(t1);
	STORE_TIME(CUBIC_TABLE,t0,t1);
	printf("timer compute cubic filter table: %d\n",timer_ticks[CUBIC_TABLE]);
#endif
#ifdef BMP_FORMAT
       sfPtr = imgOpen(sourcefilename,&img_width,&img_height,0);
#elif defined(PPM_FORMAT)
       sfPtr = imgPPMOpen(sourcefilename,&img_width,&img_height,0);
#endif
       /* compute inverse function use least multify */
       GET_TIME(t0);
       poCoefArray_x = (double *)malloc(10*sizeof(double));
       poCoefArray_y = (double *)malloc(10*sizeof(double));
       poCoefArray_x[0]=2.0434e+7;
       poCoefArray_x[1]=1.0188;
       poCoefArray_x[2]=0.21283;			       
       poCoefArray_x[3]=9.3168e-7;				
       poCoefArray_x[4]=6.3878e-7;				
       poCoefArray_x[5]=-1.1824e-6;				
       poCoefArray_x[6]=-4.1159e-12;				
       poCoefArray_x[7]=1.5281e-11;				
       poCoefArray_x[8]=2.5996e-11;				
       poCoefArray_x[9]=-1.2153e-11;

       poCoefArray_y[0]=4.4623e+6;
	poCoefArray_y[1]=-0.19111;
        poCoefArray_y[2]=1.0035;
	poCoefArray_y[3]=-4.6184e-7;				
	poCoefArray_y[4]=-2.2211e-7;				
	poCoefArray_y[5]=6.595e-7;				
	poCoefArray_y[6]=3.7652e-12;				
	poCoefArray_y[7]=-1.0826e-11;				
	poCoefArray_y[8]=-1.3948e-11;				
	poCoefArray_y[9]=8.622e-12;																
       neCoefArray_u = (double *)malloc(10*sizeof(double));
       neCoefArray_v = (double *)malloc(10*sizeof(double));
       uvMaxMin = (double *)malloc(4*sizeof(double));
       leastMult(0,img_width,-img_height,0,COARSE_X,COARSE_Y,poCoefArray_x,poCoefArray_y,neCoefArray_u,neCoefArray_v,uvMaxMin);
       GET_TIME(t1);
       STORE_TIME(LEAST_MULT,t0,t1);
       printf("least timer %d us\n",timer_ticks[LEAST_MULT]);
       
#if 0
       prtDoubleArray("neCoefArray_u",10,neCoefArray_u);
       prtDoubleArray("neCoefArray_v",10,neCoefArray_v);
#endif
       xLen = CUBIC_WIDTH;
       yLen = CUBIC_HEIGHT;	      
       stenArray = (char *)malloc(xLen*yLen);
#if 0
       for(ty = 0;ty < img_height;ty+=yLen) {
               for(tx = 0;tx < img_width ;tx+= xLen) {		      
                      imgRead(sfPtr,tx,ty,xLen,yLen,stenArray);
		      imgWrite(tfPtr,tx,ty,xLen,yLen,stenArray);
		      
	       }
       }
#endif
       /* Attention: Now I donn't know the forward
	* So I suppose the forward transform:
	*  |  cos(a)  sin(a) |
	*  | -sin(a)  cos(a) |
	*  the inverse transform:
	*  |  cos(a)  -sin(a) |
	*  |  sin(a)   cos(a) |
	*
	*/
#if 1
#define PI3   0.785398163
       max_u = -1000000000000000.;
       max_v = -1000000000000000.;
       min_u =  1000000000000000.;
       min_v =  1000000000000000.;
       tx = 0;
       for(ty = 0;ty < img_height;ty++) {
	   /*tmp_u = cos( PI3 )*tx +sin( PI3 )*(-ty);
           tmp_v = -sin( PI3 )*tx+cos( PI3 )*(-ty);*/
	   tmp_u = Calculate_u(tx,-ty,poCoefArray_x);
	   tmp_v = Calculate_v(tx,-ty,poCoefArray_y);
	   if (tmp_u >max_u) max_u = tmp_u;
	   if (tmp_u <min_u) min_u = tmp_u;
	   if (tmp_v >max_v) max_v = tmp_v;
	   if (tmp_v <min_v) min_v = tmp_v;
       }
       tx = img_width;
       for(ty = 0;ty < img_height;ty++) {
	   /*tmp_u = cos( PI3 )*tx +sin( PI3 )*(-ty);
	   tmp_v = -sin( PI3 )*tx+cos( PI3 )*(-ty);*/
	   tmp_u = Calculate_u(tx,-ty,poCoefArray_x);
	   tmp_v = Calculate_v(tx,-ty,poCoefArray_y);
	   if (tmp_u >max_u) max_u = tmp_u;					
	   if (tmp_u <min_u) min_u = tmp_u;					
	   if (tmp_v >max_v) max_v = tmp_v;					
	   if (tmp_v <min_v) min_v = tmp_v;					
       }
      ty = 0;
      for(tx = 0;tx <img_width;tx++ ) {
           /*tmp_u = cos( PI3 )*tx +sin( PI3 )*ty;
	   tmp_v = -sin( PI3 )*tx+cos( PI3 )*ty;*/
	   tmp_u = Calculate_u(tx,ty,poCoefArray_x);
	   tmp_v = Calculate_v(tx,ty,poCoefArray_y);
	   if (tmp_u >max_u) max_u = tmp_u;			          
	   if (tmp_u <min_u) min_u = tmp_u;					
	   if (tmp_v >max_v) max_v = tmp_v;
	   if (tmp_v <min_v) min_v = tmp_v;	
      } 
      ty = -img_height;
      for(tx = 0;tx <img_width;tx++ ) {
	      /*tmp_u = cos( PI3 )*tx +sin( PI3)*ty;
	      tmp_v = -sin( PI3 )*tx+cos( PI3 )*ty;*/
	      tmp_u = Calculate_u(tx,ty,poCoefArray_x);
	      tmp_v = Calculate_v(tx,ty,poCoefArray_y);
	      if (tmp_u >max_u) max_u = tmp_u;					      	      if (tmp_u <min_u) min_u = tmp_u;
	      if (tmp_v >max_v) max_v = tmp_v;
	      if (tmp_v <min_v) min_v = tmp_v;
      }
#if 1
      printf("max_u = %f max_v = %f min_u = %f min_v = %f\n",max_u,max_v,min_u,min_v);
#endif
      img_theight = (int)((max_v-min_v));
      img_twidth  = (int)((max_u-min_u));
      max_u -= min_u;
      max_u /= img_twidth;
      max_v -= min_v;
      max_v /= img_theight;
#if 1
      printf("max_u = %f max_v = %f min_u = %f min_v = %f\n",max_u,max_v,min_u,min_v);
#endif
#endif 
      /* image size of resampled image */
#ifdef BMP_FORMAT
      tfPtr = imgOpen(targetfilename,&img_twidth,&img_theight,1);
#elif defined(PPM_FORMAT)
      tfPtr = imgPPMOpen(targetfilename,&img_twidth,&img_theight,1);
#endif
			   
       sImage = (char *)malloc(img_width*img_height);
       tImage = (char *)malloc(img_twidth*img_theight);
#ifdef BMP_FORMAT
       imgRead(sfPtr,0,0,img_width,img_height,sImage);
#elif defined(PPM_FORMAT)
       imgPPMRead(sfPtr,0,0,img_width,img_height,sImage);
#endif
       
       GET_TIME(t0);
       mem_size = CUBIC_WIDTH * CUBIC_HEIGHT;
       pixel_max = 0;
       pixel_min = 255;
#if 1
       printf("Begin Resample Loop\n");
#endif
       for(ty = 0; ty < img_theight;ty++ ) {
	  for(tx = 0;tx <img_twidth;tx++ ) {
		/*printf("now pixels ty x tx: %d x %d\n",ty,tx);*/
		tmp_u =  tx * max_u + min_u;
		tmp_v =  ty * max_v + min_v;
		x = Calculate_x(tmp_u,tmp_v,neCoefArray_u,uvMaxMin);
		y = Calculate_y(tmp_u,tmp_v,neCoefArray_v,uvMaxMin);
		/*x = cos( PI3 ) * tmp_u - sin( PI3 ) * tmp_v;
		y = sin( PI3 ) * tmp_u + cos( PI3 ) * tmp_v;*/
					       
 
		x_int = (int) x;
		y_int = (int) y -1;
		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 < -2) && (y_int > (-img_height+1)) ) {
#ifdef BMP_FORMAT
			imgGet(sImage,x_int-1,img_height + y_int - 1,xLen,yLen,img_width,img_height,stenArray);
#elif defined(PPM_FORMAT)
			imgPPMGet(sImage,x_int-1,img_height + y_int - 1,xLen,yLen,img_width,img_height,stenArray);
#endif
			st = stenArray;
			pixel_d = 0.;
#if 1
			for(ch = 0;ch <yLen;ch++) {
				for(cw = 0;cw <xLen;cw++) { 
#if 0
					if( (x_int == 128)&&(y_int==128) ) {
					    printf("%f %d\n",(*ct),(*st));
					    printf("%f ----\n",(*ct)*(*st));
					}
#endif
	
					pixel_d += (*ct)*(*st);
					ct++;
					st++;
				}
			}
			if(pixel_d > 255.) pixel_d = 255.;
			if(pixel_d < 0.) pixel_d = 0.;
#endif
			pixels_num++;
			
		} else {
			pixel_d = 0.;
		}
		pixel = (unsigned char )pixel_d;
		if(pixel > pixel_max) pixel_max = pixel;
		if(pixel < pixel_min) pixel_min = pixel;
#ifdef BMP_FORMAT
		imgPut(tImage,tx,ty,1,1,img_twidth,img_theight,&pixel);
#elif defined(PPM_FORMAT)
		imgPPMPut(tImage,tx,ty,1,1,img_twidth,img_theight,&pixel);
#endif

	  } /* for(tx... */
       } /* for(ty.. */
#if 1
       printf("End Resample Loop\n");
#endif
       GET_TIME(t1);
       STORE_TIME(RESAMPLING_TIMER,t0,t1);
       printf("pixel max min %d %d\n", pixel_max,pixel_min);
#ifdef BMP_FORMAT
       imgWrite(tfPtr,0,0,img_twidth,img_theight,tImage);
#elif defined(PPM_FORMAT)
       imgPPMWrite(tfPtr,0,0,img_twidth,img_theight,tImage);
#endif
       printf("resampling time %d us\n",timer_ticks[RESAMPLING_TIMER]);
       printf("invere function time %d us\n",timer_ticks[INVERSE_FUNCTION]);
       printf("resampled pixels = %d\n",pixels_num);
       free(cubicTable);
       free(stenArray);
#ifdef BMP_FORMAT
       imgClose(sfPtr);
       imgClose(tfPtr);
#elif defined(PPM_FORMAT)
       imgPPMClose(sfPtr);
       imgPPMClose(tfPtr);
#endif
}