/*============= PROGRAM MACFAC.C ============================================
	Version 1.1: Wed 19 january 1994 14:31
	MacII (with co-processor) version of program FAC for processing images without headers
	
	This program is FREEWARE for the Academic Research Community.
	It MUST NOT BE USED (or COPIED or SOLD) by ANYONE having COMMERCIAL PURPOSES.
	The authors disclaim their responsability for any kind of trouble directly or indirectly 
	induced by the use of this program.
	
	Authors: Noel BONNET and Pierre TREBBIA, Universite de Reims, Faculte des Sciences
										B.P. 347, F51062 Reims cedex FRANCE
	Fax: (33) 26 05 32 50										E-Mail: pierre.trebbia@univ-reims.fr
	
	Questions and enquiries must be send to Pierre TREBBIA (see address & e-mail above)
	
	Those who accept to modify the overall presentation of this program in order to get full benefits
	of a more friendly interface (scrolling menus, radio buttons, windows for looking at images,...)
	are VERY WELCOME: They are encouraged to go that way and to send to Pierre TREBBIA a copy of
	both the modified source code and the compiled application.
	
	The file *.RES is created for keeping a trace of the calculations.
		Open it with any text editor (TeachText for example).
	The file *.LST is created for keeping a trace of all the files created by the program.
		This is also a text file.
	The file *.XCL is a Microsoft Excel (Registred Trade Mark) compatible file. Open it with Excel
	for editing a plot of the locations of the processed images in the factorial space.
==========================================================================
	Adjust the 3 first parameters NBIMA = NBAXE and NBPTS to your specific requirements
	(examples: 20,20,32 or 15,15,64 or 10,10,128 or 6,6,256  etc...)
==========================================================================*/

#define NBPTS 128	/* this parameter may be adjusted to your own image sizes within memory limits (<129)	*/

#define NBIMA 10	
#define NBAXE 10

/*=========================================================================*/

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#include <string.h>

#define buci for(i=0;i<NIMA;i++)
#define buck for(k=0;k<KVVP;k++)
#define buckx for(k=0;k<nombrex;k++)
#define buckr for(k=0;k<kr;k++)
#define bucn for(n=0;n<NBPTS;n++)
#define bucnbl for(nbl=1;nbl<(nblmax+1);nbl++)
#define SEEK_SET 0
#define ext_axe ".AXM"
#define ext_filtre ".RCM"
#define ext_fictif ".FCM"


	char *nomdep,*nomresult1,*nomresult2,*nomresult3,*nomfictif;
	char FILNAM[40],gener[40],nomout[NBIMA][40],inter[5];

	short IX[NBIMA][NBPTS],ITAB[NBPTS];
	short moy[NBIMA],min[NBIMA],max[NBIMA],ia[NBAXE];
	short i,j,k,l,n,nbl,nblmax,nbsig,nbcre,iti,ecranflag,NPL,NLI,NIMA,KVVP;

	long int facnor;
	long int som[NBIMA];

	double PSI[NBIMA][NBAXE],VP[NBIMA],S[NBIMA][NBIMA],vptot,vptotmem,FI[NBPTS][NBAXE];
	double PSIFIC[NBAXE],ftot,MIN,MAX,SJ[NBIMA],SI[NBPTS],w[NBIMA][NBIMA],d[NBIMA];

	short NDIM  = NBIMA;
	short KDIM  = NBAXE;

	FILE *res; FILE *res1; FILE *res2; FILE *res3;
	FILE *nom; FILE *ima; FILE *fictif;
	FILE *fnom[NBIMA];

/*========================================================================*/

void param()

