#include <stdio.h>
#include <math.h>
#include <string.h>
#include "nomouse.h"

#define MAXCOL   11
#define MAXSIZE 6553

int ScreenXs,ScreenYs;
double Mean[MAXCOL],Sdv[MAXCOL],TotVar,Var;
double XtX[MAXCOL][MAXCOL];
double Eignv[MAXCOL-1][MAXCOL-1];
unsigned char Xrow[MAXCOL][MAXSIZE/2];

double Nu[MAXCOL],Mu[MAXCOL],Lambda[MAXCOL];
double MaxVal[MAXCOL],MinVal[MAXCOL],Egn[MAXCOL][MAXCOL];
double RangeVal[MAXCOL];

unsigned char ZeroBuff[MAXSIZE],WhiteBuff[MAXSIZE];

int Cr=0xd;
long load_xtx_1();
int comp_pc();
void plot_eigen_val();
void find_max_min();

struct image_header
{
	int row,col;
	int sample;
	long header_bytes;
	int numimage;
}Hdr;
long Sample_Rate=1;
struct images
{
	char name[100];
	FILE *fp;
	int data_max,data_min;
}Img[MAXCOL-1];

main(argc,argv)

int argc;
char *argv[];

{
	int i,j,k;
	char string[100],lblname[100];
	int numimg,ret;
	char outname[100];

	if(argc<2)
	{
		printf(
"     You must give a file name without a type and there must be\n");
		printf(
"file of type '.lbl' of that name.  The minimum information needed\n");
		printf(
"in the label file is as follows:\n\n");
		printf("     FILE_TYPE          = IMAGE\n");
		printf("     IMAGE_LINES        =  400\n");
		printf("     LINE_SAMPLES       =  640\n");
		printf("     SET_POINTER        =    6\n");
		printf("       canonb1.dat  0  255\n");
		printf("       canonb2.dat  0  255\n");
		printf("       canonb3.dat  0  255\n");
		printf("       canonb4.dat  0  255\n");
		printf("       canonb5.dat  0  255\n");
		printf("       canonb7.dat  0  255\n");
		printf("     END\n\n");
		printf(
"     The 0 and 255 mark the maximum and minimum values which will\n");
		printf(
"be considered as data.  Values outside this range will be\n");
		printf(
"considered to be missing.  If any value is missing the entire\n");
		printf(
"point is missing.  You may have three values for 3 point scaling,\n");
		printf(
"in which case the first and last will be assumed to be the data\n");
		printf("range.\n");
		exit(0);
	}
	strcpy(lblname,argv[1]);
	if((ret=read_lbl(lblname))<0)
	{
		printf("Could not process label file '%s'\n",lblname);
		printf("error = %d\n",ret);
		print_error(ret);
		exit(0);
	}
	printf("How many 'Principal Component Score Images' do you want?\n");
	printf("Maximum = %d\n",Hdr.numimage);
	scanf("%d",&numimg);
	if(numimg<1)
		numimg=1;
	if(numimg>Hdr.numimage)
		numimg=Hdr.numimage;
	printf("Give file name for output set  -- 6 or less characters.\n");
	scanf("%s",outname);
	load_xtx_1();
	print_xtx("%10.0lf");
	printf("Do you want the means removed?  (y or n)\n");
	if(getch()=='y')
	{
		remove_mean();
		print_xtx("%10.3lf");
		printf("Do you want the variance removed?  (y or n)\n");
		if(getch()=='y')
		{
			remove_variance();
			print_xtx("%10.6lf");
		}
	}
	load_Egn();
	comp_pc(Hdr.numimage,numimg,outname);

}

/********************************************************************
**
**	computes principle components
**
********************************************************************* */

int comp_pc(numvar,numimg,outname)

int numvar,numimg;
char *outname;

