# AtrialFibrePlugin
# Copyright (C) 2018 <NAME>, King's College London, all rights reserved, see LICENSE file

'''
Atrial fibre generation plugin.
'''

import os
import stat
import ast
import shutil
import datetime
import zipfile
import warnings
from itertools import starmap
from collections import defaultdict

try:
    import configparser
except ImportError:
    import ConfigParser as configparser
    
try:
    from sfepy.base.conf import ProblemConf
    from sfepy.applications import solve_pde
    from sfepy.base.base import output
except ImportError:
    warnings.warn('SfePy needs to be installed or in PYTHONPATH to generate fiber directions.')

from eidolon import ( ui,
        ScenePlugin, Project, avg, vec3, successive, first, RealMatrix, IndexMatrix, StdProps, timing, ReprType,
        listToMatrix,MeshSceneObject, BoundBox, ElemType, reduceMesh, PyDataSet, taskmethod, taskroutine, printFlush
        )
import eidolon

import numpy as np
from scipy.spatial import cKDTree

plugindir= os.path.dirname(os.path.abspath(__file__)) # this file's directory

# directory/file names
uifile=os.path.join(plugindir,'AtrialFibrePlugin.ui') 
deformDir=os.path.join(plugindir,'deformetricaC')
deformExe=os.path.join(deformDir,'deformetrica')
architecture=os.path.join(plugindir,'architecture.ini') 
problemFile=os.path.join(plugindir,'problemfile.py') 

# registration file names
decimatedFile='subject.vtk'
targetFile='target.vtk'
datasetFile='data_set.xml'
modelFile='model.xml'
optimFile='optimization_parameters.xml'
registeredFile='Registration_subject_to_subject_0__t_9.vtk'
decimate='decimate-surface'

# deformetrica parameters
kernelWidthSub=5000
kernelWidthDef=8000
kernelType='cudaexact'
dataSigma=0.1
stepSize=0.000001

# field names
#regionField='regions'
#landmarkField='landmarks'
#directionField='directions'
#gradientDirField='gradientDirs'
#elemDirField='elemDirs'
#elemRegionField='elemRegions'
#elemthickness='elemThickness'
#elemGradient='elemGradient'

fieldNames=eidolon.enum(
    'regions','landmarks',
    'directions','gradientDirs',
    'elemDirs','elemRegions','elemThickness', 
    'nodeGradient'
)

objNames=eidolon.enum(
    'atlasmesh',
    'origmesh',
    'epimesh','epinodes',
    'endomesh','endonodes',
    'architecture'
)

regTypes=eidolon.enum('endo','epi')


# load the UI file into the ui namespace, this is subtyped below
ui.loadUI(open(uifile).read())


def showLines(nodes,lines,name='Lines',matname='Default'):
    mgr=eidolon.getSceneMgr()
    lineds=eidolon.LineDataSet(name+'DS',nodes,lines)
    obj=MeshSceneObject(name,lineds)
    mgr.addSceneObject(obj)
    
    rep=obj.createRepr(eidolon.ReprType._line,matname=matname)
    mgr.addSceneObjectRepr(rep)
    
    return obj,rep


def showElemDirs(obj,glyphscale,mgr):
    ds=obj.datasets[0]
    nodes=ds.getNodes()
    tets=first(i for i in ds.enumIndexSets() if i.getType()==ElemType._Tet1NL)
    elemdirfield=ds.getDataField(fieldNames._elemDirs)
    
    mid=ElemType.Tet1NL.basis(0.25,0.25,0.25)
    
    elemnodes=[ElemType.Tet1NL.applyCoeffs([nodes[e] for e in elem],mid) for elem in tets]
    elemobj=MeshSceneObject('elemobj',PyDataSet('elemobjds',elemnodes,[],[elemdirfield]))
    mgr.addSceneObject(elemobj)
    
    rep=elemobj.createRepr(eidolon.ReprType._glyph,0,externalOnly=False,drawInternal=True,glyphname='arrow',
                           glyphscale=glyphscale,dfield=elemdirfield.getName(),vecfunc=eidolon.VecFunc._Linear)
    
    mgr.addSceneObjectRepr(rep)


