; Calculates average power spectra and defocus values for a set of micrographs
;
; PURPOSE:   Calculates average power spectra for a set of micrographs.
;            First converts a scanned file to SPIDER format.
;            Places power spectra in: power/roo****
;            Estimates defocus from 2D power spectra and places in: defocus
;            Uses SPIDER operation 'TF ED'
;            Inputs: Window size, percentage of the overlap
;                    distance of the window from the micrograph border
;                    PS mask radius
;
; SOURCE:  spider/docs/techs/recon/newprogs/powdefocus.spi 
;
; Edit following parameters and filenames as needed.
;
; [deci] = Decimation factor for 'DC S' operation
;      0 = Get value from params file (key=16)
;      1 = Full sized image
;      2 = 1/2 size
;      4 = 1/4 size
[deci]      = 0      ; Decimation factor (0 = Get value from params file)

[winsiz]    = 500    ; Window size of small pieces (Sx=Sy)
[xover]     = 50     ; X overlap % 
[yover]     = 50     ; Y overlap % 
[xd]        = 500    ; X Dist. from the edge 
[yd]        = 500    ; Y Dist. from the edge 
[radius]    = 225    ; Power spectrum mask radius (A)  (0 = don't mask)

[dodefocus] = 1      ; Find defocus values 
[keepspi]   = 1      ; Keep the on-disk temp spider file (0 = discard) 
[maskr]     = 0      ; Mask output power spectrum this far from border   (0 = nomask)  

; ----------- Input files --------------

[params]   = '../params'                     ; Parameter doc file

[sel_mic]  = '../sel_micrograph'             ; Micrograph selection doc file 

[micgr]    = '../Micrographs/raw{****[mic]}' ; Micrograph images

; ----------- Output files --------------

[outdir]   = 'power'                         ; Power spectra directory

[pow]      = '[outdir]/pw_avg{****[mic]}'    ; Power spectra images

[roo]      = '[outdir]/roo{****[mic]}'       ; Rotational averagedoc file

[ctf]      = '[outdir]/ctf{****[mic]}'       ; Spectrum, envelope, noise doc file

[out]      = 'defocus'                       ; Defocus values doc file

; ----------- Temp file --------------

[spi] = '_1'                                  ; Temp Spider file

; -------------- END BATCH HEADER --------------------------

MD
TR OFF                             ; Loop info turned off
MD
VB OFF                             ; File info turned off
MD
() OFF                             ; No need for () in DO loops          
MD
SET MP                             ; Use all available OMP processors
0

VM                                 ; Make sure output dir. present
mkdir -p [outdir]

;  -----  Get zip & format flags (can params vary??)
UD 1,[zflag]                       ; Get zip flag
[params]                           ; Params file           (input)
UD 2,[fflag]                       ; Get tif flag
[params]                           ; Params file           (input)
UD 3,[nx]                          ; HiScan  X parameters 
[params]                           ; Params file           (input)
UD 4,[ny]                          ; HiScan  Y dimension        
[params]                           ; Params file           (input)
UD 5,[sp_pixsiz]                   ; Get pixel size
[params]                           ; Params file           (input)
IF ([deci].LT.1) THEN
   UD 16,[deci]                    ; Get decimation factor  
   [params]                        ; Params file           (input)
   IF ([deci].LT.1) [deci] = 1.0   ; Should not be zero
ENDIF

;;[sp_pixsiz]=[sp_pixsiz]*[deci]     ; Adjust pixel size for decimation??
;;(I think we should NOT be multiplying by the pixel size.
;; Every other usage of key #16 in PARAMS ignores it. --trs)

IF ([radius] .GT. 0) THEN 
  [circf] = 2*[sp_pixsiz]/[radius] ; Mask radius (now in 1/px)
ENDIF

IF ([dodefocus] .GT. 0) THEN
  ; Want to determine defocus parameters for each group

  UD 7,[sp_sph_abb]                ; Spherical aberration (mm)
  [params]                         ; Params file           (input)
  UD 12,[sp_acr]                   ; Amplitude contrast ratio
  [params]                         ; Params file           (input)
  UD 14,[sp_lambda]                ; Lambda
  [params]                         ; Params file           (input)

  ; Put defocus parameters title in doc file
  DE
  [out]                            ; Defocus file           (removed)

  SD /    Micrograph     Defocus     Astig.Ang    Astig.Mag  Cutoff.Freq
  [out]                            ; Title for defocus file (output)
ENDIF

[tmparr] = 'tmparray'              ; Name for temporary incore doc file

[winsizd2]=[winsiz]/2
SD IC NEW                          ; Create an in-core doc file array
[tmparr]
4,[winsizd2]                       ; Creates a 4 column doc file

UD N [nummics]                     ; Get number of micrographs
[sel_mic]

DO [key]=1,[nummics]               ; Loop over all micrographs -------------------

  UD [key], [mic]                  ; Get micrograph number
  [sel_mic]                        ; Micrograph select doc file  (input)

   VM
   echo ' 'Loading Micrograph:   {*****[mic]}  

   ; Convert micrographs and load into incore file
   @loadmic([mic],[zflag],[fflag],[deci],[nx],[ny],[keepspi])
   [micgr]                         ; Micrograph template     (input)
   [spi]                           ; Spider file             (output)
   _4                              ; Hiscan & Nikon scratch  (output)

   VM
   echo " Creating spectra:   [pow] and [roo]"  
   DE                              ; If output doc file exists, delete 
   [roo]                           ;                      (removed)

   FI H [x],[y]                    ; Get dimensions of converted image
   [spi]                           ; Micrograph image     (input)
   NSAM,NROW

   ; [winsiz] = large window size

   ;[v87]=(([x]-2*[xd])*([y]-2*[yd]))number of pieces in this micrograph

   [nox]=100/(100-[xover]); Normalization of % of the overlap in x 
   [noy]=100/(100-[yover]); Normalization of % of the overlap in y

   [npx]=INT([nox]*(([x]-2*[xd])/[winsiz]-1)+1)  ; No. of pieces horizontal dim.(X)
   [npy]=INT([noy]*(([y]-2*[yd])/[winsiz]-1)+1)  ; No. of pieces vertical dim.(Y)

   [v94] = 0
   ; -------------------- Loops over X and Y -------------- 
   DO [v12] = 1, [npy]
      [uly] = ([winsiz]/[noy])*([v12]-1) + [yd]     ; Y-direction
      DO  [v13] = 1, [npx]
         [ulx] = ([winsiz]/[nox])*([v13]-1) + [xd]  ; X-direction

         WI                       ; Window image
         [spi]                    ; Get small pieces of input image 
         _5
         [winsiz],[winsiz]        ; Size of small pieces
         [ulx],[uly]              ; Coord. of upper left corner

         RA                       ; Correct ramp effects 
         _5                       ;                  (input)
         _6                       ;                  (output)

         PW                       ; Calc power spectrum of each small piece
         _6                       ;                  (input)
         _7                       ; PS               (output)

         SQ                       ; Square image             
         _7                       ;                  (input) 
         _9                       ; Squared          (output)

         [v94] = [v94]+1

         IF ([v94].ge.2) THEN
            AD                    ; Add images
            _3                    ;                  (input)
            _9                    ;                  (input)
            _3                    ; Sum              (output)
            *
         ELSE
            CP                    ; Copy image
            _9                    ;                  (input)
            _3                    ;                  (output)
         ENDIF
      ENDDO
   ENDDO

   ; --------- End Loops. Resulting spectrum in:  _3 ------------

   AR                            ; Average over power spectra of small pieces 
   _3                            ; Inline file              (input)
   _8                            ; Inline file              (output)
   P1/[v94]

   WU                            ; Calculate the square root 
   _8                            ; Inline file              (input) 
   _2                            ; Inline file              (output)  

   IF ([circf].GT.0.0) THEN      ; Mask out center 
      FI H [nsam],[nrow]
      _2                         ; Inline file              (input)
      NSAM,NROW                  ; Get image dimensions

      [cx] = [nsam]/2 + 1        ; x center
      [cy] = [nrow]/2 + 1        ; Y center

      GP [pv]                    ; Get pixel intensity
      _2                         ; Inline file              (input)
      ([cx]+10,5)                ; Location for background   

       
      [circ] =[circf]*[nsam]     ; Mask radius
      [circo]=[maskr]            ; Outer Mask radius
      if ([maskr].GT. 0) THEN
        [circo]=[cx]-[maskr]     ; Outer Mask radius
      ENDIF

      MA                         ; Mask image
      _2                         ; Inline file              (input)
      _9                         ; Inline file              (output)
      [circo],[circ]             ; Radii of central dot. e.g. (2*5A/px/200A)/x 500 = 25 pixels
      D                          ; Disk
      E                          ; External background
      [pv]                       ; Background
      [cx],[cy]                  ; Center coordinates

      VM
      echo " Masked spectra:     {*****[circo]} and {***[circ]}"  

      CP                         ; Copy image
      _9                         ; Inline file              (input)
      [pow]                      ; Disk power spectrum      (output)
   ELSE
      CP                         ; Copy to disk
      _2                         ; Inline spectrum file     (input)
      [pow]                      ; Disk power spectrum      (output)
   ENDIF                         
 
   RO                            ; Rotational average
   _2                            ; File                     (input)
   _3                            ; 1D rotational average    (output)

   LI D                          ; List in doc file
   _3                            ; 1D rotational average    (input)
   [tmparr]                      ; Doc file                 (output)
   R                             ; Row wanted
   (1)                           ; Which row

   ; Labels for columns in output doc file
   SD /      Amplitude    X-axis     Spat.Freq.
   [roo]                                ; Doc file    (output)

   [winsizd2t] = [winsiz]/2             ; 1/2 large window size

   DO [v70] = 1,[winsizd2]
      UD IC [v70], [v72],[v73]
      [tmparr]                          ; Doc file    (input)

      ; Compute spatial freq.
      [sfreq] = [v70] / (2.0 * [winsizd2t]) 

      SD [v70],[v72],[v73],[sfreq]
      [roo]                             ; PS Doc file (output)
   ENDDO

   ; ----------- Get defocus value  -------------------

   IF ([dodefocus] .GT. 0) THEN
     ; Save defocus parameters for each group in doc file

     DE                                ; Delete group doc file
     [ctf]                             ;                    (removed)

     ; Transfer Function -- Estimation of CTF parameters

     TF ED [v12],[v13],[v14],[v15],[v16] 
     [pow]                             ; 2D spectrum        (input)
     ([sp_pixsiz], [sp_sph_abb])       ; Pixel size, Spherical aberration
     ([sp_lambda])                     ; Lambda
     ([sp_acr])                        ; Ampl. contrast ratio
     [ctf]                             ; Doc file           (output)

     ; Save parameters for this group in summary doc file
     SD [key], [mic],[v14],[v12],[v13],[v16]
     [out]                             ; Doc file           (output)

     VM
     echo ' 'Micrograph: {*****[mic]}'   'Defocus: {%f8.2%[v14]};echo ' '

   ENDIF
ENDDO

SD E
[out]                                  ; Doc file           (output)

VM
echo ' ' ; cat [out].$DATEXT ; echo ' '
 
EN
;