{
	short naxe;

	nomdep=(char *)malloc(65);
	nomresult1=(char *)malloc(65);
	nomresult2=(char *)malloc(65);
	nomresult3=(char *)malloc(65);

	printf("\nMacFAC : FACTORIAL ANALYSIS OF CORRESPONDENCE");
	printf("\n\nBEWARE that ALL the files MUST be in the same folder as this program");
	printf("\n\n\nThere should exist a file named *.LST containing,");
	printf("\nline after line, the names of the N (N<=%d) images to be processed.",NBIMA);
	printf("\n\nThese images should be short integer 16 bits without header");
	printf("\nand must contain an integer multiple of %d pixels.\n",NBPTS);
	printf("\n\n Generic name of the file containing this information [Apple Q => Abort]? ");
	fflush(stdin); scanf("%s",nomdep);
	strcat(nomdep,".LST");
	printf("\n How many pixels per line? ");
	fflush(stdin); scanf("%hd",&NPL);
	printf("\n How many lines per image? ");
	fflush(stdin); scanf("%hd",&NLI);

	NIMA=NBIMA+1;
	while(NIMA>NBIMA)
		{
			printf("\n     How many images to be processed (<=%d) in %s? ",NDIM,nomdep);
			fflush(stdin); scanf( "%hd",&NIMA);
			if(NIMA==0) NIMA=NBIMA+1;
		}

	naxe=NIMA-1; if (naxe>KDIM) naxe=KDIM;
	KVVP=naxe+1;
	while(KVVP>naxe)
		{
			printf("\n     How many factorial axes (<=%d)? ",naxe);
			fflush(stdin); scanf("%hd",&KVVP);
			if(KVVP==0) KVVP=naxe+1;
			if (KVVP>naxe) printf("\n  The number of axes must be smaller than the number of images");
		}	/* end of while */
	nbsig=KVVP;
	printf("\n Generic name of the result file (8 charact. max.)? ");
	fflush(stdin); scanf("%s",nomresult1);
	strcpy(nomresult2,nomresult1); strcpy(nomresult3,nomresult1);
	strcat(nomresult1,".RES"); strcat(nomresult2,".XCL"); strcat(nomresult3,".LST");
	printf("\n Your results will be in the files:\n %s, %s and %s",nomresult1,nomresult2,nomresult3);
	printf("\n Detailed output on the screen (Yes=0; No=1)? ");
	fflush(stdin); scanf("%hd",&ecranflag);
}	/* end of void */

/*==========================================================*/

void AFCINI()

{
	long int temp,tempint;
	short sugg;
	double tempfp;

	res=fopen(nomresult1,"w");  fclose(res); res=fopen(nomresult1,"a");
	res2=fopen(nomresult2,"w"); fclose(res2); res2=fopen(nomresult2,"a");
	res3=fopen(nomresult3,"w"); fclose(res3); res3=fopen(nomresult3,"a");
	fprintf(res3,"%s\n",nomresult1); fprintf(res3,"%s\n",nomresult2);

	fprintf(res,"\nMacFAC: FACTORIAL ANALYSIS OF CORRESPONDENCE");
	fprintf(res,"\n****************************************");
	sugg=128;
	
	tempint=(long) NPL*(long)NLI; tempfp=(double)tempint/(double)NBPTS; temp=tempint/(long)NBPTS;
	if(tempfp!=(double)temp)
			{
					while((tempint%sugg)!=0) sugg--;
					if(sugg==0) printf("\n Not readable images ==>  Create ROIs\n");
					printf("\n  ===> Modify the parameter NBPTS: %5d  I suggest: %5d\n",NBPTS,sugg); exit(0);
			}	/* end of if */
	nblmax=(short) temp;

	nom=fopen(nomdep,"r");
	buci
		{
			fflush(stdin); fscanf(nom,"%s",FILNAM);
			fprintf(res,"\n image # %d : %s",i+1,FILNAM);
			if(ecranflag==0) printf("\n image # %d : %s",i+1,FILNAM);
		}	/* end of buci */
	fclose(res);
	fclose(res2);
	fclose(res3);
	fclose(nom);
	puts("\n\n Coffee break ...");
}	/* end of void */

/*=====================================================================*/

/* Statistics set-up  */

void mima(ind)

short ind;

{
	for(k=0;k<ind;k++)
		{
			som[k]=0; min[k]=32767; max[k]=-32767;
		}	/* end of for */
}	/* end of void */

/*===================================================================*/

void lect_fich(mas)

/*   Reading NBPTS data in each image	*/

short mas;

