codebase / stimlib

Commit 52b09bd9ba821d85032c832e3a7b3c020385964d

Authored by Pavel Govyadinov
2 parents 6ada8448 15a53122

Merge branch 'pranathi_stimlib' into 'master' 

Merge Pranathi_stimlib into master

See merge request !3
Showing 1 changed file with 248 additions and 4 deletions   Show diff stats
Inline Side-by-side
stim/biomodels/network.h
  Show/Hide comments   View file @52b09bd
... ... @@ -90,6 +90,7 @@ class network{
90 90 return ss.str();
91 91 }
92 92  
  93 +
93 94 };
94 95  
95 96 protected:
... ... @@ -109,6 +110,150 @@ public:
109 110 return V.size();
110 111 }
111 112  
  113 + std::vector<vertex> operator*(T s){
  114 + for (unsigned i=0; i< vertices; i ++ ){
  115 + V[i] = V[i] * s;
  116 + }
  117 + return V;
  118 + }
  119 +
  120 + std::vector<vertex> operator*(vec<T> s){
  121 + for (unsigned i=0; i< vertices; i ++ ){
  122 + for (unsigned dim = 0 ; dim< 3; dim ++){
  123 + V[i][dim] = V[i][dim] * s[dim];
  124 + }
  125 + }
  126 + return V;
  127 + }
  128 +
  129 + // Returns an average of branching index in the network
  130 +
  131 + double BranchingIndex(){
  132 + double B=0;
  133 + for(unsigned v=0; v < V.size(); v ++){
  134 + B += ((V[v].e[0].size()) + (V[v].e[1].size()));
  135 + }
  136 + B = B / V.size();
  137 + return B;
  138 +
  139 + }
  140 +
  141 + // Returns number of branch points in thenetwork
  142 +
  143 + unsigned int BranchP(){
  144 + unsigned int B=0;
  145 + unsigned int c;
  146 + for(unsigned v=0; v < V.size(); v ++){
  147 + c = ((V[v].e[0].size()) + (V[v].e[1].size()));
  148 + if (c > 2){
  149 + B += 1;}
  150 + }
  151 + return B;
  152 +
  153 + }
  154 +
  155 + // Returns number of end points (tips) in thenetwork
  156 +
  157 + unsigned int EndP(){
  158 + unsigned int B=0;
  159 + unsigned int c;
  160 + for(unsigned v=0; v < V.size(); v ++){
  161 + c = ((V[v].e[0].size()) + (V[v].e[1].size()));
  162 + if (c == 1){
  163 + B += 1;}
  164 + }
  165 + return B;
  166 +
  167 + }
  168 +
  169 + //// Returns a dictionary with the key as the vertex
  170 + //std::map<std::vector<vertex>,unsigned int> DegreeDict(){
  171 + // std::map<std::vector<vertex>,unsigned int> dd;
  172 + // unsigned int c = 0;
  173 + // for(unsigned v=0; v < V.size(); v ++){
  174 + // c = ((V[v].e[0].size()) + (V[v].e[1].size()));
  175 + // dd[V[v]] = c;
  176 + // }
  177 + // return dd;
  178 + //}
  179 +
  180 + //// Return number of branching stems
  181 + //unsigned int Stems(){
  182 + // unsigned int s = 0;
  183 + // std::map<std::vector<vertex>,unsigned int> dd;
  184 + // dd = DegreeDict();
  185 + // //for(unsigned v=0; v < V.size(); v ++){
  186 + // // V[v].e[0].
  187 + // return s;
  188 + //}
  189 +
  190 +
  191 + // Returns an average of fiber/edge lengths in the network
  192 + double Lengths(){
  193 + stim::vec<T> L;double sumLength = 0;
  194 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  195 + L.push_back(E[e].length()); //append the edge length
  196 + sumLength = sumLength + E[e].length();
  197 + }
  198 + double avg = sumLength / E.size();
  199 + return avg;
  200 + }
  201 +
  202 +
  203 + // Returns an average of tortuosities in the network
  204 + double Tortuosities(){
  205 + stim::vec<T> t;
  206 + stim::vec<T> id1, id2; // starting and ending vertices of the edge
  207 + double distance;double tortuosity;double sumTortuosity = 0;
  208 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  209 + id1 = E[e][0]; //get the edge starting point
  210 + id2 = E[e][E[e].size() - 1]; //get the edge ending point
  211 + distance = (id1 - id2).len(); //displacement between the starting and ending points
  212 + if(distance > 0){
  213 + tortuosity = E[e].length()/ distance ; // tortuoisty = edge length / edge displacement
  214 + }
  215 + else{
  216 + tortuosity = 0;}
  217 + t.push_back(tortuosity);
  218 + sumTortuosity += tortuosity;
  219 + }
  220 + double avg = sumTortuosity / E.size();
  221 + return avg;
  222 + }
  223 +
  224 + // Returns average contraction of the network
  225 + double Contractions(){
  226 + stim::vec<T> t;
  227 + stim::vec<T> id1, id2; // starting and ending vertices of the edge
  228 + double distance;double contraction;double sumContraction = 0;
  229 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  230 + id1 = E[e][0]; //get the edge starting point
  231 + id2 = E[e][E[e].size() - 1]; //get the edge ending point
  232 + distance = (id1 - id2).len(); //displacement between the starting and ending points
  233 + contraction = distance / E[e].length(); // tortuoisty = edge length / edge displacement
  234 + t.push_back(contraction);
  235 + sumContraction += contraction;
  236 + }
  237 + double avg = sumContraction / E.size();
  238 + return avg;
  239 + }
  240 +
  241 + // returns average fractal dimension of the branches of the network
  242 + double FractalDimensions(){
  243 + stim::vec<T> t;
  244 + stim::vec<T> id1, id2; // starting and ending vertices of the edge
  245 + double distance;double fract;double sumFractDim = 0;
  246 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  247 + id1 = E[e][0]; //get the edge starting point
  248 + id2 = E[e][E[e].size() - 1]; //get the edge ending point
  249 + distance = (id1 - id2).len(); //displacement between the starting and ending points
  250 + fract = std::log(distance) / std::log(E[e].length()); // tortuoisty = edge length / edge displacement
  251 + t.push_back(sumFractDim);
  252 + sumFractDim += fract;
  253 + }
  254 + double avg = sumFractDim / E.size();
  255 + return avg;
  256 + }
