python/network.py

# -*- coding: utf-8 -*-
"""
Created on Sat Sep 16 16:34:49 2017
@author: pavel
"""
import struct
import numpy as np
import scipy as sp
import networkx as nx
import matplotlib.pyplot as plt
import math
'''
    Definition of the Node class
    Duplicate of the node class in network
    Stores the physical position, outgoing edges list and incoming edges list.
'''
class Node:
    def __init__(self, point, outgoing, incoming):
        self.p = point
        self.o = outgoing
        self.i = incoming
#    def p():
#        return self.p
'''
    Definition of the Fiber class.
    Duplicate of the Node class in network
    Stores the starting vertex, the ending vertex, the points array and the radius array
'''
class Fiber:
    
    def __init__ (self):
        self.v0 = 0
        self.v1 = 0
        self.points = []
        self.radii = []
        
        #NOTE: there is no function overloading in Python
#    def __init__ (self, p1, p2, pois, rads):
#        self.v0 = p1
#        self.v1 = p2
#        self.points = pois
#        self.radii = rads
    
    '''
        return the length of the fiber.
    '''        
    def length(self):
        length = 0
        for i in range(len(self.points)-1):
            length = length + math.sqrt(pow(self.points[i][0]- self.points[i+1][0],2) + pow(self.points[i][1]- self.points[i+1][1],2) + pow(self.points[i][2]- self.points[i+1][2],2))
        return length
        
    '''
        returns the turtuosity of the fiber.
    '''    
    def turtuosity(self):
        turtuosity = 0
        distance = math.sqrt(math.pow(self.points[0][0]- self.points[len(self.points)-1][0],2) + math.pow(self.points[0][1]- self.points[len(self.points)-1][1],2) + math.pow(self.points[0][2]- self.points[len(self.points)-1][2],2))
        turtuosity = self.length()/distance
        #print(turtuosity)
        
        return turtuosity
        
    '''
        returns the volume of the fiber.
    '''
    def volume(self):
        volume = 0
        for i in range(len(self.points)-1):
            volume = volume + 1.0/3.0 * math.pi * (math.pow(self.radii[i],2) + math.pow(self.radii[i+1],2) + self.radii[i]*self.radii[i+1]) * math.sqrt(math.pow(self.points[i][0]- self.points[i+1][0],2) + math.pow(self.points[i][1]- self.points[i+1][1],2) + math.pow(self.points[i][2]- self.points[i+1][2],2))
        #print(volume)
        return volume
    
