{
This program will evaluate all images with chosen criteria of Standard Deviation
and create Average and Standard Deviation Images in selected directory <Path>

Program require ANY image of the same size as processed images (usually some any mask image)

Program DO NOT require initialization file
After processing directories MUST be cleaned and 
all temporary files should be deleted 

File names for temporarily files in parent directory are

<"M"><Fle header><_><index>.tif
<"B"><Fle header><_><index>.tif 
Sum.tif
AB.tif


The input parameters requested from operator:

- 'Select IMAGE for size definition' --> enter the file header 
- Browse and select directory where images are stored
- 'Enter the FILE HEADER' --> enter the file header convention <FILE HEADER><_><Index>.tif
- 'Enter the FIRST frame' --> enter the index number for start
- 'Enter the LAST frame' --> enter the index number for end

- 'Enter the NUMBER of LOOPS' 	   --> number of circles of rejection 
- 'Enter the REJECTION CRITERIA'   --> rejection criteria in sigmas 
- 'Enter the REJECTION TRESHOLD'   --> threshold for saturated pixels or 
  		 	 		   			   	   pixels with anomalously high amplitude

The program (V++) has a memory leak which lead to run out of memory
It is advised to close and restart V++ after processing 3-4 directories in a row

All averaged and standard deviation images stored for each loop
Naming convention for output files:
<File header><"_"><av><loop index> 	--> averaged file
<File header><"_"><sd><loop index>	--> standard deviation file

Usually the last averaged file is used for further processing
									   
Created by Vladimir Lobastov 13.01.2002}



const
  	 	  CommonDlgsLib = 'VppCommonDlgs.dll' ;	 	{DLL for common dialog} 
		  InPath='D:\DATA\28082003\';			 	{default path for images}
		  Path2='D:\DATA\Utilities\';				{default path for initial mask image}
		  Qwst='Do you want to use Patch for saturated pixels?';
//		  InPath = 'D:\DATA\';
		  		  
		  function BrowseFolders( Path:pointer ) : integer ; external CommonDlgsLib ;
		  

var 
	CR;		  							{rejection criteria}
	ind, jnd, k;						{indices}
	Loops;
	F,N;		  						{first and last number of images to process}
	FHead;	  							{file header}
	xSize, ySize;						{sizes of the images}	
	Img;								{image for size definition}
	CH;
	Path;								{path for images}
	TR;

{************************************************************************}
{ ***	                      Procedure InitMASK	   				  ***}									   
{************************************************************************}
{This procedure will create array of mask and new data images }
procedure InitMASK;
var
	MI0, BI;
	ind;
	i;	  
  begin
  	   	  	if Query(Qwst) = id_Yes then  
			begin			   
			   Open( Path2+'SpotPatch.tif', MI0) ; {create and mask images}
			end   
			else
			   Open( Path2+'BlankBinaryInv.tif', MI0) ; {create and mask images}			   
			   {MI0:=1;	}		   			   		  		   {initiate mask imge as} 

{---->}	for ind:=F to N do	   
		   begin
			   i:='_'+Str(ind);
   			   Save(MI0, Path+'M'+FHead+i+'.tif');  		   {save mask image}			
 			   Open(Path+FHead+i+'.tif',BI);
			   Save(BI, Path+'B'+FHead+i+'.tif');
			   WriteStatus('Initiating... '+Str(ind));
			   Delete(BI);
			   Delete(i); 				   					
		   end;
			   Delete(MI0);	   		   
  end;		   					
  {-----------------------> end of procedure InitMASK}

{************************************************************************}
{ ***	                      Procedure ABI	  		   				  ***}									   
{************************************************************************}
{This procedure calculates mean image}
procedure ABI(k);
var 
   AB,AVG;
   NI,MI,BI;
   i;		  
   Zeros;
   FNAME;
begin
	 	 Open( Path2+'BlankSingle.tif', AB) ; {open blank Single image}
		
{---->}for ind:=F to N do
		begin
			   i:='_'+Str(ind);			   
 					Open(Path+'B'+FHead+i+'.tif',BI);					
			   		AB:=AB+BI;
					WriteStatus('Summed...  '+ Str(ind));
			   	   	Show(AB, 'Summing');			
			     Delete(BI);
				 Delete(i);
		end;
				 Save(AB,Path+'AB'+'.tif');
				 Delete(AB);

	 	 Open( Path2+'BlankShortInt.tif', NI) ; {open blank ShortInt image}
{---->}
	for ind:=F to N do
		begin		 
			   i:='_'+Str(ind);			   
					Open(Path+'M'+FHead+i+'.tif',MI);					
					NI:=NI+MI;
Show(NI,'NI');
			     Delete(MI);
				 Delete(i);									   
		end;			
					Zeros:= ((NI)= 0);
				Open(Path+'AB'+'.tif',AB);						
				AVG:=Single(AB/(NI+Zeros))*(not Zeros);
			    Save(AVG, Path+FHead+'_av'+Str(k)+'.tif');											
		
			Delete(AVG);						
			Delete(Zeros);
			Delete(AB);
			Delete(NI);
	end;
{------------------------------------>end of procedure --ABI--}  		
	