112 257 stim::cylinder<T> get_cylinder(unsigned f){
113 258 return E[f]; //return the specified edge (casting it to a fiber)
114 259 }
... ... @@ -147,7 +292,6 @@ public:
147 292  
148 293 //find out if the nodes for this fiber have already been created
149 294 it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node
150   - it_idx = std::distance(id2vert.begin(), it);
151 295 if(it == id2vert.end()){ //if i[0] hasn't already been used
152 296 vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position
153 297 new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing
... ... @@ -156,20 +300,21 @@ public:
156 300 id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list
157 301 }
158 302 else{ //if the vertex already exists
  303 + it_idx = std::distance(id2vert.begin(), it);
159 304 V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing
160 305 new_edge.v[0] = it_idx;
161 306 }
162 307  
163 308 it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID
164   - it_idx = std::distance(id2vert.begin(), it);
165 309 if(it == id2vert.end()){ //if i[1] hasn't already been used
166 310 vertex new_vertex = new_edge[I-1]; //create a new vertex, assign it a position
167 311 new_vertex.e[1].push_back(E.size()); //add the current edge as incoming
168   - new_edge.v[1] = V.size();
  312 + new_edge.v[1] = V.size(); //add the new vertex to the edge
169 313 V.push_back(new_vertex); //add the new vertex to the vertex list
170 314 id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list
171 315 }
172 316 else{ //if the vertex already exists
  317 + it_idx = std::distance(id2vert.begin(), it);
173 318 V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming
174 319 new_edge.v[1] = it_idx;
175 320 }
... ... @@ -308,8 +453,107 @@ public:
308 453 unsigned nmags(){
309 454 return E[0].nmags();
310 455 }
  456 + // split a string in text by the character sep
  457 + stim::vec<T> split(std::string &text, char sep)
  458 + {
  459 + stim::vec<T> tokens;
  460 + std::size_t start = 0, end = 0;
  461 + while ((end = text.find(sep, start)) != std::string::npos) {
  462 + tokens.push_back(atof(text.substr(start, end - start).c_str()));
  463 + start = end + 1;
  464 + }
  465 + tokens.push_back(atof(text.substr(start).c_str()));
  466 + return tokens;
  467 + }
  468 + // load a network in text file to a network class
  469 + void load_txt(std::string filename)
  470 + {
  471 + std::vector <std::string> file_contents;
  472 + std::ifstream file(filename);
  473 + std::string line;
  474 + std::vector<unsigned> id2vert; //this list stores the vertex ID associated with each network vertex
  475 + //for each line in the text file, store them as strings in file_contents
  476 + while (std::getline(file, line))
  477 + {
  478 + std::stringstream ss(line);
  479 + file_contents.push_back(ss.str());
  480 + }
  481 + unsigned int numEdges = atoi(file_contents[0].c_str()); //number of edges in the network
  482 + unsigned int I = atoi(file_contents[1].c_str()) ; //calculate the number of points3d on the first edge
  483 + unsigned int count = 1; unsigned int k = 2; // count is global counter through the file contents, k is for the vertices on the edges
  484 + // for each edge in the network.
  485 + for (unsigned int i = 0; i < numEdges; i ++ )
  486 + {
  487 + // pre allocate a position vector p with number of points3d on the edge p
  488 + std::vector< stim::vec<T> > p(0, I);
  489 + // for each point on the nth edge
  490 + for (unsigned int j = k; j < I + k; j++)
  491 + {
  492 + // split the points3d of floats with separator space and form a float3 position vector out of them
  493 + p.push_back(split(file_contents[j], ' '));
  494 + }
  495 + count += p.size() + 1; // increment count to point at the next edge in the network
  496 + I = atoi(file_contents[count].c_str()); // read in the points3d at the next edge and convert it to an integer
  497 + k = count + 1;
  498 + edge new_edge = p; // create an edge with a vector of points3d on the edge
  499 + E.push_back(new_edge); // push the edge into the network
  500 + }
  501 + unsigned int numVertices = atoi(file_contents[count].c_str()); // this line in the text file gives the number of distinct vertices
  502 + count = count + 1; // this line of text file gives the first verrtex
  503 + // push each vertex into V
  504 + for (unsigned int i = 0; i < numVertices; i ++)
  505 + {
  506 + vertex new_vertex = split(file_contents[count], ' ');
  507 + V.push_back(new_vertex);
  508 + count += atoi(file_contents[count + 1].c_str()) + 2; // Skip number of edge ids + 2 to point to the next vertex
  509 + }
  510 + } // end load_txt function