{

	long int position_seek;

	nom=fopen(nomdep,"r");
	bucn *(SI+n)=0.;
	buci
		{
			fflush(stdin); fscanf(nom,"%s",FILNAM);
			ima=fopen(FILNAM,"rb");
			position_seek=(2*(long)NBPTS*((long)nbl-1));
			fseek(ima,position_seek,SEEK_SET);
			fread(ITAB,NBPTS*sizeof(short),1,ima);
			bucn
				{
					iti=ITAB[n]; IX[i][n]=iti;
					if(mas==0)
						{
							min[i]=(iti<min[i])?iti:min[i];
							max[i]=(iti>max[i])?iti:max[i];
							ftot+=(double) iti;
							*(SJ+i)+=(double) iti;
						}	/* end of if */
					SI[n]+=(double) iti;
				}	/* end of bucn */
			fclose(ima);
		}	/* end of buci */
	fclose(nom);
}	/* end of void */

/*======================================================================*/

void VVPRO3()

/*  Compute the eigenvalues and eigenvectors of a W(N,N) real symmetric matrix.
	W(N,N) is destroyed and replaced by the columns of eigenvectors
	D(N) contains the eigenvalues in decreasing order	*/

{

	double w2,ep,rn,epsil,ww,w22,w2a,wwa,tteta,cteta,steta,ateta,flip;
	short ni,ki,n1,k1,kp;

	buci
		{
			for(l=0;l<NIMA;l++) w[i][l]=S[i][l];
		}	/* end of buci */
	n=NIMA; epsil=0.0000000001; rn=(double) n; w2=0.;
	for(l=0;l<n;l++)
		{
			for(k=0;k<n;k++) w2+=w[l][k]*w[l][k];
		}	/* end of for */
	ep = epsil*w2/rn;
	ni=n*(n-1)/2;
	ki=ni;
	n1=n-1;
	while(ki>0)
		{
			for(k=0;k<n1;k++)
				{
					k1=k+1;
					for(kp=k1;kp<n;kp++)
						{
							w2=0.; ww=0.;
							for(l=0;l<n;l++)
								{
									w2+=w[l][k]*w[l][kp];
									ww+=(w[l][k]+w[l][kp])*(w[l][k]-w[l][kp]);
								}	/* end of for */
							w22=w2*2.;
							w2a=fabs((double)w22);
							if((w2a>=ep)||(ww<0.)) wwa=fabs((double)ww);
							else
								{
									ki-=1;
									if(ki<=0) goto E150;
									else goto E140;
								}	/* end of else */
							if(w2a<=wwa)
								{
									tteta=w2a/wwa;
									cteta=1./sqrt((double)1.+tteta*tteta);
									steta=tteta*cteta;
								}	/* end of if */
							else
								{
									ateta=wwa/w2a;
									steta=1./sqrt((double)1.+ateta*ateta);
									cteta=ateta*steta;
								}	/* end of else */
							cteta=sqrt((double)(1.+cteta)/2.);
							steta/=2.*cteta;
							if(ww<0.)
								{
									flip=cteta;
									cteta=steta;
									steta=flip;
								}	/* end of if */
							if(w22<0.) steta=-steta;
							for(l=0;l<n;l++)
								{
									flip=w[l][k];
									w[l][k]=flip*cteta+w[l][kp]*steta;
									w[l][kp]=-flip*steta+w[l][kp]*cteta;
								}	/* end of for */
							ki=ni;
							E140:;
						}	/* end of for */
				}	/* end of for */
		}	/* end of while */
	E150: for(k=0;k<n;k++)
		{
			d[k]=0.;
			for(l=0;l<n;l++) d[k]+=w[l][k]*w[l][k];
			d[k]=sqrt((double)d[k]);
		}	/* end of for */
	for( k=0;k<n;k++)
		{
			for(l=0;l<n;l++) w[l][k]/=d[k];
		}	/* end of for */
	buci
		{
			VP[i]=(double)d[i];
			for(l=0;l<NIMA;l++) S[i][l]=w[i][l];
		}	/* end of buci */
}	/* end of void */

/*========================================================================*/

void AFCCA()

