Python EMAN類代碼示例

def main():
    (options, args) = parse_command_line()

    LOGbegin(sys.argv)

    map1 = EMAN.EMData()
    map1.readImage(args[0])

    map2 = EMAN.EMData()
    map2.readImage(args[1])
    map2t = map2.copy()

    print "Start searching the scale factor in range [%g, %g]" % (options.scale_min, options.scale_max)
    scale = brent(matching_score, brack=[options.scale_min, options.scale_max], args=(map1, map2t))
    print "Final result: scale = %g" % (scale)
    if options.save_map or options.save_fsc:
        map2.setTAlign(0, 0, 0)
        map2.setRAlign(0, 0, 0)
        map2.rotateAndTranslate(scale)
        if options.save_map:
            map2.writeImage(options.save_map)
            print "Final transformed map2 is saved to %s" % (options.save_map)
        if options.save_fsc:
            fsc = map1.fsc(map2)
            if options.apix:
                EMAN.save_data(0, 1.0 / (options.apix * 2.0 * len(fsc)), fsc, len(fsc), options.save_fsc)
            else:
                EMAN.save_data(0, 1.0, fsc, len(fsc), options.save_fsc)
            print "Final Fourier Shell Correlation curve is saved to %s" % (options.save_fsc)

    LOGend()

def reconstruction_by_make3d(imagefile, mapfile, options, mapfile2e, mapfile2o, fscfile):
	cmd_prefix = EMAN.pre_mpirun(mpilib = options.mpi, mpi_nodefile = options.mpi_nodefile, single_job_per_node = 1)
	cmd_3drec = "%s `which make3d` %s hard=90 lowmem sym=%s mode=2 apix=%g pad=%d usecenter iter=1" % \
			(cmd_prefix, imagefile, options.endSym, options.apix, options.pad)
	if options.sffile: 
		cmd_3drec += " wiener3d=%s" % ( options.sffile )
		if not options.phasecorrected: cmd_3drec += " flipphase"
	cmd = "%s out=%s" % (cmd_3drec, mapfile)
	if options.eotest and mapfile2e and mapfile2o:
		cmd += "\n%s subset=0/2 out=%s" % (cmd_3drec, mapfile2e)
		cmd += "\n%s subset=1/2 out=%s" % (cmd_3drec, mapfile2o)
	print cmd
	if options.mpi:
		os.system(cmd)
	else:
		runpar_file = "runpar.make3d"
		runparfp = open(runpar_file,"w")
		runparfp.write("%s\n" % (cmd))
		runparfp.close()
		cmd = "runpar proc=%d,%d file=%s nofs" % (options.cpus, options.cpus, runpar_file)
		print cmd
		os.system(cmd)
	if not (os.path.exists(mapfile) and os.path.getsize(mapfile)):
		print "ERROR: 3D map \"%s\" is not generated properly" % (mapfile)
		sys.exit()
	if mapfile2e and mapfile2o and os.path.exists(mapfile2e) and os.path.exists(mapfile2o):
		cmd = "proc3d %s %s fsc=%s" % (mapfile2e,mapfile2o,fscfile)
		print cmd
		os.system(cmd)
	EMAN.post_mpirun(mpilib = options.mpi)

