From:	SMTP%"pierre.trebbia@univ-reims.fr"  8-OCT-1993 04:29:29.53
To:	ZALUZEC
CC:	
Subj:	sending MacFAC1.0 (message 1)

Date: Fri, 8 Oct 93 08:12:58 +0100
To: Zaluzec@ANLEMC.MSD.ANL.GOV
From: pierre.trebbia@univ-reims.fr (Pierre TREBBIA)
X-Sender: trebbia@centre01
Subject: sending MacFAC1.0 (message 1)
X-Attachments: :Pierre No 1:5530:MACFAC.ABS:

Hi Nestor,
Well, I'm going to send you the following files:
1) MACFAC.ABS is attached to this message
2) MACFAC.DOC is attached to a 2nd message sent just after message 1
3) MACFAC.SRC is attached to a 3rd message sent just after message 2
4) MACFAC.SEA.HQX (whole package with demo and compiled program) is
attached to a 4th message sent just after message 3.
If you have any trouble (MACFAC.SRC is rather long and may be split into 2
messages...) in receiving these files, tell it to me!
Cheers,
Pierre
--------------
Pierre Trebbia, LASSI/GRSM, Universite de Reims, F51062 Reims cedex, France
Email : pierre.trebbia@univ-reims.fr
Phone : (33) 26 05 33 62
FAX : (33) 26 05 32 50
Title: MacFAC 1.0
Keywords: multivariate statistics, data analysis, processing, filtering
Computer: MacIntosh II family 
Operating system: system 7
Programming language: THINK C 5.0
Hardware requirements: math co-processor
Author(s): Noel BONNET & Pierre TREBBIA
Correspondence: Faculte des sciences, Universite de Reims, B.P. 347, F51062
Reims cedex FRANCE.
Fax: (33) 26 05 32 50, e-mail: pierre.trebbia@univ-reims.fr
Abstract:
The purpose of MacFAC is to extract orthogonal sources of information from
data (images, spectra,...). See any basic statistic book for more
information about multivariate statistics and Factorial Analysis of
Correspondence (FAC). MacFAC is given with a test data set in order to let
you be sure that it works properly on your own computer. Read the text file
READ.ME which is embedded in the file MacFAC.sea.hqx for more information
on how the test must be done.


From:	SMTP%"pierre.trebbia@univ-reims.fr"  8-OCT-1993 04:06:22.37
To:	ZALUZEC
CC:	
Subj:	sending MACFAC (message 2)

Date: Fri, 8 Oct 93 08:13:50 +0100
To: Zaluzec@ANLEMC.MSD.ANL.GOV
From: pierre.trebbia@univ-reims.fr (Pierre TREBBIA)
X-Sender: trebbia@centre01
Subject: sending MACFAC (message 2)
X-Attachments: :Pierre No 1:5530:MACFAC.DOC:

Here is MACFAC.DOC
--------------
Pierre Trebbia, LASSI/GRSM, Universite de Reims, F51062 Reims cedex, France
Email : pierre.trebbia@univ-reims.fr
Phone : (33) 26 05 33 62
FAX : (33) 26 05 32 50
Title: MacFAC 1.0
Documentation:
The purpose of MacFAC is to extract orthogonal sources of information from
data (images, spectra,...).
See any basic statistic book for more information about multivariate
statistics and Factorial Analysis of Correspondence (FAC).
See also (for example) P. Trebbia, N. Bonnet, Ultramicroscopy 34 (1990) 165,
P. Trebbia, C. Mory, Ultramicroscopy 34 (1990) 179,
N. Bonnet et al, Proc. MAS Meeting, Los Angeles (1993) p. S212

MacFAC is given with a test data set in order to let you be sure that it
works properly on your own computer. Read the text file READ.ME which is
embedded in the file MacFAC.sea.hqx for more information on how the test
must be done.


From:	SMTP%"pierre.trebbia@univ-reims.fr"  8-OCT-1993 05:56:05.88
To:	ZALUZEC
CC:	
Subj:	sending MACFAC (message 3)

Date: Fri, 8 Oct 93 08:14:38 +0100
To: Zaluzec@ANLEMC.MSD.ANL.GOV
From: pierre.trebbia@univ-reims.fr (Pierre TREBBIA)
X-Sender: trebbia@centre01
Subject: sending MACFAC (message 3)
X-Attachments: :Pierre No 1:2:MACFAC.SRC:

here is MACFAC.SRC
--------------
Pierre Trebbia, LASSI/GRSM, Universite de Reims, F51062 Reims cedex, France
Email : pierre.trebbia@univ-reims.fr
Phone : (33) 26 05 33 62
FAX : (33) 26 05 32 50
/*============= PROGRAM MACFAC.C ============================================
	Version 1.0: Wed 6 october 1993 15: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.
	French version of Microsoft Excel needs a comma as decimal point. This
translation is made at line
	473 in this listing. Discard this line if you don't need this translation.
==========================================================================
	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 NBIMA 10
#define NBAXE 10
#define NBPTS 128

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

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#include <string.h>
/*	#include <malloc.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,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\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? ");
	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);

	printf("\n     How many images to be processed (<=%d) in %s? ",NDIM,nomdep);
	fflush(stdin);
	scanf( "%hd",&NIMA);
	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>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;

	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****************************************");
	temp=((long) NPL*(long)NLI)/(long)NBPTS;
	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);
	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);
	fclose(res);
	printf("\n Sum of the eigenvalues after diagonalisation = %f +
%15.7e",VP[0],vptot);
	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;
	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.;
	buck
		{
			l=k+1;
			p=VP[l]/vptot;
			pc+=p;
			if(pc<=0.99) nbsig=l+1;
			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 are significant",nbsig);
	fprintf(res,"\n Only the %d first axes are 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]);
					for (j=0; j<strlen(chaine); j++) if (chaine[j] == '.') chaine[j] =
',';	/* french coding of the decimal point	*/
					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=nbsig; 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 */
	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 void */

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

void AFCRE()

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

/*  Choose the axes for filtering   */

	flagaxezero=0;
	nombre=KDIM; if(nombre>nbsig) nombre=nbsig;
	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 */
	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 void */

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

void AFCRF()

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

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

	nombre=KDIM; if(nombre>nbsig) nombre=nbsig;
	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 */
	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 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 */

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