{
	double sd,conv,conv1;
	short l1,temp;

/*  Set-up */

	ftot=0.;   /* sum of all the pixels in all images */
	AFCINI();
	buci SJ[i]=0.;  /* sum of all the pixels in image i */
	buci
		{
			for(l=0;l<=i;l++) S[i][l]=0.; /* covariance(imagei,imagel) */
		}	/* end of buci */
	mima(NIMA);

/*  Opening the image files */

	bucnbl
		{
			lect_fich(0);

/*		Build the matrix S(I,L) */

			buci
				{
					for(l=0;l<=i;l++)
						{
							bucn
								{
									conv=(double) IX[i][n]; conv1= (double) IX[l][n]; /* conversion floating */
									S[i][l]=(SI[n]!=0.)?S[i][l]+conv*conv1/SI[n]:S[i][l];
								}	/* end of bucn */
						}	/* end of for */
				}	/* end of buci */
		}	/* end of bucnbl */

/*  End of reading */

	res=fopen(nomresult1,"a");
	buci
		{
			conv=SJ[i]/(long)NLI/(long)NPL + 0.5;
			temp=(short) conv;
			if(ecranflag==0)
				printf("\n image # %d: Min=%d, Max=%d, Mean=%d, Tot=%1.9e, Wght=%f", i+1,min[i],max[i],temp,SJ[i],(SJ[i]/ftot));
			fprintf(res,"\n image # %d: Min=%d, Max=%d, Mean=%d, Tot=%1.9e, Wght=%f", i+1,min[i],max[i],temp,SJ[i],(SJ[i]/ftot));

			for(l=0;l<=i;l++)
				{
					sd=SJ[i]*SJ[l]; S[i][l]/=sqrt((double)sd);
					S[l][i]=S[i][l]; /* symetrisation of the covariance matrix */
				}	/* end of for */
		}	/* end of buci */

/*  Working on the matrix */

	vptot=0.;
	buci vptot+=S[i][i]; /* trace of the matrix */
	vptot-=1.;
	fprintf(res,"\n Trace of the covariance matrix before diagonalisation = 1 + %12.7e",vptot);
	printf("\n Trace of the covariance matrix before diagonalisation = 1 + %12.7e",vptot);
	vptotmem=vptot+1.0;
	/*	fprintf(res,"\n Matrix to diagonalise: ");	*/
	if(ecranflag==0) printf("\n Matrix to diagonalise: ");
	buci
		{
			/*	fprintf(res,"\n");	*/
			if(ecranflag==0) printf("\n");
			for(l=0;l<=i;l++)
				{
					/*	fprintf(res,"  %12.7f",S[i][l]);	*/
					if(ecranflag==0) printf("  %12.7f",S[i][l]);
				}	/* end of for */
		}	/* end of buci */

/*  DIAGONALISATION of the matrix  */

	VVPRO3();
	vptot=0.;
	buci vptot+=VP[i];
	vptot-=VP[0];
	fprintf(res,"\n Sum of the eigenvalues after diagonalisation = %f + %15.7e",VP[0],vptot);
	printf("\n Sum of the eigenvalues after diagonalisation = %f + %15.7e",VP[0],vptot);
	vptotmem=vptotmem-VP[0]-vptot;
	printf("\n Difference = %e",vptotmem);
	fprintf(res,"\n Difference = %e",vptotmem);
	fclose(res);
	buci
		{
			if(SJ[i]<1.0e-6)
				{
					printf("\n SJ(%d) RAISED UP TO 1E-6 !!!",i); SJ[i]=1.0e-6;
				}	/* end of if */
			SJ[i]/=ftot; /* normalisation of the sums per image */
		}	/* end of buci */

/*  COORDINATES OF THE IMAGES ON THE EIGENVECTORS  */

	buci
		{
			for(l=0;l<KVVP;l++)
				{
					l1=l+1;
					PSI[i][l]=(S[i][l1]/sqrt((double)SJ[i]))*sqrt((double)VP[l1]);
				}	/* end of for */
		}	/* end of buci */
}	/* end of void */

/*======================================================================*/

void AFCRS()

