Commit 2f365fd55d25eb9891eae475442e6e256ef02ed3

Authored by David Mayerich
2 parents 5d272e9b e06b2e0b

Merge branch 'pranathi_stimlib'

stim/biomodels/fiber.h 0 → 100644
  1 +#ifndef STIM_FIBER_H
  2 +#define STIM_FIBER_H
  3 +
  4 +#include <vector>
  5 +#include <ANN/ANN.h>
  6 +
  7 +namespace stim{
  8 +
  9 +/** This class stores information about a single fiber represented as a set of geometric points
  10 + * between two branch or end points. This class is used as a fundamental component of the stim::network
  11 + * class to describe an interconnected (often biological) network.
  12 + */
  13 +template<typename T>
  14 +class fiber{
  15 +
  16 +protected:
  17 + unsigned int N; //number of points in the fiber
  18 + double **c; //centerline (array of double pointers)
  19 +
  20 + T* r; // array of fiber radii
  21 + ANNkd_tree* kdt; //kd-tree stores all points in the fiber for fast searching
  22 +
  23 + /// Initialize an empty fiber
  24 + void init()
  25 + {
  26 + kdt = NULL;
  27 + c=NULL;
  28 + r=NULL;
  29 + N=0;
  30 + }
  31 +
  32 + /// Initialize a fiber with N centerline points (all located at [0, 0, 0] with radius 0)
  33 + void init(unsigned int n)
  34 + {
  35 +
  36 + N = n; //set the number of points
  37 + kdt = NULL;
  38 + c = (double**) malloc(sizeof(double*) * N); //allocate the array pointer
  39 +
  40 + for(unsigned int i = 0; i < N; i++) //allocate space for each point
  41 + c[i] = (double*) malloc(sizeof(double) * 3);
  42 +
  43 + r = (T*) malloc(sizeof(T) * N); //allocate space for the radii
  44 + }
  45 +
  46 + /// Copies an existing fiber to the current fiber
  47 +
  48 + /// @param cpy stores the new copy of the fiber
  49 + void copy( const stim::fiber<T>& cpy ){
  50 +
  51 + ///allocate space for the new fiber
  52 + init(cpy.N);
  53 +
  54 + ///copy the points
  55 + for(unsigned int i = 0; i < N; i++){
  56 + for(unsigned int d = 0; d < 3; d++) //for each dimension
  57 + c[i][d] = cpy.c[i][d]; //copy the coordinate
  58 +
  59 + r[i] = cpy.r[i]; //copy the radius
  60 + }
  61 +
  62 + gen_kdtree(); //generate the kd tree for the new fiber
  63 + }
  64 +
  65 + /// generate a KD tree for points on fiber
  66 + void gen_kdtree()
  67 + {
  68 + int n_data = N; //create an array of data points
  69 + ANNpointArray pts = (ANNpointArray)c; //cast the centerline list to an ANNpointArray
  70 + kdt = new ANNkd_tree(pts, n_data, 3); //build a KD tree
  71 + }
  72 +
  73 + /// find distance between two points
  74 + double dist(double* p0, double* p1){
  75 +
  76 + double sum = 0; // initialize variables
  77 + float v;
  78 + for(unsigned int d = 0; d < 3; d++)
  79 + {
  80 + v = p1[d] - p0[d];
  81 + sum +=v * v;
  82 +
  83 + }
  84 + return sqrt(sum);
  85 + }
  86 +
  87 + /// This function retreives the index for the fiber point closest to q
  88 +
  89 + /// @param q is a reference point used to find the closest point on the fiber center line
  90 + unsigned int ann( stim::vec<double> q ){
  91 +
  92 + ANNidxArray idx = new ANNidx[1]; //variable used to hold the nearest point
  93 + ANNdistArray sq_dist = new ANNdist[1]; //variable used to hold the squared distance to the nearest point
  94 +
  95 + kdt->annkSearch(q.data(), 1, idx, sq_dist); //search the KD tree for the nearest neighbor
  96 +
  97 + return *idx;
  98 + }
  99 +
  100 + /// Returns a stim::vec representing the point at index i
  101 +
  102 + /// @param i is an index of the desired centerline point
  103 + stim::vec<T> get_vec(unsigned i){
  104 + stim::vec<T> r;
  105 + r.resize(3);
  106 + r[0] = c[i][0];
  107 + r[1] = c[i][1];
  108 + r[2] = c[i][2];
  109 +
  110 + return r;
  111 + }
  112 +
  113 +
  114 +public:
  115 +
  116 + fiber(){
  117 + init();
  118 + }
  119 +
  120 + /// Copy constructor
  121 + fiber(const stim::fiber<T> &obj){
  122 +
  123 + copy(obj);
  124 +
  125 + }
  126 +
  127 + //temp constructor for graph visualization
  128 + fiber(int n)
  129 + {
  130 + init(n);
  131 + }
  132 +
  133 + /// Constructor takes a list of stim::vec points, the radius at each point is set to zero
  134 + fiber(std::vector< stim::vec<T> > p){
  135 + init(p.size()); //initialize the fiber
  136 +
  137 + //for each point, set the centerline position and radius
  138 + for(unsigned int i = 0; i < N; i++){
  139 +
  140 + //set the centerline position
  141 + for(unsigned int d = 0; d < 3; d++)
  142 + c[i][d] = (double) p[i][d];
  143 +
  144 + //set the radius
  145 + r[i] = 0;
  146 + }
  147 +
  148 + //generate a kd tree
  149 + gen_kdtree();
  150 + }
  151 +
  152 + /// constructor takes a list of points and radii
  153 + fiber(std::vector< stim::vec< T > > pos, std::vector< T > radii){
  154 + init(pos.size()); //initialize the fiber
  155 +
  156 + //for each point, set the centerline position and radius
  157 + for(unsigned int i = 0; i < N; i++){
  158 +
  159 + //set the centerline position
  160 + for(unsigned int d = 0; d < 3; d++)
  161 + c[i][d] = (double) pos[i][d];
  162 +
  163 + //set the radius
  164 + r[i] = radii[i];
  165 + }
  166 +
  167 + //generate a kd tree
  168 + gen_kdtree();
  169 + }
  170 +
  171 + /// constructor takes an array of points and radii
  172 + // this function is used when the radii are represented as a stim::vec,
  173 + // since this may be easier when importing OBJs
  174 + fiber(std::vector< stim::vec<T> > pos, std::vector< stim::vec<T> > radii){
  175 +
  176 + init(pos.size());
  177 +
  178 + //for each point, set the position and radius
  179 + for(unsigned int i = 0; i < N; i++){
  180 + //at(i) = (double*)malloc(sizeof(double) * 3);
  181 + for(unsigned int d = 0; d < 3; d++)
  182 + c[i][d] = (double) pos[i][d];
  183 +
  184 + r[i] = radii[i][(unsigned int)0];
  185 + }
  186 +
  187 + gen_kdtree();
  188 + }
  189 +
  190 + /// Assignment operation
  191 + fiber& operator=(const fiber &rhs){
  192 +
  193 + if(this == &rhs) return *this; //test for and handle self-assignment
  194 +
  195 + copy(rhs);
  196 + }
  197 +
  198 + /// Calculate the length of the fiber and return it.
  199 + double length(){
  200 +
  201 + double* p0;
  202 + double *p1;
  203 + double l = 0; //initialize the length to zero
  204 +
  205 + //for each point
  206 + //typename std::list< point<T> >::iterator i; //create a point iterator
  207 + for(unsigned int i = 0; i < N; i++){ //for each point in the fiber
  208 +
  209 + if(i == 0) //if this is the first point, just store it
  210 + p1 = c[0];
  211 + else{ //if this is any other point
  212 + p0 = p1; //shift p1->p0
  213 + p1 = c[i]; //set p1 to the new point
  214 + l += dist(p0, p1); //add the length of p1 - p0 to the running sum
  215 + }
  216 + }
  217 +
  218 + return l; //return the length
  219 + }
  220 +
  221 + /// Calculates the length and average radius of the fiber
  222 +
  223 + /// @param length is filled with the fiber length
  224 + T radius(T& length){
  225 +
  226 + double* p0; //temporary variables to store point positions
  227 + double* p1;
  228 + T r0, r1; //temporary variables to store radii at points
  229 + double l;
  230 + T r_mean; //temporary variable to store the length and average radius of a fiber segment
  231 + double length_sum = 0; //initialize the length to zero
  232 + T radius_sum = 0; //initialize the radius sum to zero
  233 +
  234 + //for each point
  235 + //typename std::list< point<T> >::iterator i; //create a point iterator
  236 + for(unsigned int i = 0; i < N; i++){ //for each point in the fiber
  237 +
  238 + if(i == 0){ //if this is the first point, just store it
  239 + p1 = c[0];
  240 + r1 = r[0];
  241 + }
  242 + else{ //if this is any other point
  243 + p0 = p1; //shift p1->p0 and r1->r0
  244 + r0 = r1;
  245 + p1 = c[i]; //set p1 to the new point
  246 + r1 = r[i];
  247 +
  248 + l = dist(p0, p1); //calculate the length of the p0-p1 segment
  249 + r_mean = (r0 + r1) / 2; //calculate the average radius of the segment
  250 +
  251 + radius_sum += r_mean * (T) l; //add the radius scaled by the length to a running sum
  252 + length_sum += l; //add the length of p1 - p0 to the running sum
  253 + }
  254 + }
  255 +
  256 + length = length_sum; //store the total length
  257 +
  258 + //if the total length is zero, store a radius of zero
  259 + if(length == 0)
  260 + return 0;
  261 + else
  262 + return (T)(radius_sum / length); //return the average radius of the fiber
  263 + }
  264 + T average_radius()
  265 + {
  266 + T r_sum = 0.;
  267 + for(unsigned int i = 0; i < N; i++)
  268 + {
  269 + r_sum = r_sum + r[i];
  270 + }
  271 + return r_sum/((T) N);
  272 + }
  273 +
  274 + /// Calculates the average radius of the fiber
  275 + T radius(){
  276 + T length;
  277 + return radius(length);
  278 + }
  279 +
  280 + /// Returns the radius at index idx.
  281 + T radius(int idx){
  282 + return r[idx];
  283 + }
  284 +
  285 + /// Return the point on the fiber closest to q
  286 + /// @param q is the query point used to locate the nearest point on the fiber centerline
  287 + stim::vec<T> nearest(stim::vec<T> q){
  288 +
  289 + stim::vec<double> temp( (double) q[0], (double) q[1], (double) q[2]);
  290 +
  291 + unsigned int idx = ann(temp); //determine the index of the nearest neighbor
  292 +
  293 + return stim::vec<T>((T) c[idx][0], (T) c[idx][1], (T) c[idx][2]); //return the nearest centerline point
  294 + }
  295 +
  296 + /// Return the point index on the fiber closest to q
  297 + /// @param q is the query point used to locate the nearest point on the fiber centerline
  298 + unsigned int nearest_idx(stim::vec<T> q){
  299 +
  300 + stim::vec<double> temp((double) q[0], (double) q[1], (double) q[2]);
  301 +
  302 + unsigned int idx = ann(temp); //determine the index of the nearest neighbor
  303 +
  304 + return idx; //return the nearest centerline point index
  305 + }
  306 +
  307 + /// Returns the fiber centerline as an array of stim::vec points
  308 + std::vector< stim::vec<T> > centerline(){
  309 +
  310 + //create an array of stim vectors
  311 + std::vector< stim::vec<T> > pts(N);
  312 +
  313 + //cast each point to a stim::vec, keeping only the position information
  314 + for(unsigned int i = 0; i < N; i++)
  315 + pts[i] = stim::vec<T>((T) c[i][0], (T) c[i][1], (T) c[i][2]);
  316 +
  317 + //return the centerline array
  318 + return pts;
  319 + }
  320 +
  321 + /// Returns the fiber centerline magnitudes as an array of stim::vec points
  322 + std::vector< stim::vec<T> > centerlinemag(){
  323 +
  324 + //create an array of stim vectors
  325 + std::vector< stim::vec<T> > pts(N);
  326 +
  327 + //cast each point to a stim::vec, keeping only the position information
  328 + for(unsigned int i = 0; i < N; i++)
  329 + pts[i] = stim::vec<T>(r[i], r[i]);;
  330 +
  331 + //return the centerline array
  332 + return pts;
  333 + }
  334 +
  335 + /// Split the fiber at the specified index. If the index is an end point, only one fiber is returned
  336 + std::vector< stim::fiber<T> > split(unsigned int idx){
  337 +
  338 + std::vector< stim::fiber<T> > fl; //create an array to store up to two fibers
  339 +
  340 + //if the index is an end point, only the existing fiber is returned
  341 + if(idx == 0 || idx == N-1){
  342 + fl.resize(1); //set the size of the fiber to 1
  343 + fl[0] = *this; //copy the current fiber
  344 + }
  345 +
  346 + //if the index is not an end point
  347 + else{
  348 +
  349 + unsigned int N1 = idx + 1; //calculate the size of both fibers
  350 + unsigned int N2 = N - idx;
  351 +
  352 + fl.resize(2); //set the array size to 2
  353 +
  354 + fl[0].init(N1); //set the size of each fiber
  355 + fl[1].init(N2);
  356 +
  357 + //copy both halves of the fiber
  358 + unsigned int i, d;
  359 +
  360 + //first half
  361 + for(i = 0; i < N1; i++){ //for each centerline point
  362 + for(d = 0; d < 3; d++)
  363 + fl[0].c[i][d] = c[i][d]; //copy each coordinate
  364 + fl[0].r[i] = r[i]; //copy the corresponding radius
  365 + }
  366 +
  367 + //second half
  368 + for(i = 0; i < N2; i++){
  369 + for(d = 0; d < 3; d++)
  370 + fl[1].c[i][d] = c[idx + i][d];
  371 + fl[1].r[i] = r[idx + i];
  372 + }
  373 + }
  374 +
  375 + return fl; //return the array
  376 +
  377 + }
  378 +
  379 + /// Calculates the set of fibers resulting from a connection between the current fiber and a fiber f
  380 +
  381 + /// @param f is the fiber that will be connected to the current fiber
  382 + std::vector< stim::fiber<T> > connect( stim::fiber<T> &f, double dist){
  383 +
  384 + double min_dist;
  385 + unsigned int idx0, idx1;
  386 +
  387 + //go through each point in the query fiber, looking for the indices for the closest points
  388 + for(unsigned int i = 0; i < f.n_pts(); i++){
  389 + //Run through all points and find the index with the closest point, then partition the fiber and return two fibers.
  390 +
  391 + }
  392 +
  393 +
  394 +
  395 + }
  396 +
  397 + /// Outputs the fiber as a string
  398 + std::string str(){
  399 + std::stringstream ss;
  400 +
  401 + //create an iterator for the point list
  402 + //typename std::list< point<T> >::iterator i;
  403 + for(unsigned int i = 0; i < N; i++){
  404 + ss<<" [ ";
  405 + for(unsigned int d = 0; d < 3; d++){
  406 + ss<<c[i][d]<<" ";
  407 + }
  408 + ss<<"] r = "<<r[i]<<std::endl;
  409 + }
  410 +
  411 + return ss.str();
  412 + }
  413 + /// Returns the number of centerline points in the fiber
  414 + unsigned int size(){
  415 + return N;
  416 + }
  417 +
  418 +
  419 + /// Bracket operator returns the element at index i
  420 +
  421 + /// @param i is the index of the element to be returned as a stim::vec
  422 + stim::vec<T> operator[](unsigned i){
  423 + return get_vec(i);
  424 + }
  425 +
  426 + /// Back method returns the last point in the fiber
  427 + stim::vec<T> back(){
  428 + return get_vec(N-1);
  429 + }
  430 + ////resample a fiber in the network
  431 + stim::fiber<T> resample(T spacing)
  432 + {
  433 +
  434 +
  435 + std::vector<T> v(3); //v-direction vector of the segment
  436 + stim::vec<T> p(3); //- intermediate point to be added
  437 + stim::vec<T> p1(3); // p1 - starting point of an segment on the fiber,
  438 + stim::vec<T> p2(3); // p2 - ending point,
  439 + double sum=0; //distance summation
  440 + std::vector<stim::vec<T> > fiberPositions = centerline();
  441 + std::vector<stim::vec<T> > newPointList; // initialize list of new resampled points on the fiber
  442 + // for each point on the centerline (skip if it is the last point on centerline)
  443 + //unsigned int N = fiberPositions.size(); // number of points on the fiber
  444 + for(unsigned int f=0; f< N-1; f++)
  445 + {
  446 +
  447 + p1 = fiberPositions[f]; p2 = fiberPositions[f + 1]; v = p2 - p1;
  448 + for(unsigned int d = 0; d < 3; d++){
  449 + sum +=v[d] * v[d];} //length of segment-distance between starting and ending point
  450 +
  451 + T lengthSegment = sqrt(sum); //find Length of the segment as distance between the starting and ending points of the segment
  452 +
  453 + if(lengthSegment >= spacing) // if length of the segment is greater than standard deviation resample
  454 + {
  455 + // repeat resampling until accumulated stepsize is equsl to length of the segment
  456 + for(T step=0.0; step<lengthSegment; step+=spacing)
  457 + {
  458 + // calculate the resampled point by travelling step size in the direction of normalized gradient vector
  459 + for(unsigned int i=0; i<3;i++)
  460 + {
  461 + p[i] = p1[i] + v[i]*(step/lengthSegment);
  462 + } //for each dimension
  463 + // add this resampled points to the new fiber list
  464 + newPointList.push_back(p);
  465 + }
  466 + }
  467 + else // length of the segment is now less than standard deviation, push the ending point of the segment and proceed to the next point in the fiber
  468 + newPointList.push_back(fiberPositions[f+1]);
  469 + }
  470 + newPointList.push_back(fiberPositions[N-1]); //add the last point on the fiber to the new fiber list
  471 + fiber newFiber(newPointList);
  472 + return newFiber;
  473 + }
  474 +
  475 +};
  476 +
  477 +
  478 +
  479 +} //end namespace stim
  480 +
  481 +
  482 +
  483 +#endif
