pyEM.py

# -*- coding: utf-8 -*-

# pyEM.py

# Copyright (c) 2020, Martin Schorb
# Copyright (c) 2020, European Molecular Biology Laboratory
# Produced at the EMBL
# All rights reserved.

"""    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <https://www.gnu.org/licenses/>.
"""

# Python module for interpreting and modifying Navigator files created using SerialEM. For detailed information,
# see  https://doi.org/10.1038/s41592-019-0396-9
# Information and documentation of the individual functions can be found in Readme.md

# dependencies
#
import os
import os.path
import sys
import numpy
import skimage.transform as tf
import copy
from skimage import io
import re
import mrcfile as mrc
import time
from operator import itemgetter
import fnmatch
from subprocess import Popen, PIPE
import xml.etree.ElementTree as ET
import xml.dom.minidom

mrcext = ('.st', '.mrc', '.map', '.rec', '.ali', '.preali')

# get python version
py_ver = sys.version_info

my_env = os.environ.copy()

# get IMOD version
p1 = Popen("imodinfo", shell=True, stdout=PIPE)
o = list()
for imodline in p1.stdout:
    o.append(imodline)

if o == []:
    print('Warning! - No IMOD installation found')
    imod_ver = [0, 0, 0]
else:
    o1 = str(o[0]).split(' ')
    imod_ver = list(map(int, o1[o1.index('Version') + 1].split('.')))


# define functions


# %%

def loadtext(fname):
    """
    Reads a text file and returns the lines.

    Parameters
    ----------
    fname : str
        input path of a text file, such as nav (adoc or XML) or adoc

    Returns
    -------
    list of str

    """

    # check if file exists
    if not os.path.exists(fname):
        raise FileNotFoundError('ERROR: ' + fname + ' does not exist! Exiting' + '\n')

    f = open(fname, "r")

    lines = list()

    for line in f.readlines():
        lines.append(line.strip())

    f.close()
    return lines


# -------------------------------
# %%

def nav_item(inlines, label):
    """
    Extracts the content block of a single navItem of given label.
    Reads and parses navigator adoc files version >2, also in XML format !!
    Returns the first item found as a dictionary and the remaining list with that item removed.
    This is useful when multiple items have the exact same label and the function is called \\
    from within a loop to retrieve them all.

    Parameters
    ----------
    inlines : list of str
        input navigator file lines

    label : str
        the navigator label to extract

    Returns
    -------
    dict

    """
    #

    lines = inlines[:]

    if lines[0] == '<?xml version="1.0" encoding="utf-8"?>':
        # load XML
        root = ET.fromstringlist(lines)
        el = root.findall('*[@name="%s"]' % label)
        if el == []:
            print('ERROR: Navigator Item ' + label + ' not found!')
            result = []
        else:
            root.remove(el[0])
            newroot = ET.Element('navigator')
            newroot.append(el[0])

            result = xmltonav(ET.tostringlist(newroot, encoding='unicode'))[0]
            lines = ET.tostring(root)
            lines = xml.dom.minidom.parseString(lines).toprettyxml().replace('\t', '').split('\n')
    else:
        if lines[-1] != '':
            lines = lines + ['']

        # search for a navigator item of given label
        searchstr = '[Item = ' + label + ']'
        if searchstr not in lines:
            print('ERROR: Navigator Item ' + label + ' not found!')
            result = []
        else:
            itemstartline = lines.index(searchstr) + 1
            itemendline = lines[itemstartline:].index('')

            item = lines[itemstartline:itemstartline + itemendline]
            result = parse_adoc(item)
            # create a dictionary entry that contains the label. Comment # is given in case it is exported directly.
            result['# Item'] = lines[itemstartline - 1][lines[itemstartline - 1].find(' = ') + 3:-1]

            lines[itemstartline - 1:itemstartline + itemendline + 1] = []

    return result, lines


# -------------------------------
# %%

def adoc_items(lines1, search, header=False):
    """
    Extracts the content block of an item of given label in an adoc file and \\
    returns it as a list of dictionaries

    Parameters
    ----------
    lines1 : list of str
        input adoc file lines
    search : str
        the search string
    header : bool
        if True, returns the header of the file instead.

    Returns
    -------
    list of dict

    """

    result = []

    if lines1[-1] != '':
        lines = lines1 + ['']
    else:
        lines = lines1

    if header:
        item = lines[:lines.index('')]
        result.append(parse_adoc(item))
    else:
        pattern = '*' + search + '*'
        matching = fnmatch.filter(lines, pattern)

        if len(matching) == 0:
            print('ERROR: String ' + search + ' not found!')

        # search for mdoc key item with the given label
        for searchstr in matching:
            if searchstr not in lines:
                print('ERROR: String ' + search + ' not found!')
            elif '=' in searchstr[:searchstr.index(search)]:
                print('ERROR: String ' + search + ' not found in labels!')
            else:
                itemstartline = lines.index(searchstr) + 1
                itemendline = lines[itemstartline:].index('')

                item = lines[itemstartline:itemstartline + itemendline]

                itemdict = parse_adoc(item)
                itemdict['# ' + searchstr] = lines[itemstartline - 1][
                                             lines[itemstartline - 1].find(' = ') + 2:-1].lstrip(' ')
                result.append(itemdict)

    return result


# -------------------------------
# %%

def mdoc_item(lines1, label, header=False):
    """
    Extracts the content block of an item of given label in a mdoc file and \\
    returns it as a dictionary

    Parameters
    ----------
    lines1 : list of str
        input mdoc file lines
    label : str
        the mdoc label to search for
    header : bool
        if True, returns the header of the file instead.

    Returns
    -------
    list of dict

    """

    if lines1[-1] != '':
        lines = lines1 + ['']
    else:
        lines = lines1

    if header:
        item = lines[:lines.index('')]
    else:
        # search for mdoc key item with the given label
        searchstr = '[' + label + ']'
        if searchstr not in lines:
            print('ERROR: Item ' + label + ' not found!')
            item = []
        else:
            itemstartline = lines.index(searchstr) + 1
            itemendline = lines[itemstartline:].index('')

            item = lines[itemstartline:itemstartline + itemendline]
    result = parse_adoc(item)

    return result


# -------------------------------
# %%

def parse_adoc(lines):
    """
    Converts an adoc-format string list into a dictionary.

    Parameters
    ----------
    lines1 : list of str
        input adoc file lines

    Returns
    -------
    dict

    """
    #

    output = {}
    for line in lines:
        entry = line.split()
        if entry:
            output.update({entry[0]: entry[2:]})

    return output


# -------------------------------
# %%

def map_file(mapitem):
    """
    Extracts map file name from navigator and checks for existance.
    Checks all sub-directories if file is not present in the root (navigator location).

    Parameters
    ----------
    mapitem : dict
        a Map Item dict from the Navigator

    Returns
    -------
    str

    """

    # get string from navigator item
    mapfile = ' '.join(mapitem['MapFile'])
    cdir = os.getcwd()

    if os.path.exists(mapfile):
        return mapfile
        # print('current map: ' + mapfile)
    else:
        #    print('Warning: ' + mapfile + ' does not exist!' + '\n')
        mapfile1 = mapfile[mapfile.rfind('\\') + 1:]
        dir1 = mapfile[:mapfile.rfind('\\')]
        dir2 = dir1[dir1.rfind('\\') + 1:]

        #  print('will try ' + mapfile1 + ' in current directory or subdirectories.' + '\n')

        # check subdirectories recursively
        mapfile2 = None

        for path, directories, files in os.walk(cdir):
            if mapfile1 in files:
                mapfile = os.path.join(path, mapfile1)

                if path == os.path.join(cdir, dir2):
                    return mapfile
                else:
                    mapfile2 = mapfile

        if mapfile2 is None:
            raise FileNotFoundError('ERROR: ' + mapfile1 + ' does not exist! Exiting.')
        else:
            return mapfile2