{

	double p,pc,plim;
	char  chaine[20];

	res=fopen(nomresult1,"a");
	res2=fopen(nomresult2,"a");
	fprintf(res,"\n\n              RESULTS\n");
	fprintf(res,"\n  Eigenvalues\tPercentage\tPercentage");
	printf("\n  Eigenvalues\tPercentage\tPercentage");
	fprintf(res,"\n                      of variance       total");
	printf("\n                of variance       total");
	fprintf(res,"\n %12.4e",VP[0]);
	printf("\n %12.4e",VP[0]);

	pc=0.;
	plim=1.0/((double) NIMA-1.);
	buck
		{
			l=k+1;
			p=VP[l]/vptot;
			pc+=p;
			if(p>=plim) nbsig=l;
			fprintf(res,"\n %12.4e\t  %10.2f\t  %10.2f",VP[l],100.*p,100.*pc);
			printf("\n %12.4e\t  %10.2f\t  %10.2f",VP[l],100.*p,100.*pc);
		}	/* end of buck */
	printf("\n Only the %d first axes seem significant",nbsig);
	fprintf(res,"\n Only the %d first axes seem significant",nbsig);

/*  COORDINATES of the eigenvectors on the images   */

	fprintf(res,"\n\n Coordinates of the axes on the images");
	if(ecranflag==0) printf("\n Coordinates of the axes on the images");
	for(k=0;k<(KVVP+1);k++)
		{
			fprintf(res,"\n axe # %d: ",k);
			if(ecranflag==0) printf("\n axe # %d: ",k);
			buci
				{
					fprintf(res,"  %10.6f",S[i][k]);
					if(ecranflag==0) printf("  %10.6f",S[i][k]);
				}	/* end of buci */			
		}	/* end of for */

/*  COORDINATES of the images on the eigenvectors   */

	fprintf(res,"\n\n Coordinates of the images on the %d axes",KVVP);
	if(ecranflag==0) printf("\n Coordinates of the images on the %d axes",KVVP);
	buci
		{
			fprintf(res,"\n image # %d: ",i+1);
			if(ecranflag==0) printf("\n image # %d: ",i+1);
			buck
				{
					fprintf(res,"  %10.6f",PSI[i][k]);
					if(ecranflag==0) printf("  %10.6f",PSI[i][k]);
					sprintf(chaine,"%10.6f	",PSI[i][k]);
					fprintf(res2,chaine);
				}	/* end of buck */
			fprintf(res2,"%d\n",i+1);
		}	/* end of buci */
	fclose(res); fclose(res2);
}	/* end of void */

/*========================================================================*/

void menu()

{
	printf("\n***************************************************************");
	printf("\n 1 : Compute the images of the factorial axes");
	printf("\n 2 : Filter the images with a reduced number of axes");
	printf("\n 3 : Build a fictitious image");
	printf("\n            10 : PROGRAM EXIT");
	printf("\n***************************************************************");
}	/* end of void */

/*==================================================================*/

void quitter()

{
	printf("\n Remember: Your results are in:\n %s, %s and %s\n",nomresult1,nomresult2,nomresult3);
	printf("\n %d images have been created",nbcre);
	res=fopen(nomresult1,"a");
	fprintf(res,"\n %d images have been created",nbcre);
	fclose(res);
	free(nomdep); free(nomresult1); free(nomresult2); free(nomresult3);
	printf("\n\n");
	exit(0);
}	/* end of void */

/*======================================================================*/

/*  Loop on FI : pixel values in the factorial images  */

void buclfi(ind)

short ind;

{
	short l1;
	double fi,iu;

	l1=k+1;
	bucn
		{
			if(*(SI+n)!=0.)
				{
					fi=0.;
					buci fi+=(double)IX[i][n]*S[i][l1]/sqrt((double)SJ[i]);
					fi/=*(SI+n);
					if(fabs((double)fi)>32767.0)
						{
							printf("\n OVERFLOW %f  @ buclfi",fi);
							quitter();
						}	/* end of if */
					if(ind !=0)
						{
							iu=0.5+fi*(double) facnor;

							min[k]=(iu<(double)min[k])?(short)iu:min[k];
							max[k]=(iu>(double)max[k])?(short)iu:max[k];
							MIN=(iu<MIN)?iu:MIN;
							MAX=(iu>MAX)?iu:MAX;

							*(som+k)+= (long) iu;
							*(ITAB+n)=(short) iu;
						}	/* end of if */
					else FI[n][k]=fi;
				}	/* end of if */
		}	/* end of bucn */
}	/* end of void */

/*======================================================================*/

void AFCPR()

