def shiftToCenter(infile,shiftfile,isEMAN=False):
	'''
	EMAN defines the rotation origin differently from other packages.
	Therefore, it needs to be recenterred according to the package
	after using EMAN proc3d rotation functions.
	'''
	# center of rotation for eman is not at length/2.
	if isEMAN:
		formatoffset = getEmanCenter()
		prefix = ''
	else:
		formatoffset = (0,0,0)
		prefix = 'non-'

	apDisplay.printMsg('Shifting map center for %sEMAN usage' % (prefix,))
	# Find center of mass of the density map
	a = mrc.read(infile)
	t = a.mean()+2*a.std()
	numpy.putmask(a,a>=t,t)
	numpy.putmask(a,a<t,0)
	center = ndimage.center_of_mass(a)
	offset = (center[0]+formatoffset[0]-a.shape[0]/2,center[1]+formatoffset[1]-a.shape[1]/2,center[2]+formatoffset[2]-a.shape[2]/2)
	offset = (-offset[0],-offset[1],-offset[2])
	apDisplay.printMsg('Shifting map center by (x,y,z)=(%.2f,%.2f,%.2f)' % (offset[2],offset[1],offset[0]))
	# shift the map
	a = mrc.read(infile)
	a = ndimage.interpolation.shift(a,offset)
	mrc.write(a,shiftfile)
	h = mrc.readHeaderFromFile(infile)
	mrc.update_file_header(shiftfile,h)

def writeMrcStack(path, stackname, mrc_files, binning=1):
        apDisplay.printMsg("Writing MRC stack file... ")
        stackname = os.path.join(path, stackname)
        im = mrc.read(mrc_files[0])
        image = imagefun.bin(im, binning)
        mrc.write(image,stackname)
        del mrc_files[0]
        for mrcfile in mrc_files:
                im = mrc.read(mrcfile)
                image = imagefun.bin(im, binning)
                mrc.append(image, stackname)

