{This program will generate a combined image (array of ratio images) 
in a grid with specified step. 



created 13/01/2002 by V.Lobastov Ver. 08}

const

{	 			   	  			 	  	 	 		   		 		 										--->>OLD
	 Path ='D:\DATA\';								//Base path
	 PathU ='D:\DATA\Utilities\'; 		 			//default Utilities directory
	 PathR ='D:\DATA\Utilities\XYarrayNEW\';		//default directory for polar images    		  
	 DM ='AppCOT.tif';								//default Mask
	 														  	  				  					  		<<---OLD
}
	 	  //P ='D:\UED3_Data_processed\Utilities\Delays\*.ini'; 		 	  //default directory for INI file
		  P ='D:\UED4_Data_raw\Utilities_UED4\Delays\*.ini'; 		 	  //default directory for INI file
  	 	 
		  Prompt = 'Select an image for calculation' ;	  
	 
		  CommonDlgsLib = 'VppCommonDlgs.dll' ;		// declaration of CommonDialog constant 

//Windows API constants
  		  wm_Close 		= 16 ;
		  wm_SysCommand = $0112 ; 
		  sc_Close 		= 61536 ;		  	  		// command to close Excel
		  

		  { Declaration of external DLL functions }
		  function OpenDialog( FileName,Title,Filter,DefaultExt:pointer ) : integer ; external CommonDlgsLib ;		  
		  function FindWindow( ClassName,WindowName:pointer ) : integer ;
		  external 'User32' ; name 'FindWindowA' ;

		  procedure SendMessage( Handle,Message,wParam,lParam:integer ) ; 
		  external 'User32' ; name 'SendMessageA' ;		  
		  
		  
		  
		  
		  
{************************************************************************}
{ ***	                      Global variable declaration			  ***}									   
{************************************************************************}

var
   cP;
   X0,Y0; 									//initial coordinates for center   
   X0_new,Y0_new; 							//output coordinates for center   
   S; 										//Search radius
   cur, Ed;									//
   i,j, m,n, q;										//indices
   Fnam;									//extended file name for .CSV file
   extFile;									//file name
   I_Mask; 			   		   				//default Mask image
   Size_x, Size_y, HALF;                     //Size of image in X and Y, and half of the size
   Img, iNam;								 //Experimental image, extended and short names
   iC;										 //Return image
   NRAT;									 //ratio of experiment/converted
   RR;			   							 //appended image name
   PAT;										 //array of variance
   Data, cellY;						 		 //Array of 1D intensities, colomn in ASCII file
   ST,k,iT;								 	 //Search step
   sR, eR;							 		 //search limits
   arSz;
   minS;									 //minimum step size
   Combo, Variance;

   Sel;
   AppMask;
   Path;
   PathU;
   PathR;
   PathC;
   RXY;
{
******************************************************************
********************** 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);				  	 	 		//Find the first ini file in the delays directory 
	 		 nName:=ExtractFileName(Name);										//Convert the file name and the extension to a string
	 		 nDIR:= ExtractFilePath(Name);										//Return the drive and directory portions of a file name including the trailing backslash
			 ix:=0;																//initialize the variable ix to value of zero
	  		 
			 while StrLen(Name) > 0 do											//If there is a file in the directory, do

	  		 	   begin
	  	   		   		Name:=ExtractFileName(FindNextFile) ;					//Finde the next ini file in the delays directory
		   				nName:=nName+';'+Name;									//Augment the nName string with the next file name
					 if ix>100 then halt('Too many Files in Directory');   		//Stop, if there are more than 100 files
					 	ix:=ix+1;			  	   		 						//Increase the counter ix by 1
	 				end; 

	 			 if SelectString( 'Select File',nName , Choice ) >= 0 then		//If a selection among nName strings is made from the dialog box, then
	 			 	begin
			 			 SELECT:=nDIR+Choice;		  		   	 	  			//Initialize SELECT to the selected ini file
	 				end	
	 			else
	 				begin 
	 	 	  			  WriteInfo('No Selection made! Restart program! ');	//Inform the user that no selection was made
		 	  			  halt;					  						 		//abort the program
	 				end; 
		end;
//----------------------> end of function SELECT <-----------------------------------			



