least.c

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

#include <least.h>
double MaxValue(double *array,int x_num,int y_num)
{
	double *temp,result;
	int i;
	temp=(double *)malloc(4*sizeof(double));
	*temp=*array;
	*(temp+1)=*(array+x_num-1);
	*(temp+2)=*(array+(y_num-1)*x_num);
	*(temp+3)=*(array+y_num*x_num-1);
	result=*temp;
	for (i=1;i<4;i++)
	{
		if (result<(*(temp+i)))
			result=*(temp+i);
	}
	free(temp);
	return result;
}

double MinValue(double *array,int x_num,int y_num)
{
	double *temp,result;
	int i;
	temp=(double *)malloc(4*sizeof(double));
	*temp=*array;
	*(temp+1)=*(array+x_num-1);
	*(temp+2)=*(array+(y_num-1)*x_num);
	*(temp+3)=*(array+y_num*x_num-1);
	result=*temp;
	for (i=1;i<4;i++)
	{
		if (result>(*(temp+i)))
			result=*(temp+i);
	}
	free(temp);
	return result;
}


void q(double X,double Y,double *result)
{
	result[0]=1;
	result[1]=X;
	result[2]=Y;
	result[3]=0.5*(3*X*X-1);
	result[4]=X*Y;
	result[5]=0.5*(3*Y*Y-1);
	result[6]=0.5*X*(5*X*X-3);
	result[7]=0.5*Y*(3*X*X-1);
	result[8]=0.5*(3*Y*Y-1)*X;
	result[9]=0.5*Y*(5*Y*Y-3);
}

