network.h 12.6 KB
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459
#ifndef STIM_NETWORK_H
#define STIM_NETWORK_H

#include <stim/math/vector.h>
#include <stim/visualization/obj.h>
#include <list>
#include <ANN/ANN.h>

namespace stim{

/** This class provides an interface for dealing with biological networks.
 *  It takes the following aspects into account:
 *  	1) Network geometry and centerlines
 *  	2) Network connectivity (a graph structure can be extracted)
 *  	3) Network surface structure (the surface is represented as a triangular mesh and referenced to the centerline)
 */

template<typename T>
class network{

	//-------------------HELPER CLASSES-----------------------
	/// Stores information about a geometric point on the network centerline (including point position and radius)
	//template<typename T>
	class point : public stim::vec<T>{

	public:
		T r;

		point() : stim::vec<T>(){}

		//casting constructor
		point(stim::vec<T> rhs) : stim::vec<T>(rhs){}
	};

	//template<typename T>
	class t_node;
	class fiber;

	//create typedefs for the iterators to simplify the network code
	typedef typename std::list< fiber >::iterator fiber_i;
	typedef typename std::list< t_node >::iterator t_node_i;

	/// Stores information about a single capillary (a length of vessel between two branch or end points)
	//template<typename T>
	class fiber : public std::list< point >{

		using std::list< point >::begin;
		using std::list< point >::end;
		using std::list< point >::size;


	public:
		//std::list< point > P;		//geometric point positions

		typename std::list< t_node >::iterator n[2];				//indices to terminal nodes
		unsigned int id;

	public:

		/// Calculate the length of the fiber and return it.
		T length(){

			point p0, p1;
			T l = 0;				//initialize the length to zero

			//for each point
			typename std::list< point >::iterator i;	//create a point iterator
			for(i = begin(); i != end(); i++){		//for each point in the fiber

				if(i == begin())						//if this is the first point, just store it
					p1 = *i;
				else{									//if this is any other point
					p0 = p1;							//shift p1->p0
					p1 = *i;							//set p1 to the new point
					l += (p1 - p0).len();				//add the length of p1 - p0 to the running sum
				}
			}

			return l;									//return the length
		}

		T radius(T& length){

			point p0, p1;				//temporary variables to store point positions
			T r0, r1;					//temporary variables to store radii at points
			T l, r;						//temporary variable to store the length and average radius of a fiber segment
			T length_sum = 0;			//initialize the length to zero
			T radius_sum = 0;			//initialize the radius sum to zero

			//for each point
			typename std::list< point >::iterator i;	//create a point iterator
			for(i = begin(); i != end(); i++){		//for each point in the fiber

				if(i == begin()){						//if this is the first point, just store it
					p1 = *i;
					r1 = i->r;
				}
				else{									//if this is any other point
					p0 = p1;							//shift p1->p0 and r1->r0
					r0 = r1;
					p1 = *i;							//set p1 to the new point
					r1 = i->r;								//and r1

					l = (p1 - p0).len();				//calculate the length of the p0-p1 segment
					r = (r0 + r1) / 2;					//calculate the average radius of the segment

					radius_sum += r * l;				//add the radius scaled by the length to a running sum
					length_sum += l;					//add the length of p1 - p0 to the running sum
				}
			}

			length = length_sum;						//store the total length

			//if the total length is zero, store a radius of zero
			if(length == 0)
				return 0;
			else
				return radius_sum / length;					//return the average radius of the fiber
		}

		std::vector< stim::vec<T> > geometry(){

			std::vector< stim::vec<T> > result;				//create an array to store the fiber geometry
			result.resize( size() );					//pre-allocate the array

			typename std::list< point >::iterator p;				//create a list iterator
			unsigned int pi = 0;						//create an index into the result array

			//for each geometric point on the fiber centerline
			for(p = begin(); p != end(); p++){
				result[pi] = *p;
				pi++;
			}

			return result;			//return the geometry array

		}

		std::string str(){
			std::stringstream ss;

			//create an iterator for the point list
			typename std::list<point>::iterator i;
			for(i = begin(); i != end(); i++){
				ss<<i->str()<<"  r = "<<i->r<<std::endl;
			}

			return ss.str();
		}
	};

	/// Terminal node for a capillary. This is analogous to a graph vertex and contains a list of edge indices.
	//template<typename T>
	class t_node{

