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

  	}
  

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

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

  		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

  	}
  

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