def testMRCImages():
        file1,file2 = sys.argv[1:3]
        print 'reading MRCs'
        image1 = mrc.read(file1)
        image2 = mrc.read(file2)
        image1 = imagefun.bin(image1, 4)
        image2 = imagefun.bin(image2, 4)
        print 'register...'
        #result = register(image1, image2, range()
        #result = register2(image1, image2, range(86,95))
        result = register2(image1, image2, range(90,91))

 def makeStack(self, tiltseries, mrc_files):
         stackname = self.getFilename(tiltseries) + '.st'
         stackname = os.path.join(self.settings['path'], stackname)
         im = mrc.read(mrc_files[0])
         image = imagefun.bin(im, int(self.settings['binning']))
         mrc.write(image,stackname)
         #shutil.copy(mrc_files[0], stackname)
         del mrc_files[0]
         for mrcfile in mrc_files:
                 im = mrc.read(mrcfile)
                 image = imagefun.bin(im, int(self.settings['binning']))
                 mrc.append(image, stackname)

 def readUploadInfo(self,info=None):
         if info is None:
                 # example
                 info = ['test.mrc','2e-10','1','1','50000','-2e-6','120000']
         self.logger.info('reading image info')
         try:
                 self.uploadedInfo = {}
                 self.uploadedInfo['original filepath'] = os.path.abspath(info[0])
                 self.uploadedInfo['unbinned pixelsize'] = float(info[1])
                 self.uploadedInfo['binning'] = {'x':int(info[2]),'y':int(info[3])}
                 self.uploadedInfo['magnification'] = int(info[4])
                 self.uploadedInfo['defocus'] = float(info[5])
                 self.uploadedInfo['high tension'] = int(info[6])
                 if len(info) > 7:
                         self.uploadedInfo['stage a'] = float(info[7])*3.14159/180.0
                 # add other items in the dictionary and set to instrument in the function
                 # setInfoToInstrument if needed
         except:
                 #self.logger.exception('Bad batch file parameters')
                 raise
         try:
                 self.uploadedInfo['image'] = mrc.read(self.uploadedInfo['original filepath'])
         except IOError, e:
                 self.logger.exception('File %s not available for upload' % self.uploadedInfo['original filepath'])
                 raise

 def averageStack(self, stackfile):
     avgfile = "avg.mrc"
     a = mrc.read(stackfile)
     a = np.sum(a, axis=0)
     a = (a - a.min()) / (a.max() - a.min())
     mrc.write(a, avgfile)
     return avgfile

 def load(self, filename=None):
         if filename is None:
                 filename = self.filename.get()
         try:
                 self.pluginpipeline.process(Image(mrc.read(filename)))
         except IOError:
                 self.logger.error('Load file "%s" failed' % filename)

	def start(self):
		#def runCTFdisplayTools(imgdata, ctfvalues, opimagedir, fftpath=None, fftfreq=None):
		### RUN CTF DISPLAY TOOLS
		imgdata = {
			'filename': os.path.abspath(self.params['imagefile']),
			'image': mrc.read(self.params['imagefile']),
		}
		ctfdata = {
			'volts': self.params['kv']*1e3,
			'cs': self.params['cs'],
			'apix': self.params['apix'],
			'defocus1': self.params['defocus1']*1e-6,
			'defocus2': self.params['defocus2']*1e-6,
			'angle_astigmatism': self.params['astigangle'],
			'amplitude_contrast': self.params['ampcontrast'],
		}
		a = ctfdisplay.CtfDisplay()
		a.debug = self.params['debug']
		ctfdisplay.ctftools.debug = self.params['debug']
		ctfdisplaydict = a.CTFpowerspec(imgdata, ctfdata, None, None, True)
		if ctfdisplaydict is None:
			raise

		ctfdata['confidence_30_10'] = ctfdisplaydict['conf3010']
		ctfdata['confidence_5_peak'] = ctfdisplaydict['conf5peak']
		ctfdata['overfocus_conf_30_10'] = ctfdisplaydict['overconf3010']
		ctfdata['overfocus_conf_5_peak'] = ctfdisplaydict['overconf5peak']
		ctfdata['resolution_80_percent'] = ctfdisplaydict['res80']
		ctfdata['resolution_50_percent'] = ctfdisplaydict['res50']
		### override the confidence
		ctfdata['confidence'] = max(ctfdisplaydict['conf5peak'], ctfdisplaydict['conf3010'])

		return ctfdata

	def processCTF(self, event):  # wxGlade: MyFrame.<event_handler>
		if self.checkCTFvalues() is False:
			event.Skip()
			return
		self.statbar.PushStatusText("Processing please wait...", 0)
		self.processButton.SetBackgroundColour(wx.Colour(255, 191, 191))
		self.Update()
		#self.statbar.SetBackgroundColour(wx.Colour(255, 191, 191))
		imgdata = {
			'filename': self.fullimagepath,
			'image': mrc.read(self.fullimagepath),
		}
		ctfdata = {
			'volts': self.voltValue.GetValue()*1e3,
			'cs': self.csValue.GetValue(),
			'apix': self.pixelSizeValue.GetValue(),
			'defocus1': self.defoc1Value.GetValue()*1e-6,
			'defocus2': self.defoc2Value.GetValue()*1e-6,
			'angle_astigmatism': self.angleValue.GetValue(),
			'amplitude_contrast': self.ampConValue.GetValue(),
		}
		a = ctfdisplay.CtfDisplay()
		ctfdisplaydict = a.CTFpowerspec(imgdata, ctfdata, None, None, True)
		self.statbar.PushStatusText("Finished", 0)
		self.processButton.SetBackgroundColour(wx.Colour(191, 255, 191))
		self.Update()
		
		self.showImages(ctfdisplaydict)

		event.Skip()

 def makeDenoisedParticleJPGFrames(self,partdatas,shiftdata,shortname):
         '''
         Denoise the boxed particles and then read in the resulting frame stack of the particle for saving as movie frames.
         '''
         if not partdatas:
                 return
         if not self.is_dd_stack:
                 apDisplay.printError('Denoising works only with ddstack for now')
         imgdata = partdatas[0]['image']
         ddstackdir,stackfile = self.getDDStackDirFile(imgdata)
         framestacks = []
         for p,partdata in enumerate(partdatas):
                 # denoise within the particle box
                 col_start,row_start = apBoxer.getBoxStartPosition(imgdata,self.half_box,partdata, shiftdata)
                 row_end = row_start + self.boxsize
                 col_end = col_start + self.boxsize
                 roi = ((row_start,row_end),(col_start,col_end))
                 paramstr = self.denoise.setupKSVDdenoise(self.frameavg,self.firstframe,self.nframe,roi)
                 apDisplay.printMsg('denoise param string: %s' % paramstr)
                 self.denoise.makeDenoisedStack(ddstackdir, stackfile)
                 outputstackfile = '%s_%s.mrc' % (stackfile[:-4], paramstr)
                 framestacks.append(outputstackfile)
         for i,start_frame in enumerate(range(self.firstframe,self.nframe-self.frameavg-self.firstframe+1,self.framestep)):
                 for p,partdata in enumerate(partdatas):
                         array = mrc.read(os.path.join(self.params['rundir'],'results/mrc',framestacks[p]),i)
                         # bin is not used for now
                         movieframe_number = start_frame
                         self.saveFrameFromArray(array,shortname,p,movieframe_number)

