codebase / stimlib

  #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/structures/kdtree.cuh>

  #include <stim/cuda/cudatools/timer.h>

  #ifdef __CUDACC__
  //device gaussian function
  __device__ float gaussianFunction(float x, float std){ return expf(-x/(2*std*std));}
  
  //compute metric in parallel
  template <typename T>
  __global__ void d_metric(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 d_findedge(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 = 0;

  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 cylinder<T>

  	{
  		public:
  		unsigned 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 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;
  		}
  

  	};
  

  	///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.

  
  			//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();

  			}

  	};
  

  protected:

  	std::vector<edge> E;       //list of edges

  	std::vector<vertex> V;	   //list of vertices.

  	//std::vector<int> NT;	   //stores a list of neuronal type for each point in the centerline (will set value only when centerline is built from swc file)

  
  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();

  	}
  

  	//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

  				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

  		}

  	}
  

  	//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
  

  		//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 = 1; l < S.numL(); l++) {		// for every node starts from second one

  			//NT.push_back(S.node[l].type);

  			stim::centerline<T> c3(2);						// every fiber contains two vertices

  			int p_idx = S.node[l].parent_idx - 1;			// get the parent node loop idx of current node
  			c3[0] = S.node[p_idx].point;					// store the begin vertex

  			c3[1] = S.node[l].point;						// store the end vertex
  		
  			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;
  			radius.push_back(S.node[p_idx].radius);
  			radius.push_back(S.node[l].radius);
  			C3.copy_r(radius);
  
  			edge new_edge(C3);								// contains two points								
  			

  			//get the first and last vertex IDs for the line

  			i[0] = S.node[p_idx].idx;						//get the SWC ID for the start point

  			i[1] = S.node[l].idx;							//get the SWC ID for the end point
  
  			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[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
  		}
  	}
  

  	/// 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

  	}
  

  	// host gaussian function
  	__host__ float gaussianFunction(float x, float std = 25){ return exp(-x/(2*std*std));}  // std default value is 25

  	/// 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];	 
  		}
  	}
  

  	/// 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