class initdict(defaultdict):
    def __init__(self,initfunc,*args,**kwargs):
        defaultdict.__init__(self,None,*args,**kwargs)
        self.initfunc=initfunc
        
    def __missing__(self,key):
        value=self.initfunc(key)
        self.__setitem__(key,value)
        return value
    

class plane(object):
    def __init__(self,center,norm):
        self.center=center
        self.norm=norm.norm()
        
    def dist(self,pt):
        return pt.planeDist(self.center,self.norm)
        
    def moveUp(self,dist):
        self.center+=self.norm*dist
        
    def numPointsAbove(self,nodes):
        return sum(1 for n in nodes if self.dist(n)>=0)
        
    def between(self,nodes,otherplane):
        numnodes=len(nodes)
        return self.numPointsAbove(nodes)==numnodes and otherplane.numPointsAbove(nodes)==numnodes
       
    def findIntersects(self,nodes,inds):
        numnodes=inds.m()
        result=[]        
        for n in range(inds.n()):
            if 0<self.numPointsAbove(nodes.mapIndexRow(inds,n))<numnodes:
                result.append(n)
                
        return result
        
        
class TriMeshGraph(object):
    def __init__(self,nodes,tris,ocdepth=3):
        self.nodes=nodes if isinstance(nodes,eidolon.Vec3Matrix) else listToMatrix(nodes,'nodes')
        self.tris=tris if isinstance(tris,eidolon.IndexMatrix) else listToMatrix(tris,'tris')
        
        self.tricenters=[avg(self.getTriNodes(r),vec3()) for r in range(self.tris.n())]
        self.adj,self.ragged=generateTriAdj(self.tris) # elem -> elems
        self.nodeelem=generateNodeElemMap(self.nodes.n(),self.tris) # node -> elems
        self.edges=generateSimplexEdgeMap(self.nodes.n(),self.tris) # node -> nodes
        self.boundbox=BoundBox(nodes)
        
#        self.octree=eidolon.Octree(ocdepth,self.boundbox.getDimensions(),self.boundbox.center)
#        self.octree.addMesh(self.nodes,self.tris)
        
        def computeDist(key):
            i,j=key
            return self.tricenters[i].distTo(self.tricenters[j])
            
        self.tridists=initdict(computeDist)
        
#    def getIntersectedTri(self,start,end):
#        '''
#        Returns the triangle index and the (t,u,v) triple if the line from `start' to `end' intersects the indexed triangle 
#        at a distance of `t' from `start' with xi coord (u,v). Returns None if no triangle intersected.
#        '''
#        startoc=self.octree.getLeaf(start)
#        endoc=self.octree.getLeaf(end)
#        inds=(startoc.leafdata if startoc is not None else []) + (endoc.leafdata if endoc is not None else []) 
#        
#        r=eidolon.Ray(start,end-start)
#        
#        for tri in inds:
#            d=r.intersectsTri(*self.getTriNodes(tri))
#            if d:# and d[0]<=start.distTo(end):
#                return tri,d
#            
#        return None
        
    def getPathSubgraph(self,starttri,endtri):
        return getAdjTo(self.adj,starttri,endtri)
        
    def getTriNodes(self,triindex):
        return self.nodes.mapIndexRow(self.tris,triindex)
    
    def getTriNorm(self,triindex):
        a,b,c=self.getTriNodes(triindex)
        return a.planeNorm(b,c)
    
    def getSharedNodeTris(self,triindex):
        tris=set()
        for n in self.tris[triindex]:
            tris.update(self.nodeelem[n])
            
        tris.remove(triindex)
        return list(sorted(tris))
    
    def getNearestTri(self,pt):
        def triDist(tri):
            norm=self.getTriNorm(tri)
            pt=self.tricenters[tri]
            return abs(pt.planeDist(pt,norm))
        
        nearestnode=min([n for n in range(self.nodes.n()) if self.nodeelem[n]],key=lambda n:self.nodes[n].distToSq(pt))
        tris=self.nodeelem[nearestnode]
        
        return min(tris,key=triDist)
        
    def getPath(self,starttri,endtri,acceptTri=None):
        return dijkstra(self.adj,starttri,endtri,lambda i,j:self.tridists[(i,j)],acceptTri)
    
    
