Skip to content

Commit

Permalink
Got rid of SPARX/EMAN2 dependency almost completely, updated CTF and …
Browse files Browse the repository at this point in the history 
…radial processing utilities
  • Loading branch information
@rdrighetto
rdrighetto committed Feb 20, 2017
1 parent 221802f commit 02c9867
Show file tree
Hide file tree
Showing 6 changed files with 240 additions and 112 deletions.
32 changes: 26 additions & 6 deletions apps/resources/config/2dx_master.cfg
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
48 changes: 31 additions & 17 deletions scripts/proc/SPR_AverageParticles.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@
import numpy as np
import matplotlib.pyplot as plt
import sys
import ioMRC
import focus_utilities as util

def main():

Expand All @@ -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]
Expand Down Expand Up @@ -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

Expand All @@ -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]
Expand Down Expand Up @@ -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()
139 changes: 89 additions & 50 deletions scripts/proc/SPR_ExtractParticles.py.new
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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:

Expand All @@ -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:

Expand Down Expand Up @@ -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:
Expand Down
Loading

0 comments on commit 02c9867

Please sign in to comment.