void leastMult(int startX,int endX,int startY,int endY,int coarseX, int coarseY, double* poCoefArray_x, double* poCoefArray_y, double* neCoefArray_u, double* neCoefArray_v,double *uvMaxMin)
{
	int *x,*y;
	double *u,*v;
	int x_num,y_num,i,j;
	long xy_num;
	x_num=(int)((endX-startX)/coarseX)+1;
	y_num=(int)((endY-startY)/coarseY)+1;
	x=(int *)malloc(x_num*sizeof(int));
	y=(int *)malloc(y_num*sizeof(int));
	xy_num=x_num*y_num;

	for (i=0;i<x_num;i++)
		x[i]=i*coarseX + startX;
	for (j=0;j<y_num;j++)
		y[j]=j*coarseY + startY;

	u=(double *)malloc(xy_num*sizeof(double));
    v=(double *)malloc(xy_num*sizeof(double));
	
	for (j=0;j<y_num;j++)
	{
		for (i=0;i<x_num;i++)
		{
		    u[j*x_num+i]=x[i]*(x[i]*(x[i]*poCoefArray_x[6]+poCoefArray_x[3])+poCoefArray_x[1])+poCoefArray_x[0];
			u[j*x_num+i]+=y[j]*(y[j]*(y[j]*poCoefArray_x[9]+poCoefArray_x[5])+poCoefArray_x[2]);
			u[j*x_num+i]+=x[i]*y[j]*(poCoefArray_x[4]+poCoefArray_x[7]*x[i]+poCoefArray_x[8]*y[j]);

			v[j*x_num+i]=x[i]*(x[i]*(x[i]*poCoefArray_y[6]+poCoefArray_y[3])+poCoefArray_y[1])+poCoefArray_y[0];
			v[j*x_num+i]+=y[j]*(y[j]*(y[j]*poCoefArray_y[9]+poCoefArray_y[5])+poCoefArray_y[2]);
			v[j*x_num+i]+=x[i]*y[j]*(poCoefArray_y[4]+poCoefArray_y[7]*x[i]+poCoefArray_y[8]*y[j]);
		}
	}
	long t;
	uvMaxMin[0]=MaxValue(u,x_num,y_num);
	uvMaxMin[1]=MaxValue(v,x_num,y_num);
	uvMaxMin[2]=MinValue(u,x_num,y_num);
    uvMaxMin[3]=MinValue(v,x_num,y_num);
    
	double *X,*Y;
	X=(double *)malloc(xy_num*sizeof(double));
	Y=(double *)malloc(xy_num*sizeof(double));
	for (t=0;t<xy_num;t++)
	{
		X[t]=(2*u[t]-(uvMaxMin[0]+uvMaxMin[2]))/(uvMaxMin[0]-uvMaxMin[2]);
		Y[t]=(2*v[t]-(uvMaxMin[1]+uvMaxMin[3]))/(uvMaxMin[1]-uvMaxMin[3]);
	}
	free(u);
	free(v);
	int k,l;
	int n=10;
	double *sum_x,*sum_y,*sum;
	sum=(double *)malloc(n*n*sizeof(double));
	sum_x=(double *)malloc(n*sizeof(double));
	sum_y=(double *)malloc(n*sizeof(double));
	for (k=0;k<n;k++)
	{
		sum_x[k]=0.0;
		sum_y[k]=0.0;
		for (l=0;l<=k;l++)
		{
			*(sum+k*n+l)=0.0;
		}
	}
	double *tep;
	tep=(double *)malloc(n*sizeof(double));

	for (t=0;t<xy_num;t++)
	{
		q(X[t],Y[t],tep);
	    for (k=0;k<n;k++)
		{
			sum_x[k]=sum_x[k]+tep[k]*x[t%x_num];
            sum_y[k]=sum_y[k]+tep[k]*y[t/x_num];
		    for (l=0;l<=k;l++)
                *(sum+k*n+l)+=tep[k]*tep[l];
		}
	}

	for (k=0;k<n-1;k++)
	{
		for (l=k+1;l<n;l++)
		{
			*(sum+k*n+l)=*(sum+l*n+k);
		}
	}
	
	free(tep);
	free(x);
	free(y);
	free(X);
	free(Y);
    
	int r;
	double temp;
/*	//the first row of U is equal to A
	//the first column of matrix L */
	for (i=1;i<n;i++)
	    *(sum+i*n+0)=*(sum+i*n+0)/(*sum);
/*    //A=LU*/
	for (r=1;r<n;r++)
	{  
	/*	//the r-th row of matrix U*/
		for (i=r;i<n;i++)
		{
		    temp=0.0;
		    for (k=0;k<r;k++)
			{
			    temp+=(*(sum+r*n+k))*(*(sum+k*n+i));
			}
		    *(sum+r*n+i)=*(sum+r*n+i)-temp;
		}
/*        //the r-th column of matrix L*/
        for (i=r+1;i<n;i++)
		{
		    temp=0.0;
		    for (k=0;k<r;k++)
			{
				temp+=(*(sum+i*n+k))*(*(sum+k*n+r));
			}
			*(sum+i*n+r)=(*(sum+i*n+r)-temp)/(*(sum+r*n+r));
		}
	}

/*	//////////Calculate neCoefArray_u///////////
	//solve Ly=b*/
	for (i=1;i<n;i++)
	{
		temp=0.0;
		for (k=0;k<i;k++)
			temp+=(*(sum+i*n+k))*(*(sum_x+k));
		sum_x[i]-=temp;
	}
	/*//solve Ux=y*/
	double * neCoefmid_u;
	neCoefmid_u=(double *)malloc(n*sizeof(double));
	neCoefmid_u[n-1]=sum_x[n-1]/(*(sum+(n-1)*n+n-1));
	for (i=n-2;i>=0;i--)
	{
		temp=0.0;
		for (k=i+1;k<n;k++)
		{
			temp+=(*(sum+i*n+k))*neCoefmid_u[k];
		}
		neCoefmid_u[i]=(sum_x[i]-temp)/(*(sum+i*n+i));
	}
	free(sum_x);
    * neCoefArray_u=neCoefmid_u[0]-0.5*(neCoefmid_u[3]+neCoefmid_u[5]);
	* (neCoefArray_u+1)=neCoefmid_u[1]-0.5*(neCoefmid_u[8]+3*neCoefmid_u[6]);
	* (neCoefArray_u+2)=neCoefmid_u[3]*1.5;
	* (neCoefArray_u+3)=neCoefmid_u[6]*2.5;
	* (neCoefArray_u+4)=neCoefmid_u[2]-0.5*(neCoefmid_u[7]+3*neCoefmid_u[9]);
	* (neCoefArray_u+5)=neCoefmid_u[5]*1.5;
	* (neCoefArray_u+6)=2.5*neCoefmid_u[9];
	* (neCoefArray_u+7)=neCoefmid_u[4];
    * (neCoefArray_u+8)=1.5*neCoefmid_u[7];
	* (neCoefArray_u+9)=1.5*neCoefmid_u[8];
	free(neCoefmid_u);
    
	/*/////////////Calculate neCoefArray_v///////////
	//solve Ly=b*/
	for (i=1;i<n;i++)
	{
		temp=0.0;
		for (k=0;k<i;k++)
			temp+=(*(sum+i*n+k))*(*(sum_y+k));
		sum_y[i]-=temp;
	}
	/*//solve Ux=y*/
	double * neCoefmid_v;
	neCoefmid_v=(double *)malloc(n*sizeof(double));
	neCoefmid_v[n-1]=sum_y[n-1]/(*(sum+(n-1)*n+n-1));
	for (i=n-2;i>=0;i--)
	{
		temp=0.0;
		for (k=i+1;k<n;k++)
		{
			temp+=(*(sum+i*n+k))*neCoefmid_v[k];
		}
		neCoefmid_v[i]=(sum_y[i]-temp)/(*(sum+i*n+i));
	}
	free(sum_y);
	free(sum);
	* neCoefArray_v=neCoefmid_v[0]-0.5*(neCoefmid_v[3]+neCoefmid_v[5]);
	* (neCoefArray_v+1)=neCoefmid_v[1]-0.5*(neCoefmid_v[8]+3*neCoefmid_v[6]);
	* (neCoefArray_v+2)=neCoefmid_v[3]*1.5;
	* (neCoefArray_v+3)=neCoefmid_v[6]*2.5;
	* (neCoefArray_v+4)=neCoefmid_v[2]-0.5*(neCoefmid_v[7]+3*neCoefmid_v[9]);
	* (neCoefArray_v+5)=neCoefmid_v[5]*1.5;
	* (neCoefArray_v+6)=2.5*neCoefmid_v[9];
	* (neCoefArray_v+7)=neCoefmid_v[4];
	* (neCoefArray_v+8)=1.5*neCoefmid_v[7];
	* (neCoefArray_v+9)=1.5*neCoefmid_v[8];
	free(neCoefmid_v);
}

