Commit dde0226c authored by Martin Schorb's avatar Martin Schorb
Browse files

mapfix application2

parent 01276be8
......@@ -41,16 +41,28 @@ target_map = 'refmap'
import os
import os.path
import numpy
from operator import itemgetter
#import matplotlib.pyplot as plt
#import tifffile as tiff
import mrcfile as mrc
import pyEM as em
# ====================================================================================
#%%
# start script
......@@ -60,12 +72,9 @@ navlines = em.loadtext(navname)
if targetitem == []:
raise Exception('ERROR! No reference map with label "'+target_map+'" specified!')
target_merge = em.mergemap(targetitem)
targetheader = target_merge['mergeheader']
t_mat = em.map_matrix(targetitem)
newnavf = navname[:-4] + '_automaps.nav'
#nnf = open(newnavf,'w')
#nnf.write("%s\n" % navlines[0])
......@@ -83,6 +92,7 @@ non_acq = [x for x in allitems if x not in acq]
non_acq.remove(targetitem)
maps = {}
maps['mapnav'] = []
newmapid = em.newID(allitems,10000)
......@@ -91,131 +101,24 @@ ntotal = len(acq)
newnav = list()
for idx,acq_item in enumerate(acq):
print('Processing navitem '+ str(idx+1) + '/' + str(ntotal) + ' (%2.0f%% done)' %(idx*100/ntotal))
newmapid = em.newID(allitems,newmapid + 1)
mapitem = em.realign_map(acq_item,allitems)
itemid = mapitem['# Item']
if not itemid in maps.keys():
maps[itemid] = em.mergemap(mapitem)
groupid = em.newID(allitems,999000000+int(mapitem['MapID'][0][-6:]))
non_acq.remove(mapitem)
# NoRealign
mapitem['Color'] = '5'
newnav.append(mapitem)
# combine rotation matrices
map_mat = maps[itemid]['matrix']
maptf = (numpy.linalg.inv(map_mat) * t_mat).T
xval = float(acq_item['StageXYZ'][0]) #(float(acq_item['PtsX'][0]))
yval = float(acq_item['StageXYZ'][1]) #(float(acq_item['PtsY'][0]))
pt = numpy.array([xval,yval])
# calculate the pixel coordinates
pt_px1 = em.get_pixel(acq_item,maps[itemid])
px_scale = targetheader['pixelsize'] /( maps[itemid]['mapheader']['pixelsize'] )
imsz1 = numpy.array([targetheader['ysize'],targetheader['xsize']])
im = numpy.array(maps[itemid]['im'])
im2, p2 = em.map_extract(im,pt_px1,pt_px1,px_scale,imsz1,maptf)
im2size = im2.shape
if min(im2.shape)<200:
print('Warning! Item ' + acq_item['# Item'] + ' is not within the map frame. Ignoring it')
else:
# pad item numbers to 4 digits
if acq_item['# Item'].isdigit(): acq_item['# Item'] = acq_item['# Item'].zfill(4)
imfile = 'virt_map_' + acq_item['# Item'] + '.mrc'
if os.path.exists(imfile): os.remove(imfile)
# tiff.imsave(imfile,im2)
im2 = numpy.rot90(im2,3)
with mrc.new(imfile) as mrcf:
mrcf.set_data(im2.T)
mrcf.close()#
cx = im2.shape[0]
cy = im2.shape[1]
a = [[0,0],[cx,0],[cx,cy],[0,cy],[0,0]]
a = numpy.matrix(a) - [cx/2 , cy/2]
t_mat_i = numpy.linalg.inv(t_mat)
c1 = a*t_mat_i.T + pt
cnx = numpy.array(numpy.transpose(c1[:,1]))
cnx = numpy.array2string(cnx,separator=' ')
cnx = cnx[2:-2]
cny = numpy.array(numpy.transpose(c1[:,0]))
cny = " ".join(list(map(str,cny)))
cny = cny[1:-2]
# fill navigator
acq_item['Acquire'] = '0'
# NoRealign
# acq_item['Color'] = '5'
outnav.append(acq_item)
newnavitem = dict(targetitem)
newnavitem['MapFile'] = [imfile]
newnavitem['StageXYZ'] = acq_item['StageXYZ']
newnavitem['RawStageXY'] = acq_item['StageXYZ'][0:2]
newnavitem['PtsY'] = cnx.split()
newnavitem['PtsX'] = cny.split()
newnavitem['NumPts'] = ['1']
newnavitem['Note'] = newnavitem['MapFile']
newnavitem['MapID'] = [str(newmapid)]
newnavitem['Acquire'] = ['1']
newnavitem['MapSection'] = ['0']
newnavitem['SamePosId'] = acq_item['MapID']
newnavitem['GroupID'] = [str(groupid)]
newnavitem['MapWidthHeight'] = [str(im2size[1]),str(im2size[0])]
newnavitem['ImageType'] = ['0']
newnavitem['MapMinMaxScale'] = [str(numpy.min(im2)),str(numpy.max(im2))]
newnavitem['NumPts'] = ['5']
newnavitem['# Item'] = 'map_' + acq_item['# Item']
outnav.append(newnavitem)
newnavitem, maps, item = em.virt_map_at_point(acq_item,idx,maps,allitems,targetitem,targetheader,outnav,outformat='mrc',numtiles=1)
outnav.append(item)
outnav.append(newnavitem)
outnav.sort(key=itemgetter('# Item'))
outnav.sort(key=itemgetter('# Item'))
#for nitem in non_acq:
# newnav.append(nitem)
for nitem in outnav:
newnav.append(nitem)
on1=em.ordernav(outnav, delim='_')
finalnav = maps['mapnav'].copy()
finalnav.extend(on1.copy())
em.write_navfile(newnavf,newnav,xml=False)
em.write_navfile(newnavf,finalnav,xml=False)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment