#ifndef STIM_NETWORK_H
#define STIM_NETWORK_H

#include <list>
#include <stdlib.h>
#include <sstream>
#include <fstream>
#include <algorithm>
#include <string.h>
#include <math.h>
#include <stim/math/vector.h>
#include <stim/visualization/obj.h>
#include <stim/visualization/fiber.h>
#include <ANN/ANN.h>
#include <boost/tuple/tuple.hpp>


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 ANN library.
*/


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 fiber<T>
	{
		public:
		int Nodes[2];		//unique id's designating the starting and ending

		///default constructor
		edge() : fiber<T>()
		{
			Nodes[1] = -1; Nodes[2] = -1;
		}

		///sized costructor with two known nodes.
		///@param startId: int storing idx assigned to Nodes[1].
		///@param endId: int storing idx assigned to Nodes[2].
		///@param n: int for the number of points in the fiber.
		edge(int startId, int endId, int n)
		 :fiber<T>(n)
		{
			Nodes[1] = startId; Nodes[2] = endId;
		}
		
		///constructor from a std::vector of stim::vecs of positions and radii.
		///@param pos: Vector of stim vecs storing the positions of the fiber.
		///@param mag: Vector of stim vecs storing the radii of the fiber.
		edge(std::vector< stim::vec<T> > pos, std::vector< stim::vec<T> > radii)
			 : fiber<T>(pos, radii)
		{
			Nodes[1] = -1; Nodes[2] = -1;
		}
		
		///constructor from an std::vector of stim::vecs of positions and radii as well as a known starting and ending node.
		///@param pos: Vector of stim vecs storing the positions of the fiber.
		///@param mag: Vector of stim vecs storing the radii of the fiber.
		///@param id1: int storing idx assigned to Nodes[1].
		///@param id2: int storing idx assigned to Nodes[2].
		edge(std::vector< stim::vec<T> > pos, std::vector< stim::vec<T> > radii, int id1, int id2)
			 : fiber<T>(pos, radii)
		{
			Nodes[1] = id1; Nodes[2] = id2;
		}


		edge(std::vector< stim::vec<T> > pos, std::vector<T> radii)
			 : fiber<T>(pos, radii)
		{
			Nodes[1] = -1; Nodes[2] = -1;
		}

		///sets the endpoints to the two int values.
		///@param int id1: index of Nodes[1].
		///@param int id2: index of Nodes[2].
		void
		setEndpoints(int id1, int id2)
		{
			Nodes[1] = id1; Nodes[2] = id2;
		}		

	};	
	
	///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 node
	{
		private:
			std::vector<int> edges;		//indices of edges connected to this node.
			stim::vec<T> p;			//position of this node in physical space.
		public:
			//no default constructor

			///@param pos: stim vec with the x, y, z position of the edge.
			///stim::vec constructure with a position but no attached edges.
			node(stim::vec<T> pos)
			{
				p = pos;
			}

			///@param pos: stim vec with the x, y, z position of the edge.
			///@param i: int i storing the index of an attached edge.
			//stim::vec constructor with a position and an attached edge.
			node(stim::vec<T> pos, int i)
			{
				p = pos;
				edges.push_back(i);
			}

			///@param x: x coordinate of the node..
			///@param y: y coordinate of the node..
			///@param z: z coordinate of the node..
			//float value constructor. 
			node(T x, T y, T z)
			{
				p = stim::vec<T>(x,y,z);
			}
			///@param x: x coordinate of the node..
			///@param y: y coordinate of the node..
			///@param z: z coordinate of the node..
			///@param i: int i storing the index of an attached edge.
			//float value constructor with an attached edge.
			node(T x, T y, T z, int i)
			{
				p = stim::vec<T>(x,y,z);
				edges.push_back(i);
			}

			///@param id: int index of the fiber to attach to this node.
			///adds the fiber to the rest of the fibers connected to this node.
			void
			addEdge(int id)
			{
			 	edges.push_back(id);
			}
	
			///returns the position of the node.
			stim::vec<T>
			getPosition()
			{	
				return p;
			}
		
			///@param id: int index of the fiber to detach to this node.
			///removes the edge from the list of the edges attached to this node.
			void
			removeEdge(int id)
			{
				for(int i = 0; i < edges.size(); i++)
				{
					if(edges[i] == id)
						edges.erase(edges.begin()+i);
				}
			}
			
			///returns and std::string with the position of this node.
			std::string
			str()
			{
				return p.str();
			}

			///returns and std::string with the list of edges attached to this node.
			std::string
			edges_to_str()
			{
				std::ostringstream ss;
				std::cout<<"here"<<std::endl;
//				ss << "[";
				for(int i = 0; i < edges.size()-1; i++)
				{
					std::cout<<"here"<<i<<std::endl;
					ss << edges[i] << ";";
					std::cout<<edges.size()-1<<std::endl;
				}
				ss << edges[edges.size()-1];
//				ss << "]";
				return ss.str();
			}
			
	};

	public:

	std::vector<edge*> E;       //list of pointers to edges.
	std::vector<node> V;	    //list of nodes.
	std::vector< stim::vec<T> > allVerticesList; // all nodes before sampling
	std::vector<stim::vec<T> > allVerticesAfterSampling ; //all vertices after sampling
	
	///Returns the number of edges in the network.
	unsigned int
	sizeE()
	{
		return E.size();
	}

	///Returns the number of nodes in the network.
	unsigned int
	sizeV()
	{
		return V.size();
	}
	unsigned int
	sizeAfterSamplingV()
	{
		return allVerticesAfterSampling.size();
	}