double Calculate_x(double u,double v,double * neCoefArray_u,double *uvMaxMin)
{
	double X,Y;
	X=(2*u-(uvMaxMin[0]+uvMaxMin[2]))/(uvMaxMin[0]-uvMaxMin[2]);
	Y=(2*v-(uvMaxMin[1]+uvMaxMin[3]))/(uvMaxMin[1]-uvMaxMin[3]);
	double a;
	a=neCoefArray_u[0]+X*(neCoefArray_u[1]+X*(neCoefArray_u[2]+X*neCoefArray_u[3]))+Y*(neCoefArray_u[4]+Y*(neCoefArray_u[5]+Y*neCoefArray_u[6]))+X*Y*(neCoefArray_u[7]+neCoefArray_u[8]*X+neCoefArray_u[9]*Y);
 	return a;
}

double Calculate_y(double u,double v,double * neCoefArray_v,double *uvMaxMin)
{
	double X,Y;
	X=(2*u-(uvMaxMin[0]+uvMaxMin[2]))/(uvMaxMin[0]-uvMaxMin[2]);
	Y=(2*v-(uvMaxMin[1]+uvMaxMin[3]))/(uvMaxMin[1]-uvMaxMin[3]);
	double a;
    a=neCoefArray_v[0]+X*(neCoefArray_v[1]+X*(neCoefArray_v[2]+X*neCoefArray_v[3]))+Y*(neCoefArray_v[4]+Y*(neCoefArray_v[5]+Y*neCoefArray_v[6]))+X*Y*(neCoefArray_v[7]+neCoefArray_v[8]*X+neCoefArray_v[9]*Y);
 	return a;
}
double Calculate_u(double x,double y,double * poCoefArray_x)
{
	double a;
	a=x*(x*(x*poCoefArray_x[6]+poCoefArray_x[3])+poCoefArray_x[1])+poCoefArray_x[0];
	a+=y*(y*(y*poCoefArray_x[9]+poCoefArray_x[5])+poCoefArray_x[2]);
	a+=x*y*(poCoefArray_x[4]+poCoefArray_x[7]*x+poCoefArray_x[8]*y);

	//a=poCoefArray_x[0]+poCoefArray_x[1]*x+poCoefArray_x[2]*y+poCoefArray_x[3]*x*x+poCoefArray_x[4]*x*y+poCoefArray_x[5]*y*y+poCoefArray_x[6]*x*x*x+poCoefArray_x[7]*x*x*y+poCoefArray_x[8]*x*y*y+poCoefArray_x[9]*y*y*y;
	return a;
}