def main():
	
	(options, input_imagefiles, output_lstfile) = parse_command_line()
	
	Particle.sym = options.sym
	Particle.max_ang_diff = options.max_ang_diff * math.pi/180.0
	Particle.max_cen_diff = options.max_cen_diff
		
	image_lists = []
	
	for f in input_imagefiles:	# build the image sets
		tmp_image_list = EMAN.image2list(f)
		
		# first check and remove duplicate images
		seen = sets.Set()
		tmp_image_list2 = []
		for i, img in enumerate(tmp_image_list): 
			imgid = "%s-%d" % (os.path.abspath(img[0]), img[1])
			if imgid in seen:
				print "%s: particle %d/%d (%s %d) is a duplicate image, ignored" % (f, i, len(tmp_image_list), img[0], img[1])
			else:
				seen.add(imgid)
				tmp_image_list2.append(img)
		
		image_list = sets.Set()
		for i, img in enumerate(tmp_image_list2): 
			p = Particle(img)
			image_list.add(p)
		print "%s: %d images" % (f, len(image_list))
		image_lists.append( image_list )

	all_images = image_lists[0]
	print "Begining with %d images in %s" % (len(all_images), input_imagefiles[0])
	for i in range(1,len(image_lists)):
		image_list = image_lists[i]
		if options.mode == "common":	# all_images & imageset
			all_images = intersection(all_images, image_list)
			print "%d common images after processing %d images in %s" % ( len(all_images), len(image_list), input_imagefiles[i] )
		elif options.mode == "union":
			all_images = union(all_images, image_list)	# all_images | imageset
			print "%d images after merging %d images in %s" % ( len(all_images), len(image_list), input_imagefiles[i] )
		elif options.mode == "diff":
			all_images = difference(all_images,image_list)	# all_images-imageset
			print "%d different images after processing %d images in %s" % ( len(all_images), len(image_list), input_imagefiles[i] )
		elif options.mode == "symdiff":
			all_images = symmetric_difference(all_images,image_list)	# all_images ^ imageset
			print "%d different images after processing %d images in %s" % ( len(all_images), len(image_list), input_imagefiles[i] )
	
	all_images=list(all_images)
	all_images.sort()
	
	all_images_output_list = [i.image for i in  all_images]
			
	if len(all_images_output_list):
		print "%d images saved to %s" % ( len(all_images_output_list), output_lstfile )
		EMAN.imagelist2lstfile(all_images_output_list, output_lstfile)
	else:
		print "No image left after the image sets operation"

def random_sample(num_maps, sf, particles):
	#subset_num = int (len(particles)*fraction)
	num_part = int(len(particles))
	num_proj = EMAN.fileCount("proj.img")[0]
	
	for i in range(num_maps):
		#print "Model %d" % (i)
		#random.shuffle(particles)
		directory = "model_%d" % (i)
		os.mkdir(directory)
		os.chdir(directory)
		for files in ["start.hed","proj.hed"]:
		#up_dir_file = "../"+file
			cmd="lstcat.py "+files+" ../"+files
			os.system(cmd)
			os.link(files, files[:-3]+"img")
		os.link("../"+sf,sf)
		for c in range(num_proj):
			clsfile="cls%04d.lst" % (c)
			fp = file(clsfile,"w")
			fp.write("#LST\n%d\tproj.hed\n" % (c) )
			fp.close()
		#for j in range(subset_num):
		for j in range(num_part):
			#p = int(particles[j][0])	
			#c = int(particles[j][1])
			partt=random.choice(particles)
			p = partt[0]
			c = int(partt[1])
			#print "Index %d/%d/%d\tParticle %d -> Class %d" % (j, subset_num, len(particles), p, c)
			clsfile="cls%04d.lst" % (c)
			fp = open(clsfile,"a")
			fp.write(p)
			fp.close()
		os.chdir("..")

def numberParticlesInStack(stackname, startnum=0, verbose=True):
        ### new faster methond
        apImagicFile.numberStackFile(stackname, startnum=0)
        return

        # store the particle number in the stack header
        # NOTE!!! CONFORMS TO EMAN CONVENTION, STARTS AT 0!!!
        import EMAN
        t0 = time.time()
        apDisplay.printMsg("saving particle numbers to stack header")
        n=EMAN.fileCount(stackname)[0]
        im=EMAN.EMData()
        i = 0
        back = "\b\b\b\b\b\b\b"
        while i < n:
                j=startnum+i
                im.readImage(stackname,i)
                im.setNImg(j)
                im.writeImage(stackname,i)
                if verbose is True and i%100 == 0:
                        sys.stderr.write(back+back+back+str(j)+" of "+str(n))
                i+=1
        sys.stderr.write("\n")
        apDisplay.printMsg("finished in "+apDisplay.timeString(time.time()-t0))
        return

