Commit 5176c0be authored by Martin Schorb's avatar Martin Schorb
Browse files

partial merge test routine

parents 7d2615b2 f9e59012
......@@ -2,4 +2,5 @@
__pycache__/
build/
dist/
pyem_tests/
*.egg-info/
......@@ -722,10 +722,10 @@ def mergemap(mapitem,crop=False,black=False,blendmont=True,bigstitch=False):
else:
tilepx = tilepx1
if mdoc_item(mdoclines,'MontSection = 0') == []: #older mdoc file format, created before SerialEM 3.7x
if not '[MontSection = 0]' in mdoclines: #older mdoc file format, created before SerialEM 3.7x
print('Warning: mrc stack without montage information. Assume pixel size is consistent for all sections.')
str1=mdoclines[0]
pixelsize = float(mdoclines[0][str1.find('=')+1:])
pixelsize = float(mdoclines[0][str1.find('=')+1:])/ 10000 # in um
else:
mont_item = mdoc_item(mdoclines,'MontSection = '+str(mapsection))
pixelsize = float(mont_item['PixelSpacing'][0])/ 10000 # in um
......@@ -1134,46 +1134,90 @@ def map_extract(im,c,p,px_scale,t_size,mat,int8=False):
if (limitsize==t_size).all():
im3=im2.copy()
shift = [0,0]
else:
im3=im2.copy()
shift = [0,0]
else:
# determine limitation of image by the borders of rotated crop
rotmat=M2[:2,:2]
sin_angle=abs(numpy.mean([-rotmat[0,1],rotmat[1,0]]))
rotmat=M2[:2,:2]
outsize = numpy.min([[realsize/px_scale*numpy.cos(numpy.arcsin(sin_angle))],[t_size]],axis=0).squeeze()
# make sure map size matches the original minimum camera pixel block limits (powers of 2)
ii=1
while (numpy.mod(t_size,2**ii)==0).all() and ii<4:
ii=ii+1
shift = numpy.mod(outsize,2**ii)/2
outsize = 2**ii*numpy.floor(outsize/(2**ii))
im3 = imcrop(im2,[t_size[1]/2+shift[1],t_size[0]/2+shift[0]],numpy.flip(outsize))
im3[im3==0]=numpy.mean(im3[:])
f_size = im3.shape
if int8:
im3=im3.astype(numpy.int8)
else:
im3=im3.astype(numpy.int16)
p4 = p1 * mat.T
p4[:,0] = f_size[1]/2 + p4[:,0]
p4[:,1] = f_size[0]/2 + p4[:,1]
o_size=numpy.max(numpy.abs([[0,realsize[1]] *mat,[realsize[0],0] *mat]),axis=0).astype(numpy.int).squeeze()
M = numpy.concatenate((mat_i,numpy.array([[0],[0]])),axis=1)
M = numpy.concatenate((M,[[0,0,1]]),axis=0)
mat1 = numpy.eye(3)
mat1[0,2] = -(o_size[1]-1)/2
mat1[1,2] = -(o_size[0]-1)/2
mat2 = numpy.eye(3)
mat2[0,2] = (realsize[1]-1)/2
mat2[1,2] = (realsize[0]-1)/2
M1 = numpy.dot(mat2,M)
M2 = numpy.dot(M1,mat1)
if int8:
im1 = numpy.round(255.0 * (im1 - im1.min()) / (im1.max() - im1.min() - 1.0)).astype(numpy.uint8)
# interpolate image
im2 = tf.warp(im1,M2,output_shape=o_size,preserve_range=True,mode='edge')#cval=numpy.median(im1[:]))
#check if output image size needs to be modified
limitsize = numpy.min([o_size,t_size],axis=0).squeeze()
if (limitsize==t_size).all():
# im3=im2.copy()
outsize = t_size
shift = [0,0]
else:
# determine limitation of image by the borders of rotated crop
rotmat=mat_i/numpy.sqrt(numpy.linalg.det(mat_i))
ca = numpy.mean(abs(numpy.diag(rotmat)))
sq2 = numpy.sqrt(2)/4
factor = 2*sq2 + sq2 *ca
outsize = numpy.min([o_size*factor,limitsize,o_size],axis=0).squeeze()
# make sure map size matches the original minimum camera pixel block limits (powers of 2)
ii=1
while (numpy.mod(t_size,2**ii)==0).all() and ii<4:
ii=ii+1
shift = numpy.mod(outsize,2**ii)/2
outsize = 2**ii*numpy.floor(outsize/(2**ii))
im3 = imcrop(im2,[o_size[1]/2+shift[1],o_size[0]/2+shift[0]],outsize)
# im3[im3==0]=numpy.median(im3[:])
f_size = im3.shape
if int8:
im3=im3.astype(numpy.int8)
else:
im3=im3.astype(numpy.int16)
p4 = p1 * mat.T
p4[:,0] = f_size[1]/2 + p4[:,0]
p4[:,1] = f_size[0]/2 + p4[:,1]
return im3, p4
return im3, p4
......@@ -1755,3 +1799,241 @@ def pointitem(label,regis=1):
point['Type'] = ['0']
return point
# ------------------------------------------------------------
def virt_map_at_point(item,idx,maps,allitems,targetitem,targetheader,outnav,numtiles=1,outformat='mrc'):
# creates a virtual map at a point item
"""
Parameters
----------
item : navitem dict
The target point item.
idx : int
acquisition index.
allitems : list of nav items
full navigator.
targetitem : navitem dict
target/reference map.
targetheader : mapheader dict
header of the target/ref map.
outnav : list
already generated virtual maps (nav items) to check for ID duplicates
numtiles : int, optional
Number of montage tiles for virtual map (squared) The default is 1.
outformat : 'tif' or 'mrc', optional
Output format of the virtual map. The default is 'mrc'.
Returns
-------
newnavitem : navitem dict
the navitem of the new virtual map.
item : navitem dict
the original navitem (with some modifications).
"""
if numtiles>1:
# montages can only be created in TIF/idoc
outformat = 'tif'
mapitem = realign_map(item,allitems)
itemid = mapitem['# Item']
delim = 100000
checkitems=allitems.copy()
checkitems.extend(outnav)
newmapid = newID(checkitems,divmod(int(targetitem['MapID'][0]),delim)[0]*delim+idx)
if not itemid in maps.keys():
maps[itemid] = mergemap(mapitem)
groupid = [str(newID(allitems,999000000+int(mapitem['MapID'][0][-6:])))]
print('add original map to navigator')
# NoRealign
mapitem['Color'] = '5'
maps['mapnav'].append(mapitem)
else:
groupid = outnav[-1]['GroupID']
map_mat = maps[itemid]['matrix']
t_mat = map_matrix(targetitem)
maptf = (numpy.linalg.inv(map_mat) * t_mat).T
xval = float(item['StageXYZ'][0]) #(float(item['PtsX'][0]))
yval = float(item['StageXYZ'][1]) #(float(item['PtsY'][0]))
pt = numpy.array([xval,yval])
# calculate the pixel coordinates
pt_px1 = get_pixel(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 = map_extract(im,pt_px1,pt_px1,px_scale,numtiles*imsz1,maptf)
im2size = im2.shape
if min(im2.shape)<200:
print('Warning! Item ' + item['# Item'] + ' is not within the map frame. Ignoring it')
newnavitem=None
else:
# pad item numbers to 4 digits
if item['# Item'].isdigit(): item['# Item'] = item['# Item'].zfill(4)
newnavitem = dict(targetitem)
if 'MapLDConSet' in targetitem.keys():
newnavitem['MapLDConSet'] = targetitem['MapLDConSet']
newmapid = newmapid + int(targetitem['MapLDConSet'][0])
if targetitem['MapLDConSet'] == ['0']:
prefix = 'V_'
newnavitem['Acquire'] = ['0']
elif targetitem['MapLDConSet'] == ['4']:
newnavitem['Acquire'] = ['0']
prefix = 'P_'
else:
prefix='m_'
newnavitem['Acquire'] = ['1']
else:
prefix=''
if numtiles <2 :
if outformat=='tif':
newnavitem['ImageType'] = ['2']
imfile = prefix + item['# Item'] + ".tif" # '.mrc'
if os.path.exists(imfile): os.remove(imfile)
io.imsave(imfile,im2.astype('uint16'), check_contrast=False)
elif outformat=='mrc':
im4 = numpy.rot90(im2,3)
imfile = prefix + item['# Item'] + '.mrc'
if os.path.exists(imfile): os.remove(imfile)
with mrc.new(imfile) as mrcf:
mrcf.set_data(im4.T)
mrcf.close()
else:
# generate virtual montage
imfile = prefix + item['# Item'] +".idoc"
newnavitem['ImageType'] = ['3']
newnavitem['MapFramesXY'] = [str(numtiles),str(numtiles)]
newnavitem['MontBinning'] = '1'
newnavitem[' MapMontage'] = '1'
MinMaxMean=' '.join(map(str,[im2.min(),im2.max(),int(round(im2.mean()))]))
idocout = ['ImageSeries = 1']
psstr = 'PixelSpacing = '+ ' %0.2f' % (targetheader['pixelsize']*10000)
idocout.append(psstr)
tilesz = (numpy.array(im2size)/numtiles).astype(int)
# tilesz=numpy.flip(tilesz)
idocout.append('ImageSize = '+str(tilesz[1])+' '+str(tilesz[0]))
idocout.append('Montage = 1')
idocout.append('DataMode = 6')
idocout.append('')
for tile in range(numtiles**2):
tilepos = list(numpy.divmod(tile,numtiles))
tilepos0= tilepos.copy()
tilepos[0]=numtiles-tilepos[0]-1
imrange = [[(tilepos[ix]*tilesz[ix]),((tilepos[ix]+1)*tilesz[ix])] for ix in range(2)]
im3 = im2[imrange[0][0]:imrange[0][1],imrange[1][0]:imrange[1][1]]
tilefile = prefix + item['# Item'] + '_' + str(tile)+".tif"
if os.path.exists(imfile): os.remove(imfile)
io.imsave(tilefile,im3.astype('uint16'), check_contrast=False)
idocout.append('[Image = '+tilefile+']')
idocout.append('PieceCoordinates = '+str(tilepos0[1]*tilesz[1])+' '+str(tilepos0[0]*tilesz[0])+' 0')
idocout.append('MinMaxMean = '+MinMaxMean)
idocout.append('')
with open(imfile,'w') as imf:
imf.writelines("%s\n" % line for line in idocout)
cx = im2.shape[1]
cy = im2.shape[0]
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
item['Acquire'] = '0'
# NoRealign
# item['Color'] = '5'
newnavitem['MapFile'] = [imfile]
newnavitem['StageXYZ'] = item['StageXYZ']
newnavitem['RawStageXY'] = 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'] = item['MapID']
newnavitem['GroupID'] = groupid
newnavitem['MapWidthHeight'] = [str(im2size[1]),str(im2size[0])]
newnavitem['MapMinMaxScale'] = [str(numpy.min(im2)),str(numpy.max(im2))]
newnavitem['NumPts'] = ['5']
newnavitem['# Item'] = prefix + item['# Item']
return newnavitem, maps, item
......@@ -6,7 +6,11 @@ with open("readme.md", "r") as fh:
long_description = fh.read()
setup(name='py-EM',
<<<<<<< HEAD
version='20201216',
=======
version='20210406',
>>>>>>> dev
py_modules=['pyEM'],
description='Tools to interact with Serial EM enabling automated transmission electron microscopy.',
long_description=long_description,
......
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 6 10:49:32 2021
@author: schorb
"""
import requests
import zipfile
import os
from subprocess import Popen, PIPE
# import hashlib
cwd = os.getcwd()
if not os.path.exists('pyem_tests'):
# download and prepare test data
print('downloading test data ...')
url = 'https://oc.embl.de/index.php/s/NDAkNzEr7zjSZ3Q/download'
myfile = requests.get(url)
open('test.zip','wb').write(myfile.content)
with zipfile.ZipFile('test.zip', 'r') as unzip:
unzip.extractall()
os.remove('test.zip')
if not os.path.exists('applications'):
# download test scripts
os.mkdir('applications')
url = 'https://git.embl.de/schorb/pyem/-/raw/master/applications/virt_anchormaps.py?inline=false'
anchfile = requests.get(url)
open('applications/virt_anchormaps.py','wb').write(anchfile.content)
url = 'https://git.embl.de/schorb/pyem/-/raw/master/applications/maps_acquire_cmd.py?inline=false'
reffile = requests.get(url)
open('applications/maps_acquire_cmd.py','wb').write(reffile.content)
#start tests
print('Testing map extraction for Tecnai...')
os.chdir(os.path.join(cwd,'pyem_tests'))
p1 = Popen('python "'+os.path.join(cwd,'applications','maps_acquire_cmd.py')+'" "'+ os.path.join(cwd,'pyem_tests','tecnai.nav')+'"', shell=True,stderr=PIPE,stdout=PIPE)
# test_hashlist = []
# print('comparing file hashes')
# for file in ['tecnai_automaps.nav','virt_map_p1.mrc','virt_map_p2.mrc','virt_map_p3.mrc','virt_map_p4.mrc']:
# test_hashlist.append(hashlib.md5(open(file,'rb').read()).hexdigest())
print('Testing map extraction for Krios...')
p2 = Popen('python "'+os.path.join(cwd,'applications','virt_anchormaps.py')+'" "'+ os.path.join(cwd,'pyem_tests','cryo_test_edge.nav')+'"', shell=True,stderr=PIPE,stdout=PIPE)
p1.wait()
p2.wait()
print('Now, please check the *.automaps.nav and linked virtual maps in SerialEM and compare them with the files in the `result` directory. ')
\ No newline at end of file
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