def writeVarianceImage(imagicfile, varmrcfile):
        imgdict = readImagic(imagicfile)
        if imgdict is None:
                return
        vararray = imgdict['images'].std(0)
        mrc.write(vararray, varmrcfile)
        return vararray

	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 correlate_template(self):
		'''
		Correlate template that is already created and configured.
		'''
		fromimage = 'original'
		if self.__results[fromimage] is None or self.__results['template'] is None:
			raise RuntimeError('need image %s and template before correlation' % (fromimage,))
		edges = self.__results[fromimage]
		edges = self.maskBlack(edges)
		template = self.__results['template']
		cortype = self.correlation_config['cortype']
		corfilt = self.correlation_config['corfilt']
		if cortype == 'cross':
			cc = correlator.cross_correlate(edges, template)
		elif cortype == 'phase':
			cc = correlator.phase_correlate(edges, template, zero=False)
		else:
			raise RuntimeError('bad correlation type: %s' % (cortype,))

		if corfilt is not None:
			kernel = convolver.gaussian_kernel(*corfilt)
			self.convolver.setKernel(kernel)
			cc = self.convolver.convolve(image=cc)
		self.__update_result('correlation', cc)
		if self.save_mrc:
			mrc.write(cc, 'correlation.mrc')

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 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 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 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 create_template(self):
		fromimage = 'original'
		if self.__results[fromimage] is None:
			raise RuntimeError('need image %s before creating template' % (fromimage,))
		self.configure_template()

		# read template file
		filename = self.template_config['template filename']
		tempim = self.read_hole_template(filename)

		# create template of proper size
		shape = self.__results[fromimage].shape
		center = (0,0)
		filediameter = self.template_config['file diameter']
		diameter = self.template_config['template diameter']
		scale = float(diameter) / filediameter

		im2 = scipy.ndimage.zoom(tempim, scale)
		origshape = im2.shape
		edgevalue = im2[0,0]
		template = edgevalue * numpy.ones(shape, im2.dtype)
		offset = ( (shape[0]-origshape[0])/2, (shape[1]-origshape[1])/2 )
		template[offset[0]:offset[0]+origshape[0], offset[1]:offset[1]+origshape[1]] = im2
		shift = (shape[0]/2, shape[1]/2)
		template = scipy.ndimage.shift(template, shift, mode='wrap')

		template = template.astype(numpy.float32)
		self.__update_result('template', template)
		if self.save_mrc:
			mrc.write(template, 'template.mrc')

        def processAndSaveFFT(self, imgdata, fftpath):
                if os.path.isfile(fftpath):
                        print "FFT file found"
                        if fftpath in self.freqdict.keys():
                                print "Freq found"
                                return False
                        print "Freq not found"
                print "creating FFT file: ", fftpath

                ### downsize and filter leginon image
                if self.params['uncorrected']:
                        imgarray = imagefilter.correctImage(imgdata, params)
                else:
                        imgarray = imgdata['image']

                ### calculate power spectra
                apix = apDatabase.getPixelSize(imgdata)
                fftarray, freq = ctfpower.power(imgarray, apix, mask_radius=0.5, fieldsize=self.params['fieldsize'])
                #fftarray = imagefun.power(fftarray, mask_radius=1)

                fftarray = ndimage.median_filter(fftarray, 2)

                ## preform a rotational average and remove peaks
                rotfftarray = ctftools.rotationalAverage2D(fftarray)
                stdev = rotfftarray.std()
                rotplus = rotfftarray + stdev*4
                fftarray = numpy.where(fftarray > rotplus, rotfftarray, fftarray)

                ### save to jpeg
                self.freqdict[fftpath] = freq
                mrc.write(fftarray, fftpath)

                self.saveFreqFile()

                return True

	def processParticles(self, imgdata, partdatas, shiftdata):
		self.shortname = apDisplay.short(imgdata['filename'])

		### if only selected points along helix,
		### fill in points with helical step
		if self.params['helicalstep']:
			apix = apDatabase.getPixelSize(imgdata)
			partdatas = self.fillWithHelicalStep(partdatas, apix)

		### run batchboxer
		self.boxedpartdatas, self.imgstackfile, self.partmeantree = self.boxParticlesFromImage(imgdata, partdatas, shiftdata)
		if self.boxedpartdatas is None:
			self.stats['lastpeaks'] = 0
			apDisplay.printWarning("no particles were boxed from "+self.shortname+"\n")
			self.badprocess = True
			return None

		self.stats['lastpeaks'] = len(self.boxedpartdatas)

		apDisplay.printMsg("do not break function now otherwise it will corrupt stack")
		#time.sleep(1.0)

		### merge image particles into big stack
		totalpart = self.mergeImageStackIntoBigStack(self.imgstackfile, imgdata)

		### create a stack average every so often
		if self.stats['lastpeaks'] > 0:
			totalPartices = self.existingParticleNumber+self.stats['peaksum']+self.stats['lastpeaks']
			logpeaks = math.log(totalPartices)
			if logpeaks > self.logpeaks:
				self.logpeaks = math.ceil(logpeaks)
				numpeaks = math.ceil(math.exp(self.logpeaks))
				apDisplay.printMsg("writing averaging stack, next average at %d particles"%(numpeaks))
				mrc.write(self.summedParticles/float(totalPartices), "average.mrc")
		return totalpart

	def correlate_template(self):
		fromimage = 'edges'

		if None in (self.__results[fromimage], self.__results['template']):
			raise RuntimeError('need image %s and template before correlation' % (fromimage,))
		edges = self.__results[fromimage]
		template = self.__results['template']
		cortype = self.correlation_config['cortype']
		corfilt = self.correlation_config['corfilt']
		if cortype == 'cross':
			cc = correlator.cross_correlate(edges, template)
		elif cortype == 'phase':
			cc = correlator.phase_correlate(edges, template, zero=False)
		else:
			raise RuntimeError('bad correlation type: %s' % (cortype,))
		cc = numpy.absolute(cc)

		if corfilt is not None:
			kernel = convolver.gaussian_kernel(*corfilt)
			self.edgefinder.setKernel(kernel)
			cc = self.edgefinder.convolve(image=cc)
		#cc = imagefun.zscore(smooth)
		#cc = imagefun.zscore(cc)
		self.__update_result('correlation', cc)
		if self.save_mrc:
			mrc.write(cc, 'correlation.mrc')

def createSingleImage(inputimages, globaloutput, outfilename, outformat, outtext, outxml=False):
	n = len(inputimages)
	tiles = []
	for i,input in enumerate(inputimages):
		print 'inserting %d of %d' % (i+1,n)
		tile = globaloutput.insertImage(input, outtext=outtext)
		tiles.append(tile)
	'''
	if targetinputs:
		for input,target in targetinputs:
			globaloutput.insertImage(input, target)
		print 'inserted %s targets' % (len(globaloutput.targets),)
	for target in globaloutput.targets:
		print 'T', target
		markTarget(globaloutput.image, target)
	'''
	if outtext is not None:
		f = open(outtext, 'w')
		if outxml:
			xmldoc = getMosaicXMLData(tiles)
			xmldoc.writexml(f, "    ", "", "\n", "UTF-8")
		else:
			lines = [ str(tile)+"\n" for tile in tiles ]
			f.writelines(lines)
		f.close()

	if outfilename is not None:
		mrc.write(globaloutput.image, outfilename)

	def create_template(self, ring_list=None,tilt_axis=None,tilt_angle=None):
		'''
		This creates the template image that will be correlated
		with the edge image.  This will fail if there is no
		existing edge image, which is necessary to determine the
		size of the template.
		'''
		fromimage = 'edges'
		if self.__results[fromimage] is None:
			raise RuntimeError('need image %s before creating template' % (fromimage,))
		self.configure_template(ring_list,tilt_axis,tilt_angle)
		shape = self.__results[fromimage].shape
		center = (0,0)
		ring_list = self.template_config['ring_list']
		template = numpy.zeros(shape, numpy.int8)
		tilt_axis = self.template_config['tilt_axis']
		tltangle = self.template_config['tilt_angle']
		tltaxis=(tilt_axis*numpy.pi/180)
		for ring in ring_list:
			temp = self.oval.get(shape, center, ring[0], ring[1],tltangle,tltaxis)
			template = template | temp
		template = template.astype(numpy.float32)
		#template = imagefun.zscore(template)
		self.__update_result('template', template)
		if self.save_mrc:
			mrc.write(template, 'template.mrc')

def writeTemp(imdata, newarray):
                path = imdata['session']['image path']
                filename = imdata['filename'] + '.mrc'
                tmppath = os.path.split(path)[:-1] + ('tmp',)
                tmppath = os.path.join(*tmppath)
                fullname = os.path.join(tmppath, filename)
                mrc.write(newarray, fullname)