{
	int i,j,k,ev;
	double num,numpts;
	int numit;
	int size=numvar+1;
	FILE *fp;
	int numev=numvar;
	int ret;

	printf("\nX'X Array\n\n");
	print_Egn("%10.4lg");
	TotVar=0.0;
	for(i=0;i<size-1;i++)
		for(j=0;j<size-1;j++)
			TotVar+=Egn[i][j];
	printf("\nTotal variance = %lf\n\n   lambda\n",TotVar);
	for(ev=0;ev<numev;ev++)
	{
		numit=comp_eigen(Egn,size,ev);
		for(i=0;i<size-1;i++)
			for(j=0;j<size-1;j++)
				Egn[i][j]-=Lambda[ev]*Nu[i]*Nu[j];
		printf("%10.5lf ev #%d ",Lambda[ev],ev+1);
		for(i=0;i<size-1;i++)
		{
			Eignv[ev][i]=Nu[i];
			printf("%8.3lf",Eignv[ev][i]);
		}
		printf("\n");
	}
	printf("\n");
	scale_images();
	make_lbl(numimg,outname);
	for(i=0;i<numimg;i++)
		if((ret=make_pcscore_img(i,outname))<0)
		{
			printf("\nError %d in making image\n\n",ret);
			exit(0);
		}
	return(1);
}

/********************************************************************
**
**	sets principle component max and min for making PC score images
**
**       -1 = hit end of file reading data
**
********************************************************************* */

int make_pcscore_img(img,name)

int img;
char *name;

{
	int i,j,k,m,n;
	long offset=Hdr.header_bytes;
	int row=Hdr.row,col=Hdr.col,valid,check=0;
	long lrow=row,lcol=col,collast,datasize=lrow*lcol,dataread=0;
	long maxread=MAXSIZE/2;
	double dval;
	FILE *fpout;
	char string[100];
	int val;

	sprintf(string,"%s%02d.egn",name,img+1);
	fpout=fopen(string,"wb");
	if(!fpout)
	{
		printf("Could not open '%s' to write.\n\n",string);
		exit(0);
	}
	lrow=(datasize+maxread-1)/maxread;
	lcol=maxread;
	for(k=0;k<Hdr.numimage;k++)
		rewind(Img[k].fp);
	if(offset>0)
		for(k=0;k<Hdr.numimage;k++)
			fseek(Img[k].fp,offset,SEEK_SET);
	for(k=0;k<Hdr.numimage;k++)
	{
		if(Img[k].data_min>0||Img[k].data_max<255)
			check=1;
	}
	printf("%cMakeing PC Score image '%s':\n       of %ld reads",Cr,string,lrow);
	for(i=0;i<lrow;i++)
	{
		printf("%c%d",Cr,i+1);
		if((datasize-dataread)<lcol)
			lcol=datasize-dataread;
		for(k=0;k<Hdr.numimage;k++)
			if(fread((char *)Xrow[k],sizeof(char),(int)lcol,Img[k].fp)<(int)lcol)
				return(-1);
		dataread+=lcol;
		for(j=0;j<(int)lcol;j++)
		{
			valid=1;
			if(check==1)
			{
				for(k=0;k<Hdr.numimage;k++)
				{
					if(Xrow[k][j]<Img[k].data_min||Xrow[k][j]>Img[k].data_max)
					{
						valid=0;
						break;
					}
				}
			}
			val=0;
			if(valid==1)
			{
				dval=0.0;
				for(n=0;n<Hdr.numimage;n++)
					dval+=(double)Xrow[n][j]*Eignv[img][n];
				val=(dval-MinVal[img])/RangeVal[img]*255.0+1.0;
				if(val<1)
					val=1;
				if(val>255)
					val=255;
			}
			Xrow[MAXCOL-1][j]=val;
		}
		if(fwrite((char *)Xrow[MAXCOL-1],sizeof(char),(int)lcol,fpout)<(int)lcol)
			return(-2);
	}
	fclose(fpout);
	return(1);
}






/***********************************************************************
**
**
**
**
************************************************************************/

int make_lbl(num,name)

int num;
char *name;

