Merge branch 'master' of git.stim.ee.uh.edu:codebase/stimlib

sberisha
2 parents 35534291 176a9b1c
Showing 22 changed files with 2712 additions and 2795 deletions Show diff stats
stim/biomodels/nwt_format.pptx → docs/nwt_format.pptx
docs/test.nwt
python/fibernet.py
python/network.py → python/network_dep.py
python/nwt.py
stim/biomodels/centerline.h
stim/biomodels/network.h
stim/biomodels/test.nwt
stim/cuda/branch_detection.cuh
stim/cuda/filter.cuh
stim/gl/gl_spider.h
stim/iVote/ivote2.cuh
stim/iVote/ivote2/local_max.cuh
stim/image/image.h
stim/math/filters/median2.cuh
stim/math/matrix.h
stim/math/matrix_sq.h
stim/math/spharmonics.h
stim/optics/scalarmie.h
stim/visualization/gl_network.h
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 19 14:45:39 2018
+
+@author: david
+"""
+
+import nwt
+import matplotlib.pyplot as plt
+import numpy
+
+
+class fibernet:
+    
+    def __init__(self, filename):
+        self.nwt(filename)
+        self.buildgraph()
+    
+    #load a network from a NWT file
+    def nwt(self, filename):
+        net = nwt.NWT(filename)
+        
+        #first insert the points where fibers are connected
+        self.P = []
+        for vi in range(0, len(net.v)):                                         #store each vertex in the point list
+            self.P.append( (net.v[vi].p[0], net.v[vi].p[1], net.v[vi].p[2]) )
+        
+        #now insert each fiber, connecting them appropriately using previously inserted points (vertices)
+        self.F = []
+        for ei in range(0, len(net.e)):
+            f = [numpy.int(net.e[ei].v[0])]                                     #create a new fiber to store point indices, initialize with the first vertex ID
+            for pi in range(1, net.e[ei].p.shape[0]):
+                self.P.append( (net.e[ei].p[pi,0], net.e[ei].p[pi,1], net.e[ei].p[pi,2]))    #append the current point to the global point list
+                f.append(len(self.P) - 1)                                       #append the current point ID
+            f.append(numpy.int(net.e[ei].v[1]))                                 #append the last vertex ID
+            self.F.append(f)
+    
+    #build a graph from the available geometry
+    def buildgraph(self):
+        p_to_v = numpy.ones((self.npoints()), numpy.int) * (-1)                            #create an array that maps a point ID to a vertex ID
+        
+        self.V = []                                                                  #create an empty vertex list
+        self.E = []                                                                  #create an empty edge list
+        for fi in range(0, self.nfibers()):                                     #for each fiber in the network
+            if p_to_v[self.F[fi][0]] == -1:                                          #if the first point hasn't been encountered
+                p_to_v[self.F[fi][0]] = len(self.V)                                       #add it to the mapping
+                self.V.append( (self.F[fi][0], []) )                                              #add a new vertex representing this point
+            if p_to_v[self.F[fi][-1]] == -1:
+                p_to_v[self.F[fi][-1]] = len(self.V)
+                self.V.append( (self.F[fi][-1], []) )
+                
+            v0 = p_to_v[self.F[fi][0]]                                          #get the two vertex indices corresponding to this edge
+            v1 = p_to_v[self.F[fi][-1]]
+            
+            self.E.append( (v0, v1) )                                                #create an edge representing this fiber and add it to the edge list
+            self.V[v0][1].append(fi)                                                 #add this edge ID to the corresponding vertices to update the adjacency list
+            self.V[v1][1].append(fi)
+            
+    #return the length of a fiber specified by fiber ID fid
+    def length(self, fid):
+        p = self.points(fid)                                                 #get the points corresponding to the fiber
+        l = 0                                                                   #initialize the length to zero
+        for pi in range(1, len(p)):                                             #for each line segment
+            v = p[pi] - p[pi-1]                                                 #calculate the vector representing the line segment
+            l = l + numpy.linalg.norm(numpy.array(v))                                  #add the length of the line segment
+        return l                                                                #return the total length   
+
+    #return a list of vertices adjacent to the vertex specified by vid
+    def adjacent(self, vid):
+        adj = []                                                                #initialize an empty adjacency list
+        for e in self.V[vid][1]:                                                #for each edge connected to vid
+            if self.E[e][0] != vid:
+                adj.append(self.E[e][0])
+            else:
+                adj.append(self.E[e][1])                
+        return adj
+    
+    #return the fiber network graph as an adjacency list
+    def adjacencylist(self):
+        A = []
+        for vi in range(0, len(self.V)):                                                        #for each vertex in the graph
+            A.append( self.adjacent(vi) )                                       #add the adjacency list for the vertex into the main adjacency list
+        return A
+
+    #return the fiber network graph as an adjacency matrix
+    # possible values for the matrix include:
+    #           binary - 1 when there is a connection, 0 otherwise
+    #           edgeID - stores the edge ID corresponding to a connection, and -1 otherwise
+    #           length - stores the length of the connecting fiber, and -1 if there is no connection
+    def adjacencymatrix(self, value="binary"):
+        if value == "binary":
+            M = numpy.ones( (len(self.V), len(self.V)), numpy.bool )
+        elif value == "edgeID":
+            M = numpy.ones( (len(self.V), len(self.V)), numpy.int ) * -1
+        elif value == "length":
+            M = numpy.ones( (len(self.V), len(self.V)) ) * -1
+        else:
+            return []
+            
+        for vi in range(0, len(self.V)):                                        #for each vertex in the graph
+            A = self.adjacent(vi)                                               #get the adjacency matrix
+            for a in A:
+                if value == "binary":
+                    M[vi,a] = M[a,vi] = 1
+                if value == "edgeID":
+                    M[vi,a] = M[a,vi] = self.findedge(a, vi)
+                if value == "length":
+                    M[vi,a] = M[a,vi] = self.length(self.findedge(a,vi))
+        return M
+    
+    #find the edge ID corresponding to two vertices
+    def findedge(self, v0, v1):
+        
+        edges = self.V[v0][1]                                                   #get the list of candidate edges connected to v0
+        for e in edges:                                                         #for each edge
+            if self.E[e][0] == v1 or self.E[e][1] == v1:
+                return e
+        
+        return -1                                                               #return -1 if the edge doesn't exist
+    
+    #get the points associated with a specified fiber ID fid    
+    def points(self, fid):                                                   #get the points corresponding to a specified fiber index
+        f = self.F[fid]                                                         #get the list of point indices
+        return numpy.array(self.P)[f]
+    
+    #get the total number of points in the network
+    def npoints(self):
+        return len(self.P)
+    
+    #get the total number of fibers in the network
+    def nfibers(self):
+        return len(self.F)
+    
+    #plot the network as a 3D line plot
+    def plot(self):
+        fig = plt.figure()
+        ax = fig.add_subplot(111, projection='3d')
+        
+        for fi in range(0, self.nfibers()):
+            test = fibers.points(fi)
+            ax.plot(test[:, 0], test[:, 1], test[:, 2])
+        return fig
 \ No newline at end of file
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 19 12:38:30 2018
+
+@author: david
+"""
+
+import numpy
+
+class Vertex:
+    def __init__(self, f):                                                      #load a vertex given a file ID at that vertex location
+        self.p = numpy.fromfile(f, dtype=numpy.float32, count=3)                                     #load the vertex coordinates
+        
+        EO = int.from_bytes(f.read(4), "little")                                #get the number of edges moving out of and into this vertex
+        EI = int.from_bytes(f.read(4), "little")
+        
+        self.Eout = numpy.fromfile(f, dtype=numpy.uint32, count=EO)             #get the incoming and outgoing edge lists
+        self.Ein = numpy.fromfile(f, dtype=numpy.uint32, count=EI)
+        
+    def __str__(self):
+        sp = "(" + str(self.p[0]) + ", " + str(self.p[1]) + ", " + str(self.p[2]) + ")"
+        
+        seo = "["
+        for e in self.Eout:
+            seo = seo + " " + str(e)
+        seo = seo + " ]"
+            
+        sei = "["
+        for e in self.Ein:
+            sei = sei + " " + str(e)
+        sei = sei + " ]"
+        
+        return sp + "  out = " + seo + "  in = " + sei
+        
+class Edge:
+    def __init__(self, f):
+        self.v = numpy.fromfile(f, dtype=numpy.uint32, count=2)
+        P = int.from_bytes(f.read(4), "little")
+        p1d = numpy.fromfile(f, dtype=numpy.float32, count=4*P)
+        self.p = numpy.reshape(p1d, (P, 4))
+        
+    def __str__(self):
+        sv = "[ " + str(self.v[0]) + "---> " + str(self.v[1]) + " ]"
+        return sv + str(self.p)
+
+class NWT:
+    def __init__(self, fname):
+        f = open(fname, "rb")                                                        #open the NWT file for binary reading
+        self.filetype = f.read(14)                                                           #this should be "nwtfileformat "
+        self.desc = f.read(58)                                                               #NWT file description
+        V = int.from_bytes(f.read(4), "little")                                         #retrieve the number of vertices (graph)
+        E = int.from_bytes(f.read(4), "little")
+              
+        self.v = []
+        for i in range(0, V):
+            self.v.append(Vertex(f))
+            
+        self.e = []
+        for i in range(1, E):
+            self.e.append(Edge(f))
 \ No newline at end of file
