Commit ee63d5f7570ef3fd3925f4b1dd0fff2bfb9679e3

Authored by Jiaming Guo
1 parent 86061a86

modified network for flow3

Showing 1 changed file with 141 additions and 1 deletions   Show diff stats
stim/biomodels/network.h
... ... @@ -473,6 +473,146 @@ public:
473 473 return E[e]; //return the specified edge (casting it to a fiber)
474 474 }
475 475  
  476 + /// subdivide current network
  477 + void subdivision() {
  478 +
  479 + std::vector<unsigned> ori_index; // original index
  480 + std::vector<unsigned> new_index; // new index
  481 + std::vector<edge> nE; // new edge
  482 + std::vector<vertex> nV; // new vector
  483 + unsigned id = 0;
  484 +
  485 + for (unsigned i = 0; i < num_edge; i++) {
  486 + if (E[i].size() == 2) { // if current edge can't be subdivided
  487 + stim::centerline<T> line(2);
  488 + for (unsigned k = 0; k < 2; k++)
  489 + line[k] = E[i][k];
  490 + line.update();
  491 +
  492 + edge new_edge(line);
  493 +
  494 + vertex new_vertex = new_edge[0];
  495 + id = E[i].v[0];
  496 + auto position = std::find(ori_index.begin(), ori_index.end(), id);
  497 + if (position == ori_index.end()) { // new vertex
  498 + ori_index.push_back(id);
  499 + new_index.push_back(nV.size());
  500 +
  501 + new_vertex.e[0].push_back(nE.size());
  502 + new_edge.v[0] = nV.size();
  503 + nV.push_back(new_vertex); // push back vertex as a new vertex
  504 + }
  505 + else { // existing vertex
  506 + int k = std::distance(ori_index.begin(), position);
  507 + new_edge.v[0] = new_index[k];
  508 + nV[new_index[k]].e[0].push_back(nE.size());
  509 + }
  510 +
  511 + new_vertex = new_edge[1];
  512 + id = E[i].v[1];
  513 + position = std::find(ori_index.begin(), ori_index.end(), id);
  514 + if (position == ori_index.end()) { // new vertex
  515 + ori_index.push_back(id);
  516 + new_index.push_back(nV.size());
  517 +
  518 + new_vertex.e[1].push_back(nE.size());
  519 + new_edge.v[1] = nV.size();
  520 + nV.push_back(new_vertex); // push back vertex as a new vertex
  521 + }
  522 + else { // existing vertex
  523 + int k = std::distance(ori_index.begin(), position);
  524 + new_edge.v[1] = new_index[k];
  525 + nV[new_index[k]].e[1].push_back(nE.size());
  526 + }
  527 +
  528 + nE.push_back(new_edge);
  529 +
  530 + nE[nE.size() - 1].cylinder<T>::set_r(0, E[i].cylinder<T>::r(0));
  531 + nE[nE.size() - 1].cylinder<T>::set_r(1, E[i].cylinder<T>::r(1));
  532 + }
  533 + else { // subdivide current edge
  534 + for (unsigned j = 0; j < E[i].size() - 1; j++) {
  535 + stim::centerline<T> line(2);
  536 + for (unsigned k = 0; k < 2; k++)
  537 + line[k] = E[i][j + k];
  538 + line.update();
  539 +
  540 + edge new_edge(line);
  541 +
  542 + if (j == 0) { // edge contains original starting point
  543 + vertex new_vertex = new_edge[0];
  544 + id = E[i].v[0];
  545 + auto position = std::find(ori_index.begin(), ori_index.end(), id);
  546 + if (position == ori_index.end()) { // new vertex
  547 + ori_index.push_back(id);
  548 + new_index.push_back(nV.size());
  549 +
  550 + new_vertex.e[0].push_back(nE.size());
  551 + new_edge.v[0] = nV.size();
  552 + nV.push_back(new_vertex); // push back vertex as a new vertex
  553 + }
  554 + else { // existing vertex
  555 + int k = std::distance(ori_index.begin(), position);
  556 + new_edge.v[0] = new_index[k];
  557 + nV[new_index[k]].e[0].push_back(nE.size());
  558 + }
  559 +
  560 + new_vertex = new_edge[1];
  561 + new_vertex.e[1].push_back(nE.size());
  562 + new_edge.v[1] = nV.size();
  563 + nV.push_back(new_vertex); // push back internal point as a new vertex
  564 +
  565 + nE.push_back(new_edge);
  566 + }
  567 +
  568 + else if (j == E[i].size() - 2) { // edge contains original ending point
  569 +
  570 + vertex new_vertex = new_edge[1];
  571 + nV[nV.size() - 1].e[0].push_back(nE.size());
  572 + new_edge.v[0] = nV.size() - 1;
  573 +
  574 + id = E[i].v[1];
  575 + auto position = std::find(ori_index.begin(), ori_index.end(), id);
  576 + if (position == ori_index.end()) { // new vertex
  577 + ori_index.push_back(id);
  578 + new_index.push_back(nV.size());
  579 +
  580 + new_vertex.e[1].push_back(nE.size());
  581 + new_edge.v[1] = nV.size();
  582 + nV.push_back(new_vertex); // push back vertex as a new vertex
  583 + }
  584 + else { // existing vertex
  585 + int k = std::distance(ori_index.begin(), position);
  586 + new_edge.v[1] = new_index[k];
  587 + nV[new_index[k]].e[1].push_back(nE.size());
  588 + }
  589 +
  590 + nE.push_back(new_edge);
  591 + }
  592 +
  593 + else {
  594 + vertex new_vertex = new_edge[1];
  595 +
  596 + nV[nV.size() - 1].e[0].push_back(nE.size());
  597 + new_vertex.e[1].push_back(nE.size());
  598 + new_edge.v[0] = nV.size() - 1;
  599 + new_edge.v[1] = nV.size();
  600 + nV.push_back(new_vertex);
  601 +
  602 + nE.push_back(new_edge);
  603 + }
  604 +
  605 + // get radii
  606 + nE[nE.size() - 1].cylinder<T>::set_r(0, E[i].cylinder<T>::r(j));
  607 + nE[nE.size() - 1].cylinder<T>::set_r(1, E[i].cylinder<T>::r(j + 1));
  608 + }
  609 + }
  610 + }
  611 +
  612 + (*this).E = nE;
  613 + (*this).V = nV;
  614 + }
  615 +
476 616 //load a network from an OBJ file
477 617 void load_obj(std::string filename){
478 618  
... ... @@ -715,7 +855,7 @@ public:
715 855 edge new_edge(C3); // new edge
716 856  
717 857 //create an edge from the given centerline
718   - unsigned int I = new_edge.size(); //calculate the number of points on the centerline
  858 + unsigned int I = (unsigned)new_edge.size(); //calculate the number of points on the centerline
719 859  
720 860 //get the first and last vertex IDs for the line
721 861 i[0] = S.E[l].front();
... ...