def getImageFiles(imgtree, rawdir, link, copy):
	#This function should replace linkImageFiles in all calls (i.e. in tomoaligner.py and everywhere else)
	filenamelist = []
	newimgtree=[]
	for imagedata in imgtree:
		#set up names
		imgpath=imagedata['session']['image path']
		presetname=imagedata['preset']['name']
		imgprefix=presetname+imagedata['filename'].split(presetname)[-1]
		imgname=imgprefix+'.mrc'
		filenamelist.append(imgprefix)
		destpath = os.path.join(rawdir,imgname)
		newimgtree.append(destpath)
		imgfullpath = os.path.join(imgpath,imagedata['filename']+'.mrc')

		if link == "True":
			#create symlinks to files
			if os.path.islink(destpath):
				os.remove(destpath)
			if not os.path.isfile(destpath):
				os.symlink(imgfullpath,destpath)
		elif copy == "True":
			shutil.copy(imgfullpath,destpath)	
			
			#Y-flip raw images, normalize them, and convert them to float32 because Protomo
			image=mrc.read(destpath)
			image=numpy.flipud(image)
			image=imagenorm.normStdev(image)
			image=numpy.float32(image)
			mrc.write(image,destpath)
		#else: just return values
	return filenamelist, newimgtree

	def applyEnvelope(self, inimage, outimage, scaleFactor=1, msg=False):
		"""
		input path to image and envelope, output amplitude-adjusted image
		"""

		if msg is True:
			apDisplay.printColor("now applying envelope function to: "+inimage, "cyan")

		if self.envamp is None:
			self.prepareEnvelope(scaleFactor)

		### read image
		im = mrc.read(inimage)

		### fourier transform
		imfft = self.real_fft2d(im)

		### mutliply real envelope function by image fft
		newfft = self.envamp * imfft

		### inverse transform
		newimg = self.inverse_real_fft2d(newfft)

		### normalize between 0 and 1
		newimg = (newimg-newimg.mean()) / newimg.std()

		### save image
		mrc.write(newimg, outimage)

		### workaround for now
		time.sleep(0.1)

		return

	def medianVolume(self):
		volpath = os.path.join(self.params['rundir'], "volumes/*a.mrc")
		mrcfiles = glob.glob(volpath)
		volumes = []
		for filename in mrcfiles:
			if os.path.isfile(filename):
				vol = mrc.read(filename)
				print filename, vol.shape
				volumes.append(vol)
		volarray = numpy.asarray(volumes, dtype=numpy.float32)
		try:
			medarray = numpy.median(volarray, axis=0)
		except:
			medarray = numpy.median(volarray)
		medfile = os.path.join(self.params['rundir'], "volumes/medianVolume.mrc")
		print medfile, medarray.shape
		mrc.write(medarray, medfile)

		apix = apStack.getStackPixelSizeFromStackId(self.params['stackid'])
		sessiondata = apStack.getSessionDataFromStackId(self.params['stackid'])

		uploadcmd = ( ("uploadModel.py --projectid=%d --session=%s --file=%s "
				+"--apix=%.3f --sym=%s --name=satmedian-recon%d.mrc --res=30 --description='%s %d'")
			%(self.params['projectid'], sessiondata['name'], medfile, 
				apix, self.params['symmname'], self.params['reconid'],
				"SAT selected median volume for recon", self.params['reconid'], ) )
		apDisplay.printColor(uploadcmd, "purple")
		f = open("upload.sh", "w")
		f.write(uploadcmd+"\n")
		f.close()

