Commit f1bb77987f4b4912b5327d148963de4487a4f52b

Authored by David Mayerich
1 parent 9832489d

moved the static network class into the netmets repository from STIMLIB, where it is no longer used

centerline.h 0 → 100644
  1 +#ifndef STIM_CENTERLINE_H
  2 +#define STIM_CENTERLINE_H
  3 +
  4 +#include <vector>
  5 +#include <stim/math/vec3.h>
  6 +#include <stim/structures/kdtree.cuh>
  7 +
  8 +namespace stim{
  9 +
  10 +/** This class stores information about a single fiber represented as a set of geometric points
  11 + * between two branch or end points. This class is used as a fundamental component of the stim::network
  12 + * class to describe an interconnected (often biological) network.
  13 + */
  14 +template<typename T>
  15 +class centerline : public std::vector< stim::vec3<T> >{
  16 +
  17 +protected:
  18 +
  19 + std::vector<T> L; //stores the integrated length along the fiber (used for parameterization)
  20 +
  21 + ///Return the normalized direction vector at point i (average of the incoming and outgoing directions)
  22 + vec3<T> d(size_t i) {
  23 + if (size() <= 1) return vec3<T>(0, 0, 0); //if there is insufficient information to calculate the direction, return a null vector
  24 + if (size() == 2) return (at(1) - at(0)).norm(); //if there are only two points, the direction vector at both is the direction of the line segment
  25 + if (i == 0) return (at(1) - at(0)).norm(); //the first direction vector is oriented towards the first line segment
  26 + if (i == size() - 1) return (at(size() - 1) - at(size() - 2)).norm(); //the last direction vector is oriented towards the last line segment
  27 +
  28 + //all other direction vectors are the average direction of the two joined line segments
  29 + vec3<T> a = at(i) - at(i - 1);
  30 + vec3<T> b = at(i + 1) - at(i);
  31 + vec3<T> ab = a.norm() + b.norm();
  32 + return ab.norm();
  33 + }
  34 +
  35 + //initializes the integrated length vector to make parameterization easier, starting with index idx (all previous indices are assumed to be correct)
  36 + void update_L(size_t start = 0) {
  37 + L.resize(size()); //allocate space for the L array
  38 + if (start == 0) {
  39 + L[0] = 0; //initialize the length value for the first point to zero (0)
  40 + start++;
  41 + }
  42 +
  43 + stim::vec3<T> d;
  44 + for (size_t i = start; i < size(); i++) { //for each line segment in the centerline
  45 + d = at(i) - at(i - 1);
  46 + L[i] = L[i - 1] + d.len(); //calculate the running length total
  47 + }
  48 + }
  49 +
  50 + void init() {
  51 + if (size() == 0) return; //return if there aren't any points
  52 + update_L();
  53 + }
  54 +
  55 + /// Returns a stim::vec representing the point at index i
  56 +
  57 + /// @param i is an index of the desired centerline point
  58 + stim::vec<T> get_vec(unsigned i){
  59 + return std::vector< stim::vec3<T> >::at(i);
  60 + }
  61 +
  62 + ///finds the index of the point closest to the length l on the lower bound.
  63 + ///binary search.
  64 + size_t findIdx(T l) {
  65 + for (size_t i = 1; i < L.size(); i++) { //for each point in the centerline
  66 + if (L[i] > l) return i - 1; //if we have passed the desired length value, return i-1
  67 + }
  68 + return L.size() - 1;
  69 + /*size_t i = L.size() / 2;
  70 + size_t max = L.size() - 1;
  71 + size_t min = 0;
  72 + while (i < L.size() - 1){
  73 + if (l < L[i]) {
  74 + max = i;
  75 + i = min + (max - min) / 2;
  76 + }
  77 + else if (L[i] <= l && L[i + 1] >= l) {
  78 + break;
  79 + }
  80 + else {
  81 + min = i;
  82 + i = min + (max - min) / 2;
  83 + }
  84 + }
  85 + return i;*/
  86 + }
  87 +
  88 + ///Returns a position vector at the given length into the fiber (based on the pvalue).
  89 + ///Interpolates the radius along the line.
  90 + ///@param l: the location of the in the cylinder.
  91 + ///@param idx: integer location of the point closest to l but prior to it.
  92 + stim::vec3<T> p(T l, int idx) {
  93 + T rat = (l - L[idx]) / (L[idx + 1] - L[idx]);
  94 + stim::vec3<T> v1 = at(idx);
  95 + stim::vec3<T> v2 = at(idx + 1);
  96 + return(v1 + (v2 - v1)*rat);
  97 + }
  98 +
  99 +
  100 +public:
  101 +
  102 + using std::vector< stim::vec3<T> >::at;
  103 + using std::vector< stim::vec3<T> >::size;
  104 +
  105 + centerline() : std::vector< stim::vec3<T> >() {
  106 + init();
  107 + }
  108 + centerline(size_t n) : std::vector< stim::vec3<T> >(n){
  109 + init();
  110 + }
  111 + centerline(std::vector<stim::vec3<T> > pos) :
  112 + std::vector<stim::vec3<T> > (pos)
  113 + {
  114 + init();
  115 + }
  116 +
  117 + //overload the push_back function to update the length vector
  118 + void push_back(stim::vec3<T> p) {
  119 + std::vector< stim::vec3<T> >::push_back(p);
  120 + update_L(size() - 1);
  121 + }
  122 +
  123 + ///Returns a position vector at the given p-value (p value ranges from 0 to 1).
  124 + ///interpolates the position along the line.
  125 + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
  126 + stim::vec3<T> p(T pvalue) {
  127 + if (pvalue <= 0.0) return at(0); //return the first element
  128 + if (pvalue >= 1.0) return back(); //return the last element
  129 +
  130 + T l = pvalue*L[L.size() - 1];
  131 + int idx = findIdx(l);
  132 + return p(l, idx);
  133 + }
  134 +
  135 + ///Update centerline internal parameters (currently the L vector)
  136 + void update() {
  137 + init();
  138 + }
  139 + ///Return the length of the entire centerline
  140 + T length() {
  141 + return L.back();
  142 + }
  143 +
  144 +
  145 + /// stitch two centerlines
  146 + ///@param c1, c2: two centerlines
  147 + ///@param sigma: sample rate
  148 + static std::vector< stim::centerline<T> > stitch(stim::centerline<T> c1, stim::centerline<T> c2 = stim::centerline<T>()) {
  149 +
  150 + std::vector< stim::centerline<T> > result;
  151 + stim::centerline<T> new_centerline;
  152 + stim::vec3<T> new_vertex;
  153 +
  154 + // if only one centerline, stitch itself!
  155 + if (c2.size() == 0) {
  156 + size_t num = c1.size();
  157 + size_t id = 100000; // store the downsample start position
  158 + T threshold;
  159 + if (num < 4) { // if the number of vertex is less than 4, do nothing
  160 + result.push_back(c1);
  161 + return result;
  162 + }
  163 + else {
  164 + // test geometry start vertex
  165 + stim::vec3<T> v1 = c1[1] - c1[0]; // vector from c1[0] to c1[1]
  166 + for (size_t p = 2; p < num; p++) { // 90° standard???
  167 + stim::vec3<T> v2 = c1[p] - c1[0];
  168 + float cosine = v2.dot(v1);
  169 + if (cosine < 0) {
  170 + id = p;
  171 + threshold = v2.len();
  172 + break;
  173 + }
  174 + }
  175 + if (id != 100000) { // find a downsample position on the centerline
  176 + T* c;
  177 + c = (T*)malloc(sizeof(T) * (num - id) * 3);
  178 + for (size_t p = id; p < num; p++) {
  179 + for (size_t d = 0; d < 3; d++) {
  180 + c[p * 3 + d] = c1[p][d];
  181 + }
  182 + }
  183 + stim::kdtree<T, 3> kdt;
  184 + kdt.create(c, num - id, 5); // create tree
  185 +
  186 + T* query = (T*)malloc(sizeof(T) * 3);
  187 + for (size_t d = 0; d < 3; d++)
  188 + query[d] = c1[0][d];
  189 + size_t index;
  190 + T dist;
  191 +
  192 + kdt.search(query, 1, &index, &dist);
  193 +
  194 + free(query);
  195 + free(c);
  196 +
  197 + if (dist > threshold) {
  198 + result.push_back(c1);
  199 + }
  200 + else {
  201 + // the loop part
  202 + new_vertex = c1[index];
  203 + new_centerline.push_back(new_vertex);
  204 + for (size_t p = 0; p < index + 1; p++) {
  205 + new_vertex = c1[p];
  206 + new_centerline.push_back(new_vertex);
  207 + }
  208 + result.push_back(new_centerline);
  209 + new_centerline.clear();
  210 +
  211 + // the tail part
  212 + for (size_t p = index; p < num; p++) {
  213 + new_vertex = c1[p];
  214 + new_centerline.push_back(new_vertex);
  215 + }
  216 + result.push_back(new_centerline);
  217 + }
  218 + }
  219 + else { // there is one potential problem that two positions have to be stitched
  220 + // test geometry end vertex
  221 + stim::vec3<T> v1 = c1[num - 2] - c1[num - 1];
  222 + for (size_t p = num - 2; p > 0; p--) { // 90° standard
  223 + stim::vec3<T> v2 = c1[p - 1] - c1[num - 1];
  224 + float cosine = v2.dot(v1);
  225 + if (cosine < 0) {
  226 + id = p;
  227 + threshold = v2.len();
  228 + break;
  229 + }
  230 + }
  231 + if (id != 100000) { // find a downsample position
  232 + T* c;
  233 + c = (T*)malloc(sizeof(T) * (id + 1) * 3);
  234 + for (size_t p = 0; p < id + 1; p++) {
  235 + for (size_t d = 0; d < 3; d++) {
  236 + c[p * 3 + d] = c1[p][d];
  237 + }
  238 + }
  239 + stim::kdtree<T, 3> kdt;
  240 + kdt.create(c, id + 1, 5); // create tree
  241 +
  242 + T* query = (T*)malloc(sizeof(T) * 1 * 3);
  243 + for (size_t d = 0; d < 3; d++)
  244 + query[d] = c1[num - 1][d];
  245 + size_t index;
  246 + T dist;
  247 +
  248 + kdt.search(query, 1, &index, &dist);
  249 +
  250 + free(query);
  251 + free(c);
  252 +
  253 + if (dist > threshold) {
  254 + result.push_back(c1);
  255 + }
  256 + else {
  257 + // the tail part
  258 + for (size_t p = 0; p < index + 1; p++) {
  259 + new_vertex = c1[p];
  260 + new_centerline.push_back(new_vertex);
  261 + }
  262 + result.push_back(new_centerline);
  263 + new_centerline.clear();
  264 +
  265 + // the loop part
  266 + for (size_t p = index; p < num; p++) {
  267 + new_vertex = c1[p];
  268 + new_centerline.push_back(new_vertex);
  269 + }
  270 + new_vertex = c1[index];
  271 + new_centerline.push_back(new_vertex);
  272 + result.push_back(new_centerline);
  273 + }
  274 + }
  275 + else { // no stitch position
  276 + result.push_back(c1);
  277 + }
  278 + }
  279 + }
  280 + }
  281 +
  282 +
  283 + // two centerlines
  284 + else {
  285 + // find stitch position based on nearest neighbors
  286 + size_t num1 = c1.size();
  287 + T* c = (T*)malloc(sizeof(T) * num1 * 3); // c1 as reference point
  288 + for (size_t p = 0; p < num1; p++) // centerline to array
  289 + for (size_t d = 0; d < 3; d++) // because right now my kdtree code is a relatively close code, it has its own structure
  290 + c[p * 3 + d] = c1[p][d]; // I will merge it into stimlib totally in the near future
  291 +
  292 + stim::kdtree<T, 3> kdt; // kdtree object
  293 + kdt.create(c, num1, 5); // create tree
  294 +
  295 + size_t num2 = c2.size();
  296 + T* query = (T*)malloc(sizeof(T) * num2 * 3); // c2 as query point
  297 + for (size_t p = 0; p < num2; p++) {
  298 + for (size_t d = 0; d < 3; d++) {
  299 + query[p * 3 + d] = c2[p][d];
  300 + }
  301 + }
  302 + std::vector<size_t> index(num2);
  303 + std::vector<T> dist(num2);
  304 +
  305 + kdt.search(query, num2, &index[0], &dist[0]); // find the nearest neighbors in c1 for c2
  306 +
  307 + // clear up
  308 + free(query);
  309 + free(c);
  310 +
  311 + // find the average vertex distance of one centerline
  312 + T sigma1 = 0;
  313 + T sigma2 = 0;
  314 + for (size_t p = 0; p < num1 - 1; p++)
  315 + sigma1 += (c1[p] - c1[p + 1]).len();
  316 + for (size_t p = 0; p < num2 - 1; p++)
  317 + sigma2 += (c2[p] - c2[p + 1]).len();
  318 + sigma1 /= (num1 - 1);
  319 + sigma2 /= (num2 - 1);
  320 + float threshold = 4 * (sigma1 + sigma2) / 2; // better way to do this?
  321 +
  322 + T min_d = *std::min_element(dist.begin(), dist.end()); // find the minimum distance between c1 and c2
  323 +
  324 + if (min_d > threshold) { // if the minimum distance is too large
  325 + result.push_back(c1);
  326 + result.push_back(c2);
  327 +
  328 +#ifdef DEBUG
  329 + std::cout << "The distance between these two centerlines is too large" << std::endl;
  330 +#endif
  331 + }
  332 + else {
  333 + // auto smallest = std::min_element(dist.begin(), dist.end());
  334 + unsigned int smallest = std::min_element(dist.begin(), dist.end());
  335 + // auto i = std::distance(dist.begin(), smallest); // find the index of min-distance in distance list
  336 + unsigned int i = std::distance(dist.begin(), smallest); // find the index of min-distance in distance list
  337 +
  338 + // stitch position in c1 and c2
  339 + int id1 = index[i];
  340 + int id2 = i;
  341 +
  342 + // actually there are two cases
  343 + // first one inacceptable
  344 + // second one acceptable
  345 + if (id1 != 0 && id1 != num1 - 1 && id2 != 0 && id2 != num2 - 1) { // only stitch one end vertex to another centerline
  346 + result.push_back(c1);
  347 + result.push_back(c2);
  348 + }
  349 + else {
  350 + if (id1 == 0 || id1 == num1 - 1) { // if the stitch vertex is the first or last vertex of c1
  351 + // for c2, consider two cases(one degenerate case)
  352 + if (id2 == 0 || id2 == num2 - 1) { // case 1, if stitch position is also on the end of c2
  353 + // we have to decide which centerline get a new vertex, based on direction
  354 + // for c1, computer the direction change angle
  355 + stim::vec3<T> v1, v2;
  356 + float alpha1, alpha2; // direction change angle
  357 + if (id1 == 0)
  358 + v1 = (c1[1] - c1[0]).norm();
  359 + else
  360 + v1 = (c1[num1 - 2] - c1[num1 - 1]).norm();
  361 + v2 = (c2[id2] - c1[id1]).norm();
  362 + alpha1 = v1.dot(v2);
  363 + if (id2 == 0)
  364 + v1 = (c2[1] - c2[0]).norm();
  365 + else
  366 + v1 = (c2[num2 - 2] - c2[num2 - 1]).norm();
  367 + v2 = (c1[id1] - c2[id2]).norm();
  368 + alpha2 = v1.dot(v2);
  369 + if (abs(alpha1) > abs(alpha2)) { // add the vertex to c1 in order to get smooth connection
  370 + // push back c1
  371 + if (id1 == 0) { // keep geometry information
  372 + new_vertex = c2[id2];
  373 + new_centerline.push_back(new_vertex);
  374 + for (size_t p = 0; p < num1; p++) { // stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
  375 + new_vertex = c1[p];
  376 + new_centerline.push_back(new_vertex);
  377 + }
  378 + }
  379 + else {
  380 + for (size_t p = 0; p < num1; p++) { // stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1
  381 + new_vertex = c1[p];
  382 + new_centerline.push_back(new_vertex);
  383 + }
  384 + new_vertex = c2[id2];
  385 + new_centerline.push_back(new_vertex);
  386 + }
  387 + result.push_back(new_centerline);
  388 + new_centerline.clear();
  389 +
  390 + // push back c2
  391 + for (size_t p = 0; p < num2; p++) {
  392 + new_vertex = c2[p];
  393 + new_centerline.push_back(new_vertex);
  394 + }
  395 + result.push_back(new_centerline);
  396 + }
  397 + else { // add the vertex to c2 in order to get smooth connection
  398 + // push back c1
  399 + for (size_t p = 0; p < num1; p++) {
  400 + new_vertex = c1[p];
  401 + new_centerline.push_back(new_vertex);
  402 + }
  403 + result.push_back(new_centerline);
  404 + new_centerline.clear();
  405 +
  406 + // push back c2
  407 + if (id2 == 0) { // keep geometry information
  408 + new_vertex = c1[id1];
  409 + new_centerline.push_back(new_vertex);
  410 + for (size_t p = 0; p < num2; p++) { // stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
  411 + new_vertex = c2[p];
  412 + new_centerline.push_back(new_vertex);
  413 + }
  414 + }
  415 + else {
  416 + for (size_t p = 0; p < num2; p++) { // stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1
  417 + new_vertex = c2[p];
  418 + new_centerline.push_back(new_vertex);
  419 + }
  420 + new_vertex = c1[id1];
  421 + new_centerline.push_back(new_vertex);
  422 + }
  423 + result.push_back(new_centerline);
  424 + }
  425 + }
  426 + else { // case 2, the stitch position is on c2
  427 + // push back c1
  428 + if (id1 == 0) { // keep geometry information
  429 + new_vertex = c2[id2];
  430 + new_centerline.push_back(new_vertex);
  431 + for (size_t p = 0; p < num1; p++) { // stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
  432 + new_vertex = c1[p];
  433 + new_centerline.push_back(new_vertex);
  434 + }
  435 + }
  436 + else {
  437 + for (size_t p = 0; p < num1; p++) { // geometry end vertex on c1 -> geometry start vertex on c1 -> stitch vertex on c2
  438 + new_vertex = c1[p];
  439 + new_centerline.push_back(new_vertex);
  440 + }
  441 + new_vertex = c2[id2];
  442 + new_centerline.push_back(new_vertex);
  443 + }
  444 + result.push_back(new_centerline);
  445 + new_centerline.clear();
  446 +
  447 + // push back c2
  448 + for (size_t p = 0; p < id2 + 1; p++) { // first part
  449 + new_vertex = c2[p];
  450 + new_centerline.push_back(new_vertex);
  451 + }
  452 + result.push_back(new_centerline);
  453 + new_centerline.clear();
  454 +
  455 + for (size_t p = id2; p < num2; p++) { // second part
  456 + new_vertex = c2[p];
  457 + new_centerline.push_back(new_vertex);
  458 + }
  459 + result.push_back(new_centerline);
  460 + }
  461 + }
  462 + else { // if the stitch vertex is the first or last vertex of c2
  463 + // push back c2
  464 + if (id2 == 0) { // keep geometry information
  465 + new_vertex = c1[id1];
  466 + new_centerline.push_back(new_vertex);
  467 + for (size_t p = 0; p < num2; p++) { // stitch vertex on c1 -> geometry start vertex on c2 -> geometry end vertex on c2
  468 + new_vertex = c2[p];
  469 + new_centerline.push_back(new_vertex);
  470 + }
  471 + }
  472 + else {
  473 + for (size_t p = 0; p < num2; p++) { // geometry end vertex on c2 -> geometry start vertex on c2 -> stitch vertex on c1
  474 + new_vertex = c2[p];
  475 + new_centerline.push_back(new_vertex);
  476 + }
  477 + new_vertex = c1[id1];
  478 + new_centerline.push_back(new_vertex);
  479 + result.push_back(new_centerline);
  480 + new_centerline.clear();
  481 +
  482 + // push back c1
  483 + for (size_t p = 0; p < id1 + 1; p++) { // first part
  484 + new_vertex = c1[p];
  485 + new_centerline.push_back(new_vertex);
  486 + }
  487 + result.push_back(new_centerline);
  488 + new_centerline.clear();
  489 +
  490 + for (size_t p = id1; p < num1; p++) { // second part
  491 + new_vertex = c1[p];
  492 + new_centerline.push_back(new_vertex);
  493 + }
  494 + result.push_back(new_centerline);
  495 + }
  496 + }
  497 + }
  498 + }
  499 + }
  500 + return result;
  501 + }
  502 +
  503 + /// Split the fiber at the specified index. If the index is an end point, only one fiber is returned
  504 + std::vector< stim::centerline<T> > split(unsigned int idx){
  505 +
  506 + std::vector< stim::centerline<T> > fl; //create an array to store up to two fibers
  507 + size_t N = size();
  508 +
  509 + //if the index is an end point, only the existing fiber is returned
  510 + if(idx == 0 || idx == N-1){
  511 + fl.resize(1); //set the size of the fiber to 1
  512 + fl[0] = *this; //copy the current fiber
  513 + }
  514 +
  515 + //if the index is not an end point
  516 + else{
  517 +
  518 + unsigned int N1 = idx + 1; //calculate the size of both fibers
  519 + unsigned int N2 = N - idx;
  520 +
  521 + fl.resize(2); //set the array size to 2
  522 +
  523 + fl[0] = stim::centerline<T>(N1); //set the size of each fiber
  524 + fl[1] = stim::centerline<T>(N2);
  525 +
  526 + //copy both halves of the fiber
  527 + unsigned int i;
  528 +
  529 + //first half
  530 + for(i = 0; i < N1; i++) //for each centerline point
  531 + fl[0][i] = std::vector< stim::vec3<T> >::at(i);
  532 + fl[0].init(); //initialize the length vector
  533 +
  534 + //second half
  535 + for(i = 0; i < N2; i++)
  536 + fl[1][i] = std::vector< stim::vec3<T> >::at(idx+i);
  537 + fl[1].init(); //initialize the length vector
  538 + }
  539 +
  540 + return fl; //return the array
  541 +
  542 + }
  543 +
  544 + /// Outputs the fiber as a string
  545 + std::string str(){
  546 + std::stringstream ss;
  547 + size_t N = std::vector< stim::vec3<T> >::size();
  548 + ss << "---------[" << N << "]---------" << std::endl;
  549 + for (size_t i = 0; i < N; i++)
  550 + ss << std::vector< stim::vec3<T> >::at(i) << std::endl;
  551 + ss << "--------------------" << std::endl;
  552 +
  553 + return ss.str();
  554 + }
  555 +
  556 + /// Back method returns the last point in the fiber
  557 + stim::vec3<T> back(){
  558 + return std::vector< stim::vec3<T> >::back();
  559 + }
  560 +
  561 + ////resample a fiber in the network
  562 + stim::centerline<T> resample(T spacing)
  563 + {
  564 + //std::cout<<"fiber::resample()"<<std::endl;
  565 +
  566 + stim::vec3<T> v; //v-direction vector of the segment
  567 + stim::vec3<T> p; //- intermediate point to be added
  568 + stim::vec3<T> p1; // p1 - starting point of an segment on the fiber,
  569 + stim::vec3<T> p2; // p2 - ending point,
  570 + //double sum=0; //distance summation
  571 +
  572 + size_t N = size();
  573 +
  574 + centerline<T> new_c; // initialize list of new resampled points on the fiber
  575 + // for each point on the centerline (skip if it is the last point on centerline)
  576 + for(unsigned int f=0; f< N-1; f++)
  577 + {
  578 + p1 = at(f);
  579 + p2 = at(f+1);
  580 + v = p2 - p1;
  581 +
  582 + T lengthSegment = v.len(); //find Length of the segment as distance between the starting and ending points of the segment
  583 +
  584 + if(lengthSegment >= spacing){ // if length of the segment is greater than standard deviation resample
  585 +
  586 + // repeat resampling until accumulated stepsize is equsl to length of the segment
  587 + for(T step=0.0; step<lengthSegment; step+=spacing){
  588 + // calculate the resampled point by travelling step size in the direction of normalized gradient vector
  589 + p = p1 + v * (step / lengthSegment);
  590 +
  591 + // add this resampled points to the new fiber list
  592 + new_c.push_back(p);
  593 + }
  594 + }
  595 + 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
  596 + new_c.push_back(at(f));
  597 + }
  598 + new_c.push_back(at(N-1)); //add the last point on the fiber to the new fiber list
  599 + //centerline newFiber(newPointList);
  600 + return new_c;
  601 + }
  602 +
  603 +};
  604 +
  605 +
  606 +
  607 +} //end namespace stim
  608 +
  609 +
  610 +
  611 +#endif