# -------------------------------
# %%
# def findfile(searchstr, searchdir):
#     # will find files that match a search string in subfolders of the provided search directory
#     output = list()
#     for rootdir, dirs, files in os.walk(searchdir):
#         for file in files:
#             if fnmatch.fnmatch(file, searchstr):
#                 output.append(os.path.join(rootdir, file))
#     return output


# -------------------------------
# %%
def map_header(m):
    """
    Extrracts the header of a map file.

    Parameters
    ----------
    m : str or dict or mrc.mrcfile.MrcFile
        A map file, map dictionary or MRC object

    Returns
    -------
    dict

    """
    header = {}

    if (type(m) is mrc.mrcfile.MrcFile) | (type(m) is mrc.mrcmemmap.MrcMemmap):
        # extracts MRC header information for a given mrc.object (legacy from reading mrc headers)

        header['xsize'] = int(m.header.nx)
        header['ysize'] = int(m.header.ny)

        header['stacksize'] = int(m.header.nz)

        # determine the scale

        header['pixelsize'] = m.voxel_size.x / 10000  # in um
    elif type(m) is dict:
        if "MapFile" in m.keys():
            merge = mergemap(m)# THERE IS A BUG!!!!!, blendmont=False)
            header = merge['mapheader']
    elif type(m) is str:
        if os.path.exists(m):
            if os.path.splitext(m)[1] in mrcext:
                header = map_header(mrc.mmap(m, permissive='True'))

    return header


# -------------------------------
# %%

def itemtonav(item, name):
    """
    Converts a dictionary autodoc item variable into text list suitable for export into navigator format.

    Parameters
    ----------
    item : dict
        Navigator Item
    name : str
        Item label name

    Returns
    -------
    list of str

    """

    dlist = list()
    dlist.append('[Item = ' + name + ']')

    # Python 2
    if py_ver[0] < 3:
        for key, value in item.iteritems():
            if not key == '# Item':
                dlist.append(key + ' = ' + " ".join(value))
    else:
        # Python 3+
        for key, value in item.items():
            if not key == '# Item':
                dlist.append(key + ' = ' + " ".join(value))
    dlist.append('')
    return dlist


# -------------------------------

# %%

def write_navfile(filename, outitems, xml=False):
    """
    Creates a new navigator file from a list of navItems (default is mdoc format)


    Parameters
    ----------
    filename : str
        File name to ouput (needs to be writable.)
    outitems : list of dict
        A list of items to be written.
    xml : bool
        If true, write Navigator in XML format.

    """

    allitems = copy.deepcopy(outitems)

    if xml:
        root = ET.Element('navigator')
        pd = ET.SubElement(root, 'PreData')
        adv = ET.SubElement(pd, 'AdocVersion')
        adv.text = '2.00'
        lsa = ET.SubElement(pd, 'LastSavedAs')
        lsa.text = filename

        for item in allitems:
            print(item)
            print(item['# Item'])

            ci = ET.SubElement(root, 'Item')
            ci.set('name', item['# Item'])

            for key, val in item.items():
                if not key == '# Item':
                    cp = ET.SubElement(ci, key)
                    cp.text = ' '.join(val)
        indent_xml(root)
        tree = ET.ElementTree(root)
        tree.write(filename, encoding="utf-8", xml_declaration=True)

    else:
        # MDOC format
        head0 = 'AdocVersion = 2.00'
        head1 = 'LastSavedAs = ' + filename

        nnf = open(filename, 'w')
        nnf.write("%s\n" % head0)
        nnf.write("%s\n" % head1)
        nnf.write("\n")

        # fill the new file
        for nitem in allitems:
            out = itemtonav(nitem, nitem['# Item'])
            for item in out:
                nnf.write("%s\n" % item)

        nnf.close()


# -------------------------------

# pretty print xml, from:
# http://effbot.org/zone/element-lib.htm#prettyprint

def indent_xml(elem, level=0):
    i = "\n" + level * "  "
    if len(elem):
        if not elem.text or not elem.text.strip():
            elem.text = i + "  "
        if not elem.tail or not elem.tail.strip():
            elem.tail = i
        for elem in elem:
            indent_xml(elem, level + 1)
        if not elem.tail or not elem.tail.strip():
            elem.tail = i
    else:
        if level and (not elem.tail or not elem.tail.strip()):
            elem.tail = i


# ------------------------------


# %%

def newID(allitems, startid):
    """
    Checks if provided item ID already exists in a navigator and gives the next unique ID.

    Parameters
    ----------
    allitems : list of dict
        list of all items in a Navigator (<-fullnav)
    startid : int
        starting ID

    Returns
    -------
    int

    """

    newid = startid

    for item in allitems:
        if 'MapID' in item:
            if str(startid) == item['MapID'][0]:
                newid = newID(allitems, startid + 1)

    return newid


# -------------------------------
# %%

def newreg(navitems):
    """
    Returns the next available registration for the input set of navigator items

    Parameters
    ----------
    navitems : list of dict
        list of Navigator items

    Returns
    -------
    int

    """
    #
    reg = list()
    for item in navitems:
        reg.append(int(item['Regis'][0]))

    return max(reg) + 1


# -------------------------------

# %%

def fullnav(inlines, header=False):
    """
    Parses a full nav file and returns a list of dictionaries.

    Parameters
    ----------
    inlines :  list of str
        input navigator file lines

    header : bool
        if True, returns the header of the file instead.

    Returns
    -------
    list of dict

    """
    navlines = inlines[:]

    if '<?xml version=' in navlines[0]:
        # load XML
        c = xmltonav(navlines)
    else:
        # mdoc format
        c = []
        for index, item in enumerate(navlines):
            if item.find('[Item') > -1:
                if header:
                    return inlines[:index - 1]

                (b, navlines) = nav_item(navlines, item[item.find(' = ') + 3:-1])
                b['# Item'] = item[item.find(' = ') + 3:-1]
                c.append(b)

    return c


# -------------------------------
# %%

def xmltonav(navlines):
    """
    Parses a XML block and returns the contents as a list of dictionaries.

    Parameters
    ----------
    inlines :  list of str
        input XML file lines

    Returns
    -------
    list of dict

    """
    if not navlines[0].lstrip(' ').startswith('<'):
        raise ValueError('Input string list not in XML format!')

    root = ET.fromstringlist(navlines)

    c = []

    for child in root:
        if child.tag == 'Item':
            item = dict({'# Item': child.attrib['name']})
            for prop in child:
                item[prop.tag] = prop.text.split(' ')
            c.append(item)

    return c


# -------------------------------
# %%

def navlabel_match(navitems, searchstr):
    """
    Identifies navigator items whose labels contain the given string

    Parameters
    ----------
    navitems : list of dict
        List of Navigator items
    searchstr : str
        String to be searched

    Returns
    -------
    list of dict

    """

    r = re.compile(r'.*' + searchstr + '.*')

    return list(filter(lambda item: r.match(item['# Item']), navitems))


# -------------------------------
# %%