def loadArchitecture(path,section):
    '''
    Load the architecture from the given file `path' and return values from the given section (endo or epi). The
    return value is a tuple containing:
        landmarks: 0-based indices of landmark nodes in the atlas
        lmlines  : 0-based index pairs defining lines between indices in landmarks
        lmregions: a list of maps, each map defining a region which are mappings from a 0-based indices into lmlines to 
                   the line's landmark index pair
        lmstim   : a per-region specifier list stating which lines (L# for index #) or atlas node (N#) defines stimulation
        lground  : a per-region specifier list stating which lines (L# for index #) or atlas node (N#) defines ground
        appendageregion: region number for the appendage
        appendagenode: node index for the appendage's division node which must be generated
    '''
    c=configparser.SafeConfigParser()
    assert len(c.read(path))>0
    
    landmarks=ast.literal_eval(c.get(section,'landmarks')) # 0-based node indices
    lines=ast.literal_eval(c.get(section,'lines')) # 1-based landmark indices
    regions=ast.literal_eval(c.get(section,'regions')) # 1-based landmark indices
    stimulus=ast.literal_eval(c.get(section,'stimulus')) # per region
    ground=ast.literal_eval(c.get(section,'ground')) # per region
    appendageregion=ast.literal_eval(c.get(section,'appendageregion'))
    appendagenode=ast.literal_eval(c.get(section,'appendagenode'))
    appendagelmindex=ast.literal_eval(c.get(section,'appendagelmindex'))
    
#    types=ast.literal_eval(c.get(section,'type')) # per region    
    
    # indices that don't exist are for landmarks that need to be calculated
    
    lmlines=[subone(l) for l in lines ]#if max(l)<=len(landmarks)] # filter for lines with existing node indices
    lmregions=[subone(r) for r in regions]
#    lmregions=[subone(r) for r in regions if all(i<=len(landmarks) for i in r)]
    
    lmstim=stimulus#[:len(lmregions)]
    lmground=ground#[:len(lmregions)]
    
    allregions=[]
    for r in lmregions:
        lr={i:(a,b) for i,(a,b) in enumerate(lmlines) if a in r and b in r}
        if len(lr)>2:
            allregions.append(lr)

    return landmarks,lmlines,allregions,lmstim,lmground, appendageregion,appendagenode,appendagelmindex
    

def writeMeshFile(filename,nodes,inds,nodegroup,indgroup,dim):
    '''Write a medit format mesh file to `filename'.'''
    with open(filename,'w') as o:
        print('MeshVersionFormatted 1',file=o)
        print('Dimension %i'%dim,file=o)
        print('Vertices',file=o)
        print(len(nodes),file=o)
        
        for n in range(len(nodes)):
            for v in tuple(nodes[n])[:dim]:
                print('%20.10f'%v,end=' ',file=o)
                
            group=0 if nodegroup is None else nodegroup[n]
            
            print(group,file=o)
            
        print('Triangles' if len(inds[0])==3 else 'Tetrahedra',file=o)
        print(len(inds),file=o)
        
        for n in range(len(inds)):
            print(*['%10i'%(t+1) for t in inds[n]],file=o,end=' ')
            group=0 if indgroup is None else indgroup[n]
            
            print(group,file=o)


def createNodeQuery(nodes):
    '''
    Create a cKDTree object from `nodes' and return a query function which accepts a position and radius value. The
    function will return the nearest point index to the given position if radius<=0 and a list of indices of points
    within the given radius of the position otherwise. The node list is also returned as a second return value.
    '''
    tree=cKDTree(np.asarray(list(map(tuple,nodes))))
    
    def _query(pos,radius=0):
        '''Query `nodes' for the nearest node to `pos' if `radius'<=0 or a list of those within `radius' otherwise.'''
        pos=tuple(pos)
        if radius<=0:
            return tree.query(pos)[1],tree.data
        else:
            return tree.query_ball_point(pos,radius),tree.data
        
    return _query