cylinder.h 0 → 100644
  1 +#ifndef STIM_CYLINDER_H
  2 +#define STIM_CYLINDER_H
  3 +#include <iostream>
  4 +#include <stim/math/circle.h>
  5 +#include "centerline.h"
  6 +#include <stim/visualization/obj.h>
  7 +
  8 +
  9 +namespace stim
  10 +{
  11 +template<typename T>
  12 +class cylinder : public centerline<T> {
  13 +protected:
  14 +
  15 + using stim::centerline<T>::d;
  16 +
  17 + std::vector< stim::vec3<T> > U; //stores the array of U vectors defining the Frenet frame
  18 + std::vector< T > R; //stores a list of magnitudes for each point in the centerline (assuming mags[0] is the radius)
  19 +
  20 + using stim::centerline<T>::findIdx;
  21 +
  22 + void init() {
  23 + U.resize(size()); //allocate space for the frenet frame vectors
  24 +// if (R.size() != 0)
  25 + R.resize(size());
  26 +
  27 + stim::circle<T> c; //create a circle
  28 + stim::vec3<T> d0, d1;
  29 + for (size_t i = 0; i < size() - 1; i++) { //for each line segment in the centerline
  30 + c.rotate(d(i)); //rotate the circle to match that normal
  31 + U[i] = c.U; //save the U vector from the circle
  32 + }
  33 + U[size() - 1] = c.U; //for the last point, duplicate the final frenet frame vector
  34 + }
  35 +
  36 + //calculates the U values for each point to initialize the frenet frame
  37 + // this function assumes that the centerline has already been set
  38 +
  39 +public:
  40 +
  41 + using stim::centerline<T>::size;
  42 + using stim::centerline<T>::at;
  43 + using stim::centerline<T>::L;
  44 + using stim::centerline<T>::length;
  45 +
  46 + cylinder() : centerline<T>(){}
  47 +
  48 + cylinder(centerline<T>c) : centerline<T>(c) {
  49 + init();
  50 + }
  51 +
  52 + cylinder(std::vector<stim::vec3<T> > p, std::vector<T> s)
  53 + : centerline<T>(p)
  54 + {
  55 + R = s;
  56 + init();
  57 + }
  58 +
  59 + cylinder(stim::centerline<T> p, std::vector<T> s)
  60 + {
  61 + //was d = s;
  62 + p = s;
  63 + init();
  64 + }
  65 +
  66 + //cylinder(centerline<T>c, T r) : centerline(c) {
  67 + // init();
  68 + // //add_mag(r);
  69 + //}
  70 +
  71 + //initialize a cylinder with a list of points and magnitude values
  72 + //cylinder(centerline<T>c, std::vector<T> r) : centerline(c){
  73 + // init();
  74 + // //add_mag(r);
  75 + //}
  76 +
  77 + //copy the original radius
  78 + void copy_r(std::vector<T> radius) {
  79 + for (unsigned i = 0; i < radius.size(); i++)
  80 + R[i] = radius[i];
  81 + }
  82 +
  83 + ///Returns magnitude i at the given length into the fiber (based on the pvalue).
  84 + ///Interpolates the position along the line.
  85 + ///@param l: the location of the in the cylinder.
  86 + ///@param idx: integer location of the point closest to l but prior to it.
  87 + T r(T l, int idx) {
  88 + T a = (l - L[idx]) / (L[idx + 1] - L[idx]);
  89 + T v2 = (R[idx] + (R[idx + 1] - R[idx])*a);
  90 +
  91 + return v2;
  92 + }
  93 +
  94 + ///Returns the ith magnitude at the given p-value (p value ranges from 0 to 1).
  95 + ///interpolates the radius along the line.
  96 + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
  97 + T rl(T pvalue) {
  98 + if (pvalue <= 0.0) return R[0];
  99 + if (pvalue >= 1.0) return R[size() - 1];
  100 +
  101 + T l = pvalue*L[L.size() - 1];
  102 + int idx = findIdx(l);
  103 + return r(l, idx);
  104 + }
  105 +
  106 + /// Returns the magnitude at the given index
  107 + /// @param i is the index of the desired point
  108 + /// @param r is the index of the magnitude value
  109 + T r(unsigned i) {
  110 + return R[i];
  111 + }
  112 +
  113 +
  114 + ///adds a magnitude to each point in the cylinder
  115 + /*void add_mag(V val = 0) {
  116 + if (M.size() == 0) M.resize(size()); //if the magnitude vector isn't initialized, resize it to match the centerline
  117 + for (size_t i = 0; i < size(); i++) //for each point
  118 + R[i].push_back(val); //add this value to the magnitude vector at each point
  119 + }*/
  120 +
  121 + //adds a magnitude based on a list of magnitudes for each point
  122 + /*void add_mag(std::vector<T> val) {
  123 + if (M.size() == 0) M.resize(size()); //if the magnitude vector isn't initialized, resize it to match the centerline
  124 + for (size_t i = 0; i < size(); i++) //for each point
  125 + R[i].push_back(val[i]); //add this value to the magnitude vector at each point
  126 + }*/
  127 +
  128 + //sets the value of magnitude m at point i
  129 + void set_r(size_t i, T r) {
  130 + R[i] = r;
  131 + }
  132 +
  133 + /*size_t nmags() {
  134 + if (M.size() == 0) return 0;
  135 + else return R[0].size();
  136 + }*/
  137 +
  138 + ///Returns a circle representing the cylinder cross section at point i
  139 + stim::circle<T> circ(size_t i) {
  140 + return stim::circle<T>(at(i), R[i], d(i), U[i]);
  141 + }
  142 +
  143 + ///Returns an OBJ object representing the cylinder with a radial tesselation value of N using magnitude m
  144 + stim::obj<T> OBJ(size_t N) {
  145 + stim::obj<T> out; //create an OBJ object
  146 + stim::circle<T> c0, c1;
  147 + std::vector< stim::vec3<T> > p0, p1; //allocate space for the point sets representing the circles bounding each cylinder segment
  148 + T u0, u1, v0, v1; //allocate variables to store running texture coordinates
  149 + for (size_t i = 1; i < size(); i++) { //for each line segment in the cylinder
  150 + c0 = circ(i - 1); //get the two circles bounding the line segment
  151 + c1 = circ(i);
  152 +
  153 + p0 = c0.points(N); //get t points for each of the end caps
  154 + p1 = c1.points(N);
  155 +
  156 + u0 = L[i - 1] / length(); //calculate the texture coordinate (u, v) where u runs along the cylinder length
  157 + u1 = L[i] / length();
  158 +
  159 + for (size_t n = 1; n < N; n++) { //for each point in the circle
  160 + v0 = (double)(n-1) / (double)(N - 1); //v texture coordinate runs around the cylinder
  161 + v1 = (double)(n) / (double)(N - 1);
  162 + out.Begin(OBJ_FACE); //start a face (quad)
  163 + out.TexCoord(u0, v0);
  164 + out.Vertex(p0[n - 1][0], p0[n - 1][1], p0[n - 1][2]); //output the points composing a strip of quads wrapping the cylinder segment
  165 + out.TexCoord(u1, v0);
  166 + out.Vertex(p1[n - 1][0], p1[n - 1][1], p1[n - 1][2]);
  167 +
  168 + out.TexCoord(u0, v1);
  169 + out.Vertex(p1[n][0], p1[n][1], p1[n][2]);
  170 + out.TexCoord(u1, v1);
  171 + out.Vertex(p0[n][0], p0[n][1], p0[n][2]);
  172 + out.End();
  173 + }
  174 + v0 = (double)(N - 2) / (double)(N - 1); //v texture coordinate runs around the cylinder
  175 + v1 = 1.0;
  176 + out.Begin(OBJ_FACE);
  177 + out.TexCoord(u0, v0);
  178 + out.Vertex(p0[N - 1][0], p0[N - 1][1], p0[N - 1][2]); //output the points composing a strip of quads wrapping the cylinder segment
  179 + out.TexCoord(u1, v0);
  180 + out.Vertex(p1[N - 1][0], p1[N - 1][1], p1[N - 1][2]);
  181 +
  182 + out.TexCoord(u0, v1);
  183 + out.Vertex(p1[0][0], p1[0][1], p1[0][2]);
  184 + out.TexCoord(u1, v1);
  185 + out.Vertex(p0[0][0], p0[0][1], p0[0][2]);
  186 + out.End();
  187 + }
  188 + return out;
  189 + }
  190 +
  191 + std::string str() {
  192 + std::stringstream ss;
  193 + size_t N = std::vector< stim::vec3<T> >::size();
  194 + ss << "---------[" << N << "]---------" << std::endl;
  195 + for (size_t i = 0; i < N; i++)
  196 + ss << std::vector< stim::vec3<T> >::at(i) << " r = " << R[i] << " u = " << U[i] << std::endl;
  197 + ss << "--------------------" << std::endl;
  198 +
  199 + return ss.str();
  200 + }
  201 +
  202 + /// Integrates a magnitude value along the cylinder.
  203 + /// @param m is the magnitude value to be integrated (this is usually the radius)
  204 + T integrate() {
  205 + T sum = 0; //initialize the integral to zero
  206 + if (L.size() == 1)
  207 + return sum;
  208 + else {
  209 + T m0, m1; //allocate space for both magnitudes in a single segment
  210 + m0 = R[0]; //initialize the first point and magnitude to the first point in the cylinder
  211 + T len = L[1]; //allocate space for the segment length
  212 +
  213 +
  214 + for (unsigned p = 1; p < size(); p++) { //for every consecutive point in the cylinder
  215 + m1 = R[p];
  216 + if (p > 1) len = (L[p] - L[p - 1]); //calculate the segment length using the L array
  217 + sum += (m0 + m1) / (T)2.0 * len; //add the average magnitude, weighted by the segment length
  218 + m0 = m1; //move to the next segment by shifting points
  219 + }
  220 + return sum; //return the integral
  221 + }
  222 + }
  223 +
  224 + /// Resamples the cylinder to provide a maximum distance of "spacing" between centerline points. All current
  225 + /// centerline points are guaranteed to exist in the new cylinder
  226 + /// @param spacing is the maximum spacing allowed between sample points
  227 + cylinder<T> resample(T spacing) {
  228 + cylinder<T> c = stim::centerline<T>::resample(spacing); //resample the centerline and use it to create a new cylinder
  229 +
  230 + //size_t nm = nmags(); //get the number of magnitude values in the current cylinder
  231 + //if (nm > 0) { //if there are magnitude values
  232 + // std::vector<T> magvec(nm, 0); //create a magnitude vector for a single point
  233 + // c.M.resize(c.size(), magvec); //allocate space for a magnitude vector at each point of the new cylinder
  234 + //}
  235 +
  236 + T l, t;
  237 + for (size_t i = 0; i < c.size(); i++) { //for each point in the new cylinder
  238 + l = c.L[i]; //get the length along the new cylinder
  239 + t = l / length(); //calculate the parameter value along the new cylinder
  240 + //for (size_t mag = 0; mag < nm; mag++) { //for each magnitude value
  241 + c.R[i] = r(t); //retrieve the interpolated magnitude from the current cylinder and store it in the new one
  242 + //}
  243 + }
  244 + return c;
  245 + }
  246 +
  247 + std::vector< stim::cylinder<T> > split(unsigned int idx) {
  248 +
  249 + unsigned N = size();
  250 + std::vector< stim::centerline<T> > LL;
  251 + LL.resize(2);
  252 + LL = (*this).centerline<T>::split(idx);
  253 + std::vector< stim::cylinder<T> > C(LL.size());
  254 + unsigned i = 0;
  255 + C[0] = LL[0];
  256 + //C[0].R.resize(idx);
  257 + for (; i < idx + 1; i++) {
  258 + //for(unsigned d = 0; d < 3; d++)
  259 + //C[0][i][d] = LL[0].c[i][d];
  260 + C[0].R[i] = R[i];
  261 + //C[0].R[i].resize(1);
  262 + }
  263 + if (C.size() == 2) {
  264 + C[1] = LL[1];
  265 + i--;
  266 + //C[1].M.resize(N - idx);
  267 + for (; i < N; i++) {
  268 + //for(unsigned d = 0; d < 3; d++)
  269 + //C[1][i - idx][d] = LL[1].c[i - idx][d];
  270 + C[1].R[i - idx] = R[i];
  271 + //C[1].M[i - idx].resize(1);
  272 + }
  273 + }
  274 +
  275 + return C;
  276 + }
  277 +
  278 +
  279 + /*
  280 + ///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and magnitudes (inM)
  281 + void
  282 + init(centerline inP, std::vector< std::vector<T> > inM) {
  283 + M = inM; //the magnitude vector can be copied directly
  284 + (*this) = inP; //the centerline can be copied to this class directly
  285 + stim::vec3<float> v1;
  286 + stim::vec3<float> v2;
  287 + e.resize(inP.size());
  288 +
  289 + norms.resize(inP.size());
  290 + Us.resize(inP.size());
  291 +
  292 + if(inP.size() < 2)
  293 + return;
  294 +
  295 + //calculate each L.
  296 + L.resize(inP.size()); //the number of precomputed lengths will equal the number of points
  297 + T temp = (T)0; //length up to that point
  298 + L[0] = temp;
  299 + for(size_t i = 1; i < L.size(); i++)
  300 + {
  301 + temp += (inP[i-1] - inP[i]).len();
  302 + L[i] = temp;
  303 + }
  304 +
  305 + stim::vec3<T> dr = (inP[1] - inP[0]).norm();
  306 + s = stim::circle<T>(inP[0], inR[0][0], dr, stim::vec3<T>(1,0,0));
  307 + norms[0] = s.N;
  308 + e[0] = s;
  309 + Us[0] = e[0].U;
  310 + for(size_t i = 1; i < inP.size()-1; i++)
  311 + {
  312 + s.center(inP[i]);
  313 + v1 = (inP[i] - inP[i-1]).norm();
  314 + v2 = (inP[i+1] - inP[i]).norm();
  315 + dr = (v1+v2).norm();
  316 + s.normal(dr);
  317 +
  318 + norms[i] = s.N;
  319 +
  320 + s.scale(inR[i][0]/inR[i-1][0]);
  321 + e[i] = s;
  322 + Us[i] = e[i].U;
  323 + }
  324 +
  325 + int j = inP.size()-1;
  326 + s.center(inP[j]);
  327 + dr = (inP[j] - inP[j-1]).norm();
  328 + s.normal(dr);
  329 +
  330 + norms[j] = s.N;
  331 +
  332 + s.scale(inR[j][0]/inR[j-1][0]);
  333 + e[j] = s;
  334 + Us[j] = e[j].U;
  335 + }
  336 +
  337 + ///returns the direction vector at point idx.
  338 + stim::vec3<T>
  339 + d(int idx)
  340 + {
  341 + if(idx == 0)
  342 + {
  343 + stim::vec3<T> temp(
  344 + c[idx+1][0]-c[idx][0],
  345 + c[idx+1][1]-c[idx][1],
  346 + c[idx+1][2]-c[idx][2]
  347 + );
  348 +// return (e[idx+1].P - e[idx].P).norm();
  349 + return (temp.norm());
  350 + }
  351 + else if(idx == N-1)
  352 + {
  353 + stim::vec3<T> temp(
  354 + c[idx][0]-c[idx+1][0],
  355 + c[idx][1]-c[idx+1][1],
  356 + c[idx][2]-c[idx+1][2]
  357 + );
  358 + // return (e[idx].P - e[idx-1].P).norm();
  359 + return (temp.norm());
  360 + }
  361 + else
  362 + {
  363 +// return (e[idx+1].P - e[idx].P).norm();
  364 +// stim::vec3<float> v1 = (e[idx].P-e[idx-1].P).norm();
  365 + stim::vec3<T> v1(
  366 + c[idx][0]-c[idx-1][0],
  367 + c[idx][1]-c[idx-1][1],
  368 + c[idx][2]-c[idx-1][2]
  369 + );
  370 +
  371 +// stim::vec3<float> v2 = (e[idx+1].P-e[idx].P).norm();
  372 + stim::vec3<T> v2(
  373 + c[idx+1][0]-c[idx][0],
  374 + c[idx+1][1]-c[idx][1],
  375 + c[idx+1][2]-c[idx][2]
  376 + );
  377 +
  378 + return (v1.norm()+v2.norm()).norm();
  379 + }
  380 + // return e[idx].N;
  381 +
  382 + }
  383 +
  384 + stim::vec3<T>
  385 + d(T l, int idx)
  386 + {
  387 + if(idx == 0 || idx == N-1)
  388 + {
  389 + return norms[idx];
  390 +// return e[idx].N;
  391 + }
  392 + else
  393 + {
  394 +
  395 + T rat = (l-L[idx])/(L[idx+1]-L[idx]);
  396 + return( norms[idx] + (norms[idx+1] - norms[idx])*rat);
  397 +// return( e[idx].N + (e[idx+1].N - e[idx].N)*rat);
  398 + }
  399 + }
  400 +
  401 +
  402 + ///finds the index of the point closest to the length l on the lower bound.
  403 + ///binary search.
  404 + int
  405 + findIdx(T l)
  406 + {
  407 + unsigned int i = L.size()/2;
  408 + unsigned int max = L.size()-1;
  409 + unsigned int min = 0;
  410 + while(i > 0 && i < L.size()-1)
  411 + {
  412 +// std::cerr << "Trying " << i << std::endl;
  413 +// std::cerr << "l is " << l << ", L[" << i << "]" << L[i] << std::endl;
  414 + if(l < L[i])
  415 + {
  416 + max = i;
  417 + i = min+(max-min)/2;
  418 + }
  419 + else if(L[i] <= l && L[i+1] >= l)
  420 + {
  421 + break;
  422 + }
  423 + else
  424 + {
  425 + min = i;
  426 + i = min+(max-min)/2;
  427 + }
  428 + }
  429 + return i;
  430 + }
  431 +
  432 + public:
  433 + ///default constructor
  434 + cylinder()
  435 + // : centerline<T>()
  436 + {
  437 +
  438 + }
  439 +
  440 + ///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder.
  441 + ///@param inP: Vector of stim vec3 composing the points of the centerline.
  442 + ///@param inM: Vector of stim vecs composing the radii of the centerline.
  443 + cylinder(std::vector<stim::vec3<T> > inP, std::vector<stim::vec<T> > inM)
  444 + : centerline<T>(inP)
  445 + {
  446 + init(inP, inM);
  447 + }
  448 +
  449 + ///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder.
  450 + ///@param inP: Vector of stim vec3 composing the points of the centerline.
  451 + ///@param inM: Vector of stim vecs composing the radii of the centerline.
  452 + cylinder(std::vector<stim::vec3<T> > inP, std::vector< T > inM)
  453 + : centerline<T>(inP)
  454 + {
  455 + std::vector<stim::vec<T> > temp;
  456 + stim::vec<T> zero(0.0,0.0);
  457 + temp.resize(inM.size(), zero);
  458 + for(int i = 0; i < inM.size(); i++)
  459 + temp[i][0] = inR[i];
  460 + init(inP, temp);
  461 + }
  462 +
  463 +
  464 + ///Constructor defines a cylinder with centerline inP and magnitudes of zero
  465 + ///@param inP: Vector of stim vec3 composing the points of the centerline
  466 + cylinder(std::vector< stim::vec3<T> > inP)
  467 + : centerline<T>(inP)
  468 + {
  469 + std::vector< stim::vec<T> > inM; //create an array of arbitrary magnitudes
  470 +
  471 + stim::vec<T> zero;
  472 + zero.push_back(0);
  473 +
  474 + inM.resize(inP.size(), zero); //initialize the magnitude values to zero
  475 + init(inP, inM);
  476 + }
  477 +
  478 + //assignment operator creates a cylinder from a centerline (default radius is zero)
  479 + cylinder& operator=(stim::centerline<T> c) {
  480 + init(c);
  481 + }
  482 +
  483 + ///Returns the number of points on the cylinder centerline
  484 +
  485 + unsigned int size(){
  486 + return N;
  487 + }
  488 +
  489 +
  490 + ///Returns a position vector at the given p-value (p value ranges from 0 to 1).
  491 + ///interpolates the position along the line.
  492 + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
  493 + stim::vec3<T>
  494 + p(T pvalue)
  495 + {
  496 + if(pvalue < 0.0 || pvalue > 1.0)
  497 + {
  498 + return stim::vec3<float>(-1,-1,-1);
  499 + }
  500 + T l = pvalue*L[L.size()-1];
  501 + int idx = findIdx(l);
  502 + return (p(l,idx));
  503 + }
  504 +
  505 + ///Returns a position vector at the given length into the fiber (based on the pvalue).
  506 + ///Interpolates the radius along the line.
  507 + ///@param l: the location of the in the cylinder.
  508 + ///@param idx: integer location of the point closest to l but prior to it.
  509 + stim::vec3<T>
  510 + p(T l, int idx)
  511 + {
  512 + T rat = (l-L[idx])/(L[idx+1]-L[idx]);
  513 + stim::vec3<T> v1(
  514 + c[idx][0],
  515 + c[idx][1],
  516 + c[idx][2]
  517 + );
  518 +
  519 + stim::vec3<T> v2(
  520 + c[idx+1][0],
  521 + c[idx+1][1],
  522 + c[idx+1][2]
  523 + );
  524 +// return( e[idx].P + (e[idx+1].P-e[idx].P)*rat);
  525 + return( v1 + (v2-v1)*rat);
  526 +// return(
  527 +// return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
  528 + }
  529 +
  530 + ///Returns a radius at the given p-value (p value ranges from 0 to 1).
  531 + ///interpolates the radius along the line.
  532 + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
  533 + T
  534 + r(T pvalue)
  535 + {
  536 + if(pvalue < 0.0 || pvalue > 1.0){
  537 + std::cerr<<"Error, value "<<pvalue<<" is outside of [0 1]."<<std::endl;
  538 + exit(1);
  539 + }
  540 + T l = pvalue*L[L.size()-1];
  541 + int idx = findIdx(l);
  542 + return (r(l,idx));
  543 + }
  544 +
  545 + ///Returns a radius at the given length into the fiber (based on the pvalue).
  546 + ///Interpolates the position along the line.
  547 + ///@param l: the location of the in the cylinder.
  548 + ///@param idx: integer location of the point closest to l but prior to it.
  549 + T
  550 + r(T l, int idx)
  551 + {
  552 + T rat = (l-L[idx])/(L[idx+1]-L[idx]);
  553 + T v1 = (e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*rat);
  554 + T v3 = (Us[idx].len() + (Us[idx+1].len() - Us[idx].len())*rat);
  555 + T v2 = (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat);
  556 +// std::cout << (float)v1 = (float) v2 << std::endl;
  557 + if(abs(v3 - v1) >= 10e-6)
  558 + {
  559 + std::cout << "-------------------------" << std::endl;
  560 + std::cout << e[idx].str() << std::endl << std::endl;
  561 + std::cout << Us[idx].str() << std::endl;
  562 + std::cout << (float)v1 - (float) v2 << std::endl;
  563 + std::cout << "failed" << std::endl;
  564 + }
  565 +// std::cout << e[idx].U.len() << " " << mags[idx][0] << std::endl;
  566 +// std::cout << v2 << std::endl;
  567 + return(v2);
  568 +// return (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat);
  569 + // (
  570 + }
  571 +
  572 + /// Returns the magnitude at the given index
  573 + /// @param i is the index of the desired point
  574 + /// @param m is the index of the magnitude value
  575 + T ri(unsigned i, unsigned m = 0){
  576 + return mags[i][m];
  577 + }
  578 +
  579 + /// Adds a new magnitude value to all points
  580 + /// @param m is the starting value for the new magnitude
  581 + void add_mag(T m = 0){
  582 + for(unsigned int p = 0; p < N; p++)
  583 + mags[p].push_back(m);
  584 + }
  585 +
  586 + /// Sets a magnitude value
  587 + /// @param val is the new value for the magnitude
  588 + /// @param p is the point index for the magnitude to be set
  589 + /// @param m is the index for the magnitude
  590 + void set_mag(T val, unsigned p, unsigned m = 0){
  591 + mags[p][m] = val;
  592 + }
  593 +
  594 + /// Returns the number of magnitude values at each point
  595 + unsigned nmags(){
  596 + return mags[0].size();
  597 + }
  598 +
  599 + ///returns the position of the point with a given pvalue and theta on the surface
  600 + ///in x, y, z coordinates. Theta is in degrees from 0 to 360.
  601 + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
  602 + ///@param theta: the angle to the point of a circle.
  603 + stim::vec3<T>
  604 + surf(T pvalue, T theta)
  605 + {
  606 + if(pvalue < 0.0 || pvalue > 1.0)
  607 + {
  608 + return stim::vec3<float>(-1,-1,-1);
  609 + } else {
  610 + T l = pvalue*L[L.size()-1];
  611 + int idx = findIdx(l);
  612 + stim::vec3<T> ps = p(l, idx);
  613 + T m = r(l, idx);
  614 + s = e[idx];
  615 + s.center(ps);
  616 + s.normal(d(l, idx));
  617 + s.scale(m/e[idx].U.len());
  618 + return(s.p(theta));
  619 + }
  620 + }
  621 +
  622 + ///returns a vector of points necessary to create a circle at every position in the fiber.
  623 + ///@param sides: the number of sides of each circle.
  624 + std::vector<std::vector<vec3<T> > >
  625 + getPoints(int sides)
  626 + {
  627 + std::vector<std::vector <vec3<T> > > points;
  628 + points.resize(N);
  629 + for(int i = 0; i < N; i++)
  630 + {
  631 + points[i] = e[i].getPoints(sides);
  632 + }
  633 + return points;
  634 + }
  635 +
  636 + ///returns the total length of the line at index j.
  637 + T
  638 + getl(int j)
  639 + {
  640 + return (L[j]);
  641 + }
  642 + /// Allows a point on the centerline to be accessed using bracket notation
  643 +
  644 + vec3<T> operator[](unsigned int i){
  645 + return e[i].P;
  646 + }
  647 +
  648 + /// Returns the total length of the cylinder centerline
  649 + T length(){
  650 + return L.back();
  651 + }
  652 +
  653 + /// Integrates a magnitude value along the cylinder.
  654 + /// @param m is the magnitude value to be integrated (this is usually the radius)
  655 + T integrate(unsigned m = 0){
  656 +
  657 + T M = 0; //initialize the integral to zero
  658 + T m0, m1; //allocate space for both magnitudes in a single segment
  659 +
  660 + //vec3<T> p0, p1; //allocate space for both points in a single segment
  661 +
  662 + m0 = mags[0][m]; //initialize the first point and magnitude to the first point in the cylinder
  663 + //p0 = pos[0];
  664 +
  665 + T len = L[0]; //allocate space for the segment length
  666 +
  667 + //for every consecutive point in the cylinder
  668 + for(unsigned p = 1; p < N; p++){
  669 +
  670 + //p1 = pos[p]; //get the position and magnitude for the next point
  671 + m1 = mags[p][m];
  672 +
  673 + if(p > 1) len = (L[p-1] - L[p-2]); //calculate the segment length using the L array
  674 +
  675 + //add the average magnitude, weighted by the segment length
  676 + M += (m0 + m1)/(T)2.0 * len;
  677 +
  678 + m0 = m1; //move to the next segment by shifting points
  679 + }
  680 + return M; //return the integral
  681 + }
  682 +
  683 + /// Averages a magnitude value across the cylinder
  684 + /// @param m is the magnitude value to be averaged (this is usually the radius)
  685 + T average(unsigned m = 0){
  686 +
  687 + //return the average magnitude
  688 + return integrate(m) / L.back();
  689 + }
  690 +
  691 + /// Resamples the cylinder to provide a maximum distance of "spacing" between centerline points. All current
  692 + /// centerline points are guaranteed to exist in the new cylinder
  693 + /// @param spacing is the maximum spacing allowed between sample points
  694 + cylinder<T> resample(T spacing){
  695 +
  696 + std::vector< vec3<T> > result;
  697 +
  698 + vec3<T> p0 = e[0].P; //initialize p0 to the first point on the centerline
  699 + vec3<T> p1;
  700 + unsigned N = size(); //number of points in the current centerline
  701 +
  702 + //for each line segment on the centerline
  703 + for(unsigned int i = 1; i < N; i++){
  704 + p1 = e[i].P; //get the second point in the line segment
  705 +
  706 + vec3<T> v = p1 - p0; //calculate the vector between these two points
  707 + T d = v.len(); //calculate the distance between these two points (length of the line segment)
  708 +
  709 + size_t nsteps = (size_t)std::ceil(d / spacing); //calculate the number of steps to take along the segment to meet the spacing criteria
  710 + T stepsize = (T)1.0 / nsteps; //calculate the parametric step size between new centerline points
  711 +
  712 + //for each step along the line segment
  713 + for(unsigned s = 0; s < nsteps; s++){
  714 + T alpha = stepsize * s; //calculate the fraction of the distance along the line segment covered
  715 + result.push_back(p0 + alpha * v); //push the point at alpha position along the line segment
  716 + }
  717 +
  718 + p0 = p1; //shift the points to move to the next line segment
  719 + }
  720 +
  721 + result.push_back(e[size() - 1].P); //push the last point in the centerline
  722 +
  723 + return cylinder<T>(result);
  724 +
  725 + }*/
  726 +
  727 +
  728 +};
  729 +
  730 +}
  731 +#endif