	public:

		unsigned int id;

		//lists of edge indices for capillaries
			//the "in" and "out" just indicate how the geometry is defined:
			//		edges in the "in" list are geometrically oriented such that the terminal node is last
			//		edges in the "out" list are geometrically oriented such that the terminal node is first
		std::list< fiber_i > in;			//edge indices for incoming capillaries
		std::list< fiber_i > out;		//edge indices for outgoing capillaries

		std::string str(){

			std::stringstream ss;

			ss<<id<<": ";						//output the node ID

			//output the IDs for both lists
			typename std::list< fiber_i >::iterator f;

			for(f = in.begin(); f != in.end(); f++){

				if(f != in.begin())
					ss<<", ";
				ss<<(*f)->n[0]->id;
			}

			//if there are nodes in both lists, separate them by a comma
			if(out.size() > 0 && in.size() > 0)
				ss<<", ";

			for(f = out.begin(); f != out.end(); f++){

				if(f != out.begin())
					ss<<", ";
				ss<<(*f)->n[1]->id;
			}


			return ss.str();



		}
	};


//---------------NETWORK CLASS-----------------------------

protected:

	//list of terminal nodes
	std::list<t_node> N;

	//list of fibers
	std::list<fiber> F;

	/// Sets a unique ID for each terminal node and fiber
	void set_names(){

		unsigned int i;

		i = 0;
		for(t_node_i ti = N.begin(); ti != N.end(); ti++)
			ti->id = i++;

		i = 0;
		for(fiber_i fi = F.begin(); fi != F.end(); fi++)
			fi->id = i++;
	}

public:

	std::string str(){


		//assign names to elements of the network
		set_names();

		//create a stringstream for output
		std::stringstream ss;

		//output the nodes
		ss<<"Nodes----------------------------"<<std::endl;
		for(t_node_i i = N.begin(); i != N.end(); i++){
			ss<<i->str()<<std::endl;
		}

		//output the fibers
		ss<<std::endl<<"Fibers---------------------------"<<std::endl;

		T length, radius;
		//output every fiber
		for(fiber_i f = F.begin(); f != F.end(); f++){

			//calculate the length and average radius
			radius = f->radius(length);

			//output the IDs of the terminal nodes
			ss<<f->n[0]->id<<" -- "<<f->n[1]->id<<": length = "<<length<<",  average radius = "<<radius<<std::endl;
		}

		return ss.str();
	}

	/// Load a network from an OBJ object

	/// @param object is the object file to be used as the basis for the network
	void load( stim::obj<T> object){

		//get the number of vertices in the object
		unsigned int nV = object.numV();

		//allocate an array of pointers to nodes, which will be used to preserve connectivity
			//initiate all values to T.end()
		std::vector< t_node_i > node_hash(nV, N.end());

		unsigned int nL = object.numL();			//get the number of lines in the OBJ

		//for each line in the OBJ structure
		for(unsigned int li = 0; li < nL; li++){

			F.push_back(fiber());				//push a new fiber onto the fiber list

			fiber_i f = --(F.end());			//get an iterator to the new fiber

			//----------Handle the terminating nodes for the fiber

			//get the indices of the line vertices
			std::vector< unsigned int > Li = object.getL_Vi(li);
			unsigned int i0 = Li.front() - 1;
			unsigned int i1 = Li.back() - 1;

			//deal with the first end point of the capillary
			if(node_hash[i0] != N.end()){			//if the node has been used before
				(*f).n[0] = node_hash[i0];			//assign the node to the new capillary
				(*node_hash[i0]).out.push_back(f);	//add an out pointer to the existing node
			}
			else{									//otherwise
				N.push_back(t_node());				//create a new node and add it to the node list
				t_node_i t = --(N.end());			//get an iterator to the new node
				node_hash[i0] = t;					//add a pointer to the new node to the hash list
				(*f).n[0] = t;						//add a pointer to the new node to the capillary
				(*t).out.push_back(f);				//add a pointer to the capillary to the new node
			}

			//deal with the last end point of the capillary
			if(node_hash[i1] != N.end()){
				(*f).n[1] = node_hash[i1];
				(*node_hash[i1]).in.push_back(f);
			}
			else{
				N.push_back(t_node());
				t_node_i t = --(N.end());
				node_hash[i1] = t;			//add the new node to the hash list
				(*f).n[1] = t;
				(*t).in.push_back(f);
			}

			//-------------Handle the geometric points for the fiber
			std::vector< vec<T> > L = object.getL_V(li);
			std::vector< vec<T> > R = object.getL_VT(li);

			unsigned int nP = L.size();				//get the number of geometric points in the fiber
			//for each vertex in the fiber
			for(unsigned int pi = 0; pi < nP; pi++){
				point p = (point)L[pi];					//move the geometric coordinates into a point structure
				p.r = R[pi][0];							//store the radius
				f->push_back(p);						//push the point onto the current fiber
			}
		}

	}	//end load()

