From 02c98675e4337b734b8a4e9aeb3a82cfabc1265d Mon Sep 17 00:00:00 2001 From: rdrighetto Date: Mon, 20 Feb 2017 20:27:38 +0100 Subject: [PATCH] Got rid of SPARX/EMAN2 dependency almost completely, updated CTF and radial processing utilities --- apps/resources/config/2dx_master.cfg | 32 +++- scripts/proc/SPR_AverageParticles.py | 48 +++--- scripts/proc/SPR_ExtractParticles.py.new | 139 +++++++++++------- scripts/proc/focus_ctf.py | 63 ++++++-- scripts/proc/focus_utilities.py | 54 +++---- .../frealign/ExportParticleStack.script | 16 +- 6 files changed, 240 insertions(+), 112 deletions(-) diff --git a/apps/resources/config/2dx_master.cfg b/apps/resources/config/2dx_master.cfg index a5056ab2..0fd79c68 100644 --- a/apps/resources/config/2dx_master.cfg +++ b/apps/resources/config/2dx_master.cfg @@ -5278,22 +5278,42 @@ set SPR_CalculateDefocusTilted = "y" # # LABEL: Save also phase-flipped stack? # LEGEND: This parameter allows to select, if a phase-flipped version of the stack should also be created. Only effective if particle-based CTF correction is selected (SPR_WhichCTF = "Particle"). -# EXAMPLE: SPR_SavePhaseFlipped = "y" +# EXAMPLE: SPR_SavePhaseFlipped = "n" # TYPE: Bool "y;n" # LOCKED: NO # SYNC_WITH_UPPER_LEVEL: YES # MODE: 1 -set SPR_SavePhaseFlipped = "y" +set SPR_SavePhaseFlipped = "n" +# +# LABEL: Save also CTF-multiplied stack? +# LEGEND: This parameter allows to select, if a CTF-multiplied version of the stack should also be created. Only effective if particle-based CTF correction is selected (SPR_WhichCTF = "Particle"). +# EXAMPLE: SPR_SaveCTFMultiplied = "y" +# TYPE: Bool "y;n" +# LOCKED: NO +# SYNC_WITH_UPPER_LEVEL: YES +# MODE: 1 +set SPR_SaveCTFMultiplied = "y" # # LABEL: Save also Wiener-filtered stack? # LEGEND: This parameter allows to select, if a Wiener-filtered version of the stack should also be created. Only effective if particle-based CTF correction is selected (SPR_WhichCTF = "Particle"). -# EXAMPLE: SPR_SaveWienerFiltered = "y" +# EXAMPLE: SPR_SaveWienerFiltered = "n" # TYPE: Bool "y;n" # LOCKED: NO # SYNC_WITH_UPPER_LEVEL: YES # MODE: 1 set SPR_SaveWienerFiltered = "n" # +# LABEL: Wiener-filter constant. Only used if Wiener-filtering is applied as CTF correction. +# LEGEND: The noise constant for the Wiener filter. Try 0.4. +# EXAMPLE: SPR_WienerConstant = "0.4" +# HELP: http://2dx.org/documentation/2dx-software/parameters/Defocus +# TYPE: Float "MIN=0.0;MAX=1.0;DEFAULT=0.4" +# LOCKED: NO +# SYNC_WITH_UPPER_LEVEL: YES +# ISWRONG: NO +# MODE: 1 +set SPR_WienerConstant = "0.4" +# # LABEL: Ignore 'lat2' images? # LEGEND: This parameter allows to select, if images from a second lattice should be ignored in the single-particle dataset. In most cases this should be set to 'y'. # EXAMPLE: SPR_Ignore_Lat2 = "y" @@ -5351,7 +5371,7 @@ set SPR_WhichCTF = "Micrograph" # LABEL: Which CTF-corrected stack to use for crystal averaging? # LEGEND: This parameter allows to select, which of the different possible CTF-corrected stacks to use for averaging particles on a per-crystal basis. "ctfcor" is the stack created from the CTF-corrected micrographs from 2dx_image (image_ctfcor.mrc files); "phase-flipped" if CTF was corrected on each particle individually by phase-flipping; "wiener-filtered" if CTF was corrected on each particle individually by Wiener-filtering. # EXAMPLE: SPR_WhichStack = "ctfcor" -# TYPE: Drop_Down_Menu "ctfcor;phase-flipped;wiener-filtered" +# TYPE: Drop_Down_Menu "ctfcor;phase-flipped;ctf-multiplied;wiener-filtered" # LOCKED: NO # SYNC_WITH_UPPER_LEVEL: YES # MODE: 1 @@ -5360,7 +5380,7 @@ set SPR_WhichStack = "ctfcor" # LABEL: Which type of stack? # LEGEND: This parameter allows to select, which of the different available stacks to use for averaging particles on a per-crystal basis. "non-corrected" is the non-CTF-CORRECTED particle stack, e.g. "particles.mrcs", which you should use if you want FREALIGN and RELION to perform CTF correction internally; "ctfcor" is the stack created from the CTF-corrected micrographs from 2dx_image (image_ctfcor.mrc files); "phase-flipped" if CTF was corrected on each particle individually by phase-flipping; "wiener-filtered" if CTF was corrected on each particle individually by Wiener-filtering. # EXAMPLE: SPR_WhichStackAll = "non-corrected" -# TYPE: Drop_Down_Menu "non-corrected;ctfcor;phase-flipped;wiener-filtered" +# TYPE: Drop_Down_Menu "non-corrected;ctfcor;phase-flipped;ctf-multiplied;wiener-filtered" # LOCKED: NO # SYNC_WITH_UPPER_LEVEL: YES # MODE: 1 @@ -5369,7 +5389,7 @@ set SPR_WhichStackAll = "non-corrected" # LABEL: Inherit pre-refinement results from which stack? # LEGEND: To prepare the files to do a full Single-Particle Refinement and/or 3D Classification, the optimized parameters obtained in the pre-refinement will be inherited. Choose here from which pre-refined stack you want to inherit the alignment parameters. "ctfcor" is the stack created from the CTF-corrected micrographs from 2dx_image (image_ctfcor.mrc files); "phase-flipped" if CTF was corrected on each particle individually by phase-flipping; "wiener-filtered" if CTF was corrected on each particle individually by Wiener-filtering. # EXAMPLE: SPR_WhichStackInherit = "ctfcor" -# TYPE: Drop_Down_Menu "ctfcor;phase-flipped;wiener-filtered" +# TYPE: Drop_Down_Menu "ctfcor;phase-flipped;ctf-multiplied;wiener-filtered" # LOCKED: NO # SYNC_WITH_UPPER_LEVEL: YES # MODE: 1 diff --git a/scripts/proc/SPR_AverageParticles.py b/scripts/proc/SPR_AverageParticles.py index 08207c58..e9cbe182 100755 --- a/scripts/proc/SPR_AverageParticles.py +++ b/scripts/proc/SPR_AverageParticles.py @@ -14,6 +14,8 @@ import numpy as np import matplotlib.pyplot as plt import sys +import ioMRC +import focus_utilities as util def main(): @@ -38,8 +40,9 @@ def main(): f = open(stack_path+stack_rootname+'_crystal-avg_1_r1-%.4d.par' % this_thread, 'w+') - first = spx.EMData() - first.read_image(stack_file, 0) + # first = spx.EMData() + # first.read_image(stack_file, 0) + header = ioMRC.readMRCHeader( stack_file ) par = np.loadtxt(stack_path+stack_rootname+'_1_r1.par', comments='C') labels = par[:,7] @@ -71,19 +74,24 @@ def main(): img_list = np.where(labels == x)[0] - avg = spx.EMData(first.get_xsize(),first.get_ysize()) + # avg = spx.EMData(first.get_xsize(),first.get_ysize()) + # ioMRC header is Z,Y,X: + avg = np.zeros( [header['dimensions'][2], header['dimensions'][1]] ) if do_frc: plt.figure() - odd = spx.EMData(first.get_xsize(),first.get_ysize()) - even = spx.EMData(first.get_xsize(),first.get_ysize()) + # odd = spx.EMData(first.get_xsize(),first.get_ysize()) + # even = spx.EMData(first.get_xsize(),first.get_ysize()) + odd = np.zeros( [header['dimensions'][2], header['dimensions'][1]] ) + even = np.zeros( [header['dimensions'][2], header['dimensions'][1]] ) for i in img_list: - img = spx.EMData() - img.read_image(stack_file, int(i)) + # img = spx.EMData() + # img.read_image(stack_file, int(i)) + img = ioMRC.readMRC( stack_file, idx=i ) avg += img @@ -97,9 +105,15 @@ def main(): even += img - avg = NormalizeStack([avg], sigma)[0] + # if normalize_box: + + # box = NormalizeStack([box], sigma)[0] + avg = util.NormalizeImg( avg, std=sigma ) + # avg = NormalizeStack([avg], sigma)[0] + + # avg.write_image(stack_path+stack_rootname+'_crystal-avg-%.4d.mrcs' % this_thread, j-1) + ioMRC.writeMRC( avg, stack_path+stack_rootname+'_crystal-avg-%.4d.mrcs' % this_thread, dtype='float32', idx=j-1 ) - avg.write_image(stack_path+stack_rootname+'_crystal-avg-%.4d.mrcs' % this_thread, j-1) # Write .par file with the parameters for each particle in the dataset: # print >>f, ' %d' % (x),' %.2f' % par[img_list[0],1],' %.2f' % par[img_list[0],2],' %.2f' % par[img_list[0],3],' %.2f' % par[img_list[0],4],' %.2f' % par[img_list[0],5],' %d' % par[img_list[0],6],' %d' % par[img_list[0],7],' %.2f' % par[img_list[0],8],' %.2f' % par[img_list[0],9],' %.2f' % par[img_list[0],10],' %.2f' % par[img_list[0],11],' %d' % par[img_list[0],12],' %.4f' % par[img_list[0],13],' %.2f' % par[img_list[0],14],' %.2f' % par[img_list[0],15] @@ -134,17 +148,17 @@ def main(): # print ':: ' -def NormalizeStack(stack, sigma): -# Normalizes a stack of images to zero-mean and standard deviation given by sigma: +# def NormalizeStack(stack, sigma): +# # Normalizes a stack of images to zero-mean and standard deviation given by sigma: - stack_norm = [] - for i in range(len(stack)): +# stack_norm = [] +# for i in range(len(stack)): - zero_float = stack[i] - stack[i]['mean'] - zero_float.mult(sigma/zero_float['sigma']) - stack_norm.append(zero_float) +# zero_float = stack[i] - stack[i]['mean'] +# zero_float.mult(sigma/zero_float['sigma']) +# stack_norm.append(zero_float) - return stack_norm +# return stack_norm if __name__ == '__main__': main() \ No newline at end of file diff --git a/scripts/proc/SPR_ExtractParticles.py.new b/scripts/proc/SPR_ExtractParticles.py.new index fb59a895..044ba70c 100644 --- a/scripts/proc/SPR_ExtractParticles.py.new +++ b/scripts/proc/SPR_ExtractParticles.py.new @@ -56,39 +56,44 @@ def main(): else: save_phase_flipped = False if sys.argv[18] == 'y': + save_ctf_multiplied = True + else: + save_ctf_multiplied = False + if sys.argv[19] == 'y': save_wiener_filtered = True else: save_wiener_filtered = False - sigma = float(sys.argv[19]) # Sigma for normalization of the windowed images (if normalize_box == True) - if sys.argv[20] == 'Defocus/Lattice': + wiener_constant = float(sys.argv[20]) + sigma = float(sys.argv[21]) # Sigma for normalization of the windowed images (if normalize_box == True) + if sys.argv[22] == 'Defocus/Lattice': tiltgeom = '' - elif sys.argv[20] == 'Defocus': + elif sys.argv[22] == 'Defocus': tiltgeom = 'DEFOCUS_' - elif sys.argv[20] == 'Lattice': + elif sys.argv[22] == 'Lattice': tiltgeom = 'LATTICE_' - elif sys.argv[20] == 'Merge': + elif sys.argv[22] == 'Merge': tiltgeom = 'MERGE_' - if sys.argv[21] == 'Micrograph': + if sys.argv[23] == 'Micrograph': ctfcor = True stack_rootname = stack_rootname + '_ctfcor' else: ctfcor = False - if sys.argv[22] == 'y': + if sys.argv[24] == 'y': save_pick_fig = True else: save_pick_fig = False - if sys.argv[23] == 'y': + if sys.argv[25] == 'y': use_masked_image = True else: use_masked_image = False - if sys.argv[24] == 'y': + if sys.argv[26] == 'y': shuffle_order = True else: shuffle_order = False - n_threads = int(sys.argv[25]) + n_threads = int(sys.argv[27]) if n_threads < 1: n_threads = 1 - this_thread = int(sys.argv[26]) + this_thread = int(sys.argv[28]) # End arguments @@ -377,33 +382,37 @@ def main(): # spx.EMNumPy.em2numpy(box).tofile(mrcf) ioMRC.writeMRC( box, stack_path+stack_rootname+'-%.4d.mrcs' % this_thread, dtype='float32', idx=idx ) - if (save_phase_flipped or save_wiener_filtered) and not ctfcor: + # if (save_phase_flipped or save_wiener_filtered) and not ctfcor: - # Convert CTF parameters to SPARX convention: - defocus = (RLDEF1+RLDEF2)/2 - ast = RLDEF1-RLDEF2 - if params['AST_ANGLE'] < 0.0: + # # Convert CTF parameters to SPARX convention: + # defocus = (RLDEF1+RLDEF2)/2 + # ast = RLDEF1-RLDEF2 + # if params['AST_ANGLE'] < 0.0: - astang = 360.0 + params['AST_ANGLE'] + # astang = 360.0 + params['AST_ANGLE'] - else: + # else: - astang = params['AST_ANGLE'] + # astang = params['AST_ANGLE'] - # Generate SPARX CTF object: - p = [defocus * 1e-4, microscope_cs, microscope_voltage, apix, 0, ampcontrast * 100, ast * 1e-4, astang] + # # Generate SPARX CTF object: + # p = [defocus * 1e-4, microscope_cs, microscope_voltage, apix, 0, ampcontrast * 100, ast * 1e-4, astang] - spx.set_ctf(box, p) + # spx.set_ctf(box, p) - ctf = spx.generate_ctf(p) + # ctf = spx.generate_ctf(p) # Phase-flip the image: if save_phase_flipped and not ctfcor: # Apply CTF correction on whole micrograph to reduce delocalization effects: - imgctfcor = spx.filt_ctf( img, ctf, binary=1 ) + # imgctfcor = spx.filt_ctf( img, ctf, binary=1 ) + + # boxctfcor = spx.Util.window( imgctfcor, int( box_size ), int( box_size ), 1, int( round( x[i] ) ), int( round( y[i] ) ) ) + imgctfcor = CTF.CorrectCTF( img, DF1=RLDEF1, DF2=RLDEF2, AST=params['AST_ANGLE'], WGH=ampcontrast, apix=apix, Cs=microscope_cs, kV=microscope_voltage, ctftype=0, return_ctf=False, invert_contrast=False ) + + boxctfcor = imgctfcor[yi-w/2-box_size/2:yi-w/2+box_size/2, xi-w/2-box_size/2:xi-w/2+box_size/2] - boxctfcor = spx.Util.window( imgctfcor, int( box_size ), int( box_size ), 1, int( round( x[i] ) ), int( round( y[i] ) ) ) if normalize_box: @@ -421,13 +430,43 @@ def main(): # spx.EMNumPy.em2numpy(boxctfcor).tofile( mrcf ) ioMRC.writeMRC( boxctfcor, stack_path+stack_rootname+'_phase-flipped-%.4d.mrcs' % this_thread, dtype='float32', idx=idx ) + # CTF-multiply the image: + if save_ctf_multiplied and not ctfcor: + + # # Apply CTF correction on whole micrograph to reduce delocalization effects: + # imgctfcor = spx.filt_ctf( img, ctf, binary=0 ) + + # boxctfcor = spx.Util.window( imgctfcor, int( box_size ), int( box_size ), 1, int( round( x[i] ) ), int( round( y[i] ) ) ) + imgctfcor = CTF.CorrectCTF( img, DF1=RLDEF1, DF2=RLDEF2, AST=params['AST_ANGLE'], WGH=ampcontrast, apix=apix, Cs=microscope_cs, kV=microscope_voltage, ctftype=1, return_ctf=False, invert_contrast=False ) + + boxctfcor = imgctfcor[yi-w/2-box_size/2:yi-w/2+box_size/2, xi-w/2-box_size/2:xi-w/2+box_size/2] + + if normalize_box: + + # boxctfcor = NormalizeStack([boxctfcor], sigma)[0] + boxctfcor = util.NormalizeImg( boxctfcor, std=sigma ) + + # Write image to the particle stack: + # if idx == 0: + # # If this is the first image, we initiate as a normal .mrcs stack. + # boxctfcor.write_image( stack_path+stack_rootname+'_wiener-filtered-%.4d.mrcs' % this_thread, idx ) + + # else: + # # Subsequent images are directly appended to the file: + # with open( stack_path+stack_rootname+'_wiener-filtered-%.4d.mrcs' % this_thread, 'ab' ) as mrcf: + # spx.EMNumPy.em2numpy(boxctfcor).tofile( mrcf ) + ioMRC.writeMRC( boxctfcor, stack_path+stack_rootname+'_ctf-multiplied-%.4d.mrcs' % this_thread, dtype='float32', idx=idx ) + # Wiener-filter the image: if save_wiener_filtered and not ctfcor: - # Apply CTF correction on whole micrograph to reduce delocalization effects: - imgctfcor = spx.filt_ctf( img, ctf, binary=0 ) + # # Apply CTF correction on whole micrograph to reduce delocalization effects: + # imgctfcor = spx.filt_ctf( img, ctf, binary=0 ) + + # boxctfcor = spx.Util.window( imgctfcor, int( box_size ), int( box_size ), 1, int( round( x[i] ) ), int( round( y[i] ) ) ) + imgctfcor = CTF.CorrectCTF( img, DF1=RLDEF1, DF2=RLDEF2, AST=params['AST_ANGLE'], WGH=ampcontrast, apix=apix, Cs=microscope_cs, kV=microscope_voltage, ctftype=2, return_ctf=False, invert_contrast=False, C=wiener_constant ) - boxctfcor = spx.Util.window( imgctfcor, int( box_size ), int( box_size ), 1, int( round( x[i] ) ), int( round( y[i] ) ) ) + boxctfcor = imgctfcor[yi-w/2-box_size/2:yi-w/2+box_size/2, xi-w/2-box_size/2:xi-w/2+box_size/2] if normalize_box: @@ -464,30 +503,30 @@ def main(): print '\nBoxed %d/%d CC peaks from micrograph %d/%d.\n' % (m, i+1, n, N) - # Update the counts in the MRC headers: - # First the normal particle stack: - header = ioMRC.readMRCHeader( stack_path+stack_rootname+'-%.4d.mrcs' % this_thread ) - header['dimensions'][0] = idx - # Now we write the header back: - with open( stack_path+stack_rootname+'-%.4d.mrcs' % this_thread, 'rb+' ) as mrcf: - ioMRC.writeMRCHeader( mrcf, header, endchar = '<' ) - # Then the phase-flipped: - if save_phase_flipped and not ctfcor: + # # Update the counts in the MRC headers: + # # First the normal particle stack: + # header = ioMRC.readMRCHeader( stack_path+stack_rootname+'-%.4d.mrcs' % this_thread ) + # header['dimensions'][0] = idx + # # Now we write the header back: + # with open( stack_path+stack_rootname+'-%.4d.mrcs' % this_thread, 'rb+' ) as mrcf: + # ioMRC.writeMRCHeader( mrcf, header, endchar = '<' ) + # # Then the phase-flipped: + # if save_phase_flipped and not ctfcor: - header = ioMRC.readMRCHeader( stack_path+stack_rootname+'_phase-flipped-%.4d.mrcs' % this_thread ) - header['dimensions'][0] = idx - # Now we write the header back: - with open( stack_path+stack_rootname+'_phase-flipped-%.4d.mrcs' % this_thread, 'rb+' ) as mrcf: - ioMRC.writeMRCHeader( mrcf, header, endchar = '<' ) - - # Then the Wiener-filtered: - if save_wiener_filtered and not ctfcor: + # header = ioMRC.readMRCHeader( stack_path+stack_rootname+'_phase-flipped-%.4d.mrcs' % this_thread ) + # header['dimensions'][0] = idx + # # Now we write the header back: + # with open( stack_path+stack_rootname+'_phase-flipped-%.4d.mrcs' % this_thread, 'rb+' ) as mrcf: + # ioMRC.writeMRCHeader( mrcf, header, endchar = '<' ) + + # # Then the Wiener-filtered: + # if save_wiener_filtered and not ctfcor: - header = ioMRC.readMRCHeader( stack_path+stack_rootname+'_wiener-filtered-%.4d.mrcs' % this_thread ) - header['dimensions'][0] = idx - # Now we write the header back: - with open( stack_path+stack_rootname+'_wiener-filtered-%.4d.mrcs' % this_thread, 'rb+' ) as mrcf: - ioMRC.writeMRCHeader( mrcf, header, endchar = '<' ) + # header = ioMRC.readMRCHeader( stack_path+stack_rootname+'_wiener-filtered-%.4d.mrcs' % this_thread ) + # header['dimensions'][0] = idx + # # Now we write the header back: + # with open( stack_path+stack_rootname+'_wiener-filtered-%.4d.mrcs' % this_thread, 'rb+' ) as mrcf: + # ioMRC.writeMRCHeader( mrcf, header, endchar = '<' ) # print '<<@progress: %d>>' % round(n*100.0/N) # Report progress to the GUI: diff --git a/scripts/proc/focus_ctf.py b/scripts/proc/focus_ctf.py index 20358220..91cab264 100644 --- a/scripts/proc/focus_ctf.py +++ b/scripts/proc/focus_ctf.py @@ -2,8 +2,11 @@ # Author: Ricardo Righetto # E-mail: ricardo.righetto@unibas.ch +# TO-DO: modify FFT-related functions to use only half the spectrum (real data, hermitian symmetry) + import numpy as np import focus_utilities +import copy try: # SciPy FFT pack is faster than NumPy's: import scipy.fftpack as fft @@ -16,32 +19,44 @@ def CTF( imsize = [100, 100], DF1 = 1000.0, DF2 = None, AST = 0.0, WGH = 0.10, C # Generates 2D CTF function # Underfocus is positive following conventions of FREALIGN and most of the packages out there (in Angstroms). # B is B-factor +# ONLY WORKS FOR SQUARE IMAGES! + + if imsize[0] != imsize[1]: + + raise RuntimeError("Error: CTF model currently only implemented for square images!") Cs *= 1e7 # Convert Cs to Angstroms + if DF2 == None: DF2 = DF1 - AST *= np.pi / 180.0 + else: + + # NOTATION FOR DEFOCUS1, DEFOCUS2, ASTIGMASTISM BELOW IS INVERTED DUE TO NUMPY CONVENTION: + DF1, DF2 = DF2, DF1 + + AST *= -np.pi / 180.0 WL = ElectronWavelength( kV ) w1 = np.sqrt( 1 - WGH*WGH ) w2 = WGH - rmesh,amesh = focus_utilities.RadialIndices( imsize, normalize=True ) + rmesh,amesh = focus_utilities.RadialIndices( imsize, normalize=True, rounding=True ) - rmesh = rmesh / apix + rmesh /= apix rmesh2 = rmesh*rmesh - # NOTATION BELOW IS INVERTED DUE TO NUMPY CONVENTION: - DF = 0.5 * (DF1 + DF2 + (DF2 - DF1) * np.cos( 2.0 * (amesh - AST) ) ) + + # From Mindell & Grigorieff, JSB 2003: + DF = 0.5 * (DF1 + DF2 + (DF1 - DF2) * np.cos( 2.0 * (amesh - AST) ) ) import warnings with warnings.catch_warnings(): warnings.filterwarnings( "ignore", category=RuntimeWarning ) - Xr = np.pi * WL * rmesh2 * ( DF - 1 / (2 * WL*WL * rmesh2 * Cs) ) + Xr = np.pi * WL * rmesh2 * ( DF - 1.0 / (2.0 * WL*WL * rmesh2 * Cs) ) Xr = np.nan_to_num( Xr ) sinXr = np.sin( Xr ) @@ -51,14 +66,18 @@ def CTF( imsize = [100, 100], DF1 = 1000.0, DF2 = None, AST = 0.0, WGH = 0.10, C # CTFim = CTFreal + CTFimag*1j CTFim = -w1 * sinXr - w2 * cosXr - CTFim = CTFim * np.exp( -B * ( rmesh2 ) / 4 ) + + if B != 0.0: # Apply B-factor only if necessary: + + CTFim = CTFim * np.exp( -B * ( rmesh2 ) / 4 ) return CTFim def ElectronWavelength( kV = 300.0 ): # Returns electorn wavelength in Angstroms - kV *= 1e3 # ensure Kilovolts for below formula - return 12.26 / np.sqrt( kV + 0.9785 * kV*kV / ( 1e6 ) ) + # kV *= 1e3 # ensure Kilovolts for below formula + # return 12.2639 / np.sqrt( kV + 0.97845 * kV*kV / ( 1e6 ) ) + return 12.2639 / np.sqrt( kV * 1e3 + 0.97845 * kV*kV ) def CorrectCTF( img, DF1 = 1000.0, DF2 = None, AST = 0.0, WGH = 0.10, invert_contrast = False, Cs = 2.7, kV = 300.0, apix = 1.0, B = 0.0, ctftype = 0, C = 1.0, return_ctf = False ): # Applies CTF correction to image @@ -99,6 +118,30 @@ def CorrectCTF( img, DF1 = 1000.0, DF2 = None, AST = 0.0, WGH = 0.10, invert_con raise ValueError( "Error: Type of CTF correction must be 0 (phase-flipping), 1 (CTF multiplication) or 2 (Wiener filtering). ctftype = %d " % ctftype ) + # if return_half: + + # # AmpHalf = np.zeros( img.shape ) + # CTFnorm = CTFim**2 + # AmpHalf = np.abs( FT )**2 + # AmpHalf = AmpHalf - AmpHalf.mean() + CTFnorm.mean() + # AmpHalf = AmpHalf*CTFnorm.std()/AmpHalf.std() + # # CTFnorm -= CTFnorm.mean() + # # CTFnorm /= CTFnorm.std() + # AmpHalf[:,img.shape[1]/2:] = CTFnorm[:,img.shape[1]/2:] + # # AmpHalf = AmpHalf**2 + + # if return_ctf and return_half: + + # return CTFcor.real, CTFim, AmpHalf + + # elif return_ctf: + + # return CTFcor.real, CTFim + + # elif return_half: + + # return CTFcor.real, AmpHalf + if return_ctf: return CTFcor.real, CTFim @@ -109,5 +152,3 @@ def CorrectCTF( img, DF1 = 1000.0, DF2 = None, AST = 0.0, WGH = 0.10, invert_con - - diff --git a/scripts/proc/focus_utilities.py b/scripts/proc/focus_utilities.py index 757c0a11..cc6f792a 100644 --- a/scripts/proc/focus_utilities.py +++ b/scripts/proc/focus_utilities.py @@ -11,19 +11,20 @@ import numpy.fft as fft -def RadialIndices( imsize = [100, 100], rounding=False, normalize=False ): +def RadialIndices( imsize = [100, 100], rounding=True, normalize=False ): # Returns radius and angles for each pixel (or voxel) in a 2D image or 3D volume of shape = imsize # For 2D returns the angle with the horizontal x- axis # For 3D returns the angle with the horizontal x,y plane # If imsize is a scalar, will default to 2D. -# Rounding is to ensure "perfect" radial symmetry, desirable for applications in real space. -# Norm will normalize the radius to values between 0 and 1. +# Rounding is to ensure "perfect" radial symmetry, desirable in most applications. +# Normalize=True will normalize the radius to values between 0.0 and 1.0 (in relation to the object's shortest dimension). +# Note: the radii values returned by this function are consistent with the sampling frequency in SciPy/NumPy if np.isscalar(imsize): imsize = [imsize, imsize] - if len(imsize) > 3: + if len( imsize ) > 3: raise ValueError ( "Object should not have dimensions larger than 3: len(imsize) = %d " % len(imsize)) @@ -31,49 +32,50 @@ def RadialIndices( imsize = [100, 100], rounding=False, normalize=False ): with warnings.catch_warnings(): warnings.filterwarnings( "ignore", category=RuntimeWarning ) - if len(imsize) == 2: + # d = np.ones( len( imsize ) ) - if normalize: + # if not normalize: - [xmesh, ymesh] = np.mgrid[-imsize[0]/2:imsize[0]/2, -imsize[1]/2:imsize[1]/2].astype(np.float) + # for i in np.arange( len( imsize ) ) - xmesh /= imsize[0] - ymesh /= imsize[1] + # d[i] = 1./imsize[i] - else: + m = np.mod(imsize, 2) # Check if dimensions are odd or even - [xmesh, ymesh] = np.mgrid[-imsize[0]/2:imsize[0]/2, -imsize[1]/2:imsize[1]/2] - # xmesh += 1 - # ymesh += 1 + if len(imsize) == 2: + + # [xmesh, ymesh] = np.mgrid[-imsize[0]/2+1:imsize[0]/2+1, -imsize[1]/2+1:imsize[1]/2+1].astype(np.float) + # The definition below is consistent with scipy/numpy fft.fftfreq: + [xmesh, ymesh] = np.mgrid[-imsize[0]//2+m[0]:(imsize[0]-1)//2+1, -imsize[1]//2+m[1]:(imsize[1]-1)//2+1].astype(np.float) rmesh = np.sqrt( xmesh*xmesh + ymesh*ymesh ) + amesh = np.arctan( ymesh / xmesh ) else: - if normalize: + # [xmesh, ymesh, zmesh] = np.mgrid[-imsize[0]/2:imsize[0]/2, -imsize[1]/2:imsize[1]/2, -imsize[2]/2:imsize[2]/2].astype(np.float) + # The definition below is consistent with scipy/numpy fft.fftfreq: + [xmesh, ymesh, zmesh] = np.mgrid[-imsize[0]//2+m[0]:(imsize[0]-1)//2+1, -imsize[1]//2+m[1]:(imsize[1]-1)//2+1, -imsize[2]//2+m[2]:(imsize[2]-1)//2+1].astype(np.float) - [xmesh, ymesh, zmesh] = np.mgrid[-imsize[0]/2:imsize[0]/2, -imsize[1]/2:imsize[1]/2, -imsize[2]/2:imsize[2]/2].astype(np.float) + rmesh = np.sqrt( xmesh*xmesh + ymesh*ymesh + zmesh*zmesh ) - xmesh /= imsize[0] - ymesh /= imsize[1] - zmesh /= imsize[2] + amesh = np.arccos( zmesh / rmesh ) - else: - [xmesh, ymesh, zmesh] = np.mgrid[-imsize[0]/2:imsize[0]/2, -imsize[1]/2:imsize[1]/2, -imsize[2]/2:imsize[2]/2] + if rounding: - rmesh = np.sqrt( xmesh*xmesh + ymesh*ymesh + zmesh*zmesh ) - amesh = np.arccos( zmesh / rmesh ) - # amesh[imsize[0]/2, imsize[1]/2, imsize[2]/2] = 0.0 + rmesh = np.round( rmesh ) - if rounding and not normalize: + if normalize: - return np.round( rmesh ), np.nan_to_num( amesh ) + rmesh /= np.min( imsize ) else: - return rmesh, np.nan_to_num( amesh ) + rmesh = rmesh.astype(np.int) + + return rmesh, np.nan_to_num( amesh ) def RotationalAverage( img ): # Compute the rotational average of a 2D image or 3D volume diff --git a/scripts/singleparticle/frealign/ExportParticleStack.script b/scripts/singleparticle/frealign/ExportParticleStack.script index 3b8a3f58..e5391664 100755 --- a/scripts/singleparticle/frealign/ExportParticleStack.script +++ b/scripts/singleparticle/frealign/ExportParticleStack.script @@ -29,7 +29,9 @@ # DISPLAY: SPR_ShuffleParticleOrder # DISPLAY: SPR_CalculateDefocusTilted # DISPLAY: SPR_SavePhaseFlipped +# DISPLAY: SPR_SaveCTFMultiplied # DISPLAY: SPR_SaveWienerFiltered +# DISPLAY: SPR_WienerConstant # DISPLAY: SPR_WhichTiltGeometry # DISPLAY: CS # DISPLAY: KV @@ -63,7 +65,9 @@ set SPR_InvertContrast = "" set SPR_NormalizeBox = "" set SPR_CalculateDefocusTilted = "" set SPR_SavePhaseFlipped = "" +set SPR_SaveCTFMultiplied = "" set SPR_SaveWienerFiltered = "" +set SPR_WienerConstant = "" set SPR_SigmaNorm = "" set SPR_WhichTiltGeometry = "" set SPR_WhichCTF = "" @@ -119,7 +123,7 @@ endif set t = 1 set pids = "" while ( $t <= ${Thread_Number} ) - ${app_python} ${proc_2dx}/SPR_ExtractParticles.py ${SPR_IMGS_DIR} ${SPR_DIR}/2dx_merge_dirfile-unique.dat ${SPR_PICKING_DIR} ${SPR_STACKS_DIR} ${SPR_STACK_ROOTNAME} ${SPR_boxsize} ${SPR_PhaseShift} ${sample_pixel} ${KV} ${CS} ${phacon} ${magnification} ${SPR_SIGCC} ${SPR_InvertContrast} ${SPR_NormalizeBox} ${SPR_CalculateDefocusTilted} ${SPR_SavePhaseFlipped} ${SPR_SaveWienerFiltered} ${SPR_SigmaNorm} ${SPR_WhichTiltGeometry} ${SPR_WhichCTF} ${SPR_SavePickFig} ${use_masked_image} ${SPR_ShuffleParticleOrder} ${Thread_Number} ${t} & + ${app_python} ${proc_2dx}/SPR_ExtractParticles.py ${SPR_IMGS_DIR} ${SPR_DIR}/2dx_merge_dirfile-unique.dat ${SPR_PICKING_DIR} ${SPR_STACKS_DIR} ${SPR_STACK_ROOTNAME} ${SPR_boxsize} ${SPR_PhaseShift} ${sample_pixel} ${KV} ${CS} ${phacon} ${magnification} ${SPR_SIGCC} ${SPR_InvertContrast} ${SPR_NormalizeBox} ${SPR_CalculateDefocusTilted} ${SPR_SavePhaseFlipped} ${SPR_SaveCTFMultiplied} ${SPR_SaveWienerFiltered} ${SPR_WienerConstant} ${SPR_SigmaNorm} ${SPR_WhichTiltGeometry} ${SPR_WhichCTF} ${SPR_SavePickFig} ${use_masked_image} ${SPR_ShuffleParticleOrder} ${Thread_Number} ${t} & set pids = "$pids $!" @ t++ end @@ -144,6 +148,9 @@ else if ( ${SPR_SavePhaseFlipped} == "y" ) then set stacks = "$stacks _phase-flipped" endif + if ( ${SPR_SaveCTFMultiplied} == "y" ) then + set stacks = "$stacks _ctf-multiplied" + endif endif # foreach s ( raw _ctfcor _phase-flipped _wiener-filtered ) @@ -193,11 +200,12 @@ foreach s ( ${stacks} ) else tail -c +1025 ${SPR_STACKS_DIR}/${SPR_STACK_ROOTNAME}${i}-${num}.mrcs >> ${SPR_STACKS_DIR}/${SPR_STACK_ROOTNAME}${i}.mrcs endif - rm ${SPR_STACKS_DIR}/${SPR_STACK_ROOTNAME}${i}-${num}.mrcs @ t++ end + echo "::Cleaning up partial ${SPR_STACK_ROOTNAME} stacks..." + rm ${SPR_STACKS_DIR}/${SPR_STACK_ROOTNAME}${i}-*.mrcs # The below script call is necessary to ensure correct number of images in the header of the final .mrcs files. If the --stats option makes it too slow, one can (almost always) safely omit it. # ${proc_2dx}/MRCheaderUpdate.py ${SPR_STACKS_DIR}/${SPR_STACK_ROOTNAME}${i}.mrcs --stats @@ -266,6 +274,10 @@ if ( ${SPR_SaveWienerFiltered} == "y" && ${SPR_WhichCTF} != "Micrograph" ) then cp ${SPR_STACKS_DIR}/${SPR_STACK_ROOTNAME}_1_r1.par ${SPR_STACKS_DIR}/${SPR_STACK_ROOTNAME}_wiener-filtered_1_r1.par endif +if ( ${SPR_SaveCTFMultiplied} == "y" && ${SPR_WhichCTF} != "Micrograph" ) then + cp ${SPR_STACKS_DIR}/${SPR_STACK_ROOTNAME}_1_r1.par ${SPR_STACKS_DIR}/${SPR_STACK_ROOTNAME}_ctf-multiplied_1_r1.par +endif + if ( ${SPR_SavePickFig} == "y" ) then foreach i (${SPR_PICKING_DIR}/*.png)