gl_network.h 0 → 100644
  1 +#ifndef STIM_GL_NETWORK
  2 +#define STIM_GL_NETWORK
  3 +
  4 +#include <GL/glut.h>
  5 +#include "network.h"
  6 +#include <stim/visualization/aaboundingbox.h>
  7 +
  8 +namespace stim{
  9 +
  10 +template <typename T>
  11 +class gl_network : public stim::network<T>{
  12 +
  13 +protected:
  14 + using stim::network<T>::E;
  15 + using stim::network<T>::V;
  16 +
  17 + GLuint dlist;
  18 +
  19 +public:
  20 +
  21 + /// Default constructor
  22 + gl_network() : stim::network<T>(){
  23 + dlist = 0;
  24 + }
  25 +
  26 + /// Constructor creates a gl_network from a stim::network
  27 + gl_network(stim::network<T> N) : stim::network<T>(N){
  28 + dlist = 0;
  29 + }
  30 +
  31 + /// Fills the parameters with the minimum and maximum spatial positions in the network,
  32 + /// specifying a bounding box for the network geometry
  33 + aaboundingbox<T> boundingbox(){
  34 +
  35 + aaboundingbox<T> bb; //create a bounding box
  36 +
  37 + //loop through every edge
  38 + for(unsigned e = 0; e < E.size(); e++){
  39 + //loop through every point
  40 + for(unsigned p = 0; p < E[e].size(); p++)
  41 + bb.expand(E[e][p]); //expand the bounding box to include the point
  42 + }
  43 +
  44 + return bb; //return the bounding box
  45 + }
  46 +
  47 + ///render cylinder based on points from the top/bottom hat
  48 + ///@param C1 set of points from one of the hat
  49 + void renderCylinder(std::vector< stim::vec3<T> > C1, std::vector< stim::vec3<T> > C2, stim::vec3<T> P1, stim::vec3<T> P2) {
  50 + glBegin(GL_QUAD_STRIP);
  51 + for (unsigned i = 0; i < C1.size(); i++) { // for every point on the circle
  52 + stim::vec3<T> n1 = C1[i] - P1; // set normal vector for every vertex on the quads
  53 + stim::vec3<T> n2 = C2[i] - P2;
  54 + n1 = n1.norm();
  55 + n2 = n2.norm();
  56 + glNormal3f(n1[0], n1[1], n1[2]);
  57 + glVertex3f(C1[i][0], C1[i][1], C1[i][2]);
  58 + glNormal3f(n2[0], n2[1], n2[2]);
  59 + glVertex3f(C2[i][0], C2[i][1], C2[i][2]);
  60 + }
  61 + glEnd();
  62 + }
  63 +
  64 + ///render the vertex as sphere based on glut build-in function
  65 + ///@param x, y, z are the three coordinates of the center point
  66 + ///@param radius is the radius of the sphere
  67 + ///@param subdivisions is the slice/stride along/around z-axis
  68 + void renderBall(T x, T y, T z, T radius, int subdivisions) {
  69 + glPushMatrix();
  70 + glTranslatef(x, y, z);
  71 + glutSolidSphere(radius, subdivisions, subdivisions);
  72 + glPopMatrix();
  73 + }
  74 +
  75 + ///render the vertex as sphere based on transformation
  76 + ///@param x, y, z are the three coordinates of the center point
  77 + ///@param radius is the radius of the sphere
  78 + ///@param slice is the number of subdivisions around the z-axis
  79 + ///@param stack is the number of subdivisions along the z-axis
  80 + void renderBall(T x, T y, T z, T radius, T slice, T stack) {
  81 + T step_z = stim::PI / slice; // step angle along z-axis
  82 + T step_xy = 2 * stim::PI / stack; // step angle in xy-plane
  83 + T xx[4], yy[4], zz[4]; // store coordinates
  84 +
  85 + T angle_z = 0.0; // start angle
  86 + T angle_xy = 0.0;
  87 +
  88 + glBegin(GL_QUADS);
  89 + for (unsigned i = 0; i < slice; i++) { // around the z-axis
  90 + angle_z = i * step_z; // step step_z each time
  91 +
  92 + for (unsigned j = 0; j < stack; j++) { // along the z-axis
  93 + angle_xy = j * step_xy; // step step_xy each time, draw floor by floor
  94 +
  95 + xx[0] = radius * std::sin(angle_z) * std::cos(angle_xy); // four vertices
  96 + yy[0] = radius * std::sin(angle_z) * std::sin(angle_xy);
  97 + zz[0] = radius * std::cos(angle_z);
  98 +
  99 + xx[1] = radius * std::sin(angle_z + step_z) * std::cos(angle_xy);
  100 + yy[1] = radius * std::sin(angle_z + step_z) * std::sin(angle_xy);
  101 + zz[1] = radius * std::cos(angle_z + step_z);
  102 +
  103 + xx[2] = radius * std::sin(angle_z + step_z) * std::cos(angle_xy + step_xy);
  104 + yy[2] = radius * std::sin(angle_z + step_z) * std::sin(angle_xy + step_xy);
  105 + zz[2] = radius * std::cos(angle_z + step_z);
  106 +
  107 + xx[3] = radius * std::sin(angle_z) * std::cos(angle_xy + step_xy);
  108 + yy[3] = radius * std::sin(angle_z) * std::sin(angle_xy + step_xy);
  109 + zz[3] = radius * std::cos(angle_z);
  110 +
  111 + for (unsigned k = 0; k < 4; k++) {
  112 + glVertex3f(x + xx[k], y + yy[k], z + zz[k]); // draw the floor plane
  113 + }
  114 + }
  115 + }
  116 + glEnd();
  117 + }
  118 +
  119 + /// Render the network centerline as a series of line strips.
  120 + /// glCenterline0 is for only one input
  121 + void glCenterline0(){
  122 + if (!glIsList(dlist)) { //if dlist isn't a display list, create it
  123 + dlist = glGenLists(1); //generate a display list
  124 + glNewList(dlist, GL_COMPILE); //start a new display list
  125 + for (unsigned e = 0; e < E.size(); e++) { //for each edge in the network
  126 + glBegin(GL_LINE_STRIP);
  127 + for (unsigned p = 0; p < E[e].size(); p++) { //for each point on that edge
  128 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); //set the vertex position based on the current point
  129 + glTexCoord1f(0); //set white color
  130 + }
  131 + glEnd();
  132 + }
  133 + glEndList(); //end the display list
  134 + }
  135 + glCallList(dlist); // render the display list
  136 + }
  137 +
  138 + ///render the network centerline as a series of line strips(when loading at least two networks, otherwise using glCenterline0())
  139 + ///colors are based on metric values
  140 + void glCenterline(){
  141 +
  142 + if(!glIsList(dlist)){ //if dlist isn't a display list, create it
  143 + dlist = glGenLists(1); //generate a display list
  144 + glNewList(dlist, GL_COMPILE); //start a new display list
  145 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  146 + //unsigned errormag_id = E[e].nmags() - 1;
  147 + glBegin(GL_LINE_STRIP);
  148 + for(unsigned p = 0; p < E[e].size(); p++){ //for each point on that edge
  149 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); //set the vertex position based on the current point
  150 + glTexCoord1f(E[e].r(p)); //set the texture coordinate based on the specified magnitude index
  151 + }
  152 + glEnd();
  153 + }
  154 + glEndList(); //end the display list
  155 + }
  156 + glCallList(dlist); //render the display list
  157 + }
  158 +
  159 + ///render the network cylinder as a series of tubes(when only one network loaded)
  160 + void glCylinder0(T scale = 1.0f, bool undo = false) {
  161 +
  162 + float r1, r2;
  163 + if (undo == true)
  164 + glDeleteLists(dlist, 1); // delete display list
  165 + if (!glIsList(dlist)) { // if dlist isn't a display list, create it
  166 + dlist = glGenLists(1); // generate a display list
  167 + glNewList(dlist, GL_COMPILE); // start a new display list
  168 + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network
  169 + for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge
  170 + stim::circle<T> C1 = E[e].circ(p - 1);
  171 + stim::circle<T> C2 = E[e].circ(p);
  172 + r1 = E[e].r(p - 1);
  173 + r2 = E[e].r(p);
  174 + C1.set_R(scale * r1); // re-scale the circle to the same
  175 + C2.set_R(scale * r2);
  176 + std::vector< stim::vec3<T> > Cp1 = C1.glpoints(20); // get 20 points on the circle plane
  177 + std::vector< stim::vec3<T> > Cp2 = C2.glpoints(20);
  178 + glBegin(GL_QUAD_STRIP);
  179 + for (unsigned i = 0; i < Cp1.size(); i++) {
  180 + glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
  181 + glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
  182 + }
  183 + glEnd();
  184 + } // set the texture coordinate based on the specified magnitude index
  185 + }
  186 + for (unsigned n = 0; n < V.size(); n++) {
  187 + for (unsigned i = 0; i < E.size(); i++) {
  188 + if (E[i].v[0] == n) {
  189 + r1 = E[i].r(0) * scale;
  190 + break;
  191 + }
  192 + else if (E[i].v[1] == n) {
  193 + r1 = E[i].r(E[i].size() - 1) * scale;
  194 + break;
  195 + }
  196 + }
  197 + renderBall(V[n][0], V[n][1], V[n][2], r1, 20);
  198 + }
  199 + glEndList(); // end the display list
  200 + }
  201 + glCallList(dlist); // render the display list
  202 + }
  203 +
  204 + ///render the network cylinder as a series of tubes
  205 + ///colors are based on metric values
  206 + void glCylinder(float sigma, float radius) {
  207 +
  208 + if (radius != sigma) // if render radius was changed by user, create a new display list
  209 + glDeleteLists(dlist, 1);
  210 + if (!glIsList(dlist)) { // if dlist isn't a display list, create it
  211 + dlist = glGenLists(1); // generate a display list
  212 + glNewList(dlist, GL_COMPILE); // start a new display list
  213 + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network
  214 + for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge
  215 + stim::circle<T> C1 = E[e].circ(p - 1);
  216 + stim::circle<T> C2 = E[e].circ(p);
  217 + C1.set_R(2*radius); // re-scale the circle to the same
  218 + C2.set_R(2*radius);
  219 + std::vector< stim::vec3<T> > Cp1 = C1.glpoints(20);// get 20 points on the circle plane
  220 + std::vector< stim::vec3<T> > Cp2 = C2.glpoints(20);
  221 + glBegin(GL_QUAD_STRIP);
  222 + for (unsigned i = 0; i < Cp1.size(); i++) {
  223 + stim::vec3<T> n1 = Cp1[i] - E[e][p - 1]; // set normal vector for every vertex on the quads
  224 + stim::vec3<T> n2 = Cp2[i] - E[e][p];
  225 + n1 = n1.norm();
  226 + n2 = n2.norm();
  227 +
  228 + glNormal3f(n1[0], n1[1], n1[2]);
  229 + glTexCoord1f(E[e].r(p - 1));
  230 + glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
  231 + glNormal3f(n2[0], n2[1], n2[2]);
  232 + glTexCoord1f(E[e].r(p));
  233 + glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
  234 + }
  235 + glEnd();
  236 + } // set the texture coordinate based on the specified magnitude index
  237 + }
  238 + for (unsigned n = 0; n < V.size(); n++) {
  239 + size_t num = V[n].e[0].size(); // store the number of outgoing edge of that vertex
  240 + if (num != 0) { // if it has outgoing edge
  241 + unsigned idx = V[n].e[0][0]; // find the index of first outgoing edge of that vertex
  242 + glTexCoord1f(E[idx].r(0)); // bind the texture as metric of first point on that edge
  243 + }
  244 + else {
  245 + unsigned idx = V[n].e[1][0]; // find the index of first incoming edge of that vertex
  246 + glTexCoord1f(E[idx].r(E[idx].size() - 1)); // bind the texture as metric of last point on that edge
  247 + }
  248 + renderBall(V[n][0], V[n][1], V[n][2], 2*radius, 20);
  249 + }
  250 + glEndList(); // end the display list
  251 + }
  252 + glCallList(dlist); // render the display list
  253 + }
  254 +
  255 + ///render a T as a adjoint network of GT in transparancy(darkgreen, overlaid)
  256 + void glAdjointCylinder(float sigma, float radius) {
  257 +
  258 + if (radius != sigma) // if render radius was changed by user, create a new display list
  259 + glDeleteLists(dlist + 4, 1);
  260 + if (!glIsList(dlist + 4)) {
  261 + glNewList(dlist + 4, GL_COMPILE);
  262 + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network
  263 + for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge
  264 + stim::circle<T> C1 = E[e].circ(p - 1);
  265 + stim::circle<T> C2 = E[e].circ(p);
  266 + C1.set_R(2 * radius); // scale the circle to the same
  267 + C2.set_R(2 * radius);
  268 + std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
  269 + std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
  270 + glBegin(GL_QUAD_STRIP);
  271 + for (unsigned i = 0; i < Cp1.size(); i++) { // for every point on the circle(+1 means closing the circle)
  272 + glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
  273 + glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
  274 + }
  275 + glEnd();
  276 + } // set the texture coordinate based on the specified magnitude index
  277 + }
  278 + for (unsigned n = 0; n < V.size(); n++) {
  279 + size_t num = V[n].e[0].size(); // store the number of outgoing edge of that vertex
  280 + if (num != 0) { // if it has outgoing edge
  281 + unsigned idx = V[n].e[0][0]; // find the index of first outgoing edge of that vertex
  282 + }
  283 + else {
  284 + unsigned idx = V[n].e[1][0]; // find the index of first incoming edge of that vertex
  285 + }
  286 + renderBall(V[n][0], V[n][1], V[n][2], 2 * radius, 20);
  287 + }
  288 + glEndList();
  289 + }
  290 + glCallList(dlist + 4);
  291 + }
  292 +
  293 + ///render the network cylinder as series of tubes
  294 + ///@param I is a indicator: 0 -> GT, 1 -> T
  295 + ///@param map is the mapping relationship between two networks
  296 + ///@param colormap is the random generated color set for render
  297 + void glRandColorCylinder(int I, std::vector<unsigned> map, std::vector<T> colormap, float sigma, float radius) {
  298 +
  299 + if (radius != sigma) // if render radius was changed by user, create a new display list
  300 + glDeleteLists(dlist + 2, 1);
  301 + if (!glIsList(dlist + 2)) { // if dlist isn't a display list, create it
  302 + glNewList(dlist + 2, GL_COMPILE); // start a new display list
  303 + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network
  304 + if (map[e] != unsigned(-1)) {
  305 + if (I == 0) { // if it is to render GT
  306 + glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
  307 + }
  308 + else { // if it is to render T
  309 + glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
  310 + }
  311 +
  312 + for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
  313 + stim::circle<T> C1 = E[e].circ(p - 1);
  314 + stim::circle<T> C2 = E[e].circ(p);
  315 + C1.set_R(2*radius); // re-scale the circle to the same
  316 + C2.set_R(2*radius);
  317 + std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
  318 + std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
  319 + renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]);
  320 + }
  321 + }
  322 + else {
  323 + glColor3f(1.0, 1.0, 1.0); // white color for the un-mapping edges
  324 + for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
  325 + stim::circle<T> C1 = E[e].circ(p - 1);
  326 + stim::circle<T> C2 = E[e].circ(p);
  327 + C1.set_R(2*radius); // scale the circle to the same
  328 + C2.set_R(2*radius);
  329 + std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
  330 + std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
  331 + renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]);
  332 + }
  333 + }
  334 + }
  335 + for (unsigned n = 0; n < V.size(); n++) {
  336 + size_t num_edge = V[n].e[0].size() + V[n].e[1].size();
  337 + if (num_edge > 1) { // if it is the joint vertex
  338 + glColor3f(0.3, 0.3, 0.3); // gray
  339 + renderBall(V[n][0], V[n][1], V[n][2], 3*radius, 20);
  340 + }
  341 + else { // if it is the terminal vertex
  342 + glColor3f(0.6, 0.6, 0.6); // more white gray
  343 + renderBall(V[n][0], V[n][1], V[n][2], 3*radius, 20);
  344 + }
  345 + }
  346 + glEndList();
  347 + }
  348 + glCallList(dlist + 2);
  349 + }
  350 +
  351 + ///render the network centerline as a series of line strips in random different color
  352 + ///@param I is a indicator: 0 -> GT, 1 -> T
  353 + ///@param map is the mapping relationship between two networks
  354 + ///@param colormap is the random generated color set for render
  355 + void glRandColorCenterline(int I, std::vector<unsigned> map, std::vector<T> colormap) {
  356 + if (!glIsList(dlist + 2)) {
  357 + glNewList(dlist + 2, GL_COMPILE);
  358 + for (unsigned e = 0; e < E.size(); e++) {
  359 + if (map[e] != unsigned(-1)) { // if it has corresponding edge in another network
  360 + if (I == 0) // if it is to render GT
  361 + glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
  362 + else // if it is to render T
  363 + glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
  364 +
  365 + glBegin(GL_LINE_STRIP);
  366 + for (unsigned p = 0; p < E[e].size(); p++) {
  367 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  368 + }
  369 + glEnd();
  370 + }
  371 + else {
  372 + glColor3f(1.0, 1.0, 1.0); // white color for the un-mapping edges
  373 + glBegin(GL_LINE_STRIP);
  374 + for (unsigned p = 0; p < E[e].size(); p++) {
  375 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  376 + }
  377 + glEnd();
  378 + }
  379 + }
  380 + glEndList();
  381 + }
  382 + glCallList(dlist + 2);
  383 + }
  384 +
  385 + void glAdjointCenterline() {
  386 + if (!glIsList(dlist + 4)) {
  387 + glNewList(dlist + 4, GL_COMPILE);
  388 + for (unsigned e = 0; e < E.size(); e++) { //for each edge in the network
  389 +
  390 + glBegin(GL_LINE_STRIP);
  391 + for (unsigned p = 0; p < E[e].size(); p++) { //for each point on that edge
  392 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); //set the vertex position based on the current point
  393 + glTexCoord1f(E[e].r(p)); //set the texture coordinate based on the specified magnitude index
  394 + }
  395 + glEnd();
  396 + }
  397 + glEndList();
  398 + }
  399 + glCallList(dlist + 4);
  400 + }
  401 +
  402 + // highlight the difference part
  403 + void glDifferenceCylinder(int I, std::vector<unsigned> map, std::vector<T> colormap, float sigma, float radius) {
  404 +
  405 + if (radius != sigma) // if render radius was changed by user, create a new display list
  406 + glDeleteLists(dlist + 6, 1);
  407 + if (!glIsList(dlist + 6)) { // if dlist isn't a display list, create it
  408 + glNewList(dlist + 6, GL_COMPILE); // start a new display list
  409 + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network
  410 + if (map[e] != unsigned(-1)) {
  411 + glEnable(GL_BLEND); //enable color blend
  412 + glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); //set blend function
  413 + glDisable(GL_DEPTH_TEST); //should disable depth
  414 +
  415 + if (I == 0) { // if it is to render GT
  416 + glColor4f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2], 0.1);
  417 + }
  418 + else { // if it is to render T
  419 + glColor4f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2], 0.1);
  420 + }
  421 +
  422 + for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
  423 + stim::circle<T> C1 = E[e].circ(p - 1);
  424 + stim::circle<T> C2 = E[e].circ(p);
  425 + C1.set_R(2 * radius); // re-scale the circle to the same
  426 + C2.set_R(2 * radius);
  427 + std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
  428 + std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
  429 + renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]);
  430 + }
  431 + glDisable(GL_BLEND);
  432 + glEnable(GL_DEPTH_TEST);
  433 + }
  434 + else {
  435 + glColor3f(1.0, 1.0, 1.0); // white color for the un-mapping edges
  436 + for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge
  437 + stim::circle<T> C1 = E[e].circ(p - 1);
  438 + stim::circle<T> C2 = E[e].circ(p);
  439 + C1.set_R(2 * radius); // scale the circle to the same
  440 + C2.set_R(2 * radius);
  441 + std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
  442 + std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
  443 + renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]);
  444 + }
  445 + }
  446 + }
  447 + for (unsigned n = 0; n < V.size(); n++) {
  448 +
  449 + size_t num_edge = V[n].e[0].size() + V[n].e[1].size();
  450 + if (num_edge > 1) { // if it is the joint vertex
  451 + glColor4f(0.3, 0.3, 0.3, 0.1); // gray
  452 + renderBall(V[n][0], V[n][1], V[n][2], 3 * radius, 20);
  453 + }
  454 + else { // if it is the terminal vertex
  455 + glColor4f(0.6, 0.6, 0.6, 0.1); // more white gray
  456 + renderBall(V[n][0], V[n][1], V[n][2], 3 * radius, 20);
  457 + }
  458 + }
  459 + glEndList();
  460 + }
  461 + glCallList(dlist + 6);
  462 + }
  463 +
  464 + //void glRandColorCenterlineGT(GLuint &dlist1, std::vector<unsigned> map, std::vector<T> colormap){
  465 + // if(!glIsList(dlist1)){
  466 + // dlist1 = glGenLists(1);
  467 + // glNewList(dlist1, GL_COMPILE);
  468 + // for(unsigned e = 0; e < E.size(); e++){
  469 + // if(map[e] != unsigned(-1)){
  470 + // glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]);
  471 + // glBegin(GL_LINE_STRIP);
  472 + // for(unsigned p = 0; p < E[e].size(); p++){
  473 + // glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  474 + // }
  475 + // glEnd();
  476 + // for (unsigned p = 0; p < E[e].size() - 1; p++) {
  477 + // renderCylinder(E[e][p][0], E[e][p][1], E[e][p][2], E[e][p + 1][0], E[e][p + 1][1], E[e][p + 1][2], 10, 20);
  478 + // }
  479 + // }
  480 + // else{
  481 + // glColor3f(1.0, 1.0, 1.0);
  482 + // glBegin(GL_LINE_STRIP);
  483 + // for(unsigned p = 0; p < E[e].size(); p++){
  484 + // glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  485 + // }
  486 + // glEnd();
  487 + // }
  488 + // }
  489 + // for (unsigned v = 0; v < V.size(); v++) {
  490 + // size_t num_edge = V[v].e[0].size() + V[v].e[1].size();
  491 + // if (num_edge > 1) {
  492 + // glColor3f(0.3, 0.3, 0.3); // gray color for vertex
  493 + // renderBall(V[v][0], V[v][1], V[v][2], 20, 20);
  494 + // }
  495 + // }
  496 + // glEndList();
  497 + // }
  498 + // glCallList(dlist1);
  499 + //}
  500 +
  501 + //void glRandColorCenterlineT(GLuint &dlist2, std::vector<unsigned> map, std::vector<T> colormap){
  502 + // if(!glIsList(dlist2)){
  503 + // dlist2 = glGenLists(1);
  504 + // glNewList(dlist2, GL_COMPILE);
  505 + // for(unsigned e = 0; e < E.size(); e++){
  506 + // if(map[e] != unsigned(-1)){
  507 + // glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]);
  508 + // glBegin(GL_LINE_STRIP);
  509 + // for(unsigned p = 0; p < E[e].size(); p++){
  510 + // glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  511 + // }
  512 + // glEnd();
  513 + // for (unsigned p = 0; p < E[e].size() - 1; p++) {
  514 + // renderCylinder(E[e][p][0], E[e][p][1], E[e][p][2], E[e][p + 1][0], E[e][p + 1][1], E[e][p + 1][2], 10, 20);
  515 + // }
  516 + // }
  517 + // else{
  518 + // glColor3f(1.0, 1.0, 1.0);
  519 + // glBegin(GL_LINE_STRIP);
  520 + // for(unsigned p = 0; p < E[e].size(); p++){
  521 + // glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  522 + // }
  523 + // glEnd();
  524 + // }
  525 + // }
  526 + // for (unsigned v = 0; v < V.size(); v++) {
  527 + // size_t num_edge = V[v].e[0].size() + V[v].e[1].size();
  528 + // if (num_edge > 1) {
  529 + // glColor3f(0.3, 0.3, 0.3); // gray color for vertex
  530 + // renderBall(V[v][0], V[v][1], V[v][2], 20, 20);
  531 + // }
  532 + // }
  533 + // glEndList();
  534 + // }
  535 + // glCallList(dlist2);
  536 + //}
  537 +
  538 +
  539 + //void renderCylinder(T x1, T y1, T z1, T x2, T y2, T z2, T radius, int subdivisions) {
  540 + // T dx = x2 - x1;
  541 + // T dy = y2 - y1;
  542 + // T dz = z2 - z1;
  543 + // /// handle the degenerate case with an approximation
  544 + // if (dz == 0)
  545 + // dz = .00000001;
  546 + // T d = sqrt(dx*dx + dy*dy + dz*dz);
  547 + // T ax = 57.2957795*acos(dz / d); // 180°/pi
  548 + // if (dz < 0.0)
  549 + // ax = -ax;
  550 + // T rx = -dy*dz;
  551 + // T ry = dx*dz;
  552 +
  553 + // glPushMatrix();
  554 + // glTranslatef(x1, y1, z1);
  555 + // glRotatef(ax, rx, ry, 0.0);
  556 +
  557 + // glutSolidCylinder(radius, d, subdivisions, 1);
  558 + // glPopMatrix();
  559 + //}
  560 +
  561 +
  562 + /// render the network centerline from swc file as a series of strips in different colors based on the neuronal type
  563 + /// glCenterline0_swc is for only one input
  564 + /*void glCenterline0_swc() {
  565 + if (!glIsList(dlist)) { // if dlist isn't a display list, create it
  566 + dlist = glGenLists(1); // generate a display list
  567 + glNewList(dlist, GL_COMPILE); // start a new display list
  568 + for (unsigned e = 0; e < E.size(); e++) {
  569 + int type = NT[e]; // get the neuronal type
  570 + switch (type) {
  571 + case 0:
  572 + glColor3f(1.0f, 1.0f, 1.0f); // white for undefined
  573 + glBegin(GL_LINE_STRIP);
  574 + for (unsigned p = 0; p < E[e].size(); p++) {
  575 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  576 + }
  577 + glEnd();
  578 + break;
  579 + case 1:
  580 + glColor3f(1.0f, 0.0f, 0.0f); // red for soma
  581 + glBegin(GL_LINE_STRIP);
  582 + for (unsigned p = 0; p < E[e].size(); p++) {
  583 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  584 + }
  585 + glEnd();
  586 + break;
  587 + case 2:
  588 + glColor3f(1.0f, 0.5f, 0.0f); // orange for axon
  589 + glBegin(GL_LINE_STRIP);
  590 + for (unsigned p = 0; p < E[e].size(); p++) {
  591 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  592 + }
  593 + glEnd();
  594 + break;
  595 + case 3:
  596 + glColor3f(1.0f, 1.0f, 0.0f); // yellow for undefined
  597 + glBegin(GL_LINE_STRIP);
  598 + for (unsigned p = 0; p < E[e].size(); p++) {
  599 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  600 + }
  601 + glEnd();
  602 + break;
  603 + case 4:
  604 + glColor3f(0.0f, 1.0f, 0.0f); // green for undefined
  605 + glBegin(GL_LINE_STRIP);
  606 + for (unsigned p = 0; p < E[e].size(); p++) {
  607 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  608 + }
  609 + glEnd();
  610 + break;
  611 + case 5:
  612 + glColor3f(0.0f, 1.0f, 1.0f); // verdant for undefined
  613 + glBegin(GL_LINE_STRIP);
  614 + for (unsigned p = 0; p < E[e].size(); p++) {
  615 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  616 + }
  617 + glEnd();
  618 + break;
  619 + case 6:
  620 + glColor3f(0.0f, 0.0f, 1.0f); // blue for undefined
  621 + glBegin(GL_LINE_STRIP);
  622 + for (unsigned p = 0; p < E[e].size(); p++) {
  623 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  624 + }
  625 + glEnd();
  626 + break;
  627 + case 7:
  628 + glColor3f(0.5f, 0.0f, 1.0f); // purple for undefined
  629 + glBegin(GL_LINE_STRIP);
  630 + for (unsigned p = 0; p < E[e].size(); p++) {
  631 + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]);
  632 + }
  633 + glEnd();
  634 + break;
  635 + }
  636 + }
  637 + glEndList(); //end the display list
  638 + }
  639 + glCallList(dlist); // render the display list
  640 + }*/
  641 +
  642 + ///render the network cylinder as a series of tubes
  643 + ///colors are based on metric values
  644 + //void glCylinder(float sigma) {
  645 + // if (!glIsList(dlist)) { //if dlist isn't a display list, create it
  646 + // dlist = glGenLists(1); //generate a display list
  647 + // glNewList(dlist, GL_COMPILE); //start a new display list
  648 + // for (unsigned e = 0; e < E.size(); e++) { //for each edge in the network
  649 + // for (unsigned p = 1; p < E[e].size() - 1; p++) { // for each point on that edge
  650 + // stim::circle<T> C1 = E[e].circ(p - 1);
  651 + // stim::circle<T> C2 = E[e].circ(p);
  652 + // C1.set_R(2.5*sigma); // scale the circle to the same
  653 + // C2.set_R(2.5*sigma);
  654 + // std::vector< stim::vec3<T> >Cp1 = C1.glpoints(20);
  655 + // std::vector< stim::vec3<T> >Cp2 = C2.glpoints(20);
  656 + // glBegin(GL_QUAD_STRIP);
  657 + // for (unsigned i = 0; i < Cp1.size(); i++) { // for every point on the circle(+1 means closing the circle)
  658 + // glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]);
  659 + // glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]);
  660 + // glTexCoord1f(E[e].r(p));
  661 + // }
  662 + // glEnd();
  663 + // } //set the texture coordinate based on the specified magnitude index
  664 + // }
  665 + // for (unsigned n = 0; n < V.size(); n++) {
  666 + // size_t num = V[n].e[0].size(); //store the number of outgoing edge of that vertex
  667 + // if (num != 0) { //if it has outgoing edge
  668 + // unsigned idx = V[n].e[0][0]; //find the index of first outgoing edge of that vertex
  669 + // glTexCoord1f(E[idx].r(0)); //bind the texture as metric of first point on that edge
  670 + // }
  671 + // else {
  672 + // unsigned idx = V[n].e[1][0]; //find the index of first incoming edge of that vertex
  673 + // glTexCoord1f(E[idx].r(E[idx].size() - 1)); //bind the texture as metric of last point on that edge
  674 + // }
  675 + // renderBall(V[n][0], V[n][1], V[n][2], 2.5*sigma, 20);
  676 + // }
  677 + // glEndList(); //end the display list
  678 + // }
  679 + // glCallList(dlist); //render the display list
  680 + //}
  681 +
  682 +}; //end stim::gl_network class
  683 +}; //end stim namespace
  684 +
  685 +
  686 +
  687 +#endif