	/// Returns an array of node positions
	std::vector< stim::vec<T> > get_node_positions(){

		std::vector< stim::vec<T> > result;				//create an array to store the result
		result.resize(N.size());						//set the array size

		t_node_i ni;									//create a terminal node iterator
		unsigned int vi = 0;							//vertex index into the result array

		//for every terminal node
		for(ni = N.begin(); ni != N.end(); ni++){

			//create a vector based on the node position

			//if the number of outgoing nodes is nonzero
			if(ni->out.size() != 0)
				result[vi] = ni->out.front()->front();
			else if(ni->in.size() != 0)
				result[vi] = ni->in.front()->back();

			vi++;										//increment the array index
		}

		//return the resulting array
		return result;
	}

	std::vector< stim::vec<T> > get_fiber_geometry( fiber_i f ){
		return f->geometry();
	}

	/// Generate an OBJ file from the network

	stim::obj<T> obj(){

		//create an OBJ object
		stim::obj<T> object;

		//name the nodes
		set_names();

		//retrieve a list of terminal node positions
		std::vector< stim::vec<T> > node_pos = get_node_positions();

		//add the nodes to the obj file
		object.addV(node_pos);

		//counter for vertex indices in the object class
		unsigned int nP;

		//for each fiber
		fiber_i fi;						//create a fiber iterator
		for(fi = F.begin(); fi != F.end(); fi++){

			//get an array of fiber points
			std::vector< stim::vec<T> > fiber_p = get_fiber_geometry(fi);

			//create a subset of this array
			typename std::vector< stim::vec<T> >::iterator start = fiber_p.begin() + 1;
			typename std::vector< stim::vec<T> >::iterator end = fiber_p.end() - 1;
			typename std::vector< stim::vec<T> > fiber_subset(start, end);

			//add this subset to the geometry object
			nP = object.addV(fiber_subset);

			//create an array to hold vertex indices for a line
			std::vector<unsigned int> line;
			line.resize(fiber_p.size());

			//add the terminal nodes to the line list (make sure to add 1 to make them compatible with the OBJ)
			line[0] = fi->n[0]->id + 1;
			line[line.size() - 1] = fi->n[1]->id + 1;

			//add the intermediate vertex indices to the line array
			for(unsigned int i = 0; i < fiber_subset.size(); i++){
				line[1 + i] = nP + i;
			}

			//add the line list to the object class
			object.addLine(line);

		}

		return object;
	}

	/// This function returns the information necessary for a simple graph-based physical (ex. fluid) simulation.

	/// @param n0 is a array which will contain the list of source nodes
	/// @param n1 is a array which will contain the list of destination nodes
	/// @param length is a array containing the lengths of fibers in the network
	/// @param radius is a array containing the average radii of fibers in the network
	void build_simgraph(std::vector<unsigned int>& n0, std::vector<unsigned int>& n1, std::vector<T>& length, std::vector<T>& radius){

		//determine the number of fibers in the network
		unsigned int nF = F.size();

		//allocate the necessary space to store the fiber information
		n0.resize(nF);
		n1.resize(nF);
		length.resize(nF);
		radius.resize(nF);

		//assign names (identifiers) to the network components
		set_names();

		//fill the arrays
		unsigned int i = 0;
		T l, r;
		for(fiber_i f = F.begin(); f != F.end(); f++){
			n0[i] = f->n[0]->id;	//get the identifiers for the first and second nodes for the current fiber
			n1[i] = f->n[1]->id;

			r = f->radius(l);		//get the length and radius of the capillary (calculated at the same time)

			radius[i] = r;			//store the radius in the output array
			length[i] = l;			//store the length in the output array

			i++;					//increment the array index
		}


	}

};

};	//end namespace stim


#endif