def parse_command_line():
	usage="Usage: %prog <image filename> <output filename> [options]"
	parser = OptionParser(usage=usage)
	parser.add_option("--mode",metavar="['defocus', 'snr', 'ctfb', 'ctfamp']",dest="mode",type="choice",choices=['defocus', 'snr', 'ctfb', 'ctfamp'], help="ctf mode to plot. default to \"defocus\"", default="defocus")
	parser.add_option("--sf",dest="sffile",type="string",help="structural factor file name", default="")
	parser.add_option("--first",dest="first",type="int",help="first image to use. default to 0", default=0)
	parser.add_option("--last",dest="last",type="int",help="last image to use. default to last image in the file", default=-1)
	
	(options, args) = parser.parse_args()
	if len(args) != 2:
		parser.print_help()
		sys.exit(-1)
	else:
		image = args[0]
		output = args[1]
	
	if options.mode in ['snr']:
		if not options.sffile:
			print "Error: --sf should be specified"
			sys.exit(1)
		elif not os.path.exists(options.sffile):
			print "Error: sffile %s does not exist" % (options.sffile)
			sys.exit(2)
	imgnum = EMAN.fileCount(image)[0]
	if options.last==-1: options.last=imgnum-1
	if options.first>0 or options.last>0:
		if not (0<=options.first<imgnum):
			print "Error: --first=%d is out of correct range [0, %d]" % (options.first, imgnum-1)
			sys.exit(3)
		if not (options.first<=options.last<imgnum):
			print "Error: --last=%d is out of correct range [%d, %d]" % (options.last, options.first, imgnum-1)
			sys.exit(4)

	return (options, image, output)

def alignParticlesInLST(lstfile,outstack):
        import EMAN
        ### create a stack of particles aligned according to an LST file
        images=EMAN.readImages(lstfile,-1,-1,0)
        for i in images:
                i.edgeNormalize()
                i.rotateAndTranslate()
                if i.isFlipped():
                        i.hFlip()
                i.writeImage(outstack,-1)

def alignParticlesInLST(lstfile,outstack):
        import EMAN
        ### create a stack of particles aligned according to an LST file
        apDisplay.printMsg("aligning particles in '%s', saving to stack: %s"%(lstfile,outstack))        
        images=EMAN.readImages(lstfile,-1,-1,0)
        for i in images:
                i.edgeNormalize()
                i.rotateAndTranslate()
                if i.isFlipped():
                        i.hFlip()
                i.writeImage(outstack,-1)

def main():
	if len(sys.argv) != 3:
		print "Usage::  %s <input particles file name> <output particles file name>"%(sys.argv[0])
		sys.exit(-1)
	input = sys.argv[1]
	output = sys.argv[2]
	
	stack=EMAN.readImages(input,-1,-1)
	i = 0
	for image in stack:
		print "Working on image %s"%(i)
		new_image = speckle_filter(image)
		new_image.writeImage(output,i)
		i=i+1

def writeStackParticlesToFile(stackname, filename):
        import EMAN
        # write out the particle numbers from imagic header to a file
        # NOTE!!! CONFORMS TO EMAN CONVENTION, STARTS AT 0!!!
        apDisplay.printMsg("saving list of saved particles to:")
        apDisplay.printMsg(filename)
        f = open(filename,'w')
        n=EMAN.fileCount(stackname)[0]
        im=EMAN.EMData()
        for i in range(n):
                im.readImage(stackname,i)
                f.write(str(im.NImg())+"\n")
        f.close()
        return

