Merge branch 'icy' into 'master'

Icy

See merge request schorb/pyem!4

applications/clem_icy.py

@@ -11,8 +11,8 @@ import xml.etree.ElementTree as ET
 
 import sys
 import os
-
-
+import glob
+import numpy as np
 
 
 
@@ -122,6 +122,108 @@ for regpt in em_regpts:
 pts2xml(mm_fm['mapfile'], fm_px)
 pts2xml(mm_em['mergefile'], em_px)
 
+
+##                ICY   runs here      
+
+
+print('----------------------------------------------------------------')
+print('        Starting Icy..................................')
+print('----------------------------------------------------------------')
+print(' Please run ec-CLEM and when done click SHOW ROIs ON ORIGINAL SOURCE IMAGE\n')
+
+
+workdir = os.getcwd()
+
+os.chdir('C:\Software\icy')
+icycmd = 'java -jar icy.jar -x plugins.tprovoost.scripteditor.main.ScriptEditorPlugin C:\Software\opener.js'
+
+os.system(icycmd +' '+ mm_fm['mapfile'] +' '+  mm_em['mergefile'] + os.path.splitext(mm_em['mapfile'])[1])
+
+os.chdir(workdir)
+
+
+#%%
+
+# import icy XMLs
+
+x_emf = mm_em['mergefile']+ os.path.splitext( mm_em['mapfile'])[1] + '_ROIsavedwhenshowonoriginaldata' + '.xml'
+x_trafo = mm_fm['mapfile'] + '_transfo' + '.xml'
+
+if not all([os.path.exists(x_emf),os.path.exists(x_trafo)]):  raise ValueError('Could not find the results of Icy registration. Please re-run.')
+    
+root_em = ET.parse(x_emf)
+root_trafo = ET.parse(x_trafo)
+
+
+# extract transformation matrices
+
+mat = []
+if root_trafo.find('MatrixTransformation') == None:
+    root_trafo = ET.parse(glob.glob(mm_fm['mapfile'] + '_transfo' + '*back-up.xml')[-1])
+
+
+for child in root_trafo.iter():
+    if child.tag == 'MatrixTransformation':
+        mat.append(child.attrib)
+
+M = np.eye(4)
+M_list = list()
+        
+for matrix in mat:
+    thismat = np.eye(4)
+    thismat[0,0] = matrix['m00']
+    thismat[0,1] = matrix['m01']
+    thismat[0,2] = matrix['m02']
+    thismat[0,3] = matrix['m03']
+    thismat[1,0] = matrix['m10']
+    thismat[1,1] = matrix['m11']
+    thismat[1,2] = matrix['m12']
+    thismat[1,3] = matrix['m13']
+    thismat[2,0] = matrix['m20']
+    thismat[2,1] = matrix['m21']
+    thismat[2,2] = matrix['m22']
+    thismat[2,3] = matrix['m23']
+    thismat[3,0] = matrix['m30']
+    thismat[3,1] = matrix['m31']
+    thismat[3,2] = matrix['m32']
+    thismat[3,3] = matrix['m33']
+    
+    M_list.append(thismat)
+    
+    M = np.dot(M.T,thismat)
+    
+
+
+# extract registration point positions
+
+
+p_idx = list()
+pts = list()
+
+for child in root_em.iter():
+    pt = np.zeros(2)
+    if child.tag == 'position':
+        
+        for coords in child.getiterator():            
+            if coords.tag == 'pos_x':
+                pt[0] = float(coords.text)
+            elif coords.tag == 'pos_y':
+                pt[1] = float(coords.text)
+            #elif coords.tag == 'pos_z':
+            #    pt[2] = float(coords.text)
+            
+        pts.append(pt)
+            
+    elif child.tag == 'name':
+        p_idx.append(int(child.text[child.text.rfind(' '):]))
+                
+    
+### TODO ##
+
+
+                
+                
+
 # create new file by copying the header of the input file
 newnavf = navfile[:-4] + '_icy.nav'
 

pyEM.py

