; Calculates average power spectra 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***
; Inputs: Window size, percentage of the overlap
; distance of the window from the micrograph border, padding size
;
; SOURCE: spider/docs/techs/recon/newprogs/power.spi
;
; Edit following parameters and filenames as needed.
[winsiz] = 256 ; Window size of small pieces (Sx=Sy)
[xover] = 12 ; % of the overlap in x
[yover] = 12 ; % of the overlap in y
[xd] = 12 ; Dist. from the edge (x) (this is ok for ccd. May need larger value for film, like 50)
[yd] = 12 ; Dist. from the edge (y) (same as above)
[mask] = 225 ; Power spectrum mask radius (Angstroms); 0 = don't mask
[padsize]= 1024 ; padding of the small pieces. Important for sampling. The bigger the better, but takes longer to calculate)
; [deci] = decimation factor for 'DC S' operation
; 0 = get value from param file (key=16)
; 1 = full sized image
; 2 = 1/2 size
; 4 = 1/4 size
[deci] = 1 ; Decimation factor (0 = get value from param file). For best results in calculating defocus values
; use a decimation factor that will give you a sampling around 3 Angstrom/pixel
; ----------- Input files --------------
[params] = '../params' ; Parameter file
[sel_mic] = '../sel_micrographs' ; Micrograph selection doc file
[micgr] = '../Micrographs/mic{*****[mic]}' ; Micrographs
; ----------- Output files --------------
[outdir] = 'power' ; Power spectra output directory
[spectrum] = '[outdir]/pw_avg{*****[mic]}' ; Power spectra
[roo] = '[outdir]/roo{*****[mic]}' ; Rotational average (doc 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
UD N [nummics] ; Get number of micrographs
[sel_mic]
VM ; Make sure output dir. present
mkdir -p [outdir]
FR G ; Name for temporary incore doc file
[tmparr]tmparray
[winsizd2]=[padsize]/2
SD IC NEW ; Create an in-core doc file array
[tmparr]
4,[winsizd2] ; Creates a 4 column doc file
DO [key]=1,[nummics] ; Loop over all micrographs -------------------
UD [key], [mic]
[sel_mic] ; Micrograph selection doc file (input)
VM
echo ' 'Micrograph: {*****[mic]} Creating spectra: [spectrum] and [roo]
; ----- Get zip & format flags from parameter file (can params vary??) ---------
UD 1,[zflag] ; Get zip flag
[params]
UD 2,[fflag] ; Get tif flag
[params]
; ----------- Checks if files are gzip compressed ------------------
IF ([zflag]*[fflag].GT.0) THEN ; Both tif & zip flags set
VM ; Unzip the file
gunzip [micgr].tif.gz
ELSEIF ([zflag].GT.0) THEN ; Zip, but not tif, flag set
VM ; Unzip the file
gunzip [micgr].$DATEXT.gz
ENDIF
; ----------- Conversion based on scanner type -------------------
IF ([fflag] .eq .0) THEN ; Spider file, need to put input in _1
CP ; Place in incore file
[micgr] ; Spider file (input)
_1 ; Spider file (output)
ELSEIF ([fflag] .EQ. 1) THEN
UD 3,[nsam] ; HiScan raw file, Get X,Y size parameters
[params] ; Params file (input)
UD 4,[nrow] ;
[params] ; Params file (input)
CP FROM RAW
16 ; Bits / pixel
[micgr].tif ; Raw file (input)
[raw] ; File (input)
[nsam],[nrow] ; Size
(342) ; Header bytes
(1)
N
_2 ; Spider file (output)
AR
_2 ; Spider file (input)
_1 ; Spider file (output)
log(p1+1)
ELSEIF ([fflag] .EQ. 3) THEN
VM ; ZI tif file. Overview should always = 1
zi2spi [micgr].tif [micgr].$DATEXT 1
CP ; Place in incore file
[micgr]
_1 ; Spider file (output)
DE ; Delete the SPIDER format file
[micgr] ;
ELSEIF ([format].EQ.4) THEN
CP FROM NIKON ; Nikon Tif Scanner file
[micgr].tif ; Nikon tif file (input)
_1 ; Spider file (output)
ENDIF
; ----------- Size reduction, if any -------------------
IF ([deci].LT.1) THEN
UD 16,[deci] ; Get reduction factor from param file
[params] ; Param file (input)
ENDIF
IF ([deci].GT.1) THEN ; Reduction
DC S ; Decimate image - Sum neighbouring pixels
_1 ; Spider file (input)
_2 ; Spider file (output)
[deci],[deci] ; Decimation factor
CP ; Copy file
_2 ; Spider file (input)
_1 ; Spider file (output)
ENDIF
; ----------- Masking ---------------------------------------
IF ([mask].NE.0) THEN
UD 5,[pxsiz] ; Get pixel size
[params] ; Params file (input)
[pxsiz]=[pxsiz]*[deci] ; Adjust pixel size for decimation
[radius]=2*[pxsiz]/[mask] ; Mask radius (now in 1/px)
ENDIF
; If the output doc files already exist, delete them
DE
[roo]
; ------------------ from: power_p1 ----------------
FI H [x],[y] ; Get dimensions of the full image
_1
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) ; Number of pieces horizontal dim.(X)
[npy]=INT([noy]*(([y]-2*[yd])/[winsiz]-1)+1) ; Number 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
_1 ; Get small pieces of the input image
_5
[winsiz],[winsiz] ; Size of small pieces
[ulx],[uly] ; Coordinates of the upper left corner
RA ; Correct ramp effects
_5
_6
PD ; Pad the small pieces into larger ones for better sampling
_6
_77
[padsize],[padsize] ; padding size
B ; padding background = border average of the small window
([padsize]-[winsiz])/2+1,([padsize]-[winsiz])/2+1 ; small window centered in the larger one.
PW ; Calc power spectrum of each small piece
_77
_7
SQ ; Square image
_7
_9
[v94] = [v94]+1
IF ([v94].ge.2) THEN
AD ; Add images
_3
_9
_3
*
ELSE
CP ; Copy image
_9
_3
ENDIF
ENDDO
ENDDO
; --------- End Loops. Resulting spectrum in _3 ------------
AR ; Average over power spectra of small pieces
_3
_8
P1/[v94]
WU ; Calculate the square root
_8
_2
IF ([radius].NE.0) THEN ; Mask out center
FI H [nsam],[nrow]
_2
NSAM,NROW ; Get image dimensions
[cx] = [nsam]/2 + 1
[cy] = [nrow]/2 + 1
GP [pv] ; Get pixel intensity
_2
([cx]+10,5)
MA ; Mask image
_2
_9
(0.0,[radius]*[nsam]) ; Radius of central dot. e.g., (2 * 5A/px/200A)/x 500 = 25 pixels
D
E
[pv]
([cx],[cy])
CP ; Copy image
_9
_2
ENDIF ; End maskout if clause. Spectrum in _2
CP ; Write the output file
_2
[spectrum]
; Create the document file, with column headings & spatial freq. column
RO
_2
_3 ; 1D rotational average
LI D
_3
[tmparr] ; Doc file (output)
R ; Row output
(1) ; Row
; Labels for columns in output file
SD / Amplitude X-axis Spat.Freq.
[roo] ; Doc file (output)
[winsizd2t] = [padsize]/2 ; 1/2 padded window size
DO [v70] = 1,[winsizd2]
UD IC [v70], [v72],[v73]
[tmparr] ; Doc file (input)
[sfreq] = [v70] / (2.0 * [winsizd2t]) ; Compute spatial freq.
SD [v70],[v72],[v73],[sfreq]
[roo] ; Doc file (output)
ENDDO
; ------------------ end: power_p1 ----------------
; ----------- Rezip if necessary -------------------
IF ([zflag]*[fflag] .GT. 0) THEN ; Both tif & zip flags set
VM ; Recompress the file
gzip [micgr].tif
ELSEIF ([zflag].eq.1) THEN ; Only zip flag
VM ; Recompress the file
gzip [micgr].$DATEXT
ENDIF
ENDDO
VM
echo ' '
EN
;