TrimAliFile

#!/usr/bin/env python

def command_line_options(parser):       """ add command line options here"""       parser.add_option("-x", "--x_alicenter", dest="c_XaliCenter", 			default="0", type="int", 			help="particle X-center (default=0)")       parser.add_option("-y", "--y_alicenter", dest="c_YaliCenter", 			default="0", type="int", 			help="particle Y-center (default=0)")	           parser.add_option("-z", "--z_center", dest="c_ZCenter", 			default="0", type="int", 			help="particle Y-center (default=0)")       parser.add_option("-w", "--thick", dest="c_Thick", 			default="300", type="int", 			help="thickness of the final volume (default=300)")       parser.add_option("-b", "--box_size", dest="c_BoxSize", 			default="300", type="int", 			help="Box Size (default=300)")       parser.add_option("-t", "--tilt", dest="c_tiltfile", default='serie.tlt', 			type="string", help="Docfile FileName tilt angle must be in the second column(default=tilt.doc)")       parser.add_option("-i", "--inputalifile", dest="c_inputalifile", default='serie.ali', 			type="string", help="Input Ali FileName (default=serie.ali)")       parser.add_option("-o", "--outputfile", dest="c_outputfile", default='trim_', 			type="string", help="Output root FileName (default=trim_)")

#path to xmipp_protocols import os,sys,shutil scriptdir=os.path.split(os.path.dirname(os.popen('which xmipp_protocols','r').read()))[0]+'/protocols' sys.path.append(scriptdir) # add default search path

#process command line import optparse parser = optparse.OptionParser("usage: %prog [options] ") command_line_options(parser)

#arg the values that do not requiere '-flags', that is, the python script (options, args) = parser.parse_args()

XCenter = options.c_XaliCenter YCenter = options.c_YaliCenter ZCenter = options.c_ZCenter Thick = options.c_Thick USize   = options.c_BoxSize docfile = options.c_tiltfile inputalifile  = options.c_inputalifile outputfile = options.c_outputfile progname = 'trimvol '

if XCenter == 0:     parser.print_help()     exit()

#read ccp4 file header import read_ccp4 nc,nr,ns = read_ccp4.read_CCP4_header(inputalifile) XTilt   = nc/2#X coordinate for tilt axis ZTilt   = Thick/2 # Z coordinate for tilt axis

print nc , nr , ns

#import docfiles,math import math doc=open(docfile,'r')

 [[DistanceTiltAxisX]] = XCenter -  XTilt [[DistanceTiltAxisZ]] = ZCenter -  ZTilt

#leave xmipp files here dirname = './IMG'  import shutil if os.path.exists(dirname):        shutil.rmtree(dirname) if not os.path.exists(dirname):     os.mkdir(dirname)

#print [[DistanceTiltAxis]] i=0 text_out =" ; Headerinfo columns: rot (1) , tilt (2), psi (3), Xoff (4), Yoff (5), Weight (6), Mirror (7)\n" out_doc = open(outputfile + '.doc','w')

for line in doc:     line = line.strip()     mylist = (line.split())     print  mylist[0] [[TiltAngleD]] = float (mylist[0]) [[TiltAngle]] = math.radians(float (mylist[0]))     #print [[TiltAngle]]     print "TiltAngle ", TiltAngle

 [[DistanceTiltAxisTiltTomographsX]] =  int (math.floor ( (DistanceTiltAxisX * math.cos(TiltAngle))+0.5))     DistanceTiltAxisTiltTomographsZ =  int (math.floor ( (DistanceTiltAxisZ * math.sin(TiltAngle))+0.5))

#print center relationship of [[DistanceTiltAxisTiltTomographs]] [[DistanceTiltAxisTiltTomographs]]  =  (DistanceTiltAxisTiltTomographsX + DistanceTiltAxisTiltTomographsZ)

#print center relationship of [[DistanceTiltAxisTiltTomographsFromCorner]] [[DistanceTiltAxisTiltTomographsFromCorner]] = (XTilt+DistanceTiltAxisTiltTomographs)

#to avoid trimvol errors

wXUpperCorner=DistanceTiltAxisTiltTomographsFromCorner - int(USize/2)     wXLowerCorner=DistanceTiltAxisTiltTomographsFromCorner + int(USize/2)     if wXUpperCorner < 1:        wXUpperCorner = 1        wXLowerCorner = USize     elif wXLowerCorner > nc:        wXUpperCorner = nc -   USize +1        wXLowerCorner = nc

wYUpperCorner=YCenter - int(USize/2)     wYLowerCorner=YCenter + int(USize/2)     if wYUpperCorner < 1:        wYUpperCorner = 1        wYLowerCorner = USize     elif wYLowerCorner > nr:        wYUpperCorner = nr -   USize +1        wYLowerCorner = nr

command` progname + ' -x ' + str(wXUpperCorner) + ',' + str(wXLowerCorner) + ' '     command +`	      ' -y ' + str(wYUpperCorner) + ',' + str(wYLowerCorner) + ' '     command += '-z ' + str(i+1) + ',' + str(i+1) + ' '     command += inputalifile + ' ' + outputfile + '_img' + str("(i+1)) + '.mrc >>' + ' ' + outputfile + '.log'

print  command     os.system(command)

# convert to xmipp format     mrc_filename   = outputfile + '_img' + str("(i+1)) + '.xmp'

command` 'xmipp_convert_spi22ccp4 -i ' +  mrc_filename + ' -o ' + xmipp_filename     print command     os.system(command)     text_out +` ' ; ' + xmipp_filename + '\n'     text_out += str ("(-1. *[[TiltAngleD]])) + '   0.00000    0.00000    0.00000   0.00000    0.00000\n'     i += 1

#create docfile out_doc.write(text_out) out_doc.close() os.system('xmipp_selfile_create "' + dirname + '/*.xmp" > ' + outputfile + '.sel') os.system('xmipp_header_assign -i '  + outputfile + '.doc')

command = 'bcat -slices -output ' + outputfile + '_ali.mrc' + ' ' + outputfile + '_img???.mrc' print command os.system(command)

#clear dir remove_seed = outputfile + '_img???.mrc' print "remove", remove_seed os.system('rm ' + remove_seed)

#leave xmipp files now here dirname2 = './' + outputfile  import shutil if os.path.exists(dirname2):        shutil.rmtree(dirname2) if not os.path.exists(dirname2):     os.mkdir(dirname2)

print "moving files..." os.system('mv ' + outputfile + '.??? '  + outputfile + '/' )    os.system('mv ' + dirname + ' ' + outputfile + '/'  )

--Main.RobertoMarabini - 05 May 2009