; NOTE : CHANGED 'PK 3D' to 'PK 3R' READ COMMENTS BELOW 'PK 3R' OPERATION ; : ADDED MASK TO RESTRICT THE PEAK SEARCH INSIDE THE MASK ONLY, LOOK ABOVE 'PK 3R' ; ; Searches for molecular signature
; sigsloop.pam                                Bimal Rath   :  JAN 2003
;                        PARALLELIZED BY      ArDean Leith
;
; SEARCHES FOR MOLECULAR SIGNATURE (REFERENCE VOLUME) INSIDE A LARGE VOLUME. 
; BE SURE BOTH VOLUMES HAVE SAME MAGNIFICATION (1 PIXEL = "N" NANOMETER) 
;
;;MD
;;VB OFF
MD
TR OFF



; get around the named variable related bug until Dean fixes it
;x72 = [Y72]
;x72

x72   ; Needed from startup line!!!!


; READ INPUT FILES AND REGISTERS
@sigs_settings[x12,x65,x66,x67,x74,x75,x76,x77,x78,x79,x83,x95,x39,x61,x64]

IF (x12 .GT. 1) THEN
   MD
   SET MP                 ; Set number of OMP processors
   x12
ENDIF

; FIND NSAM, NROW AND NSLICE OF LARGE AND REFERENCE VOLUMES
FI x20,x21,x22            ; Large volume size
[LARGE_VOLUME]
12,2,1

x27 = INT(x20/2)+1        ; Center of large volume
x28 = INT(x21/2)+1
x29 = INT(x22/2)+1
                        
FI x23,x24,x25            ; Small volume size
[PADDED_REF_VOLUME]
12,2,1

x26 = x23*x24*x25         ; Number of pixels in small volume

x33 = INT(x23/2)+1        ; Center of small volume
x34 = INT(x24/2)+1
x35 = INT(x25/2)+1 

x86 = x20-x23+1           ; big-small image size
x87 = x21-x24+1
x88 = x22-x25+1

x36 = INT((x20-x23)/2)+1  ; Corner of small image inserted into center of big
x37 = INT((x21-x24)/2)+1
x38 = INT((x22-x25)/2)+1  
 
WI           ; Window 
[MASK_PKR]   ; Peak restricting mask
_50          ; Output volume
x86,x87,x88  ; Difference size
x33,x34,x35  ; 1/2 size of padded reference volume

x97 = 0                   ; Doc. file output line number

IF (x61 .NE. 1) THEN 
   ; NONTOMOGRAPHIC VOLUME

   ; MAKE AN APPROPRIATE MASK FROM THE REFERENCE VOLUME   
   ; THRESHOLD VALUE NEEDS TO BE CHANGED FOR EACH REFERENCE VOLUME
   ; NOTE: DON'T REUSE _1

   TH M           ; Threshold to create mask
   [PADDED_REF_VOLUME]
   _1             ; Output mask volume
   B              ; Voxels below x83 set = 0
   x83            ; Threshold level
ENDIF

; EULER ANGLE SEARCH IS DONE HERE

; USED IN LOOP
x56 = x65 + (x72-1)*x74          
        