{
	int i,j,k;
	FILE *fpout;
	char string[100];

					/* make group label */
	for(k=0;k<strlen(name);k++)
		if(name[k]=='.')
			name[k]='\0';
	name[6]='\0';
	sprintf(string,"%s_e.lbl",name);
	fpout=fopen(string,"rt");
	if(fpout)
	{
		printf(
		"File '%s' already exists.  Do you want to overwrite it?  (y or n)\n\n",
			string);
		if(getch()!='y')
			exit(0);
		fclose(fpout);
	}
	fpout=fopen(string,"wt");
	if(!fpout)
	{
		printf("Could not open '%s' to write.\n\n",string);
		exit(0);
	}
	fprintf(fpout,"FILE_TYPE          = IMAGE\n");
	fprintf(fpout,"IMAGE_LINES        =  %5d\n",Hdr.row);
	fprintf(fpout,"LINE_SAMPLES       =  %5d\n",Hdr.col);
	fprintf(fpout,"SET_POINTER        =  %5d\n",num);
	for(k=0;k<num;k++)
		fprintf(fpout,"  %s%02d.egn   1  255\n",name,k+1);
	fprintf(fpout,"END\n");
	fclose(fpout);
	for(k=0;k<num;k++)
	{
		sprintf(string,"%s%02d.lbl",name,k+1);
		fpout=fopen(string,"wt");
		if(!fpout)
		{
			printf("Could not open '%s' to write.\n\n",string);
			exit(0);
		}
		fprintf(fpout,"FILE_TYPE          = IMAGE\n");
		fprintf(fpout,"IMAGE_LINES        =  %5d\n",Hdr.row);
		fprintf(fpout,"LINE_SAMPLES       =  %5d\n",Hdr.col);
		fprintf(fpout,"IMAGE_POINTER      = '%s%02d.egn'\n",name,k+1);
		fprintf(fpout,"IMAGE_NOTE:  Principal component score image\n");
		fprintf(fpout,"     input images:\n");
		for(i=0;i<Hdr.numimage;i++)
			fprintf(fpout,"     %2d -- '%s'\n",i+1,Img[i].name);
		fprintf(fpout,"     Lambda = %10.5lf\n",Lambda[k]);
		fprintf(fpout,"     Eigen Vector:\n");
		for(i=0;i<Hdr.numimage;i++)
			fprintf(fpout,"%9.5lf ",Eignv[k][i]);
		fprintf(fpout,"\n");
		fprintf(fpout,"     Principal Component Minimum = %12.5lf\n",MinVal[k]);
		fprintf(fpout,"     Principal Component Maximum = %12.5lf\n",MaxVal[k]);
		fprintf(fpout,"END_IMAGE_NOTE\n");
		fprintf(fpout,"END\n");
		fclose(fpout);
	}
}





/********************************************************************
**
**	sets principle component max and min for making PC score images
**
**       -1 = hit end of file reading data
**
********************************************************************* */

int scale_images()

{
	int i,j,k,m,n;
	long offset=Hdr.header_bytes;
	int row=Hdr.row,col=Hdr.col,valid,check=0;
	long lrow=row,lcol=col,collast,datasize=lrow*lcol,dataread=0;
	long maxread=MAXSIZE/2;
	double dval;

	lrow=(datasize+maxread-1)/maxread;
	lcol=maxread;
	for(k=0;k<Hdr.numimage;k++)
		rewind(Img[k].fp);
	if(offset>0)
		for(k=0;k<Hdr.numimage;k++)
			fseek(Img[k].fp,offset,SEEK_SET);
	for(k=0;k<Hdr.numimage;k++)
	{
		if(Img[k].data_min>0||Img[k].data_max<255)
			check=1;
		MaxVal[k]=-10000000000.0;
		MinVal[k]= 10000000000.0;
	}
	printf("Scaling Image principal components:\n       of %ld reads",lrow);
	for(i=0;i<lrow;i++)
	{
		printf("%c%d",Cr,i+1);
		if((datasize-dataread)<lcol)
			lcol=datasize-dataread;
		for(k=0;k<Hdr.numimage;k++)
			if(fread((char *)Xrow[k],sizeof(char),(int)lcol,Img[k].fp)<lcol)
				return(-1);
		dataread+=lcol;
		for(j=0;j<lcol;j++)
		{
			valid=1;
			if(check==1)
			{
				for(k=0;k<Hdr.numimage;k++)
				{
					if(Xrow[k][j]<Img[k].data_min||Xrow[k][j]>Img[k].data_max)
					{
						valid=0;
						break;
					}
				}
			}
			if(valid==1)
			{
				for(m=0;m<Hdr.numimage;m++)
				{
					dval=0.0;
					for(n=0;n<Hdr.numimage;n++)
						dval+=(double)Xrow[n][j]*Eignv[m][n];
					if(dval>MaxVal[m])
						MaxVal[m]=dval;
					if(dval<MinVal[m])
						MinVal[m]=dval;
				}
			}
		}
	}
	printf("\nEV    minimum   maximum    range\n");
	for(m=0;m<Hdr.numimage;m++)
	{
		RangeVal[m]=MaxVal[m]-MinVal[m];
		printf("%2d  %10.5lf%10.5lf%10.5lf\n",
			m+1,MinVal[m],MaxVal[m],RangeVal[m]);
	}
	return(1);
}




