# -*- 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 = [] 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 ''' Writes the header given and open file descripion, number of verticies and number of edges. ''' def writeHeader(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(open_file, vertex): open_file.write(struct.pack(' 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(nodeList, edgeList, R=(100, 100, 100)): #get a list of all node positions in the network P = [] for e in edgeList: 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 = aabb(nodeList, edgeList) #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