def duplicate_items(navitems, labels=None, prefix='', reg=True, maps=False):
    """
    Duplicates a set of Navigator items. Adds the duplicates to the original Navigator list.

    Parameters
    ----------
    navitems : list of dict
        List of Navigator items
    labels : bool
        list of labels of the items to duplicate
    prefix : str
        prefix to be added to the label of the created duplicates
    reg : bool
        determines whether the registration of duplicate items should be changed (default:yes)
    maps :bool
        if True, all maps that contain the label of the selected items or that were used to draw \\
         these are duplicated as well.

    Returns
    -------
    list of dict

    """

    if labels is None:
        labels = []
    outitems = copy.deepcopy(navitems)

    if labels == []:
        dupitems = nav_selection(navitems)
        for item in dupitems:
            item['Acquire'] = ['0']

    else:
        dupitems = nav_selection(navitems, labels, False)

    if reg:
        new_reg = [str(newreg(dupitems))]
    else:
        new_reg = dupitems[0]['Regis']

    for item in dupitems:
        newitem = copy.deepcopy(item)

        newitem['Regis'] = new_reg
        newitem['# Item'] = prefix + item['# Item']
        newitem['MapID'] = [str(newID(outitems, int(newitem['MapID'][0])))]

        if maps:
            drawnmap = realign_map(newitem, navitems)
            newmapID = newID(outitems, int(drawnmap['MapID'][0])) + 1

            if not drawnmap == []:
                dupdrawn = copy.deepcopy(drawnmap)
                dupdrawn['# Item'] = prefix + drawnmap['# Item']
                dupdrawn['Regis'] = new_reg

                existingmaps = nav_find(outitems, '# Item', val=dupdrawn['# Item'])
                existingdupmaps = nav_find(existingmaps, 'Regis', val=dupdrawn['Regis'])
                dupdrawn['MapID'] = [str(newmapID)]

                if len(existingdupmaps) == 0:
                    othermaps = navlabel_match(navitems, dupdrawn['# Item'])
                    othermaps.pop(othermaps.index(nav_selection(othermaps, sel=dupdrawn['# Item'], acquire=False)[0]))

                    outitems.append(dupdrawn)

                newitem['DrawnID'] = dupdrawn['MapID']
        outitems.append(newitem)

    return outitems


# -------------------------------------


def map_matrix(mapitem):
    """
    Extracts the matrix relating pixel and stage coordinates from a Map Item
    Parameters
    ----------
    mapitem : dict
        Navigator map item

    Returns
    -------
    numpy.array

    """

    return numpy.array(list(map(float, mapitem['MapScaleMat']))).reshape(2, 2) * (
            int(mapitem['MapBinning'][0]) / int(mapitem['MontBinning'][0]))


# -------------------------------
# %%