def boxerRotate(imgfile, parttree, outstack, boxsize):
        """
        boxes the particles with expanded size,
        applies a rotation to particle,
        reduces boxsize to requested size,
        and saves them to a imagic file
        """
        # size needed is sqrt(2)*boxsize, using 1.5 to be extra safe
        bigboxsize = int(math.ceil(1.5*boxsize))
        imgarray = mrc.read(imgfile)
        bigboxedparticles = boxerMemory(imgarray, parttree, bigboxsize)
        
        boxedparticles = []
        boxshape = (boxsize,boxsize)
        apDisplay.printMsg("Rotating particles...")
        for i in range(len(bigboxedparticles)):
                if i % 10 == 0:
                        sys.stderr.write(".")
                bigboxpart = bigboxedparticles[i]
                partdict = parttree[i]
                ### add 90 degrees because database angle is from x-axis not y-axis
                angle = partdict['angle']+90.0
                rotatepart = ndimage.rotate(bigboxpart, angle=angle, reshape=False, order=1)
                boxpart = imagefilter.frame_cut(rotatepart, boxshape)
                boxedparticles.append(boxpart)
        sys.stderr.write("done\n")
        apImagicFile.writeImagic(boxedparticles, outstack)
        return True

def boxMaskStack(bmstackf, partdatas, box, xmask, ymask, falloff, imask=None, norotate=False):
        from appionlib.apSpider import operations
        from appionlib import apEMAN
        import os

        # create blank image for mask using SPIDER
        maskfile = "boxmask.spi"
        operations.createBoxMask(maskfile,box,xmask,ymask,falloff,imask)

        # convert mask to MRC
        apEMAN.executeEmanCmd("proc2d boxmask.spi boxmask.mrc",verbose=False,showcmd=False)
        os.remove("boxmask.spi")

        maskarray = mrc.read("boxmask.mrc")

        # box particles
        maskedparts = []
        for i in range(len(partdatas)):
                if norotate is True:
                        rotatemask = maskarray
                else:
                        angle = (-partdatas[i]['angle'])-90
                        rotatemask = ndimage.rotate(maskarray, angle=angle, reshape=False, order=1)
                maskedparts.append(rotatemask)

        # write to stack
        apImagicFile.writeImagic(maskedparts, bmstackf)
        os.remove("boxmask.mrc")
        return bmstackf

        def targetTestImage(self):
                usercheck = self.settings['user check']
                self.settings['user check'] = False
                filename = self.settings['test image']
                
                try:
                        image2 = mrc.read(filename)             
                        image = numpy.asarray(image2,dtype=numpy.float)
                        
                except:
                        self.logger.error('Failed to load test image')
                        raise
                        return
                self.setImage(image, 'Image')

                if self.handle is None:
                        self.handle = pymat.open()
                pymat.put(self.handle, 'image', image)

                imdata_id = 0

                pymat.put(self.handle, 'image_id',imdata_id)
        
                self.matlabFindTargets()

                self.settings['user check'] = usercheck

	def readFile(self, oldmrcfile):
		apDisplay.printMsg('Reading %s into memory' % oldmrcfile)
		imagearray = mrc.read(oldmrcfile)
		# invert image density
		if self.params['invert'] is True:
			imagearray *= -1.0
		return imagearray