{
 ************************************************************************
 ***                                                                  ***
 ***	                      Procedure RAD							  ***
 ***																  ***
 ************************************************************************ 
This procedure generates a polar coordinate image with center passed as input parameters 
 }  

procedure RAD(X_off,Y_off);

var
  Image_R;                                      {Polar-field image (equal radius field)}
  i, j;											{indices}
  Full;
begin
   Full:='x'+Str(Word(Round(10*X_off)))+'y'+Str(Word(Round(10*Y_off)));	  		//Rounds the center coordinates, converts the expression to word and then to a string
   if FileExists(PathR+Full+'.tif') then										//Check if a polar image already exists
   begin
      WriteStatus('File Already Exists  '+ Full);								//Report that the file exists 
   end   
   else   
   begin
      WriteStatus('X_off= '+Str(X_off)+ '  Y_off= '+Str(Y_off));				//display the center coordinates 
      Image_R:= CreateImage(Double, Size_x,Size_y);								//create a radial image of type double
      for i:= 0 to (Size_x-1) do	
         for j:= 0 to (Size_y-1) do
            Image_R[i,j] := Sqrt(Sqr(i-X_off)+Sqr(j-Y_off));					//adding the value of radial distance to each pixel in the image
      Save(Image_R,PathR+Full+'.tif');											//save the radial image in the output directory 
      Free(Image_R);															//free the memory held by the variable 
   end;
end;  
  
{---------> end of procedure RAD}

   
{************************************************************************}
{ ***	                      Function GetOpenFileName 				  ***}									   
{************************************************************************}
{This function returns a full name of selected file }		  
   
function GetOpenFileName(Title,Filter,DefaultExt ) ;
{ Open the common Open dialog box and get a file name }
{ Quit if Cancel is pressed }

  const
  	   BufSize = 255 ;
  
  var
  	   FileName ;

	   begin	 	 
	   		SetDir(Path);	  	  	  	  		 	   		  			  		//Change to the path directory
  	 		FileName := StringOfChar( ' ',BufSize ) ;							//Filename is intitialized to a string of whitespaces
  			
		 if OpenDialog( FileName,Title,Filter,DefaultExt ) <> 0 then			//Open a dialog box, if a file is selected, then
     	 	GetOpenFileName := Trim( FileName )			   	  					//Trim the whitespaces around the selected filename
  		else
  			halt('File is NOT selected!!!');  									//if not file was selected, abort the program.
end;  { GetOpenFileName }



{************************************************************************}
{ ***	                      function ImgCENT(X2,Y2, RR, sR,eR)	  ***}									   
{************************************************************************}

//------>this function returns an image generated from 1D data
//------>input parameters are center coordinates (X2,Y2), file name for polar image, start and end of the caluculated radius

function ImgCENT(X2,Y2, RR, StartR, endR);

var
   R;	 									 //1D-->2D converted image image 
   R_field,F_image;							 //polar image, mask image of equal radii          
   
	begin	  							

//-----> polar image selection						  
						  RAD(X2,Y2);				  			 	   	  		//This procedure generates a polar coordinate image with center passed as input parameters
 		    		 	  Open(PathR + RR+'.tif',R_field);			  			//Open the radial coordinate image
						  WriteStatus('Selected polar image is ' + RR);			//Display a message in the status bar
						  				
//-----> create arrays				   
 			   

				   R:=CreateImage(single,Size_x,Size_y);		   				//create an image of same size as the input mask image		   				   
         		   R:=0;														//initialize all pixels to value of zero

					for i:=startR to endR do									//loop through the given radial distance
           				begin		   
               				 F_Image := ( R_field > i) and (R_field <= i+1);  	//binary ring of current radius
			   				 R:=R+F_image*Data[0,i];   	   			   			//create 2D images from 1D data
//			  				 Show(R,RR);			   
							 WriteStatus(RR+'  r= '+Str(i));					//Display a message in the status bar
		  				end;
				   			 R:=R*I_Mask*(not (R = 0));			  				//mask the 2D image with the image intensifier mask
							 WriteStatus(RR); 	   								//Display a message in the status bar
//							 Show(R,RR);
							 ImgCENT:=R;
							 							 				   
//----------->Free memory
				Delete(R); 
				Delete(F_Image); 
				Delete(R_field);

		end;		
//--------------------> end of function ImgCENT <--------------------------------- 