@@ -585,11 +585,12 @@ def mergemap(mapitem,crop=False,black=False,blendmont=True):
         imod_vercheck = (imod_ver[0]>=4 and imod_ver[1]>=10 and imod_ver[2]>=42)
 
         if blendmont:
-            mergefile = mbase + '_merged'+ '_s' + str(mapsection)
-            if not os.path.exists(mergefile+'.mrc'):                
+            mergebase = mbase + '_merged'+ '_s' + str(mapsection)
+            mergefile = mergebase+'.mrc'
+            if not os.path.exists(mergefile):                
                 if imod_vercheck:
-                    call_blendmont(mapfile,mergefile,mapsection)
-                    merge_mrc =  mrc.mmap(mergefile + '.mrc', permissive = 'True')
+                    call_blendmont(mapfile,mergebase,mapsection)
+                    merge_mrc =  mrc.mmap(mergefile, permissive = 'True')
                     im = merge_mrc.data
                     mergeheader = map_header(merge_mrc)
                 else:
@@ -598,7 +599,7 @@ def mergemap(mapitem,crop=False,black=False,blendmont=True):
                     mergeheader['ysize'] = int(tilepx[-1][1]) + mapheader['ysize']
                     mergefile = mapfile
             else:
-                merge_mrc =  mrc.mmap(mergefile + '.mrc', permissive = 'True')
+                merge_mrc =  mrc.mmap(mergefile, permissive = 'True')
                 im = merge_mrc.data
                 mergeheader = map_header(merge_mrc)
                            
@@ -716,11 +717,9 @@ def mergemap(mapitem,crop=False,black=False,blendmont=True):
 
 
     # check if map is a montage or not
-    mergefile = mapfile[:mapfile.rfind('.mrc')]
-    if mergefile == []: mergefile = mapfile
-
-    mergefile = mergefile + '_merged'+ '_s' + str(mapsection)
 
+    mergebase = os.path.splitext(mapfile)[0] + '_merged'+ '_s' + str(mapsection)
+    mergefile = mergebase + '.mrc'
 
     if mapheader['stacksize'] < 2:
         print('Single image found. No merging needed.')
@@ -742,15 +741,15 @@ def mergemap(mapitem,crop=False,black=False,blendmont=True):
 
     else:
         if blendmont:
-            if not os.path.exists(mergefile+'.mrc'):
-                call_blendmont(mapfile,mergefile,mapsection,black)
+            if not os.path.exists(mergefile):
+                call_blendmont(mapfile,mergebase,mapsection,black)
                 
-            merge_mrc =  mrc.mmap(mergefile + '.mrc', permissive = 'True')
+            merge_mrc =  mrc.mmap(mergefile, permissive = 'True')
             im = merge_mrc.data
             mergeheader = map_header(merge_mrc)            
     
                 # extract pixel coordinate of each tile
-            tilepx = list(loadtext(mergefile + '.al'))
+            tilepx = list(loadtext(mergebase + '.al'))
             for j, item in enumerate(tilepx): tilepx[j] = list(re.split(' +',item))
             
             # use original tile coordinates(pixels) from SerialEM to determine tile position in montage
@@ -811,34 +810,36 @@ def mergemap(mapitem,crop=False,black=False,blendmont=True):
   overlapx = mapheader['xsize'] - xstep
   overlapy = mapheader['ysize'] - ystep
   
+  print(mergefile)
+  print(mergebase)
+  
       
   # cropping of merged file to only include areas of interest. The user needs to create an IMOD model file and close 3dmod to proceed.
   if crop & blendmont:
-      if not os.path.exists(mergefile+'_crop.mrc'):
+      if not os.path.exists(mergebase+'_crop.mrc'):
           loopcount = 0
-          print('waiting for crop model to be created ... Please store it under this file name: \"' + mergefile + '.mod\".')
-          callcmd = '3dmod \"' +  mergefile + '.mrc\" \"' + mergefile + '.mod\"'
+          print('waiting for crop model to be created ... Please store it under this file name: \"' + mergebase + '.mod\".')
+          callcmd = '3dmod \"' +  mergefile + '\" \"' + mergebase + '.mod\"'
           os.system(callcmd)