/*	//adds an edge from two std::vectors with positions and radii.
	void
	addEdge(std::vector< stim::vec<T> > pos, std::vector<stim::vec<T> > radii, int endId)
	{
		
		edge a(pos,radii, endId);
		E.push_back(a);
	} */

	///A complicated method that adds an edge to the network.
	///Most of this functionality will be moved into fiber.
	void
	addEdge(std::vector< stim::vec<T> > pos, std::vector<stim::vec<T> > radii, int startId, int endId)
	{
		//
		if(startId == -1 && endId == -1)
		{
			//case where the edge is neither starting nor ending in a fiber.
			//i. first fiber.
			
			//Add two nodes to the nodes vector
			V.push_back(node(pos[pos.size()-1]));
			V.push_back(node(pos[0]));

			//the edge will be connected to the nodes
			edge *a = new edge(pos,radii,(V.size()-2), (V.size()-1));
	
			//add fiber to fiber list.
			E.push_back(a);

			//The last two nodes added to V will "own" the last edge added to E.
			V[V.size()-1].addEdge(E.size()-1);
			V[V.size()-2].addEdge(E.size()-1);
			
		}
	
		else if(startId != -1 && endId == -1)
		{
			//case where the edge is starting with a fiber, but not ending in one.
			
			//split the fiber behind you into two.
			unsigned int point = E[startId]->nearest_idx(pos[0]);

			//split the hit fiber at point two parts temp[0], temp[1]
			std::vector < stim::fiber <T> > temp = E[startId]->split(point);
			if(temp.size() > 1)
			{
				//add the nearest point in the behind fiber into the hitting fiber.
				pos.insert(pos.begin(), E[startId]->nearest(pos[0]));
				stim::vec<T> v(E[startId]->radius(point), E[startId]->radius(point));
				radii.insert(radii.begin(), v);

				//reset the fiber at the endId to be a shorter fiber(temp[0]).
				std::vector<stim::vec<T> > ce = temp[0].centerline();
				std::vector<stim::vec<T> > cm = temp[0].centerlinemag();

				//remake the edge, such that the starting point of this edge
				//is the same the split point, but the end edge is the old end.
				V.push_back(node(ce[ce.size()-1]));
				int tempNodeId = E[startId]->Nodes[1];
				E[startId] = new edge(ce, cm, (V.size()-1), E[startId]->Nodes[2]);
				V[V.size()-1].addEdge(startId);
				

				//add the other part of the fiber (temp[1])
				ce = temp[1].centerline();
				cm = temp[1].centerlinemag();
				E.push_back(new edge(ce, cm,tempNodeId ,(V.size()-1)));
				V[V.size()-1].addEdge(E.size()-1);

				V[tempNodeId].removeEdge(startId);
				V[tempNodeId].addEdge(E.size()-1);
//				V[V.size()-1].removeEdge(startId);

				//make the new hitting fiber..
				V.push_back(node(pos[pos.size()-1]));
				edge *a = new edge(pos, radii, (V.size()-2), (V.size()-1));
				E.push_back(a);
				V[V.size()-1].addEdge(E.size()-1);
				V[V.size()-2].addEdge(E.size()-1);
			} else {
				stim::vec<T> pz = E[startId]->nearest(pos[0]);
				if((V[E[startId]->Nodes[1]].getPosition() - pz).len() <
					(V[E[startId]->Nodes[2]].getPosition() - pz).len())
				{
					unsigned int point = E[startId]->nearest_idx(pos[0]);
					pos.insert(pos.begin(), E[startId]->nearest(pos[0]));
					stim::vec<T> v(E[startId]->radius(point), E[startId]->radius(point));
					radii.insert(radii.begin(), v);
					V.push_back(node(pos[pos.size()-1]));
					edge *a = new edge(pos, radii, E[startId]->Nodes[1], (V.size()-1));
					E.push_back(a);
					
					V[V.size()-1].addEdge(E.size()-1);
					V[ E[startId]->Nodes[1]].addEdge(E.size()-1);

				} 
				else
				{
					unsigned int point = E[startId]->nearest_idx(pos[0]);
					pos.insert(pos.begin(), E[startId]->nearest(pos[0]));
					stim::vec<T> v(E[startId]->radius(point), E[startId]->radius(point));
					radii.insert(radii.begin(), v);
					V.push_back(node(pos[pos.size()-1]));
					edge *a = new edge(pos, radii, E[startId]->Nodes[2], (V.size()-1));
					E.push_back(a);
					
					V[V.size()-1].addEdge(E.size()-1);
					V[ E[startId]->Nodes[2]].addEdge(E.size()-1);
				}
			}
		}	

		//case where the edge is starting at a seedpoint but ends in a fiber.
		if(startId == -1 && endId != -1 && endId < sizeE())
		{	
			//split the hit fiber into two.
			unsigned int point = E[endId]->nearest_idx(pos[pos.size() -1]);
			
			//split the hit fiber at point into two parts temp[0], temp[1]
			std::vector < stim::fiber <T> > temp
				 = E[endId]->split(point);
			if(temp.size() > 1)
			{
				//add the nearest point in the hit fiber into the hitting fiber.
				pos.push_back(E[endId]->nearest(pos[pos.size()-1]));
	//			pos.insert(pos.end(), E[endId].nearest(pos[pos.size()-1]));
				stim::vec<T> v(E[endId]->radius(point), E[endId]->radius(point));
				radii.push_back(v);

				//split the hit fiber at point into two parts temp[0], temp[1]
				std::vector < stim::fiber <T> > temp
					 = E[endId]->split(point);
				
				//reset the fiber at endId to be a shorter fiber (temp[0]).
				std::vector<stim::vec<T> > ce = temp[0].centerline();
				std::vector<stim::vec<T> > cm = temp[0].centerlinemag();

				//remake the edge, such that the ending point of this edge
				//is the same as before, but the starting edge is new.
				V.push_back(node(ce[ce.size()-1])); //split point.
				int tempNodeId = E[endId]->Nodes[2];
				E[endId] = new edge(ce, cm, E[endId]->Nodes[1], (V.size()-1));
				V[V.size()-1].addEdge(endId);

				//add that other part of the fiber (temp[1])
				//such that it begins with the middle node, and ends with 
				//the last node of the previous fiber.
				ce = temp[1].centerline();
				cm = temp[1].centerlinemag();
				E.push_back(new edge(ce, cm, (V.size()-1), tempNodeId));
				V[V.size()-1].addEdge(E.size()-1);
	//			V[V.size()-1].removeEdge(endId);
				
				
		
				//make the new hitting fiber..
				V.push_back(pos[0]);
				edge *a = new edge(pos,radii,(V.size()-1), (V.size()-2));
				E.push_back(a);
				V[V.size()-1].addEdge(E.size()-1);
				V[V.size()-2].addEdge(E.size()-1);

				//in the end we have added two new nodes and two new edges.
			}
			else {
				stim::vec<T> pz = E[endId]->nearest(pos[0]);
				if((V[ E[endId]->Nodes[1]].getPosition() - pz).len() <
					(V[E[endId]->Nodes[2]].getPosition() - pz).len())
				{
					unsigned int point = E[endId]->nearest_idx(pos[0]);
					pos.insert(pos.begin(), E[endId]->nearest(pos[0]));
					stim::vec<T> v(E[endId]->radius(point), E[endId]->radius(point));
					radii.insert(radii.begin(), v);
					V.push_back(node(pos[pos.size()-1]));
					edge *a = new edge(pos, radii, E[endId]->Nodes[1], (V.size()-1));
					E.push_back(a);
					
					V[V.size()-1].addEdge(E.size()-1);
					V[ E[endId]->Nodes[1]].addEdge(E.size()-1);

				} 
				else
				{
					unsigned int point = E[endId]->nearest_idx(pos[0]);
					pos.insert(pos.begin(), E[endId]->nearest(pos[0]));
					stim::vec<T> v(E[endId]->radius(point), E[endId]->radius(point));
					radii.insert(radii.begin(), v);
					V.push_back(node(pos[pos.size()-1]));
					edge *a = new edge(pos, radii, E[endId]->Nodes[2], (V.size()-1));
					E.push_back(a);
					
					V[V.size()-1].addEdge(E.size()-1);
					V[ E[endId]->Nodes[2]].addEdge(E.size()-1);
				}
			}
		}

		if(startId != -1 && endId != -1 && endId < sizeE())
		{
			//case where the edge is starting with a fiber, and ends in one.
			
			//split the fiber behind you into two.
			unsigned int point = E[startId]->nearest_idx(pos[0]);
//			std::cout << "in merge to both" << std::endl;

			//split the hit fiber at point two parts temp[0], temp[1]
			std::vector < stim::fiber <T> > temp = E[startId]->split(point);
			if(temp.size() > 1)
			{
				//extend the previous fiber (guaranteed to be added last) by one
				//and add the
				pos = E[E.size()-1]->centerline();
				radii = E[E.size()-1]->centerlinemag();
				pos.insert(pos.begin(), E[startId]->nearest(pos[0]));
				stim::vec<T> v(E[startId]->radius(point), E[startId]->radius(point));
				radii.insert(radii.begin(), v);
				V.erase(V.end());
				V.push_back(node(pos[0]));


				//something weird is going on here.
				E[E.size()-1] = new edge(pos, radii, (V.size()-2), (V.size()-1));
				V[V.size()-1].addEdge(E.size()-1);

				//reset the fiber at the endId to be a shorter fiber(temp[0]).
				std::vector<stim::vec<T> > ce = temp[0].centerline();
				std::vector<stim::vec<T> > cm = temp[0].centerlinemag();

//				std::cout << ce.size() << std::endl;
				//remake the edge, such that the starting point of this edge
				//is the same as before, but the end edge is new.
				int tempNodeId = E[startId]->Nodes[1];
				E[startId] = new edge(ce, cm, (V.size()-1), E[startId]->Nodes[2]);
				V[V.size()-1].addEdge(startId);

				//add the other part of the fiber (temp[1])
				ce = temp[1].centerline();
				cm = temp[1].centerlinemag();
				E.push_back(new edge(ce, cm,tempNodeId, (V.size()-1)));
				V[V.size()-1].addEdge(E.size()-1);
				V[tempNodeId].removeEdge(startId);
				V[tempNodeId].addEdge(E.size()-1);
//				V[V.size()-1].removeEdge(startId);
			}
			else {
				stim::vec<T> pz = E[endId]->nearest(pos[0]);
				if((V[ E[endId]->Nodes[1]].getPosition() - pz).len() <
					(V[E[endId]->Nodes[2]].getPosition() - pz).len())
				{
					unsigned int point = E[endId]->nearest_idx(pos[0]);
					pos.insert(pos.begin(), E[endId]->nearest(pos[0]));
					stim::vec<T> v(E[endId]->radius(point), E[endId]->radius(point));
					radii.insert(radii.begin(), v);
					V.push_back(node(pos[pos.size()-1]));
					edge *a = new edge(pos, radii, E[endId]->Nodes[1], (V.size()-1));
					E.push_back(a);
					
					V[V.size()-1].addEdge(E.size()-1);
					V[ E[endId]->Nodes[1]].addEdge(E.size()-1);

				} 
				else
				{
					unsigned int point = E[endId]->nearest_idx(pos[0]);
					pos.insert(pos.begin(), E[endId]->nearest(pos[0]));
					stim::vec<T> v(E[endId]->radius(point), E[endId]->radius(point));
					radii.insert(radii.begin(), v);
					V.push_back(node(pos[pos.size()-1]));
					edge *a = new edge(pos, radii, E[endId]->Nodes[2], (V.size()-1));
					E.push_back(a);
					
					V[V.size()-1].addEdge(E.size()-1);
					V[ E[endId]->Nodes[2]].addEdge(E.size()-1);
				}
			}
		}
		
	}

	///@param pos: std::vector of stim vecs with the positions of the point in the fiber.
	///@param radii: std::vector of floats with the radii of the fiber at positions in pos.
	///returns the forest as a std::string. For testing only.
	std::string
	edges_to_str()
	{
		std::stringstream ss;
		for(unsigned int i = 0; i < E.size(); i++)
		{
			ss << i << ": " << E[i]->str() << std::endl;
		}
		return(ss.str());
	}
	// total number of points on all edges!=fibers in the network
	int
	totalEdges()
	{
		int totalEdgeVertices=0;int N=0;
		for (unsigned int i=0; i < sizeE(); i ++)
		{// FOR N points on the fiber N-1 edges are possible
			N = E[i]->n_pts();
			totalEdgeVertices = totalEdgeVertices + N- 1;
		}
		return totalEdgeVertices;
	}
	// sum of all the fiber lengths
	float
	lengthOfNetwork()
	{
		float networkLength=0;float N=0;
		for (unsigned int i=0; i < sizeE(); i ++)
		{
			N = E[i]->length();
			networkLength = networkLength + N;
		}
		return networkLength;
	}
	///@param i: index of the required fiber.
	///Returns an std::vector with the centerline of the ith fiber in the edgelist.
	std::vector< stim::vec<T> >
	getEdgeCenterLine(int i)
	{
		std::vector < stim::vec<T> > a;
		return E[i]->centerline();
	}

	///@param i: index of the required fiber.
	///Returns an std::vector with the centerline radii of the ith fiber in the edgelist.
	std::vector< stim::vec<T> >
	getEdgeCenterLineMag(int i)
	{
		std::vector < stim::vec<T> > a;
		return E[i]->centerlinemag();
	}
	
	///@param i: index of the required fibers nodes..
	///Returns an std::string with the start and end points of this node..
	std::string
	nodes_to_str(int i)
	{
		std::stringstream ss;
		ss << "[from Node " << E[i] -> Nodes[1] << " to " << E[i] -> Nodes[2] << "]";
		return ss.str();
	}
	//load an obj file into a network class
	stim::network<T>
	objToNetwork(stim::obj<T> objInput)
	{
		stim::network<T> nwc;
		//network network2;
		// function to convert an obj file loaded using stim/visualization/obj.h
		// to a 3D network class using network.h methods.
		std::vector< stim::vec<T> > fiberPositions; //initialize positions on the fiber to be added to network
		objInput.getLine(1, fiberPositions); int numDim = fiberPositions[0].size();//find dimensions of the vertices in obj file
		// to verify if the nodes are already pushed into node list
		std::vector<bool> validList;
		validList.assign(objInput.numV(), false);
		// go through each of the lines "l followed by indices in obj and add all start and end vertices as the nodes
		// using addNode function for adding nodes
		// and the line as an edge(belongs to fiber class) using addEdge function
		std::vector<unsigned> vertexIndices(objInput.numV());
		std::vector< stim::vec<T> > vertices = objInput.getAllV(vertexIndices);
		nwc.addVertices(vertices);
		for (unsigned i =1; i< objInput.numL() + 1; i++) 
			{
			    // edges/fibers could be added to the network class
				std::vector< stim::vec<T> > fiberPositions;

			    objInput.getLine(i, fiberPositions);
			    // finding size to allocate radii
			    int numPointsOnFiber = fiberPositions.size();
			    // to extend it to a 3D postion if it is a 1D/2D vertex in space
				std::vector< stim::vec<T> > fiberPositions1(numPointsOnFiber);
			   // 2D case append and assign the last dimension to zero
			   if (numDim == 2)
			   {//  2D case append and assign the last dimension to zero repeating it for all the points on fiber
				   for (int j = 0; j < numPointsOnFiber; j ++)
				   {
					   fiberPositions1[j][numDim - 2] = fiberPositions[j][0];
					   fiberPositions1[j][numDim -1 ] = fiberPositions[j][1];
					   fiberPositions1[j][numDim] = 0;
				   }
			   }
			   // 1D case append and assign last two dimensions to zero
			   else if (numDim == 1)
			   {
				   for (int j = 0; j < numPointsOnFiber; j ++)
				   {
					   fiberPositions1[j][numDim - 2] = fiberPositions[j][0];
					   fiberPositions1[j][numDim -1 ] = 0;
					   fiberPositions1[j][numDim] = 0;
				   }
			   }
			   else
			   {
				   fiberPositions1 = fiberPositions;
			   }
		   		// then add edge 
			    //edge* a = new edge(fiberPositions1,radii); 
				//std::cout<<"here"<<std::endl;
			    //E.push_back(a);
			    std::vector<stim::vec<T> > newPointList;
			    newPointList  = Resample(fiberPositions1);
				int numPointsOnnewFiber = newPointList.size();
				nwc.addVerticesAfterSamplimg(newPointList);
				std::vector<T> radii(numPointsOnnewFiber); // allocating radii to zero
				nwc.addEdge(newPointList, radii);
				// add nodes
				stim::vec<T> v0(3);stim::vec<T> vN(3);
				int endId = numPointsOnnewFiber -1;
				v0[0] = newPointList[0][0];v0[1] = newPointList[0][1];v0[2] = newPointList[0][2];
				vN[0] = newPointList[endId][0];vN[1] = newPointList[endId][1];vN[2] = newPointList[endId][2];
				// VISITED INDEXES OF the nodes are set to true
				if(!validList[objInput.getIndex(vertices, v0)])
				{
					//V.push_back(node(v0));
					nwc.addNode(v0);
					validList[objInput.getIndex(vertices, v0)] = true;
				}
				if(!validList[objInput.getIndex(vertices, vN)])
				{
					//V.push_back(node(vN));
					nwc.addNode(vN);
					validList[objInput.getIndex(vertices, vN)] = true;
				}
			}
		return nwc;
		}
	// convert ground and test case to network ,kd trees
	boost::tuple<  ANNkd_tree*, ANNkd_tree*, stim::network<T>, stim::network<T> > 
	LoadNetworks(std::string gold_filename, std::string test_filename)
	{
		ANNkd_tree* kdTree1;ANNkd_tree* kdTree2;
		using namespace stim;
		network network1;network network2;
		stim::obj<T> objFile1(gold_filename);
		network1 = objToNetwork(objFile1);
		kdTree1 = generateKdTreeFromNetwork(network1);
		std::cout<<"Gold Standard:network and kdtree generated"<<std::endl;
		stim::obj<T> objFile2(test_filename);
		network2 = objToNetwork(objFile2);
		kdTree2 = generateKdTreeFromNetwork(network2);
		std::cout<<"Test Network:network and kdtree generated"<<std::endl;
		return boost::make_tuple(kdTree1, kdTree2, network1, network2);
		std::cout<<"out of this loop"<<std::endl;
	}
	//distance between two points
	double dist(std::vector<T> p0, std::vector<T> p1)
	{

		double sum = 0;
		for(unsigned int d = 0; d < 3; d++)
			sum += p0[d] * p1[d];
		return sqrt(sum);
	}
	// sum of elements in a vector
	double sum(stim::vec<T> metricList)
	{
		float sumMetricList = 0;
		for (unsigned int count=0; count<metricList.size(); count++)
		{ sumMetricList += metricList[count];}
		return sumMetricList;
	}
	//generate a kd tree from network
	ANNkd_tree* 
	generateKdTreeFromNetwork(stim::network<T> network)			//kd-tree stores all points in the fiber for fast searching
	{       
		std::cout<<"kd trees generated"<<std::endl;
		ANNkd_tree* kdt;
		double **c;						//centerline (array of double pointers)
		float* r;						   // array of fiber radii
		unsigned int n_data = network.sizeAfterSamplingV(); //set the number of points
		c = (double**) malloc(sizeof(double*) * n_data);  //allocate the array pointer
		for(unsigned int i = 0; i < n_data; i++)		 //allocate space for each point
			{c[i] = (double*) malloc(sizeof(double) * 3);}	
		r = (float*) malloc(sizeof(float) * n_data);			//allocate space for the radii
		//stim::vec<T> node;
		for(unsigned int i = 0; i < n_data; i++)
			{
				//node = network.V[i].getPosition();
				//convert_to_double
				for(unsigned int d = 0; d < 3; d++){		//for each dimension
					c[i][d] = double(network.allVerticesAfterSampling[i][d]);				//copy the coordinate
				}
			r[i] = r[i];						//copy the radius
			}
		ANNpointArray pts = (ANNpointArray)c;           //create an array of data points
		kdt = new ANNkd_tree(pts, n_data, 3);			//build a KD tree
		return kdt;
	}

	///@param pos: std::vector of stim vecs with the positions of the point in the fiber.
	///@param radii: std::vector of floats with the radii of the fiber at positions in pos.
	///adds an edge from two std::vectors with positions and radii.
	void
	addEdge(std::vector< stim::vec<T> > pos, std::vector<T> radii)
	{
		edge *a = new edge(pos,radii);
		E.push_back(a);
	}
	void
	addNode(stim::vec<T> nodes)
	{
		V.push_back(nodes);
	}
	void
	addVertices(std::vector< stim::vec<T> > vertices)
	{
		for (unsigned int i=0; i < vertices.size(); i ++)
		{
			allVerticesList.push_back(vertices[i]);
		}
	}
	void
	addVerticesAfterSamplimg(std::vector< stim::vec<T> > vertices)
	{
		for (unsigned int i=0; i < vertices.size(); i ++)
		{
			allVerticesAfterSampling.push_back(vertices[i]);
		}
	}
	// gaussian function
	float gaussianFunction(float x, float std=25)
	{
		float evaluate = 1.0f - ((exp(-x/(2*std*std))));
		return evaluate;
	}
	// compares point on a network to a kd tree for two skeletons
	boost::tuple< float, float > 
	compareSkeletons(boost::tuple<  ANNkd_tree*, ANNkd_tree*, stim::network<float>, stim::network<float> > networkKdtree)
	{
		float gFPR, gFNR;
		gFPR = CompareNetKD(networkKdtree.get<0>(), networkKdtree.get<3>());
		gFNR = CompareNetKD(networkKdtree.get<1>(), networkKdtree.get<2>());
		return boost::make_tuple(gFPR, gFNR);
	}
	// gaussian distance of points on network to Kdtree
	float 
	CompareNetKD(ANNkd_tree* kdTreeGroundtruth, stim::network<T> networkTruth)
	{
		std::vector< stim::vec<T> > fiberPoints;
		float netmetsMetric=0;
		double eps = 0; // error bound
		ANNdistArray dists1;dists1 = new ANNdist[1]; // near neighbor distances
		ANNdistArray dists2; dists2 = new ANNdist[1];// near neighbor distances
		ANNidxArray nnIdx1; nnIdx1 = new ANNidx[1]; // near neighbor indices // allocate near neigh indices
		ANNidxArray nnIdx2; nnIdx2 = new ANNidx[1]; // near neighbor indices // allocate near neigh indices
		int N; int numQueryPoints = networkTruth.totalEdges();
		float totalNetworkLength = networkTruth.lengthOfNetwork();
		stim::vec<float> fiberMetric(networkTruth.sizeE());
		//for each fiber
		for (unsigned int i=0; i < networkTruth.sizeE(); i ++)
		{
			std::vector<T> p1(3); std::vector<T> p2(3);//temporary variables to store point positions
			fiberPoints = networkTruth.E[i]->centerline();
			N = networkTruth.E[i]->n_pts();
			stim::vec<float> segmentMetric(N-1);
			// for each segment in the fiber
			for (unsigned int j = 0; j < N - 1; j++)
			{
				ANNpoint queryPt1; queryPt1 = annAllocPt(3);
				ANNpoint queryPt2; queryPt2 = annAllocPt(3);
				//for each dimension of the point
				for(unsigned int d = 0; d < 3; d++)
				{ 
					queryPt1[d] = double(fiberPoints[j][d]);
					p1[d] = double(fiberPoints[j][d]);
					queryPt2[d] = double(fiberPoints[j + 1][d]);
					p2[d] = double(fiberPoints[j + 1][d]);// for each dimension
				}
				kdTreeGroundtruth->annkSearch( queryPt1, 1, nnIdx1, dists1, eps); // error bound
				kdTreeGroundtruth->annkSearch( queryPt2, 1, nnIdx2, dists2, eps); // error bound
				std::cout<<float(dists1[0])<<"dist1-----"<<float(dists2[0])<<"dist2---"<<std::endl;
				float dist1 = gaussianFunction(float(dists1[0]));float dist2 = gaussianFunction(float(dists2[0]));
				segmentMetric[j] = (((dist1 + dist2) / 2) * dist(p1, p2)) ;
			}
			fiberMetric[i] = sum(segmentMetric)/;
		}
		netmetsMetric =  sum(fiberMetric)/totalNetworkLength ;
		return netmetsMetric;	
	}
    //resample a fiber in the network
	std::vector<stim::vec<T> > Resample(std::vector< stim::vec<T> > fiberPositions, float spacing=25.0)
	{
		std::vector<T> p1(3), p2(3), v(3);
		stim::vec<T> p(3);
		std::vector<stim::vec<T> > newPointList;
		for(unsigned int f=0; f<fiberPositions.size()-1; f++)
		{
			T lengthSegment = dist(p1,p2);
			for(int d=0; d<3;d++)
			{
				p1[d] = T(fiberPositions[f][d]);
			    p2[d] = T(fiberPositions[f + 1][d]);
			}// for each dimension
			if( lengthSegment >= spacing )
				{	for (int dim=0; dim<3;dim++) //find the direction of travel
						{v[dim] = p2[dim] - p1[dim];}
					//set the step size to the voxel size
					T step;
					for(step=0.0; step<lengthSegment; step+=spacing)
					{
						for(int i=0; i<3;i++)
							{p[i] = p1[i] + v[i]*(step/lengthSegment);}
						newPointList.push_back(p);
					}
				}
			else 
				newPointList = fiberPositions;
			}
		return newPointList;
	}
	// modulus of a vector
	T
	Length(std::vector<T> v)
	{
	 	T sum=0;
        for (int i=0; i<v.size(); i++)
			sum += v[i] * v[i];
	    return sum;
	}
	// function to remove characters from a 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 the network graph to obj
	void
	to_obj()
	{
		std::ofstream ofs;
		ofs.open("Graph.obj", std::ofstream::out | std::ofstream::app);
		int num = allVerticesList.size();
		std::string *strArray = new std::string[num];
		for(unsigned int i = 0; i < allVerticesList.size(); i++)
		{
			std::string str;
			str = allVerticesList[i].str();
			removeCharsFromString(str, "[],");
			ofs << "v " << str << "\n";
			removeCharsFromString(str," ");
			strArray[i] = str;
		}
		for(unsigned int i = 0; i < E.size(); i++)
		{
			 std::string str;
			 str = E[i]->strObj(strArray, num);
			 //removeCharsFromString(str,"0");
             ofs << "l " << str << "\n";
		} 
		ofs.close();
	}	
		///exports the graph.
	void
	to_csv()
	{
		std::ofstream ofs;
		ofs.open("Graph.csv", std::ofstream::out | std::ofstream::app);
		std::cout<<"here"<<std::endl;
		for(int i = 0; i < V.size(); i++)
		{
			std::cout<<"here"<<i<<std::endl;
			ofs << V[i].edges_to_str() << "\n";
		}
		ofs.close();
	}
	///exports the graph.
	void
	to_gdf()
	{
		std::ofstream ofs;
		ofs.open("Graph.gdf", std::ofstream::out | std::ofstream::app);
		ofs << "nodedef>name VARCHAR\n";
		for(int i = 0; i < V.size(); i++)
		{
			ofs << i << "\n";
		}
		ofs << "edgedef>Nodes[1] VARCHAR, Nodes[2] VARCHAR, weight INT, length FLOAT, av_radius FLOAT \n";
		for(int i = 0; i < E.size(); i++)
		{
			ofs << E[i]->Nodes[1] << "," << E[i]->Nodes[2] << "," <<E[i]->n_pts()
			 << ","<< E[i]->length() << "," << E[i]->average_radius() << "\n";
		} 
		ofs.close();
	}

};
};
#endif