{************************************************************************}
{ ***	                      function DAT(HALF)				  	  ***}									   
{************************************************************************}
{This procedure will read an ASCII file into the array Data}
procedure DAT(HALF);
//----->
	var
   	   WrkSht; 		   							 						  				//name for Excel worksheet
   	   vCh;									 											// DDE channels 		
		begin
				   Data:=CreateArray(single,1,HALF);					   				//create an array of type single with 1 column and HALF rows
				   
				   WrkSht := Trim(ChangeFileExt( extFile, '' )) ;						//trim away the file extension of the radial file
				   WriteStatus(extFile);

	 		   Execute('C:\Program Files\Microsoft Office\Office11\EXCEL.EXE '+Fnam);	//Run Excel from within the module to open the radial file
			  // WrkSht := Trim(ChangeFileExt( extFile, '' )) ;				 			
				   
				   vCh := DdeInitiate( 'Excel',WrkSht ) ;  	  							//Initiate a conversation with Excel
	
	  			if vCh = 0 then halt('Communication problem!!!') ;						//if communication could not be established, abort the program
  
  	  			   for i:=1 to HALF do
	  	  		   	   begin
			   		   		cellY:= 'R'+Str(i+1)+'C2';			  						//Absolute adress of the cell in Excel: Row i, Column 2
			   				Data[0,i-1] := DdeRequest( vCh, cellY ,fmt_Number ) ;		//Request a number from cellY and read it into the array Data
  		   				end;	 
  					DdeTerminate(vCh);	   			   			  			  			//Close the conversation with Excel
  					SendMessage( FindWindow( 'XLMAIN',0 ),wm_SysCommand,sc_Close,0 ) ;	//???
		// Halt('Checkpoint');
		end;			
//--------->  

{************************************************************************}
{ ***	                      function iRAT(Fst, Sec)				  ***}									   
{************************************************************************}
//------> this function will generate a ratio of FIRST image over SECOND image

function iRAT(Fst, Sec);
		 begin
		 	  iRAT:= Fst*(not (Sec = 0))*(not (Fst = 0))/(Sec+(Sec = 0)) 	  			//avoid division by 0
		 end;
//--------------------> end of function iRAT <---------------------------------

{************************************************************************}
{ ***	                      function minPAT						  ***}									   
{************************************************************************}
//------> this function return a coordinate of the minimum of an array

function minPAT;
	var
	   ind;
	   minV;	 
		 begin
		 	minV:=MinOf(PAT[2,..]);  		 
				for ind:=0 to arSz-1 do
					begin	
						if (PAT[2,ind]=minV) then minPAT:=ind;				 
					end;	

		 end;
//--------------------> end of function minPAT <---------------------------------

{************************************************************************}
{ ***	                      function VARI(iM)				          ***}									   
{************************************************************************}
//----> this function return variance over an image for all non zero pixels

function Vari(iM);
	var
	   pix;
	   MeanImg;
	   MS;
		 begin
		 	   	  		MS:= not(iM = 0); 	   	  														  
						pix:= Single(SumOf(MS));							{count all nonzero pixels}						
                        MeanImg:=Single(SumOf(iM)/pix); 	   			{mean count for selected Image}
						Vari:=SumOf(MS*sqr(iM-MeanImg));						

		
		 end;
//--------------------> end of function VARI <---------------------------------

{************************************************************************}
{ ***	                      procedure ITER(Step)			          ***}									   
{************************************************************************}

//this function return a COMBO image of a series calculated images 