-/*
-Copyright <2017> <David Mayerich>
-
-Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
-
-The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-*/
-
-#ifndef STIM_CENTERLINE_H
-#define STIM_CENTERLINE_H
-
-#include <vector>
-#include <stim/math/vec3.h>
-#include <stim/structures/kdtree.cuh>
-
-namespace stim{
-
-/**	This class stores information about a single fiber represented as a set of geometric points
- *	between two branch or end points. This class is used as a fundamental component of the stim::network
- *	class to describe an interconnected (often biological) network.
- */
-template<typename T>
-class centerline : public std::vector< stim::vec3<T> >{
-
-protected:
-
-	std::vector<T> L;										//stores the integrated length along the fiber (used for parameterization)
-
-	///Return the normalized direction vector at point i (average of the incoming and outgoing directions)
-	vec3<T> d(size_t i) {
-		if (size() <= 1) return vec3<T>(0, 0, 0);						//if there is insufficient information to calculate the direction, return a null vector
-		if (size() == 2) return (at(1) - at(0)).norm();					//if there are only two points, the direction vector at both is the direction of the line segment
-		if (i == 0) return (at(1) - at(0)).norm();						//the first direction vector is oriented towards the first line segment
-		if (i == size() - 1) return (at(size() - 1) - at(size() - 2)).norm();	//the last direction vector is oriented towards the last line segment
-
-		//all other direction vectors are the average direction of the two joined line segments
-		vec3<T> a = at(i) - at(i - 1);
-		vec3<T> b = at(i + 1) - at(i);
-		vec3<T> ab = a.norm() + b.norm();
-		return ab.norm();
-	}
-
-	//initializes the integrated length vector to make parameterization easier, starting with index idx (all previous indices are assumed to be correct)
-	void update_L(size_t start = 0) {
-		L.resize(size());									//allocate space for the L array
-		if (start == 0) {
-			L[0] = 0;											//initialize the length value for the first point to zero (0)
-			start++;
-		}
-
-		stim::vec3<T> d;
-		for (size_t i = start; i < size(); i++) {		//for each line segment in the centerline
-			d = at(i) - at(i - 1);
-			L[i] = L[i - 1] + d.len();				//calculate the running length total
-		}
-	}
-
-	void init() {
-		if (size() == 0) return;								//return if there aren't any points
-		update_L();
-	}
-
-	/// Returns a stim::vec representing the point at index i
-
-	/// @param i is an index of the desired centerline point
-	stim::vec<T> get_vec(unsigned i){
-		return std::vector< stim::vec3<T> >::at(i);
-	}
-
-	///finds the index of the point closest to the length l on the lower bound.
-	///binary search.
-	size_t findIdx(T l) {
-		for (size_t i = 1; i < L.size(); i++) {				//for each point in the centerline
-			if (L[i] > l) return i - 1;						//if we have passed the desired length value, return i-1
-		}
-		return L.size() - 1;
-		/*size_t i = L.size() / 2;
-		size_t max = L.size() - 1;
-		size_t min = 0;
-		while (i < L.size() - 1){
-			if (l < L[i]) {
-				max = i;
-				i = min + (max - min) / 2;
-			}
-			else if (L[i] <= l && L[i + 1] >= l) {
-				break;
-			}
-			else {
-				min = i;
-				i = min + (max - min) / 2;
-			}
-		}
-		return i;*/
-	}
-
-	///Returns a position vector at the given length into the fiber (based on the pvalue).
-	///Interpolates the radius along the line.
-	///@param l: the location of the in the cylinder.
-	///@param idx: integer location of the point closest to l but prior to it.
-	stim::vec3<T> p(T l, int idx) {
-		T rat = (l - L[idx]) / (L[idx + 1] - L[idx]);
-		stim::vec3<T> v1 = at(idx);
-		stim::vec3<T> v2 = at(idx + 1);
-		return(v1 + (v2 - v1)*rat);
-	}
-
-
-public:
-
-	using std::vector< stim::vec3<T> >::at;
-	using std::vector< stim::vec3<T> >::size;
-
-	centerline() : std::vector< stim::vec3<T> >() {
-		init();
-	}
-	centerline(size_t n) : std::vector< stim::vec3<T> >(n){
-		init();
-	}
-	centerline(std::vector<stim::vec3<T> > pos) :
-		std::vector<stim::vec3<T> > (pos)
-	{
-		init();
-	}
-	
-	//overload the push_back function to update the length vector
-	void push_back(stim::vec3<T> p) {
-		std::vector< stim::vec3<T> >::push_back(p);
-		update_L(size() - 1);
-	}
-
-	///Returns a position vector at the given p-value (p value ranges from 0 to 1).
-	///interpolates the position along the line.
-	///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
-	stim::vec3<T> p(T pvalue) {
-		if (pvalue <= 0.0) return at(0);			//return the first element
-		if (pvalue >= 1.0) return back();			//return the last element
-
-		T l = pvalue*L[L.size() - 1];
-		int idx = findIdx(l);
-		return p(l, idx);
-	}
-
-	///Update centerline internal parameters (currently the L vector)
-	void update() {
-		init();
-	}
-	///Return the length of the entire centerline
-	T length() {
-		return L.back();
-	}
-
-
-	/// stitch two centerlines
-	///@param c1, c2: two centerlines
-	///@param sigma: sample rate
-	static std::vector< stim::centerline<T> > stitch(stim::centerline<T> c1, stim::centerline<T> c2 = stim::centerline<T>()) {
-		
-		std::vector< stim::centerline<T> > result;
-		stim::centerline<T> new_centerline;
-		stim::vec3<T> new_vertex;
-
-		// if only one centerline, stitch itself!
-		if (c2.size() == 0) {
-			size_t num = c1.size();
-			size_t id = 100000;							// store the downsample start position
-			T threshold;
-			if (num < 4) {								// if the number of vertex is less than 4, do nothing
-				result.push_back(c1);
-				return result;
-			}
-			else {
-				// test geometry start vertex
-				stim::vec3<T> v1 = c1[1] - c1[0];		// vector from c1[0] to c1[1]
-				for (size_t p = 2; p < num; p++) {		// 90° standard???
-					stim::vec3<T> v2 = c1[p] - c1[0];
-					float cosine = v2.dot(v1);
-					if (cosine < 0) {
-						id = p;
-						threshold = v2.len();
-						break;
-					}
-				}
-				if (id != 100000) {						// find a downsample position on the centerline
-					T* c;
-					c = (T*)malloc(sizeof(T) * (num - id) * 3);
-					for (size_t p = id; p < num; p++) {
-						for (size_t d = 0; d < 3; d++) {
-							c[p * 3 + d] = c1[p][d];
-						}
-					}
-					stim::kdtree<T, 3> kdt;
-					kdt.create(c, num - id, 5);			// create tree
-
-					T* query = (T*)malloc(sizeof(T) * 3);
-					for (size_t d = 0; d < 3; d++)
-						query[d] = c1[0][d];
-					size_t index;
-					T dist;
-
-					kdt.search(query, 1, &index, &dist);
-
-					free(query);
-					free(c);
-
-					if (dist > threshold) {
-						result.push_back(c1);
-					}
-					else {
-						// the loop part
-						new_vertex = c1[index];
-						new_centerline.push_back(new_vertex);
-						for (size_t p = 0; p < index + 1; p++) {
-							new_vertex = c1[p];
-							new_centerline.push_back(new_vertex);
-						}
-						result.push_back(new_centerline);
-						new_centerline.clear();
-
-						// the tail part
-						for (size_t p = index; p < num; p++) {
-							new_vertex = c1[p];
-							new_centerline.push_back(new_vertex);
-						}
-						result.push_back(new_centerline);
-					}
-				}
-				else {	// there is one potential problem that two positions have to be stitched
-						// test geometry end vertex
-					stim::vec3<T> v1 = c1[num - 2] - c1[num - 1];
-					for (size_t p = num - 2; p > 0; p--) {		// 90° standard
-						stim::vec3<T> v2 = c1[p - 1] - c1[num - 1];
-						float cosine = v2.dot(v1);
-						if (cosine < 0) {
-							id = p;
-							threshold = v2.len();
-							break;
-						}
-					}
-					if (id != 100000) {						// find a downsample position
-						T* c;
-						c = (T*)malloc(sizeof(T) * (id + 1) * 3);
-						for (size_t p = 0; p < id + 1; p++) {
-							for (size_t d = 0; d < 3; d++) {
-								c[p * 3 + d] = c1[p][d];
-							}
-						}
-						stim::kdtree<T, 3> kdt;
-						kdt.create(c, id + 1, 5);				// create tree
-
-						T* query = (T*)malloc(sizeof(T) * 1 * 3);
-						for (size_t d = 0; d < 3; d++)
-							query[d] = c1[num - 1][d];
-						size_t index;
-						T dist;
-
-						kdt.search(query, 1, &index, &dist);
-
-						free(query);
-						free(c);
-
-						if (dist > threshold) {
-							result.push_back(c1);
-						}
-						else {
-							// the tail part
-							for (size_t p = 0; p < index + 1; p++) {
-								new_vertex = c1[p];
-								new_centerline.push_back(new_vertex);
-							}
-							result.push_back(new_centerline);
-							new_centerline.clear();
-
-							// the loop part
-							for (size_t p = index; p < num; p++) {
-								new_vertex = c1[p];
-								new_centerline.push_back(new_vertex);
-							}
-							new_vertex = c1[index];
-							new_centerline.push_back(new_vertex);
-							result.push_back(new_centerline);
-						}
-					}
-					else {	// no stitch position
-						result.push_back(c1);
-					}
-				}
-			}
-		}
-
-
-		// two centerlines
-		else {
-			// find stitch position based on nearest neighbors												
-			size_t num1 = c1.size();
-			T* c = (T*)malloc(sizeof(T) * num1 * 3);		// c1 as reference point
-			for (size_t p = 0; p < num1; p++)				// centerline to array
-				for (size_t d = 0; d < 3; d++)				// because right now my kdtree code is a relatively close code, it has its own structure
-					c[p * 3 + d] = c1[p][d];				// I will merge it into stimlib totally in the near future
-
-			stim::kdtree<T, 3> kdt;							// kdtree object
-			kdt.create(c, num1, 5);							// create tree
-
-			size_t num2 = c2.size();						
-			T* query = (T*)malloc(sizeof(T) * num2 * 3);	// c2 as query point
-			for (size_t p = 0; p < num2; p++) {
-				for (size_t d = 0; d < 3; d++) {
-					query[p * 3 + d] = c2[p][d];
-				}
-			}
-			std::vector<size_t> index(num2);
-			std::vector<T> dist(num2);
-
-			kdt.search(query, num2, &index[0], &dist[0]);	// find the nearest neighbors in c1 for c2
-
-			// clear up
-			free(query);
-			free(c);
-
-			// find the average vertex distance of one centerline
-			T sigma1 = 0;
-			T sigma2 = 0;
-			for (size_t p = 0; p < num1 - 1; p++) 
-				sigma1 += (c1[p] - c1[p + 1]).len();
-			for (size_t p = 0; p < num2 - 1; p++)
-				sigma2 += (c2[p] - c2[p + 1]).len();
-			sigma1 /= (num1 - 1);
-			sigma2 /= (num2 - 1);
-			float threshold = 4 * (sigma1 + sigma2) / 2;			// better way to do this?
-
-			T min_d = *std::min_element(dist.begin(), dist.end());	// find the minimum distance between c1 and c2
-
-			if (min_d > threshold) {								// if the minimum distance is too large
-				result.push_back(c1);
-				result.push_back(c2);
-
-#ifdef DEBUG
-				std::cout << "The distance between these two centerlines is too large" << std::endl;
-#endif
-			}
-			else {
-			//	auto smallest = std::min_element(dist.begin(), dist.end());
-				unsigned int smallest = std::min_element(dist.begin(), dist.end());
-			//	auto i = std::distance(dist.begin(), smallest);		// find the index of min-distance in distance list
-				unsigned int i = std::distance(dist.begin(), smallest);		// find the index of min-distance in distance list
-
-				// stitch position in c1 and c2
-				int id1 = index[i];
-				int id2 = i;
-
-				// actually there are two cases
-				// first one inacceptable
-				// second one acceptable
-				if (id1 != 0 && id1 != num1 - 1 && id2 != 0 && id2 != num2 - 1) {		// only stitch one end vertex to another centerline
-					result.push_back(c1);
-					result.push_back(c2);
-				}
-				else {
-					if (id1 == 0 || id1 == num1 - 1) {			// if the stitch vertex is the first or last vertex of c1
-						// for c2, consider two cases(one degenerate case)
-						if (id2 == 0 || id2 == num2 - 1) {		// case 1, if stitch position is also on the end of c2
-							// we have to decide which centerline get a new vertex, based on direction
-							// for c1, computer the direction change angle
-							stim::vec3<T> v1, v2;
-							float alpha1, alpha2;				// direction change angle
-							if (id1 == 0)
-								v1 = (c1[1] - c1[0]).norm();
-							else
-								v1 = (c1[num1 - 2] - c1[num1 - 1]).norm();
-							v2 = (c2[id2] - c1[id1]).norm();
-							alpha1 = v1.dot(v2);
-							if (id2 == 0)
-								v1 = (c2[1] - c2[0]).norm();
-							else
-								v1 = (c2[num2 - 2] - c2[num2 - 1]).norm();
-							v2 = (c1[id1] - c2[id2]).norm();
-							alpha2 = v1.dot(v2);
-							if (abs(alpha1) > abs(alpha2)) {					// add the vertex to c1 in order to get smooth connection
-								// push back c1
-								if (id1 == 0) {									// keep geometry information
-									new_vertex = c2[id2];
-									new_centerline.push_back(new_vertex);
-									for (size_t p = 0; p < num1; p++) {			// stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
-										new_vertex = c1[p];
-										new_centerline.push_back(new_vertex);
-									}
-								}
-								else {
-									for (size_t p = 0; p < num1; p++) {			// stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1
-										new_vertex = c1[p];
-										new_centerline.push_back(new_vertex);
-									}
-									new_vertex = c2[id2];
-									new_centerline.push_back(new_vertex);
-								}
-								result.push_back(new_centerline);
-								new_centerline.clear();
-
-								// push back c2
-								for (size_t p = 0; p < num2; p++) {
-									new_vertex = c2[p];
-									new_centerline.push_back(new_vertex);
-								}
-								result.push_back(new_centerline);
-							}
-							else {												// add the vertex to c2 in order to get smooth connection
-								// push back c1
-								for (size_t p = 0; p < num1; p++) {
-									new_vertex = c1[p];
-									new_centerline.push_back(new_vertex);
-								}
-								result.push_back(new_centerline);
-								new_centerline.clear();
-
-								// push back c2
-								if (id2 == 0) {									// keep geometry information
-									new_vertex = c1[id1];
-									new_centerline.push_back(new_vertex);
-									for (size_t p = 0; p < num2; p++) {			// stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
-										new_vertex = c2[p];
-										new_centerline.push_back(new_vertex);
-									}
-								}
-								else {
-									for (size_t p = 0; p < num2; p++) {			// stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1
-										new_vertex = c2[p];
-										new_centerline.push_back(new_vertex);
-									}
-									new_vertex = c1[id1];
-									new_centerline.push_back(new_vertex);
-								}
-								result.push_back(new_centerline);
-							}
-						}
-						else {												// case 2, the stitch position is on c2
-							// push back c1
-							if (id1 == 0) {									// keep geometry information
-								new_vertex = c2[id2];
-								new_centerline.push_back(new_vertex);
-								for (size_t p = 0; p < num1; p++) {			// stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
-									new_vertex = c1[p];
-									new_centerline.push_back(new_vertex);
-								}
-							}
-							else {
-								for (size_t p = 0; p < num1; p++) {			// geometry end vertex on c1 -> geometry start vertex on c1 -> stitch vertex on c2
-									new_vertex = c1[p];
-									new_centerline.push_back(new_vertex);
-								}
-								new_vertex = c2[id2];
-								new_centerline.push_back(new_vertex);
-							}
-							result.push_back(new_centerline);
-							new_centerline.clear();
-
-							// push back c2
-							for (size_t p = 0; p < id2 + 1; p++) {			// first part
-								new_vertex = c2[p];
-								new_centerline.push_back(new_vertex);
-							}
-							result.push_back(new_centerline);
-							new_centerline.clear();
-
-							for (size_t p = id2; p < num2; p++) {			// second part
-								new_vertex = c2[p];
-								new_centerline.push_back(new_vertex);
-							}
-							result.push_back(new_centerline);
-						}
-					}
-					else {							// if the stitch vertex is the first or last vertex of c2
-						// push back c2
-						if (id2 == 0) {										// keep geometry information
-							new_vertex = c1[id1];
-							new_centerline.push_back(new_vertex);
-							for (size_t p = 0; p < num2; p++) {				// stitch vertex on c1 -> geometry start vertex on c2 -> geometry end vertex on c2
-								new_vertex = c2[p];
-								new_centerline.push_back(new_vertex);
-							}
-						}
-						else {
-							for (size_t p = 0; p < num2; p++) {				// geometry end vertex on c2 -> geometry start vertex on c2 -> stitch vertex on c1
-								new_vertex = c2[p];
-								new_centerline.push_back(new_vertex);
-							}
-							new_vertex = c1[id1];
-							new_centerline.push_back(new_vertex);
-							result.push_back(new_centerline);
-							new_centerline.clear();
-
-							// push back c1
-							for (size_t p = 0; p < id1 + 1; p++) {			// first part
-								new_vertex = c1[p];
-								new_centerline.push_back(new_vertex);
-							}
-							result.push_back(new_centerline);
-							new_centerline.clear();
-
-							for (size_t p = id1; p < num1; p++) {			// second part
-								new_vertex = c1[p];
-								new_centerline.push_back(new_vertex);
-							}
-							result.push_back(new_centerline);
-						}
-					}
-				}
-			}
-		}
-		return result;
-	}
-
-	/// Split the fiber at the specified index. If the index is an end point, only one fiber is returned
-	std::vector< stim::centerline<T> > split(unsigned int idx){
-
-		std::vector< stim::centerline<T> > fl;				//create an array to store up to two fibers
-		size_t N = size();
-
-		//if the index is an end point, only the existing fiber is returned
-		if(idx == 0 || idx == N-1){
-			fl.resize(1);							//set the size of the fiber to 1
-			fl[0] = *this;							//copy the current fiber
-		}
-
-		//if the index is not an end point
-		else{
-
-			unsigned int N1 = idx + 1;					//calculate the size of both fibers
-			unsigned int N2 = N - idx;
-
-			fl.resize(2);								//set the array size to 2
-
-			fl[0] = stim::centerline<T>(N1);			//set the size of each fiber
-			fl[1] = stim::centerline<T>(N2);
-
-			//copy both halves of the fiber
-			unsigned int i;
-
-			//first half
-			for(i = 0; i < N1; i++)					//for each centerline point
-				fl[0][i] = std::vector< stim::vec3<T> >::at(i);
-			fl[0].init();							//initialize the length vector
-
-			//second half
-			for(i = 0; i < N2; i++)
-				fl[1][i] = std::vector< stim::vec3<T> >::at(idx+i);
-			fl[1].init();							//initialize the length vector
-		}
-
-		return fl;										//return the array
-
-	}
-
-	/// Outputs the fiber as a string
-	std::string str(){
-		std::stringstream ss;
-		size_t N = std::vector< stim::vec3<T> >::size();
-		ss << "---------[" << N << "]---------" << std::endl;
-		for (size_t i = 0; i < N; i++)
-			ss << std::vector< stim::vec3<T> >::at(i) << std::endl;
-		ss << "--------------------" << std::endl;
-
-		return ss.str();
-	}
-
-	/// Back method returns the last point in the fiber
-	stim::vec3<T> back(){
-		return std::vector< stim::vec3<T> >::back();
-	}
-
-		////resample a fiber in the network
-	stim::centerline<T> resample(T spacing)
-	{	
-		//std::cout<<"fiber::resample()"<<std::endl;
-
-		stim::vec3<T> v;    //v-direction vector of the segment
-		stim::vec3<T> p;      //- intermediate point to be added
-		stim::vec3<T> p1;   // p1 - starting point of an segment on the fiber,
-		stim::vec3<T> p2;   // p2 - ending point,
-		//double sum=0;  //distance summation
-
-		size_t N = size();
-
-		centerline<T> new_c; // initialize list of new resampled points on the fiber
-		// for each point on the centerline (skip if it is the last point on centerline)
-		for(unsigned int f=0; f< N-1; f++)
-		{			
-			p1 = at(f); 
-			p2 = at(f+1);
-			v = p2 - p1;
-			
-			T lengthSegment = v.len();			//find Length of the segment as distance between the starting and ending points of the segment
-
-			if(lengthSegment >= spacing){ // if length of the segment is greater than standard deviation resample
-				
-				// repeat resampling until accumulated stepsize is equsl to length of the segment
-				for(T step=0.0; step<lengthSegment; step+=spacing){
-					// calculate the resampled point by travelling step size in the direction of normalized gradient vector
-					p = p1 + v * (step / lengthSegment);
-					
-					// add this resampled points to the new fiber list
-					new_c.push_back(p);
-				}
-			}
-			else       // length of the segment is now less than standard deviation, push the ending point of the segment and proceed to the next point in the fiber
-				new_c.push_back(at(f));
-		}
-		new_c.push_back(at(N-1));   //add the last point on the fiber to the new fiber list
-		//centerline newFiber(newPointList);
-		return new_c;
-	}
-
-};
-
-
-
-}	//end namespace stim
-
-
-
-#endif
+#ifndef JACK_CENTERLINE_H
+#define JACK_CENTERLINE_H
+
+#include <stdlib.h>
+#include <vector>
+#include <stim/math/vec3.h>
+#include <stim/structures/kdtree.cuh>
+
+namespace stim {
+
+	/// we always assume that one centerline has a flow direction even it actually does not have. Also, we allow loop centerline
+	/// NOTE: centerline is not derived from std::vector<stim::vec3<T>> class!!!
+	template<typename T>
+	class centerline {
+
+	private:
+		size_t n;								// number of points on the centerline, can be zero if NULL
+
+		// update length information at each point (distance from starting point) starting from index "start"
+		void update_L(size_t start = 0) {
+			L.resize(n);
+
+			if (start == 0) {
+				L[0] = 0.0f;
+				start++;
+			}
+
+			stim::vec3<T> dir;					// temp direction vector for calculating length between two points
+			for (size_t i = start; i < n; i++) {
+				dir = C[i] - C[i - 1];			// calculate the certerline extending direction
+				L[i] = L[i - 1] + dir.len();	// addition
+			}
+		}
+
+	protected:
+		std::vector<stim::vec3<T> > C;			// points on the centerline
+		std::vector<T> L;						// stores the integrated length along the fiber
+
+	public:
+		/// constructors
+		// empty constructor
+		centerline() {
+			n = 0;
+		}
+
+		// constructor that allocate memory
+		centerline(size_t s) {
+			n = s;
+			C.resize(s);		// allocate memory for points
+			L.resize(s);		// allocate memory for lengths
+
+			update_L();
+		}
+
+		// constructor that constructs a centerline based on a list of points
+		centerline(std::vector<stim::vec3<T> > rhs) {
+			n = rhs.size();						// get the number of points
+			C.resize(n);
+			for (size_t i = 0; i < n; i++)
+				C[i] = rhs[i];					// copy data
+			update_L();
+		}
+
+
+		/// vector operations
+		// add a new point to current centerline
+		void push_back(stim::vec3<T> p) {
+			C.push_back(p);
+			n++;								// increase the number of points
+			update_L(n - 1);
+		}
+
+		// insert a new point at specific location to current centerline
+		void insert(typename std::vector<stim::vec3<T> >::iterator pos, stim::vec3<T> p) {
+			C.insert(pos, p);									// insert a new point
+			n++;
+			size_t d = std::distance(C.begin(), pos);			// get the index
+			update_L(d);
+		}
+		// insert a new point at C[idx]
+		void insert(size_t idx, stim::vec3<T> p) {
+			n++;
+			C.resize(n);
+			for (size_t i = n - 1; i > idx; i--)				// move point location
+				C[i] = C[i - 1];
+			C[idx] = p;
+			update_L(idx);
+		}
+
+		// assign a point at specific location to current centerline
+		void assign(size_t idx, stim::vec3<T> p) {
+			C[idx] = p;
+			update_L(idx);
+		}
+
+		// erase a point at specific location on current centerline
+		void erase(typename std::vector<stim::vec3<T> >::iterator pos) {
+			C.erase(pos);										// erase a point
+			n--;
+			size_t d = std::distance(C.begin(), pos);			// get the index
+			update_L(d);
+		}
+		// erase a point at C[idx]
+		void erase(size_t idx) {
+			n--;
+			for (size_t i = idx; i < n; i++)
+				C[i] = C[i + 1];
+			C.resize(n);
+			update_L(idx);
+		}
+
+		// clear up all the points
+		void clear() {
+			C.clear();			// clear list
+			L.clear();			// clear length information
+			n = 0;				// set number to zero
+		}
+
+		// reverse current centerline in terms of points order
+		centerline<T> reverse() {
+			centerline<T> result = *this;
+
+			std::reverse(result.C.begin(), result.C.end());
+			result.update_L();
+
+			return result;
+		}
+
+		/// functions for reading centerline information
+		// return the number of points on current centerline
+		size_t size() {
+			return n;
+		}
+
+		// return the length
+		T length() {
+			return L.back();
+		}
+
+		// finds the index of the point closest to the length "l" on the lower bound
+		size_t findIdx(T l) {
+			for (size_t i = 1; i < L.size(); i++) {
+				if (L[i] > l)
+					return i - 1;
+			}
+
+			return L.size() - 1;
+		}
+
+		// get a position vector at the given length into the fiber (based on the pvalue), interpolate
+		stim::vec3<T> p(T l, size_t idx) {
+			T rate = (l - L[idx]) / (L[idx + 1] - L[idx]);
+			stim::vec3<T> v1 = C[idx];
+			stim::vec3<T> v2 = C[idx + 1];
+		
+			return (v1 + (v2 - v1) * rate);
+		}
+		// get a position vector at the given pvalue(pvalue[0.0f, 1.0f])
+		stim::vec3<T> p(T pvalue) {
+			// degenerated cases
+			if (pvalue <= 0.0f) return C[0];
+			if (pvalue >= 1.0f) return C.back();
+		
+			T l = pvalue * L.back();		// get the length based on the given pvalue
+			size_t idx = findIdx(l);
+			
+			return p(l, idx);				// interpolation
+		}
+
+		// get the normalized direction vector at point idx (average of the incoming and outgoing directions)
+		stim::vec3<T> d(size_t idx) {
+			if (n <= 1) return stim::vec3<T>(0.0f, 0.0f, 0.0f);		// if there is insufficient information to calculate the direction, return null
+			if (n == 2) return (C[1] - C[0]).norm();				// if there are only two points, the direction vector at both is the direction of the line segment
+			
+			// degenerate cases at two ends
+			if (idx == 0) return (C[1] - C[0]).norm();				// the first direction vector is oriented towards the first line segment
+			if (idx == n - 1) return (C[n - 1] - C[n - 2]).norm();	// the last direction vector is oriented towards the last line segment
+
+			// all other direction vectors are the average direction of the two joined line segments
+			stim::vec3<T> a = C[idx] - C[idx - 1];
+			stim::vec3<T> b = C[idx + 1] - C[idx];
+			stim::vec3<T> ab = a.norm() + b.norm();
+			return ab.norm();
+		}
+
+
+		/// arithmetic operations
+		// '=' operation
+		centerline<T> & operator=(centerline<T> rhs) {
+			n = rhs.n;
+			C = rhs.C;
+			L = rhs.L;
+
+			return *this;
+		}
+
+		// "[]" operation
+		stim::vec3<T> & operator[](size_t idx) {
+			return C[idx];
+		}
+
+		// "==" operation
+		bool operator==(centerline<T> rhs) const {
+
+			if (n != rhs.size())
+				return false;
+			else {
+				size_t num = rhs.size();
+				stim::vec3<T> tmp;			// weird situation that I can only use tmp instead of C itself in comparison
+				for (size_t i = 0; i < num; i++) {
+					stim::vec3<T> tmp = C[i];
+					if (tmp[0] != rhs[i][0] || tmp[1] != rhs[i][1] || tmp[2] != rhs[i][2])
+						return false;
+				}
+				return true;
+			}
+		}
+
+		// "+" operation
+		centerline<T> operator+(stim::vec3<T> p) const {
+			centerline<T> result(*this);
+			result.C.push_back(p);
+			result.n++;
+			result.update_L(n - 1);
+			
+			return result;
+		}
+		centerline<T> operator+(centerline<T> c) const {
+			centerline<T> result(*this);
+			size_t num1 = result.size();
+			size_t num2 = c.size();
+			for (size_t i = 0; i < num2; i++)
+				result.push_back(c[i]);
+			result.update_L(num1);
+
+			return result;
+		}
+
+		
+		/// advanced operation
+		// stitch two centerlines if possible (mutual-stitch and self-stitch)
+		static std::vector<centerline<T> > stitch(centerline<T> c1, centerline<T> c2, size_t end = 0) {
+			std::vector<centerline<T> > result;
+			centerline<T> new_centerline;
+			stim::vec3<T> new_vertex;
+			// ********** for Pavel **********
+			// ********** JACK thinks that ultimately we want it AUTOMATEDLY! **********
+
+			// check stitch case
+			if (c1 == c2) {		// self-stitch case
+				// ***** don't know how it works *****
+			}
+			else {				// mutual-stitch case
+				size_t num1 = c1.size();	// get the numbers of two centerlines
+				size_t num2 = c2.size();
+
+				T* reference = (T*)malloc(sizeof(T) * num1 * 3);	// c1 as reference set
+				T* query = (T*)malloc(sizeof(T) * num2 * 3);		// c2 as query set
+				for (size_t p = 0; p < num1; p++)					// read points
+					for (size_t d = 0; d < 3; d++)
+						reference[p * 3 + d] = c1[p][d];			// KDTREE is stilla close code, it has its own structure
+				for (size_t p = 0; p < num2; p++)					// read points
+					for (size_t d = 0; d < 3; d++)
+						query[p * 3 + d] = c2[p][d];
+
+				stim::kdtree<T, 3> kdt;
+				kdt.create(reference, num1, 5);						// create a tree based on reference set, max_tree_level is set to 5
+				
+				std::vector<size_t> index(num2);					// stores index results
+				std::vector<float> dist(num2);
+
+				kdt.search(query, num2, &index[0], &dist[0]);		// find the nearest neighbor in c1 for c2
+
+				// clear up
+				free(reference);
+				free(query);
+
+				std::vector<float>::iterator sm = std::min_element(dist.begin(), dist.end());	// find smallest distance
+				size_t id = std::distance(dist.begin(), sm);				// find the target index
+
+				size_t id1 = index[id];
+				size_t id2 = id;
+
+				// for centerline c1
+				bool flag = false;						// flag indicates whether it has already added the bifurcation point to corresponding new centerline
+				if (id1 == 0 || id1 == num1 - 1) { 		// splitting bifurcation is on the end
+					new_centerline = c1;
+					new_vertex = c2[id2];
+					if (id1 == 0) {
+						new_centerline.insert(0, new_vertex);
+						flag = true;
+					}
+					if (id1 == num1 - 1) {
+						new_centerline.push_back(new_vertex);
+						flag = true;
+					}
+
+					result.push_back(new_centerline);
+				}
+				else {									// splitting bifurcation is on the centerline
+					std::vector<centerline<T> > tmp_centerline = c1.split(id1);
+					result = tmp_centerline;
+				}
+
+				// for centerline c2
+				if (id2 == 0 || id2 == num2 - 1) {		// splitting bidurcation is on the end
+					new_centerline = c2;
+					if (flag)
+						result.push_back(new_centerline);
+					else {			// add the bifurcation point to this centerline
+						new_vertex = c1[id1];
+						if (id2 == 0) {
+							new_centerline.insert(0, new_vertex);
+							flag = true;
+						}
+						if (id2 == num2 - 1) {
+							new_centerline.push_back(new_vertex);
+							flag = true;
+						}
+
+						result.push_back(new_centerline);
+					}
+				}
+				else {									// splitting bifurcation is on the centerline
+					std::vector<centerline<T> > tmp_centerline = c2.split(id2);
+					result.push_back(tmp_centerline[0]);
+					result.push_back(tmp_centerline[1]);
+				}
+			}
+
+			return result;
+		}
+
+		// split current centerline at specific position
+		std::vector<centerline<T> > split(size_t idx) {
+			std::vector<centerline<T> > result;
+
+			// won't split
+			if (idx <= 0 || idx >= n - 1) {
+				result.resize(1);
+				result[0] = *this;				// return current centerline
+			}
+			// do split
+			else {
+				size_t n1 = idx + 1;			// vertex idx would appear twice
+				size_t n2 = n - idx;			// in total n + 1 points
+
+				centerline<T> tmp;				// temp centerline
+
+				result.resize(2);
+
+				for (size_t i = 0; i < n1; i++)	// first half
+					tmp.push_back(C[i]);
+				tmp.update_L();
+				result[0] = tmp;
+				tmp.clear();					// clear up for next computation
+
+				for (size_t i = 0; i < n2; i++)	// second half
+					tmp.push_back(C[i + idx]);
+				tmp.update_L();
+				result[1] = tmp;
+			}
+
+			return result;
+		}
+
+		// resample current centerline
+		centerline<T> resample(T spacing) {
+			
+			stim::vec3<T> dir;				// direction vector
+			stim::vec3<T> tmp;				// intermiate point to be added
+			stim::vec3<T> p1;				// starting point
+			stim::vec3<T> p2;				// ending point
+
+			centerline<T> result;
+
+			for (size_t i = 0; i < n - 1; i++) {
+				p1 = C[i];
+				p2 = C[i + 1];
+				
+				dir = p2 - p1;				// compute the direction of current segment
+				T seg_len = dir.len();
+
+				if (seg_len > spacing) {	// current segment can be sampled
+					for (T step = 0.0f; step < seg_len; step += spacing) {
+						tmp = p1 + dir * (step / seg_len);		// add new point
+						result.push_back(tmp);
+					}
+				}
+				else
+					result.push_back(p1);	// push back starting point
+			}
+			result.push_back(p2);			// push back ending point
+
+			return result;
+		}
+	};
+}
+
+#endif
-/*
-Copyright <2017> <David Mayerich>
-
-Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
-
-The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-*/
-#ifndef STIM_NETWORK_H
-#define STIM_NETWORK_H
-
-#include <stdlib.h>
-#include <assert.h>
-#include <sstream>
-#include <fstream>
-#include <algorithm>
-#include <string.h>
-#include <math.h>
-#include <stim/math/vec3.h>
-#include <stim/visualization/obj.h>
-#include <stim/visualization/swc.h>
-#include <stim/visualization/cylinder.h>
-#include <stim/cuda/cudatools/timer.h>
-#include <stim/cuda/cudatools/callable.h>
-#include <stim/structures/kdtree.cuh>
-//********************help function********************
-// gaussian_function
-CUDA_CALLABLE float gaussianFunction(float x, float std = 25) { return exp(-x / (2 * std*std)); }  // std default sigma value is 25
-
-// compute metric in parallel
-#ifdef __CUDACC__
-template <typename T>
-__global__ void find_metric_parallel(T* M, size_t n, T* D, float sigma){
-	size_t x = blockDim.x * blockIdx.x + threadIdx.x;
-	if(x >= n) return;
-	M[x] = 1.0f - gaussianFunction(D[x], sigma);
-}
-
-//find the corresponding edge index from array index
-__global__ void find_edge_index_parallel(size_t* I, size_t n, unsigned* R, size_t* E, size_t ne){
-	size_t x = blockDim.x * blockIdx.x + threadIdx.x;
-	if(x >= n) return;
-	unsigned i = 0;
-	size_t N = 0;
-	for(unsigned e = 0; e < ne; e++){
-		N += E[e];
-		if(I[x] < N){
-			R[x] = i;
-			break;
-		}
-		i++;
-	}
-}
-#endif
-
-//hard-coded factor
-int threshold_fac;
-
-namespace stim{
-/** This is the a class that interfaces with gl_spider in order to store the currently
- *   segmented network. The following data is stored and can be extracted:
- *   1)Network geometry and centerline.
- *   2)Network connectivity (a graph of nodes and edges), reconstructed using kdtree.
-*/
-
-template<typename T>
-class network{
-
-	///Each edge is a fiber with two nodes.
-	///Each node is an in index to the endpoint of the fiber in the nodes array.
-	class edge : public cylinder<T>
-	{
-		public:
-	
-		unsigned int v[2];		//unique id's designating the starting and ending
-		// default constructor
-		edge() : cylinder<T>() {
-			v[1] = (unsigned)(-1); v[0] = (unsigned)(-1);
-		}
-		/// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
-/*
-		///@param v0, the starting index.
-		///@param v1, the ending index.
-		///@param sz, the number of point in the fiber.
-		edge(unsigned int v0, unsigned int v1, unsigned int sz) : cylinder<T>(
-		{
-
-		}
-*/
-		edge(std::vector<stim::vec3<T> > p, std::vector<T> s)
-			: cylinder<T>(p,s)
-		{
-		}
-		///@param p is an array of positions in space
-		edge(stim::centerline<T> p) : cylinder<T>(p){}
-
-		/// Copy constructor creates an edge from a cylinder
-		edge(stim::cylinder<T> f) : cylinder<T>(f) {}
-
-		/// Resamples an edge by calling the fiber resampling function
-		edge resample(T spacing){
-			edge e(cylinder<T>::resample(spacing));	//call the fiber->edge constructor
-			e.v[0] = v[0];					//copy the vertex data
-			e.v[1] = v[1];
-
-			return e;						//return the new edge
-		}
-
-		/// Output the edge information as a string
-		std::string str(){
-			std::stringstream ss;
-			ss<<"("<<cylinder<T>::size()<<")\tl = "<<this->length()<<"\t"<<v[0]<<"----"<<v[1];
-			return ss.str();
-		}
-
-		std::vector<edge> split(unsigned int idx){
-		
-			std::vector< stim::cylinder<T> > C;
-			C.resize(2);
-			C =	(*this).cylinder<T>::split(idx);
-			std::vector<edge> E(C.size());
-
-			for(unsigned e = 0; e < E.size(); e++){
-				E[e] = C[e];
-			}
-			return E;
-		}
-
-		/// operator for writing the edge information into a binary .nwt file.
-		friend std::ofstream& operator<<(std::ofstream& out, const edge& e)
-		{
-			out.write(reinterpret_cast<const char*>(&e.v[0]), sizeof(unsigned int));	///write the starting point.
-			out.write(reinterpret_cast<const char*>(&e.v[1]), sizeof(unsigned int));	///write the ending point.
-			unsigned int sz = e.size();	///write the number of point in the edge.
-			out.write(reinterpret_cast<const char*>(&sz), sizeof(unsigned int));
-			for(int i = 0; i < sz; i++)	///write each point
-			{
-				stim::vec3<T> point = e[i];
-				out.write(reinterpret_cast<const char*>(&point[0]), 3*sizeof(T));
-			//	for(int j = 0; j < nmags(); j++)	//future code for multiple mags
-			//	{
-				out.write(reinterpret_cast<const char*>(&e.R[i]), sizeof(T));	///write the radius
-				//std::cout << point.str() << " " << e.R[i] << std::endl;
-			//	}
-			}
-			return out;	//return stream
-		}
-
-		/// operator for reading an edge from a binary .nwt file.
-		friend std::ifstream& operator>>(std::ifstream& in, edge& e)
-		{
-			unsigned int v0, v1, sz;
-			in.read(reinterpret_cast<char*>(&v0), sizeof(unsigned int));	//read the staring point.
-			in.read(reinterpret_cast<char*>(&v1), sizeof(unsigned int));	//read the ending point
-			in.read(reinterpret_cast<char*>(&sz), sizeof(unsigned int));	//read the number of points in the edge
-//			stim::centerline<T> temp = stim::centerline<T>(sz);		//allocate the new edge
-//			e = edge(temp);
-			std::vector<stim::vec3<T> > p(sz);
-			std::vector<T> r(sz);
-			for(int i = 0; i < sz; i++)		//set the points and radii to the newly read values
-			{
-				stim::vec3<T> point;
-				in.read(reinterpret_cast<char*>(&point[0]), 3*sizeof(T));
-				p[i] = point;
-				T mag;
-		//				for(int j = 0; j < nmags(); j++)		///future code for mags
-		//				{
-				in.read(reinterpret_cast<char*>(&mag), sizeof(T));	
-				r[i] = mag;
-				//std::cout << point.str() << " " << mag << std::endl;
-		//				}
-			}
-			e = edge(p,r);
-			e.v[0] = v0; e.v[1] = v1;
-			return in;
-		}
-	};
-
-	///Node class that stores the physical position of the node as well as the edges it is connected to (edges that connect to it), As well as any additional data necessary.
-	class vertex : public stim::vec3<T>
-	{
-		public:
-			//std::vector<unsigned int> edges;					//indices of edges connected to this node.
-			std::vector<unsigned int> e[2];						//indices of edges going out (e[0]) and coming in (e[1])
-			//stim::vec3<T> p;							//position of this node in physical space.
-			//default constructor
-			vertex() : stim::vec3<T>()
-			{
-			}
-			//constructor takes a stim::vec
-			vertex(stim::vec3<T> p) : stim::vec3<T>(p){}
-
-			/// Output the vertex information as a string
-			std::string 
-			str(){
-				std::stringstream ss;
-				ss<<"\t(x, y, z) = "<<stim::vec3<T>::str();
-
-				if(e[0].size() > 0){
-					ss<<"\t> ";
-					for(unsigned int o = 0; o < e[0].size(); o++)
-						ss<<e[0][o]<<" ";
-				}
-				if(e[1].size() > 0){
-					ss<<"\t< ";
-					for(unsigned int i = 0; i < e[1].size(); i++)
-						ss<<e[1][i]<<" ";
-				}
-
-				return ss.str();
-			}
-			///operator for writing the vector into the stream;
-			friend std::ofstream& operator<<(std::ofstream& out, const vertex& v)
-			{
-				unsigned int s0, s1;
-				s0 = v.e[0].size();
-				s1 = v.e[1].size();
-				out.write(reinterpret_cast<const char*>(&v.ptr[0]), 3*sizeof(T));	///write physical vertex location
-				out.write(reinterpret_cast<const char*>(&s0), sizeof(unsigned int));	///write the number of "outgoing edges"
-				out.write(reinterpret_cast<const char*>(&s1), sizeof(unsigned int));	///write the number of "incoming edges"	
-				if (s0 != 0)
-					out.write(reinterpret_cast<const char*>(&v.e[0][0]), sizeof(unsigned int)*v.e[0].size());	///write the "outgoing edges"
-				if (s1 != 0)
-					out.write(reinterpret_cast<const char*>(&v.e[1][0]), sizeof(unsigned int)*v.e[1].size());	///write the "incoming edges"
-				return out;
-			}
-
-			///operator for reading the vector out of the stream;
-			friend std::ifstream& operator>>(std::ifstream& in, vertex& v)
-			{
-				in.read(reinterpret_cast<char*>(&v[0]), 3*sizeof(T));	///read the physical position
-				unsigned int s[2];					
-				in.read(reinterpret_cast<char*>(&s[0]), 2*sizeof(unsigned int));	///read the sizes of incoming and outgoing edge arrays
-
-				std::vector<unsigned int> one(s[0]);
-				std::vector<unsigned int> two(s[1]);
-				v.e[0] = one;
-				v.e[1] = two;
-				if (one.size() != 0)
-					in.read(reinterpret_cast<char*>(&v.e[0][0]), s[0] * sizeof(unsigned int));		///read the arrays of "outgoing edges"
-				if (two.size() != 0)
-					in.read(reinterpret_cast<char*>(&v.e[1][0]), s[1] * sizeof(unsigned int));		///read the arrays of "incoming edges"
-				return in;
-			}
-
-	};
-
-protected:
-
-	std::vector<edge> E;       //list of edges
-	std::vector<vertex> V;	   //list of vertices.
-
-public:
-
-	///default constructor
-	network()
-	{
-		
-	}
-
-	///constructor with a file to load.
-	network(std::string fileLocation)
-	{
-		load_obj(fileLocation);
-	}
-
-	///Returns the number of edges in the network.
-	unsigned int edges(){
-		return E.size();
-	}
-
-	///Returns the number of nodes in the network.
-	unsigned int vertices(){
-		return V.size();
-	}
-
-	///Returns the radius at specific point in the edge
-	T get_r(unsigned e, unsigned i) {
-		return E[e].r(i);
-	}
-
-	///Returns the average radius of specific edge
-	T get_average_r(unsigned e) {
-		T result = 0.0;
-		unsigned n = E[e].size();
-		for (unsigned p = 0; p < n; p++)
-			result += E[e].r(p);
-
-		return (T)result / n;
-	}
-
-	///Returns the length of current edge
-	T get_l(unsigned e) {
-		return E[e].length();
-	}
-
-	///Returns the start vertex of current edge
-	size_t get_start_vertex(unsigned e) {
-		return E[e].v[0];
-	}
-
-	///Returns the end vertex of current edge
-	size_t get_end_vertex(unsigned e) {
-		return E[e].v[1];
-	}
-
-	///Returns one vertex
-	stim::vec3<T> get_vertex(unsigned i) {
-		return V[i];
-	}
-
-	///Returns the boundary vertices' indices
-	std::vector<unsigned> get_boundary_vertex() {
-		std::vector<unsigned> result;
-
-		for (unsigned v = 0; v < V.size(); v++) {
-			if (V[v].e[0].size() + V[v].e[1].size() == 1) {	// boundary vertex
-				result.push_back(v);
-			}
-		}
-
-		return result;
-	}
-
-	///Set radius
-	void set_r(unsigned e, std::vector<T> radius) {
-		E[e].cylinder<T>::copy_r(radius);
-	}
-
-	void set_r(unsigned e, T radius) {
-		for (size_t i = 0; i < E[e].size(); i++)
-			E[e].cylinder<T>::set_r(i, radius);
-	}
-	//scale the network by some constant value
-	//	I don't think these work??????
-	/*std::vector<vertex> operator*(T s){
-		for (unsigned i=0; i< vertices; i ++ ){
-			V[i] = V[i] * s;
-		}
-		return V;
-	}
-
-	std::vector<vertex> operator*(vec<T> s){
-		for (unsigned i=0; i< vertices; i ++ ){
-			for (unsigned dim = 0 ; dim< 3; dim ++){
-				V[i][dim] = V[i][dim] * s[dim];
-			}
-		}
-		return V;
-	}*/
-
-	// Returns an average of branching index in the network
-	double BranchingIndex(){
-		double B=0;
-		for(unsigned v=0; v < V.size(); v ++){
-			B += ((V[v].e[0].size()) + (V[v].e[1].size()));
-		}
-		B = B / V.size();
-		return B;
-
-	}
-
-	// Returns number of branch points in thenetwork
-	unsigned int BranchP(){
-		unsigned int B=0;
-		unsigned int c;
-		for(unsigned v=0; v < V.size(); v ++){
-			c = ((V[v].e[0].size()) + (V[v].e[1].size()));
-			if (c > 2){
-			B += 1;}
-		}		
-		return B;
-
-	}
-
-	// Returns number of end points (tips) in thenetwork
-	unsigned int EndP(){
-		unsigned int B=0;
-		unsigned int c;
-		for(unsigned v=0; v < V.size(); v ++){
-			c = ((V[v].e[0].size()) + (V[v].e[1].size()));
-			if (c == 1){
-			B += 1;}
-		}		
-		return B;
-
-	}
-
-	//// Returns a dictionary with the key as the vertex
-	//std::map<std::vector<vertex>,unsigned int> DegreeDict(){
-	//	std::map<std::vector<vertex>,unsigned int> dd;
-	//	unsigned int c = 0;
-	//	for(unsigned v=0; v < V.size(); v ++){
-	//		c = ((V[v].e[0].size()) + (V[v].e[1].size()));
-	//		dd[V[v]] = c;
-	//	}
-	//	return dd;
-	//}
-
-	//// Return number of branching stems
-	//unsigned int Stems(){
-	//	unsigned int s = 0;
-	//	std::map<std::vector<vertex>,unsigned int> dd;
-	//	dd = DegreeDict();
-	//	//for(unsigned v=0; v < V.size(); v ++){
-	//	//	V[v].e[0].
-	//	return s;
-	//}
-
-	//Calculate Metrics---------------------------------------------------
-	// Returns an average of fiber/edge lengths in the network
-	double Lengths(){
-		stim::vec<T> L;
-		double sumLength = 0;
-		for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
-			L.push_back(E[e].length());						//append the edge length
-			sumLength = sumLength + E[e].length();
-		}
-		double avg = sumLength / E.size();
-		return avg;
-	}
-
-
-	// Returns an average of tortuosities in the network
-	double Tortuosities(){
-		stim::vec<T> t;
-		stim::vec<T> id1, id2;                        // starting and ending vertices of the edge
-		double distance;double tortuosity;double sumTortuosity = 0;
-		for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
-			id1 = E[e][0];									//get the edge starting point
-			id2 = E[e][E[e].size() - 1];					//get the edge ending point
-			distance = (id1 - id2).len();                   //displacement between the starting and ending points
-			if(distance > 0){
-				tortuosity = E[e].length()/	distance	;		// tortuoisty = edge length / edge displacement
-			}
-			else{
-				tortuosity = 0;}
-			t.push_back(tortuosity);
-			sumTortuosity += tortuosity;
-		}
-		double avg = sumTortuosity / E.size();
-		return avg;
-	}
-
-	// Returns average contraction of the network
-	double Contractions(){
-	stim::vec<T> t;
-	stim::vec<T> id1, id2;                        // starting and ending vertices of the edge
-	double distance;double contraction;double sumContraction = 0;
-	for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
-		id1 = E[e][0];									//get the edge starting point
-		id2 = E[e][E[e].size() - 1];					//get the edge ending point
-		distance = (id1 - id2).len();                   //displacement between the starting and ending points
-		contraction = distance / E[e].length();		// tortuoisty = edge length / edge displacement
-		t.push_back(contraction);
-		sumContraction += contraction;
-	}
-	double avg = sumContraction / E.size();
-	return avg;
-	}
-
-	// returns average fractal dimension of the branches of the network
-	double FractalDimensions(){
-	stim::vec<T> t;
-	stim::vec<T> id1, id2;                        // starting and ending vertices of the edge
-	double distance;double fract;double sumFractDim = 0;
-	for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
-		id1 = E[e][0];									//get the edge starting point
-		id2 = E[e][E[e].size() - 1];					//get the edge ending point
-		distance = (id1 - id2).len();                   //displacement between the starting and ending points
-		fract = std::log(distance) / std::log(E[e].length());		// tortuoisty = edge length / edge displacement
-		t.push_back(sumFractDim);
-		sumFractDim += fract;
-	}
-	double avg = sumFractDim / E.size();
-	return avg;
-	}
-
-	//returns a cylinder represented a given fiber (based on edge index)
-	stim::cylinder<T> get_cylinder(unsigned e){
-		return E[e];									//return the specified edge (casting it to a fiber)
-	}
-
-	//load a network from an OBJ file
-	void load_obj(std::string filename){
-
-		stim::obj<T> O;									//create an OBJ object
-		O.load(filename);								//load the OBJ file as an object
-
-		std::vector<unsigned> id2vert;							//this list stores the OBJ vertex ID associated with each network vertex
-
-		unsigned i[2];									//temporary, IDs associated with the first and last points in an OBJ line
-
-		//for each line in the OBJ object
-		for(unsigned int l = 1; l <= O.numL(); l++){
-
-			std::vector< stim::vec<T> > c;						//allocate an array of points for the vessel centerline
-			O.getLine(l, c);							//get the fiber centerline
-
-			stim::centerline<T> c3(c.size());
-			for(size_t j = 0; j < c.size(); j++)
-				c3[j] = c[j];
-			c3.update();
-
-	//		edge new_edge = c3;		///This is dangerous.
-			edge new_edge(c3);
-					
-			//create an edge from the given centerline
-			unsigned int I = new_edge.size();					//calculate the number of points on the centerline
-
-			//get the first and last vertex IDs for the line
-			std::vector< unsigned > id;						//create an array to store the centerline point IDs
-			O.getLinei(l, id);							//get the list of point IDs for the line
-			i[0] = id.front();							//get the OBJ ID for the first element of the line
-			i[1] = id.back();							//get the OBJ ID for the last element of the line
-
-			std::vector<unsigned>::iterator it;					//create an iterator for searching the id2vert array
-			unsigned it_idx;							//create an integer for the id2vert entry index
-
-			//find out if the nodes for this fiber have already been created
-			it = find(id2vert.begin(), id2vert.end(), i[0]);	//look for the first node
-			if(it == id2vert.end()){							//if i[0] hasn't already been used
-				vertex new_vertex = new_edge[0];				//create a new vertex, assign it a position
-				bool flag = false;
-				unsigned j = 0;
-				for (; j < V.size(); j++) {						// check whether current vertex is already exist
-					if (new_vertex == V[j]) {
-						flag = true;
-						break;
-					}
-				}
-				if (!flag) {									// unique one
-					new_vertex.e[0].push_back(E.size());				//add the current edge as outgoing
-					new_edge.v[0] = V.size();					//add the new edge to the edge
-					V.push_back(new_vertex);					//add the new vertex to the vertex list
-					id2vert.push_back(i[0]);					//add the ID to the ID->vertex conversion list
-				}
-				else {
-					V[j].e[0].push_back(E.size());
-					new_edge.v[0] = j;
-				}
-			}
-			else{									//if the vertex already exists
-				it_idx = std::distance(id2vert.begin(), it);
-				V[it_idx].e[0].push_back(E.size());				//add the current edge as outgoing
-				new_edge.v[0] = it_idx;
-			}
-
-			it = find(id2vert.begin(), id2vert.end(), i[1]);			//look for the second ID
-			if(it == id2vert.end()){						//if i[1] hasn't already been used
-				vertex new_vertex = new_edge[I-1];				//create a new vertex, assign it a position
-				bool flag = false;
-				unsigned j = 0;
-				for (; j < V.size(); j++) {					// check whether current vertex is already exist
-					if (new_vertex == V[j]) {
-						flag = true;
-						break;
-					}
-				}
-				if (!flag) {
-					new_vertex.e[1].push_back(E.size());				//add the current edge as incoming
-					new_edge.v[1] = V.size();                                  	//add the new vertex to the edge
-					V.push_back(new_vertex);					//add the new vertex to the vertex list
-					id2vert.push_back(i[1]);					//add the ID to the ID->vertex conversion list
-				}
-				else {
-					V[j].e[1].push_back(E.size());
-					new_edge.v[1] = j;
-				}
-			}
-			else{									//if the vertex already exists
-				it_idx = std::distance(id2vert.begin(), it);
-				V[it_idx].e[1].push_back(E.size());				//add the current edge as incoming
-				new_edge.v[1] = it_idx;
-			}
-
-			E.push_back(new_edge);							//push the edge to the list
-
-		}
-
-		// copy the radii information from OBJ
-		/*if (O.numVT()) {
-			unsigned k = 0;
-			for (unsigned i = 0; i < E.size(); i++) {
-				for (unsigned j = 0; j < E[i].size(); j++) {
-					E[i].cylinder<T>::set_r(j, O.getVT(k)[0] / 2);
-					k++;
-				}
-			}
-		}*/
-		// OBJ class assumes that in L the two values are equal
-		if (O.numVT()) {
-			std::vector< unsigned > id;						//create an array to store the centerline point IDs
-			for (unsigned i = 0; i < O.numL(); i++) {
-				id.clear();
-				O.getLinei(i + 1, id);							//get the list of point IDs for the line
-				for (unsigned j = 0; j < id.size(); j++)
-					E[i].cylinder<T>::set_r(j, O.getVT(id[j] - 1)[0] / 2);
-			}
-		}
-	}
-
-	///loads a .nwt file. Reads the header and loads the data into the network according to the header.
-	void
-	loadNwt(std::string filename)
-	{
-		int dims[2];		///number of vertex, number of edges
-		readHeader(filename, &dims[0]);		//read header
-		std::ifstream file;
-		file.open(filename.c_str(), std::ios::in | std::ios::binary);		///skip header information.
-		file.seekg(14+58+4+4, file.beg);
-		vertex v;
-		for(int i = 0; i < dims[0]; i++)		///for every vertex, read vertex, add to network.
-		{
-			file >> v;
-			V.push_back(v);
-//			std::cout << i << " " << v.str() << std::endl;
-		}
-
-		std::cout << std::endl;
-		for(int i = 0; i < dims[1]; i++)		///for every edge, read edge, add to network.
-		{
-			edge e;
-			file >> e;
-			E.push_back(e);
-			//std::cout << i << " " << E[i].str() << std::endl;		// not necessary?
-		}
-		file.close();
-	}
-
-	///saves a .nwt file. Writes the header in raw text format, then saves the network as a binary file.
-	void
-	saveNwt(std::string filename)
-	{
-		writeHeader(filename);
-		std::ofstream file;
-		file.open(filename.c_str(), std::ios::out | std::ios::binary | std::ios::app);	///since we have written the header we are not appending.
-		for(int i = 0; i < V.size(); i++)	///look through the Vertices and write each one.
-		{
-//			std::cout << i << " " << V[i].str() << std::endl;
-			file << V[i];
-		}
-		for(int i = 0; i < E.size(); i++)	///loop through the Edges and write each one.
-		{
-			//std::cout << i << " " << E[i].str() << std::endl;		// not necesarry?
-			file << E[i];
-		}
-		file.close();
-	}
-
-
-	///Writes the header information to a .nwt file.
-	void
-	writeHeader(std::string filename)
-	{
-		std::string magicString = "nwtFileFormat ";		///identifier for the file.
-		std::string desc = "fileid(14B), desc(58B), #vertices(4B), #edges(4B): bindata";
-		int hNumVertices = V.size();		///int byte header storing the number of vertices in the file
-		int hNumEdges = E.size();		///int byte header storing the number of edges.
-		std::ofstream file;
-		file.open(filename.c_str(), std::ios::out | std::ios::binary);
-		std::cout << hNumVertices << " " << hNumEdges << std::endl;
-		file.write(reinterpret_cast<const char*>(&magicString.c_str()[0]), 14);	//write the file id
-		file.write(reinterpret_cast<const char*>(&desc.c_str()[0]), 58);	//write the description
-		file.write(reinterpret_cast<const char*>(&hNumVertices), sizeof(int));	//write #vert.
-		file.write(reinterpret_cast<const char*>(&hNumEdges), sizeof(int));	//write #edges
-//		file << magicString.c_str() << desc.c_str() << hNumVertices << hNumEdges;
-		file.close();
-		
-	}
-
-	///Reads the header information from a .nwt file.
-	void
-	readHeader(std::string filename, int *dims)
-	{
-		char magicString[14];		///id
-		char desc[58];			///description
-		int hNumVertices;		///#vert
-		int hNumEdges;			///#edges
-		std::ifstream file;		////create stream
-		file.open(filename.c_str(), std::ios::in | std::ios::binary);
-		file.read(reinterpret_cast<char*>(&magicString[0]), 14);	///read the file id.
-		file.read(reinterpret_cast<char*>(&desc[0]), 58);		///read the description
-		file.read(reinterpret_cast<char*>(&hNumVertices), sizeof(int));	///read the number of vertices
-		file.read(reinterpret_cast<char*>(&hNumEdges), sizeof(int));	///read the number of edges
-//		std::cout << magicString << desc << hNumVertices << " " <<  hNumEdges << std::endl;
-		file.close();							///close the file.
-		dims[0] = hNumVertices;						///fill the returned reference.
-		dims[1] = hNumEdges;
-	}
-
-	//load a network from an SWC file
-	void load_swc(std::string filename) {
-		stim::swc<T> S;										// create swc variable
-		S.load(filename);									// load the node information
-		S.create_tree();									// link those node according to their linking relationships as a tree
-		S.resample();
-
-		//NT.push_back(S.node[0].type);						// set the neuronal_type value to the first vertex in the network
-		std::vector<unsigned> id2vert;						// this list stores the SWC vertex ID associated with each network vertex
-		unsigned i[2];										// temporary, IDs associated with the first and last points
-
-		for (unsigned int l = 0; l < S.numE(); l++) {		// for every edge
-			//NT.push_back(S.node[l].type);
-
-			std::vector< stim::vec3<T> > c;
-			S.get_points(l, c);
-
-			stim::centerline<T> c3(c.size());				// new fiber
-			
-			for (unsigned j = 0; j < c.size(); j++)
-				c3[j] = c[j];								// copy the points
-		
-			c3.update();									// update the L information
-			
-			stim::cylinder<T> C3(c3);						// create a new cylinder in order to copy the origin radius information
-			// upadate the R information
-			std::vector<T> radius;
-			S.get_radius(l, radius);
-
-			C3.copy_r(radius);
-
-			edge new_edge(C3);								// new edge	
-
-			//create an edge from the given centerline
-			unsigned int I = new_edge.size();				//calculate the number of points on the centerline
-			
-			//get the first and last vertex IDs for the line
-			i[0] = S.E[l].front();
-			i[1] = S.E[l].back();
-
-			std::vector<unsigned>::iterator it;				//create an iterator for searching the id2vert array
-			unsigned it_idx;								//create an integer for the id2vert entry index
-
-			//find out if the nodes for this fiber have already been created
-			it = find(id2vert.begin(), id2vert.end(), i[0]);	//look for the first node
-			if (it == id2vert.end()) {							//if i[0] hasn't already been used
-				vertex new_vertex = new_edge[0];				//create a new vertex, assign it a position
-				new_vertex.e[0].push_back(E.size());			//add the current edge as outgoing
-				new_edge.v[0] = V.size();						//add the new edge to the edge
-				V.push_back(new_vertex);						//add the new vertex to the vertex list
-				id2vert.push_back(i[0]);						//add the ID to the ID->vertex conversion list
-			}
-			else {									//if the vertex already exists
-				it_idx = std::distance(id2vert.begin(), it);
-				V[it_idx].e[0].push_back(E.size());				//add the current edge as outgoing
-				new_edge.v[0] = it_idx;
-			}
-
-			it = find(id2vert.begin(), id2vert.end(), i[1]);	//look for the second ID
-			if (it == id2vert.end()) {							//if i[1] hasn't already been used
-				vertex new_vertex = new_edge[I - 1];			//create a new vertex, assign it a position
-				new_vertex.e[1].push_back(E.size());			//add the current edge as incoming
-				new_edge.v[1] = V.size();                       //add the new vertex to the edge
-				V.push_back(new_vertex);						//add the new vertex to the vertex list
-				id2vert.push_back(i[1]);						//add the ID to the ID->vertex conversion list
-			}
-			else {									//if the vertex already exists
-				it_idx = std::distance(id2vert.begin(), it);
-				V[it_idx].e[1].push_back(E.size());				//add the current edge as incoming
-				new_edge.v[1] = it_idx;
-			}
-
-			E.push_back(new_edge);								//push the edge to the list
-		}
-	}
-
-	/// Get adjacency matrix of the network
-	std::vector< typename std::vector<int> > get_adj_mat() {
-		
-		unsigned n = V.size();		// get the number of vertices in the networks
-
-		std::vector< typename std::vector<int> > result(n, std::vector<int>(n, 0));	// initialize every entry in the matrix to be 0
-		result.resize(n);			// resize rows
-		for (unsigned i = 0; i < n; i++)
-			result[i].resize(n);	// resize columns
-		
-		for (unsigned i = 0; i < n; i++) {			// for every vertex
-			unsigned num_out = V[i].e[0].size();	// number of outgoing edges of current vertex
-			if (num_out != 0) {
-				for (unsigned j = 0; j < num_out; j++) {
-					int edge_idx = V[i].e[0][j];		// get the jth out-going edge index of current vertex
-					int vertex_idx = E[edge_idx].v[1];	// get the ending vertex of specific out-going edge
-					result[i][vertex_idx] = 1;			// can simply set to 1 if it is simple-graph
-					result[vertex_idx][i] = 1;			// symmetric
-				}
-			}
-		}
-
-		return result;
-	}
-
-	/// Output the network as a string
-	std::string str(){
-
-		std::stringstream ss;
-		ss<<"Nodes ("<<V.size()<<")--------"<<std::endl;
-		for(unsigned int v = 0; v < V.size(); v++){
-			ss<<"\t"<<v<<V[v].str()<<std::endl;
-		}
-
-		ss<<"Edges ("<<E.size()<<")--------"<<std::endl;
-		for(unsigned e = 0; e < E.size(); e++){
-			ss<<"\t"<<e<<E[e].str()<<std::endl;
-		}
-
-		return ss.str();
-	}
-
-	/// This function resamples all fibers in a network given a desired minimum spacing
-	/// @param spacing is the minimum distance between two points on the network
-	stim::network<T> resample(T spacing){
-		stim::network<T> n;								//create a new network that will be an exact copy, with resampled fibers
-		n.V = V;									//copy all vertices
-		//n.NT = NT;										//copy all the neuronal type information
-		n.E.resize(edges());								//allocate space for the edge list
-
-		//copy all fibers, resampling them in the process
-		for(unsigned e = 0; e < edges(); e++){						//for each edge in the edge list
-			n.E[e] = E[e].resample(spacing);					//resample the edge and copy it to the new network
-		}
-
-		return n;							              	//return the resampled network
-	}
-
-	/// Calculate the total number of points on all edges.
-	unsigned total_points(){
-		unsigned n = 0;
-		for(unsigned e = 0; e < E.size(); e++)
-			n += E[e].size();
-		return n;
-	}
-
-	//Copy the point cloud representing the centerline for the network into an array
-	void centerline_cloud(T* dst) {
-		size_t p;										//stores the current edge point
-		size_t P;										//stores the number of points in an edge
-		size_t i = 0;									//index into the output array of points
-		for (size_t e = 0; e < E.size(); e++) {			//for each edge in the network
-			P = E[e].size();							//get the number of points in this edge
-			for (p = 0; p < P; p++) {
-				dst[i * 3 + 0] = E[e][p][0];		
-				dst[i * 3 + 1] = E[e][p][1];
-				dst[i * 3 + 2] = E[e][p][2];
-				i++;
-			}
-		}
-	}
-
-    // convert vec3 to array
-	void stim2array(float *a, stim::vec3<T> b){
-		a[0] = b[0];
-		a[1] = b[1];
-		a[2] = b[2];
-	}
-
-	// convert vec3 to array in bunch
-	void edge2array(T* a, edge b){
-		size_t n = b.size();
-		for(size_t i = 0; i < n; i++){
-			a[i * 3 + 0] = b[i][0];
-			a[i * 3 + 1] = b[i][1];
-			a[i * 3 + 2] = b[i][2];	 
-		}
-	}
-
-	// get list of metric
-	std::vector<T> metric() {
-		std::vector<T> result;
-		T m;
-		for (size_t e = 0; e < E.size(); e++) {
-			for (size_t p = 0; p < E[e].size(); p++) {
-				m = E[e].r(p);
-				result.push_back(m);
-			}
-		}
-		return result;
-	}
-
-	/// Calculate the average magnitude across the entire network.
-	/// @param m is the magnitude value to use. The default is 0 (usually radius).
-	T average(unsigned m = 0){
-
-		T M, L;										//allocate space for the total magnitude and length
-		M = L = 0;									//initialize both the initial magnitude and length to zero
-		for(unsigned e = 0; e < E.size(); e++){						//for each edge in the network
-			M += E[e].integrate();							//get the integrated magnitude
-			L += E[e].length();							//get the edge length
-		}
-
-		return M / L;									//return the average magnitude
-	}
-
-	/// This function compares two networks and returns the percentage of the current network that is missing from A.
-	/// @param A is the network to compare to - the field is generated for A
-	/// @param sigma is the user-defined tolerance value - smaller values provide a stricter comparison
-	stim::network<T> compare(stim::network<T> A, float sigma, int device = -1){
-
-		stim::network<T> R;										//generate a network storing the result of the comparison
-		R = (*this);											//initialize the result with the current network
-
-		T *c;						                 			// centerline (array of double pointers) - points on kdtree must be double
-		size_t n_data = A.total_points();          				// set the number of points
-		c = (T*) malloc(sizeof(T) * n_data * 3);				// allocate an array to store all points in the data set				
-
-		unsigned t = 0;
-		for(unsigned e = 0; e < A.E.size(); e++){				//for each edge in the network
-			for(unsigned p = 0; p < A.E[e].size(); p++){		//for each point in the edge
-				for(unsigned d = 0; d < 3; d++){				//for each coordinate
-
-					c[t * 3 + d] = A.E[e][p][d];				//copy the point into the array c
-				}
-				t++;
-			}
-		}
-
-		//generate a KD-tree for network A
-		size_t MaxTreeLevels = 3;								// max tree level
-		
-#ifdef __CUDACC__
-		cudaSetDevice(device);
-		stim::kdtree<T, 3> kdt;								// initialize a pointer to a kd tree
-
-		kdt.create(c, n_data, MaxTreeLevels);				// build a KD tree
-
-		for(unsigned e = 0; e < R.E.size(); e++){					//for each edge in A
-			//R.E[e].add_mag(0);							//add a new magnitude for the metric
-			//size_t errormag_id = R.E[e].nmags() - 1;		//get the id for the new magnitude
-			
-			size_t n = R.E[e].size();						// the number of points in current edge
-			T* queryPt = new T[3 * n];
-			T* m1 = new T[n];
-			T* dists = new T[n];
-			size_t* nnIdx = new size_t[n];
-
-			T* d_dists;										
-			T* d_m1;										
-			cudaMalloc((void**)&d_dists, n * sizeof(T));
-			cudaMalloc((void**)&d_m1, n * sizeof(T));
-
-			edge2array(queryPt, R.E[e]);
-			kdt.search(queryPt, n, nnIdx, dists);		
-
-			cudaMemcpy(d_dists, dists, n * sizeof(T), cudaMemcpyHostToDevice);					// copy dists from host to device
-
-			// configuration parameters
-			size_t threads = (1024>n)?n:1024;
-			size_t blocks = n/threads + (n%threads)?1:0;
-
-			find_metric_parallel<<<blocks, threads>>>(d_m1, n, d_dists, sigma);					//calculate the metric value based on the distance
-
-			cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost);
-
-			for(unsigned p = 0; p < n; p++){
-				R.E[e].set_r(p, m1[p]);
-			}
-
-			//d_set_mag<<<blocks, threads>>>(R.E[e].M, errormag_id, n, m1);
-		}
-
-#else
-		stim::kdtree<T, 3> kdt;
-		kdt.create(c, n_data, MaxTreeLevels);
-	
-		for(unsigned e = 0; e < R.E.size(); e++){			//for each edge in A
-
-			size_t n = R.E[e].size();						// the number of points in current edge
-			T* query = new T[3 * n];
-			T* m1 = new T[n];
-			T* dists = new T[n];
-			size_t* nnIdx = new size_t[n];
-
-			edge2array(query, R.E[e]);
-
-			kdt.cpu_search(query, n, nnIdx, dists);			//find the distance between A and the current network
-
-			for(unsigned p = 0; p < R.E[e].size(); p++){
-				m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma);	//calculate the metric value based on the distance
-				R.E[e].set_r(p, m1[p]);					//set the error for the second point in the segment
-			}
-		}
-#endif
-		return R;		//return the resulting network
-	}
-
-	/// This function compares two networks and split the current one according to the nearest neighbor of each point in each edge
-	/// @param A is the network to split
-	/// @param B is the corresponding mapping network
-	/// @param sigma is the user-defined tolerance value - smaller values provide a stricter comparison
-	/// @param device is the device that user want to use
-	void split(stim::network<T> A, stim::network<T> B, float sigma, int device, float threshold){
-
-		T *c;						                 	
-		size_t n_data = B.total_points();          				
-		c = (T*) malloc(sizeof(T) * n_data * 3); 				
-
-		size_t NB = B.E.size();								// the number of edges in B
-		unsigned t = 0;
-		for(unsigned e = 0; e < NB; e++){					// for every edge in B			
-			for(unsigned p = 0; p < B.E[e].size(); p++){	// for every points in B.E[e]
-				for(unsigned d = 0; d < 3; d++){				
-
-					c[t * 3 + d] = B.E[e][p][d];			// convert to array
-				}
-				t++;
-			}
-		}
-		size_t MaxTreeLevels = 3;							// max tree level
-
-#ifdef __CUDACC__
-		cudaSetDevice(device);
-		stim::kdtree<T, 3> kdt;								// initialize a pointer to a kd tree
-	
-		//compare each point in the current network to the field produced by A
-		kdt.create(c, n_data, MaxTreeLevels);				// build a KD tree
-
-		std::vector<std::vector<unsigned> > relation;		// the relationship between GT and T corresponding to NN
-		relation.resize(A.E.size());										
-
-		for(unsigned e = 0; e < A.E.size(); e++){			//for each edge in A
-			//A.E[e].add_mag(0);								//add a new magnitude for the metric
-			//size_t errormag_id = A.E[e].nmags() - 1;
-			
-			size_t n = A.E[e].size();						// the number of edges in A
-
-			T* queryPt = new T[3 * n];							// set of all the points in current edge
-			T* m1 = new T[n];								// array of metrics for every point in current edge
-			T* dists = new T[n];							// store the distances for every point in current edge
-			size_t* nnIdx = new size_t[n];					// store the indices for every point in current edge
-			
-			// define pointers in device
-			T* d_dists;														
-			T* d_m1;
-			size_t* d_nnIdx;
-
-			// allocate memory for defined pointers
-			cudaMalloc((void**)&d_dists, n * sizeof(T));
-			cudaMalloc((void**)&d_m1, n * sizeof(T));
-			cudaMalloc((void**)&d_nnIdx, n * sizeof(size_t));
-
-			edge2array(queryPt, A.E[e]);						// convert edge to array
-			kdt.search(queryPt, n, nnIdx, dists);				// search the tree to find the NN for every point in current edge
-
-			cudaMemcpy(d_dists, dists, n * sizeof(T), cudaMemcpyHostToDevice);					// copy dists from host to device
-			cudaMemcpy(d_nnIdx, nnIdx, n * sizeof(size_t), cudaMemcpyHostToDevice);				// copy Idx from host to device
-
-			// configuration parameters
-			size_t threads = (1024>n)?n:1024;													// test to see whether the number of point in current edge is more than 1024
-			size_t blocks = n/threads + (n%threads)?1:0;
-
-			find_metric_parallel<<<blocks, threads>>>(d_m1, n, d_dists, sigma);								// calculate the metrics in parallel
-
-			cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost);
-
-			for(unsigned p = 0; p < n; p++){
-				A.E[e].set_r(p, m1[p]);											// set the error(radius) value to every point in current edge
-			}
-
-			relation[e].resize(n);																// resize every edge relation size
-
-			unsigned* d_relation;
-			cudaMalloc((void**)&d_relation, n * sizeof(unsigned));								// allocate memory
-
-			std::vector<size_t> edge_point_num(NB);												// %th element is the number of points that %th edge has
-			for(unsigned ee = 0; ee < NB; ee++)
-				edge_point_num[ee] = B.E[ee].size();
-
-			size_t* d_edge_point_num;
-			cudaMalloc((void**)&d_edge_point_num, NB * sizeof(size_t));
-			cudaMemcpy(d_edge_point_num, &edge_point_num[0], NB * sizeof(size_t), cudaMemcpyHostToDevice);
-
-			find_edge_index_parallel<<<blocks, threads>>>(d_nnIdx, n, d_relation, d_edge_point_num, NB);			// find the edge corresponding to the array index in parallel
-
-			cudaMemcpy(&relation[e][0], d_relation, n * sizeof(unsigned), cudaMemcpyDeviceToHost);	//copy relationship from device to host
-		}
-#else
-		stim::kdtree<T, 3> kdt;
-		kdt.create(c, n_data, MaxTreeLevels);
-	
-		std::vector<std::vector<unsigned>> relation;		// the mapping relationship between two networks
-		relation.resize(A.E.size());										
-		for(unsigned i = 0; i < A.E.size(); i++)
-			relation[i].resize(A.E[i].size());
-
-		std::vector<size_t> edge_point_num(NB);				//%th element is the number of points that %th edge has
-		for(unsigned ee = 0; ee < NB; ee++)
-			edge_point_num[ee] = B.E[ee].size();
-
-		for(unsigned e = 0; e < A.E.size(); e++){			//for each edge in A
-			
-			size_t n = A.E[e].size();						//the number of edges in A
-
-			T* queryPt = new T[3 * n];
-			T* m1 = new T[n];
-			T* dists = new T[n];							//store the distances
-			size_t* nnIdx = new size_t[n];					//store the indices
-			
-			edge2array(queryPt, A.E[e]);
-			kdt.search(queryPt, n, nnIdx, dists);		
-
-			for(unsigned p = 0; p < A.E[e].size(); p++){
-				m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma);	//calculate the metric value based on the distance
-				A.E[e].set_r(p, m1[p]);									//set the error for the second point in the segment
-				
-				unsigned id = 0;																	//mapping edge's idx
-				size_t num = 0;																		//total number of points before #th edge
-				for(unsigned i = 0; i < NB; i++){
-					num += B.E[i].size();
-					if(nnIdx[p] < num){																//find the edge it belongs to
-						relation[e][p] = id;
-						break;
-					}
-					id++;																			//current edge won't be the one, move to next edge
-				}
-			}
-		}
-#endif
-		E = A.E;
-		V = A.V;
-
-		unsigned int id = 0;									// split value								
-		for(unsigned e = 0; e < E.size(); e++){					// for every edge
-			for(unsigned p = 0; p < E[e].size() - 1; p++){		// for every point in each edge
-				int t = (int)(E[e].length() / sigma * 2);
-				if (t <= 20)
-					threshold_fac = E[e].size();
-				else
-					threshold_fac = (E[e].length() / sigma * 2)/10;
-				if(relation[e][p] != relation[e][p + 1]){		// find the nearest edge changing point
-					id = p + 1;									// if there is no change in NN
-					if(id < threshold_fac || (E[e].size() - id) < threshold_fac)				
-						id = E[e].size() - 1;																			// extreme situation is not acceptable
-					else
-						break;
-				}
-				if(p == E[e].size() - 2)																// if there is no splitting index, set the id to the last point index of current edge
-					id = E[e].size() - 1;
-			}
-			//unsigned errormag_id = E[e].nmags() - 1;
-			T G = 0;																					// test to see whether it has its nearest neighbor
-			for(unsigned i = 0; i < E[e].size(); i++)
-				G += E[e].r(i);																			// won't split special edges
-			if(G / E[e].size() > threshold)															// should based on the color map
-				id = E[e].size() - 1;																	// set split idx to outgoing direction vertex
-
-			std::vector<edge> tmpe;
-			tmpe.resize(2);
-			tmpe = E[e].split(id);
-			vertex tmpv = stim::vec3<T>(-1, -1, 0);														// store the split point as vertex
-			if(tmpe.size() == 2){
-				relation.resize(relation.size() + 1);
-				for(unsigned d = id; d < E[e].size(); d++)
-					relation[relation.size() - 1].push_back(relation[e][d]);
-				tmpe[0].v[0] = E[e].v[0];																// begining vertex of first half edge -> original begining vertex
-				tmpe[1].v[1] = E[e].v[1];																// ending vertex of second half edge -> original ending vertex
-				tmpv = E[e][id];
-				V.push_back(tmpv);
-				tmpe[0].v[1] = (unsigned)V.size() - 1;													// ending vertex of first half edge -> new vertex
-				tmpe[1].v[0] = (unsigned)V.size() - 1;													// begining vertex of second half edge -> new vertex
-				edge tmp(E[e]);
-				E[e] = tmpe[0];																			// replace original edge by first half edge
-				E.push_back(tmpe[1]);																	// push second half edge to the last
-				V[V.size() - 1].e[1].push_back(e);														// push first half edge to the incoming of new vertex
-				V[V.size() - 1].e[0].push_back((unsigned)E.size() - 1);									// push second half edge to the outgoing of new vertex
-				for(unsigned i = 0; i < V[tmp.v[1]].e[1].size(); i++)									// find the incoming edge of original ending vertex
-					if(V[tmp.v[1]].e[1][i] == e)
-						V[tmp.v[1]].e[1][i] = (unsigned)E.size() - 1;									// set to new edge
-			}
-		}
-	}
-
-	/// This function compares two splitted networks and yields a mapping relationship between them according to NN
-	/// @param B is the network that the current network is going to map to
-	/// @param C is the mapping relationship: C[e1] = _e1 means e1 edge in current network is mapping the _e1 edge in B
-	/// @param device is the device that user want to use
-	void mapping(stim::network<T> B, std::vector<unsigned> &C, int device, float threshold){
-		stim::network<T> A;								//generate a network storing the result of the comparison
-		A = (*this);
-
-		size_t n = A.E.size();							// the number of edges in A
-		size_t NB = B.E.size();							// the number of edges in B
-
-		C.resize(A.E.size());	
-
-		T *c;						                 	// centerline (array of double pointers) - points on kdtree must be double
-		size_t n_data = B.total_points();          		// set the number of points
-		c = (T*) malloc(sizeof(T) * n_data * 3); 				
-
-		unsigned t = 0;
-		for(unsigned e = 0; e < NB; e++){					// for each edge in the network
-			for(unsigned p = 0; p < B.E[e].size(); p++){	// for each point in the edge
-				for(unsigned d = 0; d < 3; d++){			// for each coordinate
-
-					c[t * 3 + d] = B.E[e][p][d];
-				}
-				t++;
-			}
-		}
-
-		//generate a KD-tree for network A
-		//float metric = 0.0;                               		// initialize metric to be returned after comparing the network
-		size_t MaxTreeLevels = 3;									// max tree level
-		
-#ifdef __CUDACC__
-		cudaSetDevice(device);
-		stim::kdtree<T, 3> kdt;								// initialize a pointer to a kd tree
-	
-		kdt.create(c, n_data, MaxTreeLevels);				// build a KD tree
-
-		for(unsigned e = 0; e < n; e++){					//for each edge in A
-			//size_t errormag_id = A.E[e].nmags() - 1;		//get the id for the new magnitude
-			
-			//pre-judge to get rid of impossibly mapping edges
-			T M = 0;
-			for(unsigned p = 0; p < A.E[e].size(); p++)
-				M += A.E[e].r(p);
-			M = M / A.E[e].size();
-			if(M > threshold)
-				C[e] = (unsigned)-1;						//set the nearest edge of impossibly mapping edges to maximum of unsigned
-			else{
-				T* queryPt = new T[3];
-				T* dists = new T[1];
-				size_t* nnIdx = new size_t[1];
-
-				stim2array(queryPt, A.E[e][A.E[e].size()/2]);
-				kdt.search(queryPt, 1, nnIdx, dists);
-				
-				unsigned id = 0;							//mapping edge's idx
-				size_t num = 0;								//total number of points before #th edge
-				for(unsigned i = 0; i < NB; i++){
-					num += B.E[i].size();
-					if(nnIdx[0] < num){
-						C[e] = id;
-						break;
-					}
-					id++;
-				}
-			}
-		}
-#else
-		stim::kdtree<T, 3> kdt;
-		kdt.create(c, n_data, MaxTreeLevels);
-		T *dists = new T[1];								// near neighbor distances
-		size_t *nnIdx = new size_t[1];						// near neighbor indices // allocate near neigh indices
-
-		stim::vec3<T> p0, p1;
-		T* queryPt = new T[3];
-
-		for(unsigned int e = 0; e < R.E.size(); e++){			// for each edge in A
-			T M;											// the sum of metrics of current edge
-			for(unsigned p = 0; p < R.E[e].size(); p++)
-				M += A.E[e].r(p);
-			M = M / A.E[e].size();
-			if(M > threshold)								
-				C[e] = (unsigned)-1;
-			else{											// if it should have corresponding edge in B, then...
-				p1 = R.E[e][R.E[e].size()/2];							
-				stim2array(queryPt, p1);
-				kdt.cpu_search(queryPt, 1, nnIdx, dists);	// search the tree		
-				
-				unsigned id = 0;							//mapping edge's idx
-				size_t num = 0;								//total number of points before #th edge
-				for(unsigned i = 0; i < NB; i++){
-					num += B.E[i].size();
-					if(nnIdx[0] < num){
-						C[e] = id;
-						break;
-					}
-					id++;
-				}
-			}
-		}
-#endif
-	}
-
-	/// Returns the number of magnitude values stored in each edge. This should be uniform across the network.
-	//unsigned nmags(){
-	//	return E[0].nmags();
-	//}
-	// split a string in text by the character sep
-	stim::vec<T> split(std::string &text, char sep) 
-	{
-		stim::vec<T> tokens;
-		std::size_t start = 0, end = 0;
-		while ((end = text.find(sep, start)) != std::string::npos) {
-		tokens.push_back(atof(text.substr(start, end - start).c_str()));
-		start = end + 1;
-		}
-		tokens.push_back(atof(text.substr(start).c_str()));
-		return tokens;
-	}
-	// load a network in text file to a network class
-	void load_txt(std::string filename)
-	{
-		std::vector <std::string> file_contents;
-		std::ifstream file(filename.c_str());
-		std::string line;
-		std::vector<unsigned> id2vert;	//this list stores the vertex ID associated with each network vertex
-		//for each line in the text file, store them as strings in file_contents
-		while (std::getline(file, line))
-		{
-			std::stringstream ss(line);
-			file_contents.push_back(ss.str());
-		}
-		unsigned int numEdges = atoi(file_contents[0].c_str()); //number of edges in the network
-		unsigned int I = atoi(file_contents[1].c_str()) ;		//calculate the number of points3d on the first edge
-		unsigned int count = 1; unsigned int k = 2; // count is global counter through the file contents, k is for the vertices on the edges
-		// for each edge in the network.
-		for (unsigned int i = 0; i < numEdges; i ++ )
-		{
-			// pre allocate a position vector p with number of points3d on the edge p
-			std::vector< stim::vec<T> > p(0, I);
-			// for each point on the nth edge
-		  for (unsigned int j = k; j < I + k; j++)
-		  {
-			 // split the points3d of floats with separator space and form a float3 position vector out of them
-			  p.push_back(split(file_contents[j], ' '));
-		  }
-			count += p.size() + 1; // increment count to point at the next edge in the network
-			I = atoi(file_contents[count].c_str()); // read in the points3d at the next edge and convert it to an integer
-			k = count + 1;
-			edge new_edge = p; // create an edge with a vector of points3d  on the edge
-			E.push_back(new_edge); // push the edge into the network
-		}
-		unsigned int numVertices = atoi(file_contents[count].c_str()); // this line in the text file gives the number of distinct vertices
-		count = count + 1; // this line of text file gives the first verrtex
-		// push each vertex into V
-		for (unsigned int i = 0; i < numVertices; i ++)
-		{
-			vertex new_vertex = split(file_contents[count], ' ');
-			V.push_back(new_vertex);
-			count += atoi(file_contents[count + 1].c_str()) + 2; // Skip number of edge ids + 2 to point to the next vertex
-		}
-	} // end load_txt function
-
-	// strTxt returns a string of edges
-	std::string
-	strTxt(std::vector< stim::vec<T> > p)
-	{
-		std::stringstream ss;
-		std::stringstream oss;
-		for(unsigned int i = 0; i < p.size(); i++){
-			ss.str(std::string());
-			for(unsigned int d = 0; d < 3; d++){
-				ss<<p[i][d];
-			}
-			ss << "\n";
-		}
-		return ss.str();
-	}
-	// removes specified character from string
-	void removeCharsFromString(std::string &str, char* charsToRemove ) {
-	   for ( unsigned int i = 0; i < strlen(charsToRemove); ++i ) {
-		  str.erase( remove(str.begin(), str.end(), charsToRemove[i]), str.end() );
-	   }
-	}
-	//exports network to txt file
-	void
-	to_txt(std::string filename)
-	{
-		std::ofstream ofs(filename.c_str(), std::ofstream::out | std::ofstream::app);
-		//int num;
-		ofs << (E.size()).str() << "\n";
-		for(unsigned int i = 0; i < E.size(); i++)
-		{
-			 std::string str;
-			 ofs << (E[i].size()).str() << "\n";
-			 str = E[i].strTxt();
-             ofs << str << "\n";
-		} 
-		for(int i = 0; i < V.size(); i++)
-		{
-			std::string str;
-			str = V[i].str();
-			char temp[4] = "[],";
-			removeCharsFromString(str, temp);
-			ofs << str << "\n";
-		}
-		ofs.close();
-	}
-};		//end stim::network class
-};		//end stim namespace
-#endif
+#ifndef JACK_NETWORK_H
+#define JACK_NETWORK_H
+
+#include "centerline.h"
+#include<stim/visualization/obj.h>
+#include<stim/visualization/swc.h>
+#include <stim/math/circle.h>
+#include <stim/structures/kdtree.cuh>
+
+
+// *****help function*****
+
+template<typename T>
+CUDA_CALLABLE T gaussian(T x, T std = 25.0f) {
+	return exp(-x / (2.0f * std * std));
+}
+
+#ifdef __CUDACC__
+template<typename T>
+__global__ void find_metric(T* M, size_t n, T* D, T sigma) {
+	size_t x = blockDim.x * blockIdx.x + threadIdx.x;
+	if (x >= n) return;		// segfault
+	M[x] = 1.0f - gaussian<T>(D[x], sigma);
+}
+#endif
+
+namespace stim{
+
+	template<typename T>
+	class network {
+	public:
+		// define edge class that extends centerline class with radius information
+		class edge : public stim::centerline<T> {
+		protected:
+			std::vector<T> R;	// radius at each point on current edge
+
+		public:
+			unsigned int v[2];		// unique idx for starting and ending point
+			using stim::centerline<T>::d;
+			using stim::centerline<T>::C;
+
+
+			/// constructors
+			// empty constructor
+			edge() : stim::centerline<T>() {
+				v[0] = UINT_MAX;		// set to default value, risky!
+				v[1] = UINT_MAX;
+			}
+
+			// constructor that contructs an edge based on a centerline
+			edge(stim::centerline<T> c) : stim::centerline<T>(c) {
+				size_t num = c.size();
+				R.resize(num);
+			}
+
+			// constructor that constructs an edge based on a centerline and a list of radii
+			edge(stim::centerline<T> c, std::vector<T> r) : stim::centerline<T>(c) {
+				R = r;			// copy radii
+			}
+
+			// constructor that constructs an edge based on a list of points and a list of radii
+			edge(std::vector<stim::vec3<T> > c, std::vector<T> r) : stim::centerline<T>(c) {
+				R = r;			// copy radii
+			}
+
+
+			/// basic operations
+			// get radius
+			T r(size_t idx) {
+				return R[idx];
+			}
+			
+			// set radius
+			void set_r(T value) {
+				size_t num = R.size();
+				for (size_t i = 0; i < num; i++)
+					R[i] = value;
+			}
+			void set_r(size_t idx, T value) {
+				R[idx] = value;
+			}
+			void set_r(std::vector<T> value) {
+				size_t num = value.size();
+				for (size_t i = 0; i < num; i++)
+					R[i] = value[i];
+			}
+
+			// push back a new radius
+			void push_back_r(T value) {
+				R.push_back(value);
+			}
+
+
+			/// vector operation
+			// insert a new point and radius at specific location
+			void insert(size_t idx, stim::vec3<T> p, T r) {
+				centerline<T>::insert(idx, p);		// insert a new point on current centerline
+
+				R.insert(R.begin() + idx, r);		// insert a new radius for that point
+			}
+
+			// reverse the order of an edge
+			edge reverse() {
+				centerline<T> new_centerline = (*this).centerline<T>::reverse();
+				std::vector<T> new_radius = R;
+				std::reverse(new_radius.begin(), new_radius.end());
+
+				edge result(new_centerline, new_radius);
+
+				return result;
+			}
+
+			// copy edge to array
+			void edge_to_array(T* a) {
+				edge b = (*this);
+				size_t n = b.size();
+				for (size_t i = 0; i < n; i++) {
+					a[i * 3 + 0] = b[i][0];
+					a[i * 3 + 1] = b[i][1];
+					a[i * 3 + 2] = b[i][2];
+				}
+			}
+
+
+			/// arithmetic operations
+			// '+' operation
+			edge operator+(edge e) const {
+				edge result(*this);
+				size_t num = e.size();
+				for (size_t i = 0; i < num; i++) {
+					result.push_back(e[i]);
+					result.push_back_r(e.R[i]);
+				}
+
+				return result;
+			}
+			
+
+			/// advanced operations
+			// concatenate two edges from specific point, The deference between function "+" and "concatenate" is that this requires that new edges should share a joint point
+			edge concatenate(edge e2, size_t p1, size_t p2) {
+				edge e1 = *this;
+				size_t num1 = e1.size();		// get the number of points
+				size_t num2 = e2.size();
+				centerline<T> new_centerline;
+				std::vector<T> new_radius;
+
+				// four situations
+				if (p1 == 0) {
+					if (p2 == 0) 
+						e2 = e2.reverse();
+					for (size_t i = 0; i < num2 - 1; i++) {
+						new_centerline.push_back(e2[i]);
+						new_radius.push_back(e2.R[i]);
+					}
+					for (size_t i = 0; i < num1; i++) {
+						new_centerline.push_back(e1[i]);
+						new_radius.push_back(e1.R[i]);
+					}
+				}
+				else {
+					if (p2 != 0)
+						e2 = e2.reverse();
+					for (size_t i = 0; i < num1 - 1; i++) {
+						new_centerline.push_back(e1[i]);
+						new_radius.push_back(e1.R[i]);
+					}
+					for (size_t i = 0; i < num2; i++) {
+						new_centerline.push_back(e2[i]);
+						new_radius.push_back(e2.R[i]);
+					}
+				}
+
+				edge result(new_centerline, new_radius);
+
+				return result;
+			}
+
+			// split current edge at specific position
+			std::vector<edge> split(size_t idx) {
+
+				// can't update v!!!
+				std::vector<centerline<T> > tmp;
+				tmp = (*this).centerline<T>::split(idx);			// split current edge in terms of centerline
+				size_t num = tmp.size();
+				std::vector<edge> result(num);						// construct a list of edges
+				for (size_t i = 0; i < num; i++) {
+					edge new_edge(tmp[i]);							// construct new edge based on one centerline
+					result[i] = new_edge;
+				}
+		
+				for (size_t i = 0; i < num; i++) {					// for every edge
+					for (size_t j = 0; j < result[i].size(); j++) {	// for every point on that edge
+						result[i].R[j] = R[j + i * idx];			// update radius information
+					}
+				}
+
+				return result;
+			}
+
+			// resample current edge
+			edge resample(T spacing) {
+				edge result(centerline<T>::resample(spacing));		// resample current edge and output as a new edge
+				result.v[0] = v[0];									// updates unique indices
+				result.v[1] = v[1];
+
+				return result;
+			}
+
+			// compute a circle that represent the shape of cylinder cross section at point idx, (TGC -> truncated generalized cones)
+			stim::circle<T> circ(size_t idx) {
+				
+				stim::circle<T> c;			// create a circle to orient for finding the circle plane at point idx
+				c.rotate(d(idx));			// rotate the circle
+				stim::vec3<T> U = c.U;		// copy the frenet frame vector
+
+				return stim::circle<T>(C[idx], R[idx], d(idx), U);
+			}
+
+			/// output operation
+			// output the edge information as a string
+			std::string str() {
+				std::stringstream ss;
+				ss << "(" << centerline<T>::size() << ")\tl = " << this->length() << "\t" << v[0] << "----" << v[1];
+				return ss.str();
+			}
+
+			/// operator for writing the edge information into a binary .nwt file.
+			friend std::ofstream& operator<<(std::ofstream& out, edge& e) {
+				out.write(reinterpret_cast<const char*>(&e.v[0]), sizeof(unsigned int));	// write the starting point.
+				out.write(reinterpret_cast<const char*>(&e.v[1]), sizeof(unsigned int));	// write the ending point.
+				size_t sz = e.size();				// write the number of point in the edge.
+				out.write(reinterpret_cast<const char*>(&sz), sizeof(unsigned int));
+				for (size_t i = 0; i < sz; i++) {	// write each point
+					stim::vec3<T> point = e[i];
+					out.write(reinterpret_cast<const char*>(&point[0]), 3 * sizeof(T));
+					out.write(reinterpret_cast<const char*>(&e.R[i]), sizeof(T));			// write the radius
+				}
+				return out;	// return stream
+			}
+
+			/// operator for reading an edge from a binary .nwt file.
+			friend std::ifstream& operator>>(std::ifstream& in, edge& e) {
+				unsigned int v0, v1, sz;
+				in.read(reinterpret_cast<char*>(&v0), sizeof(unsigned int));	// read the staring point.
+				in.read(reinterpret_cast<char*>(&v1), sizeof(unsigned int));	// read the ending point
+				in.read(reinterpret_cast<char*>(&sz), sizeof(unsigned int));	// read the number of points in the edge
+
+				std::vector<stim::vec3<T> > p(sz);
+				std::vector<T> r(sz);
+				for (size_t i = 0; i < sz; i++) {	// set the points and radii to the newly read values
+					stim::vec3<T> point;
+					in.read(reinterpret_cast<char*>(&point[0]), 3 * sizeof(T));
+					p[i] = point;
+					T mag;
+					in.read(reinterpret_cast<char*>(&mag), sizeof(T));
+					r[i] = mag;
+				}
+				e = edge(p, r);
+				e.v[0] = v0; e.v[1] = v1;
+				return in;
+			}
+		};
+
+		// define vertex class that extends vec3 class with connectivity information
+		class vertex : public stim::vec3<T> {
+		public:
+			std::vector<unsigned int> e[2];	// incoming and outgoing edges of that vertex
+			using stim::vec3<T>::ptr;
+
+			/// constructors
+			// empty constructor
+			vertex() : stim::vec3<T>() {
+			}
+
+			// constructor that constructs a vertex based on a vec3 vector
+			vertex(stim::vec3<T> v) : stim::vec3<T>(v) {
+			}
+
+			stim::vec3<T>
+			getPosition()
+			{
+				return stim::vec3<T>(ptr[0], ptr[1], ptr[2]);
+			}
+
+			/// output operation
+			// output the vertex information as a string
+			std::string str() {
+				std::stringstream ss;
+				ss << "\t(x, y, z) = " << stim::vec3<T>::str();
+
+				if (e[0].size() > 0) {
+					ss << "\t> ";
+					for (size_t i = 0; i < e[0].size(); i++)
+						ss << e[0][i] << " ";
+				}
+				if (e[1].size() > 0) {
+					ss << "\t< ";
+					for (size_t i = 0; i < e[1].size(); i++)
+						ss << e[1][i] << " ";
+				}
+			
+				return ss.str();
+			}
+
+			/// operator for writing the vector into the stream;
+			friend std::ofstream& operator<<(std::ofstream& out, const vertex& v) {
+				unsigned int s0, s1;
+				s0 = v.e[0].size();
+				s1 = v.e[1].size();
+				out.write(reinterpret_cast<const char*>(&v.ptr[0]), 3 * sizeof(T));		// write physical vertex location
+				out.write(reinterpret_cast<const char*>(&s0), sizeof(unsigned int));			// write the number of "outgoing edges"
+				out.write(reinterpret_cast<const char*>(&s1), sizeof(unsigned int));			// write the number of "incoming edges"	
+				if (s0 != 0)
+					out.write(reinterpret_cast<const char*>(&v.e[0][0]), sizeof(unsigned int)*v.e[0].size());	// write the "outgoing edges"
+				if (s1 != 0)
+					out.write(reinterpret_cast<const char*>(&v.e[1][0]), sizeof(unsigned int)*v.e[1].size());	// write the "incoming edges"
+				return out;
+			}
+
+			/// operator for reading the vector out of the stream;
+			friend std::ifstream& operator >> (std::ifstream& in, vertex& v) {
+				in.read(reinterpret_cast<char*>(&v[0]), 3 * sizeof(T));				// read the physical position
+				unsigned int s[2];
+				in.read(reinterpret_cast<char*>(&s[0]), 2 * sizeof(unsigned int));	// read the sizes of incoming and outgoing edge arrays
+
+				std::vector<unsigned int> one(s[0]);
+				std::vector<unsigned int> two(s[1]);
+				v.e[0] = one;
+				v.e[1] = two;
+				if (one.size() != 0)
+					in.read(reinterpret_cast<char*>(&v.e[0][0]), s[0] * sizeof(unsigned int));		// read the arrays of "outgoing edges"
+				if (two.size() != 0)
+					in.read(reinterpret_cast<char*>(&v.e[1][0]), s[1] * sizeof(unsigned int));		// read the arrays of "incoming edges"
+				return in;
+			}
+		};
+
+	public:
+
+		std::vector<edge> E;	// list of edges
+		std::vector<vertex> V;	// list of vertices
+
+		/// constructors
+		// empty constructor
+		network() {
+		}
+
+		// constructor with a file to load
+		network(std::string fileLocation) {
+			load_obj(fileLocation);
+		}
+
+		// constructor that constructs a network based on lists of vertices and edges
+		network(std::vector<edge> nE, std::vector<vertex> nV) {
+			E = nE;
+			V = nV;
+		}
+
+
+		/// basic operations
+		// get the number of edges
+		size_t edges() {
+			return E.size();
+		}
+
+		// get the number of vertices
+		size_t vertices() {
+			return V.size();
+		}
+
+		// get the radius at specific point
+		T r(size_t f, size_t p) {		// edge f, point p
+			return E[f].r(p);
+		}
+		T r(size_t c) {					// vertex c
+			T result;
+			if (V[c].e[0].size()) {				// if this vertex has outgoing edges
+				size_t f = V[c].e[0][0];		// get the index of first outgoing edge of this vertex
+				result = r(f, 0);				// this vertex should be the starting point of that edge
+			}
+			else {								// if this vertex only has incoming edges
+				size_t f = V[c].e[1][0];		// get the index of first incoming edge of this vertex
+				result = r(f, E[f].size() - 1);	// this vertex should be the ending point of that edge
+			}
+
+			return result;
+		}
+
+		// get the average radius of one specific edge
+		T ar(size_t f) {
+			T result = 0.0f;
+			size_t num = E[f].size();
+			for (size_t i = 0; i < num; i++)
+				result += E[f].R[i];
+			result = result / num;
+
+			return result;
+		}
+
+		// get the length of edge "f"
+		T length(size_t f) {
+			return E[f].length();
+		}
+
+		// copy specific edge
+		edge get_edge(size_t f) {
+			return E[f];
+		}
+
+		// copy specific vertex
+		vertex get_vertex(size_t c) {
+			return V[c];
+		}
+
+		// get boundary/pendant vertices
+		std::vector<size_t> pendant() {
+			std::vector<size_t> result;
+
+			for (size_t i = 0; i < V.size(); i++)
+				if (V[i].e[0].size() + V[i].e[1].size() == 1)
+					result.push_back(i);
+			
+			return result;
+		}
+
+		// set radius for specific point on edge "f"
+		void set_r(size_t f, size_t p, T value) {
+			E[f].set_r(p, value);
+		}
+		void set_r(size_t f, T value) {
+			E[f].set_r(value);
+		}
+		void set_r(size_t f, std::vector<T> value) {
+			E[f].set_r(value);
+		}
+
+		// copy all points (coordinates) to 1D array
+		void copy_to_array(T* dst) {
+			size_t t = 0;								// indicator for points
+			for (size_t i = 0; i < E.size(); i++) {
+				for (size_t j = 0; j < E[i].size(); j++) {
+					for (size_t k = 0; k < 3; k++) {
+						dst[t * 3 + k] = E[i][j][k];
+					}
+					t++;			// next point
+				}
+			}
+		}
+
+		// get an average of branching index in the network
+		T BranchingIndex() {
+			T B = 0.0f;
+			size_t num = V.size();
+			for (size_t i = 0; i < num; i++)
+				B += (T)(V[i].e[0].size() + V[i].e[1].size());
+
+			B = B / (T)num;
+
+			return B;
+		}
+
+		// get the number of branching points in the network
+		size_t BranchP() {
+			size_t B = 0;
+			size_t c;
+			size_t num = V.size();
+			for (size_t i = 0; i < num; i++) {
+				c = (V[i].e[0].size() + V[i].e[1].size());
+				if (c > 2)
+					B += 1;
+			}
+
+			return B;
+		}
+
+		// get the number of starting or ending points in the network
+		size_t EndP() {
+			size_t B = 0;
+			size_t c;
+			size_t num = V.size();
+			for (size_t i = 0; i < num; i++) {
+				c = (V[i].e[0].size() + V[i].e[1].size());
+				if (c == 1)
+					B += 1;
+			}
+
+			return B;
+		}
+
+		// get an average of fiber length in the network
+		T Lengths() {
+			std::vector<T> L;
+			T sumLength = 0.0f;
+			size_t num = E.size();
+			for (size_t i = 0; i < num; i++) {
+				L.push_back(E[i].length());
+				sumLength += E(i).length();
+			}
+			T avg = sumLength / (T)num;
+			
+			return avg;
+		}
+
+		// get the total number of points in the network
+		size_t total_points() {
+			size_t n = 0;
+			size_t num = E.size();
+			for (size_t i = 0; i < num; i++)
+				n += num;
+
+			return n;
+		}
+
+		// get an average of tortuosities in the network
+		T Tortuosities() {
+			std::vector<T> t;
+			std::vector<T> id1, id2;
+			T distance;
+			T tortuosity;
+			T sumTortuosity = 0.0f;
+			size_t num = E.size();
+			for (size_t i = 0; i < num; i++) {
+				id1 = E[i][0];						// get the starting point
+				id2 = E[i][num - 1];				// get the ending point
+				distance = (id1 - id2).len();
+				if (distance > 0)
+					tortuosity = E[i].length() / distance;	// tortuosity = edge length / edge displacement
+				else
+					tortuosity = 0.0f;
+				t.push_back(tortuosity);
+				sumTortuosity += tortuosity;
+			}
+			T avg = sumTortuosity / (T)num;
+
+			return avg;
+		}
+
+		// get an average contraction of the network
+		T Contraction() {
+			std::vector<T> t;
+			std::vector<T> id1, id2;							// starting and ending vertices of the edge
+			T distance;
+			T contraction;
+			T sumContraction = 0.0f;
+			size_t num = E.size();
+			for (size_t i = 0; i < num; i++) {					// for each edge in the network
+				id1 = E[i][0];									// get the edge starting point
+				id2 = E[i][num - 1];							// get the edge ending point
+				distance = (id1 - id2).len();                   // displacement between the starting and ending points
+				contraction = distance / E[i].length();			// contraction = edge displacement / edge length
+				t.push_back(contraction);
+				sumContraction += contraction;
+			}
+			T avg = sumContraction / (T)num;
+
+			return avg;
+		}
+
+		// get an average fractal dimension of the branches of the network
+		T FractalDimensions() {
+			std::vector<T> t;
+			std::vector<T> id1, id2;							// starting and ending vertices of the edge
+			T distance;
+			T fract;
+			T sumFractDim = 0.0f;
+			size_t num = E.size();
+			for (size_t i = 0; i < num; i++) {							// for each edge in the network
+				id1 = E[i][0];											// get the edge starting point
+				id2 = E[i][num - 1];									// get the edge ending point
+				distance = (id1 - id2).len();							// displacement between the starting and ending points
+				fract = std::log(distance) / std::log(E[i].length());	// fractal dimension = log(edge displacement) / log(edge length)
+				t.push_back(sumFractDim);
+				sumFractDim += fract;
+			}
+			T avg = sumFractDim / (T)num;
+			
+			return avg;
+		}
+
+		
+		/// construct network from files
+		// load network from OBJ files
+		void load_obj(std::string filename) {
+			stim::obj<T> O;				// create OBJ object
+			O.load(filename);			// load OBJ file to an object
+
+			size_t ii[2];				// starting/ending point index of one centerline/edge
+			std::vector<size_t> index;	// added vertex index
+			std::vector<size_t>::iterator it;	// iterator for searching
+			size_t pos;					// position of added vertex
+
+			// get the points
+			for (size_t l = 1; l <= O.numL(); l++) {	// for every line of points
+				std::vector<stim::vec<T> > tmp;			// temp centerline
+				O.getLine(l, tmp);						// get points
+				size_t n = tmp.size();
+				std::vector<stim::vec3<T> > c(n);
+				for (size_t i = 0; i < n; i++) {		// switch from vec to vec3
+					for (size_t j = 0; j < 3; j++) {
+						c[i][j] = tmp[i][j];
+					}
+				}
+
+				centerline<T> C(c);				// construct centerline
+				edge new_edge(C);				// construct edge without radii
+
+				std::vector<unsigned> id;		// temp point index
+				O.getLinei(l, id);				// get point index
+				
+				ii[0] = (size_t)id.front();
+				ii[1] = (size_t)id.back();
+
+				size_t num = new_edge.size();	// get the number of point on current edge
+
+				// for starting point
+				it = std::find(index.begin(), index.end(), ii[0]);
+				if (it == index.end()) {		// new vertex
+					vertex new_vertex = new_edge[0];
+					new_vertex.e[0].push_back(E.size());	// push back the outgoing edge index
+					new_edge.v[0] = V.size();				// save the starting vertex index
+					V.push_back(new_vertex);
+					index.push_back(ii[0]);					// save the point index
+				}
+				else {							// already added vertex
+					pos = std::distance(index.begin(), it);	// find the added vertex position
+					V[pos].e[0].push_back(E.size());		// push back the outgoing edge index
+					new_edge.v[0] = pos;
+				}
+
+				// for ending point
+				it = std::find(index.begin(), index.end(), ii[1]);
+				if (it == index.end()) {		// new vertex
+					vertex new_vertex = new_edge[num - 1];
+					new_vertex.e[1].push_back(E.size());	// push back the incoming edge index
+					new_edge.v[1] = V.size();				// save the ending vertex index
+					V.push_back(new_vertex);
+					index.push_back(ii[1]);					// save the point index
+				}
+				else {							// already added vertex
+					pos = std::distance(index.begin(), it);	// find the added vertex position
+					V[pos].e[1].push_back(E.size());		// push back the incoming edge index
+					new_edge.v[1] = pos;
+				}
+
+				E.push_back(new_edge);
+			}
+
+			// get the radii
+			if (O.numVT()) {		// copy radii information if provided
+				std::vector<unsigned> id;		// a list stores the indices of points
+				for (size_t i = 1; i <= O.numL(); i++) {
+					id.clear();		// clear up temp for current round computation
+					O.getLinei(i, id);
+					size_t num = id.size();	// get the number of points
+					T radius;
+					for (size_t j = 0; j < num; j++) {
+						radius = O.getVT(id[j] - 1)[0] / 2;		// hard-coded: radius = diameter / 2
+						set_r(i - 1, j, radius);				// copy the radius
+					}
+				}
+			}
+		}
+
+		// load network from SWC files (neuron)
+		void load_swc(std::string filename) {
+			stim::swc<T> S;				// create a SWC object
+			S.load(filename);			// load data from SWC file to an object
+			S.create_tree();			// link nodes according to connectivity as a tree
+			S.resample();
+
+			size_t i[2];				// starting/ending point index of one centerline/edge
+			std::vector<size_t> index;	// added vertex index
+			std::vector<size_t>::iterator it;	// iterator for searching
+			size_t pos;					// position of added vertex
+
+			for (size_t l = 0; l < S.numE; l++) {
+				std::vector<stim::vec<T> > c;	// temp centerline
+				S.get_points(l, c);
+
+				centerline<T> C(c);				// construct centerline
+
+				std::vector<T> radius;
+				S.get_radius(l, radius);		// get radius
+
+				edge new_edge(C, radius);		// construct edge
+				size_t num = new_edge.size();	// get the number of point on current edge
+
+				i[0] = S.E[l].front();
+				i[1] = S.E[l].back();
+
+				// for starting point
+				it = std::find(index.begin(), index.end(), i[0]);
+				if (it == index.end()) {		// new vertex
+					vertex new_vertex = new_edge[0];
+					new_vertex.e[0].push_back(E.size());	// push back the outgoing edge index
+					new_edge.v[0] = V.size();				// save the starting vertex index
+					V.push_back(new_vertex);
+					index.push_back(i[0]);					// save the point index
+				}
+				else {							// already added vertex
+					pos = std::distance(index.begin(), it);	// find the added vertex position
+					V[pos].e[0].push_back(E.size());		// push back the outgoing edge index
+					new_edge.v[0] = pos;
+				}
+
+				// for ending point
+				it = std::find(index.begin(), index.end(), i[1]);
+				if (it == index.end()) {		// new vertex
+					vertex new_vertex = new_edge[num - 1];
+					new_vertex.e[1].push_back(E.size());	// push back the incoming edge index
+					new_edge.v[1] = V.size();				// save the ending vertex index
+					V.push_back(new_vertex);
+					index.push_back(i[1]);					// save the point index
+				}
+				else {							// already added vertex
+					pos = std::distance(index.begin(), it);	// find the added vertex position
+					V[pos].e[1].push_back(E.size());		// push back the incoming edge index
+					new_edge.v[1] = pos;
+				}
+
+				E.push_back(new_edge);
+			}
+		}
+
+/*
+		// load a network in text file to a network class
+		void load_txt(std::string filename) {
+			std::vector <std::string> file_contents;
+			std::ifstream file(filename.c_str());
+			std::string line;
+			std::vector<size_t> id2vert;	// this list stores the vertex ID associated with each network vertex
+											// for each line in the text file, store them as strings in file_contents
+			while (std::getline(file, line)) {
+				std::stringstream ss(line);
+				file_contents.push_back(ss.str());
+			}
+			size_t numEdges = atoi(file_contents[0].c_str());	// number of edges in the network
+			size_t I = atoi(file_contents[1].c_str());			// calculate the number of points3d on the first edge
+			size_t count = 1; size_t k = 2;						// count is global counter through the file contents, k is for the vertices on the edges
+
+			for (size_t i = 0; i < numEdges; i++) {
+				// pre allocate a position vector p with number of points3d on the edge p
+				std::vector<stim::vec<T> > p(0, I);
+				// for each point on the nth edge
+				for (size_t j = k; j < I + k; j++) {
+					// split the points3d of floats with separator space and form a float3 position vector out of them
+					p.push_back(std::split(file_contents[j], ' '));
+				}
+				count += p.size() + 1;	// increment count to point at the next edge in the network
+				I = atoi(file_contents[count].c_str()); // read in the points3d at the next edge and convert it to an integer
+				k = count + 1;
+				edge new_edge = p;		// create an edge with a vector of points3d  on the edge
+				E.push_back(new_edge);	// push the edge into the network
+			}
+			size_t numVertices = atoi(file_contents[count].c_str()); // this line in the text file gives the number of distinct vertices
+			count = count + 1;			// this line of text file gives the first verrtex
+	
+			for (size_t i = 0; i < numVertices; i++) {
+				vertex new_vertex = std::split(file_contents[count], ' ');
+				V.push_back(new_vertex);
+				count += atoi(file_contents[count + 1].c_str()) + 2; // Skip number of edge ids + 2 to point to the next vertex
+			}
+		}
+*/
+		// load network from NWT files
+		void load_nwt(std::string filename) {
+			int dims[2];						// number of vertex, number of edges
+			read_nwt_header(filename, &dims[0]);		// read header
+			std::ifstream file;
+			file.open(filename.c_str(), std::ios::in | std::ios::binary);		// skip header information.
+			file.seekg(14 + 58 + 4 + 4, file.beg);
+			vertex v;
+			for (int i = 0; i < dims[0]; i++) {	// for every vertex, read vertex, add to network.
+				file >> v;
+				std::cerr << v.str() << std::endl;
+				V.push_back(v);
+			}
+
+			std::cout << std::endl;
+			for (int i = 0; i < dims[1]; i++) {	// for every edge, read edge, add to network.
+				edge e;
+				file >> e;
+				std::cerr << e.str() << std::endl;
+				E.push_back(e);
+			}
+			file.close();
+		}
+
+		// save network to NWT files
+		void save_nwt(std::string filename) {
+			write_nwt_header(filename);
+			std::ofstream file;
+			file.open(filename.c_str(), std::ios::out | std::ios::binary | std::ios::app);	///since we have written the header we are not appending.
+			for (int i = 0; i < V.size(); i++) {	// look through the Vertices and write each one.
+				file << V[i];
+			}
+			for (int i = 0; i < E.size(); i++) {	// loop through the Edges and write each one.
+				file << E[i];
+			}
+			file.close();
+		}
+
+		/// NWT format functions
+		void read_nwt_header(std::string filename, int *dims) {
+			char magicString[14];		// id
+			char desc[58];				// description
+			int hNumVertices;			// #vert
+			int hNumEdges;				// #edges
+			std::ifstream file;			// create stream
+			file.open(filename.c_str(), std::ios::in | std::ios::binary);
+			file.read(reinterpret_cast<char*>(&magicString[0]), 14);		// read the file id.
+			file.read(reinterpret_cast<char*>(&desc[0]), 58);				// read the description
+			file.read(reinterpret_cast<char*>(&hNumVertices), sizeof(int));	// read the number of vertices
+			file.read(reinterpret_cast<char*>(&hNumEdges), sizeof(int));	// read the number of edges
+			file.close();								// close the file.
+			dims[0] = hNumVertices;						// fill the returned reference.
+			dims[1] = hNumEdges;
+		}
+		void write_nwt_header(std::string filename) {
+			std::string magicString = "nwtFileFormat ";				// identifier for the file.
+			std::string desc = "fileid(14B), desc(58B), #vertices(4B), #edges(4B): bindata";
+			int hNumVertices = V.size();							// int byte header storing the number of vertices in the file
+			int hNumEdges = E.size();								// int byte header storing the number of edges.
+			std::ofstream file;
+			file.open(filename.c_str(), std::ios::out | std::ios::binary);
+			std::cout << hNumVertices << " " << hNumEdges << std::endl;
+			file.write(reinterpret_cast<const char*>(&magicString.c_str()[0]), 14);	// write the file id
+			file.write(reinterpret_cast<const char*>(&desc.c_str()[0]), 58);		// write the description
+			file.write(reinterpret_cast<const char*>(&hNumVertices), sizeof(int));	// write #vert.
+			file.write(reinterpret_cast<const char*>(&hNumEdges), sizeof(int));		// write #edges
+			file.close();
+		}
+
+		// output the network as a string
+		std::string str() {
+			std::stringstream ss;
+			size_t nv = V.size();
+			size_t ne = E.size();
+			ss << "Node (" << nv << ")--------" << std::endl;
+			for (size_t i = 0; i < nv; i++)
+				ss << "\t" << i << V[i].str() << std::endl;
+			ss << "Edge (" << ne << ")--------" << std::endl;
+			for (size_t i = 0; i < ne; i++)
+				ss << "\t" << i << E[i].str() << std::endl;
+		
+			return ss.str();
+		}
+
+		// get a string of edges
+		std::string strTxt(std::vector<std::vector<T> > p) {
+			std::stringstream ss;
+			std::stringstream oss;
+			size_t num = p.size();
+			for (size_t i = 0; i < p; i++) {
+				ss.str(std::string());
+				for (size_t j = 0; j < 3; j++)
+					ss << p[i][j];
+				ss << "\n";
+			}
+
+			return ss.str();
+		}
+
+		// removes specified character from string
+		void removeCharsFromString(std::string &str, char* charsToRemove) {
+			for (size_t i = 0; i < strlen(charsToRemove); i++)
+				str.erase((remove(str.begin(), str.end(), charsToRemove[i])), str.end());
+		}
+
+		// exports network to txt file
+		void to_txt(std::string filename) {
+			std::ofstream ofs(filename.c_str(), std::ofstream::out | std::ofstream::app);
+
+			ofs << (E.size()).str() << "\n";
+			for (size_t i = 0; i < E.size(); i++) {
+				std::string str;
+				ofs << (E[i].size()).str() << "\n";
+				str = E[i].strTxt();
+				ofs << str << "\n";
+			}
+			for (size_t i = 0; i < V.size(); i++) {
+				std::string str;
+				str = V[i].str();
+				char temp[4] = "[],";
+				removeCharsFromString(str, temp);
+				ofs << str << "\n";
+			}
+			ofs.close();
+		}
+
+
+		/// advanced operations
+		// adding a fiber to current network
+		// prior information: attaching fiber "f" and point "p" information, if "order" = 0 means the first point on fiber "e" is the attaching one (others mean the last point)
+		// "add" = 0 means replace the point on fiber "e" which is to be attached to the point on current fiber, "add" = 1 means add that one. Default "add" = 0
+		// ********** we don't accept that one fiber has only 3 points on it, reorder if this happens **********
+		void add_fiber(edge e, size_t f, size_t p, size_t order, size_t add = 0) {
+			size_t num = E[f].size();		// get the number of points on this fiber
+			size_t num1 = p + 1;			// first "half"
+			size_t num2 = num - p;			// second "half"
+			size_t id = p;					// split point on the fiber that attach to
+			size_t ne = e.size();
+
+			// if a new fiber only has points that less than 4, either add or change (default) an attaching point
+			if (num1 < 4)
+				id = 0;
+			else if (num2 < 4)
+				id = num - 1;
+			
+			// check to see whether it needs to replace/add the joint point
+			if (add == 0) {				// change the point
+				if (order == 0) {
+					e[0] = E[f][id];
+					e.set_r(0, r(f, id));
+				}
+				else {
+					e[ne - 1] = E[f][id];
+					e.set_r(ne - 1, r(f, id));
+				}
+			}
+			else {						// add a point if necessary
+				if (order == 0) {
+					if (e[0] == E[f][id])	// if they share the same joint point, make sure radii are consistent
+						e.set_r(0, r(f, id));
+					else
+						e.insert(0, E[f][id], r(f, id));
+				}
+				else {
+					if (e[ne - 1] == E[f][id])
+						e.set_r(ne - 1, r(f, id));
+					else
+						e.insert(ne - 1, E[f][id], r(f, id));
+				}
+			}
+
+			std::vector<edge> tmp_edge = E[f].split(id);			// check and split
+			// current one hasn't been splitted
+			if (tmp_edge.size() == 1) {
+				if (id == 0) {		// stitch location is the starting point of current edge
+					if (V[E[f].v[0]].e[0].size() + V[E[f].v[0]].e[1].size() > 1) {	// branching point
+						if (order == 0) {
+							V[E[f].v[0]].e[0].push_back(E.size());
+							edge new_edge(e);
+							new_edge.v[0] = E[f].v[0];	// set the starting and ending points for the new edge
+							new_edge.v[1] = V.size();
+							E.push_back(new_edge);
+							vertex new_vertex = e[e.size() - 1];
+							new_vertex.e[1].push_back(V.size());	// set the incoming edge for the new point
+							V.push_back(new_vertex);
+						}
+						else {
+							V[E[f].v[0]].e[1].push_back(E.size());
+							edge new_edge(e);
+							new_edge.v[0] = V.size();	// set the starting and ending points for the new edge
+							new_edge.v[1] = E[f].v[0];
+							E.push_back(new_edge);
+							vertex new_vertex = e[e.size() - 1];
+							new_vertex.e[0].push_back(V.size());	// set the outgoing edge for the new point
+							V.push_back(new_vertex);
+						}
+					}
+					else {								// not branching point
+						size_t k = E[f].v[0];			// get the index of the starting point on current edge
+						vertex new_vertex;
+						edge new_edge;
+						if (order == 0) {
+							new_vertex = e[e.size() - 1];
+							new_edge = E[f].concatenate(e, 0, 0);
+						}
+						else {
+							new_vertex = e[0];
+							new_edge = E[f].concatenate(e, 0, e.size() - 1);
+						}
+						new_vertex.e[0].push_back(f);
+						new_edge.v[1] = E[f].v[1];		// set starting and ending points for the new concatenated edge
+						new_edge.v[0] = E[f].v[0];
+						V[k] = new_vertex;
+						E[f] = new_edge;
+					}
+				}
+				else {			// stitch location is the ending point of current edge
+					if (V[E[f].v[1]].e[0].size() + V[E[f].v[1]].e[1].size() > 1) {	// branching point
+						if (order == 0) {
+							V[E[f].v[1]].e[0].push_back(E.size());
+							edge new_edge(e);
+							new_edge.v[0] = E[f].v[1];	// set the starting and ending points for the new edge
+							new_edge.v[1] = V.size();
+							E.push_back(new_edge);
+							vertex new_vertex = e[e.size() - 1];
+							new_vertex.e[1].push_back(V.size());	// set the incoming edge for the new point
+							V.push_back(new_vertex);
+						}
+						else {
+							V[E[f].v[1]].e[1].push_back(E.size());
+							edge new_edge(e);
+							new_edge.v[0] = V.size();	// set the starting and ending points for the new edge
+							new_edge.v[1] = E[f].v[1];
+							E.push_back(new_edge);
+							vertex new_vertex = e[e.size() - 1];
+							new_vertex.e[0].push_back(V.size());	// set the outgoing edge for the new point
+							V.push_back(new_vertex);
+						}
+					}
+					else {								// not branching point
+						size_t k = E[f].v[1];			// get the index of the ending point on current edge
+						vertex new_vertex;
+						edge new_edge;
+						if (order == 0) {
+							new_vertex = e[e.size() - 1];
+							new_edge = E[f].concatenate(e, num - 1, 0);
+						}
+						else {
+							new_vertex = e[0];
+							new_edge = E[f].concatenate(e, num - 1, e.size() - 1);
+						}
+						new_vertex.e[1].push_back(f);
+						new_edge.v[1] = E[f].v[1];	// set starting and ending points for the new concatenated edge
+						new_edge.v[0] = E[f].v[0];
+						V[k] = new_vertex;
+						E[f] = new_edge;
+					}
+				}
+			}
+			// current one has been splitted
+			else {
+				vertex new_vertex = E[f][id];
+				V.push_back(new_vertex);
+				tmp_edge[0].v[0] = E[f].v[0];
+				tmp_edge[0].v[1] = V.size() - 1;		// set the ending point of the first half edge
+				tmp_edge[1].v[0] = V.size() - 1;		// set the starting point of the second half edge
+				tmp_edge[1].v[1] = E[f].v[1];
+				edge tmp(E[f]);
+				E[f] = tmp_edge[0];						// replace current edge by the first half edge
+				E.push_back(tmp_edge[1]);
+				V[V.size() - 1].e[0].push_back(E.size() - 1);			// set the incoming and outgoing edges for the splitted point
+				V[V.size() - 1].e[1].push_back(f);						// push "f" fiber as an incoming edge for the splitted point
+				for (size_t i = 0; i < V[tmp.v[1]].e[1].size(); i++)	// set the incoming edge for the original ending vertex
+					if (V[tmp.v[1]].e[1][i] == f)
+						V[tmp.v[1]].e[1][i] = E.size() - 1;
+
+				if (order == 0) {
+					e.v[0] = V.size() - 1;				// set the starting and ending points for the new edge
+					e.v[1] = V.size();
+					V[V.size() - 1].e[0].push_back(E.size());	// we assume "flow" flows from starting point to ending point!
+					new_vertex = e[e.size() - 1];		// get the ending point on the new edge
+					E.push_back(e);
+					V.push_back(new_vertex);
+					V[V.size() - 1].e[1].push_back(E.size() - 1);
+				}
+				else {
+					e.v[0] = V.size();					// set the starting and ending points for the new edge
+					e.v[1] = V.size() - 1;
+					V[V.size() - 1].e[1].push_back(E.size());
+					new_vertex = e[0];					// get the ending point on the new edge
+					E.push_back(e);
+					V.push_back(new_vertex);
+					V[V.size() - 1].e[0].push_back(E.size() - 1);
+				}
+			}
+		}
+
+		// THIS IS FOR PAVEL
+		// @param "e" is the edge that is to be stitched
+		// @param "f" is the index of edge that is to be stiched to
+		void stitch(edge e, size_t f) {
+			network<T> A = (*this);			// make a copy of current network
+
+			T* query_point = new T[3];
+			for (size_t k = 0; k < 3; k++)
+				query_point[k] = e[0][k];			// we assume the first one is the one to be stitched
+
+			size_t num = A.E[f].size();		// get the number of points on edge "f"
+			T* reference_point = (T*)malloc(sizeof(T) * num * 3);
+			A.E[f].edge_to_array(reference_point);
+			size_t max_tree_level = 3;
+
+			stim::kdtree<T, 3> kdt;					// initialize a tree
+			kdt.create(reference_point, num, max_tree_level);		// build a tree
+			
+			T* dist = new T[1];
+			size_t* nnIdx = new size_t[1];
+
+#ifdef __CUDACC__
+			kdt.search(query_point, 1, nnIdx, dist);	// search for nearest neighbor
+#else
+			kdt.cpu_search(query_point, 1, nnIdx, dist);// search for nearest neighbor
+#endif
+			add_fiber(e, f, nnIdx[0], 0, 1);
+
+			free(reference_point);
+			delete(dist);
+			delete(nnIdx);
+		}
+		
+		// split current network at "idx" location on edge "f"
+		network<T> split(size_t f, size_t idx) {
+			size_t num = E.size();
+
+			if (f >= num) {				// outside vector size
+			}
+			else {
+				if (idx <= 0 || idx >= num - 1) {	// can't split at this position
+				}
+				else {
+					std::vector<edge> list;			// a list of edges
+					list = E[f].split(idx);			// split in tems of edges
+					// first segment replaces original one
+					edge new_edge = list[0];		// new edge
+					edge tmp(E[f]);					// temp edge
+					new_edge.v[0] = E[f].v[0];		// copy starting point
+					new_edge.v[1] = V.size();		// set ending point
+					E[f] = new_edge;				// replacement
+					vertex new_vertex(new_edge[idx]);
+					new_vertex.e[1].push_back(f);			// incoming edge for the new vertex
+					new_vertex.e[0].push_back(E.size());	// outgoing edge for the new vertex
+
+					// second segment gets newly push back
+					new_edge = list[1];				// new edge
+					new_edge.v[1] = tmp.v[1];		// copy ending point
+					new_edge.v[0] = V.size();		// set starting point
+					size_t n = V[tmp.v[1]].e[1].size();
+					for (size_t i = 0; i < n; i++) {
+						if (V[tmp.v[1]].e[1][i] == f) {
+							V[tmp.v[1]].e[1][i] = E.size();
+							break;
+						}
+					}
+					
+					V.push_back(new_vertex);
+					E.push_back(new_edge);
+				}
+			}
+			
+			return (*this);
+		}
+
+		// resample current network
+		network<T> resample(T spacing) {
+			stim::network<T> result;
+
+			result.V = V;	// copy vertices
+			size_t num = E.size();	// get the number of edges
+			result.E.resize(num);
+
+			for (size_t i = 0; i < num; i++)
+				result.E[i] = E[i].resample(spacing);
+
+			return result;
+		}
+
+		// copy the point cload representing the centerline for the network into an array
+		void centerline_cloud(T* dst) {
+			size_t p;				// store the current edge point
+			size_t P;				// store the number of points in an edge
+			size_t t = 0;			// index in to the output array of points
+			size_t num = E.size();
+			for (size_t i = 0; i < num; i++) {
+				P = E[i].size();
+				for (p = 0; p < P; p++) {
+					dst[t * 3 + 0] = E[i][p][0];
+					dst[t * 3 + 1] = E[i][p][1];
+					dst[t * 3 + 2] = E[i][p][2];
+					t++;
+				}
+			}
+		}
+
+		// subdivide current network and represent as a explicit undirected graph
+		network<T> to_graph() {
+			std::vector<size_t> OI;			// a list of original vertex index
+			std::vector<size_t> NI;			// a list of new vertex index
+			std::vector<edge> nE;			// a list of new edge
+			std::vector<vertex> nV;			// a list of new vector
+			size_t id = 0;					// temp vertex index
+			size_t n = E.size();
+
+			for (size_t i = 0; i < n; i++) {
+				if (E[i].size() == 2) {		// if current edge only contain 2 points, can't be subdivided
+					stim::centerline<T> line;
+					for (size_t k = 0; k < 2; k++)		// copy points on current network
+						line.push_back(E[i][k]);
+
+					edge new_edge(line);	// construct edge based on current centerline
+
+					for (size_t k = 0; k < 2; k++) {
+						vertex new_vertex = new_edge[k];
+						id = E[i].v[k];			// get the starting/ending point index
+						std::vector<size_t>::iterator pos = std::find(OI.begin(), OI.end(), id);		// search this index through the list of new vertex index and see whether it appears
+						if (pos == OI.end()) {			// current vertex hasn't been pushed back
+							OI.push_back(id);			// add current vertex to the original list
+							NI.push_back(nV.size());	// add new index
+
+							new_vertex.e[k].push_back(nE.size());	// set outgoing/incoming edge for the new vertex
+							new_edge.v[k] = nV.size();				// set the starting/ending vertex for the new edge
+							nV.push_back(new_vertex);
+						}
+						else {							// current vertex has been pushed back before
+							int d = std::distance(OI.begin(), pos);
+							new_edge.v[k] = NI[d];
+							nV[NI[d]].e[k].push_back(nE.size());
+						}
+					}
+
+					nE.push_back(new_edge);
+
+					// copy radii information
+					nE[nE.size() - 1].set_r(0, E[i].r(0));
+					nE[nE.size() - 1].set_r(1, E[i].r(1));
+				}
+				else {						// if current edge contain at least 3 points, can be subdivided
+					size_t num = E[i].size();
+					for (size_t j = 0; j < num - 1; j++) {
+						stim::centerline<T> line;
+						for (size_t k = 0; k < 2; k++)		// copy points on current network
+							line.push_back(E[i][k]);
+
+						edge new_edge(line);
+						if (j == 0) {					// for edge that contains original starting point
+							vertex new_vertex = new_edge[0];
+							id = E[i].v[0];
+							std::vector<size_t>::iterator pos = std::find(OI.begin(), OI.end(), id);
+							if (pos == OI.end()) {		// new vertex
+								OI.push_back(id);
+								NI.push_back(nV.size());
+
+								new_vertex.e[0].push_back(nE.size());
+								new_edge.v[0] = nV.size();
+								nV.push_back(new_vertex);
+							}
+							else {
+								int d = std::distance(OI.begin(), pos);
+								new_edge.v[0] = NI[d];
+								nV[NI[d]].e[0].push_back(nE.size());
+							}
+
+							new_vertex = new_edge[1];
+							new_vertex.e[1].push_back(nE.size());
+							new_edge.v[1] = nV.size();
+							nV.push_back(new_vertex);
+							nE.push_back(new_edge);
+						}
+
+						else if (j == E[i].size() - 2) {// for edge that contains original ending point
+							vertex new_vertex = new_edge[1];
+							nV[nV.size() - 1].e[0].push_back(nE.size());
+							new_edge.v[0] = nV.size() - 1;
+
+							id = E[i].v[1];				// get ending vertex index
+							std::vector<size_t>::iterator pos = std::find(OI.begin(), OI.end(), id);
+							if (pos == OI.end()) {
+								OI.push_back(id);
+								NI.push_back(nV.size());
+
+								new_vertex.e[1].push_back(nE.size());
+								new_edge.v[1] = nV.size();
+								nV.push_back(new_vertex);
+							}
+							else {
+								int d = std::distance(OI.begin(), pos);
+								new_edge.v[1] = NI[d];
+								nV[NI[d]].e[1].push_back(nE.size());
+							}
+
+							nE.push_back(new_edge);
+						}
+
+						else {							// for edge that doesn't contain original ending point
+							vertex new_vertex = new_edge[1];
+
+							// push back the latter one on that segment
+							nV[nV.size() - 1].e[0].push_back(nE.size());
+							new_vertex.e[1].push_back(nE.size());
+							new_edge.v[0] = nV.size() - 1;
+							new_edge.v[1] = nV.size();
+							nV.push_back(new_vertex);
+							nE.push_back(new_edge);
+						}
+
+						// copy radii information
+						nE[nE.size() - 1].set_r(0, E[i].r(j));
+						nE[nE.size() - 1].set_r(1, E[i].r(j + 1));
+					}
+				}
+			}
+
+			stim::network<T> result(nE, nV);
+
+			return result;
+		}
+
+		// this function compares two networks and returns the percentage of the current network that is missing from "A"
+		// @param "A" is the network to compare to - the field is generated for A
+		// @param "sigma" is the user-defined tolerance value - smaller values provide a stricter comparison
+		// @param "device" is the GPU device to use - default no GPU provided
+		network<T> compare(network<T> A, T sigma, int device = -1) {
+			network<T> R;		// generate a network storing the result of comparison
+			R = (*this);		// initialize to current network
+
+			size_t num = A.total_points();
+			T* c = (T*)malloc(sizeof(T) * num * 3);			// allocate memory for the centerline of A
+
+			A.copy_to_array(c);				// copy points in A to a 1D array
+
+			size_t max_tree_level = 3;		// set max tree level parameter to 3
+
+#ifdef __CUDACC__
+			cudaSetDevice(device);
+			stim::kdtree<T, 3> kdt;					// initialize a tree object
+
+			kdt.create(c, num, max_tree_level);		// build tree
+
+			for (size_t i = 0; i < R.E.size(); i++) {
+				size_t n = R.E[i].size();			// the number of points in current edge
+				T* query_point = new T[3 * n];		// allocate memory for points
+				T* m1 = new T[n];					// allocate memory for metrics
+				T* dists = new T[n];				// allocate memory for distances
+				size_t* nnIdx = new size_t[n];		// allocate memory for indices
+
+				T* d_dists;
+				T* d_m1;
+				cudaMalloc((void**)&d_dists, n * sizeof(T));
+				cudaMalloc((void**)&d_m1, n * sizeof(T));
+
+				edge_to_array(query_point, R.E[i]);
+				kdt.search(query_point, n, nnIdx, dists);
+
+				cudaMemcpy(d_dists, dists, n * sizeof(T), cudaMemcpyHostToDevice);			// copy dists from host to device
+
+				// configuration parameters
+				size_t threads = (1024 > n) ? n : 1024;
+				size_t blocks = n / threads + (n % threads) ? 1 : 0;
+
+				find_metric<< <blocks, threads >> >(d_m1, n, d_dists, sigma);		// calculate the metric value based on the distance
+
+				cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost);
+
+				for (size_t j  = 0; j < n; j++) {
+					R.E[i].set_r(j, m1[j]);
+				}
+			}
+
+#else		// if there is any GPU device, use CPU - much slower
+			stim::kdtree<T, 3> kdt;
+			kdt.create(c, num, max_tree_level);
+
+			for (size_t i = 0; i < R.E.size(); i++) {			// for each edge in A
+
+				size_t n = R.E[i].size();						// the number of points in current edge
+				T* query = new T[3 * n];
+				T* m1 = new T[n];
+				T* dists = new T[n];
+				size_t* nnIdx = new size_t[n];
+
+				edge_to_array(query, R.E[i]);
+
+				kdt.cpu_search(query, n, nnIdx, dists);			// find the distance between A and the current network
+
+				for (size_t j = 0; j < n; j++) {
+					m1[j] = 1.0f - gaussian(dists[j], sigma);	// calculate the metric value based on the distance
+					R.E[i].set_r(j, m1[j]);						// set the error for the second point in the segment
+				}
+			}
+#endif
+			return R;
+		}
+
+		// this function compares two splitted networks to yield a mapping relationship between them according to nearest neighbor principle
+		// @param "B" is the template network
+		// @param "C" is the mapping relationship: C[e1] = _e1 means edge "e1" in current network is mapped to edge "_e1" in "B"
+		// @param "device" is the GPU device that user want to use
+		// @param "threshold" is to control mapping tolerance (threshold)
+		void mapping(network<T> B, std::vector<size_t> &C, T threshold, int device = -1) {
+			network<T> A;			// generate a network storing the result of the comparison
+			A = (*this);
+
+			size_t nA = A.E.size();	// the number of edges in A
+			size_t nB = B.E.size();	// the number of edges in B
+
+			C.resize(A.E.size());
+
+			size_t num = B.total_points();			// set the number of points
+			T* c = (T*)malloc(sizeof(T) * num * 3);
+			
+			B.copy_to_array(c);
+
+			size_t max_tree_level = 3;
+
+#ifdef __CUDACC__
+			cudaSetDevice(device);
+			stim::kdtree<T, 3> kdt;				// initialize a tree
+			
+			kdt.create(c, num, max_tree_level);		// build a tree
+
+			T M = 0.0f;
+			for (size_t i = 0; i < nA; i++) {
+				M = A.ar(i);			// get the average metric of edge "i"
+				if (M > threshold)
+					C[i] = UINT_MAX;	// set to MAX
+				else {
+					T* query_point = new T[3];
+					T* dist = new T[1];
+					size_t* nnIdx = new size_t[1];
+
+					for (size_t k = 0; k < 3; k++)
+						query_point[k] = A.E[i][A.E[i].size() / 2][k];		// search by the middle one, risky?
+					kdt.search(query_point, 1, nnIdx, dist);
+
+					size_t id = 0;
+					size_t sum = 0;
+					for (size_t j = 0; j < nB; j++) {		// low_band
+						sum += B.E[j].size();
+						if (nnIdx[0] < sum) {
+							C[i] = id;
+							break;
+						}
+						id++;
+					}
+				}
+			}
+
+#else
+			stim::kdtree<T, 3> kdt;
+			kdt.create(c, num, max_tree_level);
+			T* dist = new T[1];
+			size_t* nnIdx = new size_t[1];
+
+			stim::vec3<T> p;
+			T* query_point = new T[3];
+
+			T M = 0.0f;
+			for (size_t i = 0; i < nA; i++) {
+				M = A.ar(i);
+				if (M > threshold)
+					C[i] = UINT_MAX;
+				else {
+					p = A.E[i][A.E[i].size() / 2];
+					for (size_t k = 0; k < 3; k++)
+						query_point[k] = p[k];
+					kdt.cpu_search(query_point, 1, nnIdx, dist);
+
+					size_t id = 0;
+					size_t sum = 0;
+					for (size_t j = 0; j < nB; j++) {
+						sum += B.E[j].size();
+						if (nnIdx[0] < sum) {
+							C[i] = id;
+							break;
+						}
+						id++;
+					}
+				}
+			}
+#endif
+		}
+	};
+
+}
+
+
+#endif
@@ -15,14 +15,14 @@ find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y, i
 {
 	float		sigma		= 2.0;
 	unsigned int	conn		= 7;
-	float		threshold	= 40.0;
+//	float		threshold	= 40.0;
 	float*		cpuCenters	= (float*) malloc(x*y*sizeof(float));
 	int		sizek		= 7;
  
-	stringstream name;
+	std::stringstream name;
  
  
-	cpuCenters = stim::cuda::get_centers(texbufferID, texType, x, y, sizek, sigma, conn, threshold, iter);
+	cpuCenters = stim::cuda::get_centers(texbufferID, texType, x, y, sizek, sigma, conn, iter);
 	cudaDeviceSynchronize();
  
  
@@ -11,7 +11,7 @@
 #include <stim/cuda/cudatools/devices.h>
 #include <stim/cuda/cudatools/threads.h>
 #include <stim/cuda/cuda_texture.cuh>
-#include <stim/cuda/ivote.cuh>
+#include <stim/iVote/ivote2.cuh>
 #include <stim/cuda/arraymath.cuh>
  
 #define IMAD(a,b,c) ( __mul24((a), (b)) + (c) )
@@ -156,7 +156,7 @@ namespace stim
  
 	extern "C"
 	float *
-	get_centers(GLint texbufferID, GLenum texType, int DIM_X, int DIM_Y, int sizeK, float sigma, float conn, float threshold, int iter = 0)
+	get_centers(GLint texbufferID, GLenum texType, int DIM_X, int DIM_Y, int sizeK, float sigma, float conn, int iter = 0)
 	{
 		tx.SetTextureCoordinates(1);
 		tx.SetAddressMode(1, 3);
@@ -180,7 +180,7 @@ namespace stim
 		#endif
  
  
-		stim::cuda::gpu_local_max<float>(centers, res, threshold, conn, DIM_X, DIM_Y);
+		stim::cuda::gpu_local_max<float>(centers, res, conn, DIM_X, DIM_Y);
  
 		cudaDeviceSynchronize();
  
@@ -2,33 +2,46 @@
 #define STIM_GL_SPIDER_H
  
 //#include <GL/glew.h>
+///basic GL and Cuda libs
 #include <GL/glut.h>
 #include <cuda.h>
 #include <cuda_gl_interop.h>
 #include <cudaGL.h>
-#include <math.h>
+
+///stim .*h for gl tracking and visualization
 #include <stim/gl/gl_texture.h>
-#include <stim/visualization/camera.h>
 #include <stim/gl/error.h>
+#include <stim/visualization/camera.h>
+#include <stim/math/matrix_sq.h>
+#include <stim/cuda/spider_cost.cuh>
+#include <stim/visualization/glObj.h>
+#include <stim/cuda/cudatools/glbind.h>
+
+///stim *.h for CUDA
+#include <stim/cuda/cuda_texture.cuh>
+#include <stim/cuda/cudatools.h>
+
+///stim *.h for branch detection
+#include <stim/visualization/cylinder.h>
+#include <stim/cuda/branch_detection.cuh>
+#include <stim/visualization/gl_network.h>
+
+/// *.h for math
+#include <math.h>
 #include <stim/math/vector.h>
 #include <stim/math/vec3.h>
 #include <stim/math/rect.h>
-#include <stim/math/matrix.h>
 #include <stim/math/constants.h>
-#include <stim/cuda/spider_cost.cuh>
-#include <stim/cuda/cudatools/glbind.h>
 #include <stim/cuda/arraymath.cuh>
-#include <stim/cuda/cuda_texture.cuh>
-#include <stim/cuda/cudatools.h>
-#include <stim/cuda/ivote.cuh>
-#include <stim/visualization/glObj.h>
+#include <stim/math/random.h>
+
+
+///c++ libs for everything
 #include <vector>
 #include <stack>
-#include <stim/cuda/branch_detection.cuh>
-#include "../../../volume-spider/glnetwork.h"
-#include <stim/visualization/cylinder.h>
 #include <iostream>
 #include <fstream>
+
 #ifdef TIMING
 	#include <ctime>
 	#include <cstdio>	
@@ -72,7 +85,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 		std::vector<float> mV;	//A list of all the size vectors.
 		std::vector<float> lV;	//A list of all the size vectors.
  
-		stim::matrix<float, 4> cT;			//current Transformation matrix (tissue)->(texture)
+		stim::matrix_sq<float, 4> cT;			//current Transformation matrix (tissue)->(texture)
 		GLuint texID;						//OpenGL ID for the texture to be traced
 		stim::vec3<float> S;					//Size of a voxel in the volume.
 		stim::vec3<float> R;					//Dimensions of the volume.
@@ -121,7 +134,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 		std::vector< float > cM;		//radius up to the current point
 		std::vector< float > cLen;		//radius up to the current point
  
-		stim::glnetwork<float> nt;				//network object holding the currently traced centerlines
+		stim::gl_network<float> nt;				//network object holding the currently traced centerlines
 		stim::glObj<float> sk;					//OBJ file storing the network (identical to above)
  
 		//consider replacing with two seed points facing opposite directions
@@ -1318,25 +1331,25 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 			}
 		}
  