311 511  
312   -
  512 + // strTxt returns a string of edges
  513 + std::string
  514 + strTxt(std::vector< stim::vec<T> > p)
  515 + {
  516 + std::stringstream ss;
  517 + std::stringstream oss;
  518 + for(unsigned int i = 0; i < p.size(); i++){
  519 + ss.str(std::string());
  520 + for(unsigned int d = 0; d < 3; d++){
  521 + ss<<c[i][d];
  522 + }
  523 + ss < "\n"
  524 + }
  525 + return ss.str();
  526 + }
  527 + // removes specified character from string
  528 + void removeCharsFromString(std::string &str, char* charsToRemove ) {
  529 + for ( unsigned int i = 0; i < strlen(charsToRemove); ++i ) {
  530 + str.erase( remove(str.begin(), str.end(), charsToRemove[i]), str.end() );
  531 + }
  532 + }
  533 + //exports network to txt file
  534 + void
  535 + to_txt(std::string filename)
  536 + {
  537 + std::ofstream ofs;
  538 + ofs.open(filename, std::ofstream::out | std::ofstream::app);
  539 + int num;
  540 + ofs << (E.size()).str() << "\n";
  541 + for(unsigned int i = 0; i < E.size(); i++)
  542 + {
  543 + std::string str;
  544 + ofs << (E[i].size()).str() << "\n";
  545 + str = E[i].strTxt();
  546 + ofs << str << "\n";
  547 + }
  548 + for(int i = 0; i < V.size(); i++)
  549 + {
  550 + std::string str;
  551 + str = V[i].str();
  552 + removeCharsFromString(str, "[],");
  553 + ofs << str << "\n";
  554 + }
  555 + ofs.close();
  556 + }
313 557 }; //end stim::network class
314 558 }; //end stim namespace
315 559 #endif
... ...