function ITER(Step,X1,Y1);
		  var
		  	 jX, jY; 		   				 //indices
   			 X, Y;							 //current coordinates
			 iRR;							 //current Name
   			 NROI;							 //Cut ROI
			 v,w;							 //indices
		begin	

				v:=0; 			 						
				for jX:=-S to S do						   		   		  		//step through the grid parameter limits
					begin	
							X:=X1+jX*Step;										//X is the shifted center x coordinate  
				w:=0;
						for jY:=-S to S do										//step through the grid parameter limits
							begin	

									Y:=Y1+jY*Step;								//Y is the shifted center y coordinate 
						  			iRR:='x'+Str(Word(Round(10*X)))+'y'+Str(Word(Round(10*Y)));	  	//Returns a string of the currently used center coordinates
									iC:=ImgCENT(X,Y,iRR,sR,eR);					//this function returns an image generated from 1D data
									NRAT:=iRAT(Img,iC);							//this function will generate a ratio of ratio image over the 2D image
																	
									NROI:=NRAT[(X-eR)..(X+eR),(Y-eR)..(Y+eR)];	//Define the outermost region of the NRAT image						
				Combo[v*(2*eR)..(v+1)*(2*eR),w*(2*eR)..(w+1)*(2*eR)]:=NROI;		//Paste the NRAT image into the Combo image
				Show(Combo,'x'+Str(Word(Round(10*X0)))+'y'+Str(Word(Round(10*Y0)))+'s'+Str(Word(Round(Step*10)))+'x'+Str(2*(iT-k)+1));			//Show the combo image
				SetDisplayZoom(Combo,25);										//reduce the size of the image on the desktop					
//			WriteInfo(RR);	
//						 			Show(NROI, 'ROI'+RR);
//									SetDisplayZoom(NRAT,50);									
						 			SetDisplayRange( Combo , 0.9, 1.1 );		//Set the display range for the combo image
									PAT[0,cur]:=X;	 	   	 	  	  			//Write the shifted x coordinate into the first column of the array
									PAT[1,cur]:=Y;								//Write the shifted y coordinate into the second column of the array
									PAT[2,cur]:=Vari(NRAT);						//this function returns the variance over an image for all non zero pixels											
									cur:=cur+1;									//increase the line index by one 
								Free(NRAT);										//Free the memory held by NRAT
								Free(NROI);										//Free the memory held by NROI
					w:=w+1;														//increase the index by one
						 	end;			  
					v:=v+1;					  									//increase the index by one	
					end;	 
				ITER:=PAT;	 													//The function returns the updated array PAT
			end;
//--------------------> end of procedure ITER <---------------------------------		 
		 
{
 ************************************************************************
 ************************************************************************
 ***                                                                  ***
 ***	                      Main program							  ***
 ***																  ***
 ************************************************************************ 
 ************************************************************************
}		  

		  begin		   
		  
//-----------------------> Initialization parameters<---------------------------------------
	 		Sel:=SELECT;												  			// selector for initialization file
			
	 			 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				 	
				 PathC:=ReadPrivateINIString( Sel, 'Paths', 'Center');	 			//read a path for output file				 
				 

				 					  
//--------> Get image file
			  	
			   iNam:=GetOpenFileName( 'Open','Image Files (*.tif)|*.tif','' );		//select a ratio image file from the browser
			   Open(iNam, Img);			  				  							//open the ratio image
			 

//--------> Get 1D file			    			   
	
			   FNam:=GetOpenFileName( 'Open','Radial Files (*.csv)|*.csv','' );		//select the 1d radial csv file from the browser
			   extFile:=ExtractFileName(Fnam);			   					 		//Return the name and extension portions of a file name


//-------->get apparatus mask				 
				
				 AppMask:=ReadPrivateINIString( Sel, 'Images', 'AppMask'); 			//read the filename of the apparatus mask	 
 		    	 Open(PathU + AppMask,I_Mask);							 			//read the apparatus mask from disk
				 				 
                   Size_x:=GetXSize(I_Mask);				  						//Return the x-dimension of the size definition imag
                   Size_y:=GetYSize(I_Mask);					  					//Return the y-dimension of the size definition imag
				   HALF:=Integer(Size_x/2);											//Assign the center point in the x-dimension
				   DAT(HALF);						  								//This procedure will read an ASCII file into the array Data				   

//------> Default center

				 X0:=Val(ReadPrivateINIString( Sel, 'Center', 'cX'));	 			//read a center coordinates for X
				 Y0:=Val(ReadPrivateINIString( Sel, 'Center', 'cY'));	 			//read a center coordinates for Y

					GetNumber('Enter the X0 coordinate',X0);  						//Display a dialog box and prompt the user to enter a number for the x center coordinate
					GetNumber('Enter the Y0 coordinate',Y0);						//Display a dialog box and prompt the user to enter a number for the x center coordinate