stim/biomodels/network.h
1 #ifndef STIM_NETWORK_H 1 #ifndef STIM_NETWORK_H
2 #define STIM_NETWORK_H 2 #define STIM_NETWORK_H
3 3
  4 +#include <stdlib.h>
  5 +#include <assert.h>
  6 +#include <sstream>
  7 +#include <fstream>
  8 +#include <algorithm>
  9 +#include <string.h>
  10 +#include <math.h>
4 #include <stim/math/vector.h> 11 #include <stim/math/vector.h>
5 #include <stim/visualization/obj.h> 12 #include <stim/visualization/obj.h>
6 -#include <list> 13 +#include <stim/biomodels/fiber.h>
7 #include <ANN/ANN.h> 14 #include <ANN/ANN.h>
  15 +#include <boost/tuple/tuple.hpp>
  16 +
8 17
9 namespace stim{ 18 namespace stim{
  19 +/** This is the a class that interfaces with gl_spider in order to store the currently
  20 + * segmented network. The following data is stored and can be extracted:
  21 + * 1)Network geometry and centerline.
  22 + * 2)Network connectivity (a graph of nodes and edges), reconstructed using ANN library.
  23 +*/
10 24
11 -/** This class provides an interface for dealing with biological networks.  
12 - * It takes the following aspects into account:  
13 - * 1) Network geometry and centerlines  
14 - * 2) Network connectivity (a graph structure can be extracted)  
15 - * 3) Network surface structure (the surface is represented as a triangular mesh and referenced to the centerline)  
16 - */  
17 25
18 template<typename T> 26 template<typename T>
19 class network{ 27 class network{
20 28
21 - //-------------------HELPER CLASSES-----------------------  
22 - /// Stores information about a geometric point on the network centerline (including point position and radius)  
23 - //template<typename T>  
24 - class point : public stim::vec<T>{  
25 -  
26 - public:  
27 - T r;  
28 -  
29 - point() : stim::vec<T>(){}  
30 -  
31 - //casting constructor  
32 - point(stim::vec<T> rhs) : stim::vec<T>(rhs){}  
33 - };  
34 -  
35 - //template<typename T>  
36 - class t_node;  
37 - class fiber;  
38 -  
39 - //create typedefs for the iterators to simplify the network code  
40 - typedef typename std::list< fiber >::iterator fiber_i;  
41 - typedef typename std::list< t_node >::iterator t_node_i;  
42 -  
43 - /// Stores information about a single capillary (a length of vessel between two branch or end points)  
44 - //template<typename T>  
45 - class fiber : public std::list< point >{  
46 -  
47 - using std::list< point >::begin;  
48 - using std::list< point >::end;  
49 - using std::list< point >::size;  
50 -  
51 -  
52 - public:  
53 - //std::list< point > P; //geometric point positions  
54 -  
55 - typename std::list< t_node >::iterator n[2]; //indices to terminal nodes  
56 - unsigned int id;  
57 -  
58 - public:  
59 -  
60 - /// Calculate the length of the fiber and return it.  
61 - T length(){  
62 -  
63 - point p0, p1;  
64 - T l = 0; //initialize the length to zero  
65 -  
66 - //for each point  
67 - typename std::list< point >::iterator i; //create a point iterator  
68 - for(i = begin(); i != end(); i++){ //for each point in the fiber  
69 -  
70 - if(i == begin()) //if this is the first point, just store it  
71 - p1 = *i;  
72 - else{ //if this is any other point  
73 - p0 = p1; //shift p1->p0  
74 - p1 = *i; //set p1 to the new point  
75 - l += (p1 - p0).len(); //add the length of p1 - p0 to the running sum  
76 - }  
77 - }  
78 -  
79 - return l; //return the length  
80 - }  
81 -  
82 - T radius(T& length){  
83 -  
84 - point p0, p1; //temporary variables to store point positions  
85 - T r0, r1; //temporary variables to store radii at points  
86 - T l, r; //temporary variable to store the length and average radius of a fiber segment  
87 - T length_sum = 0; //initialize the length to zero  
88 - T radius_sum = 0; //initialize the radius sum to zero  
89 -  
90 - //for each point  
91 - typename std::list< point >::iterator i; //create a point iterator  
92 - for(i = begin(); i != end(); i++){ //for each point in the fiber  
93 -  
94 - if(i == begin()){ //if this is the first point, just store it  
95 - p1 = *i;  
96 - r1 = i->r;  
97 - }  
98 - else{ //if this is any other point  
99 - p0 = p1; //shift p1->p0 and r1->r0  
100 - r0 = r1;  
101 - p1 = *i; //set p1 to the new point  
102 - r1 = i->r; //and r1  
103 -  
104 - l = (p1 - p0).len(); //calculate the length of the p0-p1 segment  
105 - r = (r0 + r1) / 2; //calculate the average radius of the segment  
106 -  
107 - radius_sum += r * l; //add the radius scaled by the length to a running sum  
108 - length_sum += l; //add the length of p1 - p0 to the running sum  
109 - }  
110 - }  
111 -  
112 - length = length_sum; //store the total length  
113 -  
114 - //if the total length is zero, store a radius of zero  
115 - if(length == 0)  
116 - return 0;  
117 - else  
118 - return radius_sum / length; //return the average radius of the fiber  
119 - }  
120 -  
121 - std::vector< stim::vec<T> > geometry(){  
122 -  
123 - std::vector< stim::vec<T> > result; //create an array to store the fiber geometry  
124 - result.resize( size() ); //pre-allocate the array  
125 -  
126 - typename std::list< point >::iterator p; //create a list iterator  
127 - unsigned int pi = 0; //create an index into the result array  
128 -  
129 - //for each geometric point on the fiber centerline  
130 - for(p = begin(); p != end(); p++){  
131 - result[pi] = *p;  
132 - pi++;  
133 - }  
134 -  
135 - return result; //return the geometry array  
136 - 29 + ///Each edge is a fiber with two nodes.
  30 + ///Each node is an in index to the endpoint of the fiber in the nodes array.
  31 + class edge : public fiber<T>
  32 + {
  33 + public:
  34 + unsigned v[2]; //unique id's designating the starting and ending
  35 + // default constructor
  36 + edge() : fiber<T>(){v[1] = -1; v[0] = -1;}
  37 + /// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
  38 +
  39 + ///@param p is an array of positions in space
  40 + edge(std::vector< stim::vec<T> > p) : fiber<T>(p){}
  41 +
  42 + /// Copy constructor creates an edge from a fiber
  43 + edge(stim::fiber<T> f) : fiber<T>(f) {}
  44 +
  45 + /// Resamples an edge by calling the fiber resampling function
  46 + edge resample(T spacing){
  47 + edge e(fiber<T>::resample(spacing)); //call the fiber->edge constructor
  48 + e.v[0] = v[0]; //copy the vertex data
  49 + e.v[1] = v[1];
  50 +
  51 + return e; //return the new edge
137 } 52 }
138 - 53 +
  54 + /// Output the edge information as a string
139 std::string str(){ 55 std::string str(){
140 std::stringstream ss; 56 std::stringstream ss;
141 -  
142 - //create an iterator for the point list  
143 - typename std::list<point>::iterator i;  
144 - for(i = begin(); i != end(); i++){  
145 - ss<<i->str()<<" r = "<<i->r<<std::endl;  
146 - }  
147 - 57 + ss<<"("<<fiber<T>::N<<")\tl = "<<length()<<"\t"<<v[0]<<"----"<<v[1];
148 return ss.str(); 58 return ss.str();
149 } 59 }
150 - };  
151 -  
152 - /// Terminal node for a capillary. This is analogous to a graph vertex and contains a list of edge indices.  
153 - //template<typename T>  
154 - class t_node{  
155 -  
156 - public:  
157 60
158 - unsigned int id;  
159 -  
160 - //lists of edge indices for capillaries  
161 - //the "in" and "out" just indicate how the geometry is defined:  
162 - // edges in the "in" list are geometrically oriented such that the terminal node is last  
163 - // edges in the "out" list are geometrically oriented such that the terminal node is first  
164 - std::list< fiber_i > in; //edge indices for incoming capillaries  
165 - std::list< fiber_i > out; //edge indices for outgoing capillaries  
166 -  
167 - std::string str(){  
168 -  
169 - std::stringstream ss;  
170 -  
171 - ss<<id<<": "; //output the node ID  
172 -  
173 - //output the IDs for both lists  
174 - typename std::list< fiber_i >::iterator f;  
175 -  
176 - for(f = in.begin(); f != in.end(); f++){  
177 -  
178 - if(f != in.begin())  
179 - ss<<", ";  
180 - ss<<(*f)->n[0]->id;  
181 - }  
182 -  
183 - //if there are nodes in both lists, separate them by a comma  
184 - if(out.size() > 0 && in.size() > 0)  
185 - ss<<", ";  
186 -  
187 - for(f = out.begin(); f != out.end(); f++){ 61 + };
  62 +
  63 + ///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.
  64 + class vertex : public stim::vec<T>
  65 + {
  66 + public:
  67 + //std::vector<unsigned int> edges; //indices of edges connected to this node.
  68 + std::vector<unsigned int> e[2]; //indices of edges going out (e[0]) and coming in (e[1])
  69 + //stim::vec<T> p; //position of this node in physical space.
  70 +
  71 + //constructor takes a stim::vec
  72 + vertex(stim::vec<T> p) : stim::vec<T>(p){}
  73 +
  74 + /// Output the vertex information as a string
  75 + std::string str(){
  76 + std::stringstream ss;
  77 + ss<<"\t(x, y, z) = "<<stim::vec<T>::str();
  78 +
  79 + if(e[0].size() > 0){
  80 + ss<<"\t> ";
  81 + for(unsigned int o = 0; o < e[0].size(); o++)
  82 + ss<<e[0][o]<<" ";
  83 + }
  84 + if(e[1].size() > 0){
  85 + ss<<"\t< ";
  86 + for(unsigned int i = 0; i < e[1].size(); i++)
  87 + ss<<e[1][i]<<" ";
  88 + }
188 89
189 - if(f != out.begin())  
190 - ss<<", ";  
191 - ss<<(*f)->n[1]->id; 90 + return ss.str();
192 } 91 }
193 -  
194 -  
195 - return ss.str();  
196 -  
197 -  
198 -  
199 - } 92 +
200 }; 93 };
201 94
  95 + private:
202 96
203 -//---------------NETWORK CLASS-----------------------------  
204 -  
205 -protected:  
206 -  
207 - //list of terminal nodes  
208 - std::list<t_node> N;  
209 -  
210 - //list of fibers  
211 - std::list<fiber> F;  
212 -  
213 - /// Sets a unique ID for each terminal node and fiber  
214 - void set_names(){  
215 -  
216 - unsigned int i;  
217 -  
218 - i = 0;  
219 - for(t_node_i ti = N.begin(); ti != N.end(); ti++)  
220 - ti->id = i++; 97 + std::vector<edge> E; //list of edges
  98 + std::vector<vertex> V; //list of vertices.
  99 +
  100 + public:
221 101
222 - i = 0;  
223 - for(fiber_i fi = F.begin(); fi != F.end(); fi++)  
224 - fi->id = i++; 102 + ///Returns the number of edges in the network.
  103 + unsigned int edges(){
  104 + return E.size();
225 } 105 }
226 106
227 -public:  
228 -  
229 - std::string str(){  
230 -  
231 -  
232 - //assign names to elements of the network  
233 - set_names();  
234 -  
235 - //create a stringstream for output  
236 - std::stringstream ss;  
237 -  
238 - //output the nodes  
239 - ss<<"Nodes-----------------------------------------------"<<std::endl;  
240 - for(t_node_i i = N.begin(); i != N.end(); i++){  
241 - ss<<i->str()<<std::endl;  
242 - }  
243 -  
244 - //output the fibers  
245 - ss<<std::endl<<"Fibers----------------------------------------------"<<std::endl;  
246 -  
247 - T length, radius;  
248 - //output every fiber  
249 - for(fiber_i f = F.begin(); f != F.end(); f++){  
250 -  
251 - //calculate the length and average radius  
252 - radius = f->radius(length);  
253 -  
254 - //output the IDs of the terminal nodes  
255 - ss<<f->n[0]->id<<" -- "<<f->n[1]->id<<": length = "<<length<<", average radius = "<<radius<<std::endl;  
256 - }  
257 -  
258 - return ss.str(); 107 + ///Returns the number of nodes in the network.
  108 + unsigned int vertices(){
  109 + return V.size();
259 } 110 }
260 111
261 - /// Load a network from an OBJ object 112 + stim::fiber<T> get_fiber(unsigned f){
  113 + return E[f]; //return the specified edge (casting it to a fiber)
  114 + }
262 115
263 - /// @param object is the object file to be used as the basis for the network  
264 - void load( stim::obj<T> object){ 116 + //load a network from an OBJ file
  117 + void load_obj(std::string filename){
265 118
266 - //get the number of vertices in the object  
267 - unsigned int nV = object.numV(); 119 + stim::obj<T> O; //create an OBJ object
  120 + O.load(filename); //load the OBJ file as an object
268 121
269 - //allocate an array of pointers to nodes, which will be used to preserve connectivity  
270 - //initiate all values to T.end()  
271 - std::vector< t_node_i > node_hash(nV, N.end()); 122 + std::vector<unsigned> id2vert; //this list stores the OBJ vertex ID associated with each network vertex
272 123
273 - unsigned int nL = object.numL(); //get the number of lines in the OBJ 124 + unsigned i[2]; //temporary, IDs associated with the first and last points in an OBJ line
274 125
275 - //for each line in the OBJ structure  
276 - for(unsigned int li = 0; li < nL; li++){ 126 + //for each line in the OBJ object
  127 + for(unsigned int l = 1; l <= O.numL(); l++){
277 128
278 - F.push_back(fiber()); //push a new fiber onto the fiber list 129 + std::vector< stim::vec<T> > c; //allocate an array of points for the vessel centerline
  130 + O.getLine(l, c); //get the fiber centerline
279 131
280 - fiber_i f = --(F.end()); //get an iterator to the new fiber 132 + edge new_edge = c; //create an edge from the given centerline
281 133
282 - //----------Handle the terminating nodes for the fiber 134 + //get the first and last vertex IDs for the line
  135 + std::vector< unsigned > id; //create an array to store the centerline point IDs
  136 + O.getLinei(l, id); //get the list of point IDs for the line
  137 + i[0] = id.front(); //get the OBJ ID for the first element of the line
  138 + i[1] = id.back(); //get the OBJ ID for the last element of the line
283 139
284 - //get the indices of the line vertices  
285 - std::vector< unsigned int > Li = object.getL_Vi(li);  
286 - unsigned int i0 = Li.front() - 1;  
287 - unsigned int i1 = Li.back() - 1; 140 + std::vector<unsigned>::iterator it; //create an iterator for searching the id2vert array
  141 + unsigned it_idx; //create an integer for the id2vert entry index
288 142
289 - //deal with the first end point of the capillary  
290 - if(node_hash[i0] != N.end()){ //if the node has been used before  
291 - (*f).n[0] = node_hash[i0]; //assign the node to the new capillary  
292 - (*node_hash[i0]).out.push_back(f); //add an out pointer to the existing node 143 + //find out if the nodes for this fiber have already been created
  144 + it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node
  145 + it_idx = std::distance(id2vert.begin(), it);
  146 + if(it == id2vert.end()){ //if i[0] hasn't already been used
  147 + vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position
  148 + new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing
  149 + new_edge.v[0] = V.size(); //add the new vertex to the edge
  150 + V.push_back(new_vertex); //add the new vertex to the vertex list
  151 + id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list
293 } 152 }
294 - else{ //otherwise  
295 - N.push_back(t_node()); //create a new node and add it to the node list  
296 - t_node_i t = --(N.end()); //get an iterator to the new node  
297 - node_hash[i0] = t; //add a pointer to the new node to the hash list  
298 - (*f).n[0] = t; //add a pointer to the new node to the capillary  
299 - (*t).out.push_back(f); //add a pointer to the capillary to the new node 153 + else{ //if the vertex already exists
  154 + V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing
  155 + new_edge.v[0] = it_idx;
300 } 156 }
301 157
302 - //deal with the last end point of the capillary  
303 - if(node_hash[i1] != N.end()){  
304 - (*f).n[1] = node_hash[i1];  
305 - (*node_hash[i1]).in.push_back(f); 158 + it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID
  159 + it_idx = std::distance(id2vert.begin(), it);
  160 + if(it == id2vert.end()){ //if i[1] hasn't already been used
  161 + vertex new_vertex = new_edge.back(); //create a new vertex, assign it a position
  162 + new_vertex.e[1].push_back(E.size()); //add the current edge as incoming
  163 + new_edge.v[1] = V.size();
  164 + V.push_back(new_vertex); //add the new vertex to the vertex list
  165 + id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list
306 } 166 }
307 - else{  
308 - N.push_back(t_node());  
309 - t_node_i t = --(N.end());  
310 - node_hash[i1] = t; //add the new node to the hash list  
311 - (*f).n[1] = t;  
312 - (*t).in.push_back(f); 167 + else{ //if the vertex already exists
  168 + V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming
  169 + new_edge.v[1] = it_idx;
313 } 170 }
  171 +
  172 + E.push_back(new_edge); //push the edge to the list
314 173
315 - //-------------Handle the geometric points for the fiber  
316 - std::vector< vec<T> > L = object.getL_V(li);  
317 - std::vector< vec<T> > R = object.getL_VT(li);  
318 -  
319 - unsigned int nP = L.size(); //get the number of geometric points in the fiber  
320 - //for each vertex in the fiber  
321 - for(unsigned int pi = 0; pi < nP; pi++){  
322 - point p = (point)L[pi]; //move the geometric coordinates into a point structure  
323 - p.r = R[pi][0]; //store the radius  
324 - f->push_back(p); //push the point onto the current fiber  
325 - }  
326 } 174 }
  175 + }
327 176
328 - } //end load()  
329 -  
330 - /// Returns an array of node positions  
331 - std::vector< stim::vec<T> > get_node_positions(){  
332 -  
333 - std::vector< stim::vec<T> > result; //create an array to store the result  
334 - result.resize(N.size()); //set the array size  
335 -  
336 - t_node_i ni; //create a terminal node iterator  
337 - unsigned int vi = 0; //vertex index into the result array  
338 -  
339 - //for every terminal node  
340 - for(ni = N.begin(); ni != N.end(); ni++){  
341 -  
342 - //create a vector based on the node position  
343 -  
344 - //if the number of outgoing nodes is nonzero  
345 - if(ni->out.size() != 0)  
346 - result[vi] = ni->out.front()->front();  
347 - else if(ni->in.size() != 0)  
348 - result[vi] = ni->in.front()->back(); 177 + /// Output the network as a string
  178 + std::string str(){
349 179
350 - vi++; //increment the array index 180 + std::stringstream ss;
  181 + ss<<"Nodes ("<<V.size()<<")--------"<<std::endl;
  182 + for(unsigned int v = 0; v < V.size(); v++){
  183 + ss<<"\t"<<v<<V[v].str()<<std::endl;
351 } 184 }
352 185
353 - //return the resulting array  
354 - return result;  
355 - } 186 + ss<<"Edges ("<<E.size()<<")--------"<<std::endl;
  187 + for(unsigned e = 0; e < E.size(); e++){
  188 + ss<<"\t"<<e<<E[e].str()<<std::endl;
  189 + }
356 190
357 - std::vector< stim::vec<T> > get_fiber_geometry( fiber_i f ){  
358 - return f->geometry(); 191 + return ss.str();
359 } 192 }
360 193
361 - /// Generate an OBJ file from the network  
362 -  
363 - stim::obj<T> obj(){  
364 -  
365 - //create an OBJ object  
366 - stim::obj<T> object;  
367 -  
368 - //name the nodes  
369 - set_names();  
370 -  
371 - //retrieve a list of terminal node positions  
372 - std::vector< stim::vec<T> > node_pos = get_node_positions();  
373 -  
374 - //add the nodes to the obj file  
375 - object.addV(node_pos);  
376 -  
377 - //counter for vertex indices in the object class  
378 - unsigned int nP;  
379 -  
380 - //for each fiber  
381 - fiber_i fi; //create a fiber iterator  
382 - for(fi = F.begin(); fi != F.end(); fi++){  
383 -  
384 - //get an array of fiber points  
385 - std::vector< stim::vec<T> > fiber_p = get_fiber_geometry(fi);  
386 -  
387 - //create a subset of this array  
388 - typename std::vector< stim::vec<T> >::iterator start = fiber_p.begin() + 1;  
389 - typename std::vector< stim::vec<T> >::iterator end = fiber_p.end() - 1;  
390 - typename std::vector< stim::vec<T> > fiber_subset(start, end);  
391 -  
392 - //add this subset to the geometry object  
393 - nP = object.addV(fiber_subset); 194 + /// This function resamples all fibers in a network given a desired minimum spacing
  195 + stim::network<T> resample(T spacing){
  196 + stim::network<T> n; //create a new network that will be an exact copy, with resampled fibers
  197 + n.V = V; //copy all vertices
394 198
395 - //create an array to hold vertex indices for a line  
396 - std::vector<unsigned int> line;  
397 - line.resize(fiber_p.size());  
398 -  
399 - //add the terminal nodes to the line list (make sure to add 1 to make them compatible with the OBJ)  
400 - line[0] = fi->n[0]->id + 1;  
401 - line[line.size() - 1] = fi->n[1]->id + 1;  
402 -  
403 - //add the intermediate vertex indices to the line array  
404 - for(unsigned int i = 0; i < fiber_subset.size(); i++){  
405 - line[1 + i] = nP + i;  
406 - }  
407 -  
408 - //add the line list to the object class  
409 - object.addLine(line); 199 + n.E.resize(edges()); //allocate space for the edge list
410 200
  201 + //copy all fibers, resampling them in the process
  202 + for(unsigned e = 0; e < edges(); e++){ //for each edge in the edge list
  203 + n.E[e] = E[e].resample(spacing); //resample the edge and copy it to the new network
411 } 204 }
412 205
413 - return object; 206 + return n; //return the resampled network
414 } 207 }
415 208
416 - /// This function returns the information necessary for a simple graph-based physical (ex. fluid) simulation.  
417 209
418 - /// @param n0 is a array which will contain the list of source nodes  
419 - /// @param n1 is a array which will contain the list of destination nodes  
420 - /// @param length is a array containing the lengths of fibers in the network  
421 - /// @param radius is a array containing the average radii of fibers in the network  
422 - void build_simgraph(std::vector<unsigned int>& n0, std::vector<unsigned int>& n1, std::vector<T>& length, std::vector<T>& radius){  
423 210
424 - //determine the number of fibers in the network  
425 - unsigned int nF = F.size();  
426 -  
427 - //allocate the necessary space to store the fiber information  
428 - n0.resize(nF);  
429 - n1.resize(nF);  
430 - length.resize(nF);  
431 - radius.resize(nF);  
432 -  
433 - //assign names (identifiers) to the network components  
434 - set_names(); 211 + // total number of points on all fibers after resampling -- function used only on a resampled network
  212 + unsigned total_points(){
  213 + unsigned n = 0;
  214 + for(unsigned e = 0; e < E.size(); e++)
  215 + n += E[e].size();
  216 + return n;
  217 + }
435 218
436 - //fill the arrays  
437 - unsigned int i = 0;  
438 - T l, r;  
439 - for(fiber_i f = F.begin(); f != F.end(); f++){  
440 - n0[i] = f->n[0]->id; //get the identifiers for the first and second nodes for the current fiber  
441 - n1[i] = f->n[1]->id; 219 + // gaussian function
  220 + float gaussianFunction(float x, float std=25){ return exp(-x/(2*std*std));} // by default std = 25
442 221
443 - r = f->radius(l); //get the length and radius of the capillary (calculated at the same time) 222 + // stim 3d vector to annpoint of 3 dimensions
  223 + void stim2ann(ANNpoint &a, stim::vec<T> b){
  224 + a[0] = b[0];
  225 + a[1] = b[1];
  226 + a[2] = b[2];
  227 + }
444 228
445 - radius[i] = r; //store the radius in the output array  
446 - length[i] = l; //store the length in the output array  
447 229
448 - i++; //increment the array index 230 + /// This function compares two networks and returns a metric
  231 + float compare(stim::network<T> A, float sigma){
  232 + float metric = 0.0; // initialize metric to be returned after comparing the networks
  233 + ANNkd_tree* kdt; // initialize a pointer to a kd tree
  234 + double **c; // centerline (array of double pointers) - points on kdtree must be double
  235 + unsigned int n_data = total_points(); // set the number of points
  236 + c = (double**) malloc(sizeof(double*) * n_data); // allocate the array pointer
  237 + for(unsigned int i = 0; i < n_data; i++) // allocate space for each point of 3 dimensions
  238 + c[i] = (double*) malloc(sizeof(double) * 3);
  239 + //std::vector<stim::vec<T> > pointsVector = points_afterResampling(); //get vertices in the network after resampling
  240 + unsigned t = 0;
  241 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  242 + for(unsigned p = 0; p < E[e].size(); p++){ //for each point in the edge
  243 + for(unsigned d = 0; d < 3; d++){ //for each coordinate
  244 +
  245 + c[t][d] = E[e][p][d];
  246 + }
  247 + t++;
  248 + }
449 } 249 }
450 250
  251 + ANNpointArray pts = (ANNpointArray)c; // create an array of data points of type double
  252 + kdt = new ANNkd_tree(pts, n_data, 3); // build a KD tree using the annpointarray
  253 + double eps = 0; // error bound
  254 + ANNdistArray dists = new ANNdist[1]; // near neighbor distances
  255 + ANNidxArray nnIdx = new ANNidx[1]; // near neighbor indices // allocate near neigh indices
  256 +
  257 + stim::vec<T> p0, p1;
  258 + float m0, m1;
  259 + float M = 0; //stores the total metric value
  260 + float l; //stores the segment length
  261 + float L = 0; //stores the total network length
  262 + ANNpoint queryPt = annAllocPt(3);
  263 + for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in A
  264 + M = 0; //total metric value for the fiber
  265 + p0 = A.E[e][0]; //get the first point in the edge
  266 + stim2ann(queryPt, p0);
  267 + kdt->annkSearch( queryPt, 1, nnIdx, dists, eps); //find the distance between A and the current network
  268 + m0 = 1.0f - gaussianFunction(dists[0], sigma); //calculate the metric value based on the distance
  269 +
  270 + for(unsigned p = 1; p < A.E[e].size(); p++){ //for each point in the edge
  271 +
  272 + p1 = A.E[e][p]; //get the next point in the edge
  273 + stim2ann(queryPt, p1);
  274 + kdt->annkSearch( queryPt, 1, nnIdx, dists, eps); //find the distance between A and the current network
  275 + m1 = 1.0f - gaussianFunction(dists[0], sigma); //calculate the metric value based on the distance
  276 +
  277 + //average the metric and scale
  278 + l = (p1 - p0).len(); //calculate the length of the segment
  279 + L += l; //add the length of this segment to the total network length
  280 + M += (m0 + m1)/2 * l; //scale the metric based on segment length
  281 +
  282 + //swap points
  283 + p0 = p1;
  284 + m0 = m1;
  285 + }
  286 + }
451 287
  288 + return M/L;
452 } 289 }
453 -  
454 -};  
455 -  
456 -}; //end namespace stim  
457 -  
458 - 290 +}; //end stim::network class
  291 +}; //end stim namespace