def register(image1, image2):
        trim = 8
        image1 = image1[trim:-trim, trim:-trim]
        image2 = image2[trim:-trim, trim:-trim]
        fft1 = scipy.fftpack.fft2(image1)
        fft2 = scipy.fftpack.fft2(image2)

        fft1 = scipy.fftpack.fftshift(fft1, axes=[0])
        fft2 = scipy.fftpack.fftshift(fft2, axes=[0])

        c = int(fft1.shape[0] / 2.0)
        fft1 = fft1[:,:c+1]
        fft2 = fft2[:,:c+1]

        mag1 = numpy.abs(fft1)
        mag2 = numpy.abs(fft2)
        mrc.write(mag1, 'mag1.mrc')
        mrc.write(mag2, 'mag2.mrc')

        center = c,0
        output_shape = c,c
        print 'P1'
        p1 = polar_transform(mag1, output_shape, center)
        #scipy.misc.imsave('p1.jpg', p1)
        mrc.write(p1, 'p1.mrc')
        print 'P2'
        p2 = polar_transform(mag2, output_shape, center)
        #scipy.misc.imsave('p2.jpg', p2)
        mrc.write(p2, 'p2.mrc')

        pc = correlator.phase_correlate(p1, p2, zero=False)
        #pc = correlator.cross_correlate(p1, p2)
        mrc.write(pc, 'pc.mrc')

 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 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)

def emanMrcToStack(partlist):
        apFile.removeStack("emanmrc.hed", warn=False)
        for part in partlist:
                mrc.write(part, "temp.mrc")
                emancmd = "proc2d temp.mrc emanmrc.hed"
                apEMAN.executeEmanCmd(emancmd, verbose=False, showcmd=False)
                apFile.removeFile("temp.mrc")
        return

	def find_edges(self):
		'''
		find edges on the original image
		'''
		if self.__results['original'] is None:
			raise RuntimeError('no original image to find edges on')

		sourceim = self.__results['original']
		filt = self.edges_config['filter']
		sigma = self.edges_config['sigma']
		ab = self.edges_config['abs']
		lpsig = self.edges_config['lpsig']
		edgethresh = self.edges_config['thresh']
		edgesflag = self.edges_config['edges']

		kernel = convolver.gaussian_kernel(lpsig)
		n = len(kernel)
		self.edgefinder.setKernel(kernel)
		smooth = self.edgefinder.convolve(image=sourceim)

		if not edgesflag:
			edges = smooth
		elif filt == 'laplacian3':
			kernel = convolver.laplacian_kernel3
			self.edgefinder.setKernel(kernel)
			edges = self.edgefinder.convolve(image=smooth)
		elif filt == 'laplacian5':
			kernel = convolver.laplacian_kernel5
			self.edgefinder.setKernel(kernel)
			edges = self.edgefinder.convolve(image=smooth)
		elif filt == 'laplacian-gaussian':
			kernel = convolver.laplacian_of_gaussian_kernel(n,sigma)
			self.edgefinder.setKernel(kernel)
			edges = self.edgefinder.convolve(image=smooth)
		elif filt == 'sobel':
			self.edgefinder.setImage(smooth)
			kernel1 = convolver.sobel_row_kernel
			kernel2 = convolver.sobel_col_kernel
			edger = self.edgefinder.convolve(kernel=kernel1)
			edgec = self.edgefinder.convolve(kernel=kernel2)
			edges = numpy.hypot(edger,edgec)
			## zero the image edge effects
			edges[:n] = 0
			edges[:,:n] = 0
			edges[:,-n:] = 0
			edges[-n:] = 0
		else:
			raise RuntimeError('no such filter type: %s' % (filt,))

		if ab and edgesflag:
			edges = numpy.absolute(edges)

		if edgethresh and edgesflag:
			edges = imagefun.threshold(edges, edgethresh)

		self.__update_result('edges', edges)
		if self.save_mrc:
			mrc.write(edges, 'edges.mrc')

 def mergeResults(self):
     allout = "all_out.mrc"
     for proc in range(self.params["nproc"]):
         infilepath = "proc%03d/out.mrc" % (proc)
         a = mrc.read(infilepath)
         if proc == 0:
             mrc.write(a, allout)
         else:
             mrc.append(a, allout)
     return allout