def checkStackNumbering(stackname):
        import EMAN
        # check that the numbering is stored in the NImg parameter
        apDisplay.printMsg("checking that original stack is numbered")
        n=EMAN.fileCount(stackname)[0]
        im=EMAN.EMData()
        im.readImage(stackname,n-1)

        # if last particle is not numbered with same value as # of particles,
        # renumber the entire stack
        if n-1 != im.NImg():
                apDisplay.printWarning("Original stack is not numbered! numbering now...")
                # check that the stack exists in a writable directory and is not read-only
                if os.access(stackname, os.W_OK) is True:
                        numberParticlesInStack(stackname, startnum=0, verbose=True)
                else:
                        apDisplay.printWarning("old stack header exists in a READ-only directory and cannot be numbered according to EMAN")
        return

def getClassInfo(classes):
        import EMAN
        # read a classes.*.img file, get # of images
        imgnum, imgtype = EMAN.fileCount(classes)
        img = EMAN.EMData()
        img.readImage(classes, 0, 1)

        # for projection images, get eulers
        projeulers=[]
        for i in range(imgnum):
                img.readImage(classes, i, 1)
                e = img.getEuler()
                alt = e.thetaMRC()*180./math.pi
                az = e.phiMRC()*180./math.pi
                phi = e.omegaMRC()*180./math.pi
                eulers=[alt,az,phi]
                if i%2==0:
                        projeulers.append(eulers)
        return projeulers

def makeClassAverages(lst, outputstack,e,mask):
        import EMAN
        #align images in class
        print "creating class average from",lst,"to",outputstack
        images=EMAN.readImages(lst,-1,-1,0)
        for image in images:
                image.rotateAndTranslate()
                if image.isFlipped():
                        image.hFlip()

        #make class average
        avg=EMAN.EMData()
        avg.makeMedian(images)

        #write class average
        avg.setRAlign(e)
        avg.setNImg(len(images))
        avg.applyMask(mask,0)
        avg.writeImage(outputstack,-1)

def makeClassAverages(lst, outputstack, classdata, params):
    # align images in class
    images = EMAN.readImages(lst, -1, -1, 0)
    for image in images:
        image.rotateAndTranslate()
        if image.isFlipped():
            image.hFlip()

    # make class average
    avg = EMAN.EMData()
    avg.makeMedian(images)

    # write class average
    e = EMAN.Euler()
    alt = classdata["euler1"] * math.pi / 180
    az = classdata["euler2"] * math.pi / 180
    phi = 0.0
    e.setAngle(alt, az, phi)
    avg.setRAlign(e)
    avg.setNImg(len(images))
    avg.applyMask(params["mask"], 0)
    avg.writeImage(outputstack, -1)

def follow_lst_link(image, index, lstfiles, imagefiles):
    if not (
        imagefiles.has_key(image) or lstfiles.has_key(image)
    ):  # it has not seen this image before and will find out
        imgnum, imgtype = EMAN.fileCount(image)
        if imgtype == "lst":
            tmplst = []
            fp = open(image, "r")
            lines = fp.readlines()
            fp.close()
            if lines[0].startswith("#LST"):
                line0 = 1
            elif lines[0].startswith("#LSX"):
                line0 = 3

            for l in lines[line0:]:
                tokens = l.split()
                tmplst.append((tokens[1], int(tokens[0])))
            lstfiles[image] = tmplst
        else:
            imagefiles[image] = imgnum

    if imagefiles.has_key(image):  # it has already known this image is a binary format image
        if index >= imagefiles[image]:
            print "Error: image file %s only has %d images, not enough for %d" % (image, imagefiles[image], index)
            sys.exit(-1)
        else:
            return (image, index)
    elif lstfiles.has_key(image):  # it has already known this image is a list format image
        if index >= len(lstfiles[image]):
            print "Error: lst file %s only has images, not enough for %d" % (image, len(lstfiles[image]), index)
            sys.exit(-1)
        else:
            return follow_lst_link(lstfiles[image][index][0], lstfiles[image][index][1], lstfiles, imagefiles)
    else:
        print "Something is wrong, need to debug this program and the image file"
        sys.exit(-1)