def makeProjection(filename, xsize=512):
    mrcpath = filename
    dirpath = os.path.dirname(mrcpath)
    apDisplay.printMsg("Reading 3D recon %s" % mrcpath)
    array = mrc.read(mrcpath)
    shape = array.shape
    xsize = min(xsize, shape[2])
    # default for full tomogram is XZY
    if shape[0] > shape[1]:
        renders = {
            "a": {"axis": 0, "axisname": "z"},
            "b": {"axis": 1, "axisname": "y"},
            "c": {"axis": 2, "axisname": "x"},
        }
    else:
        renders = {
            "a": {"axis": 1, "axisname": "y"},
            "b": {"axis": 0, "axisname": "z"},
            "c": {"axis": 2, "axisname": "x"},
        }
    keys = renders.keys()
    keys.sort()
    for key in keys:
        apDisplay.printMsg("project to axis %s" % renders[key]["axisname"])
        pictpath = os.path.join(dirpath, "projection" + key)
        axis = renders[key]["axis"]
        slice = numpy.sum(array[:, :, :], axis=axis) / (shape[axis])
        mrc.write(slice, pictpath + ".mrc")
        # adjust and shrink each image
        array2jpg(pictpath, slice, size=xsize)

        def start(self):
                #get aligned stack id
                aligndata = appiondata.ApAlignStackData.direct_query(self.params['alignstackid'])
                xSizeVoxel = aligndata['boxsize']
                #get averaged image
                avgmrc = os.path.join(aligndata['path']['path'], aligndata["avgmrcfile"])
                avg = mrc.read(avgmrc)
                tmpSpiderFile="average.xmp"
                spider.write(avg, tmpSpiderFile)
                self.runFindCenter(tmpSpiderFile,xSizeVoxel)
                #get aligned stack
                alignStack = os.path.join(aligndata['path']['path'], aligndata["imagicfile"])
                tempFileNameforSpectra=self.runMakeSpectra(alignStack)

                #kerdensom will work with spectra output
                self.runKerdenSOM(tempFileNameforSpectra)
                self.createMontageInMemory()
                self.insertRotKerDenSOM()

                #apFile.removeFile(outdata)
                apFile.removeFilePattern("*.cod")
                apFile.removeFilePattern("*.err")
                apFile.removeFilePattern("*.his")
                apFile.removeFilePattern("*.inf")
                apFile.removeFilePattern("*.vs")
                apFile.removeFilePattern("*.xmp")
                apFile.removeFile(tempFileNameforSpectra)

def filterAndChimera(density, res=30, apix=None, box=None, chimtype='snapshot',
		contour=None, zoom=1.0, sym='c1', color=None, silhouette=True, mass=None):
	"""
	filter volume and then create a few snapshots for viewing on the web
	"""
	if isValidVolume(density) is False:
		apDisplay.printError("Volume file %s is not valid"%(density))
	if box is None:
		boxdims = apFile.getBoxSize(density)
		box = boxdims[0]
	### if eotest failed, filter to 30
	if not res or str(res) == 'nan':
		res = 30
	### low pass filter the volume to 60% of reported res
	tmpf = os.path.abspath(density+'.tmp.mrc')
	density = os.path.abspath(density)
	filtres = 0.6*res
	shrinkby = 1
	if box is not None and box > 250:
		shrinkby = int(math.ceil(box/160.0))
		if box % (2*shrinkby) == 0:
			### box is divisible by shrink by
			lpcmd = ('proc3d %s %s apix=%.3f tlp=%.2f shrink=%d origin=0,0,0 norm=0,1'
				% (density, tmpf, apix, filtres, shrinkby))
		else:
			### box not divisible by shrink by, need a clip
			clip = math.floor(box/shrinkby/2.0)*2*shrinkby
			lpcmd = ('proc3d %s %s apix=%.3f tlp=%.2f shrink=%d origin=0,0,0 norm=0,1 clip=%d,%d,%d'
				% (density, tmpf, apix, filtres, shrinkby, clip, clip, clip))
	else:
		lpcmd = ('proc3d %s %s apix=%.3f tlp=%.2f origin=0,0,0 norm=0,1'
			% (density, tmpf, apix, filtres))
	apDisplay.printMsg("Low pass filtering model for images")
	proc = subprocess.Popen(lpcmd, shell=True)
	proc.wait()

	### flatten solvent
	vol = mrc.read(tmpf)
	numpy.where(vol < 0, 0.0, vol)
	mrc.write(vol, tmpf)
	del vol

	### contour volume to mass
	if mass is not None:
		setVolumeMass(tmpf, apix*shrinkby, mass)
		contour = 1.0

	### set pixelsize and origin
	recmd = "proc3d %s %s apix=%.3f origin=0,0,0"%(tmpf, tmpf, apix)
	proc = subprocess.Popen(recmd, shell=True)
	proc.wait()

	### render images
	renderSlice(density, box=box, tmpfile=tmpf, sym=sym)
	if chimtype != 'snapshot':
		renderAnimation(tmpf, contour, zoom, sym, color, silhouette, name=density)
	elif chimtype != 'animate':
		renderSnapshots(tmpf, contour, zoom, sym, color, silhouette, name=density)
	apFile.removeFile(tmpf)