DO LB5 x71 = 1, x78
   x57 = x66+(x71-1)*x75
      
   DO LB6 x70 = 1, x77  
      x58 = x67+(x70-1)*x76         

      ; DO LOOP NUMBER
      x90 = (x71-1)*x77 + x70

      IF (x61 .EQ. 1) THEN
         ; TOMO VOLUME
	 ; ROTATE THE PADDED REFERENCE VOLUME
         RT 3D
         [PADDED_REF_VOLUME]
         _55
         x56,x57,x58

         x69 = x33 - 2 ; Radius of the object minus 2 pixels

         ; CREATE PROJECTIONS FROM THE ROTATED PADDED REFERENCE VOLUME
         PJ 3Q	       ; Create projection series
         _55	       ; Input vol
         x69           ; Radius of the object
         [SEL_DOC]     ; EDIT desired selection file if necessary
         [ANG_DOC]     ; EDIT desired angles file if necessary
         _56@***       ; Inline stack template
         x64           ; # of images to be stored in the inline stack        
   
         x68 = 0.3     ; Frequency cut-off for Parzen Filter
         x41 = 1       ; Begining of the Y-axis for reconstruction
   
         ; DO THE RECONSTRUCTION ON THE FLY OF THE OBJECT FROM THE PROJECTIONS
         BP W2
         _56@***       ; EDIT file prefix of aligned images if necessary
         [SEL_DOC]     ; EDIT desired selection file if necessary
         [ANG_DOC]     ; EDIT desired angles file if necessary
         x69,x23       ; EDIT (image height/2 then minus 2), (desired recon. depth)
         x41,x23       ; EDIT beginning and ending y-axis slices (max is image height)
         x68           ; Frequency cut-off for Parzen Filter
         _57           ; EDIT output filename if desired
   
      ENDIF

      IF (x39 .EQ. 1) THEN
         ; ASYMMETRIC MASK IS USED

         IF (x61 .EQ. 1) THEN
            ; TOMOGRAPHIC VOLUME
	    
	    FS x80,x81           ; Find max and min of pixel values
	    _57                  ; Reconstructed refrence volume
	    
	    x82 = (x81 - x80)/2   
	    x82 = b82 + x80      ; Thershold is halfway between max and min of pixel values

            ; THRESHOLD THE RECONSTRUCTED VOLUME
            ; NOTE: DON'T REUSE _4
            TH M           ; Threshold to create mask
            _57            ; Reconstructed refrence volume
            _4             ; Thresholded rotated mask output
            B              ; Below = 0
            x82            ; Threshold value

         ELSE
            ; NON-TOMOGRAPHIC VOLUME

            ; ROTATE THE MASK
            RT 3D          ; Rotate mask volume
            _1             ; Mask voume input
            _2             ; Rotated mask output
            x56,x57,x58    ; Rotation angle

            ; THRESHOLD THE ROTATED MASK
            ; NOTE: DON'T REUSE _4
            TH M           ; Threshold to create mask
            _2             ; Rotated mask input
            _4             ; Thresholded rotated mask output
            B              ; Below = 0
            (.5)           ; Threshold value

         ENDIF

      ELSE
         ; ROTATIONALLY INVARIANT MASK IS USED

         IF (x90 .EQ. 1) THEN
   
            x53 = (x23/2) - 2   ; Radius of the mask sphere
            x42 = 0
            ; JUST DO IT ONCE. (FIRST LOOP ONLY)
            MO 3
            _4
            x23,x24,x25
            SP
            N
            (1.0)
            x53,x42
            x33,x34,x35
            (0,0,0)
         ENDIF

      ENDIF

      ; IN ASYMMETRIC CASE SET X55 =1 FOR ALL LOOPS
      ; IN SYMMETRIC  CASE SET X55 =1 ONLY FOR THE FIRST LOOP
      IF (x39 .EQ. 1) THEN
         x55 = 1
      ELSE 
         x55 =0
      ENDIF

      IF (x90 .EQ. 1) THEN
         x55 = 1
      ENDIF

      IF ( x55 .EQ. 1) THEN
         ; FIND NUMBER OF NON-ZERO PIXELS INSIDE MASK
         ; (SAME AS SUM OF ALL PIXELS IN MASK)
         FS x11,x11,x50,x11
         _4             ; Thresholded rotated mask

         x50 = x50*x26
         IF (x50.LE.0)THEN
            VM
            echo 'NO NON-ZERO PIXELS INSIDE MASK'
            EN
         ENDIF
         x51 = 1.0 / x50 ; For speed

         ; CREATE A BLANK VOLUME, SAME SIZE AS LARGE VOLUME
         BL             ; Create blank volume
         _25            ; Blank output volume
         x20,x21,x22    ; Size
         N              ; Not average background
         (0)            ; Background

         ; INSERT ROTATED MASK INSIDE THE BLANK VOLUME
         IN              ; Insert
         _4              ; Thresholded rotated mask
         _25             ; Padded mask output volume
         (1,1,1)         ; Position

         ; FT ON PADDED MASK,  NOTE: DON'T REUSE _36
         FT              ; Fourier transform
         _25             ; Padded mask
         _36             ; Fourier of padded mask
	 
