#define MAXORDER  30  //max order of	TRRC
#define pi  3.141593		
#define E 0.01		  //max deficiency of zero	
#include <math.h>
#include <stdio.h>
#include <string.h>
#define INFILEPATH1 "E:\\peak\\bmp\\"//read path
#define INFILEPATH2 "E:\\peak\\bmp\\"//read path
#define OUTFILEPATH "E:\\peak\\rbuildbmp\\cg1.bmp"//write path
#define BMP ".bmp"
#define TXT ".txt"
//#define  maxsize    64




//////////////////MAIN//////////////////////////////////////////////
//int   readtxtTO_RC(int RC[][64]);
double CUR(double y[3]);
void   maket(double p[2],int RC[][64],int m,int n,int maxsize);
void   tran2(double  b[MAXORDER][MAXORDER][2],int RC[][64],int maxm,int maxn,int maxsize);
void printb(double b[MAXORDER][MAXORDER][2],int maxm,int maxn);
void RbuildMAP(int RC[][64],double b[][MAXORDER][2],int maxm,int maxn,int maxsize);
void stander(int RC[][64]);
int WRITEbmp(int RC[][64]);
int ReadPicPe(int Rc[64][64]);


void main()
{
int    RC[64][64];
double  b[MAXORDER][MAXORDER][2];
int maxm,a,i,j,t=0;
int maxn;
int maxsize=64;

        printf("输入maxm,maxn：");
		scanf("%d,%d",&maxm,&maxn);
		ReadPicPe(RC); 
		
		tran2(b,RC,maxm,maxn,maxsize);
		printb(b,maxm, maxn);
        RbuildMAP(RC,b,maxm,maxn,maxsize);
	    stander(RC);
		WRITEbmp(RC);

}




///////////////////////////end MAIN///////////////////////////




//1//////求阶层的函数///////////
double JC(int c)
{int i;
double w;
w=1.0;
for(i=1;i<=c;i++)
{w=w*i;}
return(w);
}
/////////end JC//////////




//2////////求径向函数中Bns的值///////////
               
double Bns(int a,int b,int p)
{double g;
double c,d,e,f,k;
int h,i,j;
double s;


h=a-p;
c=JC(h);
i=(a+fabs(b))/2-p;
d=JC(i);
e=JC(p);
j=(a-fabs(b))/2-p;
f=JC(j);
k=pow(-1,p)*c;
s=d*e*f;
g=double(k)/double(s);

return(g);
}
///////////end Bns///////////






//3////////////RED函数  ////////////////////////


double CUR(double y[3])
{
	int n,s,m,A;
	double x;   //x相当于r//
	double  B,C,D;	
	double F,G,H;
	
	n=int(y[0]);
	x=y[1];
    m=int(y[2]);
	
	
	F=0.0;
if(fabs(x)<E)
{x=E;}
	A=(n-fabs(m))/2;
	
	for(s=0;s<=A;s++)
	{
		C=Bns(n,m,s);
		D=C*pow(x,n-2*s);
	 //  printf("%f   ",D);
		F=F+D;
  //printf("%f\n",F);
	}
 
    return(F);
}
/////////////end///////////////////////////////	









//4//////////////计算图像m,n阶moment////////////////////

void   maket(double p[2],int RC[][64],int m,int n,int maxsize)

{
	 
	
	int maxtheta=int(maxsize*pi);
	int numr=maxsize/2;
	int i,r;
	int cx=maxsize/2,cy=maxsize/2;
	double theta,cthe,sthe;
	double x[3],s;
         x[0]=n;
		 x[2]=m;
		 p[0]=0;
		 p[1]=0;
		 for(i=0;i<maxtheta;i++)
		 {
			 theta=2*i/double(maxsize);
			 sthe=sin(theta);   cthe=cos(theta);
			 s=0;
			 for(r=0;r<numr;r++)
			 {
				 x[1]=r/double(numr);
				 s+=RC[cx+int(r*cthe)][cy+int(r*sthe)]*CUR(x)*x[1]*2*pi;
			 }
		 
          
		 p[0]+=s*cos(m*theta)/double(maxtheta)/numr;
         p[1]+=-s*sin(m*theta)/double(maxtheta)/numr;
		 
		 }

          p[0]=p[0]*(x[0]+1)/pi;
          p[1]=p[1]*(x[0]+1)/pi;
}
////////////////////////////////end3////////////////////*/





//5//////////把moment存到数组b[][][]/////////////////////

void tran2(double  b[MAXORDER][MAXORDER][2],int RC[][64],int maxm,int maxn,int maxsize)
		                  
{
	int i,j;
	double p[2];																		


	              for(i=0;i<maxm;i++)
					  for(j=i;j<maxn;j++)
					  {
						  if(  (j-i)%2==0 )
						  { 
						  maket(p,RC,i,j,maxsize);
					      b[i][j][0]=p[0];
                          b[i][j][1]=p[1];
						  }
					}
}
//////////////////////end4/////////////////////////////////*/