/***********************************************************************
**
**
**
**
************************************************************************/

int load_Egn()

{
	int i,j,k;

	for(i=0;i<Hdr.numimage;i++)
		for(j=0;j<Hdr.numimage;j++)
			Egn[i][j]=XtX[i+1][j+1];

}

/***********************************************************************
**
**
**
**
************************************************************************/

int print_Egn(format)

char *format;

{
	int i,j,k;

	for(i=0;i<Hdr.numimage;i++)
	{
		for(j=0;j<Hdr.numimage;j++)
			printf(format,Egn[i][j]);
		printf("\n");
	}
}



/***********************************************************************
**
**
**
**
************************************************************************/

int remove_variance()

{
	int i,j,k;
	double sd[MAXCOL];

	for(i=1;i<=Hdr.numimage;i++)
		sd[i]=sqrt(XtX[i][i]);
	for(i=1;i<=Hdr.numimage;i++)
		for(j=0;j<=Hdr.numimage;j++)
			XtX[i][j]=XtX[i][j]/(sd[i]*sd[j]);
}





/***********************************************************************
**
**
**
**
************************************************************************/

int remove_mean()

{
	int i,j,k;
	double mean[MAXCOL],num=XtX[0][0];

	for(i=0;i<=Hdr.numimage;i++)
		for(j=0;j<=Hdr.numimage;j++)
			XtX[i][j]/=num;
	for(i=0;i<=Hdr.numimage;i++)
		mean[i]=XtX[i][0];
	for(i=0;i<=Hdr.numimage;i++)
		for(j=0;j<=Hdr.numimage;j++)
			XtX[i][j]=XtX[i][j]-mean[i]*mean[j];
}




/***********************************************************************
**
**
**
**
************************************************************************/

int print_xtx(format)

char *format;

{
	int i,j,k;

	for(i=0;i<=Hdr.numimage;i++)
	{
		for(j=0;j<=Hdr.numimage;j++)
			printf(format,XtX[i][j]);
		printf("\n");
	}
}



/********************************************************************
**
**	loads X'X array  for principle components
**
**       -1 = hit end of file reading data
**
********************************************************************* */

long load_xtx_1()