def registerSubjectToTarget(subjectObj,targetObj,targetTrans,outdir,decimpath,VTK):
    '''
    Register the `subjectObj' mesh object to `targetObj' mesh object putting data into directory `outdir'. The subject 
    will be decimated to have roughly the same number of nodes as the target and then stored as subject.vtk in `outdir'. 
    Registration is done with Deformetrica and result stored as 'Registration_subject_to_subject_0__t_9.vtk' in `outdir'.
    If `targetTrans' must be None or a transform which is applied to the nodes of `targetObj' before registration.
    '''
    dpath=os.path.join(outdir,decimatedFile)
    tmpfile=os.path.join(outdir,'tmp.vtk')
    
    # if a transform is given, apply that transform to the target mesh when saving otherwise do no transformation
    vecfunc=(lambda i: tuple(targetTrans*i)) if targetTrans else None

    shutil.copy(os.path.join(deformDir,datasetFile),os.path.join(outdir,datasetFile)) # copy dataset file unchanged
 
    model=open(os.path.join(deformDir,modelFile)).read()
    model=model.replace('%1',str(dataSigma))
    model=model.replace('%2',str(kernelWidthSub))
    model=model.replace('%3',str(kernelType))
    model=model.replace('%4',str(kernelWidthDef))

    with open(os.path.join(outdir,modelFile),'w') as o: # save modified model file
        o.write(model)
        
    optim=open(os.path.join(deformDir,optimFile)).read()
    optim=optim.replace('%1',str(stepSize))
    
    with open(os.path.join(outdir,optimFile),'w') as o: # save modified optimization file
        o.write(optim)

    VTK.saveLegacyFile(tmpfile,subjectObj,datasettype='POLYDATA')
    VTK.saveLegacyFile(os.path.join(outdir,targetFile),targetObj,datasettype='POLYDATA',vecfunc=vecfunc)
        
    snodes=subjectObj.datasets[0].getNodes()
    tnodes=targetObj.datasets[0].getNodes()
    
    sizeratio=float(tnodes.n())/snodes.n()
    sizepercent=str(100*(1-sizeratio))[:6] # percent to decimate by
    
    # decimate the mesh most of the way towards having the same number of nodes as the target
    ret,output=eidolon.execBatchProgram(decimpath,tmpfile,dpath,'-reduceby',sizepercent,'-ascii',logcmd=True)
    assert ret==0,output
    
    env={'LD_LIBRARY_PATH':deformDir}
    ret,output=eidolon.execBatchProgram(deformExe,"registration", "3D", modelFile, datasetFile, optimFile, "--output-dir=.",cwd=outdir,env=env,logcmd=True)
    assert ret==0,output
    
    return output


def transferLandmarks(archFilename,fieldname,targetObj,targetTrans,subjectObj,outdir,VTK):
    '''
    Register the landmarks defined as node indices on `targetObj' to equivalent node indices on `subjectObj' via the
    decimated and registered intermediary stored in `outdir'. The result is a list of index pairs associating a node
    index in `subjectObj' for every landmark index in `targetObj'.
    '''
    decimated=os.path.join(outdir,decimatedFile)
    registered=os.path.join(outdir,registeredFile)
    
    arch=loadArchitecture(archFilename,fieldname)
    lmarks,lines=arch[:2]
    appendagelmindex=arch[-1]
    
    # append the index for the estimated appendage node, this will have to be adjusted manually after registration
    lmarks.append(appendagelmindex) 
    
    reg=VTK.loadObject(registered) # mesh registered to target
    dec=VTK.loadObject(decimated) # decimated unregistered mesh
    
    tnodes=targetObj.datasets[0].getNodes() # target points
    rnodes=reg.datasets[0].getNodes() # registered decimated points
    dnodes=dec.datasets[0].getNodes() # unregistered decimated points
    snodes=subjectObj.datasets[0].getNodes() # original subject points
    
    targetTrans=targetTrans or eidolon.transform()
    lmpoints=[(targetTrans*tnodes[m],m) for m in lmarks] # (transformed landmark node, index) pairs
    
    # find the points in the registered mesh closest to the landmark points in the target object
    query=createNodeQuery(rnodes)
    rpoints=[(query(pt)[0],m) for pt,m in lmpoints] 
    
    # find the subject nodes closes to landmark points in the decimated mesh (which are at the same indices as in the registered mesh)
    query=createNodeQuery(snodes)
    spoints=[(query(dnodes[i])[0],m) for i,m in rpoints]
        
    assert len(spoints)==len(lmpoints)
    assert all(p[0] is not None for p in spoints)
    
    slines=[l for l in lines if max(l)<len(spoints)]
    
    return spoints,slines # return list (i,m) pairs where node index i in the subject mesh is landmark m