//6///////////////////输出数组B////////////////////////////
void printb(double b[MAXORDER][MAXORDER][2],int maxm,int maxn)
{
	
	int i,j;
	for(i=0;i<maxm;i++)
	{
		for(j=i;j<maxn;j++)
		{   if(  (j-i)%2==0 )
		    	{
			     printf("%-10.5f   ",b[i][j][0]);
                 printf("%-10.5f   ",b[i][j][1]);
				 printf("\n");}

		}
        
	}
}

//7///////////////////end输出数组B///////////////////////////*/






//8////////////////////图象重建//////////////////////////////


void RbuildMAP(int RC[][64],double b[][MAXORDER][2],int maxm,int maxn,int maxsize)
			   

{
	int cx,cy;
	int i,j,m,n;
	double theta,s,p;
	double x[3];

	         
	          cx=maxsize/2;    cy=maxsize/2;


			  for(i=0;i<pi*maxsize;i++)
			  {
				  theta=i*2/double(maxsize);

				  for(j=0;j<maxsize/2;j++)
				  {
					  p=0;
					  for(n=0;n<maxn;n++)
					  {
						  s=0;
						  for(m=1;m<maxm;m++)
						  {
							if(n>=m && (n-m)%2==0 )
							{s+=2*(b[m][n][0]*cos(m*theta)-b[m][n][1]*sin(m*theta));
						     x[1]=2*j/double(maxsize);  
							 x[0]=n;
						     x[2]=m;
						     p+=CUR(x)*s;
							}
						 }

							if(n>=m && (n-m)%2==0 )
						  
							{ x[1]=2*j/double(maxsize);  
						      x[0]=n;
						      x[2]=0;
						      p+=CUR(x)*b[0][n][0];
							}
					  
					  
					  }

					  RC[cx+int(j*cos(theta))][cy+int(j*sin(theta))]=int(p);
				  }
			  }
}
////////////////////////////////end5////////////////////////////////////////*/










//9/////////////////归一化数据//////////////////////////////////////////
void stander(int RC[][64])
{
	int i,j;
	float max=0;

	for(i=0;i<64;i++)
		for(j=0;j<64;j++)
			if(RC[i][j]>max)  max=RC[i][j]; 


	for(i=0;i<64;i++)
		for(j=0;j<64;j++)
			if(RC[i][j]<0)  RC[i][j]=0;
			else RC[i][j]=RC[i][j]/max*255;
}



////////////////////////归一化数据end//////////////////////////////////////////*/




//10/////////////////write to rbuild-bmp////////////////////////////////
int WRITEbmp(int RRC[][64])
{
	int i,j;
	char ch, anyone[100],header[1078];
	FILE *fp;
		
    strcpy(anyone,OUTFILEPATH);
    
	fp=fopen(anyone,"rb");
	
	if(fp==NULL)//打开失败
	{
		printf("cannot open file\n");
		fclose(fp);
        return(0);
	}
    
    for(i=0;i<1078;i++)
	{fscanf(fp,"%c",&header[i]);}
    fp=fopen(anyone,"wb");
	if(fp==NULL)//打开失败
	{
	printf("cannot open file\n");
	fclose(fp);
	return(0);	
	}
   
	for(i=0;i<1078;i++)
	{fprintf(fp,"%c",header[i]);}
    
    for(i=63;i>=0;i--)
		for(j=0;j<64;j++)
		{
			ch=RRC[i][j];
			fprintf(fp,"%c",ch);
		}
		fclose(fp);
		return(1);
    
}
/////////////////////write to rbuild-bmp end///////////////////////////////*/




//11////////////////////////read picture pels ///////////////////
int ReadPicPe(int Rc[][64])
{FILE *fp;
 char FileName[100],ch;
 int i,j,t=0;
 long size;
 strcpy(FileName,INFILEPATH1);
 strcat(FileName,"hde.bmp");
 fp=fopen(FileName,"rb");
 if(fp==NULL)
 {	printf("cannot open file\n");
    return(0);
 
 }
 else
    printf("open file\n");

 //fseek(fp,10L,0);
 //fscanf(fp,"%l",&size);
 //fseek(fp,size,0);
 fseek(fp,1078L,1);
 for(i=63;i>=0;i--)
  for(j=0;j<64;j++)
  {    
	  fscanf(fp,"%c",&ch);
      Rc[i][j]=ch;
	  if(Rc[i][j]==-1)
		  Rc[i][j]=255;
  }
/*///print rc data////////////////// 
for(i=0;i<64;i++)
for(j=0;j<64;j++)
{    
	printf("%d ",Rc[i][j]);
    t++;//数据个数
}*/

  printf("\n");
  //printf("T=%d\n",t);
 return(1);	  
}