([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-grps],[task],[qsub],[num-angs],[incore-yn])
; [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
;