-		///Saves the network to a file.
-		void
-		saveNetwork(std::string name)
-		{
-			stim::glObj<float> sk1;
-			for(int i = 0; i < nt.sizeE(); i++)
-			{
-				std::vector<float> cm = nt.getEdgeCenterLineMag(i);
-                 		std::vector<stim::vec3< float > > ce = nt.getEdgeCenterLine(i);
-				sk1.Begin(stim::OBJ_LINE);
-				for(int j = 0; j < ce.size(); j++)
-				{
-					sk1.TexCoord(cm[j]);
-					sk1.Vertex(ce[j][0], ce[j][1], ce[j][2]);
-				}
-				sk1.End();
-			}	
-			sk1.save(name);
-		}
+                ///Saves the network to a file.
+                void
+                saveNetwork(std::string name, int xoffset = 0, int yoffset = 0, int zoffset = 0)
+                {
+                        stim::glObj<float> sk1; 
+                        for(int i = 0; i < nt.sizeE(); i++) 
+                        {
+                                std::vector<float> cm = nt.getEdgeCenterLineMag(i);
+                                std::vector<stim::vec3< float > > ce = nt.getEdgeCenterLine(i);
+                                sk1.Begin(stim::OBJ_LINE);
+                                for(int j = 0; j < ce.size(); j++) 
+                                {
+                                        sk1.TexCoord(cm[j]);
+                                        sk1.Vertex(ce[j][0]+xoffset, ce[j][1]+yoffset, ce[j][2]+zoffset);
+                                }
+                                sk1.End();
+                        }     
+                        sk1.save(name);
+                }
  
 		///Depreciated, but might be reused later()
 		///returns a COPY of the entire stim::glObj object.