class NWT:   
    
    '''
        Writes the header given and open file descripion, number of verticies and number of edges.
    '''
    def writeHeader(self, open_file, numVerts, numEdges):
        txt = "nwtFileFormat fileid(14B), desc(58B), #vertices(4B), #edges(4B): bindata"
        b = bytearray()
        b.extend(txt.encode())
        open_file.write(b)
        open_file.write(struct.pack('i', numVerts))
        open_file.write(struct.pack('i', numEdges))
        
    
    '''
        Writes a single vertex to a file.
    '''
    def writeVertex(self, open_file, vertex):
        open_file.write(struct.pack('<f',vertex.p[0]))
        open_file.write(struct.pack('<f',vertex.p[1]))
        open_file.write(struct.pack('<f',vertex.p[2]))
        open_file.write(struct.pack('i', len(vertex.o)))
        open_file.write(struct.pack('i', len(vertex.i)))
        for j in range(len(vertex.o)):
            open_file.write(struct.pack('i',vertex.o[j]))
            
        for j in range(len(vertex.i)):
            open_file.write(struct.pack('i', vertex.i[j]))    
            
        return
    
    '''
        Writes a single fiber to a file.
    '''
    def writeFiber(self, open_file, edge):
        open_file.write(struct.pack('i',edge.v0))
        open_file.write(struct.pack('i',edge.v1))
        open_file.write(struct.pack('i',len(edge.points)))
        for j in range(len(edge.points)):
            open_file.write(struct.pack('<f', edge.points[j][0]))
            open_file.write(struct.pack('<f', edge.points[j][1]))
            open_file.write(struct.pack('<f', edge.points[j][2]))
            open_file.write(struct.pack('<f', edge.radii[j]))
            
        return
    
    '''
        Writes the entire network to a file in str given the vertices array and the edges array.
    '''
    def exportNWT(self, str, vertices, edges):
        with open(str, "wb") as file:
            self.writeHeader(file, len(vertices), len(edges))
            for i in range(len(vertices)):
                self.writeVertex(file, vertices[i])
                
            for i in range(len(edges)):
                self.writeFiber(file, edges[i])
                
        return
    
    
    '''
        Reads a single vertex from an open file and returns a node Object.
    '''
    def readVertex(self, open_file):
        points = np.tile(0., 3)
        bytes = open_file.read(4)
        points[0] = struct.unpack('f', bytes)[0]
        bytes = open_file.read(4)
        points[1] = struct.unpack('f', bytes)[0]
        bytes = open_file.read(4)
        points[2] = struct.unpack('f', bytes)[0]
        bytes = open_file.read(4)
        
        numO = int.from_bytes(bytes, byteorder='little')
        outgoing = np.tile(0, numO)
        bts = open_file.read(4)
        numI = int.from_bytes(bts, byteorder='little')
        incoming = np.tile(0, numI)
        for j in range(numO):
            bytes = open_file.read(4)
            outgoing[j] = int.from_bytes(bytes, byteorder='little')
            
        for j in range(numI):
            bytes = open_file.read(4)
            incoming[j] = int.from_bytes(bytes, byteorder='little')
            
        node = Node(points, outgoing, incoming)    
        return node
        
        
    '''
        Reads a single fiber from an open file and returns a Fiber object .   
    '''
    def readFiber(self, open_file):
        bytes = open_file.read(4)
        vtx0 = int.from_bytes(bytes, byteorder = 'little')
        bytes = open_file.read(4)
        vtx1 = int.from_bytes(bytes, byteorder = 'little')
        bytes = open_file.read(4)
        numVerts = int.from_bytes(bytes, byteorder = 'little')
        pts = []
        rads = []
        
        for j in range(numVerts):
            point = np.tile(0., 3)
            bytes = open_file.read(4)
            point[0] = struct.unpack('f', bytes)[0]
            bytes = open_file.read(4)
            point[1] = struct.unpack('f', bytes)[0]
            bytes = open_file.read(4)
            point[2] = struct.unpack('f', bytes)[0]
            bytes = open_file.read(4)
            radius = struct.unpack('f', bytes)[0]
            pts.append(point)
            rads.append(radius)
            
        F = Fiber(vtx0, vtx1, pts, rads)
            
        return F
    
    '''
        Imports a NWT file at location str.
        Returns a list of Nodes objects and a list of Fiber objects.
    '''
class Network:
    
    def __init__(self, str):
        with open(str, "rb") as file:
            header = file.read(72)
            bytes = file.read(4)
            numVertex = int.from_bytes(bytes, byteorder='little')
            bytes = file.read(4)
            numEdges = int.from_bytes(bytes, byteorder='little')
            
            self.N = []
            self.F = []
            for i in range(numVertex):
                node = NWT.readVertex(file)
                self.N.append(node)
    
            for i in range(numEdges):
                edge = NWT.readFiber(file)
                self.F.append(edge)
    
    '''
    Creates a graph from a list of nodes and a list of edges.
    Uses edge length as weight.
    Returns a NetworkX Object.
    '''
#    def createLengthGraph(self):
#        G = nx.Graph()
#        for i in range(len(self.nodeList)):
#            G.add_node(i, p=V[i].p)
#        for i in range(len(self.edgeList)):
#            G.add_edge(self.edgeList[i].v0, self.edgeList[i].v1, weight = E[i].length())
#            
#        return G
#    '''
#    Creates a graph from a list of nodes and a list of edges.
#    Uses edge turtuosity as weight.
#    Returns a NetworkX Object.
#    '''    
#    def createTortuosityGraph(nodeList, edgeList):
#        G = nx.Graph()
#        for i in range(len(nodeList)):
#            G.add_node(i, p=V[i].p)
#        for i in range(len(edgeList)):
#            G.add_edge(edgeList[i].v0, edgeList[i].v1, weight = E[i].turtuosity())
#            
#        return G
    
