;
;  Example of the batch program implementing reference-free alignment
;  of the series of images using AP SA and AP RA commands. The program
;  uses following strategy: initially all the particles are centered
;  approximately using shift alignment (AP SA) with 180 degrees
;  rotation check.  This is followed by a number (up to 9 in this version)
;  steps of shift and rotational alignment of the input data series.
;  After each such step the resulting average is centered using
;  phase approximation of the center of gravity command (CG PH).
;  This program works for consecutively numbered images only, but the 
;  numbering does not have to start with one.
;

;  b02.sup  -  ALIGNMENT  01/09/98
;  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
;
;  Important notes:
;
;  (1) To start the program one has to modify the input file names
;  and specify input parameters: number of steps required, number
;  of particles and ring numbers used by APRA (determined by the
;  expected radius of the particle).
;  
;  (2) The shift alignment step does not enforce centering of the average
;  by itself.  The average produced by AP SA (AVS*** in the following
;  program) has to be centered explicitly.  In the program command
;  CG PH (phase approximation of the center of gravity) is use for
;  this purpose.  It is not guaranteed to work for all the possible
;  shapes and it is advisable to verify whether the average was
;  shifted correctly.  If not, the command can be easily replaced
;  by the other mean of centration (for example croscorrelation
;  with low-pass filtered disk).
;
;  (3) Please note that CG PH will NOT work correctly for images
;  with inverted contrast (dark particle on bright background).
;  In such case either contrast of the particles has to be inverted,
;  or the contrast of the average before application of CG PH.
;
;  (4) VERY IMPORTANT:  if you want to start the program again
;  from the beginning you have to either delete previously produced
;  document files or change the names of the existing file.
;
;  (5) The program can be easily restarted to make additional steps
;  of alignment. If initially program was ran for 5 steps of alignment
;  than it is enough to change register X61 in the input block
;  to the new iterations limit, remove part of the program between "End
;  of the input parameters" block to the main do-loop LB9 and change
;  in this line "2" to "6" (in this case, in general to the next step
;  number).
;
;  (6) The program makes it possible to calculate the current resolution
;  of 2D average. The code has to be added by the user in the place marked.
;
;
; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;
;  DESCRIPTION OF THE FILE NAMES USED:
;
;  The input data is the series of images
;  data/PAR****    1-1210
;
;  The output data is the separate series of images calculated after
;  each step of alignment
;  ali/PUR*****
;  first digit stands for the number of alignment step, the remaining
;  four digits are reserved for the particle number.
;
;  Program uses the following document files:
;
;  FALS - dummy file to initialize first shift alignment containing
;         particle numbers
;
;  APSA*** - ending with the number indicating the number of alignment
;            step. This is the document file produced by the shift
;            alignment (AP SA). Details are in the manual chapter APSA.MAN.
;
;  APRA*** - ending with the number indicating the number of alignment
;            step. This is the document file produced by the rotational
;            alignment (AP RA). Details are in the manual chapter APRA.MAN.
;
;  PRM***  - ending with the number indicating the number of alignment
;            step. This is the document file containing the combined alignment
;            parameters for each particle after all alignment steps made
;            up to this one. When applied to the initial (input) data
;            they will produce the aligned and centered series of particles.
;            The angles (after reversing the sing) can be used in 3D
;            reconstruction.
;            This file contains four columns:
;              1                      2        3        4
;            particle number        angle    X shift   Y shift
;
;  SHIFTS -  contains X,Y shifts needed to center the average calculated
;            after each shift alignment.
;
;
;  The program requires specification of the following input: 
;  X61 - number of alignment steps 
;        First is the initialization step - shift only (APSA)
;        Next X61-1 double steps are made - (rotation + shift)
;
;  X62 - number of particles
;  X63 - first particle number
;  X41 - first ring number (required by AP RA)
;  X42 - last ring number (required by AP RA)
;
;
; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;
;
;  Beginning of the program.
;  -------------------------; 
;                                                            
; Specify input parameters:
;-------------------------------------------------------------------; 
;  X61 - number of alignment steps 
;        First is the initialization step - shift only (APSA)
;        Next X61-1 double steps are made - (rotation + shift)
X61=5
;  X62 - number of particles
X62=2530
;  X63 - first particle number
X63=1
;
;  X41,X42 - specification of first and last ring numbers
X41=5
X42=29
;
; End of the input parameters block.
;-------------------------------------------------------------------; 
;
;  X71 - current step number
X71=1
; FALS is a dummy document file containing particle numbers in the first
;      column. In general they do not have to be in consecutive order
;      neither begin from one.
X31=0.0
X32=1.0
DO LB1 I=1,X62
SD X0,X0,X31,X31,X32
RESULTs/FALS
LB1
SD E
RESULTs/FALS
;
; First shift alignment with 180 degs check to center all the particles.
AP SA
data/PAR001
RESULTs/FALS
1
Y
RESULTs/APSA{***X71}
RESULTs/AVS{***X71}
;
;
;   Calculate center of gravity using CG PH to center the average.
;
CG PH,X21,X22
RESULTs/AVS{***X71}
SD X71,X21,X22
SHIFTS
;
;  X21 and X22 contain integer shift values which will be used to
;  center all the images.
;
;  
;  Shift and rotate all the input data (using center correction
;  shifts found by CG PH).  The resulting file series will be used
;  in the next alignment step.
;
DO LB2 I=1,X62
;               For short image series use UD IC instead on remove ";" signs
;               from UD ICE lines in the program.
UD S,X0,X31,X32,X33,X34,X35,X36
RESULTs/APSA{***X71}
X34=X34-X21
X35=X35-X22
X91=X71*10000+X31
RT SQ
data/PAR{***X31}       
ali/PUR{*****X91}
X32
X34,X35
SD X0,X31,X32,X34,X35
RESULTs/PRM{***X71}
LB2
SD E
RESULTs/PRM{***X71}
;
;UD ICE
;RESULTs/APSAX71
;
;
;  Main loop LB9 - make X61-1 steps of alignment.
;  ---------------------------------------------; 
;
DO LB9 I=2,X61
;
X71=X0
X70=X71-1
X93=X70*10000+X63
X92=X70*10000+X62
;
;                     Rotational alignment
AP RA
ali/PUR*****
X93-X92
X41,X42
(1)
F
RESULTs/APRA{***X71}
;
;                     Shift alignment
AP SA
ali/PUR{*****X93}
RESULTs/APRA{***X71}
(1)
N
RESULTs/APSA{***X71}
RESULTs/AVS{***X71}
;
;   Calculate center of gravity using CG PH to center the average.
;
CG PH,X82,X84,X21,X22
RESULTs/AVS{***X71}
SD X71,X21,X22
SHIFTS
;  
;  At this stage it is advisable to remove the aligned image series
;  aligned during the provious step for the sake of the disk space.
;
;DE A
;ali/PUR00X93
;
;
;  Rewrite the whole image series combining alignment parameters
;  and applying them to the input series PAR****.  A new, aligned
;  and centered series PUR***** is produced.
;
DO LB8 I=1,X62
;               For short image series use UD IC instead on remove ";" signs
;               from UD ICE lines in the program.
UD S,X0,X81,X82,X84,X85
RESULTs/PRM{***X70}
UD S,X0,X31,X32,X33,X34,X35,X36
RESULTs/APSA{***X71}
SA P,X82,X84,X85,X32,X34,X35,X52,X54,X55
X54=X54-X21
X55=X55-X22
X91=X71*10000+X81
RT SQ
data/PAR{***X81}     
ali/PUR{*****X91}
X52
X54,X55
SD X0,X81,X52,X54,X55
RESULTs/PRM{***X71}
LB8
;UD ICE
;RESULTs/PRMX70
;
;
;   The following block can be used to calculate current averages
;   and estimate the resolution.
;X93=X71*10000+X63
;X92=X71*10000+X62
;AS R
;ali/PUR00X93
;X93-X92
;E
;AASX71
;AAVX71
;ABSX71
;ABVX71
;
;  Calculate the resolution using AASX71 and ABSX71.
;  Use your own procedure ..... 
;
LB9
EN
;
;  End of the program.
;  -------------------;