//------> Search limits		  

		  		 	   		sR:=40;
							eR:=250;							
					GetNumber('Enter the search RADIUS (pixel)',eR);				//Display a dialog box and prompt the user to enter a number

							S:=2;
							minS:=0.5;	 
					
					GetNumber('Enter the GRID parameter',S);						//Display a dialog box and prompt the user to enter a number
					GetNumber('Enter the MINIMUM step (pixel)',minS);				//Display a dialog box and prompt the user to enter a number
					
					RXY:='x'+Str(Word(Round(10*X0)))+'y'+Str(Word(Round(10*Y0)))+'s'+Str(Word(Round(10*minS)));  //CenterCoordinates and Step Size

									
//------>Iterations
	
				   			iT:=1;
					GetNumber('Enter the number of ITERATIONS ',iT);				//Display a dialog box and prompt the user to enter a number
												
							arSz:=Integer(iT*sqr(2*s+1));	  						//Define the Array Size
						PAT:=CreateArray(single,3,arSz);							//Create an array of type single with 3 columns and arSz rows
						PAT:=1111111111;											//Initialize the array ???
						cur:=0;														//intialize the variable cur (line index of the array PAT

						
												
				for k:=1 to iT do													//step through the number of iterations	
					begin
Combo:=CreateImage(single, Integer((2*S+1)*2*eR+1),Integer((2*S+1)*2*eR)+1);		//Create an image of type single with size specified by the grid parameter and the search radius
Combo:=1;																			//Initialize all pixels to value 1

					ST:=minS*(2*(iT-k)+1);											//Initialize the variable search step
					WriteStatus(ST);												//Display the value of ST in the status bar
								
								PAT:=ITER(ST,X0,Y0);								//this function returns a COMBO image of a series calculated images and updates the array PAT
								cP:=minPAT;											//this function returns a coordinate of the minimum of the variance values in the array PAT
								X0_new:=PAT[0,cP]; 									//assing the x coordinate belonging to the minimum variance to X0
								Y0_new:=PAT[1,cP];									//assing the y coordinate belonging to the minimum variance to Y0
				
				    //Save(Combo,Path+PathC+RXY+'.tif');								//Save the combo image
					Free(Combo);													//Free memory
					end;	

	
					
Variance:=CreateImage(single, Integer((2*S+1)),Integer((2*S+1)));				//Create an image containing variance information
q:=0;
							for m:=0 to Integer((2*S)) do
									 begin
									 for n:=0 to Integer((2*S)) do
									 	 	  begin
											  Variance[m,n]:= PAT[2,q];
											  q:=q+1;
											  end;
									 end;
					
					Show(Variance,'VarianceImage_'+RXY);				   			//Show the Variance image	
					SetDisplayRange(Variance, 140, 300 );							//Set the display range for the combo image
					SetDisplayZoom(Variance,2000);								//increase the size of the image on the desktop	
					
					
					Ed:=CreateEditor('Variance');
							for i:=0 to arSz-1 do
								begin	
									WriteLn(Ed, PAT[0,i], ', ',PAT[1,i], ', ',PAT[2,i])				 
								end;


{
						  			RR:='x'+Str(Word(Round(10*X0)))+'y'+Str(Word(Round(10*Y0)));
									iC:=ImgCENT(X0,Y0,RR,0,HALF-1);
				
						 	 		NRAT:=iRAT(Img,iC);
									Save(NRAT,Path+PathC+'F'+RXY+'.tif');
						 			Show(NRAT, 'F'+RR);
//									SetDisplayZoom(NRAT,50);
						 			SetDisplayRange( NRAT , 0.9, 1.1 ) ;
						WriteInfo('x=  '+ Str(PAT[0,cP])+'  y= '+ Str(PAT[1,cP]) + '  var= '+Str(PAT[2,cP]));			
}													
				Delete(Data);
				Delete(I_Mask);				
				Delete(Img);
				Free(iC);
				Free(PAT);   
				Free(NRAT);
			        PlaySound( 'C:\Program Files\Digital Optics\V++\Library\Mac.wav' ) ;	
					WriteInfo('Calculation completed!                                         '+ 
							'The coordinates corresponding to the mimum variance are:         '+ 
					       	'X = '+Str(X0_new)+', Y = '+Str(Y0_new) + 
						  	'The old coordinates were:'+       
						   	'X = '+Str(X0)+', Y = '+Str(Y0));							
		  end
