  {
This program produce radial distribution of the ED image
from Ratio images 
Before start of this program need to run General_RATany4.v

Output files: 2D arrays (coloumn corresponding to delay, and row to pixel number)
<Path\Radial_1D\TRED> - combined files of intensities at set of delays <FULL(LEFT,RIGHT)>.tif
<Path\Radial_1D\TRED> - corresponding combined files of 3x standard of average 
					    <sigmaFULL(sigmaLEFT,sigmaRIGHT)>.tif

Output files: 1D arrays
<Path\Radial_1D\RAW\> - files <F(L,R)><delay>.csv
<Path\Radial_1D\TRED> - combined files of intensities at set of delays <FULL(LEFT,RIGHT)>.csv
<Path\Radial_1D\TRED> - corresponding combined files of 3x standard of average 
					    <sigmaFULL(sigmaLEFT,sigmaRIGHT)>.csv

created 12/02/2001 by V.Lobastov Ver. 01.
}

const
	 P ='D:\DATA\Utilities\Delays\*.ini'; 		 	  //default directory for INI file
//  	 P = 'D:\DATA\04072001\*.ini';	 
	 	 
//------------> Global variables
var
	Data;
    F, L, R, T, B;					 			//field image variables for full, left, and right halves CCD}
    M, M_L, M_R, M_T, M_B;						// MASKS   
    Size_x, Size_y; 		  	 	  		 	//image size	 
    SamName;									//name of the file header   
    Fuss;										//image of the apparatus imperfections
	X0,Y0, RR;											
	im, imm; 	 	  	  	   	   		
	cur, curIm, curMask;						//index, current Image, current MASK
	FULL, LEFT, RIGHT;							//combined RAW 1D data
	sigmaFULL, sigmaLEFT, sigmaRIGHT;			//standard deviations for RAW 1D data
	difFULL, difLEFT, difRIGHT;		 			//combined difference 
	sdifFULL, sdifLEFT, sdifRIGHT;	 			//standard deviations for difference

//   Mask_II;                                     //Image Intensifier Mask Image}
//   Image_I;                                     //Input ED Image}

	 NS,NF, Combo;
	 Sel, delay, FL;
	 Sub, Out, Raw, InIMG;
	 Path, PathR, PathU, PathOUT;
	 Start, Loops, Crit;
	 AppMask;
{
******************************************************************
********************** function Select ***************************
******************************************************************
This function return a full file name for a file 
}	 	 
function SELECT;											
		var
   		   nName; Name, ix;
   		   Choice, nDIR;
		begin

	 		 Name:=FindFirstFile(P, fa_Archive);
	 		 nName:=ExtractFileName(Name);			
	 		 nDIR:= ExtractFilePath(Name);			//extract the path
			 ix:=0;
 		 
			 while StrLen(Name) > 0 do

	  		 	   begin
	  	   		   		Name:=ExtractFileName(FindNextFile) ;
		   				nName:=nName+';'+Name;
					 if ix>100 then halt('Too many Files in Directory');   
					 	ix:=ix+1;
	 				end; 

	 			 if SelectString( 'Select File',nName , Choice ) >= 0 then
	 			 	begin
			 			 SELECT:=nDIR+Choice;
	 				end	
	 			else
	 				begin 
	 	 	  			  WriteInfo('No Selection made! Restart program! ');	
		 	  			  halt;
	 				end; 
		end;
//----------------------> end of function SELECT <-----------------------------------			

{
******************************************************************
********************** procedure HEADER **************************
******************************************************************
This procedure will generate header for output ASCII file
}

procedure HEADER(com, Head, RF);
   var
   	  im;
	  delay;
   begin
			Write(com,'pix'+RF);
   		for im:=NS to NF do
   	   		begin
				   delay:='delay'+Str(im);   
					Write(com, ',', Head + ReadPrivateINIString( Sel, 'FileNames', delay));
	   		end;
					WriteLn(com);
	end;	
//------------------> end of procedure HEADER <------------------------------------

{
******************************************************************
********************** procedure OUTRAT **************************
******************************************************************
This procedure will create output ASCII files in TXT format from
combined image contianing the array of 1D radial intensity data
}

procedure OUTRAT(FULLX,FNAME);