{
	short nombre,nombrex;
	double temp;

	nombre=KVVP; nombrex=nombre+1;
	while(nombrex>nombre)
		{
			printf("\n How many axes do you want to see (<=%d)? ",nombre);
			fflush(stdin); scanf("%hd",&nombrex);
			if(nombrex>nombre) printf("\n Too many ...");
		}	/* end of while */
	if(nombrex>0)
		{
			printf("\n Generic name for the files of these factorial images (4 char. max.): ");
			fflush(stdin);scanf("%s",gener);strcat(gener,"_");
			buckx
				{
					strcpy(nomout[k],gener);
					if(k<9)
						{
							sprintf(inter,"%1d",k+1); strcat(nomout[k],"00");
						}	/* end of if */
					else
						{
							if(k<99)
								{
									sprintf(inter,"%2d",k+1); strcat(nomout[k],"0");
								}	/* end of if */
							else  sprintf(inter,"%3d",k+1);
						}	/* end of else */
					strcat(nomout[k],inter);strcat(nomout[k],ext_axe);
					printf("\n Opening the file of axis %d : %s",k+1,nomout[k]);
					fnom[k]=fopen(nomout[k],"wb");
					fclose(fnom[k]);
					nbcre+=1;
				}	/* end of buckx */

			puts("\n Silence, please ! I'm working hard for you ...");

			/* First loop with a multiplicative factor of 1000 */

			mima(nombrex);
			MIN=32767.0;
			MAX=-32767.0;
			facnor=1000;
			bucnbl
				{
					lect_fich(1);
					buckx	buclfi(1);
				}	/* end of bucnbl */

			/*	if(ecranflag==0) buckx printf("\n axis  %d min = %d max = %d",(k+1),min[k],max[k]);	*/

			/* Compute the maximum multiplicative factor */

			MAX=(fabs((double)MAX)>fabs((double)MIN))?fabs((double)MAX):fabs((double)MIN);
			temp=32767000.0/(MAX+1.0);
			if(ecranflag==0) printf("\n MIN = %f, MAX = %f, Amplification = %f",(MIN/1000.0),(MAX/1000.0),temp);

			/* Second loop with the optimised multiplicative factor */

			mima(nombrex);
			facnor=(long) temp;
			bucnbl
				{
					lect_fich(1);
					buckx
						{
							buclfi(1);
							fnom[k]=fopen(nomout[k],"ab");
							fwrite(ITAB,NBPTS*sizeof(short),1,fnom[k]);
							fclose(fnom[k]);
						}	/* end of buckx */
				}	/* end of bucnbl */

			res=fopen(nomresult1,"a"); res3=fopen(nomresult3,"a");
			fprintf(res,"\n\n Building the factorial images:");
			printf("\n All the images have been multiplied by %f",temp);
			fprintf(res,"\n All the images have been multiplied by %f",temp);
			buckx
				{
					/* Statistics and exit */

					temp=((double)*(som+k))/((double)NLI*(double)NPL);
					if(temp>=0.) temp+=0.5; else temp-=0.5;
					*(moy+k)=(short) temp;
					printf("\n  Statistics of image %s:",nomout[k]);
					printf("\n  Min = %d ,  Max = %d ,  Mean = %d",min[k],max[k],moy[k]);
					fprintf(res,"\n  Statistics of image %s: ",nomout[k]);
					fprintf(res,"\n  Min = %d ,  Max = %d ,  Mean = %d",min[k],max[k],moy[k]);
					fprintf(res3,"%s\n",nomout[k]);
				}	/* end of buckx */
			fclose(res); fclose(res3);
	}	/* end of if	*/
}	/* end of void */

/*======================================================================*/

void AFCRE()