{************************************************************************}
{ ***	                      Procedure STD	  		   				  ***}									   
{************************************************************************}
{This procedure calculates standard deviation image}	
procedure STD(k);
var
		  A;
		  i;
		  NI,BI,MI,SD;
		  Sum;
		  Zeros;		 
  begin

  	    Open( Path2+'BlankSingle.tif', Sum) ; {open blank Single image}
    	Open(Path+FHead+'_AV'+Str(k)+'.tif',A);										
{---->}								   				  
	for ind:=F to N do
		begin
			   i:='_'+Str(ind);
			   WriteStatus(i);
			   	     		  Open(Path+'B'+FHead+i+'.tif',BI);
			   	     		  Sum:=Sum+(BI-A)*(BI-A);
					Delete(BI);
					Delete(i);				
		end;
					Delete(A);
				 Save(Sum,Path+'Sum'+'.tif');					
				 Delete(Sum);		  


	 	 Open( Path2+'BlankShortInt.tif', NI) ; {open blank ShortInt image}
{---->}
	for ind:=F to N do		
		begin
			   i:='_'+Str(ind);					
				   	 		  Open(Path+'M'+FHead+i+'.tif',MI);					
				   			  NI:=NI+MI;	   
					 Delete(MI);
					 Delete(i);				
		end;	
				Open(Path+'Sum'+'.tif',Sum);						
			Zeros:= ((NI-1) <= 0);
			CH:=((NI+Zeros-1) <=0);
			SD:=Single(sqrt(Sum/((NI+Zeros-1+CH)))*(not Zeros)); {floating point error}
			    Save(SD, Path+FHead+'_sd'+Str(k)+'.tif');							
			Delete(SD);
			Delete(Sum);
			Delete(Zeros);
			Delete(NI);
			Delete(CH);
  end;
{------------------------------------>end of procedure --STD--}  	

{************************************************************************}
{ ***	                      Procedure Reject 		   				  ***}									   
{************************************************************************}
{This procedure will reject pixels not satisfieng chosen criteria of 3 sigmas}

procedure Reject(k);
var
		  A,S;
		  j;
		  MI,BI,RJ;
begin
{---->}	
	for ind:=F to N do	
		begin
			j:='_'+Str(ind);  
   		 		Open(Path+'B'+FHead+j+'.tif',BI);
			    Open(Path+FHead+'_SD'+Str(k)+'.tif',S);							
			    Open(Path+FHead+'_AV'+Str(k)+'.tif',A);																   
			   		   RJ:= ((A-CR*S) < BI) and (BI< (A+CR*S)) and ( BI< TR );
					Delete(A);
					Delete(S);						   
			   		   WriteStatus('Frame= ' + Str(ind)+'   Total rejected pixels=  '+ Str(SumOf(not RJ)));
			   		BI:=BI*RJ;
			   			   Save(BI, Path+'B'+FHead+j+'.tif'); 				 	  	  {update masked images} 
					Delete(BI);
			   				Open(Path+'M'+FHead+j+'.tif',MI);					
			   		MI:=MI and RJ;		
			   			   Save(MI, Path+'M'+FHead+j+'.tif');							{update mask images}   
					Delete(MI);
					Delete(RJ);		
					Delete(j);
		end;
end;
{------------------------------------>end of procedure --Reject--}  	


{************************************************************************}
{ ***	                      Function GetFolder	   				  ***}									   
{************************************************************************}

  
function GetFolder( InitPath ) ;
{ Opens the folder browse dialog and returns a folder path }
{ InitPath may be used to set the initial folder (set to '' otherwise) }
{ Returns a null string if Cancel is pressed }

const
  BufSize = 255 ;
  
var
  FolderName ;

begin
  FolderName := InitPath + StringOfChar( ' ',BufSize - Length( InitPath ) ) ;
  if BrowseFolders( FolderName ) <> 0 then
    GetFolder := Trim( FolderName )
  else
    GetFolder := '' ;	
end;  { GetFolder }

{------------------------------------>end of Function --GetFolder--}  	

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

	
begin
		CR:=2.5;							//rejection criteria in sigma's
		F:=1;								//the first frame to process
		N:=100;								//the last frame to process
		Loops:=2;							//selection of number rejection runs
		FHead:='C_';							//file header
		TR:=60000;							//upper threshold for rejection
		xSize:=512;												  
		ySize:=512;
				
                if GetNumImages = 0 then halt( 'This program requires a desktop image' ) ;		
                   SelectImage( 'Select IMAGE for size definition',Img ) ;
                if IsNull( Img ) then halt( 'No image was selected' ) ;

				   Path:=GetFolder( InPath )+'\';
                if ( Path = '\' ) then halt( 'Directory not selected' ) ;
		WriteInfo(Path);					   
                          xSize:=GetXSize(Img);
                          ySize:=GetYSize(Img);
						  Free(Img);
						  GetString('Enter the FILE HEADER',FHead);
						  GetNumber('Enter the FIRST frame',F);						  
						  GetNumber('Enter the LAST frame',N);				   
						  GetNumber('Enter the NUMBER of LOOPS',Loops);
						  GetNumber('Enter the REJECTION CRITERIA',CR);						  						  
						  GetNumber('Enter the REJECTION TRESHOLD',TR);				
						  InitMASK;		   	   			  			{Initiate script}
		for k:=0 to Loops do
			begin
				ABI(k);				
				STD(k);
			if not (k=Loops	) then
				Reject(k);
			end;
			        PlaySound( 'C:\Program Files\Digital Optics\V++\Library\Mac.wav' ) ;	
					WriteInfo('Processing completed!!!');
end