#    '''
#    Creates a graph from a list of nodes and a list of edges.
#    Uses edge volume as weight.
#    Returns a NetworkX Object.
#    '''    
#    def createVolumeGraph(nodeList, edgeList):
#        G = nx.Graph()
#        for i in range(len(nodeList)):
#            G.add_node(i, p=V[i].p)
#        for i in range(len(edgeList)):
#            G.add_edge(edgeList[i].v0, edgeList[i].v1, weight = E[i].volume())
#            
#        return G
#'''
#Returns the positions dictionary for the Circular layout.
#'''    
#def getCircularLayout(graph, dim, radius):
#    return nx.circular_layout(graph, dim, radius)
#
#'''
#Return the positions dictionary for the Spring layout.
#'''    
#def getSpringLayout(graph, pos, iterations, scale):
#    return nx.spring_layout(graph, 2, None, pos, iterations, 'weight', scale, None)
#        
#'''
#Draws the graph.
#'''        
#def drawGraph(graph, pos):
#    nx.draw(graph, pos)
#    return
    def aabb(self):
    
        lower = self.N[0].p.copy()
        upper = lower.copy()
        for i in self.nodeList:
            for c in range(len(lower)):
                if lower[c] > i.p[c]:
                    lower[c] = i.p[c]
                if upper[c] < i.p[c]:
                    upper[c] = i.p[c]
        return lower, upper
    
    #calculate the distance field at a given resolution
    #   R (triple) resolution along each dimension
    def distancefield(self, R=(100, 100, 100)):
        
        #get a list of all node positions in the network
        P = []
        for e in self.F:
            for p in e.points:
                P.append(p)
                
        #turn that list into a Numpy array so that we can create a KD tree
        P = np.array(P)
        
        #generate a KD-Tree out of the network point array
        tree = sp.spatial.cKDTree(P)
        
        plt.scatter(P[:, 0], P[:, 1])
        
        #specify the resolution of the ouput grid
        R = (200, 200, 200)
        
        #generate a meshgrid of the appropriate size and resolution to surround the network
        lower, upper = self.aabb(self.N, self.F)           #get the space occupied by the network
        x = np.linspace(lower[0], upper[0], R[0])   #get the grid points for uniform sampling of this space
        y = np.linspace(lower[1], upper[1], R[1])
        z = np.linspace(lower[2], upper[2], R[2])
        X, Y, Z = np.meshgrid(x, y, z)
        #Z = 150 * numpy.ones(X.shape)
        
        
        Q = np.stack((X, Y, Z), 3)
        
        
        D, I = tree.query(Q)
        
        return D
    
    #returns the number of points in the network
    def npoints(self):
        n = 0
        for f in self.F:
            n = n + len(f.points) - 2
        n = n + len(self.N)
        return n
    
    #returns the number of linear segments in the network
    def nsegments(self):
        n = 0
        for f in self.F:
            n = n + len(f.points) - 1
        return n
    
    def vectorize(self):
        #initialize three coordinate vectors ()
        n = self.nsegments(self.N, self.F)
        X = np.zeros((n))
        Y = np.zeros((n))
        Z = np.zeros((n))
        idx = 0
        for i in range(0, len(self.F)):
            for j in range(0, len(self.F[i].points)-1):
                X[idx] = self.F[i].points[j][0]-self.F[i].points[j+1][0]
                Y[idx] = self.F[i].points[j][1]-self.F[i].points[j+1][1]
                Z[idx] = self.F[i].points[j][2]-self.F[i].points[j+1][2]
                idx = idx + 1
        return X, Y, Z