{
	short kr,ja,nombre,flagaxezero;
	double a,temp;

/*  Choose the axes for filtering   */

	flagaxezero=0;
	nombre=KVVP;
	kr=nombre+1;
	while(kr>nombre)
		{
			printf("\n How many significant axes are kept for filtering (<=%d)? ",nombre);
			fflush(stdin); scanf("%hd",&kr);
			if(kr>nombre) printf("\n Too many ...");

		}	/* end of while */
	if(kr>0)
		{
			buckr
				{
					*(ia+k)=nbsig+1;
					while(*(ia+k)>nbsig)
						{
							printf(" Which axis is axis %d  ",k+1);
							fflush(stdin); scanf("%hd",&ia[k]);
							if(ia[k]>nbsig) printf("\n Not significant ...");
						}	/* end of while */
					ia[k]-=1;
				}	/* end of buckr */
			printf("\n Keep axis 0 (Yes=0; No=1)? ");
			fflush(stdin); scanf("%hd",&flagaxezero);


			/*   Opening the files of filtered images  */

			printf("\n Generic name of the files for filtered images (4 char. max.): ");
			fflush(stdin);scanf("%s",gener);strcat(gener,"_");
			buci
				{
					strcpy(nomout[i],gener);
					if(i<9)
						{
							sprintf(inter,"%1d",i+1); strcat(nomout[i],"00");
						}	/* end of if */
					else
						{
							if(i<99)
								{
									sprintf(inter,"%2d",i+1); strcat(nomout[i],"0");
								}	/* end of if */
							else  sprintf(inter,"%3d",i+1);
						}	/* end of else */
					strcat(nomout[i],inter);strcat(nomout[i],ext_filtre);
					printf("\n Creation of filtered image %d: %s",i+1,nomout[i]);
					fnom[i]=fopen(nomout[i],"wb");
					fclose(fnom[i]);
					nbcre+=1;
				}	/* end of buci */

			puts("\n Silence, please ! I'm working VERY hard for you ...");

			mima(NIMA);
			bucnbl
				{
					lect_fich(1);
					buck buclfi(0);
					buci
						{
							bucn
								{
									a=0.;
									buckr
										{
											j=*(ia+k); l=j+1;
											a+=FI[n][j]*PSI[i][j]/sqrt((double)VP[l]);
										}	/* end of buckr */
									if(flagaxezero==0) a++;
									temp=a*SI[n]*SJ[i];
									if(fabs((double)temp)>32767.0)
										{
											printf("\n OVERFLOW %f  @ AFCRE",temp); quitter();
										}	/* end of if */
									ja=(short) (temp+0.5);
									ITAB[n]=ja;
									min[i]=(ja<min[i])?ja:min[i];
									max[i]=(ja>max[i])?ja:max[i];
									som[i]+=(long) ITAB[n];
								}	/* end of bucn */
							fnom[i]=fopen(nomout[i],"ab");
							fwrite(ITAB,NBPTS*sizeof(short),1,fnom[i]);
							fclose(fnom[i]);
						}	/* end of buci */
				}	/* end of bucnbl */

			/* Statistics and exit */

			res=fopen(nomresult1,"a"); res3=fopen(nomresult3,"a");
			fprintf(res,"\n\n Building filtered images:");
			fprintf(res,"\n The following axes are kept:");
			buckr fprintf(res," axis %d",(ia[k]+1));
			if(flagaxezero==0) fprintf(res," and also axis 0:"); else fprintf(res," without axis 0:");
			buci
				{
					temp=((double)*(som+i))/((double)NLI*(double)NPL);
					if(temp>=0.) temp+=0.5; else temp-=0.5;
					*(moy+i)=(short) temp;
					printf("\n  Statistics of image %s:",nomout[i]);
					printf("\n  Min = %d ,  Max = %d ,  Mean = %d",min[i],max[i],moy[i]);
					fprintf(res,"\n  Statistics of image %s:",nomout[i]);
					fprintf(res,"\n  Min = %d ,  Max = %d ,  Mean = %d",min[i],max[i],moy[i]);
					fprintf(res3,"%s\n",nomout[i]);
				}	/* end of buci */
			fclose(res); fclose(res3);
	}	/* end of if */
}	/* end of void */

/*=====================================================================*/

void AFCRF()