@@ -1347,7 +1360,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 		}
  
 		///returns a COPY of the entire stim::glnetwork object.
-		stim::glnetwork<T>
+		stim::gl_network<T>
 		getGLNetwork()
 		{
 			return nt;
@@ -152,13 +152,17 @@ namespace stim {
  
  
 	template<typename T>
-	void cpu_ivote2(T* cpuI, unsigned int rmax, size_t x, size_t y, bool invert = false, T t = 0, int iter = 8, T phi = 15.0f * (float)stim::PI / 180, int conn = 8, bool debug = false) {
+	void cpu_ivote2(T* cpuI, unsigned int rmax, size_t x, size_t y, float &gpu_time, bool invert = false, T t = 0, int iter = 8, T phi = 15.0f * (float)stim::PI / 180, int conn = 8, bool debug = false) {
 		size_t bytes = x*y * sizeof(T);
 		T* gpuI;						//allocate space on the gpu to save the input image
+
+		gpuTimer_start();
 		HANDLE_ERROR(cudaMalloc(&gpuI, bytes));
 		HANDLE_ERROR(cudaMemcpy(gpuI, cpuI, bytes, cudaMemcpyHostToDevice));		//copy the image to the gpu
 		stim::gpu_ivote2<T>(gpuI, rmax, x, y, invert, t, iter, phi, conn, debug);				//call the gpu version of the ivote
 		HANDLE_ERROR(cudaMemcpy(cpuI, gpuI, bytes, cudaMemcpyDeviceToHost));		//copy the output to the cpu
+
+		gpu_time = gpuTimer_end();
 	}
 }
 #endif
 \ No newline at end of file
@@ -10,7 +10,7 @@ namespace stim{
  
 		// this kernel calculates the local maximum for finding the cell centers
 		template<typename T>
-		__global__ void cuda_local_max(T* gpuCenters, T* gpuVote, int conn, int x, int y){
+		__global__ void cuda_local_max(T* gpuCenters, T* gpuVote, int conn, int x,  int y){
  
 			// calculate the 2D coordinates for this current thread.
 			int xi = blockIdx.x * blockDim.x + threadIdx.x;
@@ -20,16 +20,16 @@ namespace stim{
 				return;
  
 			// convert 2D coordinates to 1D
-			int i = yi * x + xi;
+		    int i = yi * x + xi;
  
 			gpuCenters[i] = 0;		//initialize the value at this location to zero
  
 			T val = gpuVote[i];
  
-			for(int xl = xi - conn; xl < xi + conn; xl++){
-				for(int yl = yi - conn; yl < yi + conn; yl++){
+			for(unsigned int xl = xi - conn; xl < xi + conn; xl++){
+				for(unsigned int yl = yi - conn; yl < yi + conn; yl++){
 					if(xl >= 0 && xl < x && yl >= 0 && yl < y){
-						int il = yl * x + xl;
+						unsigned int il = yl * x + xl;
 						if(gpuVote[il] > val){							
 							return;
 							}
@@ -47,7 +47,7 @@ namespace stim{
 		}
  
 		template<typename T>
-		void gpu_local_max(T* gpuCenters, T* gpuVote,  unsigned int conn, size_t x, size_t y){
+		void gpu_local_max(T* gpuCenters, T* gpuVote,  unsigned int conn, unsigned int x, unsigned int y){
  
 			unsigned int max_threads = stim::maxThreadsPerBlock();
 			/*dim3 threads(max_threads, 1);
@@ -275,10 +275,10 @@ public:
 		int channels = cvImage.channels();
 		allocate(cols, rows, channels);			//allocate space for the image
 		size_t img_bytes = bytes();
-		unsigned char* cv_ptr = (unsigned char*)cvImage.data;
-		//if (C() == 1)														//if this is a single-color image, just copy the data
-		//	type_copy<T, T>(cv_ptr, img, size());
-		memcpy(img, cv_ptr, bytes());
+		unsigned char* cv_ptr = (unsigned char*)cvImage.data;		
+		if (C() == 1)														//if this is a single-color image, just copy the data
+			type_copy<unsigned char, T>(cv_ptr, img, size());
+		//memcpy(img, cv_ptr, bytes());
 		if(C() == 3)														//if this is a 3-color image, OpenCV uses BGR interleaving
 			from_opencv(cv_ptr, X(), Y());
 #else
+#ifndef STIM_CUDA_MEDIAN2_H
+#define STIM_CUDA_MEDIAN2_H
+
+#include <iostream>
+#include <cuda.h>
+#include <cmath>
+#include <algorithm>
+
+#ifdef USING_CUDA
+#include "cuda_runtime.h"
+#include "device_launch_parameters.h"
+#include <stim/cuda/cudatools.h>
+#endif
+
+
+//#include <thrust/sort.h>
+//#include <thrust/execution_policy.h>
+//#include <thrust/binary_search.h>
+//#include <thrust/device_ptr.h>
+
+template <typename T> 
+__device__ void cuswap ( T& a, T& b ){
+	T c(a);
+	a=b;
+	b=c;
+}
+
+namespace stim{
+	namespace cuda{
+
+		template<typename T>
+		__global__ void cuda_median2(T* in, T* out, T* kernel, size_t X, size_t Y, size_t K){
+
+			int xi = blockIdx.x * blockDim.x + threadIdx.x;
+			int yi = blockIdx.y * blockDim.y + threadIdx.y;			
+
+			size_t Xout = X - K + 1;
+			size_t Yout = Y - K + 1;
+			int i = yi * Xout + xi;
+
+			//return if (i,j) is outside the matrix
+			if(xi >= Xout || yi >= Yout)   return;		
+
+			//scan the image, copy data from input to 2D window
+			size_t kxi, kyi;
+			for(kyi = 0; kyi < K; kyi++){
+				for(kxi = 0; kxi < K; kxi++){
+					kernel[i * K * K + kyi * K + kxi] = in[(yi + kyi)* X + xi + kxi];				      
+				}
+			}					
+					
+			//calculate kernel radius 4
+            int r = (K*K)/2;
+			
+			//sort the smallest half pixel values inside the window, calculate the middle one
+			size_t Ki = i * K * K;
+		
+			//sort the smallest half pixel values inside the window, calculate the middle one
+            for (int  p = 0; p < r+1; p++){                             
+			    for(int kk = p + 1; kk < K*K; kk++){
+                       if (kernel[Ki + kk] < kernel[Ki + p]){
+						   cuswap<T>(kernel[Ki + kk], kernel[Ki + p]);
+					   }
+                  }
+             }
+            
+			//copy the middle pixel value inside the window to output
+            out[i] = kernel[Ki + r]; 
+				     
+        }
+		 
+		
+		template<typename T>
+		void gpu_median2(T* gpu_in, T* gpu_out, T* gpu_kernel, size_t X, size_t Y, size_t K){
+               
+			//get the maximum number of threads per block for the CUDA device
+            int threads_total = stim::maxThreadsPerBlock();
+
+			//set threads in each block
+			dim3 threads(sqrt(threads_total), sqrt(threads_total));
+
+			//calculate the number of blocks
+			dim3 blocks(( (X - K + 1) / threads.x) + 1, ((Y - K + 1) / threads.y) + 1);
+
+			//call the kernel to perform median filter function
+		    cuda_median2 <<< blocks, threads >>>( gpu_in, gpu_out, gpu_kernel, X, Y, K);
+		
+		}
+
+		template<typename T>
+		void cpu_median2(T* cpu_in, T* cpu_out, size_t X, size_t Y, size_t K){
+
+		#ifndef USING_CUDA
+			
+			//output image width and height
+			 size_t X_out = X + K -1;
+			 size_t Y_out = Y + K -1;
+
+			 float* window = (float*)malloc(K * k *sizeof(float));
+
+			 for(int i = 0; i< Y; i++){
+ 		       for(int j =0; j< X; j++){
+
+				//Pick up window elements		
+				int k = 0; 			
+				for (int m=0; m< kernel_width; m++){
+					for(int n = 0; n < kernel_width; n++){
+						window[k] = in_image.at<float>(m+i, n+j);
+						k++;					
+					 }
+				}
+
+				//calculate kernel radius 4
+				r_ker = K * K/2;
+			    //Order elements (only half of them)
+				for(int p = 0; p<r_ker+1; p++){			
+				  //Find position of minimum element			       
+			       for (int l = p + 1; l < size_kernel; l++){
+     
+					 if(window[l] < window[p]){
+					     float t  = window[p];
+					    window[p] = window[l];
+					    window[l] = t;          }		
+		     	}
+					}
+					
+			  //Get result - the middle element	
+				cpu_out[i * X_out + j] =  window[r_ker];	
+			   }
+			 }
+		#else
+			//calculate input and out image pixels & calculate kernel size
+			size_t N_in = X * Y;									//number of pixels in the input image
+			size_t N_out = (X - K + 1) * (Y - K + 1);								//number of pixels in the output image
+            //size_t kernel_size = kernel_width*kernel_width;				//total number of pixels in the kernel
+            
+			//allocate memory on the GPU for the array
+			T* gpu_in;																	//allocate device memory for the input image
+			HANDLE_ERROR( cudaMalloc( &gpu_in, N_in * sizeof(T) ) );
+			T* gpu_kernel; 															//allocate device memory for the kernel
+			HANDLE_ERROR( cudaMalloc( &gpu_kernel, K * K * N_out * sizeof(T) ) );
+			T* gpu_out;																	//allocate device memory for the output image
+			HANDLE_ERROR( cudaMalloc( &gpu_out, N_out * sizeof(T) ) );
+
+			//copy the array to the GPU
+			HANDLE_ERROR( cudaMemcpy( gpu_in, cpu_in, N_in * sizeof(T), cudaMemcpyHostToDevice) );
+			
+			//call the GPU version of this function
+			gpu_median2<T>(gpu_in, gpu_out, gpu_kernel, X, Y, K);
+
+			//copy the array back to the CPU
+			HANDLE_ERROR( cudaMemcpy( cpu_out, gpu_out, N_out * sizeof(T), cudaMemcpyDeviceToHost) );
+
+			//free allocated memory
+			cudaFree(gpu_in);
+			cudaFree(gpu_kernel);
+			cudaFree(gpu_out);
+
+		}
+		
+	}
+	#endif
+}
+
+
+#endif
@@ -40,7 +40,7 @@ namespace stim{
 	public:
 		/// Constructor opens a mat4 file for writing
 		mat4file(std::string filename) {
-			matfile.open(filename, std::ios::binary);
+			matfile.open(filename.c_str(), std::ios::binary);
 		}
  
 		bool is_open() {
@@ -74,6 +74,15 @@ struct matrix_sq
 					M[i*N+j] -= rhs;
 		return *this;
 	}
+	
+	// element wise divide
+	CUDA_CALLABLE matrix_sq<T, N> operator/(T rhs)
+	{
+		for(int i=0; i<N; i++)
+			for(int j=0; j<N; j++)
+				M[i*N+j] = M[i*N+j]/rhs;
+		return *this;
+	}
  
 	template<typename Y>
 	vec<Y> operator*(vec<Y> rhs){
@@ -116,6 +125,21 @@ struct matrix_sq
 		return ss.str();
 	}
  
+        std::string toTensor()
+        {
+                std::stringstream ss; 
+
+                ss << (*this)(0, 0) << " ";
+                ss << (*this)(0, 1) << " ";
+                ss << (*this)(1, 1) << " ";
+                ss << (*this)(2, 0) << " ";
+                ss << (*this)(2, 1) << " ";
+                ss << (*this)(2, 2); 
+
+                return ss.str();
+        }
+
+
 	static matrix_sq<T, N> identity() {
 		matrix_sq<T, N> I;
 		I = 0;
@@ -43,7 +43,7 @@ namespace stim {
 		}
  
 		void setc(unsigned int l, int m, T value) {
-			unsigned int c = coeff_1d(l, m);
+			u44nsigned int c = coeff_1d(l, m);
 			C[c] = value;
 		}
  
@@ -173,57 +173,57 @@ namespace stim {
 			project(sph_pts, l, m);
 		}
  
-		/// Generates a PDF describing the density distribution of points on a sphere
-		/// @param sph_pts is a list of points in cartesian coordinates 
-		/// @param l is the maximum degree of the spherical harmonic function
-		/// @param m is the maximum order
-		/// @param c is the centroid of the points in sph_pts. DEFAULT 0,0,0
-		/// @param n is the number of points of the surface of the sphere used to create the PDF. DEFAULT 1000
-		/// @param norm, a boolean that sets where the output vectors will be normalized between 0 and 1.
-		/// @param 
-		/*void pdf(std::vector<stim::vec3<T> > sph_pts, unsigned int l, int m, stim::vec3<T> c = stim::vec3<T>(0, 0, 0), unsigned int n = 1000, bool norm = true, std::vector<T> w = std::vector<T>())
-		{
-			std::vector<double> weights;		///the weight at each point on the surface of the sphere.
-												//		weights.resize(n);
-			unsigned int nP = sph_pts.size();
-			std::vector<stim::vec3<T> > sphere = stim::Random<T>::sample_sphere(n, 1.0, stim::TAU);
-			if (w.size() < nP)
-				w = std::vector<T>(nP, 1.0);
-
-			for (int i = 0; i < n; i++)
-			{
-				T val = 0;
-				for (int j = 0; j < nP; j++)
-				{
-					stim::vec3<T> temp = sph_pts[j] - c;
-					if (temp.dot(sphere[i]) > 0)
-						val += pow(temp.dot(sphere[i]), 4)*w[j];
-				}
-				weights.push_back(val);
-			}
-
-			mcBegin(l, m);		//begin spherical harmonic sampling
-
-			if (norm)
-			{
-				T min = *std::min_element(weights.begin(), weights.end());
-				T max = *std::max_element(weights.begin(), weights.end());
-				for (unsigned int i = 0; i < n; i++)
-				{
-					stim::vec3<T> sph = sphere[i].cart2sph();
-					mcSample(sph[1], sph[2], (weights[i] - min) / (max - min));
-				}
-
-			}
-			else {
-				for (unsigned int i = 0; i < n; i++)
-				{
-					stim::vec3<T> sph = sphere[i].cart2sph();
-					mcSample(sph[1], sph[2], weights[i]);
-				}
-			}
-			mcEnd();
-		}*/
+                /// Generates a PDF describing the density distribution of points on a sphere
+                /// @param sph_pts is a list of points in cartesian coordinates 
+                /// @param l is the maximum degree of the spherical harmonic function
+                /// @param m is the maximum order
+                /// @param c is the centroid of the points in sph_pts. DEFAULT 0,0,0
+                /// @param n is the number of points of the surface of the sphere used to create the PDF. DEFAULT 1000
+                /// @param norm, a boolean that sets where the output vectors will be normalized between 0 and 1.
+                /// @param 
+                void pdf(std::vector<stim::vec3<T> > sph_pts, unsigned int l, int m, stim::vec3<T> c = stim::vec3<T>(0, 0, 0),  unsigned int n = 1000, bool norm = true, std::vector<T> w = std::vector<T>(), int normfactor = 1)
+                {
+                        std::vector<double> weights;            ///the weight at each point on the surface of the sphere.
+                                                                                                //              weights.resize(n);
+                        unsigned int nP = sph_pts.size();       ///sph_pts is the vectors we want to fit a spehrical harmonic to.
+                        std::vector<stim::vec3<T> > sphere = stim::Random<T>::sample_sphere(n, 1.0, stim::TAU);         ///randomsample a sphere of radius 1 and return the sampled vectors.
+                        if (w.size() < nP)
+                                w = std::vector<T>(nP, 1.0);    //1 weight per each m.
+
+                        for (int i = 0; i < n; i++)             //for each weight
+                        {
+                                T val = 0;
+                                for (int j = 0; j < nP; j++)    //for each vector
+                                {
+                                        stim::vec3<T> temp = sph_pts[j] - c;
+                                        if (temp.dot(sphere[i]) > 0)
+                                                val += pow(temp.dot(sphere[i]), normfactor)*w[j];
+                                }
+                                weights.push_back(val);
+                        }
+
+                        mcBegin(l, m);          //begin spherical harmonic sampling
+
+                        if (norm)
+                        {
+                                T min = *std::min_element(weights.begin(), weights.end());
+                                T max = *std::max_element(weights.begin(), weights.end());
+                                for (unsigned int i = 0; i < n; i++)
+                                {
+                                        stim::vec3<T> sph = sphere[i].cart2sph();
+                                        mcSample(sph[1], sph[2], (weights[i] - min) / (max - min));
+                                }
+
+                        }
+                        else {
+                                for (unsigned int i = 0; i < n; i++)
+                                {
+                                        stim::vec3<T> sph = sphere[i].cart2sph();
+                                        mcSample(sph[1], sph[2], weights[i]);
+                                }
+                        }
+                        mcEnd();
+                }
  
 		std::string str() {
  
@@ -390,6 +390,29 @@ namespace stim {
 			}
 		}
  
+                ///imperpolates a spherical harmonic linearly from itself to the spherical harmonic passed as an input.
+                ///@param target, spharmonic class to which we are interpolating.
+                ///@param dist float value representing the distance we are interpolating [0,1]
+                stim::spharmonics<T>
+                interp(stim::spharmonics<T> target, float dist)
+                {
+                        stim::spharmonics<T> ret(std::max(target.C.size(), C.size()));
+                        for(int i = 0; i < target.C.size(); i++)
+                        {
+                                if(i < C.size())
+                                {
+                                        T slope = target.C[i] - C[i];
+                                        ret.C[i] = C[i] + dist * slope;
+                                }
+                                else
+                                {
+                                        T slope = target.C[i];
+                                        ret.C[i] = dist * slope;
+                                }
+                        }
+                        return ret;
+                }
+
 		/// Generate spherical harmonic coefficients based on a set of N samples
 		/*void fit(std::vector<stim::vec3<T> > sph_pts, unsigned int L, bool norm = true)
 		{
@@ -389,7 +389,7 @@ void cpu_scalar_mie_scatter(stim::complex&lt;T&gt;* E, size_t N, T* x, T* y, T* z, std
 		r = p.len();
 		if(r >= a){
 			for(size_t w = 0; w < W.size(); w++){
-				Ew = W[w].E() * exp(stim::complex<double>(0, W[w].kvec().dot(c)));
+				Ew = W[w].E() * exp(stim::complex<float>(0, W[w].kvec().dot(c)));
 				kr = p.len() * W[w].kmag();											//calculate k*r
 				stim::bessjyv_sph<double>(Nl, kr, vm, j_kr, y_kr, dj_kr, dy_kr);
 				cos_phi = p.norm().dot(W[w].kvec().norm());							//calculate the cosine of the angle from the propagating direction
@@ -649,7 +649,7 @@ void cpu_scalar_mie_internal(stim::complex&lt;T&gt;* E, size_t N, T* x, T* y, T* z, st
 		if(r < a){
 			E[i] = 0;
 			for(size_t w = 0; w < W.size(); w++){
-				Ew = W[w].E() * exp(stim::complex<double>(0, W[w].kvec().dot(c)));
+				Ew = W[w].E() * exp(stim::complex<float>(0, W[w].kvec().dot(c)));
 				knr = (stim::complex<double>)n * p.len() * W[w].kmag();							//calculate k*n*r
  
 				stim::cbessjyva_sph<double>(Nl, knr, vm, j_knr, y_knr, dj_knr, dy_knr);
-#ifndef STIM_GL_NETWORK
-#define STIM_GL_NETWORK
-
-#include <GL/glut.h>
-#include <stim/biomodels/network.h>
-#include <stim/visualization/aaboundingbox.h>
-
-namespace stim{
-
-template <typename T>
-class gl_network : public stim::network<T>{
-
-protected:
-	using stim::network<T>::E;
-	using stim::network<T>::V;
-
-	GLuint dlist;
-
-public:
-
-	/// Default constructor
-	gl_network() : stim::network<T>(){
-		dlist = 0;
-	}
-
-	/// Constructor creates a gl_network from a stim::network
-	gl_network(stim::network<T> N) : stim::network<T>(N){
-		dlist = 0;
-	}
-
-	/// Fills the parameters with the minimum and maximum spatial positions in the network,
-	///     specifying a bounding box for the network geometry
-	aaboundingbox<T> boundingbox(){
-
-		aaboundingbox<T> bb;								//create a bounding box
-
-		//loop through every edge
-		for(unsigned e = 0; e < E.size(); e++){
-			//loop through every point
-			for(unsigned p = 0; p < E[e].size(); p++)
-				bb.expand(E[e][p]);						//expand the bounding box to include the point
-		}
-
-		return bb;								//return the bounding box
-	}
-
-	///render cylinder based on points from the top/bottom hat
-	///@param C1 set of points from one of the hat
-	void renderCylinder(std::vector< stim::vec3<T> > C1, std::vector< stim::vec3<T> > C2, stim::vec3<T> P1, stim::vec3<T> P2) {
-		glBegin(GL_QUAD_STRIP);
-		for (unsigned i = 0; i < C1.size(); i++) {				// for every point on the circle
-			stim::vec3<T> n1 = C1[i] - P1;						// set normal vector for every vertex on the quads
-			stim::vec3<T> n2 = C2[i] - P2;
-			n1 = n1.norm();
-			n2 = n2.norm();
-			glNormal3f(n1[0], n1[1], n1[2]);
-			glVertex3f(C1[i][0], C1[i][1], C1[i][2]);
-			glNormal3f(n2[0], n2[1], n2[2]);
-			glVertex3f(C2[i][0], C2[i][1], C2[i][2]);											
-		}	
-		glEnd();
-	}
-
-	///render the vertex as sphere based on glut build-in function
-	///@param x, y, z are the three coordinates of the center point
-	///@param radius is the radius of the sphere
-	///@param subdivisions is the slice/stride along/around z-axis
-	void renderBall(T x, T y, T z, T radius, int subdivisions) {
-		glPushMatrix();
-		glTranslatef(x, y, z);
-		glutSolidSphere(radius, subdivisions, subdivisions);
-		glPopMatrix();
-	}
-
-	///render the vertex as sphere based on transformation
-	///@param x, y, z are the three coordinates of the center point
-	///@param radius is the radius of the sphere
-	///@param slice is the number of subdivisions around the z-axis
-	///@param stack is the number of subdivisions along the z-axis
-	void renderBall(T x, T y, T z, T radius, T slice, T stack) {
-		T step_z = stim::PI / slice;					// step angle along z-axis
-		T step_xy = 2 * stim::PI / stack;				// step angle in xy-plane
-		T xx[4], yy[4], zz[4];							// store coordinates
-
-		T angle_z = 0.0;								// start angle
-		T angle_xy = 0.0;
-
-		glBegin(GL_QUADS);
-		for (unsigned i = 0; i < slice; i++) {			// around the z-axis
-			angle_z = i * step_z;						// step step_z each time
-
-			for (unsigned j = 0; j < stack; j++) {		// along the z-axis
-				angle_xy = j * step_xy;					// step step_xy each time, draw floor by floor
-
-				xx[0] = radius * std::sin(angle_z) * std::cos(angle_xy);	// four vertices
-				yy[0] = radius * std::sin(angle_z) * std::sin(angle_xy);
-				zz[0] = radius * std::cos(angle_z);
-
-				xx[1] = radius * std::sin(angle_z + step_z) * std::cos(angle_xy);
-				yy[1] = radius * std::sin(angle_z + step_z) * std::sin(angle_xy);
-				zz[1] = radius * std::cos(angle_z + step_z);
-
-				xx[2] = radius * std::sin(angle_z + step_z) * std::cos(angle_xy + step_xy);
-				yy[2] = radius * std::sin(angle_z + step_z) * std::sin(angle_xy + step_xy);
-				zz[2] = radius * std::cos(angle_z + step_z);
-
-				xx[3] = radius * std::sin(angle_z) * std::cos(angle_xy + step_xy);
-				yy[3] = radius * std::sin(angle_z) * std::sin(angle_xy + step_xy);
-				zz[3] = radius * std::cos(angle_z);
-
-				for (unsigned k = 0; k < 4; k++) {
-					glVertex3f(x + xx[k], y + yy[k], z + zz[k]);			// draw the floor plane
-				}
-			}
-		}
-		glEnd();
-	}
-
-	/// Render the network centerline as a series of line strips.
-	/// glCenterline0 is for only one input
-	void glCenterline0(){
-		if (!glIsList(dlist)) {					//if dlist isn't a display list, create it
-			dlist = glGenLists(1);				//generate a display list
-			glNewList(dlist, GL_COMPILE);		//start a new display list
-			for (unsigned e = 0; e < E.size(); e++) {				//for each edge in the network
-				glBegin(GL_LINE_STRIP);
-				for (unsigned p = 0; p < E[e].size(); p++) {			//for each point on that edge
-					glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);		//set the vertex position based on the current point
-					glTexCoord1f(0);									//set white color
-				}
-				glEnd();
-			}
-			glEndList();						//end the display list
-		}
-		glCallList(dlist);					// render the display list
-	}
-
-	///render the network centerline as a series of line strips(when loading at least two networks, otherwise using glCenterline0())
-	///colors are based on metric values
-	void glCenterline(){
-
-		if(!glIsList(dlist)){					//if dlist isn't a display list, create it
-			dlist = glGenLists(1);				//generate a display list
-			glNewList(dlist, GL_COMPILE);		//start a new display list
-			for(unsigned e = 0; e < E.size(); e++){				//for each edge in the network
-				//unsigned errormag_id = E[e].nmags() - 1;
-				glBegin(GL_LINE_STRIP);
-				for(unsigned p = 0; p < E[e].size(); p++){				//for each point on that edge
-					glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);		//set the vertex position based on the current point
-					glTexCoord1f(E[e].r(p));							//set the texture coordinate based on the specified magnitude index
-				}
-				glEnd();
-			}
-			glEndList();						//end the display list
-		}		
-		glCallList(dlist);						//render the display list
-	}
-
-	///render the network cylinder as a series of tubes(when only one network loaded) 
-	void glCylinder0() {
-
-		float r1, r2;
-		if (!glIsList(dlist)) {										// if dlist isn't a display list, create it
-			dlist = glGenLists(1);									// generate a display list
-			glNewList(dlist, GL_COMPILE);							// start a new display list
-			for (unsigned e = 0; e < E.size(); e++) {				// for each edge in the network
-				for (unsigned p = 1; p < E[e].size(); p++) {		// for each point on that edge
-					stim::circle<T> C1 = E[e].circ(p - 1);
-					stim::circle<T> C2 = E[e].circ(p);
-					r1 = E[e].r(p - 1);
-					r2 = E[e].r(p);	
-					C1.set_R(r1);									// re-scale the circle to the same
-					C2.set_R(r2);
-					std::vector< stim::vec3<T> > Cp1 = C1.glpoints(20);// get 20 points on the circle plane
-					std::vector< stim::vec3<T> > Cp2 = C2.glpoints(20);
-					glBegin(GL_QUAD_STRIP);
-					for (unsigned i = 0; i < Cp1.size(); i++) {
-						glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
-						glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
-					}
-					glEnd();
-				}													// set the texture coordinate based on the specified magnitude index
-			}
-			for (unsigned n = 0; n < V.size(); n++) {
-				for (unsigned i = 0; i < E.size(); i++) {
-					if (E[i].v[0] == n) {
-						r1 = E[i].r(0);
-						break;
-					}
-					else if (E[i].v[1] == n) {
-						r1 = E[i].r(E[i].size() - 1);
-						break;
-					}
-				}
-				renderBall(V[n][0], V[n][1], V[n][2], r1, 20);
-			}
-			glEndList();						// end the display list
-		}
-		glCallList(dlist);						// render the display list
-	}
-
-	///render the network cylinder as a series of tubes
-	///colors are based on metric values 
-	void glCylinder(float sigma, float radius) {
-
-		if (radius != sigma)					// if render radius was changed by user, create a new display list
-			glDeleteLists(dlist, 1);
-		if (!glIsList(dlist)) {										// if dlist isn't a display list, create it
-			dlist = glGenLists(1);									// generate a display list
-			glNewList(dlist, GL_COMPILE);							// start a new display list
-			for (unsigned e = 0; e < E.size(); e++) {				// for each edge in the network
-				for (unsigned p = 1; p < E[e].size(); p++) {		// for each point on that edge
-					stim::circle<T> C1 = E[e].circ(p - 1);
-					stim::circle<T> C2 = E[e].circ(p);
-					C1.set_R(2*radius);								// re-scale the circle to the same
-					C2.set_R(2*radius);
-					std::vector< stim::vec3<T> > Cp1 = C1.glpoints(20);// get 20 points on the circle plane
-					std::vector< stim::vec3<T> > Cp2 = C2.glpoints(20);
-					glBegin(GL_QUAD_STRIP);
-					for (unsigned i = 0; i < Cp1.size(); i++) {
-						stim::vec3<T> n1 = Cp1[i] - E[e][p - 1];	// set normal vector for every vertex on the quads
-						stim::vec3<T> n2 = Cp2[i] - E[e][p];
-						n1 = n1.norm();
-						n2 = n2.norm();
-
-						glNormal3f(n1[0], n1[1], n1[2]);
-						glTexCoord1f(E[e].r(p - 1));
-						glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
-						glNormal3f(n2[0], n2[1], n2[2]);
-						glTexCoord1f(E[e].r(p));
-						glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
-					}
-					glEnd();
-				}													// set the texture coordinate based on the specified magnitude index
-			}
-			for (unsigned n = 0; n < V.size(); n++) {
-				size_t num = V[n].e[0].size();					// store the number of outgoing edge of that vertex
-				if (num != 0) {									// if it has outgoing edge
-					unsigned idx = V[n].e[0][0];				// find the index of first outgoing edge of that vertex
-					glTexCoord1f(E[idx].r(0));					// bind the texture as metric of first point on that edge
-				}					
-				else {
-					unsigned idx = V[n].e[1][0];				// find the index of first incoming edge of that vertex
-					glTexCoord1f(E[idx].r(E[idx].size() - 1));	// bind the texture as metric of last point on that edge
-				}
-				renderBall(V[n][0], V[n][1], V[n][2], 2*radius, 20);
-			}
-			glEndList();						// end the display list
-		}
-		glCallList(dlist);						// render the display list
-	}
-
-	///render a T as a adjoint network of GT in transparancy(darkgreen, overlaid)
-	void glAdjointCylinder(float sigma, float radius) {
-		
-		if (radius != sigma)					// if render radius was changed by user, create a new display list
-			glDeleteLists(dlist + 4, 1);
-		if (!glIsList(dlist + 4)) {
-			glNewList(dlist + 4, GL_COMPILE);
-			for (unsigned e = 0; e < E.size(); e++) {				// for each edge in the network
-				for (unsigned p = 1; p < E[e].size(); p++) {		// for each point on that edge
-					stim::circle<T> C1 = E[e].circ(p - 1);
-					stim::circle<T> C2 = E[e].circ(p);
-					C1.set_R(2 * radius);							// scale the circle to the same
-					C2.set_R(2 * radius);
-					std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
-					std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
-					glBegin(GL_QUAD_STRIP);
-					for (unsigned i = 0; i < Cp1.size(); i++) {		// for every point on the circle(+1 means closing the circle)
-						glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
-						glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
-					}
-					glEnd();
-				}								// set the texture coordinate based on the specified magnitude index
-			}
-			for (unsigned n = 0; n < V.size(); n++) {
-				size_t num = V[n].e[0].size();					// store the number of outgoing edge of that vertex
-				if (num != 0) {									// if it has outgoing edge
-					unsigned idx = V[n].e[0][0];				// find the index of first outgoing edge of that vertex 
-				}
-				else {
-					unsigned idx = V[n].e[1][0];				// find the index of first incoming edge of that vertex
-				}
-				renderBall(V[n][0], V[n][1], V[n][2], 2 * radius, 20);
-			}
-			glEndList();
-		}
-		glCallList(dlist + 4);
-	}
-
-	///render the network cylinder as series of tubes
-	///@param I is a indicator: 0 -> GT, 1 -> T
-	///@param map is the mapping relationship between two networks
-	///@param colormap is the random generated color set for render
-	void glRandColorCylinder(int I, std::vector<unsigned> map, std::vector<T> colormap, float sigma, float radius) {
-		
-		if (radius != sigma)									// if render radius was changed by user, create a new display list
-			glDeleteLists(dlist + 2, 1);
-		if (!glIsList(dlist + 2)) {								// if dlist isn't a display list, create it
-			glNewList(dlist + 2, GL_COMPILE);					// start a new display list
-			for (unsigned e = 0; e < E.size(); e++) {			// for each edge in the network
-				if (map[e] != unsigned(-1)) {
-					if (I == 0) {								// if it is to render GT
-						glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
-					}									
-					else {										// if it is to render T
-						glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
-					}
-
-					for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
-						stim::circle<T> C1 = E[e].circ(p - 1);
-						stim::circle<T> C2 = E[e].circ(p);
-						C1.set_R(2*radius);						// re-scale the circle to the same
-						C2.set_R(2*radius);
-						std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
-						std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
-						renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]);
-					}
-				}
-				else {
-					glColor3f(1.0, 1.0, 1.0);					// white color for the un-mapping edges
-					for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
-						stim::circle<T> C1 = E[e].circ(p - 1);
-						stim::circle<T> C2 = E[e].circ(p);
-						C1.set_R(2*radius);						// scale the circle to the same
-						C2.set_R(2*radius);
-						std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
-						std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
-						renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]);
-					}
-				}
-			}
-			for (unsigned n = 0; n < V.size(); n++) {
-				size_t num_edge = V[n].e[0].size() + V[n].e[1].size();
-				if (num_edge > 1) {					// if it is the joint vertex
-					glColor3f(0.3, 0.3, 0.3);		// gray
-					renderBall(V[n][0], V[n][1], V[n][2], 3*radius, 20);
-				}
-				else {								// if it is the terminal vertex
-					glColor3f(0.6, 0.6, 0.6);		// more white gray
-					renderBall(V[n][0], V[n][1], V[n][2], 3*radius, 20);
-				}
-			}
-			glEndList();
-		}
-		glCallList(dlist + 2);
-	}
-
-	///render the network centerline as a series of line strips in random different color
-	///@param I is a indicator: 0 -> GT, 1 -> T
-	///@param map is the mapping relationship between two networks
-	///@param colormap is the random generated color set for render
-	void glRandColorCenterline(int I, std::vector<unsigned> map, std::vector<T> colormap) {
-		if (!glIsList(dlist + 2)) {
-			glNewList(dlist + 2, GL_COMPILE);
-			for (unsigned e = 0; e < E.size(); e++) {
-				if (map[e] != unsigned(-1)) {					// if it has corresponding edge in another network
-					if (I == 0)									// if it is to render GT
-						glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
-					else										// if it is to render T
-						glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
-
-					glBegin(GL_LINE_STRIP);
-					for (unsigned p = 0; p < E[e].size(); p++) {
-						glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
-					}
-					glEnd();
-				}
-				else {
-					glColor3f(1.0, 1.0, 1.0);					// white color for the un-mapping edges
-					glBegin(GL_LINE_STRIP);
-					for (unsigned p = 0; p < E[e].size(); p++) {
-						glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
-					}
-					glEnd();
-				}
-			}
-			glEndList();
-		}
-		glCallList(dlist + 2);
-	}
-
-	void glAdjointCenterline() {
-		if (!glIsList(dlist + 4)) {
-			glNewList(dlist + 4, GL_COMPILE);
-			for (unsigned e = 0; e < E.size(); e++) {				//for each edge in the network
-																
-				glBegin(GL_LINE_STRIP);
-				for (unsigned p = 0; p < E[e].size(); p++) {			//for each point on that edge
-					glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);		//set the vertex position based on the current point
-					glTexCoord1f(E[e].r(p));							//set the texture coordinate based on the specified magnitude index
-				}
-				glEnd();
-			}
-			glEndList();
-		}
-		glCallList(dlist + 4);
-	}
-
-	// highlight the difference part
-	void glDifferenceCylinder(int I, std::vector<unsigned> map, std::vector<T> colormap, float sigma, float radius) {
-
-		if (radius != sigma)									// if render radius was changed by user, create a new display list
-			glDeleteLists(dlist + 6, 1);
-		if (!glIsList(dlist + 6)) {								// if dlist isn't a display list, create it
-			glNewList(dlist + 6, GL_COMPILE);					// start a new display list
-			for (unsigned e = 0; e < E.size(); e++) {			// for each edge in the network
-				if (map[e] != unsigned(-1)) {
-					glEnable(GL_BLEND);									//enable color blend
-					glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);	//set blend function
-					glDisable(GL_DEPTH_TEST);							//should disable depth
-
-					if (I == 0) {								// if it is to render GT
-						glColor4f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2], 0.1);
-					}
-					else {										// if it is to render T
-						glColor4f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2], 0.1);
-					}
-
-					for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
-						stim::circle<T> C1 = E[e].circ(p - 1);
-						stim::circle<T> C2 = E[e].circ(p);
-						C1.set_R(2 * radius);						// re-scale the circle to the same
-						C2.set_R(2 * radius);
-						std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
-						std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
-						renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]);
-					}
-					glDisable(GL_BLEND);
-					glEnable(GL_DEPTH_TEST);
-				}
-				else {
-					glColor3f(1.0, 1.0, 1.0);					// white color for the un-mapping edges
-					for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
-						stim::circle<T> C1 = E[e].circ(p - 1);
-						stim::circle<T> C2 = E[e].circ(p);
-						C1.set_R(2 * radius);						// scale the circle to the same
-						C2.set_R(2 * radius);
-						std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
-						std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
-						renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]);
-					}
-				}
-			}
-			for (unsigned n = 0; n < V.size(); n++) {
-
-				size_t num_edge = V[n].e[0].size() + V[n].e[1].size();
-				if (num_edge > 1) {					// if it is the joint vertex
-					glColor4f(0.3, 0.3, 0.3, 0.1);		// gray
-					renderBall(V[n][0], V[n][1], V[n][2], 3 * radius, 20);
-				}
-				else {								// if it is the terminal vertex
-					glColor4f(0.6, 0.6, 0.6, 0.1);		// more white gray
-					renderBall(V[n][0], V[n][1], V[n][2], 3 * radius, 20);
-				}
-			}
-			glEndList();
-		}
-		glCallList(dlist + 6);
-	}
-
-	//void glRandColorCenterlineGT(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap){
-	//	if(!glIsList(dlist1)){
-	//		dlist1 = glGenLists(1);
-	//		glNewList(dlist1, GL_COMPILE);
-	//		for(unsigned e = 0; e < E.size(); e++){
-	//			if(map[e] != unsigned(-1)){
-	//				glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
-	//				glBegin(GL_LINE_STRIP);
-	//				for(unsigned p = 0; p < E[e].size(); p++){
-	//					glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
-	//				}
-	//				glEnd();
-	//				for (unsigned p = 0; p < E[e].size() - 1; p++) {
-	//					renderCylinder(E[e][p][0], E[e][p][1], E[e][p][2], E[e][p + 1][0], E[e][p + 1][1], E[e][p + 1][2], 10, 20);
-	//				}
-	//			}
-	//			else{
-	//				glColor3f(1.0, 1.0, 1.0);
-	//				glBegin(GL_LINE_STRIP);
-	//				for(unsigned p = 0; p < E[e].size(); p++){
-	//					glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
-	//				}
-	//				glEnd();
-	//			}
-	//		}
-	//		for (unsigned v = 0; v < V.size(); v++) {
-	//			size_t num_edge = V[v].e[0].size() + V[v].e[1].size();
-	//			if (num_edge > 1) {
-	//				glColor3f(0.3, 0.3, 0.3);		// gray color for vertex
-	//				renderBall(V[v][0], V[v][1], V[v][2], 20, 20);
-	//			}
-	//		}
-	//		glEndList();
-	//	}
-	//	glCallList(dlist1);
-	//}
-
-	//void glRandColorCenterlineT(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap){
-	//	if(!glIsList(dlist2)){
-	//		dlist2 = glGenLists(1);
-	//		glNewList(dlist2, GL_COMPILE);
-	//		for(unsigned e = 0; e < E.size(); e++){
-	//			if(map[e] != unsigned(-1)){
-	//				glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
-	//				glBegin(GL_LINE_STRIP);
-	//				for(unsigned p = 0; p < E[e].size(); p++){
-	//					glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
-	//				}
-	//				glEnd();
-	//				for (unsigned p = 0; p < E[e].size() - 1; p++) {
-	//					renderCylinder(E[e][p][0], E[e][p][1], E[e][p][2], E[e][p + 1][0], E[e][p + 1][1], E[e][p + 1][2], 10, 20);
-	//				}
-	//			}
-	//			else{
-	//				glColor3f(1.0, 1.0, 1.0);
-	//				glBegin(GL_LINE_STRIP);
-	//				for(unsigned p = 0; p < E[e].size(); p++){
-	//					glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
-	//				}
-	//				glEnd();
-	//			}
-	//		}
-	//		for (unsigned v = 0; v < V.size(); v++) {
-	//			size_t num_edge = V[v].e[0].size() + V[v].e[1].size();
-	//			if (num_edge > 1) {
-	//				glColor3f(0.3, 0.3, 0.3);		// gray color for vertex
-	//				renderBall(V[v][0], V[v][1], V[v][2], 20, 20);
-	//			}
-	//		}
-	//		glEndList();
-	//	}
-	//	glCallList(dlist2);
-	//}
-
-
-	//void renderCylinder(T x1, T y1, T z1, T x2, T y2, T z2, T radius, int subdivisions) {
-	//	T dx = x2 - x1;
-	//	T dy = y2 - y1;
-	//	T dz = z2 - z1;
-	//	/// handle the degenerate case with an approximation
-	//	if (dz == 0)
-	//		dz = .00000001;
-	//	T d = sqrt(dx*dx + dy*dy + dz*dz);					
-	//	T ax = 57.2957795*acos(dz / d);						// 180°/pi
-	//	if (dz < 0.0)
-	//		ax = -ax;
-	//	T rx = -dy*dz;
-	//	T ry = dx*dz;
-
-	//	glPushMatrix();
-	//	glTranslatef(x1, y1, z1);
-	//	glRotatef(ax, rx, ry, 0.0);
-
-	//	glutSolidCylinder(radius, d, subdivisions, 1);
-	//	glPopMatrix();
-	//}
-
-
-	/// render the network centerline from swc file as a series of strips in different colors based on the neuronal type
-	/// glCenterline0_swc is for only one input
-	/*void glCenterline0_swc() {
-	if (!glIsList(dlist)) {						// if dlist isn't a display list, create it
-	dlist = glGenLists(1);					// generate a display list
-	glNewList(dlist, GL_COMPILE);			// start a new display list
-	for (unsigned e = 0; e < E.size(); e++) {
-	int type = NT[e];					// get the neuronal type
-	switch (type) {
-	case 0:
-	glColor3f(1.0f, 1.0f, 1.0f);	// white for undefined
-	glBegin(GL_LINE_STRIP);
-	for (unsigned p = 0; p < E[e].size(); p++) {
-	glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
-	}
-	glEnd();
-	break;
-	case 1:
-	glColor3f(1.0f, 0.0f, 0.0f);	// red for soma
-	glBegin(GL_LINE_STRIP);
-	for (unsigned p = 0; p < E[e].size(); p++) {
-	glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
-	}
-	glEnd();
-	break;
-	case 2:
-	glColor3f(1.0f, 0.5f, 0.0f);	// orange for axon
-	glBegin(GL_LINE_STRIP);
-	for (unsigned p = 0; p < E[e].size(); p++) {
-	glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
-	}
-	glEnd();
-	break;
-	case 3:
-	glColor3f(1.0f, 1.0f, 0.0f);	// yellow for undefined
-	glBegin(GL_LINE_STRIP);
-	for (unsigned p = 0; p < E[e].size(); p++) {
-	glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
-	}
-	glEnd();
-	break;
-	case 4:
-	glColor3f(0.0f, 1.0f, 0.0f);	// green for undefined
-	glBegin(GL_LINE_STRIP);
-	for (unsigned p = 0; p < E[e].size(); p++) {
-	glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
-	}
-	glEnd();
-	break;
-	case 5:
-	glColor3f(0.0f, 1.0f, 1.0f);	// verdant for undefined
-	glBegin(GL_LINE_STRIP);
-	for (unsigned p = 0; p < E[e].size(); p++) {
-	glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
-	}
-	glEnd();
-	break;
-	case 6:
-	glColor3f(0.0f, 0.0f, 1.0f);	// blue for undefined
-	glBegin(GL_LINE_STRIP);
-	for (unsigned p = 0; p < E[e].size(); p++) {
-	glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
-	}
-	glEnd();
-	break;
-	case 7:
-	glColor3f(0.5f, 0.0f, 1.0f);	// purple for undefined
-	glBegin(GL_LINE_STRIP);
-	for (unsigned p = 0; p < E[e].size(); p++) {
-	glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
-	}
-	glEnd();
-	break;
-	}
-	}
-	glEndList();						//end the display list
-	}
-	glCallList(dlist);					// render the display list
-	}*/
-
-	///render the network cylinder as a series of tubes
-	///colors are based on metric values 
-	//void glCylinder(float sigma) {
-	//	if (!glIsList(dlist)) {					//if dlist isn't a display list, create it
-	//		dlist = glGenLists(1);				//generate a display list
-	//		glNewList(dlist, GL_COMPILE);		//start a new display list
-	//		for (unsigned e = 0; e < E.size(); e++) {				//for each edge in the network
-	//			for (unsigned p = 1; p < E[e].size() - 1; p++) {	// for each point on that edge
-	//				stim::circle<T> C1 = E[e].circ(p - 1);
-	//				stim::circle<T> C2 = E[e].circ(p);
-	//				C1.set_R(2.5*sigma);							// scale the circle to the same
-	//				C2.set_R(2.5*sigma);
-	//				std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
-	//				std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
-	//				glBegin(GL_QUAD_STRIP);
-	//				for (unsigned i = 0; i < Cp1.size(); i++) {		// for every point on the circle(+1 means closing the circle)
-	//					glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
-	//					glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
-	//					glTexCoord1f(E[e].r(p));
-	//				}
-	//				glEnd();
-	//			}								//set the texture coordinate based on the specified magnitude index
-	//		}
-	//		for (unsigned n = 0; n < V.size(); n++) {
-	//			size_t num = V[n].e[0].size();					//store the number of outgoing edge of that vertex
-	//			if (num != 0) {									//if it has outgoing edge
-	//				unsigned idx = V[n].e[0][0];				//find the index of first outgoing edge of that vertex 
-	//				glTexCoord1f(E[idx].r(0));					//bind the texture as metric of first point on that edge
-	//			}
-	//			else {
-	//				unsigned idx = V[n].e[1][0];				//find the index of first incoming edge of that vertex
-	//				glTexCoord1f(E[idx].r(E[idx].size() - 1));	//bind the texture as metric of last point on that edge
-	//			}
-	//			renderBall(V[n][0], V[n][1], V[n][2], 2.5*sigma, 20);
-	//		}
-	//		glEndList();						//end the display list
-	//	}
-	//	glCallList(dlist);						//render the display list
-	//}
-
-};		//end stim::gl_network class
-};		//end stim namespace
-
-
-
+#ifndef JACK_GL_NETWORK
+#define JACK_GL_NETWORK
+
+#include"network.h"
+#include<GL/glut.h>
+#include<stim/visualization/aaboundingbox.h>
+
+namespace jack {
+	
+	template<typename T>
+	class gl_network : public jack::network<T> {
+	private:
+		std::vector<T> ecolor;				// colormap for edges
+		std::vector<T> vcolor;				// colormap for vertices
+		GLint subdivision;					// rendering subdivision
+
+		/// basic geometry rendering
+		// main sphere rendering function
+		void sphere(T x, T y, T z, T radius) {
+			glPushMatrix();
+			glTranslatef((GLfloat)x, (GLfloat)y, (GLfloat)z);
+			glutSolidSphere((double)radius, subdivision, subdivision);
+			glPopMatrix();
+		}
+		// rendering sphere using glut function
+		void draw_sphere(T x, T y, T z, T radius, T alpha = 1.0f) {
+			if (alpha != 1.0f) {
+				glEnable(GL_BLEND);									// enable color blend
+				glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);	// set blend function
+				glDisable(GL_DEPTH_TEST);							// disable depth buffer
+				glColor4f(vcolor[0], vcolor[1], vcolor[2], alpha);	// set color
+				sphere(x, y, z, radius, subdivision);
+				glDisable(GL_BLEND);								// disable color blend
+				glEnable(GL_DEPTH_TEST);							// enbale depth buffer again
+			}
+			else {
+				glColor3f(vcolor[0], vcolor[1], vcolor[2]);
+				sphere(x, y, z, radius, subdivision);
+			}
+		}
+		// rendering sphere using quads
+		void draw_sphere(T x, T y, T z, T radius, T alpha = 1.0f) {
+			GLint stack = subdivision;
+			GLint slice = subdivision;
+			
+			if (alpha != 1.0f) {
+				glEnable(GL_BLEND);									// enable color blend
+				glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);	// set blend function
+				glDisable(GL_DEPTH_TEST);							// disable depth buffer
+			}
+			glColor4f(vcolor[0], vcolor[1], vcolor[2], alpha);		// set color
+
+			T step_z = stim::PI / slice;					// step angle along z-axis
+			T step_xy = 2 * stim::PI / stack;				// step angle in xy-plane
+			T xx[4], yy[4], zz[4];							// store coordinates
+
+			T angle_z = 0.0;								// start angle
+			T angle_xy = 0.0;
+
+			glBegin(GL_QUADS);
+			for (unsigned i = 0; i < slice; i++) {			// around the z-axis
+				angle_z = i * step_z;						// step step_z each time
+
+				for (unsigned j = 0; j < stack; j++) {		// along the z-axis
+					angle_xy = j * step_xy;					// step step_xy each time, draw floor by floor
+
+					xx[0] = radius * std::sin(angle_z) * std::cos(angle_xy);	// four vertices
+					yy[0] = radius * std::sin(angle_z) * std::sin(angle_xy);
+					zz[0] = radius * std::cos(angle_z);
+
+					xx[1] = radius * std::sin(angle_z + step_z) * std::cos(angle_xy);
+					yy[1] = radius * std::sin(angle_z + step_z) * std::sin(angle_xy);
+					zz[1] = radius * std::cos(angle_z + step_z);
+
+					xx[2] = radius * std::sin(angle_z + step_z) * std::cos(angle_xy + step_xy);
+					yy[2] = radius * std::sin(angle_z + step_z) * std::sin(angle_xy + step_xy);
+					zz[2] = radius * std::cos(angle_z + step_z);
+
+					xx[3] = radius * std::sin(angle_z) * std::cos(angle_xy + step_xy);
+					yy[3] = radius * std::sin(angle_z) * std::sin(angle_xy + step_xy);
+					zz[3] = radius * std::cos(angle_z);
+
+					for (unsigned k = 0; k < 4; k++) {
+						glVertex3f(x + xx[k], y + yy[k], z + zz[k]);			// draw the floor plane
+					}
+				}
+			}
+			glEnd();
+
+			if (alpha != 1.0f) {
+				glDisable(GL_BLEND);								// disable color blend
+				glEnable(GL_DEPTH_TEST);							// enbale depth buffer again
+			}
+		}
+
+		// render edge as cylinders
+		void draw_cylinder(T scale = 1.0f) {
+			stim::circle<T> C1, C2;					// temp circles
+			std::vector<stim::vec3<T> > Cp1, Cp2;	// temp lists for storing points on those circles
+			T r1, r2;								// temp radii
+
+			size_t num_edge = E.size();
+			size_t num_vertex = V.size();
+
+			for (size_t i = 0; i < num_edge; i++) {			// for every edge
+				glColor3f(ecolor[i * 3 + 0], ecolor[i * 3 + 1], ecolor[i * 3 + 2]);
+				for (size_t j = 1; j < num_vertex; j++) {	// for every vertex except first one
+					C1 = E[i].circ(j - 1);			// get the first circle plane of this segment
+					C2 = E[i].circ(j);				// get the second circle plane of this segment
+					r1 = E[i].r(j - 1);				// get the radius at first point
+					r2 = E[i].r(j);					// get the radius at second point
+					C1.set_R(scale * r1);			// rescale
+					C2.set_R(scale * r2);
+					Cp1 = C1.glpoints((unsigned)subdivision);	// get 20 points on first circle plane
+					Cp2 = C2.glpoints((unsigned)subdivision);	// get 20 points on second circle plane
+
+					glBegin(GL_QUAD_STRIP);
+					for (size_t k = 0; k < subdivision + 1; k++) {
+						glVertex3f(Cp1[k][0], Cp1[k][1], Cp1[k][2]);
+						glVertex3f(Cp2[k][0], Cp2[k][1], Cp2[k][2]);
+					}
+					glEnd();
+				}
+			}
+		}
+
+
+	protected:
+		using jack::network<T>::E;
+		using jack::network<T>::T;
+
+		GLuint dlist;
+
+	public:
+
+		/// constructors
+		// empty constructor
+		gl_network() : jack::network<T>() {
+			dlist = 0;
+			subdivision = 20;
+		}
+
+		gl_network(jack::network<T> N) : jack::network<T>(N) {
+			dlist = 0;
+			subdivision = 20;
+			ecolor.resize(N.edges * 3, 0.0f);		// default black color
+			vcolor.resize(N.vertices * 3, 0.0f);
+		}
+
+		// compute the smallest bounding box of current network
+		aabboundingbox<T> boundingbox() {
+			aaboundingbox<T> bb;			// create a bounding box object
+
+			for (size_t i = 0; i < E.size(); i++) 			// for every edge
+				for (size_t j = 0; j < E[i].size(); j++) 	// for every point on that edge
+					bb.expand(E[i][j]);		// expand the bounding box to include that point
+
+			return bb;
+		}
+
+		// change subdivision
+		void change_subdivision(GLint value) {
+			subdivision = value;
+		}
+
+
+		/// rendering functions
+		// render centerline
+		void centerline() {		
+			size_t num = E.size();				// get the number of edges
+			for (size_t i = 0; i < num; i++) {
+				glColor3f(colormap[i * 3 + 0], colormap[i * 3 + 1], colormap[i * 3 + 2]);
+				glBegin(GL_LINE_STRIP);
+				for (size_t j = 0; j < E[i].size(); j++)
+					glVertex3f(E[i][j][0], E[i][j][1], E[i][j][2]);
+				glEnd();
+			}
+		}
+
+
+		// render network as bunch of spheres and cylinders
+		void network(T scale = 1.0f) {
+			stim::vec3<T> v1, v2;			// temp vertices for rendering
+			T r1, r2;						// temp radii for rendering
+
+			if (!glIsList(dlist)) {					// if dlist is not a display list, create one
+				dlist = glGenLists(1);				// generate a display list
+				glNewList(dlist, GL_COMPILE);		// start a new display list
+			
+				// render every vertex as a sphere
+				size_t num = V.size();
+				for (size_t i = 0; i < num; i++) {
+					v1 = stim::vec3<T>(V[i][0], V[i][1], V[i][2]);	// get the vertex for rendering
+					r1 = (*this).network<T>::r(i);
+					draw_sphere(v1[0], v1[1], v1[2], r1 * scale);
+				}
+
+				// render every edge as a cylinder
+				draw_cylinder(scale);
+
+				glEndList();						// end the display list
+			}
+			glCallList(dlist);
+		}
+	};
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 #endif
 \ No newline at end of file
@@ -62,6 +62,95 @@ namespace stim {
 			if (dlist) glDeleteLists(dlist, 1);		//delete the display list when the object is destroyed
 		}
  
+                ///saves the coeffients of a spherical harmonics in a text documents.
+                void save_coeffs(std::string filename)
+                {
+                        std::ofstream myfile;
+                        myfile.open(filename.c_str());
+                        for(int i = 0; i < Sd.C.size()-1; i++)
+                        {
+                                myfile << Sd.C[i] << ",";
+                        }
+                        myfile << Sd.C[Sd.C.size()-1] << std::endl;
+                        for(int i = 0; i < Sc.C.size()-1; i++)
+                        {
+                                myfile << Sc.C[i] << ",";
+                        }
+                        myfile << Sc.C[Sc.C.size()-1] << std::endl;
+                        myfile.close();
+                }
+
+                void save(std::string filename)
+                {
+
+                        if (!tex) {
+                                init_tex();
+                        }
+
+                        stim::obj<T> object;
+                        size_t theta_i, phi_i;
+                        T d_theta = (T)stim::TAU / (T)N;
+                        T d_phi = (T)stim::PI / (T)(N-1);
+                        object.matKd("sfunc.jpg");
+                        for (phi_i = 1; phi_i < N; phi_i++) {
+                                T phi = phi_i * d_phi;
+                                object.Begin(OBJ_TRIANGLE_STRIP);
+                                for (theta_i = 0; theta_i <= N; theta_i++) {
+                                        T theta = (N - theta_i) * d_theta;
+                                        float theta_t = 1 - (float)theta_i / (float)N;
+
+                                        T r;
+                                        if (!displacement) r = 1;                                               //if no displacement, set the r value to 1 (renders a unit sphere)
+                                        else r = Sd(theta, phi);                                                //otherwise calculate the displacement value
+
+                                        glColor3f(1.0f, 1.0f, 1.0f);
+                                        if (!colormap) {                                                                                //if no colormap is being rendered
+                                                if (r < 0) glColor3f(1.0, 0.0, 0.0);                            //if r is negative, render it red
+                                                else glColor3f(0.0, 1.0, 0.0);                                          //otherwise render in green
+                                        }
+                                        if (magnitude) {                                                                                //if the magnitude is being displaced, calculate the magnitude of r
+                                                if (r < 0) r = -r;
+                                        }
+                                        stim::vec3<T> s(r, theta, phi);
+                                        stim::vec3<T> c = s.sph2cart();
+                                        stim::vec3<T> n;                                                                                                                //allocate a value to store the normal
+                                        if (!displacement) n = c;                                                                                               //if there is no displacement, the normal is spherical
+                                        else n = Sd.dphi(theta, phi).cross(Sd.dtheta(theta, phi));                              //otherwise calculate the normal as the cross product of derivatives
+
+                                        object.TexCoord(theta_t, (float)phi_i / (float)N);
+                                        //std::cout << theta_t <<"        "<<(float)phi_i / (float)N << "----------------";
+                                        object.Normal(n[0], n[1], n[2]);
+                                        object.Vertex(c[0], c[1], c[2]);
+
+                                        T r1;
+                                        if (!displacement) r1 = 1;
+                                        else r1 = Sd(theta, phi - d_phi);
+                                        if (!colormap) {                                                                                //if no colormap is being rendered
+                                                if (r1 < 0)     glColor3f(1.0, 0.0, 0.0);                               //if r1 is negative, render it red
+                                                else glColor3f(0.0, 1.0, 0.0);                                          //otherwise render in green
+                                        }
+                                        if (magnitude) {                                                                                //if the magnitude is being rendered, calculate the magnitude of r
+                                                if (r1 < 0) r1 = -r1;
+                                        }
+                                        stim::vec3<T> s1(r1, theta, phi - d_phi);
+                                        stim::vec3<T> c1 = s1.sph2cart();
+                                        stim::vec3<T> n1;
+                                        if (!displacement) n1 = c1;
+                                        else n1 = Sd.dphi(theta, phi - d_phi).cross(Sd.dtheta(theta, phi - d_phi));
+
+                                        //std::cout << theta_t << "        " << (float)(phi_i - 1) / (float)N << std::endl;
+                                        object.TexCoord(theta_t, 1.0/(2*(N)) + (float)(phi_i-1) / (float)N);
+                                        object.Normal(n1[0], n1[1], n1[2]);
+                                        object.Vertex(c1[0], c1[1], c1[2]);
+                                }
+                                object.End();
+                        }
+                        object.matKd();
+                        object.save(filename);
+                }
+
+
+
 		/// This function renders the spherical harmonic to the current OpenGL context
 		void render() {
 			//glShadeModel(GL_FLAT);
@@ -570,6 +570,7 @@ public:
 		for (size_t i = 0; i < M.size(); i++) {
 			ss << M[i].str() << std::endl;
 		}
+		return ss.str();
 	}
  
 	obj(){