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/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>(
  		{
  
  		}

  */

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

  			e.v[0] = v0; e.v[1] = v1;

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

  				e[i] = point;
  				T mag;

  		//				for(int j = 0; j < nmags(); j++)		///future code for mags
  		//				{
  				in.read(reinterpret_cast<char*>(&mag), sizeof(T));	

  				e.set_r(0, mag);
  				std::cout << point.str() << " " << mag << std::endl;

  		//				}

  			}
  			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 edge 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 "incoming edges"
  				out.write(reinterpret_cast<const char*>(&s1), sizeof(unsigned int));	///write the number of "outgoing edges"	
  				out.write(reinterpret_cast<const char*>(&v.e[0][0]), sizeof(unsigned int)*v.e[0].size());	///write the "incoming edges"
  				out.write(reinterpret_cast<const char*>(&v.e[1][0]), sizeof(unsigned int)*v.e[1].size());	///write the "outgoing edges"

  				return out;
  			}
  
  			///operator for reading the edge 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;

  				in.read(reinterpret_cast<char*>(&v.e[0][0]),s[0]*sizeof(unsigned int));		///read the arrays of "incoming edges"
  				in.read(reinterpret_cast<char*>(&v.e[1][0]),s[1]*sizeof(unsigned int));		///read the arrays of "outgoing 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

  				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

  		}

  	}

  	///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, 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;
  		}
  		file.close();