([pixsiz],[iter],[ampenhance],[r2],[sphdecon],[gold-std],[qsub]) ; ; Merge groups, finds FSC, filters volumes
 ;
 ; SOURCE: spider/docs/techs/recon1/Procs/merge-fsc-filt.spi
 ;         New                                ArDean Leith  Nov 2000
 ;         Amplitude enhancement              ArDean Leith  Apr 2005
 ;         CofG centering now done here       ArDean Leith  Sep 2010
 ;         FSC, pixsiz                        ArDean Leith  Aug 2012
 ;         Defocus group removal              ArDean Leith  Oct 2013
 ;         For CTF corrected images           ArDean Leith  Oct 2013
 ;         For gold standard reconstruction   ArDean Leith  May 2014
 ;         Rewritten                          ArDean Leith  Mar 2016
 ;         Added FSC mask creation here       ArDean Leith  Apr 2016
 ;
 ; PURPOSE: Merge ('AS S')  group volumes into a single set of volumes, 
 ;          find resolution ('FSC'),
 ;          filter ('FD C') volumes to limit resolution, and 
 ;          prepare next reference volumes.    (This is not done in parallel).
 ;
 ; CALLED FROM: pub-refine pub-refine.spi and
 ;              refine     refine.spi
 ;
 ; INPUT REGISTERS:
 ;   [pixsiz]             Pixel size        
 ;   [iter]               Refinement step iteration counter  
 ;   [ampenhance]         Use amplitude enhancement         (>0 == True) 
 ;   [r2]                 Radius of object 
 ;   [sphdecon]           Spherical deconvolution sigma     (0 == no deconvolution) 
 ;   [gold-std]           Gold standard reconstruction      (1 == gold-standard workflow)
 ;   [qsub]               Parallel deconvolution  (<0 == serial, 0 == PubSub, 1 == PBS) 
 ;
 ; FILES: ([win] denotes input directory, [out] denotes output directory 
 ;               '***' denotes group,     '%'   denotes subset,
 ;               '##'  denotes iteration, '##+' denotes next iteration)
 ;
 ; INPUT FILES:
 ;   [sel_group]          [win]/sel_group         List of groups                  (one)          
 ;   [next_group_vol]_s@  work/vol_##+_***_s@     Group subset volume             (one/group)
 ;
 ; OUTPUT FILES: ('[out] denotes output directory)
 ;   [next_vol]_s@_raw    [out]/vol##+_s@_raw     Unfiltered subset vol                (two)
 ;   [next_vol]_s@_unfilt [out]/vol##+_s@_unfilt  Unfiltered (Deconvolved?) subset vol (two) 
 ;   [next_vol]_s@        [out]/vol##+_s@         Filtered   (Deconvolved?) subset vol (two)
 ;
 ;   [next_vol]_unfilt    [out]/vol##+_unfilt     Unfiltered overall vol          (one)
 ;   [next_vol]           [out]/vol##+_all        Filtered   overall vol          (one)
 ;   [next_vol]]_cent     [out]/vol##+            Centered filtered  overall  vol (one)
 ;   [fsc_mask]           [out]/fsc_mask          Soft FSC mask volume            (one)
 ;   [next_u_fsc]         [out]/fscdoc_u_##+      Unmasked FSC curve doc file     (one)
 ;   [next_m_fsc]         [out]/fscdoc_m_##+      Masked   FSC curve doc file     (one)
 ;   [enhance_doc]        work/enhance_doc_##+    OPTIONAL Amplitude filter file  (one)
 ;
 ;-------------------------- END BATCH HEADER ------------------------------------------------

 ; All groups must be present here to use: 'AS S'!

 [next-iter] = [iter] + 1        ; Next iteration

 SYS
   echo ; echo "  Merging group volumes for iteration: {%I0%[iter]}"

 DO [s] = 1,2                    ; Loop over two volume subsets ------------- _raw

   ; Combine all subset group volumes into two subset volumes
   AS S                          ; Average group volumes for this subset
     [next_group_vol_stem]       ; Group volume template     (input)
     [sel_group]                 ; Group selection file
     A                           ; Find average
     [next_vol_s]_raw            ; Combined subset volume    (output)

 ENDDO                           ; End of loop over volume subsets ----------

 ; Create: [next_vol_s]_unfilt using deconvolution if wanted ---------------- _unfilt
 IF ( [sphdecon] > 0 ) THEN
   ; Apply Spherical deconvolution to combined subset volumes
   ; Puts de-convolved subset volumes into: [next_vol_s]_unfilt

   IF ( [qsub] >= 0) THEN
     ; Parallel deconvolution
     [task]     = 3              ; Task 3 --> pub-refine-start starts: sphdecon
     [num-jobs] = 2              ; Number of parallel groups (Here == 2 subsets)
     [script]   = './spider $PRJEXT/$DATEXT @pub-refine-start'
     @pub-submit([iter],[num-jobs],[task],[qsub])
        [script]
   ELSE
     ; Serial deconvolution
     [s] = 1
     @sphdecon([iter],[sph],[s])
     [s] = 2
     @sphdecon([iter],[sph],[s])
   ENDIF                           ; End of: IF ( [qsub] >= 0
 
   SYS
     echo ; echo "  Deconvolution finished for iteration: {%I0%[iter]}"

 ELSE
   ; Just copy subset volumes into: [next_vol_s]_unfilt (unfilt same as raw)
   DO [s] = 1,2                   ; Loop over two volume subsets 
     CP 
       [next_vol_s]_raw           ; Subset volume             (input)
       [next_vol_s]_unfilt        ; Unfiltered subset volume  (output)
   ENDDO
 ENDIF                            ; End of: IF [sphdecon] 
         
 ; Create FSC mask volume the first time it is used
 IQ FI [exists]                   ; Inquire if file exists
    [fsc_mask]                    ; FSC mask
 IF ( [exists] <= 0.0 ) THEN
   ; Must create FSC mask volume

   [falloff] = 10                 ; Falloff half width

   SYS
     echo ; echo "  Creating FSC mask volume: [fsc_mask]   Falloff: {%I0%[falloff]}"

   [s] = 1
   FI H [nx]                      ; Find volume dimensions
     [next_vol_s]_unfilt          ; Subset volume            (input)
     NX                           ; Cube dimension

   BL                             ; Creage blank volume
     [fsc_mask]_jnk               ; Blank volume             (output)
     [nx],[nx],[nx]               ; Size
     No                           ; Do not use average
     1.0                          ; Sphere density

   MA                             ; Circular soft mask 
     [fsc_mask]_jnk               ; Blank volume             (input)
     [fsc_mask]                   ; Soft mask volume         (output)
     [r2], 0                      ; Outer & inner radii
     Cos                          ; Cosine soft mask
     E                            ; Get background
     0.0                          ; Background
                                  ; Default mask center
     [falloff]                    ; Falloff half width

   DE                             ; Delete blank volume
     [fsc_mask]_jnk               ; Blank volume            (deleted)

   ; Display mask central slice 
   SYS
     echo "  Displaying FSC mask volume slice: [fsc_mask]_mid"

   [zcen] = [nx] / 2
   PS Z                           ; Pick slices
     [fsc_mask]                   ; Soft mask volume        (output)
     [fsc_mask]_mid               ; Slice image             (output)
     [zcen]                       ; Slice number
   DIS                            ; Display image
     [fsc_mask]_mid               ; Slice image             (input)
                                  ; Options
 ENDIF

 ; Determine masked reconstruction resolution
 FSC NEW [mvg],[mspfrg],[mresg], [mvl],[mspfrl],[mresl]  
   [next_vol]_s1_unfilt          ; Volume - subset 1         (input)
   [next_vol]_s2_unfilt          ; Volume - subset 2         (input)
   2                             ; Shell width (voxels) 
   [pixsiz]                      ; Pixel size 
   [fsc_mask]                    ; FSC mask                  (input)
   [next_m_fsc]                  ; Masked FSC doc file       (output)
   *                             ; No Gnuplot file wanted

 ; Determine unmasked reconstruction resolution
 FSC NEW [uvg],[uspfrg],[uresg], [uvl],[uspfrl],[uresl]   
   [next_vol]_s1_unfilt          ; Volume - subset 1         (input)
   [next_vol]_s2_unfilt          ; Volume - subset 2         (input)
   2                             ; Shell width (voxels)
   [pixsiz]                      ; Pixel size 
   *                             ; No FSC mask
   [next_u_fsc]                  ; Unmasked FSC doc file     (output)
   *                             ; No Gnuplot file wanted


 [key] = [iter] + 1              ; Key is output file number
 IF ( [gold-std] > 0 ) THEN
     ; Report the gold-std resolution
   SYS
     echo "  Iteration: {%I0%[iter]}  Gold Resolution, masked:{%f6.2%[mresg]}  unmasked:{%f6.2%[uresg]}  Spfrq:{%f7.2%[mspfrg]}" ; echo

   SD [key],[iter],[mresg],[uresg] ; Save resolutions in doc file
     [iter_resol]                  ; Resolution doc file      (output)

   [mspfr] = [mspfrg]
 ELSE
     ; Report the legacy resolution
   SYS
     echo "  Iteration: {%I0%[iter]}  Legacy Resolution, masked:{%f6.2%[mresl]}  unmasked:{%f6.2%[uresl]}  Spfrq:{%f7.2%[mspfrl]}" ; echo

   SD [key],[iter],[mresl],[uresl] ; Save resolutions in doc file
     [iter_resol]                  ; Resolution doc file      (output)

   [mspfr] = [mspfrl]

 ENDIF 
 SD END
   [iter_resol]                  ; Resolution doc file       (closed)

 ; Filter volume subsets to limit resolution 
 DO [s] = 1,2                    ; Loop over two volume subsets ------------ 

   IF ( [ampenhance] > 0 ) THEN     
     ; Apply amplitude enhancement filter to subset volume

     @enhance([mspfr],[pixsiz],[next-iter],[mspfr])  
       [next_vol_s]_unfilt       ; Reconstructed volume      (input)
       [next_vol_s]              ; Filtered volume           (output)
    
   ELSEIF ( [gold-std] == 0 ) THEN

     ; Using legacy resolution cuttoff,  set  pass-band & stop-band frequency 
     [pass-band] = [mspfr] 
     IF ( [mspfr] > 0.35 ) [pass-band] = 0.4  
     [stop-band] = [mspfrg] + 0.15           

     IF ( [s] == 1) THEN
       SYS
         echo "  Filtering subset volumes with Pass-band: {%f6.2%[pass-band]}  Stop-band:{%f6.2%[stop-band]} "; echo
     ENDIF

     ; Filter subset output volume to limit resolution
     FQ                          ; Quick filter          
       [next_vol_s]_unfilt       ; CTF'd volume              (input)
       [next_vol_s]              ; Filtered volume           (output)
       7                         ; Butterworth low pass filter
       [pass-band],[stop-band]   ; Pass-band, stop-band freq.

   ELSE 
     
     ; Using gold-standard  resolution cuttoff, filter subset volume using FSC doc file
     IF ( [s] == 1) THEN
       SYS
         echo "  Filtering subset volumes using masked FSC doc file: [next_m_fsc]" ; echo
     ENDIF

     FD C                        ; Fourier filter using FSC doc file 
       [next_vol_s]_unfilt       ; Reconstructed volume      (input)
       [next_vol_s]              ; Filtered volume           (output)
       [next_m_fsc]              ; Masked FSC doc file       (input)
       4                         ; Register col. for filtration values

   ENDIF

   ; Report center of gravity of subset volume (Not used for anything)
   CG [xcg],[ycg],[zcg]          ; Center of gravity ('CG PH' sometimes fails!) 
     [next_vol_s]                ; Filtered volume           (input)
     0                           ; Threshold
 
   [xsh] = -[xcg]
   [ysh] = -[ycg]
   [zsh] = -[zcg]

   SYS
     echo  "  Volume: [next_vol_s]  Off center: ({%F4.1%[xsh]}, {%F4.1%[ysh]}, {%F4.1%[zsh]})" 
   
 ENDDO                           ; End loop over resolution subsets -----------


 ; Combine the two subset volumes into overall unfiltered volume 
 SYS
   echo ; echo "  Combining two unfiltered subset volumes into overall volume: [next_vol]_unfilt"

 AD S                            ; Add two subset volumes
   [next_vol]_s*_unfilt          ; Subset volume template     (input)
   1-2                           ; Subset numbers
   [next_vol]_unfilt             ; Overall unfiltered volume  (output)

 ; Combine the two subset volumes into overall filtered volume 
 SYS
   echo "  Combining two filtered subset volumes into overall volume:   [next_vol]"

 AD S                            ; Add two subset volumes
   [next_vol]_s*                 ; Subset volume template     (input)
   1-2                           ; Subset numbers
   [next_vol]                    ; Overall filtered volume    (output)


 ; Find center of gravity of overall volume 
 CG [xcg],[ycg],[zcg]            ; Center of gravity ('CG PH' sometimes fails!)  
   [next_vol]                    ; Overall filtered volume    (input)
   0                             ; Threshold

 [xsh] = -[xcg]
 [ysh] = -[ycg]
 [zsh] = -[zcg]
 IF ( ABS([xsh]) < 0.5 ) [xsh] = 0.0
 IF ( ABS([ysh]) < 0.5 ) [ysh] = 0.0
 IF ( ABS([zsh]) < 0.5 ) [zsh] = 0.0

 ; Shift overall volume to its center of gravity 
 SH                              ; Shift volume
   [next_vol]                    ; Overall filtered volume    (input)
   [next_vol]_cent               ; Centered filtered volume   (output) 
   [xsh],[ysh],[zsh]             ; Shift distances for CofG
 SYS
   echo  "  Volume: [next_vol]_cent    Shifted: ({%F4.1%[xsh]}, {%F4.1%[ysh]}, {%F4.1%[zsh]})" ; echo

 IF ( [gold-std] == 0 ) THEN
   ; Replace subset volumes with overall volume to serve as next reference volumes

   DO [s] = 1,2                  ; Loop over two volume subsets ------------ 
     CP                          ; Copy
       [next_vol_s]              ; Subset volume              (input)
       [next_vol_s]_notgold      ; Preserved subset volume    (output)
   
     CP
       [next_vol]_cent           ; Overall centered volume    (input)
       [next_vol_s]              ; Filtered volume            (output)
   ENDDO
 ENDIF

 MY FL
 RE

 ;