;bimal testing begin
cp
_25
_90
	 
CP
_36
_91
;bimal testing end

         ; MULTIPLY FT OF LARGE VOLUME WITH COMPLEX CONJUGATE OF PADDED MASK
         MU M            ; Complex multiplication
         [LARGE_FT]      ; First input
         _36             ; Fourier of padded mask
         _35             ; Output
         *


;bimal testing begin
cp
[LARGE_FT]
_92
	 
CP
_35
_93
;bimal testing end


         ; INVERSE FT
         FT              ; Inverse Fourier transform
         _35             ; Input
         _27             ; Output


;bimal testing begin
cp
_27
_94
;bimal testing end

         x52=x51*x51

         ; NORMALIZE 
         AR              ; Arithmetic operation
         _27             ; Input
         _27             ; Output
         P1*P1*x52       ; P2 = (P1 / (No. OF NON-ZERO PIXELS INSIDE MASK))**2 

         ; MULTIPLY FT OF SQUARE OF LARGE VOLUME WITH COMPLEX CONJUGATE 
         ; OF FT OF BLANK VOLUME
         MU M           ; Complex multiplication
         [LARGE_SQ_FT]  ; Input
         _36            ; Multiplier
         _35            ; Output
         *

;bimal testing begin
cp
[LARGE_SQ_FT]
_95
	 
CP
_35
_96
;bimal testing end

	 
         ; DO INVERSE FT
         FT           ; Inverse Fourier transform
         _35          ; Input
         _25          ; Output

;bimal testing begin
cp
_25
 _97