def generateTriAdj(tris):
    '''
    Generates a table (n,3) giving the indices of adjacent triangles for each triangle, with a value of `n' indicating a 
    free edge. The indices in each row are in sorted order rather than per triangle edge. The result is the dual of the 
    triangle mesh represented as the (n,3) array and a map relating the mesh's ragged edges to their triangle.
    '''
    edgemap = {} # maps edges to the first triangle having that edge
    result=IndexMatrix(tris.getName()+'Adj',tris.n(),3)
    result.fill(tris.n())
    
    # Find adjacent triangles by constructing a map from edges defined by points (a,b) to the triangle having that edge,
    # when that edge is encountered twice then the current triangle is adjacent to the one that originally added the edge.
    for t1,tri in enumerate(tris): # iterate over each triangle t1
        for a,b in successive(tri,2,True): # iterate over each edge (a,b) of t1
            k=(min(a,b),max(a,b)) # key has uniform edge order
            t2=edgemap.pop(k,None) # attempt to find edge k in the map, None indicates edge not found
            
            if t2 is not None: # an edge is shared if already encountered, thus t1 is adjacent to t2
                result[t1]=sorted(set(result[t1]+(t2,)))
                result[t2]=sorted(set(result[t2]+(t1,)))
            else:
                edgemap[k]=t1 # first time edge is encountered, associate this triangle with it

    return result,edgemap


@timing
def getAdjTo(adj,start,end):
    '''
    Returns a subgraph of `adj',represented as a node->[neighbours] dict, which includes nodes `start' and `end'. 
    If `end' is None or an index not appearing in the mesh, the result will be the submesh contiguous with `start'.
    '''
    visiting=set([start])
    found={}
    numnodes=adj.n()
    
    while visiting and end not in found:
        visit=visiting.pop()
        neighbours=[n for n in adj[visit] if n<numnodes]
        found[visit]=neighbours
        visiting.update(n for n in neighbours if n not in found)
                
    return found
        
   
def generateNodeElemMap(numnodes,tris):
    '''Returns a list relating each node index to the set of element indices using that node.'''
    nodemap=[set() for _ in range(numnodes)]

    for i,tri in enumerate(tris):
        for n in tri:
            nodemap[n].add(i)
            
    #assert all(len(s)>0 for s in nodemap), 'Unused nodes in triangle topology'
            
    return nodemap
    

def generateSimplexEdgeMap(numnodes,simplices):
    '''
    Returns a list relating each node index to the set of node indices joined to it by graph edges. This assumes the mesh
    has `numnodes' number of nodes and simplex topology `simplices'.
    '''
    nodemap=[set() for _ in range(numnodes)]

    for simplex in simplices:
        simplex=set(simplex)
        for s in simplex:
            nodemap[s].update(simplex.difference((s,)))
            
    return nodemap
    

@timing
def dijkstra(adj, start, end,distFunc,acceptTri=None):
    #http://benalexkeen.com/implementing-djikstras-shortest-path-algorithm-with-python/
    # shortest paths is a dict of nodes to previous node and distance
    paths = {start: (None,0)}
    curnode = start
    visited = set()
    # consider only subgraph containing start and end, this expands geometrically so should contain the minimal path
    adj=getAdjTo(adj,start,end) 
    
    eidolon.printFlush(len(adj))
    
    if acceptTri is not None:
        accept=lambda a: (a in adj and acceptTri(a))
    else:
        accept=lambda a: a in adj
    
    while curnode != end:
        visited.add(curnode)
        destinations = list(filter(accept,adj[curnode]))
        curweight = paths[curnode][1]

        for dest in destinations:
            weight = curweight+distFunc(curnode,dest)
            
            if dest not in paths or weight < paths[dest][1]:
                paths[dest] = (curnode, weight)
        
        nextnodes = {node: paths[node] for node in paths if node not in visited}
        
        if not nextnodes:
            raise ValueError("Route %i -> %i not possible"%(start,end))
            
        # next node is the destination with the lowest weight
        curnode = min(nextnodes, key=lambda k:nextnodes[k][1])
    
    # collect path from end node back to the start
    path = []
    while curnode is not None:
        path.insert(0,curnode)
        curnode = paths[curnode][0]
        
    return path
    
    