def mergemap(mapitem, crop=False, black=False,
             blendmont=True, bigstitch=False, mapfile=None, mapsection=0):
    # %%
    # processes a map item and merges the mosaic using IMOD
    # generates a dictionary with metadata for this procedure
    # if crop is selected, a 3dmod session will be opened and the user needs to draw a model of the desired region.
    # The script continues after saving the model file and closing 3dmod.
    # black option will fill the empty spaces between tiles with 0
    # blendmont will use blendmont to merge the montage.
    # If disabled individual tile information is stored in the merge dict nevertheless.

    blend_in = blendmont
    m = dict()
    m['Sloppy'] = False

    # check if function is called with only a map file without navigator

    if mapfile is None:

        # extract map properties

        mat = map_matrix(mapitem)

        # find map file
        mapfile = map_file(mapitem)
        # print('processing mapitem '+mapitem['# Item']+' - file: '+mapfile)
        mapsection = list(map(int, mapitem['MapSection']))[0]

        m['frames'] = list(map(int, mapitem['MapFramesXY']))

    else:

        mat = numpy.eye(2)

        callcmd = 'extractpieces ' + mapfile + ' -ou ' + mapfile + '.pcs'
        os.system(callcmd)
        pxpos = loadtext(mapfile + '.pcs')
        tilepx = numpy.array([re.split(' +', line) for line in pxpos])
        frames = numpy.unique(tilepx[:, 0:2], axis=0)
        m['frames'] = [frames.shape[0], 1]

    montage_tiles = numpy.prod(m['frames'])

    tileidx_offset = 0
    mbase = os.path.splitext(mapfile)[0]
    mergebase = mbase + '_merged' + '_s' + str(mapsection)

    if os.path.splitext(mapfile)[1] not in mrcext:
        # not an mrc file

        mergeheader = {}
        # mergeheader['stacksize'] = 1

        if '.idoc' in mapfile:
            # List of tif files with additional metadata

            idoctxt = loadtext(mapfile)

            tilepos = list()
            tilepx = list()
            tilepx1 = list()

            # find the stack size from the last index of the file list (idoc)

            testlast = idoctxt.copy()
            testlast.reverse()
            lastitem = testlast[0]

            for index, item in enumerate(testlast):
                if item.strip() == '':
                    lastitem = testlast[index - 1]
                    if 'Image = ' in lastitem:
                        break

            prefix = os.path.basename(mbase)  # [:mbase.find('.idoc')]
            stacksize = int(lastitem[lastitem.find(prefix) + len(prefix):-5]) + 1
            mergeheader['stacksize'] = stacksize

            maphead0 = mdoc_item(idoctxt, [], header=True)
            im = list()

            if montage_tiles == 0:
                # stack of single images
                mergefile = mapfile[:mapfile.find('.idoc')] + '{:04d}'.format(mapsection) + '.tif'
                mergeheader['stacksize'] = 0
                tilepos = mapitem['StageXYZ'][0:2]
                tilepx = '0'
                tilepx = numpy.array([[tilepx, tilepx, mapsection], [tilepx, tilepx, mapsection]])
                tilepx1 = tilepx
                m['Sloppy'] = 'NoMont'
                im = io.imread(mergefile)
                mergeheader['xsize'] = numpy.array(im.shape)[0]
                mergeheader['ysize'] = numpy.array(im.shape)[1]
                mapheader = mergeheader.copy()
                pixelsize = float(maphead0['PixelSpacing'][0]) / 10000
            else:
                if mdoc_item(idoctxt, 'MontSection = 0') == []:  # older mdoc file format, created before SerialEM 3.7x
                    print('Warning - item' + mapitem[
                        '# Item'] + ': Series of tif images without montage information.'
                                    ' Assume pixel size is consistent for all sections.')

                    pixelsize = float(maphead0['PixelSpacing'][0]) / 10000  # in um

                else:
                    mont_item = mdoc_item(idoctxt, 'MontSection = ' + str(mapsection))
                    pixelsize = float(mont_item['PixelSpacing'][0]) / 10000  # in um

                for i in range(0, numpy.min([montage_tiles, stacksize])):
                    tilefile = mapfile[:mapfile.find('.idoc')] + '{:04d}'.format(mapsection + i) + '.tif'
                    im.append(tilefile)
                    tile = mdoc_item(idoctxt, 'Image = ' + os.path.basename(tilefile))
                    tilepos.append(tile['StagePosition'])
                    tilepx1.append(tile['PieceCoordinates'])

                    if 'AlignedPieceCoordsVS' in tile:
                        m['Sloppy'] = True
                        tilepx.append(tile['AlignedPieceCoordsVS'])
                    else:
                        tilepx.append(tile['AlignedPieceCoords'])

                mergeheader['pixelsize'] = pixelsize
                mapheader = mergeheader.copy()
                mapheader['ysize'] = int(maphead0['ImageSize'][1])
                mapheader['xsize'] = int(maphead0['ImageSize'][0])

                imsz_x = int(maphead0['ImageSize'][0])
                imsz_y = int(maphead0['ImageSize'][1])
                #
                # overlapx = imsz_x - mapheader['xsize']
                # overlapy = imsz_y - mapheader['ysize']

                # check if idoc is supported in IMOD (blendmont)
                imod_vercheck = (imod_ver[0] >= 4 or (imod_ver[0] == 4 and imod_ver[1] >= 10))

                if blendmont:
                    mergebase = mbase + '_merged' + '_s' + str(mapsection)
                    mergefile = mergebase + '.mrc'
                    if not os.path.exists(mergefile):
                        if imod_vercheck:
                            call_blendmont(mapfile, mergebase, mapsection)
                            merge_mrc = mrc.mmap(mergefile, permissive='True')
                            im = merge_mrc.data
                            mergeheader = map_header(merge_mrc)
                        else:
                            print('Please update IMOD to > 4.10.42 for merging idoc montages!')
                            mergeheader['xsize'] = int(tilepx[-1][0]) + mapheader['xsize']
                            mergeheader['ysize'] = int(tilepx[-1][1]) + mapheader['ysize']
                            mergefile = mapfile
                    else:
                        merge_mrc = mrc.mmap(mergefile, permissive='True')
                        im = merge_mrc.data
                        mergeheader = map_header(merge_mrc)

                else:
                    mergeheader['xsize'] = int(tilepx[-1][0]) + mapheader['xsize']
                    mergeheader['ysize'] = int(tilepx[-1][1]) + mapheader['ysize']
                    mergefile = mapfile
        else:
            print('Warning: ' + mapfile + ' is not an MRC file!' + '\n')
            print('Assuming it is a single tif file or a stitched montage.' + '\n')
            mergefile = mapfile
            pixelsize = 1. / numpy.sqrt(abs(numpy.linalg.det(mat)))
            mergeheader['stacksize'] = 0
            tilepos = mapitem['StageXYZ'][0:2]
            tilepx = '0'
            tilepx = numpy.array([[tilepx, tilepx, tilepx], [tilepx, tilepx, tilepx]])
            tilepx1 = tilepx
            m['Sloppy'] = 'NoMont'
            im = io.imread(mergefile)
            mergeheader['xsize'] = numpy.array(im.shape)[0]
            mergeheader['ysize'] = numpy.array(im.shape)[1]
            mapheader = mergeheader.copy()

        mappxcenter = numpy.array([mergeheader['ysize'], mergeheader['xsize']]) / 2

    else:
        # map is mrc file

        mf = mrc.mmap(mapfile, permissive='True')

        mapheader = map_header(mf)
        pixelsize = mapheader['pixelsize']
        mergeheader = mapheader.copy()

        # determine if file contains multiple montages stored in one MRC stack

        mappxcenter = [mapheader['xsize'] / 2, mapheader['ysize'] / 2]

        mdocname = mapfile + '.mdoc'

        if m['frames'] == [0, 0]:
            mapheader['stacksize'] = 0
            tileidx_offset = 0
            # tilepos = [0, 0]
            if os.path.exists(mdocname):
                mdoclines = loadtext(mdocname)
                pixelsize = float(
                    mdoc_item(mdoclines, 'ZValue = ' + str(mapsection))['PixelSpacing'][0]) / 10000  # in um

        # extract center positions of individual map tiles
        if mapheader['stacksize'] > 1:

            if os.path.exists(mdocname):
                mdoclines = loadtext(mdocname)
                if mapheader['stacksize'] > montage_tiles:
                    # multiple same-dimension montages in one MRC stack
                    tileidx_offset = mapsection * montage_tiles
                elif mapheader['stacksize'] == montage_tiles:
                    # single montage without empty tiles (polygon)
                    tileidx_offset = 0
                elif mapheader['stacksize'] < montage_tiles:
                    # polygon fit montage with empty tiles
                    tileidx_offset = 0

                tilepos = list()
                tilepx = list()
                tilepx1 = list()
                for i in range(0, numpy.min([montage_tiles, mapheader['stacksize']])):
                    tile = mdoc_item(mdoclines, 'ZValue = ' + str(tileidx_offset + i))
                    tilepos.append(tile['StagePosition'])
                    tilepx1.append(tile['PieceCoordinates'])

                    if 'AlignedPieceCoordsVS' in tile:
                        m['Sloppy'] = True
                        tilepx.append(tile['AlignedPieceCoordsVS'])
                    elif 'AlignedPieceCoords' in tile:
                        tilepx.append(tile['AlignedPieceCoords'])
                    else:
                        tilepx = tilepx1

                if '[MontSection = 0]' not 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:]) / 10000  # in um
                else:
                    mont_item = mdoc_item(mdoclines, 'MontSection = ' + str(mapsection))
                    pixelsize = float(mont_item['PixelSpacing'][0]) / 10000  # in um

            # rotation = float(mdoc_item(mdoclines,'ZValue = '+str(mapsection))['RotationAngle'][0])

            else:
                if mapheader['stacksize'] > montage_tiles:
                    print(
                        'WARNING: Multiple maps stored in an MRC stack without mdoc file for metadata.'
                        ' I will guess the pixel size.')

                pixelsize = 1. / numpy.sqrt(abs(numpy.linalg.det(mat)))
                callcmd = 'extracttilts ' + mapfile + ' -stage -all'
                # os.system(callcmd)
                p2 = Popen(callcmd, shell=True, stdout=PIPE)
                o1 = list()
                for line in p2.stdout:
                    o1.append(line)

                tilepos1 = o1[20:-1]
                tilepos = numpy.array([numpy.fromstring(tilepos1[0], dtype=float, sep=' '),
                                       numpy.fromstring(tilepos1[1], dtype=float, sep=' ')])
                for item in tilepos1[2:]:
                    tilepos = numpy.append(tilepos,
                                           [numpy.fromstring(item, dtype=float, sep=' ')], axis=0)
                tilepx = '0'
                tilepx = numpy.array([[tilepx, tilepx, tilepx], [tilepx, tilepx, tilepx]])
                # tilepx1 = tilepx
                blendmont = True

        else:
            tilepos1 = list(map(float, mapitem['StageXYZ'][0:2]))
            tilepos = numpy.array([tilepos1, tilepos1])
            tilepx = '0'
            tilepx = numpy.array([[tilepx, tilepx, tilepx], [tilepx, tilepx, tilepx]])
            # tilepx1 = tilepx

        # check if map is a montage or not

        if mapheader['stacksize'] < 2:
            print('Single image found. No merging needed.')
            tilepx = '0'
            tilepx = numpy.array([[tilepx, tilepx, mapsection], [tilepx, tilepx, mapsection]])
            tilepx1 = tilepx
            # os.system(callcmd)
            merge_mrc = mf
            mergeheader = mapheader.copy()
            # overlapx = 0
            # overlapy = 0
            # tileloc = [0, 0]
            imd = merge_mrc.data
            mergefile = mapfile

            if len(imd.shape) == 3:
                im = imd[mapsection, :, :]
            elif len(imd.shape) == 2:
                im = imd

        else:

            if blendmont:
                mergefile = mergebase + '.mrc'
                if not os.path.exists(mergefile):
                    call_blendmont(mapfile, mergebase, mapsection, black)
                mergebase = os.path.splitext(mapfile)[0] + '_merged' + '_s' + str(mapsection)
                mergefile = mergebase + '.mrc'

                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(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
                tilepx1 = list(loadtext(mapfile + '.pcs'))

                for j, item in enumerate(tilepx1):
                    tilepx1[j] = list(re.split(' +', item))
                merge_mrc.close()

            else:
                mergefile = mapfile
                im = mf.data
                tilepx = '0'
                tilepx = numpy.array([[tilepx, tilepx, tilepx], [tilepx, tilepx, tilepx]])

                if not os.path.exists(mdocname):
                    callcmd = 'extractpieces ' + mapfile + ' -ou ' + mapfile + '.pcs'
                    os.system(callcmd)
                    pxpos = loadtext(mapfile + '.pcs')
                    tilepx = [re.split(' +', line) for line in pxpos]

                tilepx1 = tilepx

                mergeheader['xsize'] = int(tilepx[-1][0]) + im.shape[0]  # int(mapitem['MapWidthHeight'][0])
                mergeheader['ysize'] = int(tilepx[-1][1]) + im.shape[1]  # int(mapitem['MapWidthHeight'][1])

        if not blend_in:
            im = mf.data

        im = numpy.rot90(numpy.transpose(im), axes=(0, 1))
        mf.close()
        # end MRC section

    tilepos = numpy.array(tilepos, float)

    tilepx = numpy.array(tilepx)
    tilepx = tilepx[tilepx[:, 2] == str(mapsection), 0:2]
    tilepx = tilepx.astype(float)

    tilepx1 = numpy.array(tilepx1)
    tilepx1 = tilepx1[tilepx1[:, 2] == str(mapsection), 0:2]
    tilepx1 = tilepx1.astype(float)

    tpx = tilepx1[:, 0]
    tpy = tilepx1[:, 1]

    if numpy.abs(tpx).max() > 0:
        xstep = tpx[tpx > 0].min()
    else:
        xstep = 1
    if numpy.abs(tpy).max() > 0:
        ystep = tpy[tpy > 0].min()
    else:
        ystep = 1

    tileloc = numpy.array([tpx / xstep, tpy / ystep]).T

    if not blendmont and bigstitch:
        # prepare coordinate list for Big Stitcher
        print('preparing coordinate list for BigStitcher for map item ' + mapitem['# Item'] + '.')
        outpx = tilepx.copy()

        stitchname = mapfile + '.stitch.csv'
        stitchfile = open(stitchname, 'w')

        stitchfile.write('dim=2\n')
        # stitchfile.write('ViewSetupID;TimePointID;(position_x, position_y, position_z)\n')

        for j, item in enumerate(outpx):
            stitchfile.write(str(j) + ";;" + "(%s, %s" % (int(item[0]), -int(item[1])) + ")\n")

        stitchfile.close()

    m['sections'] = numpy.array(list(map(int, tileloc[:, 0] * m['frames'][1] + tileloc[:, 1])))

    overlapx = mapheader['xsize'] - xstep
    overlapy = mapheader['ysize'] - ystep

    # 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(mergebase + '_crop.mrc'):
            loopcount = 0
            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(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: \"' + mergebase + '.mod\".')
                time.sleep(20)
                loopcount = loopcount + 1

            if loopcount < 21:
                print('Model file found. Will now generate the cropped map image.')
                callcmd = 'imodmop \"' + mergebase + '.mod\" \"' + mergefile + '\" \"' + mergebase + '_crop.mrc\"'
                os.system(callcmd)

        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

        merge_mrc.close()

    mergeheader['pixelsize'] = pixelsize
    mapheader['pixelsize'] = pixelsize

    # generate output

    m['mapfile'] = mapfile
    m['mergefile'] = mergefile
    m['matrix'] = mat
    m['tilepos'] = tilepos
    m['im'] = im
    m['mappxcenter'] = mappxcenter
    m['mapheader'] = mapheader
    m['mergeheader'] = mergeheader
    m['tilepx'] = tilepx
    m['tilepx1'] = tilepx1
    m['overlap'] = [overlapx, overlapy]
    m['tileloc'] = tileloc
    m['tileidx_offset'] = tileidx_offset

    return m


# -------------------------------
# %%
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 == 4 & imod_ver[1] < 10):
        raise EnvironmentError('ERROR: IMOD version needs to be at least 4.11 ! Please update. Exiting' + '\n')

    # merge the montage to a single file
    callcmd = 'extractpieces ' + '\"' + mapfile + '\" \"' + mapfile + '.pcs\"'
    os.system(callcmd)
    # print(callcmd)

    print('----------------------------------------------------\n')
    print('Merging the map montage into a single image....' + '\n')
    print('----------------------------------------------------\n')

    callcmd = 'blendmont -v -rob 1 -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'

    os.system(callcmd)


# -------------------------------

# %%

def realign_map(item, allitems):
    """
    Determines which map to align to for given navigator item.

    Parameters
    ----------
    item : dict
        Navigator Item
    allitems : list of dict
        list of all Navigator Items
    Returns
    -------
    dict

    """
    #

    if item['Type'] in [['0'], ['1']]:
        # point or polygon
        if 'DrawnID' not in item.keys():
            if 'SamePosId' not in item.keys():
                print('No map found to realign item ' + item['# Item'] + ' to, skipping it...')
                result = []
            else:
                mapitem = list(filter(lambda c_item: c_item['MapID'] == item['SamePosId'], allitems))[0]
                result = realign_map(mapitem, allitems)
        else:
            mapID = item['DrawnID']
            result = list(filter(lambda item: item['MapID'] == mapID, allitems))

    else:

        if 'RealignedID' not in item.keys():
            print('No map found to realign item ' + item['# Item'] + ' to, skipping it...')
            return []

        else:
            mapID = item['RealignedID']

        result = list(filter(lambda item: item['MapID'] == mapID, allitems))

    return result[0]


# -------------------------------------


def imcrop(im1, c, sz):
    """
    Crops an image of a given size (2 element numpy array) around a pixel coordinate (2 element numpy array)
    in case the coordinate is close to an edge, the cropped image will have the maximum possible width/height
    and centering of the image.
    
    Parameters
    ----------
    im1 : numpy.ndarray or list
        input image (2D).
    c : numpy.ndarray or list
        center around which to crop (2x1).
    sz : numpy.ndarray or list
        (2x1) size of the crop region. (maximum, depending on image boundaries)
    
    Returns
    -------
    im2 : numpy.ndarray
         (2D) output cropped image.
    
    """

    if type(im1) is list:
        im1 = numpy.array(im1)

    if type(c) is list:
        c = numpy.array(c)

    if type(sz) is list:
        sz = numpy.array(sz)

    if any([type(p) is not numpy.ndarray for p in [im1, c, sz]]):
        raise TypeError('Inputs need to be list or numpy array.')

    sz_x = sz[0]
    sz_y = sz[1]

    ximsz = im1.shape[0]
    yimsz = im1.shape[1]

    xllim = max([0, c[1] - sz_x / 2])
    xulim = min([ximsz, c[1] + sz_x / 2])

    x_range = min([c[1] - xllim, xulim - c[1]])
    xel = range(int(c[1] - x_range), int(c[1] + x_range))

    yllim = max([0, c[0] - sz_y / 2])
    yulim = min([yimsz, c[0] + sz_y / 2])

    y_range = min([c[0] - yllim, yulim - c[0]])

    yel = range(int(c[0] - y_range), int(c[0] + y_range))

    im2 = im1[xel, :]
    im2 = im2[:, yel]

    return im2


# --------------------------------------
# %%

def cart2pol(c):
    """
    Converts Cartesian into polar coordinates (2D)

    Parameters
    ----------
    c : numpy.ndarray or list

    Returns
    -------
    numpy.ndarray

    """
    #
    rho = numpy.sqrt(c[:, 0] ** 2 + c[:, 1] ** 2)
    phi = numpy.arctan2(c[:, 1], c[:, 0])
    return numpy.transpose([phi, rho])


# --------------------------------------

def pol2cart(rho, phi):
    """
    Converts  polar into cartesian coordinates (2D)

    Parameters
    ----------
    rho : float or int
    phi : float or int

    Returns
    -------
    numpy.ndarray

    """

    x = rho * numpy.cos(phi)
    y = rho * numpy.sin(phi)
    return numpy.transpose([x, y])


# --------------------------------------
# %%

def img2polygon(img, n_poly, center, radius):
    """
    Converts a binary image into a polygon (list of points) describing its outline


    Parameters
    ----------
    img : numpy.ndarray
        a binary image
    n_poly : int
        number of points for polygon
    center : list or numpy.ndarray
        center around which to draw the polygon
    radius : int or float
        maximum radius of polygon (pixel)

    Returns
    -------
    numpy.ndarray

    """

    if type(img) is not numpy.ndarray:
        raise TypeError('Input image needs to be numpy array.')

    # define threshold based on image type
    if img.dtype.kind == 'b':
        thresh = 1
    elif img.dtype.kind == 'i':
        thresh = 2 ** (8 * img.dtype.itemsize - 1) - 1
        if img.max() < thresh:
            thresh = img.max() / 2
    else:
        thresh = img.max() / 2

    # close the polygon
    n_poly = n_poly + 1

    xs, ys = img.shape

    polypt = numpy.empty((0, 2))

    # draw equi-angular lines from the center to determine the intersection with the binary image
    polyphi = numpy.linspace(0, 2 * numpy.pi, n_poly)
    endpts = pol2cart(radius, polyphi)
    endpts = numpy.array(center) + endpts

    # find the intersection points and mark the coordinates
    for pt in endpts:
        x, y = numpy.linspace(center[0], pt[0], radius), numpy.linspace(center[1], pt[1], radius)
        a, b = x.astype(int), y.astype(int)
        a[a > (ys - 1)] = ys - 1
        a[a < 0] = 0
        b[b > (xs - 1)] = xs - 1
        b[b < 0] = 0
        lpx = img[b, a].astype(int)
        lpd = numpy.diff(lpx)
        maxdiff = numpy.max(numpy.abs(lpd))
        maxdiff_ix = numpy.argmax(numpy.abs(lpd))
        if maxdiff < thresh:
            maxdiff_ix = radius - 1

        polypt = numpy.append(polypt, [(a[maxdiff_ix], b[maxdiff_ix])], axis=0)

    return polypt


# --------------------------------------

def map_extract(im, c, p, px_scale, t_size, mat, int8=False):
    """
    Extracts an image from a given position in an existing map and links positions inside.

    Parameters
    ----------
    im : numpy.ndarray
        input image.
    c : numpy.ndarray
        (2x1) - center of target region.
    p : list of numpy.ndarray
        (2xn) - points in the image that are transformed.
    px_scale : float
        relative scale of input and output image
    t_size : numpy.ndarray
        (2x1) - target size of output image (maximum, limited by original region boundaries).
    mat : numpy.mat or numpy.ndarray
        (2x2) - transformation matrix.
    int8 : bool
        if conversion of the outoput image to 8bit is desired. The default is False.

    Returns
    -------
    numpy.ndarray, list of numpy.ndarray
     img3 (2D) output extracted image.

     p4 (2xn) output pixel coordinates of input points

    """

    if type(t_size) is list:
        t_size = numpy.array(t_size)

    if type(p) is list:
        p = numpy.array(p)
        if len(p.shape) < 2:
            p = [p]

    if type(c) is list:
        c = numpy.array(c)

    imsz1 = t_size * px_scale

    # extract image (1.42x to enable rotation)
    cropsize1 = imsz1 * 1.42

    cropsize = numpy.array([cropsize1.max(), cropsize1.max()])

    im1 = imcrop(im, c, cropsize)
    realsize = numpy.array(im1.shape)

    # center to origin
    p1 = p - c

    # create homogenous matrices
    mat_i = numpy.linalg.inv(mat)

    o_size = numpy.max(numpy.abs([numpy.dot([0, realsize[1]], mat), numpy.dot([realsize[0], 0], mat)]
                                 ), axis=0).astype(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(abs(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))).astype(int)

    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 = numpy.dot(p1, mat.T)

    p4[:, 0] = f_size[1] / 2 + p4[:, 0]
    p4[:, 1] = f_size[0] / 2 + p4[:, 1]

    return im3, p4


# --------------------------------------

def get_pixel(navitem, mergedmap, tile=False, outline=False):
    """
    Determines the pixel coordinates of a navigator item in its associated map.

    Parameters
    ----------
    navitem : dict
        Navigator Item
    mergedmap : dict
        Navigator Item  as resulting from mergemap
    tile : bool
        if True: output is pixel and tile index of closest map tile
    outline : bool
        if True returns coordinates of outline (map/polygon) instead of center point

    Returns
    -------
    numpy.ndarray
        array containing the pixel coordinates and optionally the tile index
    """

    xval = float(navitem['StageXYZ'][0])  # (float(acq_item['PtsX'][0]))
    yval = float(navitem['StageXYZ'][1])  # (float(acq_item['PtsY'][0]))

    pt0 = numpy.array([xval, yval])

    pt = pt0.copy()
    # calculate the pixel coordinates

    imsz = [mergedmap['mergeheader']['ysize'], mergedmap['mergeheader']['xsize']]

    if outline:
        xval1 = numpy.array(navitem['PtsX']).astype(float)
        yval1 = numpy.array(navitem['PtsY']).astype(float)

        pt = numpy.vstack((xval1, yval1)).T

    if 'XYinPc' in navitem:
        tileid = int(navitem['PieceOn'][0])
        pt_px0 = list(map(float, navitem['XYinPc']))
        pt_px = numpy.array(pt_px0)
        pt_px = numpy.reshape(pt_px, (1, 2))
        tileidx = numpy.argwhere(mergedmap['sections'] == tileid)[0][0]

    else:

        tilepos = mergedmap['tilepos']

        if (numpy.diff(tilepos, axis=0)[0].max() == 0) and (type(mergedmap['Sloppy']) == bool) and (
                mergedmap['mapheader']['stacksize'] > 1):
            print(
                'Montage created using image shift! Problems in identifying'
                ' the positions of clicked points accurately possible!')

        if len(tilepos.shape) < 2:
            tileidx = 0
        else:
            if 'PieceOn' in navitem:
                tileid = int(navitem['PieceOn'][0])
                tileidx = numpy.argwhere(mergedmap['sections'] == tileid)[0][0]
            else:
                tiledist = numpy.sum((tilepos - pt0) ** 2, axis=1)
                tileidx = numpy.argmin(tiledist)

        # normalize coordinates

        ptn = numpy.array(pt - tilepos[tileidx])

        pt_px = numpy.dot(ptn, mergedmap['matrix'].T)

        if len(pt_px.shape) > 1:
            pt_px[:, 0] = (mergedmap['mappxcenter'][0]) + pt_px[:, 0]
            pt_px[:, 1] = (mergedmap['mappxcenter'][1]) - pt_px[:, 1]
        else:
            pt_px[0] = (mergedmap['mappxcenter'][0]) + pt_px[0]
            pt_px[1] = (mergedmap['mappxcenter'][1]) - pt_px[1]

    # output
    if not outline:
        pt_px = pt_px.squeeze()

    if tile:
        return numpy.array([pt_px, tileidx])
    else:
        if len(pt_px.shape) > 1:
            pt_px[:, 1] = mergedmap['mapheader']['xsize'] - pt_px[:, 1]
        else:
            pt_px[1] = mergedmap['mapheader']['xsize'] - pt_px[1]

        pt_px1 = pt_px + mergedmap['tilepx'][tileidx]

        return pt_px1


# ------------------------------------------------------------

def ptsonmap(mapitem, pts, nav):
    # adds a point or polygon (pixel position) to a map and creates a navigator item from
    # p is a list of numpy

    outnav = list()

    mapid = mapitem['MapID']

    for idx, p in enumerate(pts):

        p = numpy.array(p)

        polynav = dict()

        merge = mergemap(mapitem)

        imsz1 = numpy.array([merge['mergeheader']['xsize'], merge['mergeheader']['ysize']])

        p[:, 1] = imsz1[0] - p[:, 1]

        c = numpy.average(p, axis=0)

        px = numpy.array(numpy.transpose(p[:, 0]))
        px = numpy.array2string(px.astype(int), separator=' ')
        px = px.strip('"[]"')

        py = numpy.array(numpy.transpose(p[:, 1]))
        py = numpy.array2string(py.astype(int), separator=' ')
        py = py.strip('"[]"')

        if numpy.shape(p)[0] == 1:
            polynav['Type'] = ['0']
            polynav['Color'] = ['0']
            polynav['NumPts'] = ['1']

        else:
            polynav['Type'] = ['1']
            polynav['Color'] = ['1']
            polynav['NumPts'] = [str(p.shape[0])]

        label = mapitem['# Item'] + '_' + str(idx).zfill(3)

        #  map corner points
        # fill navigator

        # Polygon
        polynav['# Item'] = label
        polynav['Acquire'] = ['1']

        polynav['Draw'] = ['1']
        polynav['Regis'] = mapitem['Regis']
        polynav['DrawnID'] = mapid

        if mapitem['MapFramesXY'] == ['0', '0']:
            polynav['CoordsInMap'] = [str(int(round(c[0]))), str(int(round(c[1]))), mapitem['StageXYZ'][2]]
        else:
            tilecenters = merge['tilepx'] + numpy.array(
                [merge['mapheader']['xsize'] / 2, merge['mapheader']['ysize'] / 2])
            c_out = c.copy()
            c_out[1] = imsz1[0] - c_out[1]
            tiledist = numpy.sum((tilecenters - c_out) ** 2, axis=1)

            tileidx = numpy.argmin(tiledist)
            c_out = c_out - merge['tilepx'][tileidx]
            polynav['CoordsInPiece'] = [str(c_out[0]), str(c_out[1]), mapitem['StageXYZ'][2]]
            polynav['PieceOn'] = [str(merge['sections'][tileidx])]

        polynav['PtsX'] = px.split()
        polynav['PtsY'] = py.split()
        polynav['GroupID'] = [str(newID(nav, int(mapid[0]) + 70000))]

        outnav.append(polynav)

    return outnav


# ---------------------------------------------------------------------


def pts2nav(im, pts, cntrs, curr_map, targetitem, nav, sloppy=False, maps=False):
    # takes pixel coordinates on an input map, generates virtual maps at their positions,
    # creates polygons matching their shape.
    # input parameters:
    # - im, original image for feature detection
    # - pts, image coordinates, list of arrays (for polygons)
    # - cntrs, center of the feature (bounding box)
    # - curr_map, map item to process (emtools dict)
    # - targetitem, target map
    # - nav, all navigator items (emtools dict)
    # - sloppy (optional) if merging of original map was sloppy
    # - maps (optional) if the virtual maps instead of the points/polygons should be selected for acquisition

    # parse input data

    if type(im) != numpy.ndarray:
        raise Exception('Wrong input format of image.')

    if type(pts) != list:
        if type(pts) == numpy.ndarray:
            pts = [pts]
        else:
            raise Exception('Wrong input format of point coordinates.')

    if type(cntrs) != list:
        if type(cntrs) == numpy.ndarray:
            cntrs = [cntrs]
        else:
            raise Exception('Wrong input format of center coordinates.')
    else:
        if len(cntrs) == 2:
            if type(cntrs[0]) == int:
                cntrs = [cntrs]
            elif (type(cntrs[0]) != numpy.ndarray) and (type(cntrs[0]) != list):
                raise Exception('Wrong input format of center coordinates.')

    # generate output
    outnav = list()
    nav_pol = list()
    nav_maps = list()

    # read information from maps

    merge = mergemap(curr_map)

    # mapfile = merge['mapfile']
    # map_mrc = mrc.mmap(mapfile, permissive = 'True')
    mapheader = merge['mapheader']

    pixelsize = mapheader['pixelsize']

    imsz = numpy.array(im.shape)

    map_mat = merge['matrix']

    # target reference

    targetheader = map_header(targetitem)

    t_mat = map_matrix(targetitem)

    # combine rotation matrices
    maptf = (numpy.linalg.inv(map_mat) * t_mat).T

    px_scale = targetheader['pixelsize'] / pixelsize

    imsz1 = numpy.array([targetheader['ysize'], targetheader['xsize']])

    ntotal = len(cntrs)

    outnav.append(curr_map)

    curr_id = int(curr_map['MapID'][0])

    delim = 100000

    startid = newID(nav, divmod(curr_id, delim)[0] * delim)

    mapid = startid + 1
    idx = 0

    for c in cntrs:

        polynav = dict()
        newnavitem = dict(targetitem)
        mapid = newID(nav, mapid + 1)

        print('Processing object ' + str(idx + 1) + '/' + str(ntotal) + ' (%2.0f%%)' % (
                idx * 100 / ntotal) + ' at position %5u , %5u' % (c[0], c[1]))

        p = pts[idx]

        im4, p4 = map_extract(im, c, p, px_scale, imsz1, maptf)

        if min(im4.shape) < 400:
            print('Item is too close to border of map, skipping it.')
            ntotal = ntotal - 1
            continue

        p4[:, 1] = imsz1[0] - p4[:, 1]

        px = numpy.array(numpy.transpose(p4[:, 0]))
        px = numpy.array2string(px, separator=' ')
        px = px.strip('"[]"')

        py = numpy.array(numpy.transpose(p4[:, 1]))
        py = numpy.array2string(py, separator=' ')
        py = py.strip('"[]"')

        if numpy.shape(p4)[0] == 1:
            polynav['Type'] = ['0']
            polynav['Color'] = ['0']
            polynav['NumPts'] = ['1']

        else:
            polynav['Type'] = ['1']
            polynav['Color'] = ['1']
            polynav['NumPts'] = [str(p.shape[0])]

        label = curr_map['# Item'] + '_' + str(idx).zfill(3)
        idx = idx + 1

        imfile = 'virt_' + label + '.mrc'

        if os.path.exists(imfile):
            os.remove(imfile)

        im5 = numpy.rot90(im4, 3)

        with mrc.new(imfile) as mrcf:
            mrcf.set_data(im5.T)
            mrcf.close()  #

        #  map corner points

        cx = imsz1[1]
        cy = imsz1[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(maptf)

        c1 = a * t_mat_i

        c_out = c.copy()
        c_out[1] = imsz[0] - c_out[1]

        # c1 = c1 + [c_out[0],c_out[1]]

        # fill navigator

        # map for realignment

        newnavitem['MapFile'] = [imfile]
        newnavitem.pop('StageXYZ', '')
        newnavitem.pop('RawStageXY', '')
        if curr_map['MapFramesXY'] == ['0', '0']:
            newnavitem['CoordsInMap'] = [str(c_out[0]), str(c_out[1]), curr_map['StageXYZ'][2]]

        else:
            #  newnavitem['CoordsInAliMont'] = [str(c_out[0]),str(c_out[1]),curr_map['StageXYZ'][2]]
            #  convert aligned pixel coordinates into piece coordinates to ensure map has a proper bounding box

            tilecenters = merge['tilepx'] + numpy.array(
                [merge['mapheader']['xsize'] / 2, merge['mapheader']['ysize'] / 2])
            tiledist = numpy.sum((tilecenters - c_out) ** 2, axis=1)
            tileidx = numpy.argmin(tiledist)
            c_out = c_out - merge['tilepx'][tileidx]
            c1 = numpy.fliplr(c1 + c_out)
            newnavitem['CoordsInPiece'] = [str(c_out[0]), str(c_out[1]), curr_map['StageXYZ'][2]]
            newnavitem['PieceOn'] = [str(merge['sections'][tileidx])]

        cnx = numpy.array(numpy.transpose(c1[:, 1]))
        cnx = numpy.array2string(cnx, separator=' ')
        cnx = cnx[2:-2]

        cny = numpy.array(numpy.transpose(c1[:, 0]))
        cny = numpy.array2string(cny, separator=' ')
        cny = cny[2:-2]

        newnavitem['PtsX'] = cnx.split()
        newnavitem['PtsY'] = cny.split()
        newnavitem['Note'] = newnavitem['MapFile']
        newnavitem['MapID'] = [str(mapid)]
        newnavitem['DrawnID'] = curr_map['MapID']

        if maps:
            newnavitem['Acquire'] = ['1']
        else:
            newnavitem['Acquire'] = ['0']

        newnavitem['MapSection'] = ['0']
        newnavitem.pop('SamePosId', '')
        # newnavitem['MapWidthHeight'] = [str(im2size[0]),str(im2size[1])]
        newnavitem['ImageType'] = ['0']
        newnavitem['MapMinMaxScale'] = [str(numpy.min(im4)), str(numpy.max(im4))]
        newnavitem['NumPts'] = ['5']
        newnavitem['# Item'] = 'm_' + label
        newnavitem['GroupID'] = [str(newID(nav, startid + 50000))]
        curr_map['Acquire'] = ['0']

        # Polygon
        nav_maps.append(newnavitem)

        # do not export points/polygons when maps are inteded output.

        if not maps:
            polynav['# Item'] = label
            polynav['Acquire'] = ['1']

            polynav['Draw'] = ['1']
            polynav['Regis'] = curr_map['Regis']
            polynav['DrawnID'] = [str(mapid)]
            polynav['CoordsInMap'] = [str(int(cx / 2)), str(int(cy / 2)), curr_map['StageXYZ'][2]]
            polynav['PtsX'] = px.split()
            polynav['PtsY'] = py.split()
            polynav['GroupID'] = [str(newID(nav, startid + 70000))]

            nav_pol.append(polynav)

    outnav.extend(nav_maps)
    outnav.extend(nav_pol)

    return outnav


# ------------------------------------------------------------

def nav_find(allitems, key, val=None):
    # returns the navigator/mdoc entries with the given key/value pair
    # takes an input navigator (list of dicts), a target key (item property) and the desired values to match
    # output is the list of navigator items

    if val is None:
        val = []
    filtered = list(filter(lambda item: item.get(key), allitems))

    if val == []:
        found = filtered
    else:
        # deals with integer, string or list entries
        if type(val) == int:
            val = str(val)
        if key != '# Item':
            if type(val) == str:
                val = [val]

        found = list(filter(lambda item: item[key] == val, filtered))

    if not found == []:
        newnav = found
    else:
        newnav = []

    return newnav


# ------------------------------------------------------------

def nav_selection(allitems, sel=None, acquire=True, maps=False):
    # extracts a selection of navigator items into a new navigator
    # takes an input navigator (list of dicts), an optional list of item labels (one line each)
    # and optional whether to include items selected for acquisition, optionally,
    # all maps the selected points were clicked on are also returned
    # if only the input nav is given, it will extract all "Acquire" items
    # output is list of navigator items

    if sel is None:
        sel = []
    newnav = list()
    select = sel[:]

    if acquire:
        newnav = nav_find(allitems, 'Acquire', '1')

    if not (select == []):
        if isinstance(select, str):
            select = [select]
        for listitem in select:
            selitem = nav_find(allitems, '# Item', listitem)

            newnav.extend(selitem)

    if maps:
        for item in newnav:
            if 'DrawnID' in item.keys():
                drawnmap = nav_find(allitems, 'MapID', item['DrawnID'])
                if nav_find(newnav, 'MapID', item['DrawnID']) == []:
                    newnav.extend(drawnmap)

    return newnav


# ------------------------------------------------------------

def outline2mod(im, namebase, z=0, binning=1):
    # takes an input image of label outlines (single pixel thickness)
    # creates an IMOD model file with these outlines as contours.
    #

    filename = namebase + '.txt'
    f = open(filename, 'w')

    for label in numpy.unique(im[im > 0]):
        a = numpy.argwhere(im == label)

        # initialize geometrical sorting, working vars
        am = a
        an = a[0, :]
        pt = an
        am = numpy.delete(am, 0, 0)
        jump = 0

        while am.shape[0] > 2:
            # find neighbouring point from remaining ones
            dist = (am[:, 0] - pt[0]) ** 2 + (am[:, 1] - pt[1]) ** 2
            ix = numpy.argmin(dist)

            if numpy.min(dist) > 50:
                if jump:
                    an = numpy.vstack((an, anchor))
                    pt = anchor
                    jump = 0
                else:
                    anchor = pt
                    jump = 1
                    an = numpy.vstack((an, am[ix, :]))
                    pt = am[ix, :]
                    am = numpy.delete(am, ix, 0)
            else:
                # add to sorted list
                an = numpy.vstack((an, am[ix, :]))
                pt = am[ix, :]

                # remove this point from input
                am = numpy.delete(am, ix, 0)

        # add last 2 points to sorted list
        ix = numpy.argmin((am[:, 0] - pt[0]) ** 2 + (am[:, 1] - pt[1]) ** 2)
        an = numpy.vstack((an, am[ix, :]))
        am = numpy.delete(am, ix, 0)
        an = numpy.vstack((an, am[0, :]))

        points = numpy.vstack((an[:, 1], im.shape[0] - an[:, 0])).transpose() * binning

        # write output file
        for j, item in enumerate(points):
            f.write(" 1  " + str(label) + "  %s %s" % (item[0], item[1]) + " " + str(z) + "\n")

    f.close()

    # convert into IMOD model
    callcmd = 'point2model \"' + filename + '\" \"' + namebase + '.mod\"'
    os.system(callcmd)


# ------------------------------------------------------------

def ordernav(nav, delim=''):
    # re-orders a navigator by its label.
    # It considers the indexing after a delimiter in the string.
    # example: s01_cell-1,s02_cell-1,s01_cell-02, ...    is sorted by cells instead of s...
    # when no delimiter is given, the navigator is sorted by its label.

    nav1 = nav.copy()

    non_idx = 0

    for item in nav1:
        if delim == '':
            item['# Sorting'] = item['# Item']

        elif item['# Item'].find(delim) == -1:
            item['# Sorting'] = str(non_idx).zfill(5)
            non_idx = non_idx + 1
        else:
            item['# Sorting'] = item['# Item'][item['# Item'].find(delim) + 1:]

    newnav = sorted(nav1.copy(), key=itemgetter('# Sorting'))

    for item in newnav:
        item.pop('# Sorting')

    return newnav


# ------------------------------------------------------------

def pointitem(label, regis=1):
    # creates a default point item containing all necessary item parameters.
    # The corresponding MapID and Position need to be defined externally!

    point = dict()
    point['# Item'] = label
    point['Color'] = ['0']

    if type(regis) == list:
        regis = ''.join(map(str, regis))

    point['Regis'] = [str(regis)]
    point['NumPts'] = ['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.
    maps: dict
        maps dictionary
    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.
    maps: dict
        maps dictionary
    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 itemid not 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