{
	int i,j,k,m,n;
	long offset=Hdr.header_bytes;
	int row=Hdr.row,col=Hdr.col,valid,check=0;
	long lrow=row,lcol=col,collast,datasize=lrow*lcol,dataread=0;
	long maxread=MAXSIZE/2;

	lcol=maxread;
	lrow=(datasize+maxread-1)/maxread;
	for(k=0;k<Hdr.numimage;k++)
		rewind(Img[k].fp);
	if(offset>0)
		for(k=0;k<Hdr.numimage;k++)
			fseek(Img[k].fp,offset,SEEK_SET);
	for(k=0;k<Hdr.numimage;k++)
		if(Img[k].data_min>0||Img[k].data_max<255)
			check=1;
	printf("Loading X'X matrix:\n       of %ld reads",lrow);
	for(j=0;j<lcol;j++)
		Xrow[0][j]=1.0;
	for(i=0;i<lrow;i++)
	{
		printf("%c%d",Cr,i+1);
		if((datasize-dataread)<lcol)
			lcol=datasize-dataread;
		for(k=0;k<Hdr.numimage;k++)
			if(fread((char *)Xrow[k+1],sizeof(char),(int)lcol,Img[k].fp)<lcol)
				return(-1);
		dataread+=lcol;
		for(j=0;j<lcol;j++)
		{
			valid=1;
			if(check==1)
			{
				for(k=0;k<Hdr.numimage;k++)
				{
					if(Xrow[k+1][j]<Img[k].data_min||Xrow[k+1][j]>Img[k].data_max)
					{
						valid=0;
						break;
					}
				}
			}
			if(valid==1)
			{
				for(m=0;m<=Hdr.numimage;m++)
					for(n=0;n<=Hdr.numimage;n++)
						XtX[m][n]+=(double)Xrow[m][j]*(double)Xrow[n][j];
			}
		}
	}
	printf("\n");
	return(dataread);
}



/***********************************************************************
**
**
**
**
************************************************************************/

int print_error(err)

int err;

{
	if(err==-1)
		printf("'.lbl' file not found\n");
	if(err==-2)
		printf("row < 1\n");
	if(err==-3)
		printf("col < 1 or col >= MAXSIZE\n");
	if(err==-4)
		printf("numimage < 2 || numimage >= MAXCOL\n");
	if(err==-5)
		printf("hit end of file while reading image names\n");
	if(err==-6)
		printf("could not open image file\n");
	if(err==-7)
		printf("bad data_max/min values\n");
}




/***********************************************************************
**
**  return:
**
**       -1 = '.lbl' file not found
**       -2 = row < 1
**       -3 = col < 1 or col >= MAXSIZE
**       -4 = numimage < 2 || numimage >= MAXCOL
**       -5 = hit end of file while reading image names
**       -6 = could not open image file
**       -7 = bad data_max/min values
**
************************************************************************/

int read_lbl(name)

char *name;

{
	int i,j,k;
	char string[100],str1[100],str2[100];
	int i1,i2,i3,ret;
	FILE *fp;

	strcpy(string,name);
	for(k=0;k<strlen(string);k++)
		if(string[k]=='.')
			string[k]='\0';
	strcat(string,".lbl");
	fp=fopen(string,"rt");
	if(!fp)
		return(-1);
	Hdr.row=Hdr.col=Hdr.sample=Hdr.numimage=Hdr.header_bytes=0;
	while(fgets(string,99,fp))
	{
		sscanf(string,"%s",str1);
		{
			if(strcmpi(str1,"IMAGE_LINES")==0)
			{
				sscanf(string,"%s%s%d",str1,str2,&Hdr.row);
				if(Hdr.row<1)
					return(-2);
			}
			if(strcmpi(str1,"LINE_SAMPLES")==0)
			{
				sscanf(string,"%s%s%d",str1,str2,&Hdr.col);
				if(Hdr.col<1||Hdr.col>MAXSIZE)
					return(-3);
			}
			if(strcmpi(str1,"HEADER_BYTES")==0)
			{
				sscanf(string,"%s%s%d",str1,str2,&Hdr.header_bytes);
			}
			if(strcmpi(str1,"SET_POINTER")==0)
			{
				sscanf(string,"%s%s%d",str1,str2,&Hdr.numimage);
				if(Hdr.numimage<1||Hdr.numimage>=MAXCOL)
					return(-4);
				for(i=0;i<Hdr.numimage;i++)
				{
					if(!fgets(string,99,fp))
						return(-5);
					i3=-1;
					ret=sscanf(string,"%s%d%d%d",Img[i].name,&i1,&i2,&i3);
					if(ret>2)
					{
						Img[i].data_min=i1;
						if(ret==4)
							Img[i].data_max=i3;
						else
							Img[i].data_max=i2;
					}
					else
					{
						Img[i].data_max=255;
						Img[i].data_min=0;
					}
					if(Img[i].data_max<=Img[i].data_min||Img[i].data_max<1)
						return(-7);
					Img[i].fp=fopen(Img[i].name,"rb");
					if(!Img[i].fp)
						return(-6);
				}
			}
		}
	}	
	fclose(fp);
	return(1);
}