-          while not os.path.exists(mergefile+'.mod'):
+          while not os.path.exists(mergebase+'.mod'):
               if loopcount > 20: # wait for 6.5 minutes for the model file to be created.
                   print('Timeout - will continue without cropping!')
                   break
-              print('waiting for crop model to be created in IMOD... Please store it under this file name: \"' + mergefile + '.mod\".')
+              print('waiting for crop model to be created in IMOD... Please store it under this file name: \"' + mergebase + '.mod\".')
               time.sleep(20)
               loopcount = loopcount + 1
               
           if loopcount < 21:
               print('Model file found. Will now generate the cropped map image.')
-              callcmd = 'imodmop \"' +  mergefile + '.mod\" \"'+ mergefile + '.mrc\" \"' + mergefile + '_crop.mrc\"'
+              callcmd = 'imodmop \"' +  mergebase + '.mod\" \"'+ mergefile + '\" \"' + mergebase + '_crop.mrc\"'
               os.system(callcmd)
               
               merge_mrc.close()
   
-      merge_mrc = mrc.mmap(mergefile + '_crop.mrc', permissive = 'True')
+      merge_mrc = mrc.mmap(mergebase + '_crop.mrc', permissive = 'True')
       im_cropped = merge_mrc.data
       im_cropped = numpy.rot90(numpy.transpose(im_cropped))
-      m['im_cropped'] = im_cropped        
-      mergefile = mergefile+'.mrc'  
+      m['im_cropped'] = im_cropped       
       
   
   mergeheader['pixelsize'] = pixelsize
@@ -866,7 +867,28 @@ def mergemap(mapitem,crop=False,black=False,blendmont=True):
 
 # -------------------------------
 #%%
-def call_blendmont(mapfile,mergefile,mapsection,black=False):
+def call_blendmont(mapfile,mergebase,mapsection,black=False):
+    """
+    
+
+    Parameters
+    ----------
+    mapfile : string
+        Location of the source map file.
+    mergebase : string
+        Target file name for merged map (without extension!).
+    mapsection : int
+        which slice of the stack to process.
+    black : TYPE, bool
+        whether the background of the merged map hav value 0. The default is False.
+
+    Returns
+    -------
+    None.
+
+    """
+    
+    
     # check IMOD version
     if imod_ver[0]<4 | imod_ver[1]<10 | imod_ver[2]<29 :
         print('ERROR: IMOD version needs to be > 4.10.29! Please update. Exiting' + '\n')
@@ -881,7 +903,7 @@ def call_blendmont(mapfile,mergefile,mapsection,black=False):
     print('Merging the map montage into a single image....' + '\n')
     print('----------------------------------------------------\n')
     
-    callcmd = 'blendmont -imi ' + '\"' + mapfile + '\"' + ' -imo \"' + mergefile + '.mrc\" -pli \"' + mapfile + '.pcs\" -roo \"' + mergefile  + '.mrc\" -se ' + str(mapsection) + ' -al \"'+ mergefile + '.al\" -sloppy -nofft '    #os.system(callcmd)
+    callcmd = 'blendmont -imi ' + '\"' + mapfile + '\"' + ' -imo \"' + mergebase + '.mrc\" -pli \"' + mapfile + '.pcs\" -roo \"' + mergebase  + '.mrc\" -se ' + str(mapsection) + ' -al \"'+ mergebase + '.al\" -sloppy -nofft '    #os.system(callcmd)
     if black:
         callcmd = callcmd + ' -fill 0'
     

setup.py

@@ -6,7 +6,7 @@ with open("readme.md", "r") as fh:
     long_description = fh.read()
 
 setup(name='py-EM',
-      version='20200609',
+      version='20200618',
       py_modules=['pyEM'],
       description='Tools to interact with Serial EM enabling automated transmission electron microscopy.',
       long_description=long_description,