var
   Combo;
   jx, jjx, SX, NFX;

	begin
			SX:=GetXSize(FULLX);
			NFX:=GetYSize(FULLX);

				Combo:=CreateEditor('Out'+FNAME);        		 		
				HEADER(Combo,'','');
		for jjx:=0 to NFX-1 do
     	 	begin		 
         				 	  Write(Combo, Start+jjx) ;
     				for jx:=0 to Integer(SX-1) do	
	 				 	 begin
         				 	  Write(Combo,',', FULLX[jx,jjx]:20:15) ;
 						 end;
         					  WriteLn(Combo);						 
			end;			 	 
				 Save(Combo, PathOUT + FNAME+FL+'.csv');
				 Delete(Combo);
				 Delete(FULLX);
	end;
//----------------------> end of procedure OUTRAT <-----------------



{
*************************************************************************************************
********************** procedure Radial(Image_I,Mask_II,R_Field, Flag) **************************
*************************************************************************************************
{This program produce 1D radial distribution from ED image with rejection of bad points}

procedure Radial(Image_I,Mask_II,R_Field, Flag);				 //Radial(curIm, curMask, F, 'F', cur});

var


//   Flag;   			 	 			 		// flag for FULL, LEFT, RIGHT
   D_Image, F_Image;                  			//Equal radii Mask images}
   Masked_F;
   Masked_FF;                        			//Masked ED Image}
   Num_pix;                                     //number of pixels at certain radius}
   F_pix;                                		//left and right halves}
   Rad_F;                       				//Output Sum of all pixels with equal radius}
   F_sum;                                       //variance image for left and rigth halves}
   FS;										 	//sums}
   F_std;                                 		//standard deviation for left and rigth halves}
   F_rej;								 		//Rejection filter}
   FR_pix;								 		//Number of pixels satisfying chosen criteria}
   FF; 										 	//Number of rejected pixels}
   Rad_FF;								 		//Corrected output sum for Radial distribution}
   Fmin, Fmax;
   Rad_sum;
   MM;
   iLoop;										//number of loops for rejection, loop index
   F_keep;
   Ed1;                                         //Output ASCII file}
   i, j;										//indices   

  begin
  
						 Data := CreateArray( single, 4, Size_x/2-Start) ;
						 if IsNull( Data ) then
				Begin
						 WriteError( 'Insufficient memory to create array!' ) ;
						 Halt;
				end
				else

          for i:=0 to Integer((Size_x)/2-1-Start) do
           begin
		   
               F_Image := ( R_field > Start+i)  and (R_field <= Start+i+1) and (Mask_II) and (Fuss);  {binary ring of current radius}
{-------> Loop begins}						 
						 
		  for iLoop:=1 to Loops do;
		     begin						 
                         F_pix:= Longint(SumOf(F_Image)); 	   	   			   {number of pixels}
                         Masked_FF:=Single(Image_I*F_Image);				   {intensity ring of current radius}
	

                 if (F_pix = 0) then
                    begin
                         Rad_F:=0;
                         F_std:=0;					
                    end
                 else
					begin

                        Rad_F:=Single(SumOf(Masked_FF)/F_pix); 	   			{mean count for current radius}
                        F_sum:=Double(Masked_FF-F_Image*Rad_F);
                        F_sum:=F_sum*F_sum;
						FS:=SumOf(F_sum);					
                        F_std:= Single(sqrt(FS/F_pix));						{STD for a current radius}		

		 			end;	 

{----> calculation of "mean" and STD completed}

					   Fmin:=Rad_F-Crit*F_std;
					if Fmin<0 then Fmin:=0   ;
					   Fmax:=Rad_F+Crit*F_std;

					   F_keep:=(Masked_FF< Fmax) and (Masked_FF> Fmin);		{binary image to keep}
					   F_rej:= not F_keep; 		 	 			 			{binary image to reject}

					   FR_pix:= Word(SumOf(F_keep));						{new number of pixels for current radius}
					   FF:=F_pix-FR_pix;                                    {number of rejected pixels at current loop path}
				   
{-----------> Calculate corrected radial distribution}

						Masked_FF:=Masked_FF*F_keep;
						F_Image:=F_Image*F_keep;
						F_pix:=FR_pix;

                 if (FR_pix=0) then
                    begin
                         Rad_F:=0;
                         F_std:=0;											 
                    end
                    else

					begin
										
                        Rad_F:=Single(SumOf(Masked_FF)/F_pix);				{corrected mean for current radius}
                        F_sum:=Double(Masked_FF-F_Image*Rad_F);
                        F_sum:=F_sum*F_sum;
						FS:=SumOf(F_sum);					
                        F_std:= Single(sqrt(FS/F_pix));						{STD for a current radius}		
					end;	 
			
																														
                        WriteStatus('pix= ', Start+i, ' Rad_F= ', Rad_F, ' F_std= ', F_std, '   Rejected= ',FF);
						
			end;		{----> Loop end}

				Data[0,i]:=Start+i;
				Data[1,i]:=Rad_F;
//------> check for division by zero and calculate standard of average multiplied by factor of 3				
			if ((F_pix-1)= 0) then
			  begin 
				Data[2,i]:=F_std/sqrt(F_pix);
			  end
			  else				
				Data[2,i]:=3*F_std/sqrt(F_pix-1);
				Data[3,i]:=F_pix;		 		
				
Free(F_sum);
              end; 						 

{------> Free memory}

	Delete(D_Image);
	Delete(Masked_FF);
	Delete(Masked_F);
	Delete(R_field);
	Delete(F_Image);

	
{---------> Create Output ASCII file}
    	 Ed1:=CreateEditor(Flag+SamName);
          WriteLn(Ed1,'pix', ',', 'Mean', ',', 'MSEx3', ',', 'Npix') ;
		  
     for j:=0 to Integer(Size_x/2-1-Start) do
	 begin
         Write(Ed1, ShortInt(Data[0,j]), ',', Data[1,j]:20:15, ',', Data[2,j]:20:15, ',', Integer(Data[3,j])) ;
         WriteLn(Ed1);
 	end;	 
		 Save(Ed1, Path+Sub+Raw + Flag + SamName + FL + '.csv');

		 Delete(Ed1);
end;
//---------------> End of procedure Radial <------------------------------

{
 ************************************************************************
 ************************************************************************
 ***                                                                  ***
 ***	                      Main program							  ***
 ***																  ***
 ************************************************************************ 
 ************************************************************************
}

begin

    {-------> Get Input Data Images}
//                   SelectImage( 'select APPARATUS MASK image',Fuss) ;
//                   if IsNull( Fuss ) then halt( 'No "APPARATUS MASK" was selected' ) ;


				  
{------> Default center}
                          	X0:=250.2;
                          	Y0:=256.5;
							

//-----------------------> Initialization parameters<---------------------------------------

	 		Sel:=SELECT;												 // selector for initialization file

{------> Set Center coordinates}
						      				 	  						  
				 X0:=Val(ReadPrivateINIString( Sel, 'Center', 'cX'));	 //read a center coordinates for X
				 Y0:=Val(ReadPrivateINIString( Sel, 'Center', 'cY'));	 //read a center coordinates for Y
			
	 			 Path:= ReadPrivateINIString( Sel, 'Paths', 'Path');	 //read a base path
	 			 PathU:= ReadPrivateINIString( Sel, 'Paths', 'PathU');	 //read a path for Masks
				 PathR:=ReadPrivateINIString( Sel, 'Paths', 'PathR');	 //read a path for polar image				 	

//-------->get apparatus mask				 
				 AppMask:=ReadPrivateINIString( Sel, 'Images', 'AppMask');
 		    	 Open(PathU + AppMask,Fuss);							 
//Show(Fuss,AppMask);
WriteInfo('Xcenter= '+Str(X0)+ '  Ycenter= '+Str(Y0));
						  

                   Size_x:=GetXSize(Fuss);
                   Size_y:=GetYSize(Fuss);

				 
				 Crit:=Val(ReadPrivateINIString( Sel, 'Reject', 'Crit'));	 //read a rejection criteria
				 Loops:=Val(ReadPrivateINIString( Sel, 'Reject', 'Loops'));	 //read a number of loops

				 Start:=Val(ReadPrivateINIString( Sel, 'PixelStart', 'Start'));		//first pixel number 
				 NS:=Val(ReadPrivateINIString( Sel, 'TotalDelays', 'First'));	   //first frame to process
				 NF:=Val(ReadPrivateINIString( Sel, 'TotalDelays', 'Last'));	   //last frame to process
				 
				 				 		 
				 Sub:=ReadPrivateINIString( Sel, 'Paths', 'Sub');		 
				 if (StrLen(FindFirstFile( Path+Sub+'*.csv', fa_Archive)) = 0) then
				 CreateDir(Path+Sub);			 		  				  	   	   //create subdirectory for output
				 	 			 
				 Out:=ReadPrivateINIString( Sel, 'Paths', 'Out');		 
				 PathOUT:=Path+Sub+Out;
				 if (StrLen(FindFirstFile( PathOUT+'*.csv', fa_Archive)) = 0) then				 
				 CreateDir(PathOUT);		 		  				 	   	  		//create directory for output

				 Raw:=ReadPrivateINIString( Sel, 'Paths', 'RAW');
				 if (StrLen(FindFirstFile( Path+Sub+Raw+'*.csv', fa_Archive)) = 0) then				 
				 CreateDir(Path+Sub+Raw);		 		  				 	   	  		
				 				 
				 InIMG:=ReadPrivateINIString( Sel, 'Paths', 'OutIMG');				 				 

//-----------------------> end of Initialization <------------------------------------------

//-----> polar image selection						  
						  RR:='X'+Str(Word(Round(10*X0)))+'y'+Str(Word(Round(10*Y0)));						  

 		    		 	  Open(PathR + RR+'.tif',F);						  
						  WriteStatus('Selected polar image is ' + RR);
//-----> create binary masks for LEFT, RIGHT, TOP, BOTTOM
						  
					 M:= CreateImage( Binary, Size_x, Size_y);
			         M:=0;
						  M_R:= M;
						  M_R[Integer(X0)..Size_x-1,0..Size_y-1]:=1;
						  M_L:= not M_R;
				   
						  L:=F*M_L;
						  R:=F*M_R;
Delete(M);						  
Delete(M_L);
Delete(M_R);

//--------> create output array of 1D ratios 

 				FULL:=CreateArray(Single,NF-NS+1,size_x/2-Start);
 				LEFT:=CreateArray(Single,NF-NS+1,size_x/2-Start);
 				RIGHT:=CreateArray(Single,NF-NS+1,size_x/2-Start);
				
 				sigmaFULL:=CreateArray(Single,NF-NS+1,size_x/2-Start);
 				sigmaLEFT:=CreateArray(Single,NF-NS+1,size_x/2-Start);
 				sigmaRIGHT:=CreateArray(Single,NF-NS+1,size_x/2-Start);				
				
cur:=0;
	   //------> Browse trough files in directory defined in Path+InIMG
   for im:=NS to NF do
   	   begin   
			   delay:='delay'+Str(im);   
	   		   FL:=ReadPrivateINIString( Sel, 'FileNames', delay);
		 	   		
					Open(Path + InIMG + 'NR' + FL +'.tif',curIm); 			  //
		 			Open(Path + InIMG +'NMR'+ FL +'.tif',curMask);
					SamName:=FL;

//-------------> create combined <1D radial in colomns>x<sequential number of processed data>
 
							Radial(curIm, curMask, F, 'F');							
				   					FULL[cur,..]:=Data[1,..]; 
							   sigmaFULL[cur,..]:=Data[2,..];
		 							Delete(Data) ;

							Radial(curIm, curMask, L, 'L');
				   					LEFT[cur,..]:=Data[1,..]; 
							   sigmaLEFT[cur,..]:=Data[2,..];
		 							Delete(Data) ;							

							Radial(curIm, curMask, R, 'R');
				   					RIGHT[cur,..]:=Data[1,..]; 
							   sigmaRIGHT[cur,..]:=Data[2,..];
		 							Delete(Data) ;							
					cur:=cur+1;															
			end;

			Save(FULL, PathOUT+'FULL'+FL+'.tif');
			Save(LEFT, PathOUT+'LEFT'+FL+'.tif');
			Save(RIGHT, PathOUT+'RIGHT'+FL+'.tif');							
			
			Save(sigmaFULL, PathOUT+'sigmaFULL'+FL+'.tif');
			Save(sigmaLEFT, PathOUT+'sigmaLEFT'+FL+'.tif');
			Save(sigmaRIGHT, PathOUT+'sigmaRIGHT'+FL+'.tif');							

//-------------> create ASCII output from combined images 			

		 		OUTRAT(FULL,'FULL');
		 		OUTRAT(LEFT,'LEFT');		 
		 		OUTRAT(RIGHT,'RIGHT');

		 		OUTRAT(sigmaFULL,'sigmaFULL');
		 		OUTRAT(sigmaLEFT,'sigmaLEFT');		 
		 		OUTRAT(sigmaRIGHT,'sigmaRIGHT');


			        PlaySound( 'C:\Lib_V++\Mac.wav' ) ;							
					WriteInfo('Calculation completed!');

//------> free memory					
			Free(Fuss);						
			Free(FULL);
			Free(LEFT);
			Free(RIGHT);	
			Free(sigmaFULL);
			Free(sigmaLEFT);
			Free(sigmaRIGHT);	
end