/********************************************************************
**
**	finds max and min eigen values
**
********************************************************************* */

void find_max_min(ev,numvar,buffer,xc,xsize,yc,ysize)

int ev;
int numvar,xc,xsize,yc,ysize;
unsigned char buffer[12][512];

{
	int i,j,k,v,p,m,n,x,y;
	int ix,iy;
	int use;
	double dx;
	double val;

	/*load image buffer */

	MaxVal[ev]=-1.000e050;
	MinVal[ev]=1.000e050;

	for(i=0;i<512;i++)
	{
		plotrow(3,xc,xc+xsize,i,WhiteBuff);
		for(v=0;v<3;v++)
			getrow(v,0,511,i,buffer[v]);
		for(j=0;j<512;j++)
		{
			use=1;
			val=0.0;
			for(p=0;p<3;p++)	/* load xrow */
			{
				if((dx=buffer[p][j])==0.0)
				{
					use=0;
					break;
				}
				val+=dx*Eignv[ev][p];
			}
			if(use==1)
			{
				if(val>MaxVal[ev])MaxVal[ev]=val;
				if(val<MinVal[ev])MinVal[ev]=val;
			}
			else buffer[0][j]=0;
		}
		plotrow(3,xc,xc+xsize,i,ZeroBuff);
	}
}


/********************************************************************
**
** plots eigen values scaled from 1 - 255
**
********************************************************************* */

void plot_eigen_val(ev,numvar,buffer,xc,xsize,yc,ysize)

int ev;
int numvar,xc,xsize,yc,ysize;
unsigned char buffer[12][512];

{
	int i,j,k,v,p,m,n,x,y;
	int ix,iy;
	int use;
	double dx;
	double val;

	/*load image buffer */

	for(i=0;i<512;i++)
	{
		for(v=0;v<numvar;v++)
			getrow(v,0,511,i,buffer[v]);
		for(j=0;j<512;j++)
		{
			val=0.0;
			{
				use=1;
				for(p=0;p<numvar;p++)	/* load xrow */
				{
					if((dx=buffer[p][j])==0.0)
					{
						use=0;
						break;
					}
					val+=dx*Eignv[ev][p];
				}
			}

			val=(val-MinVal[ev])/RangeVal[ev]*254.0+1.5;
			if(val>0.0&&val<256.0&&use==1)
				buffer[ev][j]=val;
			else
				buffer[ev][j]=0;
		}

	/* plot ev */
		plotrow(ev,512,1023,i,buffer[ev]);
	}
}


/********************************************************************
**
**	Mu[]/Lambda ~= nu[]
**
********************************************************************* */

int comp_eigen(egn,size,ev)

double egn[MAXCOL][MAXCOL];
int size,ev;

{
	double error,lam,sumerr,maxerr=size;
	int i,j,k;
	int esize=size-1;
	int numit=0,maxit=50;

	if(esize<2)return;

	maxerr=1.00e-015/sqrt(maxerr);

	/* set Nu[] */

	Nu[0]=1.0;
	for(i=1;i<esize;i++)
		Nu[i]=0.0;

	/* iterate */

	do
	{
		numit+=1;
		lam=0.0;
		for(i=0;i<esize;i++)
		{
			Mu[i]=0;
			for(j=0;j<esize;j++)
			{
				Mu[i]+=Nu[j]*egn[i][j];
			}
			lam+=Mu[i]*Mu[i];
		}
		lam=sqrt(lam);
		for(i=0;i<esize;i++)
			Mu[i]/=lam;
		sumerr=0.0;
		for(i=0;i<esize;i++)
		{
			sumerr+=(Nu[i]-Mu[i])*(Nu[i]-Mu[i]);
			Nu[i]=Mu[i];
		}
	}while(sumerr>maxerr&&numit<maxit);
	Lambda[ev]=lam;
	return(numit);
}



/***********************************************************************
**
**
**
**
************************************************************************/

int set_lut()

{

}