0 \ No newline at end of file 688 \ No newline at end of file
@@ -7,8 +7,8 @@ @@ -7,8 +7,8 @@
7 #include <stim/parser/arguments.h> 7 #include <stim/parser/arguments.h>
8 #include <stim/visualization/camera.h> 8 #include <stim/visualization/camera.h>
9 #include <stim/gl/gl_texture.h> 9 #include <stim/gl/gl_texture.h>
10 -#include <stim/visualization/gl_network.h>  
11 -#include <stim/biomodels/network.h> 10 +#include "gl_network.h"
  11 +#include "network.h"
12 #include <stim/visualization/gl_aaboundingbox.h> 12 #include <stim/visualization/gl_aaboundingbox.h>
13 13
14 // OpenGL includes 14 // OpenGL includes
@@ -929,7 +929,7 @@ int main(int argc, char* argv[]) @@ -929,7 +929,7 @@ int main(int argc, char* argv[])
929 args.add("help", "prints this help"); 929 args.add("help", "prints this help");
930 args.add("sigma", "force a sigma value to specify the tolerance of the network comparison", "3"); 930 args.add("sigma", "force a sigma value to specify the tolerance of the network comparison", "3");
931 args.add("gui", "display the network or network comparison using OpenGL"); 931 args.add("gui", "display the network or network comparison using OpenGL");
932 - args.add("device", "choose specific device to run", "0"); 932 + args.add("device", "choose specific device to run (use -1 for CPU only)", "0");
933 args.add("features", "save features to a CSV file, specify file name"); 933 args.add("features", "save features to a CSV file, specify file name");
934 args.add("mapping", "mapping input according to similarity"); 934 args.add("mapping", "mapping input according to similarity");
935 args.add("stack", "load the image stacks"); 935 args.add("stack", "load the image stacks");
network.h 0 → 100644
  1 +#ifndef STIM_NETWORK_H
  2 +#define STIM_NETWORK_H
  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>
  11 +#include <stim/math/vec3.h>
  12 +#include <stim/visualization/obj.h>
  13 +#include <stim/visualization/swc.h>
  14 +#include "cylinder.h"
  15 +#include <stim/cuda/cudatools/timer.h>
  16 +#include <stim/cuda/cudatools/callable.h>
  17 +#include <stim/structures/kdtree.cuh>
  18 +//********************help function********************
  19 +// gaussian_function
  20 +CUDA_CALLABLE float gaussianFunction(float x, float std = 25) { return exp(-x / (2 * std*std)); } // std default sigma value is 25
  21 +
  22 +// compute metric in parallel
  23 +#ifdef __CUDACC__
  24 +template <typename T>
  25 +__global__ void find_metric_parallel(T* M, size_t n, T* D, float sigma){
  26 + size_t x = blockDim.x * blockIdx.x + threadIdx.x;
  27 + if(x >= n) return;
  28 + M[x] = 1.0f - gaussianFunction(D[x], sigma);
  29 +}
  30 +
  31 +//find the corresponding edge index from array index
  32 +__global__ void find_edge_index_parallel(size_t* I, size_t n, unsigned* R, size_t* E, size_t ne){
  33 + size_t x = blockDim.x * blockIdx.x + threadIdx.x;
  34 + if(x >= n) return;
  35 + unsigned i = 0;
  36 + size_t N = 0;
  37 + for(unsigned e = 0; e < ne; e++){
  38 + N += E[e];
  39 + if(I[x] < N){
  40 + R[x] = i;
  41 + break;
  42 + }
  43 + i++;
  44 + }
  45 +}
  46 +#endif
  47 +
  48 +//hard-coded factor
  49 +int threshold_fac;
  50 +
  51 +namespace stim{
  52 +/** This is the a class that interfaces with gl_spider in order to store the currently
  53 + * segmented network. The following data is stored and can be extracted:
  54 + * 1)Network geometry and centerline.
  55 + * 2)Network connectivity (a graph of nodes and edges), reconstructed using kdtree.
  56 +*/
  57 +
  58 +template<typename T>
  59 +class network{
  60 +
  61 + ///Each edge is a fiber with two nodes.
  62 + ///Each node is an in index to the endpoint of the fiber in the nodes array.
  63 + class edge : public cylinder<T>
  64 + {
  65 + public:
  66 +
  67 + unsigned int v[2]; //unique id's designating the starting and ending
  68 + // default constructor
  69 + edge() : cylinder<T>() {
  70 + v[1] = (unsigned)(-1); v[0] = (unsigned)(-1);
  71 + }
  72 + /// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
  73 +/*
  74 + ///@param v0, the starting index.
  75 + ///@param v1, the ending index.
  76 + ///@param sz, the number of point in the fiber.
  77 + edge(unsigned int v0, unsigned int v1, unsigned int sz) : cylinder<T>(
  78 + {
  79 +
  80 + }
  81 +*/
  82 + edge(std::vector<stim::vec3<T> > p, std::vector<T> s)
  83 + : cylinder<T>(p,s)
  84 + {
  85 + }
  86 + ///@param p is an array of positions in space
  87 + edge(stim::centerline<T> p) : cylinder<T>(p){}
  88 +
  89 + /// Copy constructor creates an edge from a cylinder
  90 + edge(stim::cylinder<T> f) : cylinder<T>(f) {}
  91 +
  92 + /// Resamples an edge by calling the fiber resampling function
  93 + edge resample(T spacing){
  94 + edge e(cylinder<T>::resample(spacing)); //call the fiber->edge constructor
  95 + e.v[0] = v[0]; //copy the vertex data
  96 + e.v[1] = v[1];
  97 +
  98 + return e; //return the new edge
  99 + }
  100 +
  101 + /// Output the edge information as a string
  102 + std::string str(){
  103 + std::stringstream ss;
  104 + ss<<"("<<cylinder<T>::size()<<")\tl = "<<this->length()<<"\t"<<v[0]<<"----"<<v[1];
  105 + return ss.str();
  106 + }
  107 +
  108 + std::vector<edge> split(unsigned int idx){
  109 +
  110 + std::vector< stim::cylinder<T> > C;
  111 + C.resize(2);
  112 + C = (*this).cylinder<T>::split(idx);
  113 + std::vector<edge> E(C.size());
  114 +
  115 + for(unsigned e = 0; e < E.size(); e++){
  116 + E[e] = C[e];
  117 + }
  118 + return E;
  119 + }
  120 +
  121 + /// operator for writing the edge information into a binary .nwt file.
  122 + friend std::ofstream& operator<<(std::ofstream& out, const edge& e)
  123 + {
  124 + out.write(reinterpret_cast<const char*>(&e.v[0]), sizeof(unsigned int)); ///write the starting point.
  125 + out.write(reinterpret_cast<const char*>(&e.v[1]), sizeof(unsigned int)); ///write the ending point.
  126 + unsigned int sz = e.size(); ///write the number of point in the edge.
  127 + out.write(reinterpret_cast<const char*>(&sz), sizeof(unsigned int));
  128 + for(int i = 0; i < sz; i++) ///write each point
  129 + {
  130 + stim::vec3<T> point = e[i];
  131 + out.write(reinterpret_cast<const char*>(&point[0]), 3*sizeof(T));
  132 + // for(int j = 0; j < nmags(); j++) //future code for multiple mags
  133 + // {
  134 + out.write(reinterpret_cast<const char*>(&e.R[i]), sizeof(T)); ///write the radius
  135 + //std::cout << point.str() << " " << e.R[i] << std::endl;
  136 + // }
  137 + }
  138 + return out; //return stream
  139 + }
  140 +
  141 + /// operator for reading an edge from a binary .nwt file.
  142 + friend std::ifstream& operator>>(std::ifstream& in, edge& e)
  143 + {
  144 + unsigned int v0, v1, sz;
  145 + in.read(reinterpret_cast<char*>(&v0), sizeof(unsigned int)); //read the staring point.
  146 + in.read(reinterpret_cast<char*>(&v1), sizeof(unsigned int)); //read the ending point
  147 + in.read(reinterpret_cast<char*>(&sz), sizeof(unsigned int)); //read the number of points in the edge
  148 +// stim::centerline<T> temp = stim::centerline<T>(sz); //allocate the new edge
  149 +// e = edge(temp);
  150 + std::vector<stim::vec3<T> > p(sz);
  151 + std::vector<T> r(sz);
  152 + for(int i = 0; i < sz; i++) //set the points and radii to the newly read values
  153 + {
  154 + stim::vec3<T> point;
  155 + in.read(reinterpret_cast<char*>(&point[0]), 3*sizeof(T));
  156 + p[i] = point;
  157 + T mag;
  158 + // for(int j = 0; j < nmags(); j++) ///future code for mags
  159 + // {
  160 + in.read(reinterpret_cast<char*>(&mag), sizeof(T));
  161 + r[i] = mag;
  162 + //std::cout << point.str() << " " << mag << std::endl;
  163 + // }
  164 + }
  165 + e = edge(p,r);
  166 + e.v[0] = v0; e.v[1] = v1;
  167 + return in;
  168 + }
  169 + };
  170 +
  171 + ///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.
  172 + class vertex : public stim::vec3<T>
  173 + {
  174 + public:
  175 + //std::vector<unsigned int> edges; //indices of edges connected to this node.
  176 + std::vector<unsigned int> e[2]; //indices of edges going out (e[0]) and coming in (e[1])
  177 + //stim::vec3<T> p; //position of this node in physical space.
  178 + //default constructor
  179 + vertex() : stim::vec3<T>()
  180 + {
  181 + }
  182 + //constructor takes a stim::vec
  183 + vertex(stim::vec3<T> p) : stim::vec3<T>(p){}
  184 +
  185 + /// Output the vertex information as a string
  186 + std::string
  187 + str(){
  188 + std::stringstream ss;
  189 + ss<<"\t(x, y, z) = "<<stim::vec3<T>::str();
  190 +
  191 + if(e[0].size() > 0){
  192 + ss<<"\t> ";
  193 + for(unsigned int o = 0; o < e[0].size(); o++)
  194 + ss<<e[0][o]<<" ";
  195 + }
  196 + if(e[1].size() > 0){
  197 + ss<<"\t< ";
  198 + for(unsigned int i = 0; i < e[1].size(); i++)
  199 + ss<<e[1][i]<<" ";
  200 + }
  201 +
  202 + return ss.str();
  203 + }
  204 + ///operator for writing the vector into the stream;
  205 + friend std::ofstream& operator<<(std::ofstream& out, const vertex& v)
  206 + {
  207 + unsigned int s0, s1;
  208 + s0 = v.e[0].size();
  209 + s1 = v.e[1].size();
  210 + out.write(reinterpret_cast<const char*>(&v.ptr[0]), 3*sizeof(T)); ///write physical vertex location
  211 + out.write(reinterpret_cast<const char*>(&s0), sizeof(unsigned int)); ///write the number of "outgoing edges"
  212 + out.write(reinterpret_cast<const char*>(&s1), sizeof(unsigned int)); ///write the number of "incoming edges"
  213 + if (s0 != 0)
  214 + out.write(reinterpret_cast<const char*>(&v.e[0][0]), sizeof(unsigned int)*v.e[0].size()); ///write the "outgoing edges"
  215 + if (s1 != 0)
  216 + out.write(reinterpret_cast<const char*>(&v.e[1][0]), sizeof(unsigned int)*v.e[1].size()); ///write the "incoming edges"
  217 + return out;
  218 + }
  219 +
  220 + ///operator for reading the vector out of the stream;
  221 + friend std::ifstream& operator>>(std::ifstream& in, vertex& v)
  222 + {
  223 + in.read(reinterpret_cast<char*>(&v[0]), 3*sizeof(T)); ///read the physical position
  224 + unsigned int s[2];
  225 + in.read(reinterpret_cast<char*>(&s[0]), 2*sizeof(unsigned int)); ///read the sizes of incoming and outgoing edge arrays
  226 +
  227 + std::vector<unsigned int> one(s[0]);
  228 + std::vector<unsigned int> two(s[1]);
  229 + v.e[0] = one;
  230 + v.e[1] = two;
  231 + if (one.size() != 0)
  232 + in.read(reinterpret_cast<char*>(&v.e[0][0]), s[0] * sizeof(unsigned int)); ///read the arrays of "outgoing edges"
  233 + if (two.size() != 0)
  234 + in.read(reinterpret_cast<char*>(&v.e[1][0]), s[1] * sizeof(unsigned int)); ///read the arrays of "incoming edges"
  235 + return in;
  236 + }
  237 +
  238 + };
  239 +
  240 +protected:
  241 +
  242 + std::vector<edge> E; //list of edges
  243 + std::vector<vertex> V; //list of vertices.
  244 +
  245 +public:
  246 +
  247 + ///default constructor
  248 + network()
  249 + {
  250 +
  251 + }
  252 +
  253 + ///constructor with a file to load.
  254 + network(std::string fileLocation)
  255 + {
  256 + load_obj(fileLocation);
  257 + }
  258 +
  259 + ///Returns the number of edges in the network.
  260 + unsigned int edges(){
  261 + return E.size();
  262 + }
  263 +
  264 + ///Returns the number of nodes in the network.
  265 + unsigned int vertices(){
  266 + return V.size();
  267 + }
  268 +
  269 + ///Returns the radius at specific point in the edge
  270 + T get_r(unsigned e, unsigned i) {
  271 + return E[e].r(i);
  272 + }
  273 +
  274 + ///Returns the average radius of specific edge
  275 + T get_average_r(unsigned e) {
  276 + T result = 0.0;
  277 + unsigned n = E[e].size();
  278 + for (unsigned p = 0; p < n; p++)
  279 + result += E[e].r(p);
  280 +
  281 + return (T)result / n;
  282 + }
  283 +
  284 + ///Returns the length of current edge
  285 + T get_l(unsigned e) {
  286 + return E[e].length();
  287 + }
  288 +
  289 + ///Returns the start vertex of current edge
  290 + size_t get_start_vertex(unsigned e) {
  291 + return E[e].v[0];
  292 + }
  293 +
  294 + ///Returns the end vertex of current edge
  295 + size_t get_end_vertex(unsigned e) {
  296 + return E[e].v[1];
  297 + }
  298 +
  299 + ///Returns one vertex
  300 + stim::vec3<T> get_vertex(unsigned i) {
  301 + return V[i];
  302 + }
  303 +
  304 + ///Returns the boundary vertices' indices
  305 + std::vector<unsigned> get_boundary_vertex() {
  306 + std::vector<unsigned> result;
  307 +
  308 + for (unsigned v = 0; v < V.size(); v++) {
  309 + if (V[v].e[0].size() + V[v].e[1].size() == 1) { // boundary vertex
  310 + result.push_back(v);
  311 + }
  312 + }
  313 +
  314 + return result;
  315 + }
  316 +
  317 + ///Set radius
  318 + void set_r(unsigned e, std::vector<T> radius) {
  319 + E[e].cylinder<T>::copy_r(radius);
  320 + }
  321 +
  322 + void set_r(unsigned e, T radius) {
  323 + for (size_t i = 0; i < E[e].size(); i++)
  324 + E[e].cylinder<T>::set_r(i, radius);
  325 + }
  326 + //scale the network by some constant value
  327 + // I don't think these work??????
  328 + /*std::vector<vertex> operator*(T s){
  329 + for (unsigned i=0; i< vertices; i ++ ){
  330 + V[i] = V[i] * s;
  331 + }
  332 + return V;
  333 + }
  334 +
  335 + std::vector<vertex> operator*(vec<T> s){
  336 + for (unsigned i=0; i< vertices; i ++ ){
  337 + for (unsigned dim = 0 ; dim< 3; dim ++){
  338 + V[i][dim] = V[i][dim] * s[dim];
  339 + }
  340 + }
  341 + return V;
  342 + }*/
  343 +
  344 + // Returns an average of branching index in the network
  345 + double BranchingIndex(){
  346 + double B=0;
  347 + for(unsigned v=0; v < V.size(); v ++){
  348 + B += ((V[v].e[0].size()) + (V[v].e[1].size()));
  349 + }
  350 + B = B / V.size();
  351 + return B;
  352 +
  353 + }
  354 +
  355 + // Returns number of branch points in thenetwork
  356 + unsigned int BranchP(){
  357 + unsigned int B=0;
  358 + unsigned int c;
  359 + for(unsigned v=0; v < V.size(); v ++){
  360 + c = ((V[v].e[0].size()) + (V[v].e[1].size()));
  361 + if (c > 2){
  362 + B += 1;}
  363 + }
  364 + return B;
  365 +
  366 + }
  367 +
  368 + // Returns number of end points (tips) in thenetwork
  369 + unsigned int EndP(){
  370 + unsigned int B=0;
  371 + unsigned int c;
  372 + for(unsigned v=0; v < V.size(); v ++){
  373 + c = ((V[v].e[0].size()) + (V[v].e[1].size()));
  374 + if (c == 1){
  375 + B += 1;}
  376 + }
  377 + return B;
  378 +
  379 + }
  380 +
  381 + //// Returns a dictionary with the key as the vertex
  382 + //std::map<std::vector<vertex>,unsigned int> DegreeDict(){
  383 + // std::map<std::vector<vertex>,unsigned int> dd;
  384 + // unsigned int c = 0;
  385 + // for(unsigned v=0; v < V.size(); v ++){
  386 + // c = ((V[v].e[0].size()) + (V[v].e[1].size()));
  387 + // dd[V[v]] = c;
  388 + // }
  389 + // return dd;
  390 + //}
  391 +
  392 + //// Return number of branching stems
  393 + //unsigned int Stems(){
  394 + // unsigned int s = 0;
  395 + // std::map<std::vector<vertex>,unsigned int> dd;
  396 + // dd = DegreeDict();
  397 + // //for(unsigned v=0; v < V.size(); v ++){
  398 + // // V[v].e[0].
  399 + // return s;
  400 + //}
  401 +
  402 + //Calculate Metrics---------------------------------------------------
  403 + // Returns an average of fiber/edge lengths in the network
  404 + double Lengths(){
  405 + stim::vec<T> L;
  406 + double sumLength = 0;
  407 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  408 + L.push_back(E[e].length()); //append the edge length
  409 + sumLength = sumLength + E[e].length();
  410 + }
  411 + double avg = sumLength / E.size();
  412 + return avg;
  413 + }
  414 +
  415 +
  416 + // Returns an average of tortuosities in the network
  417 + double Tortuosities(){
  418 + stim::vec<T> t;
  419 + stim::vec<T> id1, id2; // starting and ending vertices of the edge
  420 + double distance;double tortuosity;double sumTortuosity = 0;
  421 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  422 + id1 = E[e][0]; //get the edge starting point
  423 + id2 = E[e][E[e].size() - 1]; //get the edge ending point
  424 + distance = (id1 - id2).len(); //displacement between the starting and ending points
  425 + if(distance > 0){
  426 + tortuosity = E[e].length()/ distance ; // tortuoisty = edge length / edge displacement
  427 + }
  428 + else{
  429 + tortuosity = 0;}
  430 + t.push_back(tortuosity);
  431 + sumTortuosity += tortuosity;
  432 + }
  433 + double avg = sumTortuosity / E.size();
  434 + return avg;
  435 + }
  436 +
  437 + // Returns average contraction of the network
  438 + double Contractions(){
  439 + stim::vec<T> t;
  440 + stim::vec<T> id1, id2; // starting and ending vertices of the edge
  441 + double distance;double contraction;double sumContraction = 0;
  442 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  443 + id1 = E[e][0]; //get the edge starting point
  444 + id2 = E[e][E[e].size() - 1]; //get the edge ending point
  445 + distance = (id1 - id2).len(); //displacement between the starting and ending points
  446 + contraction = distance / E[e].length(); // tortuoisty = edge length / edge displacement
  447 + t.push_back(contraction);
  448 + sumContraction += contraction;
  449 + }
  450 + double avg = sumContraction / E.size();
  451 + return avg;
  452 + }
  453 +
  454 + // returns average fractal dimension of the branches of the network
  455 + double FractalDimensions(){
  456 + stim::vec<T> t;
  457 + stim::vec<T> id1, id2; // starting and ending vertices of the edge
  458 + double distance;double fract;double sumFractDim = 0;
  459 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  460 + id1 = E[e][0]; //get the edge starting point
  461 + id2 = E[e][E[e].size() - 1]; //get the edge ending point
  462 + distance = (id1 - id2).len(); //displacement between the starting and ending points
  463 + fract = std::log(distance) / std::log(E[e].length()); // tortuoisty = edge length / edge displacement
  464 + t.push_back(sumFractDim);
  465 + sumFractDim += fract;
  466 + }
  467 + double avg = sumFractDim / E.size();
  468 + return avg;
  469 + }
  470 +
  471 + //returns a cylinder represented a given fiber (based on edge index)
  472 + stim::cylinder<T> get_cylinder(unsigned e){
  473 + return E[e]; //return the specified edge (casting it to a fiber)
  474 + }
  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 + unsigned num_edge = (*this).E.size();
  485 +
  486 + for (unsigned i = 0; i < num_edge; i++) {
  487 + if (E[i].size() == 2) { // if current edge can't be subdivided
  488 + stim::centerline<T> line(2);
  489 + for (unsigned k = 0; k < 2; k++)
  490 + line[k] = E[i][k];
  491 + line.update();
  492 +
  493 + edge new_edge(line);
  494 +
  495 + vertex new_vertex = new_edge[0];
  496 + id = E[i].v[0];
  497 + auto position = std::find(ori_index.begin(), ori_index.end(), id);
  498 + if (position == ori_index.end()) { // new vertex
  499 + ori_index.push_back(id);
  500 + new_index.push_back(nV.size());
  501 +
  502 + new_vertex.e[0].push_back(nE.size());
  503 + new_edge.v[0] = nV.size();
  504 + nV.push_back(new_vertex); // push back vertex as a new vertex
  505 + }
  506 + else { // existing vertex
  507 + int k = std::distance(ori_index.begin(), position);
  508 + new_edge.v[0] = new_index[k];
  509 + nV[new_index[k]].e[0].push_back(nE.size());
  510 + }
  511 +
  512 + new_vertex = new_edge[1];
  513 + id = E[i].v[1];
  514 + position = std::find(ori_index.begin(), ori_index.end(), id);
  515 + if (position == ori_index.end()) { // new vertex
  516 + ori_index.push_back(id);
  517 + new_index.push_back(nV.size());
  518 +
  519 + new_vertex.e[1].push_back(nE.size());
  520 + new_edge.v[1] = nV.size();
  521 + nV.push_back(new_vertex); // push back vertex as a new vertex
  522 + }
  523 + else { // existing vertex
  524 + int k = std::distance(ori_index.begin(), position);
  525 + new_edge.v[1] = new_index[k];
  526 + nV[new_index[k]].e[1].push_back(nE.size());
  527 + }
  528 +
  529 + nE.push_back(new_edge);
  530 +
  531 + nE[nE.size() - 1].cylinder<T>::set_r(0, E[i].cylinder<T>::r(0));
  532 + nE[nE.size() - 1].cylinder<T>::set_r(1, E[i].cylinder<T>::r(1));
  533 + }
  534 + else { // subdivide current edge
  535 + for (unsigned j = 0; j < E[i].size() - 1; j++) {
  536 + stim::centerline<T> line(2);
  537 + for (unsigned k = 0; k < 2; k++)
  538 + line[k] = E[i][j + k];
  539 + line.update();
  540 +
  541 + edge new_edge(line);
  542 +
  543 + if (j == 0) { // edge contains original starting point
  544 + vertex new_vertex = new_edge[0];
  545 + id = E[i].v[0];
  546 + auto position = std::find(ori_index.begin(), ori_index.end(), id);
  547 + if (position == ori_index.end()) { // new vertex
  548 + ori_index.push_back(id);
  549 + new_index.push_back(nV.size());
  550 +
  551 + new_vertex.e[0].push_back(nE.size());
  552 + new_edge.v[0] = nV.size();
  553 + nV.push_back(new_vertex); // push back vertex as a new vertex
  554 + }
  555 + else { // existing vertex
  556 + int k = std::distance(ori_index.begin(), position);
  557 + new_edge.v[0] = new_index[k];
  558 + nV[new_index[k]].e[0].push_back(nE.size());
  559 + }
  560 +
  561 + new_vertex = new_edge[1];
  562 + new_vertex.e[1].push_back(nE.size());
  563 + new_edge.v[1] = nV.size();
  564 + nV.push_back(new_vertex); // push back internal point as a new vertex
  565 +
  566 + nE.push_back(new_edge);
  567 + }
  568 +
  569 + else if (j == E[i].size() - 2) { // edge contains original ending point
  570 +
  571 + vertex new_vertex = new_edge[1];
  572 + nV[nV.size() - 1].e[0].push_back(nE.size());
  573 + new_edge.v[0] = nV.size() - 1;
  574 +
  575 + id = E[i].v[1];
  576 + auto position = std::find(ori_index.begin(), ori_index.end(), id);
  577 + if (position == ori_index.end()) { // new vertex
  578 + ori_index.push_back(id);
  579 + new_index.push_back(nV.size());
  580 +
  581 + new_vertex.e[1].push_back(nE.size());
  582 + new_edge.v[1] = nV.size();
  583 + nV.push_back(new_vertex); // push back vertex as a new vertex
  584 + }
  585 + else { // existing vertex
  586 + int k = std::distance(ori_index.begin(), position);
  587 + new_edge.v[1] = new_index[k];
  588 + nV[new_index[k]].e[1].push_back(nE.size());
  589 + }
  590 +
  591 + nE.push_back(new_edge);
  592 + }
  593 +
  594 + else {
  595 + vertex new_vertex = new_edge[1];
  596 +
  597 + nV[nV.size() - 1].e[0].push_back(nE.size());
  598 + new_vertex.e[1].push_back(nE.size());
  599 + new_edge.v[0] = nV.size() - 1;
  600 + new_edge.v[1] = nV.size();
  601 + nV.push_back(new_vertex);
  602 +
  603 + nE.push_back(new_edge);
  604 + }
  605 +
  606 + // get radii
  607 + nE[nE.size() - 1].cylinder<T>::set_r(0, E[i].cylinder<T>::r(j));
  608 + nE[nE.size() - 1].cylinder<T>::set_r(1, E[i].cylinder<T>::r(j + 1));
  609 + }
  610 + }
  611 + }
  612 +
  613 + (*this).E = nE;
  614 + (*this).V = nV;
  615 + }
  616 +
  617 + //load a network from an OBJ file
  618 + void load_obj(std::string filename){
  619 +
  620 + stim::obj<T> O; //create an OBJ object
  621 + O.load(filename); //load the OBJ file as an object
  622 +
  623 + std::vector<unsigned> id2vert; //this list stores the OBJ vertex ID associated with each network vertex
  624 +
  625 + unsigned i[2]; //temporary, IDs associated with the first and last points in an OBJ line
  626 +
  627 + //for each line in the OBJ object
  628 + for(unsigned int l = 1; l <= O.numL(); l++){
  629 +
  630 + std::vector< stim::vec<T> > c; //allocate an array of points for the vessel centerline
  631 + O.getLine(l, c); //get the fiber centerline
  632 +
  633 + stim::centerline<T> c3(c.size());
  634 + for(size_t j = 0; j < c.size(); j++)
  635 + c3[j] = c[j];
  636 + c3.update();
  637 +
  638 + // edge new_edge = c3; ///This is dangerous.
  639 + edge new_edge(c3);
  640 +
  641 + //create an edge from the given centerline
  642 + unsigned int I = new_edge.size(); //calculate the number of points on the centerline
  643 +
  644 + //get the first and last vertex IDs for the line
  645 + std::vector< unsigned > id; //create an array to store the centerline point IDs
  646 + O.getLinei(l, id); //get the list of point IDs for the line
  647 + i[0] = id.front(); //get the OBJ ID for the first element of the line
  648 + i[1] = id.back(); //get the OBJ ID for the last element of the line
  649 +
  650 + std::vector<unsigned>::iterator it; //create an iterator for searching the id2vert array
  651 + unsigned it_idx; //create an integer for the id2vert entry index
  652 +
  653 + //find out if the nodes for this fiber have already been created
  654 + it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node
  655 + if(it == id2vert.end()){ //if i[0] hasn't already been used
  656 + vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position
  657 + bool flag = false;
  658 + unsigned j = 0;
  659 + for (; j < V.size(); j++) { // check whether current vertex is already exist
  660 + if (new_vertex == V[j]) {
  661 + flag = true;
  662 + break;
  663 + }
  664 + }
  665 + if (!flag) { // unique one
  666 + new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing
  667 + new_edge.v[0] = V.size(); //add the new edge to the edge
  668 + V.push_back(new_vertex); //add the new vertex to the vertex list
  669 + id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list
  670 + }
  671 + else {
  672 + V[j].e[0].push_back(E.size());
  673 + new_edge.v[0] = j;
  674 + }
  675 + }
  676 + else{ //if the vertex already exists
  677 + it_idx = std::distance(id2vert.begin(), it);
  678 + V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing
  679 + new_edge.v[0] = it_idx;
  680 + }
  681 +
  682 + it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID
  683 + if(it == id2vert.end()){ //if i[1] hasn't already been used
  684 + vertex new_vertex = new_edge[I-1]; //create a new vertex, assign it a position
  685 + bool flag = false;
  686 + unsigned j = 0;
  687 + for (; j < V.size(); j++) { // check whether current vertex is already exist
  688 + if (new_vertex == V[j]) {
  689 + flag = true;
  690 + break;
  691 + }
  692 + }
  693 + if (!flag) {
  694 + new_vertex.e[1].push_back(E.size()); //add the current edge as incoming
  695 + new_edge.v[1] = V.size(); //add the new vertex to the edge
  696 + V.push_back(new_vertex); //add the new vertex to the vertex list
  697 + id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list
  698 + }
  699 + else {
  700 + V[j].e[1].push_back(E.size());
  701 + new_edge.v[1] = j;
  702 + }
  703 + }
  704 + else{ //if the vertex already exists
  705 + it_idx = std::distance(id2vert.begin(), it);
  706 + V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming
  707 + new_edge.v[1] = it_idx;
  708 + }
  709 +
  710 + E.push_back(new_edge); //push the edge to the list
  711 +
  712 + }
  713 +
  714 + // copy the radii information from OBJ
  715 + /*if (O.numVT()) {
  716 + unsigned k = 0;
  717 + for (unsigned i = 0; i < E.size(); i++) {
  718 + for (unsigned j = 0; j < E[i].size(); j++) {
  719 + E[i].cylinder<T>::set_r(j, O.getVT(k)[0] / 2);
  720 + k++;
  721 + }
  722 + }
  723 + }*/
  724 + // OBJ class assumes that in L the two values are equal
  725 + if (O.numVT()) {
  726 + std::vector< unsigned > id; //create an array to store the centerline point IDs
  727 + for (unsigned i = 0; i < O.numL(); i++) {
  728 + id.clear();
  729 + O.getLinei(i + 1, id); //get the list of point IDs for the line
  730 + for (unsigned j = 0; j < id.size(); j++)
  731 + E[i].cylinder<T>::set_r(j, O.getVT(id[j] - 1)[0] / 2);
  732 + }
  733 + }
  734 + }
  735 +
  736 + ///loads a .nwt file. Reads the header and loads the data into the network according to the header.
  737 + void
  738 + loadNwt(std::string filename)
  739 + {
  740 + int dims[2]; ///number of vertex, number of edges
  741 + readHeader(filename, &dims[0]); //read header
  742 + std::ifstream file;
  743 + file.open(filename.c_str(), std::ios::in | std::ios::binary); ///skip header information.
  744 + file.seekg(14+58+4+4, file.beg);
  745 + vertex v;
  746 + for(int i = 0; i < dims[0]; i++) ///for every vertex, read vertex, add to network.
  747 + {
  748 + file >> v;
  749 + V.push_back(v);
  750 +// std::cout << i << " " << v.str() << std::endl;
  751 + }
  752 +
  753 + std::cout << std::endl;
  754 + for(int i = 0; i < dims[1]; i++) ///for every edge, read edge, add to network.
  755 + {
  756 + edge e;
  757 + file >> e;
  758 + E.push_back(e);
  759 + //std::cout << i << " " << E[i].str() << std::endl; // not necessary?
  760 + }
  761 + file.close();
  762 + }
  763 +
  764 + ///saves a .nwt file. Writes the header in raw text format, then saves the network as a binary file.
  765 + void
  766 + saveNwt(std::string filename)
  767 + {
  768 + writeHeader(filename);
  769 + std::ofstream file;
  770 + file.open(filename.c_str(), std::ios::out | std::ios::binary | std::ios::app); ///since we have written the header we are not appending.
  771 + for(int i = 0; i < V.size(); i++) ///look through the Vertices and write each one.
  772 + {
  773 +// std::cout << i << " " << V[i].str() << std::endl;
  774 + file << V[i];
  775 + }
  776 + for(int i = 0; i < E.size(); i++) ///loop through the Edges and write each one.
  777 + {
  778 + //std::cout << i << " " << E[i].str() << std::endl; // not necesarry?
  779 + file << E[i];
  780 + }
  781 + file.close();
  782 + }
  783 +
  784 +
  785 + ///Writes the header information to a .nwt file.
  786 + void
  787 + writeHeader(std::string filename)
  788 + {
  789 + std::string magicString = "nwtFileFormat "; ///identifier for the file.
  790 + std::string desc = "fileid(14B), desc(58B), #vertices(4B), #edges(4B): bindata";
  791 + int hNumVertices = V.size(); ///int byte header storing the number of vertices in the file
  792 + int hNumEdges = E.size(); ///int byte header storing the number of edges.
  793 + std::ofstream file;
  794 + file.open(filename.c_str(), std::ios::out | std::ios::binary);
  795 + std::cout << hNumVertices << " " << hNumEdges << std::endl;
  796 + file.write(reinterpret_cast<const char*>(&magicString.c_str()[0]), 14); //write the file id
  797 + file.write(reinterpret_cast<const char*>(&desc.c_str()[0]), 58); //write the description
  798 + file.write(reinterpret_cast<const char*>(&hNumVertices), sizeof(int)); //write #vert.
  799 + file.write(reinterpret_cast<const char*>(&hNumEdges), sizeof(int)); //write #edges
  800 +// file << magicString.c_str() << desc.c_str() << hNumVertices << hNumEdges;
  801 + file.close();
  802 +
  803 + }
  804 +
  805 + ///Reads the header information from a .nwt file.
  806 + void
  807 + readHeader(std::string filename, int *dims)
  808 + {
  809 + char magicString[14]; ///id
  810 + char desc[58]; ///description
  811 + int hNumVertices; ///#vert
  812 + int hNumEdges; ///#edges
  813 + std::ifstream file; ////create stream
  814 + file.open(filename.c_str(), std::ios::in | std::ios::binary);
  815 + file.read(reinterpret_cast<char*>(&magicString[0]), 14); ///read the file id.
  816 + file.read(reinterpret_cast<char*>(&desc[0]), 58); ///read the description
  817 + file.read(reinterpret_cast<char*>(&hNumVertices), sizeof(int)); ///read the number of vertices
  818 + file.read(reinterpret_cast<char*>(&hNumEdges), sizeof(int)); ///read the number of edges
  819 +// std::cout << magicString << desc << hNumVertices << " " << hNumEdges << std::endl;
  820 + file.close(); ///close the file.
  821 + dims[0] = hNumVertices; ///fill the returned reference.
  822 + dims[1] = hNumEdges;
  823 + }
  824 +
  825 + //load a network from an SWC file
  826 + void load_swc(std::string filename) {
  827 + stim::swc<T> S; // create swc variable
  828 + S.load(filename); // load the node information
  829 + S.create_tree(); // link those node according to their linking relationships as a tree
  830 + S.resample();
  831 +
  832 + //NT.push_back(S.node[0].type); // set the neuronal_type value to the first vertex in the network
  833 + std::vector<unsigned> id2vert; // this list stores the SWC vertex ID associated with each network vertex
  834 + unsigned i[2]; // temporary, IDs associated with the first and last points
  835 +
  836 + for (unsigned int l = 0; l < S.numE(); l++) { // for every edge
  837 + //NT.push_back(S.node[l].type);
  838 +
  839 + std::vector< stim::vec3<T> > c;
  840 + S.get_points(l, c);
  841 +
  842 + stim::centerline<T> c3(c.size()); // new fiber
  843 +
  844 + for (unsigned j = 0; j < c.size(); j++)
  845 + c3[j] = c[j]; // copy the points
  846 +
  847 + c3.update(); // update the L information
  848 +
  849 + stim::cylinder<T> C3(c3); // create a new cylinder in order to copy the origin radius information
  850 + // upadate the R information
  851 + std::vector<T> radius;
  852 + S.get_radius(l, radius);
  853 +
  854 + C3.copy_r(radius);
  855 +
  856 + edge new_edge(C3); // new edge
  857 +
  858 + //create an edge from the given centerline
  859 + unsigned int I = (unsigned)new_edge.size(); //calculate the number of points on the centerline
  860 +
  861 + //get the first and last vertex IDs for the line
  862 + i[0] = S.E[l].front();
  863 + i[1] = S.E[l].back();
  864 +
  865 + std::vector<unsigned>::iterator it; //create an iterator for searching the id2vert array
  866 + unsigned it_idx; //create an integer for the id2vert entry index
  867 +
  868 + //find out if the nodes for this fiber have already been created
  869 + it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node
  870 + if (it == id2vert.end()) { //if i[0] hasn't already been used
  871 + vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position
  872 + new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing
  873 + new_edge.v[0] = V.size(); //add the new edge to the edge
  874 + V.push_back(new_vertex); //add the new vertex to the vertex list
  875 + id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list
  876 + }
  877 + else { //if the vertex already exists
  878 + it_idx = std::distance(id2vert.begin(), it);
  879 + V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing
  880 + new_edge.v[0] = it_idx;
  881 + }
  882 +
  883 + it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID
  884 + if (it == id2vert.end()) { //if i[1] hasn't already been used
  885 + vertex new_vertex = new_edge[I - 1]; //create a new vertex, assign it a position
  886 + new_vertex.e[1].push_back(E.size()); //add the current edge as incoming
  887 + new_edge.v[1] = V.size(); //add the new vertex to the edge
  888 + V.push_back(new_vertex); //add the new vertex to the vertex list
  889 + id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list
  890 + }
  891 + else { //if the vertex already exists
  892 + it_idx = std::distance(id2vert.begin(), it);
  893 + V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming
  894 + new_edge.v[1] = it_idx;
  895 + }
  896 +
  897 + E.push_back(new_edge); //push the edge to the list
  898 + }
  899 + }
  900 +
  901 + /// Get adjacency matrix of the network
  902 + std::vector< typename std::vector<int> > get_adj_mat() {
  903 +
  904 + unsigned n = V.size(); // get the number of vertices in the networks
  905 +
  906 + std::vector< typename std::vector<int> > result(n, std::vector<int>(n, 0)); // initialize every entry in the matrix to be 0
  907 + result.resize(n); // resize rows
  908 + for (unsigned i = 0; i < n; i++)
  909 + result[i].resize(n); // resize columns
  910 +
  911 + for (unsigned i = 0; i < n; i++) { // for every vertex
  912 + unsigned num_out = V[i].e[0].size(); // number of outgoing edges of current vertex
  913 + if (num_out != 0) {
  914 + for (unsigned j = 0; j < num_out; j++) {
  915 + int edge_idx = V[i].e[0][j]; // get the jth out-going edge index of current vertex
  916 + int vertex_idx = E[edge_idx].v[1]; // get the ending vertex of specific out-going edge
  917 + result[i][vertex_idx] = 1; // can simply set to 1 if it is simple-graph
  918 + result[vertex_idx][i] = 1; // symmetric
  919 + }
  920 + }
  921 + }
  922 +
  923 + return result;
  924 + }
  925 +
  926 + /// Output the network as a string
  927 + std::string str(){
  928 +
  929 + std::stringstream ss;
  930 + ss<<"Nodes ("<<V.size()<<")--------"<<std::endl;
  931 + for(unsigned int v = 0; v < V.size(); v++){
  932 + ss<<"\t"<<v<<V[v].str()<<std::endl;
  933 + }
  934 +
  935 + ss<<"Edges ("<<E.size()<<")--------"<<std::endl;
  936 + for(unsigned e = 0; e < E.size(); e++){
  937 + ss<<"\t"<<e<<E[e].str()<<std::endl;
  938 + }
  939 +
  940 + return ss.str();
  941 + }
  942 +
  943 + /// This function resamples all fibers in a network given a desired minimum spacing
  944 + /// @param spacing is the minimum distance between two points on the network
  945 + stim::network<T> resample(T spacing){
  946 + stim::network<T> n; //create a new network that will be an exact copy, with resampled fibers
  947 + n.V = V; //copy all vertices
  948 + //n.NT = NT; //copy all the neuronal type information
  949 + n.E.resize(edges()); //allocate space for the edge list
  950 +
  951 + //copy all fibers, resampling them in the process
  952 + for(unsigned e = 0; e < edges(); e++){ //for each edge in the edge list
  953 + n.E[e] = E[e].resample(spacing); //resample the edge and copy it to the new network
  954 + }
  955 +
  956 + return n; //return the resampled network
  957 + }
  958 +
  959 + /// Calculate the total number of points on all edges.
  960 + unsigned total_points(){
  961 + unsigned n = 0;
  962 + for(unsigned e = 0; e < E.size(); e++)
  963 + n += E[e].size();
  964 + return n;
  965 + }
  966 +
  967 + //Copy the point cloud representing the centerline for the network into an array
  968 + void centerline_cloud(T* dst) {
  969 + size_t p; //stores the current edge point
  970 + size_t P; //stores the number of points in an edge
  971 + size_t i = 0; //index into the output array of points
  972 + for (size_t e = 0; e < E.size(); e++) { //for each edge in the network
  973 + P = E[e].size(); //get the number of points in this edge
  974 + for (p = 0; p < P; p++) {
  975 + dst[i * 3 + 0] = E[e][p][0];
  976 + dst[i * 3 + 1] = E[e][p][1];
  977 + dst[i * 3 + 2] = E[e][p][2];
  978 + i++;
  979 + }
  980 + }
  981 + }
  982 +
  983 + // convert vec3 to array
  984 + void stim2array(float *a, stim::vec3<T> b){
  985 + a[0] = b[0];
  986 + a[1] = b[1];
  987 + a[2] = b[2];
  988 + }
  989 +
  990 + // convert vec3 to array in bunch
  991 + void edge2array(T* a, edge b){
  992 + size_t n = b.size();
  993 + for(size_t i = 0; i < n; i++){
  994 + a[i * 3 + 0] = b[i][0];
  995 + a[i * 3 + 1] = b[i][1];
  996 + a[i * 3 + 2] = b[i][2];
  997 + }
  998 + }
  999 +
  1000 + // get list of metric
  1001 + std::vector<T> metric() {
  1002 + std::vector<T> result;
  1003 + T m;
  1004 + for (size_t e = 0; e < E.size(); e++) {
  1005 + for (size_t p = 0; p < E[e].size(); p++) {
  1006 + m = E[e].r(p);
  1007 + result.push_back(m);
  1008 + }
  1009 + }
  1010 + return result;
  1011 + }
  1012 +
  1013 + /// Calculate the average magnitude across the entire network.
  1014 + /// @param m is the magnitude value to use. The default is 0 (usually radius).
  1015 + T average(unsigned m = 0){
  1016 +
  1017 + T M, L; //allocate space for the total magnitude and length
  1018 + M = L = 0; //initialize both the initial magnitude and length to zero
  1019 + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
  1020 + M += E[e].integrate(); //get the integrated magnitude
  1021 + L += E[e].length(); //get the edge length
  1022 + }
  1023 +
  1024 + return M / L; //return the average magnitude
  1025 + }
  1026 +
  1027 + /// This function compares two networks and returns the percentage of the current network that is missing from A.
  1028 + /// @param A is the network to compare to - the field is generated for A
  1029 + /// @param sigma is the user-defined tolerance value - smaller values provide a stricter comparison
  1030 + stim::network<T> compare(stim::network<T> A, float sigma, int device = -1){
  1031 +
  1032 + stim::network<T> R; //generate a network storing the result of the comparison
  1033 + R = (*this); //initialize the result with the current network
  1034 +
  1035 + T *c; // centerline (array of double pointers) - points on kdtree must be double
  1036 + size_t n_data = A.total_points(); // set the number of points
  1037 + c = (T*) malloc(sizeof(T) * n_data * 3); // allocate an array to store all points in the data set
  1038 +
  1039 + unsigned t = 0;
  1040 + for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in the network
  1041 + for(unsigned p = 0; p < A.E[e].size(); p++){ //for each point in the edge
  1042 + for(unsigned d = 0; d < 3; d++){ //for each coordinate
  1043 +
  1044 + c[t * 3 + d] = A.E[e][p][d]; //copy the point into the array c
  1045 + }
  1046 + t++;
  1047 + }
  1048 + }
  1049 +
  1050 + //generate a KD-tree for network A
  1051 + size_t MaxTreeLevels = 3; // max tree level
  1052 +
  1053 +#ifdef __CUDACC__
  1054 + cudaSetDevice(device);
  1055 + int current_device;
  1056 + if (cudaGetDevice(&current_device) == device) {
  1057 + std::cout << "Using CUDA device " << device << " for calculations..." << std::endl;
  1058 + }
  1059 + stim::kdtree<T, 3> kdt; // initialize a pointer to a kd tree
  1060 +
  1061 + kdt.create(c, n_data, MaxTreeLevels); // build a KD tree
  1062 +
  1063 + for(unsigned e = 0; e < R.E.size(); e++){ //for each edge in A
  1064 + //R.E[e].add_mag(0); //add a new magnitude for the metric
  1065 + //size_t errormag_id = R.E[e].nmags() - 1; //get the id for the new magnitude
  1066 +
  1067 + size_t n = R.E[e].size(); // the number of points in current edge
  1068 + T* queryPt = new T[3 * n];
  1069 + T* m1 = new T[n];
  1070 + T* dists = new T[n];
  1071 + size_t* nnIdx = new size_t[n];
  1072 +
  1073 + T* d_dists;
  1074 + T* d_m1;
  1075 + cudaMalloc((void**)&d_dists, n * sizeof(T));
  1076 + cudaMalloc((void**)&d_m1, n * sizeof(T));
  1077 +
  1078 + edge2array(queryPt, R.E[e]);
  1079 + kdt.search(queryPt, n, nnIdx, dists);
  1080 +
  1081 + cudaMemcpy(d_dists, dists, n * sizeof(T), cudaMemcpyHostToDevice); // copy dists from host to device
  1082 +
  1083 + // configuration parameters
  1084 + size_t threads = (1024>n)?n:1024;
  1085 + size_t blocks = n/threads + (n%threads)?1:0;
  1086 +
  1087 + find_metric_parallel<<<blocks, threads>>>(d_m1, n, d_dists, sigma); //calculate the metric value based on the distance
  1088 +
  1089 + cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost);
  1090 +
  1091 + for(unsigned p = 0; p < n; p++){
  1092 + R.E[e].set_r(p, m1[p]);
  1093 + }
  1094 +
  1095 + //d_set_mag<<<blocks, threads>>>(R.E[e].M, errormag_id, n, m1);
  1096 + }
  1097 +
  1098 +#else
  1099 + stim::kdtree<T, 3> kdt;
  1100 + kdt.create(c, n_data, MaxTreeLevels);
  1101 +
  1102 + for(unsigned e = 0; e < R.E.size(); e++){ //for each edge in A
  1103 +
  1104 + size_t n = R.E[e].size(); // the number of points in current edge
  1105 + T* query = new T[3 * n];
  1106 + T* m1 = new T[n];
  1107 + T* dists = new T[n];
  1108 + size_t* nnIdx = new size_t[n];
  1109 +
  1110 + edge2array(query, R.E[e]);
  1111 +
  1112 + kdt.cpu_search(query, n, nnIdx, dists); //find the distance between A and the current network
  1113 +
  1114 + for(unsigned p = 0; p < R.E[e].size(); p++){
  1115 + m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma); //calculate the metric value based on the distance
  1116 + R.E[e].set_r(p, m1[p]); //set the error for the second point in the segment
  1117 + }
  1118 + }
  1119 +#endif
  1120 + return R; //return the resulting network
  1121 + }
  1122 +
  1123 + /// This function compares two networks and split the current one according to the nearest neighbor of each point in each edge
  1124 + /// @param A is the network to split
  1125 + /// @param B is the corresponding mapping network
  1126 + /// @param sigma is the user-defined tolerance value - smaller values provide a stricter comparison
  1127 + /// @param device is the device that user want to use
  1128 + void split(stim::network<T> A, stim::network<T> B, float sigma, int device, float threshold){
  1129 +
  1130 + T *c;
  1131 + size_t n_data = B.total_points();
  1132 + c = (T*) malloc(sizeof(T) * n_data * 3);
  1133 +
  1134 + size_t NB = B.E.size(); // the number of edges in B
  1135 + unsigned t = 0;
  1136 + for(unsigned e = 0; e < NB; e++){ // for every edge in B
  1137 + for(unsigned p = 0; p < B.E[e].size(); p++){ // for every points in B.E[e]
  1138 + for(unsigned d = 0; d < 3; d++){
  1139 +
  1140 + c[t * 3 + d] = B.E[e][p][d]; // convert to array
  1141 + }
  1142 + t++;
  1143 + }
  1144 + }
  1145 + size_t MaxTreeLevels = 3; // max tree level
  1146 +
  1147 +#ifdef __CUDACC__
  1148 + cudaSetDevice(device);
  1149 + stim::kdtree<T, 3> kdt; // initialize a pointer to a kd tree
  1150 +
  1151 + //compare each point in the current network to the field produced by A
  1152 + kdt.create(c, n_data, MaxTreeLevels); // build a KD tree
  1153 +
  1154 + std::vector<std::vector<unsigned> > relation; // the relationship between GT and T corresponding to NN
  1155 + relation.resize(A.E.size());
  1156 +
  1157 + for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in A
  1158 + //A.E[e].add_mag(0); //add a new magnitude for the metric
  1159 + //size_t errormag_id = A.E[e].nmags() - 1;
  1160 +
  1161 + size_t n = A.E[e].size(); // the number of edges in A
  1162 +
  1163 + T* queryPt = new T[3 * n]; // set of all the points in current edge
  1164 + T* m1 = new T[n]; // array of metrics for every point in current edge
  1165 + T* dists = new T[n]; // store the distances for every point in current edge
  1166 + size_t* nnIdx = new size_t[n]; // store the indices for every point in current edge
  1167 +
  1168 + // define pointers in device
  1169 + T* d_dists;
  1170 + T* d_m1;
  1171 + size_t* d_nnIdx;
  1172 +
  1173 + // allocate memory for defined pointers
  1174 + cudaMalloc((void**)&d_dists, n * sizeof(T));
  1175 + cudaMalloc((void**)&d_m1, n * sizeof(T));
  1176 + cudaMalloc((void**)&d_nnIdx, n * sizeof(size_t));
  1177 +
  1178 + edge2array(queryPt, A.E[e]); // convert edge to array
  1179 + kdt.search(queryPt, n, nnIdx, dists); // search the tree to find the NN for every point in current edge
  1180 +
  1181 + cudaMemcpy(d_dists, dists, n * sizeof(T), cudaMemcpyHostToDevice); // copy dists from host to device
  1182 + cudaMemcpy(d_nnIdx, nnIdx, n * sizeof(size_t), cudaMemcpyHostToDevice); // copy Idx from host to device
  1183 +
  1184 + // configuration parameters
  1185 + size_t threads = (1024>n)?n:1024; // test to see whether the number of point in current edge is more than 1024
  1186 + size_t blocks = n/threads + (n%threads)?1:0;
  1187 +
  1188 + find_metric_parallel<<<blocks, threads>>>(d_m1, n, d_dists, sigma); // calculate the metrics in parallel
  1189 +
  1190 + cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost);
  1191 +
  1192 + for(unsigned p = 0; p < n; p++){
  1193 + A.E[e].set_r(p, m1[p]); // set the error(radius) value to every point in current edge
  1194 + }
  1195 +
  1196 + relation[e].resize(n); // resize every edge relation size
  1197 +
  1198 + unsigned* d_relation;
  1199 + cudaMalloc((void**)&d_relation, n * sizeof(unsigned)); // allocate memory
  1200 +
  1201 + std::vector<size_t> edge_point_num(NB); // %th element is the number of points that %th edge has
  1202 + for(unsigned ee = 0; ee < NB; ee++)
  1203 + edge_point_num[ee] = B.E[ee].size();
  1204 +
  1205 + size_t* d_edge_point_num;
  1206 + cudaMalloc((void**)&d_edge_point_num, NB * sizeof(size_t));
  1207 + cudaMemcpy(d_edge_point_num, &edge_point_num[0], NB * sizeof(size_t), cudaMemcpyHostToDevice);
  1208 +
  1209 + find_edge_index_parallel<<<blocks, threads>>>(d_nnIdx, n, d_relation, d_edge_point_num, NB); // find the edge corresponding to the array index in parallel
  1210 +
  1211 + cudaMemcpy(&relation[e][0], d_relation, n * sizeof(unsigned), cudaMemcpyDeviceToHost); //copy relationship from device to host
  1212 + }
  1213 +#else
  1214 + stim::kdtree<T, 3> kdt;
  1215 + kdt.create(c, n_data, MaxTreeLevels);
  1216 +
  1217 + std::vector<std::vector<unsigned>> relation; // the mapping relationship between two networks
  1218 + relation.resize(A.E.size());
  1219 + for(unsigned i = 0; i < A.E.size(); i++)
  1220 + relation[i].resize(A.E[i].size());
  1221 +
  1222 + std::vector<size_t> edge_point_num(NB); //%th element is the number of points that %th edge has
  1223 + for(unsigned ee = 0; ee < NB; ee++)
  1224 + edge_point_num[ee] = B.E[ee].size();
  1225 +
  1226 + for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in A
  1227 +
  1228 + size_t n = A.E[e].size(); //the number of edges in A
  1229 +
  1230 + T* queryPt = new T[3 * n];
  1231 + T* m1 = new T[n];
  1232 + T* dists = new T[n]; //store the distances
  1233 + size_t* nnIdx = new size_t[n]; //store the indices
  1234 +
  1235 + edge2array(queryPt, A.E[e]);
  1236 + kdt.search(queryPt, n, nnIdx, dists);
  1237 +
  1238 + for(unsigned p = 0; p < A.E[e].size(); p++){
  1239 + m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma); //calculate the metric value based on the distance
  1240 + A.E[e].set_r(p, m1[p]); //set the error for the second point in the segment
  1241 +
  1242 + unsigned id = 0; //mapping edge's idx
  1243 + size_t num = 0; //total number of points before #th edge
  1244 + for(unsigned i = 0; i < NB; i++){
  1245 + num += B.E[i].size();
  1246 + if(nnIdx[p] < num){ //find the edge it belongs to
  1247 + relation[e][p] = id;
  1248 + break;
  1249 + }
  1250 + id++; //current edge won't be the one, move to next edge
  1251 + }
  1252 + }
  1253 + }
  1254 +#endif
  1255 + E = A.E;
  1256 + V = A.V;
  1257 +
  1258 + unsigned int id = 0; // split value
  1259 + for(unsigned e = 0; e < E.size(); e++){ // for every edge
  1260 + for(unsigned p = 0; p < E[e].size() - 1; p++){ // for every point in each edge
  1261 + int t = (int)(E[e].length() / sigma * 2);
  1262 + if (t <= 20)
  1263 + threshold_fac = E[e].size();
  1264 + else
  1265 + threshold_fac = (E[e].length() / sigma * 2)/10;
  1266 + if(relation[e][p] != relation[e][p + 1]){ // find the nearest edge changing point
  1267 + id = p + 1; // if there is no change in NN
  1268 + if(id < threshold_fac || (E[e].size() - id) < threshold_fac)
  1269 + id = E[e].size() - 1; // extreme situation is not acceptable
  1270 + else
  1271 + break;
  1272 + }
  1273 + if(p == E[e].size() - 2) // if there is no splitting index, set the id to the last point index of current edge
  1274 + id = E[e].size() - 1;
  1275 + }
  1276 + //unsigned errormag_id = E[e].nmags() - 1;
  1277 + T G = 0; // test to see whether it has its nearest neighbor
  1278 + for(unsigned i = 0; i < E[e].size(); i++)
  1279 + G += E[e].r(i); // won't split special edges
  1280 + if(G / E[e].size() > threshold) // should based on the color map
  1281 + id = E[e].size() - 1; // set split idx to outgoing direction vertex
  1282 +
  1283 + std::vector<edge> tmpe;
  1284 + tmpe.resize(2);
  1285 + tmpe = E[e].split(id);
  1286 + vertex tmpv = stim::vec3<T>(-1, -1, 0); // store the split point as vertex
  1287 + if(tmpe.size() == 2){
  1288 + relation.resize(relation.size() + 1);
  1289 + for(unsigned d = id; d < E[e].size(); d++)
  1290 + relation[relation.size() - 1].push_back(relation[e][d]);
  1291 + tmpe[0].v[0] = E[e].v[0]; // begining vertex of first half edge -> original begining vertex
  1292 + tmpe[1].v[1] = E[e].v[1]; // ending vertex of second half edge -> original ending vertex
  1293 + tmpv = E[e][id];
  1294 + V.push_back(tmpv);
  1295 + tmpe[0].v[1] = (unsigned)V.size() - 1; // ending vertex of first half edge -> new vertex
  1296 + tmpe[1].v[0] = (unsigned)V.size() - 1; // begining vertex of second half edge -> new vertex
  1297 + edge tmp(E[e]);
  1298 + E[e] = tmpe[0]; // replace original edge by first half edge
  1299 + E.push_back(tmpe[1]); // push second half edge to the last
  1300 + V[V.size() - 1].e[1].push_back(e); // push first half edge to the incoming of new vertex
  1301 + V[V.size() - 1].e[0].push_back((unsigned)E.size() - 1); // push second half edge to the outgoing of new vertex
  1302 + for(unsigned i = 0; i < V[tmp.v[1]].e[1].size(); i++) // find the incoming edge of original ending vertex
  1303 + if(V[tmp.v[1]].e[1][i] == e)
  1304 + V[tmp.v[1]].e[1][i] = (unsigned)E.size() - 1; // set to new edge
  1305 + }
  1306 + }
  1307 + }
  1308 +
  1309 + /// This function compares two splitted networks and yields a mapping relationship between them according to NN
  1310 + /// @param B is the network that the current network is going to map to
  1311 + /// @param C is the mapping relationship: C[e1] = _e1 means e1 edge in current network is mapping the _e1 edge in B
  1312 + /// @param device is the device that user want to use
  1313 + void mapping(stim::network<T> B, std::vector<unsigned> &C, int device, float threshold){
  1314 + stim::network<T> A; //generate a network storing the result of the comparison
  1315 + A = (*this);
  1316 +
  1317 + size_t n = A.E.size(); // the number of edges in A
  1318 + size_t NB = B.E.size(); // the number of edges in B
  1319 +
  1320 + C.resize(A.E.size());
  1321 +
  1322 + T *c; // centerline (array of double pointers) - points on kdtree must be double
  1323 + size_t n_data = B.total_points(); // set the number of points
  1324 + c = (T*) malloc(sizeof(T) * n_data * 3);
  1325 +
  1326 + unsigned t = 0;
  1327 + for(unsigned e = 0; e < NB; e++){ // for each edge in the network
  1328 + for(unsigned p = 0; p < B.E[e].size(); p++){ // for each point in the edge
  1329 + for(unsigned d = 0; d < 3; d++){ // for each coordinate
  1330 +
  1331 + c[t * 3 + d] = B.E[e][p][d];
  1332 + }
  1333 + t++;
  1334 + }
  1335 + }
  1336 +
  1337 + //generate a KD-tree for network A
  1338 + //float metric = 0.0; // initialize metric to be returned after comparing the network
  1339 + size_t MaxTreeLevels = 3; // max tree level
  1340 +
  1341 +#ifdef __CUDACC__
  1342 + cudaSetDevice(device);
  1343 + stim::kdtree<T, 3> kdt; // initialize a pointer to a kd tree
  1344 +
  1345 + kdt.create(c, n_data, MaxTreeLevels); // build a KD tree
  1346 +
  1347 + for(unsigned e = 0; e < n; e++){ //for each edge in A
  1348 + //size_t errormag_id = A.E[e].nmags() - 1; //get the id for the new magnitude
  1349 +
  1350 + //pre-judge to get rid of impossibly mapping edges
  1351 + T M = 0;
  1352 + for(unsigned p = 0; p < A.E[e].size(); p++)
  1353 + M += A.E[e].r(p);
  1354 + M = M / A.E[e].size();
  1355 + if(M > threshold)
  1356 + C[e] = (unsigned)-1; //set the nearest edge of impossibly mapping edges to maximum of unsigned
  1357 + else{
  1358 + T* queryPt = new T[3];
  1359 + T* dists = new T[1];
  1360 + size_t* nnIdx = new size_t[1];
  1361 +
  1362 + stim2array(queryPt, A.E[e][A.E[e].size()/2]);
  1363 + kdt.search(queryPt, 1, nnIdx, dists);
  1364 +
  1365 + unsigned id = 0; //mapping edge's idx
  1366 + size_t num = 0; //total number of points before #th edge
  1367 + for(unsigned i = 0; i < NB; i++){
  1368 + num += B.E[i].size();
  1369 + if(nnIdx[0] < num){
  1370 + C[e] = id;
  1371 + break;
  1372 + }
  1373 + id++;
  1374 + }
  1375 + }
  1376 + }
  1377 +#else
  1378 + stim::kdtree<T, 3> kdt;
  1379 + kdt.create(c, n_data, MaxTreeLevels);
  1380 + T *dists = new T[1]; // near neighbor distances
  1381 + size_t *nnIdx = new size_t[1]; // near neighbor indices // allocate near neigh indices
  1382 +
  1383 + stim::vec3<T> p0, p1;
  1384 + T* queryPt = new T[3];
  1385 +
  1386 + for(unsigned int e = 0; e < R.E.size(); e++){ // for each edge in A
  1387 + T M; // the sum of metrics of current edge
  1388 + for(unsigned p = 0; p < R.E[e].size(); p++)
  1389 + M += A.E[e].r(p);
  1390 + M = M / A.E[e].size();
  1391 + if(M > threshold)
  1392 + C[e] = (unsigned)-1;
  1393 + else{ // if it should have corresponding edge in B, then...
  1394 + p1 = R.E[e][R.E[e].size()/2];
  1395 + stim2array(queryPt, p1);
  1396 + kdt.cpu_search(queryPt, 1, nnIdx, dists); // search the tree
  1397 +
  1398 + unsigned id = 0; //mapping edge's idx
  1399 + size_t num = 0; //total number of points before #th edge
  1400 + for(unsigned i = 0; i < NB; i++){
  1401 + num += B.E[i].size();
  1402 + if(nnIdx[0] < num){
  1403 + C[e] = id;
  1404 + break;
  1405 + }
  1406 + id++;
  1407 + }
  1408 + }
  1409 + }
  1410 +#endif
  1411 + }
  1412 +
  1413 + /// Returns the number of magnitude values stored in each edge. This should be uniform across the network.
  1414 + //unsigned nmags(){
  1415 + // return E[0].nmags();
  1416 + //}
  1417 + // split a string in text by the character sep
  1418 + stim::vec<T> split(std::string &text, char sep)
  1419 + {
  1420 + stim::vec<T> tokens;
  1421 + std::size_t start = 0, end = 0;
  1422 + while ((end = text.find(sep, start)) != std::string::npos) {
  1423 + tokens.push_back(atof(text.substr(start, end - start).c_str()));
  1424 + start = end + 1;
  1425 + }
  1426 + tokens.push_back(atof(text.substr(start).c_str()));
  1427 + return tokens;
  1428 + }
  1429 + // load a network in text file to a network class
  1430 + void load_txt(std::string filename)
  1431 + {
  1432 + std::vector <std::string> file_contents;
  1433 + std::ifstream file(filename.c_str());
  1434 + std::string line;
  1435 + std::vector<unsigned> id2vert; //this list stores the vertex ID associated with each network vertex
  1436 + //for each line in the text file, store them as strings in file_contents
  1437 + while (std::getline(file, line))
  1438 + {
  1439 + std::stringstream ss(line);
  1440 + file_contents.push_back(ss.str());
  1441 + }
  1442 + unsigned int numEdges = atoi(file_contents[0].c_str()); //number of edges in the network
  1443 + unsigned int I = atoi(file_contents[1].c_str()) ; //calculate the number of points3d on the first edge
  1444 + unsigned int count = 1; unsigned int k = 2; // count is global counter through the file contents, k is for the vertices on the edges
  1445 + // for each edge in the network.
  1446 + for (unsigned int i = 0; i < numEdges; i ++ )
  1447 + {
  1448 + // pre allocate a position vector p with number of points3d on the edge p
  1449 + std::vector< stim::vec<T> > p(0, I);
  1450 + // for each point on the nth edge
  1451 + for (unsigned int j = k; j < I + k; j++)
  1452 + {
  1453 + // split the points3d of floats with separator space and form a float3 position vector out of them
  1454 + p.push_back(split(file_contents[j], ' '));
  1455 + }
  1456 + count += p.size() + 1; // increment count to point at the next edge in the network
  1457 + I = atoi(file_contents[count].c_str()); // read in the points3d at the next edge and convert it to an integer
  1458 + k = count + 1;
  1459 + edge new_edge = p; // create an edge with a vector of points3d on the edge
  1460 + E.push_back(new_edge); // push the edge into the network
  1461 + }
  1462 + unsigned int numVertices = atoi(file_contents[count].c_str()); // this line in the text file gives the number of distinct vertices
  1463 + count = count + 1; // this line of text file gives the first verrtex
  1464 + // push each vertex into V
  1465 + for (unsigned int i = 0; i < numVertices; i ++)
  1466 + {
  1467 + vertex new_vertex = split(file_contents[count], ' ');
  1468 + V.push_back(new_vertex);
  1469 + count += atoi(file_contents[count + 1].c_str()) + 2; // Skip number of edge ids + 2 to point to the next vertex
  1470 + }
  1471 + } // end load_txt function
  1472 +
  1473 + // strTxt returns a string of edges
  1474 + std::string
  1475 + strTxt(std::vector< stim::vec<T> > p)
  1476 + {
  1477 + std::stringstream ss;
  1478 + std::stringstream oss;
  1479 + for(unsigned int i = 0; i < p.size(); i++){
  1480 + ss.str(std::string());
  1481 + for(unsigned int d = 0; d < 3; d++){
  1482 + ss<<p[i][d];
  1483 + }
  1484 + ss << "\n";
  1485 + }
  1486 + return ss.str();
  1487 + }
  1488 + // removes specified character from string
  1489 + void removeCharsFromString(std::string &str, char* charsToRemove ) {
  1490 + for ( unsigned int i = 0; i < strlen(charsToRemove); ++i ) {
  1491 + str.erase( remove(str.begin(), str.end(), charsToRemove[i]), str.end() );
  1492 + }
  1493 + }
  1494 + //exports network to txt file
  1495 + void
  1496 + to_txt(std::string filename)
  1497 + {
  1498 + std::ofstream ofs(filename.c_str(), std::ofstream::out | std::ofstream::app);
  1499 + //int num;
  1500 + ofs << (E.size()).str() << "\n";
  1501 + for(unsigned int i = 0; i < E.size(); i++)
  1502 + {
  1503 + std::string str;
  1504 + ofs << (E[i].size()).str() << "\n";
  1505 + str = E[i].strTxt();
  1506 + ofs << str << "\n";
  1507 + }
  1508 + for(int i = 0; i < V.size(); i++)
  1509 + {
  1510 + std::string str;
  1511 + str = V[i].str();
  1512 + char temp[4] = "[],";
  1513 + removeCharsFromString(str, temp);
  1514 + ofs << str << "\n";
  1515 + }
  1516 + ofs.close();
  1517 + }
  1518 +}; //end stim::network class
  1519 +}; //end stim namespace
  1520 +#endif