{
	short kr,ja,nombre,flagaxezero;
	double a,SJFIC,temp;

/*  Choose tthe axes and the coordinates on these axes  */

	nombre=KVVP;
	kr=nombre+1;
	flagaxezero=0;
	while(kr>nombre)
		{
			printf("\n How many significant axes for this new image (<=%d)? ",nombre);
			fflush(stdin);scanf("%hd",&kr);
			if(kr>nombre) printf("\n Too many ...");
		}	/* end of while */
	if(kr>0)
		{
			res=fopen(nomresult1,"a");
			res3=fopen(nomresult3,"a");
			fprintf(res,"\n\n Building a fictitious image:");
			buckr
				{
					*(ia+k)=nbsig+1;
					while(*(ia+k)>nbsig)
						{
							printf(" Which axis is axis %d  ",k+1);
							fflush(stdin); scanf("%hd",&ia[k]);
							if(ia[k]>nbsig) printf("\n Not significant ...");
						}	/* end of while */
					ia[k]-=1;

					printf("\n Remember the coordinates on this axis of the original images\n");
					j=ia[k]; buci printf("%f ",PSI[i][j]);
					printf("\n Coordinate on this axis of the fictitious image: ");
					fflush(stdin); scanf("\n%lf",&PSIFIC[k]);
					fprintf(res,"\n Coordinate on axis %d: %lf",j+1,PSIFIC[k]);
				}	/* end of buckr */

			printf("\n Do you keep also axis 0 (Yes=0; No=1)? ");
			fflush(stdin);scanf("%hd",&flagaxezero);
			if(flagaxezero==0) fprintf(res,"\n Axis 0 is included"); else fprintf(res,"\n Axis 0 is NOT included");

			printf("\n Remember the intensity normalisation factors of the initial images\n");
			buci  printf("%f ",SJ[i]);
			printf("\n Intensity normalisation factor for this fictitious image: ");
			fflush(stdin); scanf("%lf",&SJFIC);
			fprintf(res,"\n Normalisation: %f",SJFIC);

			/*  Define the output file  */

			nomfictif= (char *) malloc (65);
			printf("\n  Generic name for the file of the fictitious image (8 char. max.): ");
			fflush(stdin); scanf("%s",nomfictif); strcat(nomfictif,ext_fictif);
			printf("\n Creating the file %s...",nomfictif);
			puts("\n Keep cool babe, I'm less relax than you!");
			fprintf(res3,"%s\n",nomfictif);
			fclose(res3);
			fictif=fopen(nomfictif,"wb");
			fclose(fictif);
			nbcre+=1;

			/* Compute and save the fictitious image*/

			mima(1);
			bucnbl
				{
					lect_fich(1);
					buck buclfi(0);
					bucn
						{
							a=0.;
							buckr
								{
									j=*(ia+k); l=j+1;
									a+=FI[n][j]*PSIFIC[k]/sqrt((double)VP[l]);
								}	/* end of buckr */
							if(flagaxezero==0) a++;
							temp=a*SI[n]*SJFIC;
							if(fabs((double)temp)>32767.0)
								{
									printf("\n OVERFLOW %f  @ AFCRF",temp); quitter();
								}	/* end of if */
							ja=(short) (temp+0.5);
							ITAB[n]=ja;
							min[0]=(ja<min[0])?ja:min[0];
							max[0]=(ja>max[0])?ja:max[0];
							som[0]+=(long) ITAB[n];
						}	/* end of bucn */
					fictif=fopen(nomfictif,"ab");
					fwrite(ITAB,sizeof(short),NBPTS,fictif);
					fclose(fictif);
				}	/* end of bucnbl */

			/* Statistics and exit */

			temp=((double)*(som))/((double)NLI*(double)NPL);
			if(temp>=0.) temp+=0.5; else temp-=0.5;
			*(moy)=(short) temp;
			printf("\n  Statistics for this image %s:",nomfictif);
			printf("\n  Min = %d ,  Max = %d ,  Mean = %d",min[0],max[0],moy[0]);
			fprintf(res,"\n  Statistics for this image %s:",nomfictif);
			fprintf(res,"\n  Min = %d ,  Max = %d ,  Mean = %d",min[0],max[0],moy[0]);
			fclose(res);
			fclose(fictif);
			free(nomfictif);
	}	/* end of if	*/
}	/* end of void */

/*========================================================================*/

main()

{
	short MI;

	param();
	AFCCA();
	AFCRS();
	nbcre=0;
	MI=10;
	while(MI!=0)
		{
			menu();
			printf("\n Your choice? ");
			fflush(stdin); scanf("%hd",&MI);
			switch(MI)
				{
					case 1:  AFCPR(); break;
					case 2:  AFCRE(); break;
					case 3:  AFCRF(); break;
					case 10: quitter(); break;
					default: exit(0);
				}	/* end of switch */
		}	/* end of while */
		return(0);
}	/* end of main */

/*========================================================================*/