Commit b2f48bca486336b45a348492262b2e3b6bf821aa

Authored by Pavel Govyadinov
2 parents 2e0f052a 52b09bd9

fixed merge conflicts

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