double Calculate_v(double x,double y,double * poCoefArray_y)
{
	double a;
	a=x*(x*(x*poCoefArray_y[6]+poCoefArray_y[3])+poCoefArray_y[1])+poCoefArray_y[0];
	a+=y*(y*(y*poCoefArray_y[9]+poCoefArray_y[5])+poCoefArray_y[2]);
	a+=x*y*(poCoefArray_y[4]+poCoefArray_y[7]*x+poCoefArray_y[8]*y);

//	a=poCoefArray_y[0]+poCoefArray_y[1]*x+poCoefArray_y[2]*y+poCoefArray_y[3]*x*x+poCoefArray_y[4]*x*y+poCoefArray_y[5]*y*y+poCoefArray_y[6]*x*x*x+poCoefArray_y[7]*x*x*y+poCoefArray_y[8]*x*y*y+poCoefArray_y[9]*y*y*y;
	return a;
}
/*
void main()
{
	int totalX,totalY;//x,y方向总象素的个数
	int coarseX,coarseY;//x,y方向的分辨率（取样间隔）
	int x_num,y_num;//x,y方向的取样点数
	//分别赋值
	totalX=totalY=30000;
	coarseX=coarseY=100;
	x_num=int(totalX/coarseX)+1;
	y_num=int(totalY/coarseY)+1;
	
	double *poCoefArray_x,*poCoefArray_y;
	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;

	double *neCoefArray_u,*neCoefArray_v;
	neCoefArray_u=(double *)malloc(10*sizeof(double));
	neCoefArray_v=(double *)malloc(10*sizeof(double));

	//最小二乘法求多项式的逆
	double *uvMaxMin;
	uvMaxMin=(double *)malloc(4*sizeof(double));
	//计算起始时间
    struct _timeb timebuffer1,timebuffer2;
    _ftime( &timebuffer1);
    int t1=timebuffer1.time;

	leastMult(totalX,totalY,coarseX,coarseY,poCoefArray_x,poCoefArray_y,neCoefArray_u,neCoefArray_v,uvMaxMin);

	//计算结束时间
	_ftime( &timebuffer2);
    int t2=timebuffer2.time;
    double usedTime=(t2-t1)+(timebuffer2.millitm-timebuffer1.millitm)*1e-6; 
    printf("The used time is %f\n",usedTime);
	//输出逆多项式的系数
	int i;
	printf("         +++++++Output Coefficient++++++\n");
	for (i=0;i<10;i++)
		printf("neCoefArray_u[%d]:  %e;\tneCoefArray_v[%d]:  %e\n",i,neCoefArray_u[i],i,neCoefArray_v[i]);
	
	//Calculate the error
	int *x,*y;
	x=(int *)malloc(x_num*sizeof(int));
	y=(int *)malloc(y_num*sizeof(int));
	long xy_num=x_num*y_num;

	int j;
    //给x,y赋值
	for (i=0;i<x_num;i++)
		x[i]=i*coarseX;
	for (j=0;j<y_num;j++)
		y[j]=j*coarseY;

	double *u,*v;
	u=(double *)malloc(xy_num*sizeof(double));
	v=(double *)malloc(xy_num*sizeof(double));
	double x_err,y_err,errorAverage;
    double errorSum=0;
	for (j=0;j<y_num;j++)
	{
		for (i=0;i<x_num;i++)
		{
			u[j*x_num+i]=Calculate_u(x[i],y[j],poCoefArray_x);
			v[j*x_num+i]=Calculate_u(x[i],y[j],poCoefArray_y);
			x_err=Calculate_x(u[j*x_num+i],v[j*x_num+i],neCoefArray_u,uvMaxMin)-x[i];
			y_err=Calculate_y(u[j*x_num+i],v[j*x_num+i],neCoefArray_v,uvMaxMin)-y[j];
			errorSum+=sqrt(x_err*x_err+y_err*y_err);
		}
	}
	free(poCoefArray_x);
    free(poCoefArray_y);
	free(neCoefArray_u);
	free(neCoefArray_v);
	free(uvMaxMin);
	free(x);
	free(y);
	free(u);
	free(v);
	errorAverage=errorSum/xy_num;
    printf("             +++++++The Error+++++++\n");
    printf("           errorAverage=%e\n\n",errorAverage);
}*/