;bimal testing end


         ; NORMALIZE
         AD F
         _25          ; Input (From FT OF SQUARE OF LARGE VOLUME)
         _27          ; Input (From FT OF LARGE VOLUME)
         x51,-1.0     ; p_25 = p_25 * x51 - p_27
         _25          ; Output
          
         ; GET RID OF SQRT OF -VE # AND DIVISION BY ZERO(WHILE DIVIDING THE 
         ; CCF BY LOCAL STANDARD DEVIATION)
	 ; IF YOU FIND CCC > 1.0 AND THE LOCATIONS OF THE MOTIF OUTSIDE THE OUTLINE OF THE SEARCHED
	 ; LARGE VOLUME THEN LOOK AT THE FILE _25 AND YOU WILL FIND PIXEL VALUES
	 ; QUITE SMALL ~ < 1E-10 BUT NOT EQUAL TO ZERO. IN THIS CASE, CHANGE THE MASK THRESHOLD
	 ; IN THE FOLLOWING OPERATION FROM ZERO (0) TO SOMETHING LIKE 1E-10. THIS MAY SOLVE
	 ; THE PROBLEM. HOWEVER, YOU MAY NEED TO CHANGE THIS PARAMETER TO FIND THE CORRECT ONE
	 ; TO BE USED. THIS PROBLEM OF GETTING -VE VARIANCE OR GETTING INCORRECT VARIANCE WHEN
	 ; PIXELS UNDER THE MASK HAVE SAME/VERY_CLOSE VALUES IS DUE TO THE WAY VARIANCE IS 
	 ; CALCULATED USING FOURIER TRANSFORM.	 
         TH M
         _25
         _80
         B
         (0)
        
         MM
         _80
         _25
         (9E+20)

         ; LOCAL STANDARD DEVIATION
         WU           ; Square root
         _25          ; Input
         _25          ; Output
           
         ; NOTE: DON'T REUSE _15  
         WI           ; Window
         _25          ; Input
         _15          ; Output                       
         x86,x87,x88  ; Size of difference
         (1,1,1)      ; Position

      ENDIF

      ; ROTATE THE REFERENCE , DON'T REUSE _2 ------------------ step 2
      
      IF (x61 .NE. 1) THEN 
         ; NON-TOMOGRAPHIC VOLUME
         
         RT 3D               ; Rotate volume
         [PADDED_REF_VOLUME] ; Input
         _2                  ; Output
         x56,x57,x58         ; Angle
      ELSE
         ; TOMOGRAPHIC VOLUME
         CP
         _57
         _2
      ENDIF

        
      ; PREPARE THE REFERENCE VOLUME SUCH THAT THE AVERAGE INSIDE 
      ; THE MASK = 0 AND THE STANDARD DEVIATION INSIDE THE MASK = 1
      MM           ; Mask multiplication
      _4           ; Rotated mask input
      _2           ; Masked rotated reference input/output volume
      (0)          ; Background for voxels < 0.5

      ; FIND AVERAGE of rotated, masked ref. volume
      FS x11,x11,x62,x63
      _2           ; Masked volume
        
      ; SUM OF ALL PIXELS IN ROTATED, masked ref. volume
      x40 = x62*x26 ; Average * No. of pixels

      SQ           ; Square the rotated, masked volume
      _2           ; Masked volume
      _7           ; Squared masked volume
               
      ; FIND AVERAGE OF SQUARED ROTATED MASK
      FS x11,x11,x62,x11
      _7           ; Squared rotated masked volume

      ; SUM OF ALL PIXELS SQUARED IN ROTATED MASK.
      x45 = x62*x26

      ; SD INSIDE MASK        
      x46 = SQRT( (x45 - (x40*x40*x51)) / (x50-1) )

      ; AVG INSIDE MASK         
      x47  =  x40/x50    
      x48  = 1.0 / x46  ; For speed

      ; NORMALIZE
      AR          ; Arithmetic operation
      _2          ; Masked volume
      _7          ; Normalized masked file
      (P1-x47) * x48

      MM          ; Mask multiply operation
      _4          ; Mask 
      _7          ; Masked input/output file
      (0)

      ; CREATE AN EMPTY VOLUME OF DIMENSION = LARGE VOLUME 
      ; PASTE THE PREPARED REFERNCE VOLUME AT THE MIDDLE OF 
      ; THIS EMPTY VOLUME            
      BL          ; Create blank volume
      _25         ; Large blank volume
      x20,x21,x22 ; Size
      N           ; Not average background
      (0)         ; Background

      IN          ; Insert
      _7          ; Small volume
      _25         ; Into large blank volume
      1,1,1       ; put small image inside big at the corner

      ; FIND CROSS CORRELATION FUNCTION OF THE ABOVE VOLUME WITH
      ; WITH THE LARGE VOLUME.                           

      FT
      _25
      _60
      
;bimal testing begin
cp
_25
_98
	 
CP
_60
_99
;bimal testing end      
      
      

      ; SET F(0,0) ELEMENT = ZERO. DONE TO DO SIMILAR NORMALIZATION
      ; AS DONE IN REAL SPACE 
      RP
      _60
      (1,1,1)
      (0)

      RP
      _60 
      (2,1,1)
      (0)

      FT
      [LARGE_VOLUME]
      _61


;bimal testing begin
cp
[LARGE_VOLUME]
_81
	 
CP
_61
_82
;bimal testing end      
      

      ; DON'T CHANGE ORDER OF INPUT IN THE FOLLOWING OPERATION
      MU M
      _61
      _60
      _62
      *

;bimal testing begin
cp
_60
_83
	 
CP
_62
_84
;bimal testing end      
      

      ; DO INVERSE FT
      FT
      _62
      _27