def subone(v):
    return tuple(i-1 for i in v)
    

def findNearestIndex(pt,nodelist):
    return min(range(len(nodelist)),key=lambda i:pt.distToSq(nodelist[i]))
        

def findFarthestIndex(pt,nodelist):
    return max(range(len(nodelist)),key=lambda i:pt.distToSq(nodelist[i]))


def getContiguousTris(graph,starttri,acceptTri):
    accepted=[starttri]
    
    adjacent=first(i for i in graph.getSharedNodeTris(starttri) if i not in accepted and acceptTri(i))
    
    while adjacent is not None:
        accepted.append(adjacent)
        
        for a in accepted[::-1]:
            allneighbours=graph.getSharedNodeTris(a)
            adjacent=first(i for i in allneighbours if i not in accepted and acceptTri(i))
            if adjacent:
                break
        
    return accepted


@timing
def findTrisBetweenNodes(start,end,landmarks,graph):
    eidolon.printFlush('Nodes:',start,end)
    start=landmarks[start]
    end=landmarks[end]

    assert 0<=start<len(graph.nodeelem)
    assert 0<=end<len(graph.nodeelem)

    starttri=first(graph.nodeelem[start])
    endtri=first(graph.nodeelem[end])
    
    assert starttri is not None
    assert endtri is not None
    
    nodes=graph.nodes
    startnode=nodes[start]
    endnode=nodes[end]
    
    easypath= graph.getPath(starttri,endtri)
    
    midnode=graph.tricenters[easypath[len(easypath)//2]]
    
    # define planes to bound the areas to search for triangles to within the space of the line
    splane=plane(startnode,midnode-startnode)
    eplane=plane(endnode,midnode-endnode)


    # adjust the plane's positions to account for numeric error
    adjustdist=1e1
    splane.moveUp(-adjustdist)
    eplane.moveUp(-adjustdist)
    
    
    assert starttri is not None
    assert endtri is not None
    
    # TODO: plane normal determination still needs work
    #linenorm=midnode.planeNorm(startnode,endnode)
    #linenorm=graph.getTriNorm(easypath[len(easypath)//2]).cross(midnode-startnode)
    linenorm=eidolon.avg(graph.getTriNorm(e) for e in easypath).cross(midnode-startnode)
    
    lineplane=plane(splane.center,linenorm)
    
    indices=set([starttri,endtri]) # list of element indices on lineplane between splane and eplane
    
    for i in range(graph.tris.n()):
        trinodes=graph.getTriNodes(i)
        
        numabove=lineplane.numPointsAbove(trinodes)
        if numabove in (1,2) and splane.between(trinodes,eplane):    
            indices.add(i)
    
    accepted=getContiguousTris(graph,starttri,lambda i:i in indices)
    
    if endtri not in accepted or len(easypath)<len(accepted):
        eidolon.printFlush('---Resorting to easypath')
        accepted=easypath 
        
    return accepted
    
    
@timing
def assignRegion(region,index,assignmat,landmarks,linemap,graph):
    
    def getEnclosedGraph(adj,excludes,start):
        visiting=set([start])
        found=set()
        numnodes=adj.n()
        
        assert start is not None
        
        while visiting:
            visit=visiting.pop()
            neighbours=[n for n in adj.getRow(visit) if n<numnodes and n not in excludes]
            found.add(visit)
            visiting.update(n for n in neighbours if n not in found)
                    
        return found
    
    # collect all tri indices on the border of this region
    bordertris=set()
    for lineindex,(a,b) in region.items():
        if (a,b) in linemap:
            line=linemap[(a,b)]
        else:
            line=findTrisBetweenNodes(a,b,landmarks,graph)
            linemap[(a,b)]=line
            linemap[(b,a)]=line
            
            # assign line ID to triangles on the line
            for tri in line:
                assignmat[tri,0]=lineindex
            
        bordertris.update(line)
        
    bordertri=graph.tricenters[first(bordertris)]
    
    farthest=max(range(len(graph.tris)),key=lambda i:graph.tricenters[i].distToSq(bordertri))
    
    maxgraph=getEnclosedGraph(graph.adj,bordertris,farthest)
    
    for tri in range(len(graph.tris)):
        if tri in bordertris or tri not in maxgraph:
            if assignmat[tri,1]<0:
                assignmat[tri,1]=index
            elif assignmat[tri,2]<0:
                assignmat[tri,2]=index
            elif assignmat[tri,3]<0:
                assignmat[tri,3]=index
    
    
@timing
def generateRegionField(obj,landmarkObj,regions,appendageregion,appendagenode,task=None):
    ds=obj.datasets[0]
    nodes=ds.getNodes()
    tris=first(ind for ind in ds.enumIndexSets() if ind.m()==3 and bool(ind.meta(StdProps._isspatial)))
    lmnodes=landmarkObj.datasets[0].getNodes()
    linemap={}
    
    landmarks={i:nodes.indexOf(lm)[0] for i,lm in enumerate(lmnodes)}
    
    assert all(0<=l<nodes.n() for l in landmarks)
    
    graph=TriMeshGraph(nodes,tris)
    
    edgenodeinds=set(eidolon.listSum(graph.ragged)) # list of all node indices on the ragged edge
    
    filledregions=RealMatrix(fieldNames._regions,tris.n(),4)
    filledregions.meta(StdProps._elemdata,'True')
    filledregions.fill(-10)
    
    
    #landmarks[appendagenode]=0 # TODO: skipping appendage node for now
    
    for region in regions:
        for a,b in region.values():
            if appendagenode not in (a,b):
                if a in landmarks and b not in landmarks:
                    oldlmnode=nodes[landmarks[a]]
                    newlm=b
                elif b in landmarks and a not in landmarks:
                    oldlmnode=nodes[landmarks[b]]
                    newlm=a
                else:
                    continue
                
                newlmnode=min(edgenodeinds,key=lambda i:nodes[i].distToSq(oldlmnode)) # ragged edge node closest to landmark
                landmarks[newlm]=newlmnode
                
#                eidolon.printFlush(newlm,newlmnode,graph.getPath(min(a,b),newlmnode),'\n')
                
#                line=findTrisBetweenNodes(a,b,landmarks,graph)
#                for tri in line:
#                    filledregions[tri,0]=max(a,b)
                
    
    if task:
        task.setMaxProgress(len(regions))

    for rindex,region in enumerate(regions):
        eidolon.printFlush('Region',rindex,'of',len(regions),region)
        
        allnodes=set(eidolon.listSum(region.values()))
        if all(a in landmarks for a in allnodes):
            assignRegion(region,rindex,filledregions,landmarks,linemap,graph)   
        else:
            eidolon.printFlush('Skipping',rindex,[a for a in allnodes if a not in landmarks])
            
        if task:
            task.setProgress(rindex+1)
            
    return filledregions,linemap


def extractTriRegion(nodes,tris,acceptFunc):
    '''
    Extract the region from the mesh (nodes,tris) as defined by the triangle acceptance function `acceptFunc'. The return
    value is a tuple containing the list of new nodes, a list of new tris, a map from old node indices in `nodes' to new 
    indices in the returned node list, and a map from triangle indices in `tris' to new ones in the returned triangle list.
    '''
    #old -> new
    newnodes=[] # new node set
    newtris=[] # new triangle set
    nodemap={} # maps old node indices to new
    trimap={} # maps old triangle indices to new
    
    for tri in range(len(tris)):
        if acceptFunc(tri):
            newtri=list(tris[tri])
            
            for i,n in enumerate(newtri):
                if n not in nodemap:
                    nodemap[n]=len(newnodes)
                    newnodes.append(nodes[n])
                    
                newtri[i]=nodemap[n]
                
            trimap[tri]=len(newtris)
            newtris.append(newtri)
            
    return newnodes,newtris,nodemap,trimap
    

def calculateMeshGradient(prefix,nodes,elems,groups,VTK):
    '''Calculate the laplace gradient for the mesh given as (nodes,elems,groups) using sfepy.'''
    tempdir=os.path.dirname(prefix)
    infile=prefix+'.mesh'
    logfile=prefix+'.log'
    outfile=prefix+'.vtk'
    probfile=prefix+'.py'
    
    writeMeshFile(infile,nodes,elems,groups,None,3)
    
    with open(problemFile) as p:
        with open(probfile,'w') as o:
            o.write(p.read()%{'inputfile':infile,'outdir':tempdir})
        
    p=
ProblemConf.from_file(probfile)