459 #endif 292 #endif
stim/biomodels/network_dep.h 0 → 100644
  1 +#ifndef STIM_NETWORK_H
  2 +#define STIM_NETWORK_H
  3 +
  4 +#include <stim/math/vector.h>
  5 +#include <stim/visualization/obj.h>
  6 +#include <list>
  7 +#include <ANN/ANN.h>
  8 +
  9 +namespace stim{
  10 +
  11 +/** This class provides an interface for dealing with biological networks.
  12 + * It takes the following aspects into account:
  13 + * 1) Network geometry and centerlines
  14 + * 2) Network connectivity (a graph structure can be extracted)
  15 + * 3) Network surface structure (the surface is represented as a triangular mesh and referenced to the centerline)
  16 + */
  17 +
  18 +template<typename T>
  19 +class network{
  20 +
  21 + //-------------------HELPER CLASSES-----------------------
  22 + /// Stores information about a geometric point on the network centerline (including point position and radius)
  23 + //template<typename T>
  24 + class point : public stim::vec<T>{
  25 +
  26 + public:
  27 + T r;
  28 +
  29 + point() : stim::vec<T>(){}
  30 +
  31 + //casting constructor
  32 + point(stim::vec<T> rhs) : stim::vec<T>(rhs){}
  33 + };
  34 +
  35 + //template<typename T>
  36 + class t_node;
  37 + class fiber;
  38 +
  39 + //create typedefs for the iterators to simplify the network code
  40 + typedef typename std::list< fiber >::iterator fiber_i;
  41 + typedef typename std::list< t_node >::iterator t_node_i;
  42 +
  43 + /// Stores information about a single capillary (a length of vessel between two branch or end points)
  44 + //template<typename T>
  45 + class fiber : public std::list< point >{
  46 +
  47 + using std::list< point >::begin;
  48 + using std::list< point >::end;
  49 + using std::list< point >::size;
  50 +
  51 +
  52 + public:
  53 + //std::list< point > P; //geometric point positions
  54 +
  55 + typename std::list< t_node >::iterator n[2]; //indices to terminal nodes
  56 + unsigned int id;
  57 +
  58 + public:
  59 +
  60 + /// Calculate the length of the fiber and return it.
  61 + T length(){
  62 +
  63 + point p0, p1;
  64 + T l = 0; //initialize the length to zero
  65 +
  66 + //for each point
  67 + typename std::list< point >::iterator i; //create a point iterator
  68 + for(i = begin(); i != end(); i++){ //for each point in the fiber
  69 +
  70 + if(i == begin()) //if this is the first point, just store it
  71 + p1 = *i;
  72 + else{ //if this is any other point
  73 + p0 = p1; //shift p1->p0
  74 + p1 = *i; //set p1 to the new point
  75 + l += (p1 - p0).len(); //add the length of p1 - p0 to the running sum
  76 + }
  77 + }
  78 +
  79 + return l; //return the length
  80 + }
  81 +
  82 + T radius(T& length){
  83 +
  84 + point p0, p1; //temporary variables to store point positions
  85 + T r0, r1; //temporary variables to store radii at points
  86 + T l, r; //temporary variable to store the length and average radius of a fiber segment
  87 + T length_sum = 0; //initialize the length to zero
  88 + T radius_sum = 0; //initialize the radius sum to zero
  89 +
  90 + //for each point
  91 + typename std::list< point >::iterator i; //create a point iterator
  92 + for(i = begin(); i != end(); i++){ //for each point in the fiber
  93 +
  94 + if(i == begin()){ //if this is the first point, just store it
  95 + p1 = *i;
  96 + r1 = i->r;
  97 + }
  98 + else{ //if this is any other point
  99 + p0 = p1; //shift p1->p0 and r1->r0
  100 + r0 = r1;
  101 + p1 = *i; //set p1 to the new point
  102 + r1 = i->r; //and r1
  103 +
  104 + l = (p1 - p0).len(); //calculate the length of the p0-p1 segment
  105 + r = (r0 + r1) / 2; //calculate the average radius of the segment
  106 +
  107 + radius_sum += r * l; //add the radius scaled by the length to a running sum
  108 + length_sum += l; //add the length of p1 - p0 to the running sum
  109 + }
  110 + }
  111 +
  112 + length = length_sum; //store the total length
  113 +
  114 + //if the total length is zero, store a radius of zero
  115 + if(length == 0)
  116 + return 0;
  117 + else
  118 + return radius_sum / length; //return the average radius of the fiber
  119 + }
  120 +
  121 + std::vector< stim::vec<T> > geometry(){
  122 +
  123 + std::vector< stim::vec<T> > result; //create an array to store the fiber geometry
  124 + result.resize( size() ); //pre-allocate the array
  125 +
  126 + typename std::list< point >::iterator p; //create a list iterator
  127 + unsigned int pi = 0; //create an index into the result array
  128 +
  129 + //for each geometric point on the fiber centerline
  130 + for(p = begin(); p != end(); p++){
  131 + result[pi] = *p;
  132 + pi++;
  133 + }
  134 +
  135 + return result; //return the geometry array
  136 +