;bimal testing begin
cp
_27
_85
;bimal testing end      
      


      WI           ; Window
      _27          ; CC volume
      _13          ; Output
      x86,x87,x88  ; Difference size
      (1,1,1)      ; left corner of small volume    

      ; DIVIDE THE CC FUNCTION WITH TOTAL NUMBER OF NON-ZERO PIXELS
      ; INSIDE THE MASK         
      AR          ; Arithmetic operation
      _13         ; Input
      _17         ; Output
      P1 * x51

      ; DIVIDE THE  ABOVE RESULT WITH CORRESPONDING ELEMENT OF 
      ; THE LOCAL STANDARD DEVIATION ARRAY 
      MU D        ; Divide
      _17         ; Input
      _15         ; Input
      _13         ; Output
      *

      DE                         ; This doc file created by each 'PK 3D'
      [DOC_FILE_DEL]_{****x72}

      ; RESTRICT THE PEAK SEARCH INSIDE A GIVEN MASK
      MM
      _50
      _13
      (0)


      ; PEAK SEARCH                 
      PK 3D                      ; Peak search
      _13                        ; Input file
      +                          ; Maxima
      (x95,1)                    ; Number of peaks, origin override flag
      N                          ; No center of gravity calc.
      (1,1)                      ; xy co-ordinate
      (1)                        ; z co-ordinate
      (1)                        ; peak no. for ratio
      N                          ; No box
      [DOC_FILE_DEL]_{****x72}   ; Output doc file


      ; WRITE EULER ANGLES, XYZ POSITIONS, & PEAK HEIGHT TO FINAL DOC FILE 
      DO LB10 x96 = 1,x95
         UD S,x96,x11,x12,x13,x14,x16,x17,x18
         [DOC_FILE_DEL]_{****x72}
         ; CORRECT FOR THE CENTER OF THE PEAK WITH RESPECT TO LARGE VOLUME.
         ; THE PEAK HEIGHT DETERMINED IN PEAK SEARCH STEP IS WITH RESPECT TO THE
         ; VOLUME CREATED BY SUBTRACTING THE DIMENSION OF REFERENCE VOLUME
         ; FROM THE LARGE VOLUME. FACTOR OF NSAM/2+1, NROW/2+1 AND NSLICE/2+1
         ; ARE ADDDED
         x30 = x11+x33
         x31 = x12+x34
         x32 = x13+x35
         x97 = x97 + 1


         ; DO NOT WRITE CROSS-CORRELATION COEFFICIENTS (CCC) IN SCIENTIFIC FORMAT. UNIX 
         ; "SORT" COMMAND WILL HAVE TROUBLE TO SORT THE FILE.
         ; CCC LESS THAN 0.05 DOES NOT MEAN MUCH
         IF (x18 .LT. (.05)) x18 = -99
	 
; bimal testing begin
IF (x18 .GT. (1.0)) then
   cp
   _81
   file_81
	 
   cp
   _82
   file_82
	 
   cp
   _83
   file_83
	 
   cp
   _84
   file_84
	 
   cp
   _85
   file_85
	 
   cp
   _90
   file_90
	 
   cp
   _91
   file_91
	 
   cp
   _92
   file_92
	 
   cp
   _93
   file_93
	 
   cp
   _94
   file_94
	 
   cp
   _95
   file_95
	 
   cp
   _96
   file_96
	 
   cp
   _97
   file_97
	 
   cp
   _98
   file_98
	 	 
	 
   cp
   _99
   file_99
	 	 	 	 	 
	 
   vm
   echo "\n\n*************ERROR ***************\n\n"
	    
endif	 
; bimal testing end
	 
        

         SD x97,x56,x57,x58,x30,x31,x32,x18
         [DOC_FILE_OUT]_{****x72}       
      LB10
         
      UD E                                                               
      DE                         ; This doc file created by each 'PK 3D'
      [DOC_FILE_DEL]_{****x72}
   LB6
      
   MY FL
      
LB5

SD E
[DOC_FILE_OUT]_{****x72}       

@signal_pub(x72)             ; Signal finished          
EN

;