From f1bb77987f4b4912b5327d148963de4487a4f52b Mon Sep 17 00:00:00 2001 From: David Mayerich Date: Tue, 10 Jul 2018 14:11:27 -0500 Subject: [PATCH] moved the static network class into the netmets repository from STIMLIB, where it is no longer used --- centerline.h | 611 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ cylinder.h | 731 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ gl_network.h | 687 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ main.cu | 6 +++--- network.h | 1520 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 3552 insertions(+), 3 deletions(-) create mode 100644 centerline.h create mode 100644 cylinder.h create mode 100644 gl_network.h create mode 100644 network.h diff --git a/centerline.h b/centerline.h new file mode 100644 index 0000000..b4b3d3f --- /dev/null +++ b/centerline.h @@ -0,0 +1,611 @@ +#ifndef STIM_CENTERLINE_H +#define STIM_CENTERLINE_H + +#include +#include +#include + +namespace stim{ + +/** This class stores information about a single fiber represented as a set of geometric points + * between two branch or end points. This class is used as a fundamental component of the stim::network + * class to describe an interconnected (often biological) network. + */ +template +class centerline : public std::vector< stim::vec3 >{ + +protected: + + std::vector L; //stores the integrated length along the fiber (used for parameterization) + + ///Return the normalized direction vector at point i (average of the incoming and outgoing directions) + vec3 d(size_t i) { + if (size() <= 1) return vec3(0, 0, 0); //if there is insufficient information to calculate the direction, return a null vector + 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 + if (i == 0) return (at(1) - at(0)).norm(); //the first direction vector is oriented towards the first line segment + if (i == size() - 1) return (at(size() - 1) - at(size() - 2)).norm(); //the last direction vector is oriented towards the last line segment + + //all other direction vectors are the average direction of the two joined line segments + vec3 a = at(i) - at(i - 1); + vec3 b = at(i + 1) - at(i); + vec3 ab = a.norm() + b.norm(); + return ab.norm(); + } + + //initializes the integrated length vector to make parameterization easier, starting with index idx (all previous indices are assumed to be correct) + void update_L(size_t start = 0) { + L.resize(size()); //allocate space for the L array + if (start == 0) { + L[0] = 0; //initialize the length value for the first point to zero (0) + start++; + } + + stim::vec3 d; + for (size_t i = start; i < size(); i++) { //for each line segment in the centerline + d = at(i) - at(i - 1); + L[i] = L[i - 1] + d.len(); //calculate the running length total + } + } + + void init() { + if (size() == 0) return; //return if there aren't any points + update_L(); + } + + /// Returns a stim::vec representing the point at index i + + /// @param i is an index of the desired centerline point + stim::vec get_vec(unsigned i){ + return std::vector< stim::vec3 >::at(i); + } + + ///finds the index of the point closest to the length l on the lower bound. + ///binary search. + size_t findIdx(T l) { + for (size_t i = 1; i < L.size(); i++) { //for each point in the centerline + if (L[i] > l) return i - 1; //if we have passed the desired length value, return i-1 + } + return L.size() - 1; + /*size_t i = L.size() / 2; + size_t max = L.size() - 1; + size_t min = 0; + while (i < L.size() - 1){ + if (l < L[i]) { + max = i; + i = min + (max - min) / 2; + } + else if (L[i] <= l && L[i + 1] >= l) { + break; + } + else { + min = i; + i = min + (max - min) / 2; + } + } + return i;*/ + } + + ///Returns a position vector at the given length into the fiber (based on the pvalue). + ///Interpolates the radius along the line. + ///@param l: the location of the in the cylinder. + ///@param idx: integer location of the point closest to l but prior to it. + stim::vec3 p(T l, int idx) { + T rat = (l - L[idx]) / (L[idx + 1] - L[idx]); + stim::vec3 v1 = at(idx); + stim::vec3 v2 = at(idx + 1); + return(v1 + (v2 - v1)*rat); + } + + +public: + + using std::vector< stim::vec3 >::at; + using std::vector< stim::vec3 >::size; + + centerline() : std::vector< stim::vec3 >() { + init(); + } + centerline(size_t n) : std::vector< stim::vec3 >(n){ + init(); + } + centerline(std::vector > pos) : + std::vector > (pos) + { + init(); + } + + //overload the push_back function to update the length vector + void push_back(stim::vec3 p) { + std::vector< stim::vec3 >::push_back(p); + update_L(size() - 1); + } + + ///Returns a position vector at the given p-value (p value ranges from 0 to 1). + ///interpolates the position along the line. + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1). + stim::vec3 p(T pvalue) { + if (pvalue <= 0.0) return at(0); //return the first element + if (pvalue >= 1.0) return back(); //return the last element + + T l = pvalue*L[L.size() - 1]; + int idx = findIdx(l); + return p(l, idx); + } + + ///Update centerline internal parameters (currently the L vector) + void update() { + init(); + } + ///Return the length of the entire centerline + T length() { + return L.back(); + } + + + /// stitch two centerlines + ///@param c1, c2: two centerlines + ///@param sigma: sample rate + static std::vector< stim::centerline > stitch(stim::centerline c1, stim::centerline c2 = stim::centerline()) { + + std::vector< stim::centerline > result; + stim::centerline new_centerline; + stim::vec3 new_vertex; + + // if only one centerline, stitch itself! + if (c2.size() == 0) { + size_t num = c1.size(); + size_t id = 100000; // store the downsample start position + T threshold; + if (num < 4) { // if the number of vertex is less than 4, do nothing + result.push_back(c1); + return result; + } + else { + // test geometry start vertex + stim::vec3 v1 = c1[1] - c1[0]; // vector from c1[0] to c1[1] + for (size_t p = 2; p < num; p++) { // 90° standard??? + stim::vec3 v2 = c1[p] - c1[0]; + float cosine = v2.dot(v1); + if (cosine < 0) { + id = p; + threshold = v2.len(); + break; + } + } + if (id != 100000) { // find a downsample position on the centerline + T* c; + c = (T*)malloc(sizeof(T) * (num - id) * 3); + for (size_t p = id; p < num; p++) { + for (size_t d = 0; d < 3; d++) { + c[p * 3 + d] = c1[p][d]; + } + } + stim::kdtree kdt; + kdt.create(c, num - id, 5); // create tree + + T* query = (T*)malloc(sizeof(T) * 3); + for (size_t d = 0; d < 3; d++) + query[d] = c1[0][d]; + size_t index; + T dist; + + kdt.search(query, 1, &index, &dist); + + free(query); + free(c); + + if (dist > threshold) { + result.push_back(c1); + } + else { + // the loop part + new_vertex = c1[index]; + new_centerline.push_back(new_vertex); + for (size_t p = 0; p < index + 1; p++) { + new_vertex = c1[p]; + new_centerline.push_back(new_vertex); + } + result.push_back(new_centerline); + new_centerline.clear(); + + // the tail part + for (size_t p = index; p < num; p++) { + new_vertex = c1[p]; + new_centerline.push_back(new_vertex); + } + result.push_back(new_centerline); + } + } + else { // there is one potential problem that two positions have to be stitched + // test geometry end vertex + stim::vec3 v1 = c1[num - 2] - c1[num - 1]; + for (size_t p = num - 2; p > 0; p--) { // 90° standard + stim::vec3 v2 = c1[p - 1] - c1[num - 1]; + float cosine = v2.dot(v1); + if (cosine < 0) { + id = p; + threshold = v2.len(); + break; + } + } + if (id != 100000) { // find a downsample position + T* c; + c = (T*)malloc(sizeof(T) * (id + 1) * 3); + for (size_t p = 0; p < id + 1; p++) { + for (size_t d = 0; d < 3; d++) { + c[p * 3 + d] = c1[p][d]; + } + } + stim::kdtree kdt; + kdt.create(c, id + 1, 5); // create tree + + T* query = (T*)malloc(sizeof(T) * 1 * 3); + for (size_t d = 0; d < 3; d++) + query[d] = c1[num - 1][d]; + size_t index; + T dist; + + kdt.search(query, 1, &index, &dist); + + free(query); + free(c); + + if (dist > threshold) { + result.push_back(c1); + } + else { + // the tail part + for (size_t p = 0; p < index + 1; p++) { + new_vertex = c1[p]; + new_centerline.push_back(new_vertex); + } + result.push_back(new_centerline); + new_centerline.clear(); + + // the loop part + for (size_t p = index; p < num; p++) { + new_vertex = c1[p]; + new_centerline.push_back(new_vertex); + } + new_vertex = c1[index]; + new_centerline.push_back(new_vertex); + result.push_back(new_centerline); + } + } + else { // no stitch position + result.push_back(c1); + } + } + } + } + + + // two centerlines + else { + // find stitch position based on nearest neighbors + size_t num1 = c1.size(); + T* c = (T*)malloc(sizeof(T) * num1 * 3); // c1 as reference point + for (size_t p = 0; p < num1; p++) // centerline to array + for (size_t d = 0; d < 3; d++) // because right now my kdtree code is a relatively close code, it has its own structure + c[p * 3 + d] = c1[p][d]; // I will merge it into stimlib totally in the near future + + stim::kdtree kdt; // kdtree object + kdt.create(c, num1, 5); // create tree + + size_t num2 = c2.size(); + T* query = (T*)malloc(sizeof(T) * num2 * 3); // c2 as query point + for (size_t p = 0; p < num2; p++) { + for (size_t d = 0; d < 3; d++) { + query[p * 3 + d] = c2[p][d]; + } + } + std::vector index(num2); + std::vector dist(num2); + + kdt.search(query, num2, &index[0], &dist[0]); // find the nearest neighbors in c1 for c2 + + // clear up + free(query); + free(c); + + // find the average vertex distance of one centerline + T sigma1 = 0; + T sigma2 = 0; + for (size_t p = 0; p < num1 - 1; p++) + sigma1 += (c1[p] - c1[p + 1]).len(); + for (size_t p = 0; p < num2 - 1; p++) + sigma2 += (c2[p] - c2[p + 1]).len(); + sigma1 /= (num1 - 1); + sigma2 /= (num2 - 1); + float threshold = 4 * (sigma1 + sigma2) / 2; // better way to do this? + + T min_d = *std::min_element(dist.begin(), dist.end()); // find the minimum distance between c1 and c2 + + if (min_d > threshold) { // if the minimum distance is too large + result.push_back(c1); + result.push_back(c2); + +#ifdef DEBUG + std::cout << "The distance between these two centerlines is too large" << std::endl; +#endif + } + else { + // auto smallest = std::min_element(dist.begin(), dist.end()); + unsigned int smallest = std::min_element(dist.begin(), dist.end()); + // auto i = std::distance(dist.begin(), smallest); // find the index of min-distance in distance list + unsigned int i = std::distance(dist.begin(), smallest); // find the index of min-distance in distance list + + // stitch position in c1 and c2 + int id1 = index[i]; + int id2 = i; + + // actually there are two cases + // first one inacceptable + // second one acceptable + if (id1 != 0 && id1 != num1 - 1 && id2 != 0 && id2 != num2 - 1) { // only stitch one end vertex to another centerline + result.push_back(c1); + result.push_back(c2); + } + else { + if (id1 == 0 || id1 == num1 - 1) { // if the stitch vertex is the first or last vertex of c1 + // for c2, consider two cases(one degenerate case) + if (id2 == 0 || id2 == num2 - 1) { // case 1, if stitch position is also on the end of c2 + // we have to decide which centerline get a new vertex, based on direction + // for c1, computer the direction change angle + stim::vec3 v1, v2; + float alpha1, alpha2; // direction change angle + if (id1 == 0) + v1 = (c1[1] - c1[0]).norm(); + else + v1 = (c1[num1 - 2] - c1[num1 - 1]).norm(); + v2 = (c2[id2] - c1[id1]).norm(); + alpha1 = v1.dot(v2); + if (id2 == 0) + v1 = (c2[1] - c2[0]).norm(); + else + v1 = (c2[num2 - 2] - c2[num2 - 1]).norm(); + v2 = (c1[id1] - c2[id2]).norm(); + alpha2 = v1.dot(v2); + if (abs(alpha1) > abs(alpha2)) { // add the vertex to c1 in order to get smooth connection + // push back c1 + if (id1 == 0) { // keep geometry information + new_vertex = c2[id2]; + new_centerline.push_back(new_vertex); + for (size_t p = 0; p < num1; p++) { // stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1 + new_vertex = c1[p]; + new_centerline.push_back(new_vertex); + } + } + else { + for (size_t p = 0; p < num1; p++) { // stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1 + new_vertex = c1[p]; + new_centerline.push_back(new_vertex); + } + new_vertex = c2[id2]; + new_centerline.push_back(new_vertex); + } + result.push_back(new_centerline); + new_centerline.clear(); + + // push back c2 + for (size_t p = 0; p < num2; p++) { + new_vertex = c2[p]; + new_centerline.push_back(new_vertex); + } + result.push_back(new_centerline); + } + else { // add the vertex to c2 in order to get smooth connection + // push back c1 + for (size_t p = 0; p < num1; p++) { + new_vertex = c1[p]; + new_centerline.push_back(new_vertex); + } + result.push_back(new_centerline); + new_centerline.clear(); + + // push back c2 + if (id2 == 0) { // keep geometry information + new_vertex = c1[id1]; + new_centerline.push_back(new_vertex); + for (size_t p = 0; p < num2; p++) { // stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1 + new_vertex = c2[p]; + new_centerline.push_back(new_vertex); + } + } + else { + for (size_t p = 0; p < num2; p++) { // stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1 + new_vertex = c2[p]; + new_centerline.push_back(new_vertex); + } + new_vertex = c1[id1]; + new_centerline.push_back(new_vertex); + } + result.push_back(new_centerline); + } + } + else { // case 2, the stitch position is on c2 + // push back c1 + if (id1 == 0) { // keep geometry information + new_vertex = c2[id2]; + new_centerline.push_back(new_vertex); + for (size_t p = 0; p < num1; p++) { // stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1 + new_vertex = c1[p]; + new_centerline.push_back(new_vertex); + } + } + else { + for (size_t p = 0; p < num1; p++) { // geometry end vertex on c1 -> geometry start vertex on c1 -> stitch vertex on c2 + new_vertex = c1[p]; + new_centerline.push_back(new_vertex); + } + new_vertex = c2[id2]; + new_centerline.push_back(new_vertex); + } + result.push_back(new_centerline); + new_centerline.clear(); + + // push back c2 + for (size_t p = 0; p < id2 + 1; p++) { // first part + new_vertex = c2[p]; + new_centerline.push_back(new_vertex); + } + result.push_back(new_centerline); + new_centerline.clear(); + + for (size_t p = id2; p < num2; p++) { // second part + new_vertex = c2[p]; + new_centerline.push_back(new_vertex); + } + result.push_back(new_centerline); + } + } + else { // if the stitch vertex is the first or last vertex of c2 + // push back c2 + if (id2 == 0) { // keep geometry information + new_vertex = c1[id1]; + new_centerline.push_back(new_vertex); + for (size_t p = 0; p < num2; p++) { // stitch vertex on c1 -> geometry start vertex on c2 -> geometry end vertex on c2 + new_vertex = c2[p]; + new_centerline.push_back(new_vertex); + } + } + else { + for (size_t p = 0; p < num2; p++) { // geometry end vertex on c2 -> geometry start vertex on c2 -> stitch vertex on c1 + new_vertex = c2[p]; + new_centerline.push_back(new_vertex); + } + new_vertex = c1[id1]; + new_centerline.push_back(new_vertex); + result.push_back(new_centerline); + new_centerline.clear(); + + // push back c1 + for (size_t p = 0; p < id1 + 1; p++) { // first part + new_vertex = c1[p]; + new_centerline.push_back(new_vertex); + } + result.push_back(new_centerline); + new_centerline.clear(); + + for (size_t p = id1; p < num1; p++) { // second part + new_vertex = c1[p]; + new_centerline.push_back(new_vertex); + } + result.push_back(new_centerline); + } + } + } + } + } + return result; + } + + /// Split the fiber at the specified index. If the index is an end point, only one fiber is returned + std::vector< stim::centerline > split(unsigned int idx){ + + std::vector< stim::centerline > fl; //create an array to store up to two fibers + size_t N = size(); + + //if the index is an end point, only the existing fiber is returned + if(idx == 0 || idx == N-1){ + fl.resize(1); //set the size of the fiber to 1 + fl[0] = *this; //copy the current fiber + } + + //if the index is not an end point + else{ + + unsigned int N1 = idx + 1; //calculate the size of both fibers + unsigned int N2 = N - idx; + + fl.resize(2); //set the array size to 2 + + fl[0] = stim::centerline(N1); //set the size of each fiber + fl[1] = stim::centerline(N2); + + //copy both halves of the fiber + unsigned int i; + + //first half + for(i = 0; i < N1; i++) //for each centerline point + fl[0][i] = std::vector< stim::vec3 >::at(i); + fl[0].init(); //initialize the length vector + + //second half + for(i = 0; i < N2; i++) + fl[1][i] = std::vector< stim::vec3 >::at(idx+i); + fl[1].init(); //initialize the length vector + } + + return fl; //return the array + + } + + /// Outputs the fiber as a string + std::string str(){ + std::stringstream ss; + size_t N = std::vector< stim::vec3 >::size(); + ss << "---------[" << N << "]---------" << std::endl; + for (size_t i = 0; i < N; i++) + ss << std::vector< stim::vec3 >::at(i) << std::endl; + ss << "--------------------" << std::endl; + + return ss.str(); + } + + /// Back method returns the last point in the fiber + stim::vec3 back(){ + return std::vector< stim::vec3 >::back(); + } + + ////resample a fiber in the network + stim::centerline resample(T spacing) + { + //std::cout<<"fiber::resample()"< v; //v-direction vector of the segment + stim::vec3 p; //- intermediate point to be added + stim::vec3 p1; // p1 - starting point of an segment on the fiber, + stim::vec3 p2; // p2 - ending point, + //double sum=0; //distance summation + + size_t N = size(); + + centerline new_c; // initialize list of new resampled points on the fiber + // for each point on the centerline (skip if it is the last point on centerline) + for(unsigned int f=0; f< N-1; f++) + { + p1 = at(f); + p2 = at(f+1); + v = p2 - p1; + + T lengthSegment = v.len(); //find Length of the segment as distance between the starting and ending points of the segment + + if(lengthSegment >= spacing){ // if length of the segment is greater than standard deviation resample + + // repeat resampling until accumulated stepsize is equsl to length of the segment + for(T step=0.0; step +#include +#include "centerline.h" +#include + + +namespace stim +{ +template +class cylinder : public centerline { +protected: + + using stim::centerline::d; + + std::vector< stim::vec3 > U; //stores the array of U vectors defining the Frenet frame + std::vector< T > R; //stores a list of magnitudes for each point in the centerline (assuming mags[0] is the radius) + + using stim::centerline::findIdx; + + void init() { + U.resize(size()); //allocate space for the frenet frame vectors +// if (R.size() != 0) + R.resize(size()); + + stim::circle c; //create a circle + stim::vec3 d0, d1; + for (size_t i = 0; i < size() - 1; i++) { //for each line segment in the centerline + c.rotate(d(i)); //rotate the circle to match that normal + U[i] = c.U; //save the U vector from the circle + } + U[size() - 1] = c.U; //for the last point, duplicate the final frenet frame vector + } + + //calculates the U values for each point to initialize the frenet frame + // this function assumes that the centerline has already been set + +public: + + using stim::centerline::size; + using stim::centerline::at; + using stim::centerline::L; + using stim::centerline::length; + + cylinder() : centerline(){} + + cylinder(centerlinec) : centerline(c) { + init(); + } + + cylinder(std::vector > p, std::vector s) + : centerline(p) + { + R = s; + init(); + } + + cylinder(stim::centerline p, std::vector s) + { + //was d = s; + p = s; + init(); + } + + //cylinder(centerlinec, T r) : centerline(c) { + // init(); + // //add_mag(r); + //} + + //initialize a cylinder with a list of points and magnitude values + //cylinder(centerlinec, std::vector r) : centerline(c){ + // init(); + // //add_mag(r); + //} + + //copy the original radius + void copy_r(std::vector radius) { + for (unsigned i = 0; i < radius.size(); i++) + R[i] = radius[i]; + } + + ///Returns magnitude i at the given length into the fiber (based on the pvalue). + ///Interpolates the position along the line. + ///@param l: the location of the in the cylinder. + ///@param idx: integer location of the point closest to l but prior to it. + T r(T l, int idx) { + T a = (l - L[idx]) / (L[idx + 1] - L[idx]); + T v2 = (R[idx] + (R[idx + 1] - R[idx])*a); + + return v2; + } + + ///Returns the ith magnitude at the given p-value (p value ranges from 0 to 1). + ///interpolates the radius along the line. + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1). + T rl(T pvalue) { + if (pvalue <= 0.0) return R[0]; + if (pvalue >= 1.0) return R[size() - 1]; + + T l = pvalue*L[L.size() - 1]; + int idx = findIdx(l); + return r(l, idx); + } + + /// Returns the magnitude at the given index + /// @param i is the index of the desired point + /// @param r is the index of the magnitude value + T r(unsigned i) { + return R[i]; + } + + + ///adds a magnitude to each point in the cylinder + /*void add_mag(V val = 0) { + if (M.size() == 0) M.resize(size()); //if the magnitude vector isn't initialized, resize it to match the centerline + for (size_t i = 0; i < size(); i++) //for each point + R[i].push_back(val); //add this value to the magnitude vector at each point + }*/ + + //adds a magnitude based on a list of magnitudes for each point + /*void add_mag(std::vector val) { + if (M.size() == 0) M.resize(size()); //if the magnitude vector isn't initialized, resize it to match the centerline + for (size_t i = 0; i < size(); i++) //for each point + R[i].push_back(val[i]); //add this value to the magnitude vector at each point + }*/ + + //sets the value of magnitude m at point i + void set_r(size_t i, T r) { + R[i] = r; + } + + /*size_t nmags() { + if (M.size() == 0) return 0; + else return R[0].size(); + }*/ + + ///Returns a circle representing the cylinder cross section at point i + stim::circle circ(size_t i) { + return stim::circle(at(i), R[i], d(i), U[i]); + } + + ///Returns an OBJ object representing the cylinder with a radial tesselation value of N using magnitude m + stim::obj OBJ(size_t N) { + stim::obj out; //create an OBJ object + stim::circle c0, c1; + std::vector< stim::vec3 > p0, p1; //allocate space for the point sets representing the circles bounding each cylinder segment + T u0, u1, v0, v1; //allocate variables to store running texture coordinates + for (size_t i = 1; i < size(); i++) { //for each line segment in the cylinder + c0 = circ(i - 1); //get the two circles bounding the line segment + c1 = circ(i); + + p0 = c0.points(N); //get t points for each of the end caps + p1 = c1.points(N); + + u0 = L[i - 1] / length(); //calculate the texture coordinate (u, v) where u runs along the cylinder length + u1 = L[i] / length(); + + for (size_t n = 1; n < N; n++) { //for each point in the circle + v0 = (double)(n-1) / (double)(N - 1); //v texture coordinate runs around the cylinder + v1 = (double)(n) / (double)(N - 1); + out.Begin(OBJ_FACE); //start a face (quad) + out.TexCoord(u0, v0); + 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 + out.TexCoord(u1, v0); + out.Vertex(p1[n - 1][0], p1[n - 1][1], p1[n - 1][2]); + + out.TexCoord(u0, v1); + out.Vertex(p1[n][0], p1[n][1], p1[n][2]); + out.TexCoord(u1, v1); + out.Vertex(p0[n][0], p0[n][1], p0[n][2]); + out.End(); + } + v0 = (double)(N - 2) / (double)(N - 1); //v texture coordinate runs around the cylinder + v1 = 1.0; + out.Begin(OBJ_FACE); + out.TexCoord(u0, v0); + 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 + out.TexCoord(u1, v0); + out.Vertex(p1[N - 1][0], p1[N - 1][1], p1[N - 1][2]); + + out.TexCoord(u0, v1); + out.Vertex(p1[0][0], p1[0][1], p1[0][2]); + out.TexCoord(u1, v1); + out.Vertex(p0[0][0], p0[0][1], p0[0][2]); + out.End(); + } + return out; + } + + std::string str() { + std::stringstream ss; + size_t N = std::vector< stim::vec3 >::size(); + ss << "---------[" << N << "]---------" << std::endl; + for (size_t i = 0; i < N; i++) + ss << std::vector< stim::vec3 >::at(i) << " r = " << R[i] << " u = " << U[i] << std::endl; + ss << "--------------------" << std::endl; + + return ss.str(); + } + + /// Integrates a magnitude value along the cylinder. + /// @param m is the magnitude value to be integrated (this is usually the radius) + T integrate() { + T sum = 0; //initialize the integral to zero + if (L.size() == 1) + return sum; + else { + T m0, m1; //allocate space for both magnitudes in a single segment + m0 = R[0]; //initialize the first point and magnitude to the first point in the cylinder + T len = L[1]; //allocate space for the segment length + + + for (unsigned p = 1; p < size(); p++) { //for every consecutive point in the cylinder + m1 = R[p]; + if (p > 1) len = (L[p] - L[p - 1]); //calculate the segment length using the L array + sum += (m0 + m1) / (T)2.0 * len; //add the average magnitude, weighted by the segment length + m0 = m1; //move to the next segment by shifting points + } + return sum; //return the integral + } + } + + /// Resamples the cylinder to provide a maximum distance of "spacing" between centerline points. All current + /// centerline points are guaranteed to exist in the new cylinder + /// @param spacing is the maximum spacing allowed between sample points + cylinder resample(T spacing) { + cylinder c = stim::centerline::resample(spacing); //resample the centerline and use it to create a new cylinder + + //size_t nm = nmags(); //get the number of magnitude values in the current cylinder + //if (nm > 0) { //if there are magnitude values + // std::vector magvec(nm, 0); //create a magnitude vector for a single point + // c.M.resize(c.size(), magvec); //allocate space for a magnitude vector at each point of the new cylinder + //} + + T l, t; + for (size_t i = 0; i < c.size(); i++) { //for each point in the new cylinder + l = c.L[i]; //get the length along the new cylinder + t = l / length(); //calculate the parameter value along the new cylinder + //for (size_t mag = 0; mag < nm; mag++) { //for each magnitude value + c.R[i] = r(t); //retrieve the interpolated magnitude from the current cylinder and store it in the new one + //} + } + return c; + } + + std::vector< stim::cylinder > split(unsigned int idx) { + + unsigned N = size(); + std::vector< stim::centerline > LL; + LL.resize(2); + LL = (*this).centerline::split(idx); + std::vector< stim::cylinder > C(LL.size()); + unsigned i = 0; + C[0] = LL[0]; + //C[0].R.resize(idx); + for (; i < idx + 1; i++) { + //for(unsigned d = 0; d < 3; d++) + //C[0][i][d] = LL[0].c[i][d]; + C[0].R[i] = R[i]; + //C[0].R[i].resize(1); + } + if (C.size() == 2) { + C[1] = LL[1]; + i--; + //C[1].M.resize(N - idx); + for (; i < N; i++) { + //for(unsigned d = 0; d < 3; d++) + //C[1][i - idx][d] = LL[1].c[i - idx][d]; + C[1].R[i - idx] = R[i]; + //C[1].M[i - idx].resize(1); + } + } + + return C; + } + + + /* + ///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and magnitudes (inM) + void + init(centerline inP, std::vector< std::vector > inM) { + M = inM; //the magnitude vector can be copied directly + (*this) = inP; //the centerline can be copied to this class directly + stim::vec3 v1; + stim::vec3 v2; + e.resize(inP.size()); + + norms.resize(inP.size()); + Us.resize(inP.size()); + + if(inP.size() < 2) + return; + + //calculate each L. + L.resize(inP.size()); //the number of precomputed lengths will equal the number of points + T temp = (T)0; //length up to that point + L[0] = temp; + for(size_t i = 1; i < L.size(); i++) + { + temp += (inP[i-1] - inP[i]).len(); + L[i] = temp; + } + + stim::vec3 dr = (inP[1] - inP[0]).norm(); + s = stim::circle(inP[0], inR[0][0], dr, stim::vec3(1,0,0)); + norms[0] = s.N; + e[0] = s; + Us[0] = e[0].U; + for(size_t i = 1; i < inP.size()-1; i++) + { + s.center(inP[i]); + v1 = (inP[i] - inP[i-1]).norm(); + v2 = (inP[i+1] - inP[i]).norm(); + dr = (v1+v2).norm(); + s.normal(dr); + + norms[i] = s.N; + + s.scale(inR[i][0]/inR[i-1][0]); + e[i] = s; + Us[i] = e[i].U; + } + + int j = inP.size()-1; + s.center(inP[j]); + dr = (inP[j] - inP[j-1]).norm(); + s.normal(dr); + + norms[j] = s.N; + + s.scale(inR[j][0]/inR[j-1][0]); + e[j] = s; + Us[j] = e[j].U; + } + + ///returns the direction vector at point idx. + stim::vec3 + d(int idx) + { + if(idx == 0) + { + stim::vec3 temp( + c[idx+1][0]-c[idx][0], + c[idx+1][1]-c[idx][1], + c[idx+1][2]-c[idx][2] + ); +// return (e[idx+1].P - e[idx].P).norm(); + return (temp.norm()); + } + else if(idx == N-1) + { + stim::vec3 temp( + c[idx][0]-c[idx+1][0], + c[idx][1]-c[idx+1][1], + c[idx][2]-c[idx+1][2] + ); + // return (e[idx].P - e[idx-1].P).norm(); + return (temp.norm()); + } + else + { +// return (e[idx+1].P - e[idx].P).norm(); +// stim::vec3 v1 = (e[idx].P-e[idx-1].P).norm(); + stim::vec3 v1( + c[idx][0]-c[idx-1][0], + c[idx][1]-c[idx-1][1], + c[idx][2]-c[idx-1][2] + ); + +// stim::vec3 v2 = (e[idx+1].P-e[idx].P).norm(); + stim::vec3 v2( + c[idx+1][0]-c[idx][0], + c[idx+1][1]-c[idx][1], + c[idx+1][2]-c[idx][2] + ); + + return (v1.norm()+v2.norm()).norm(); + } + // return e[idx].N; + + } + + stim::vec3 + d(T l, int idx) + { + if(idx == 0 || idx == N-1) + { + return norms[idx]; +// return e[idx].N; + } + else + { + + T rat = (l-L[idx])/(L[idx+1]-L[idx]); + return( norms[idx] + (norms[idx+1] - norms[idx])*rat); +// return( e[idx].N + (e[idx+1].N - e[idx].N)*rat); + } + } + + + ///finds the index of the point closest to the length l on the lower bound. + ///binary search. + int + findIdx(T l) + { + unsigned int i = L.size()/2; + unsigned int max = L.size()-1; + unsigned int min = 0; + while(i > 0 && i < L.size()-1) + { +// std::cerr << "Trying " << i << std::endl; +// std::cerr << "l is " << l << ", L[" << i << "]" << L[i] << std::endl; + if(l < L[i]) + { + max = i; + i = min+(max-min)/2; + } + else if(L[i] <= l && L[i+1] >= l) + { + break; + } + else + { + min = i; + i = min+(max-min)/2; + } + } + return i; + } + + public: + ///default constructor + cylinder() + // : centerline() + { + + } + + ///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder. + ///@param inP: Vector of stim vec3 composing the points of the centerline. + ///@param inM: Vector of stim vecs composing the radii of the centerline. + cylinder(std::vector > inP, std::vector > inM) + : centerline(inP) + { + init(inP, inM); + } + + ///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder. + ///@param inP: Vector of stim vec3 composing the points of the centerline. + ///@param inM: Vector of stim vecs composing the radii of the centerline. + cylinder(std::vector > inP, std::vector< T > inM) + : centerline(inP) + { + std::vector > temp; + stim::vec zero(0.0,0.0); + temp.resize(inM.size(), zero); + for(int i = 0; i < inM.size(); i++) + temp[i][0] = inR[i]; + init(inP, temp); + } + + + ///Constructor defines a cylinder with centerline inP and magnitudes of zero + ///@param inP: Vector of stim vec3 composing the points of the centerline + cylinder(std::vector< stim::vec3 > inP) + : centerline(inP) + { + std::vector< stim::vec > inM; //create an array of arbitrary magnitudes + + stim::vec zero; + zero.push_back(0); + + inM.resize(inP.size(), zero); //initialize the magnitude values to zero + init(inP, inM); + } + + //assignment operator creates a cylinder from a centerline (default radius is zero) + cylinder& operator=(stim::centerline c) { + init(c); + } + + ///Returns the number of points on the cylinder centerline + + unsigned int size(){ + return N; + } + + + ///Returns a position vector at the given p-value (p value ranges from 0 to 1). + ///interpolates the position along the line. + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1). + stim::vec3 + p(T pvalue) + { + if(pvalue < 0.0 || pvalue > 1.0) + { + return stim::vec3(-1,-1,-1); + } + T l = pvalue*L[L.size()-1]; + int idx = findIdx(l); + return (p(l,idx)); + } + + ///Returns a position vector at the given length into the fiber (based on the pvalue). + ///Interpolates the radius along the line. + ///@param l: the location of the in the cylinder. + ///@param idx: integer location of the point closest to l but prior to it. + stim::vec3 + p(T l, int idx) + { + T rat = (l-L[idx])/(L[idx+1]-L[idx]); + stim::vec3 v1( + c[idx][0], + c[idx][1], + c[idx][2] + ); + + stim::vec3 v2( + c[idx+1][0], + c[idx+1][1], + c[idx+1][2] + ); +// return( e[idx].P + (e[idx+1].P-e[idx].P)*rat); + return( v1 + (v2-v1)*rat); +// return( +// return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx]))); + } + + ///Returns a radius at the given p-value (p value ranges from 0 to 1). + ///interpolates the radius along the line. + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1). + T + r(T pvalue) + { + if(pvalue < 0.0 || pvalue > 1.0){ + std::cerr<<"Error, value "<= 10e-6) + { + std::cout << "-------------------------" << std::endl; + std::cout << e[idx].str() << std::endl << std::endl; + std::cout << Us[idx].str() << std::endl; + std::cout << (float)v1 - (float) v2 << std::endl; + std::cout << "failed" << std::endl; + } +// std::cout << e[idx].U.len() << " " << mags[idx][0] << std::endl; +// std::cout << v2 << std::endl; + return(v2); +// return (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat); + // ( + } + + /// Returns the magnitude at the given index + /// @param i is the index of the desired point + /// @param m is the index of the magnitude value + T ri(unsigned i, unsigned m = 0){ + return mags[i][m]; + } + + /// Adds a new magnitude value to all points + /// @param m is the starting value for the new magnitude + void add_mag(T m = 0){ + for(unsigned int p = 0; p < N; p++) + mags[p].push_back(m); + } + + /// Sets a magnitude value + /// @param val is the new value for the magnitude + /// @param p is the point index for the magnitude to be set + /// @param m is the index for the magnitude + void set_mag(T val, unsigned p, unsigned m = 0){ + mags[p][m] = val; + } + + /// Returns the number of magnitude values at each point + unsigned nmags(){ + return mags[0].size(); + } + + ///returns the position of the point with a given pvalue and theta on the surface + ///in x, y, z coordinates. Theta is in degrees from 0 to 360. + ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1). + ///@param theta: the angle to the point of a circle. + stim::vec3 + surf(T pvalue, T theta) + { + if(pvalue < 0.0 || pvalue > 1.0) + { + return stim::vec3(-1,-1,-1); + } else { + T l = pvalue*L[L.size()-1]; + int idx = findIdx(l); + stim::vec3 ps = p(l, idx); + T m = r(l, idx); + s = e[idx]; + s.center(ps); + s.normal(d(l, idx)); + s.scale(m/e[idx].U.len()); + return(s.p(theta)); + } + } + + ///returns a vector of points necessary to create a circle at every position in the fiber. + ///@param sides: the number of sides of each circle. + std::vector > > + getPoints(int sides) + { + std::vector > > points; + points.resize(N); + for(int i = 0; i < N; i++) + { + points[i] = e[i].getPoints(sides); + } + return points; + } + + ///returns the total length of the line at index j. + T + getl(int j) + { + return (L[j]); + } + /// Allows a point on the centerline to be accessed using bracket notation + + vec3 operator[](unsigned int i){ + return e[i].P; + } + + /// Returns the total length of the cylinder centerline + T length(){ + return L.back(); + } + + /// Integrates a magnitude value along the cylinder. + /// @param m is the magnitude value to be integrated (this is usually the radius) + T integrate(unsigned m = 0){ + + T M = 0; //initialize the integral to zero + T m0, m1; //allocate space for both magnitudes in a single segment + + //vec3 p0, p1; //allocate space for both points in a single segment + + m0 = mags[0][m]; //initialize the first point and magnitude to the first point in the cylinder + //p0 = pos[0]; + + T len = L[0]; //allocate space for the segment length + + //for every consecutive point in the cylinder + for(unsigned p = 1; p < N; p++){ + + //p1 = pos[p]; //get the position and magnitude for the next point + m1 = mags[p][m]; + + if(p > 1) len = (L[p-1] - L[p-2]); //calculate the segment length using the L array + + //add the average magnitude, weighted by the segment length + M += (m0 + m1)/(T)2.0 * len; + + m0 = m1; //move to the next segment by shifting points + } + return M; //return the integral + } + + /// Averages a magnitude value across the cylinder + /// @param m is the magnitude value to be averaged (this is usually the radius) + T average(unsigned m = 0){ + + //return the average magnitude + return integrate(m) / L.back(); + } + + /// Resamples the cylinder to provide a maximum distance of "spacing" between centerline points. All current + /// centerline points are guaranteed to exist in the new cylinder + /// @param spacing is the maximum spacing allowed between sample points + cylinder resample(T spacing){ + + std::vector< vec3 > result; + + vec3 p0 = e[0].P; //initialize p0 to the first point on the centerline + vec3 p1; + unsigned N = size(); //number of points in the current centerline + + //for each line segment on the centerline + for(unsigned int i = 1; i < N; i++){ + p1 = e[i].P; //get the second point in the line segment + + vec3 v = p1 - p0; //calculate the vector between these two points + T d = v.len(); //calculate the distance between these two points (length of the line segment) + + size_t nsteps = (size_t)std::ceil(d / spacing); //calculate the number of steps to take along the segment to meet the spacing criteria + T stepsize = (T)1.0 / nsteps; //calculate the parametric step size between new centerline points + + //for each step along the line segment + for(unsigned s = 0; s < nsteps; s++){ + T alpha = stepsize * s; //calculate the fraction of the distance along the line segment covered + result.push_back(p0 + alpha * v); //push the point at alpha position along the line segment + } + + p0 = p1; //shift the points to move to the next line segment + } + + result.push_back(e[size() - 1].P); //push the last point in the centerline + + return cylinder(result); + + }*/ + + +}; + +} +#endif diff --git a/gl_network.h b/gl_network.h new file mode 100644 index 0000000..adb95ba --- /dev/null +++ b/gl_network.h @@ -0,0 +1,687 @@ +#ifndef STIM_GL_NETWORK +#define STIM_GL_NETWORK + +#include +#include "network.h" +#include + +namespace stim{ + +template +class gl_network : public stim::network{ + +protected: + using stim::network::E; + using stim::network::V; + + GLuint dlist; + +public: + + /// Default constructor + gl_network() : stim::network(){ + dlist = 0; + } + + /// Constructor creates a gl_network from a stim::network + gl_network(stim::network N) : stim::network(N){ + dlist = 0; + } + + /// Fills the parameters with the minimum and maximum spatial positions in the network, + /// specifying a bounding box for the network geometry + aaboundingbox boundingbox(){ + + aaboundingbox bb; //create a bounding box + + //loop through every edge + for(unsigned e = 0; e < E.size(); e++){ + //loop through every point + for(unsigned p = 0; p < E[e].size(); p++) + bb.expand(E[e][p]); //expand the bounding box to include the point + } + + return bb; //return the bounding box + } + + ///render cylinder based on points from the top/bottom hat + ///@param C1 set of points from one of the hat + void renderCylinder(std::vector< stim::vec3 > C1, std::vector< stim::vec3 > C2, stim::vec3 P1, stim::vec3 P2) { + glBegin(GL_QUAD_STRIP); + for (unsigned i = 0; i < C1.size(); i++) { // for every point on the circle + stim::vec3 n1 = C1[i] - P1; // set normal vector for every vertex on the quads + stim::vec3 n2 = C2[i] - P2; + n1 = n1.norm(); + n2 = n2.norm(); + glNormal3f(n1[0], n1[1], n1[2]); + glVertex3f(C1[i][0], C1[i][1], C1[i][2]); + glNormal3f(n2[0], n2[1], n2[2]); + glVertex3f(C2[i][0], C2[i][1], C2[i][2]); + } + glEnd(); + } + + ///render the vertex as sphere based on glut build-in function + ///@param x, y, z are the three coordinates of the center point + ///@param radius is the radius of the sphere + ///@param subdivisions is the slice/stride along/around z-axis + void renderBall(T x, T y, T z, T radius, int subdivisions) { + glPushMatrix(); + glTranslatef(x, y, z); + glutSolidSphere(radius, subdivisions, subdivisions); + glPopMatrix(); + } + + ///render the vertex as sphere based on transformation + ///@param x, y, z are the three coordinates of the center point + ///@param radius is the radius of the sphere + ///@param slice is the number of subdivisions around the z-axis + ///@param stack is the number of subdivisions along the z-axis + void renderBall(T x, T y, T z, T radius, T slice, T stack) { + T step_z = stim::PI / slice; // step angle along z-axis + T step_xy = 2 * stim::PI / stack; // step angle in xy-plane + T xx[4], yy[4], zz[4]; // store coordinates + + T angle_z = 0.0; // start angle + T angle_xy = 0.0; + + glBegin(GL_QUADS); + for (unsigned i = 0; i < slice; i++) { // around the z-axis + angle_z = i * step_z; // step step_z each time + + for (unsigned j = 0; j < stack; j++) { // along the z-axis + angle_xy = j * step_xy; // step step_xy each time, draw floor by floor + + xx[0] = radius * std::sin(angle_z) * std::cos(angle_xy); // four vertices + yy[0] = radius * std::sin(angle_z) * std::sin(angle_xy); + zz[0] = radius * std::cos(angle_z); + + xx[1] = radius * std::sin(angle_z + step_z) * std::cos(angle_xy); + yy[1] = radius * std::sin(angle_z + step_z) * std::sin(angle_xy); + zz[1] = radius * std::cos(angle_z + step_z); + + xx[2] = radius * std::sin(angle_z + step_z) * std::cos(angle_xy + step_xy); + yy[2] = radius * std::sin(angle_z + step_z) * std::sin(angle_xy + step_xy); + zz[2] = radius * std::cos(angle_z + step_z); + + xx[3] = radius * std::sin(angle_z) * std::cos(angle_xy + step_xy); + yy[3] = radius * std::sin(angle_z) * std::sin(angle_xy + step_xy); + zz[3] = radius * std::cos(angle_z); + + for (unsigned k = 0; k < 4; k++) { + glVertex3f(x + xx[k], y + yy[k], z + zz[k]); // draw the floor plane + } + } + } + glEnd(); + } + + /// Render the network centerline as a series of line strips. + /// glCenterline0 is for only one input + void glCenterline0(){ + if (!glIsList(dlist)) { //if dlist isn't a display list, create it + dlist = glGenLists(1); //generate a display list + glNewList(dlist, GL_COMPILE); //start a new display list + for (unsigned e = 0; e < E.size(); e++) { //for each edge in the network + glBegin(GL_LINE_STRIP); + for (unsigned p = 0; p < E[e].size(); p++) { //for each point on that edge + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); //set the vertex position based on the current point + glTexCoord1f(0); //set white color + } + glEnd(); + } + glEndList(); //end the display list + } + glCallList(dlist); // render the display list + } + + ///render the network centerline as a series of line strips(when loading at least two networks, otherwise using glCenterline0()) + ///colors are based on metric values + void glCenterline(){ + + if(!glIsList(dlist)){ //if dlist isn't a display list, create it + dlist = glGenLists(1); //generate a display list + glNewList(dlist, GL_COMPILE); //start a new display list + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network + //unsigned errormag_id = E[e].nmags() - 1; + glBegin(GL_LINE_STRIP); + for(unsigned p = 0; p < E[e].size(); p++){ //for each point on that edge + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); //set the vertex position based on the current point + glTexCoord1f(E[e].r(p)); //set the texture coordinate based on the specified magnitude index + } + glEnd(); + } + glEndList(); //end the display list + } + glCallList(dlist); //render the display list + } + + ///render the network cylinder as a series of tubes(when only one network loaded) + void glCylinder0(T scale = 1.0f, bool undo = false) { + + float r1, r2; + if (undo == true) + glDeleteLists(dlist, 1); // delete display list + if (!glIsList(dlist)) { // if dlist isn't a display list, create it + dlist = glGenLists(1); // generate a display list + glNewList(dlist, GL_COMPILE); // start a new display list + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network + for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge + stim::circle C1 = E[e].circ(p - 1); + stim::circle C2 = E[e].circ(p); + r1 = E[e].r(p - 1); + r2 = E[e].r(p); + C1.set_R(scale * r1); // re-scale the circle to the same + C2.set_R(scale * r2); + std::vector< stim::vec3 > Cp1 = C1.glpoints(20); // get 20 points on the circle plane + std::vector< stim::vec3 > Cp2 = C2.glpoints(20); + glBegin(GL_QUAD_STRIP); + for (unsigned i = 0; i < Cp1.size(); i++) { + glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]); + glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]); + } + glEnd(); + } // set the texture coordinate based on the specified magnitude index + } + for (unsigned n = 0; n < V.size(); n++) { + for (unsigned i = 0; i < E.size(); i++) { + if (E[i].v[0] == n) { + r1 = E[i].r(0) * scale; + break; + } + else if (E[i].v[1] == n) { + r1 = E[i].r(E[i].size() - 1) * scale; + break; + } + } + renderBall(V[n][0], V[n][1], V[n][2], r1, 20); + } + glEndList(); // end the display list + } + glCallList(dlist); // render the display list + } + + ///render the network cylinder as a series of tubes + ///colors are based on metric values + void glCylinder(float sigma, float radius) { + + if (radius != sigma) // if render radius was changed by user, create a new display list + glDeleteLists(dlist, 1); + if (!glIsList(dlist)) { // if dlist isn't a display list, create it + dlist = glGenLists(1); // generate a display list + glNewList(dlist, GL_COMPILE); // start a new display list + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network + for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge + stim::circle C1 = E[e].circ(p - 1); + stim::circle C2 = E[e].circ(p); + C1.set_R(2*radius); // re-scale the circle to the same + C2.set_R(2*radius); + std::vector< stim::vec3 > Cp1 = C1.glpoints(20);// get 20 points on the circle plane + std::vector< stim::vec3 > Cp2 = C2.glpoints(20); + glBegin(GL_QUAD_STRIP); + for (unsigned i = 0; i < Cp1.size(); i++) { + stim::vec3 n1 = Cp1[i] - E[e][p - 1]; // set normal vector for every vertex on the quads + stim::vec3 n2 = Cp2[i] - E[e][p]; + n1 = n1.norm(); + n2 = n2.norm(); + + glNormal3f(n1[0], n1[1], n1[2]); + glTexCoord1f(E[e].r(p - 1)); + glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]); + glNormal3f(n2[0], n2[1], n2[2]); + glTexCoord1f(E[e].r(p)); + glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]); + } + glEnd(); + } // set the texture coordinate based on the specified magnitude index + } + for (unsigned n = 0; n < V.size(); n++) { + size_t num = V[n].e[0].size(); // store the number of outgoing edge of that vertex + if (num != 0) { // if it has outgoing edge + unsigned idx = V[n].e[0][0]; // find the index of first outgoing edge of that vertex + glTexCoord1f(E[idx].r(0)); // bind the texture as metric of first point on that edge + } + else { + unsigned idx = V[n].e[1][0]; // find the index of first incoming edge of that vertex + glTexCoord1f(E[idx].r(E[idx].size() - 1)); // bind the texture as metric of last point on that edge + } + renderBall(V[n][0], V[n][1], V[n][2], 2*radius, 20); + } + glEndList(); // end the display list + } + glCallList(dlist); // render the display list + } + + ///render a T as a adjoint network of GT in transparancy(darkgreen, overlaid) + void glAdjointCylinder(float sigma, float radius) { + + if (radius != sigma) // if render radius was changed by user, create a new display list + glDeleteLists(dlist + 4, 1); + if (!glIsList(dlist + 4)) { + glNewList(dlist + 4, GL_COMPILE); + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network + for (unsigned p = 1; p < E[e].size(); p++) { // for each point on that edge + stim::circle C1 = E[e].circ(p - 1); + stim::circle C2 = E[e].circ(p); + C1.set_R(2 * radius); // scale the circle to the same + C2.set_R(2 * radius); + std::vector< stim::vec3 >Cp1 = C1.glpoints(20); + std::vector< stim::vec3 >Cp2 = C2.glpoints(20); + glBegin(GL_QUAD_STRIP); + for (unsigned i = 0; i < Cp1.size(); i++) { // for every point on the circle(+1 means closing the circle) + glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]); + glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]); + } + glEnd(); + } // set the texture coordinate based on the specified magnitude index + } + for (unsigned n = 0; n < V.size(); n++) { + size_t num = V[n].e[0].size(); // store the number of outgoing edge of that vertex + if (num != 0) { // if it has outgoing edge + unsigned idx = V[n].e[0][0]; // find the index of first outgoing edge of that vertex + } + else { + unsigned idx = V[n].e[1][0]; // find the index of first incoming edge of that vertex + } + renderBall(V[n][0], V[n][1], V[n][2], 2 * radius, 20); + } + glEndList(); + } + glCallList(dlist + 4); + } + + ///render the network cylinder as series of tubes + ///@param I is a indicator: 0 -> GT, 1 -> T + ///@param map is the mapping relationship between two networks + ///@param colormap is the random generated color set for render + void glRandColorCylinder(int I, std::vector map, std::vector colormap, float sigma, float radius) { + + if (radius != sigma) // if render radius was changed by user, create a new display list + glDeleteLists(dlist + 2, 1); + if (!glIsList(dlist + 2)) { // if dlist isn't a display list, create it + glNewList(dlist + 2, GL_COMPILE); // start a new display list + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network + if (map[e] != unsigned(-1)) { + if (I == 0) { // if it is to render GT + glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]); + } + else { // if it is to render T + glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]); + } + + for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge + stim::circle C1 = E[e].circ(p - 1); + stim::circle C2 = E[e].circ(p); + C1.set_R(2*radius); // re-scale the circle to the same + C2.set_R(2*radius); + std::vector< stim::vec3 >Cp1 = C1.glpoints(20); + std::vector< stim::vec3 >Cp2 = C2.glpoints(20); + renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]); + } + } + else { + glColor3f(1.0, 1.0, 1.0); // white color for the un-mapping edges + for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge + stim::circle C1 = E[e].circ(p - 1); + stim::circle C2 = E[e].circ(p); + C1.set_R(2*radius); // scale the circle to the same + C2.set_R(2*radius); + std::vector< stim::vec3 >Cp1 = C1.glpoints(20); + std::vector< stim::vec3 >Cp2 = C2.glpoints(20); + renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]); + } + } + } + for (unsigned n = 0; n < V.size(); n++) { + size_t num_edge = V[n].e[0].size() + V[n].e[1].size(); + if (num_edge > 1) { // if it is the joint vertex + glColor3f(0.3, 0.3, 0.3); // gray + renderBall(V[n][0], V[n][1], V[n][2], 3*radius, 20); + } + else { // if it is the terminal vertex + glColor3f(0.6, 0.6, 0.6); // more white gray + renderBall(V[n][0], V[n][1], V[n][2], 3*radius, 20); + } + } + glEndList(); + } + glCallList(dlist + 2); + } + + ///render the network centerline as a series of line strips in random different color + ///@param I is a indicator: 0 -> GT, 1 -> T + ///@param map is the mapping relationship between two networks + ///@param colormap is the random generated color set for render + void glRandColorCenterline(int I, std::vector map, std::vector colormap) { + if (!glIsList(dlist + 2)) { + glNewList(dlist + 2, GL_COMPILE); + for (unsigned e = 0; e < E.size(); e++) { + if (map[e] != unsigned(-1)) { // if it has corresponding edge in another network + if (I == 0) // if it is to render GT + glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]); + else // if it is to render T + glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]); + + glBegin(GL_LINE_STRIP); + for (unsigned p = 0; p < E[e].size(); p++) { + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); + } + glEnd(); + } + else { + glColor3f(1.0, 1.0, 1.0); // white color for the un-mapping edges + glBegin(GL_LINE_STRIP); + for (unsigned p = 0; p < E[e].size(); p++) { + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); + } + glEnd(); + } + } + glEndList(); + } + glCallList(dlist + 2); + } + + void glAdjointCenterline() { + if (!glIsList(dlist + 4)) { + glNewList(dlist + 4, GL_COMPILE); + for (unsigned e = 0; e < E.size(); e++) { //for each edge in the network + + glBegin(GL_LINE_STRIP); + for (unsigned p = 0; p < E[e].size(); p++) { //for each point on that edge + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); //set the vertex position based on the current point + glTexCoord1f(E[e].r(p)); //set the texture coordinate based on the specified magnitude index + } + glEnd(); + } + glEndList(); + } + glCallList(dlist + 4); + } + + // highlight the difference part + void glDifferenceCylinder(int I, std::vector map, std::vector colormap, float sigma, float radius) { + + if (radius != sigma) // if render radius was changed by user, create a new display list + glDeleteLists(dlist + 6, 1); + if (!glIsList(dlist + 6)) { // if dlist isn't a display list, create it + glNewList(dlist + 6, GL_COMPILE); // start a new display list + for (unsigned e = 0; e < E.size(); e++) { // for each edge in the network + if (map[e] != unsigned(-1)) { + glEnable(GL_BLEND); //enable color blend + glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); //set blend function + glDisable(GL_DEPTH_TEST); //should disable depth + + if (I == 0) { // if it is to render GT + glColor4f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2], 0.1); + } + else { // if it is to render T + glColor4f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2], 0.1); + } + + for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge + stim::circle C1 = E[e].circ(p - 1); + stim::circle C2 = E[e].circ(p); + C1.set_R(2 * radius); // re-scale the circle to the same + C2.set_R(2 * radius); + std::vector< stim::vec3 >Cp1 = C1.glpoints(20); + std::vector< stim::vec3 >Cp2 = C2.glpoints(20); + renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]); + } + glDisable(GL_BLEND); + glEnable(GL_DEPTH_TEST); + } + else { + glColor3f(1.0, 1.0, 1.0); // white color for the un-mapping edges + for (unsigned p = 1; p < E[e].size(); p++) {// for each point on that edge + stim::circle C1 = E[e].circ(p - 1); + stim::circle C2 = E[e].circ(p); + C1.set_R(2 * radius); // scale the circle to the same + C2.set_R(2 * radius); + std::vector< stim::vec3 >Cp1 = C1.glpoints(20); + std::vector< stim::vec3 >Cp2 = C2.glpoints(20); + renderCylinder(Cp1, Cp2, E[e][p - 1], E[e][p]); + } + } + } + for (unsigned n = 0; n < V.size(); n++) { + + size_t num_edge = V[n].e[0].size() + V[n].e[1].size(); + if (num_edge > 1) { // if it is the joint vertex + glColor4f(0.3, 0.3, 0.3, 0.1); // gray + renderBall(V[n][0], V[n][1], V[n][2], 3 * radius, 20); + } + else { // if it is the terminal vertex + glColor4f(0.6, 0.6, 0.6, 0.1); // more white gray + renderBall(V[n][0], V[n][1], V[n][2], 3 * radius, 20); + } + } + glEndList(); + } + glCallList(dlist + 6); + } + + //void glRandColorCenterlineGT(GLuint &dlist1, std::vector map, std::vector colormap){ + // if(!glIsList(dlist1)){ + // dlist1 = glGenLists(1); + // glNewList(dlist1, GL_COMPILE); + // for(unsigned e = 0; e < E.size(); e++){ + // if(map[e] != unsigned(-1)){ + // glColor3f(colormap[e * 3 + 0], colormap[e * 3 + 1], colormap[e * 3 + 2]); + // glBegin(GL_LINE_STRIP); + // for(unsigned p = 0; p < E[e].size(); p++){ + // glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); + // } + // glEnd(); + // for (unsigned p = 0; p < E[e].size() - 1; p++) { + // 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); + // } + // } + // else{ + // glColor3f(1.0, 1.0, 1.0); + // glBegin(GL_LINE_STRIP); + // for(unsigned p = 0; p < E[e].size(); p++){ + // glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); + // } + // glEnd(); + // } + // } + // for (unsigned v = 0; v < V.size(); v++) { + // size_t num_edge = V[v].e[0].size() + V[v].e[1].size(); + // if (num_edge > 1) { + // glColor3f(0.3, 0.3, 0.3); // gray color for vertex + // renderBall(V[v][0], V[v][1], V[v][2], 20, 20); + // } + // } + // glEndList(); + // } + // glCallList(dlist1); + //} + + //void glRandColorCenterlineT(GLuint &dlist2, std::vector map, std::vector colormap){ + // if(!glIsList(dlist2)){ + // dlist2 = glGenLists(1); + // glNewList(dlist2, GL_COMPILE); + // for(unsigned e = 0; e < E.size(); e++){ + // if(map[e] != unsigned(-1)){ + // glColor3f(colormap[map[e] * 3 + 0], colormap[map[e] * 3 + 1], colormap[map[e] * 3 + 2]); + // glBegin(GL_LINE_STRIP); + // for(unsigned p = 0; p < E[e].size(); p++){ + // glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); + // } + // glEnd(); + // for (unsigned p = 0; p < E[e].size() - 1; p++) { + // 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); + // } + // } + // else{ + // glColor3f(1.0, 1.0, 1.0); + // glBegin(GL_LINE_STRIP); + // for(unsigned p = 0; p < E[e].size(); p++){ + // glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); + // } + // glEnd(); + // } + // } + // for (unsigned v = 0; v < V.size(); v++) { + // size_t num_edge = V[v].e[0].size() + V[v].e[1].size(); + // if (num_edge > 1) { + // glColor3f(0.3, 0.3, 0.3); // gray color for vertex + // renderBall(V[v][0], V[v][1], V[v][2], 20, 20); + // } + // } + // glEndList(); + // } + // glCallList(dlist2); + //} + + + //void renderCylinder(T x1, T y1, T z1, T x2, T y2, T z2, T radius, int subdivisions) { + // T dx = x2 - x1; + // T dy = y2 - y1; + // T dz = z2 - z1; + // /// handle the degenerate case with an approximation + // if (dz == 0) + // dz = .00000001; + // T d = sqrt(dx*dx + dy*dy + dz*dz); + // T ax = 57.2957795*acos(dz / d); // 180°/pi + // if (dz < 0.0) + // ax = -ax; + // T rx = -dy*dz; + // T ry = dx*dz; + + // glPushMatrix(); + // glTranslatef(x1, y1, z1); + // glRotatef(ax, rx, ry, 0.0); + + // glutSolidCylinder(radius, d, subdivisions, 1); + // glPopMatrix(); + //} + + + /// render the network centerline from swc file as a series of strips in different colors based on the neuronal type + /// glCenterline0_swc is for only one input + /*void glCenterline0_swc() { + if (!glIsList(dlist)) { // if dlist isn't a display list, create it + dlist = glGenLists(1); // generate a display list + glNewList(dlist, GL_COMPILE); // start a new display list + for (unsigned e = 0; e < E.size(); e++) { + int type = NT[e]; // get the neuronal type + switch (type) { + case 0: + glColor3f(1.0f, 1.0f, 1.0f); // white for undefined + glBegin(GL_LINE_STRIP); + for (unsigned p = 0; p < E[e].size(); p++) { + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); + } + glEnd(); + break; + case 1: + glColor3f(1.0f, 0.0f, 0.0f); // red for soma + glBegin(GL_LINE_STRIP); + for (unsigned p = 0; p < E[e].size(); p++) { + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); + } + glEnd(); + break; + case 2: + glColor3f(1.0f, 0.5f, 0.0f); // orange for axon + glBegin(GL_LINE_STRIP); + for (unsigned p = 0; p < E[e].size(); p++) { + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); + } + glEnd(); + break; + case 3: + glColor3f(1.0f, 1.0f, 0.0f); // yellow for undefined + glBegin(GL_LINE_STRIP); + for (unsigned p = 0; p < E[e].size(); p++) { + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); + } + glEnd(); + break; + case 4: + glColor3f(0.0f, 1.0f, 0.0f); // green for undefined + glBegin(GL_LINE_STRIP); + for (unsigned p = 0; p < E[e].size(); p++) { + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); + } + glEnd(); + break; + case 5: + glColor3f(0.0f, 1.0f, 1.0f); // verdant for undefined + glBegin(GL_LINE_STRIP); + for (unsigned p = 0; p < E[e].size(); p++) { + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); + } + glEnd(); + break; + case 6: + glColor3f(0.0f, 0.0f, 1.0f); // blue for undefined + glBegin(GL_LINE_STRIP); + for (unsigned p = 0; p < E[e].size(); p++) { + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); + } + glEnd(); + break; + case 7: + glColor3f(0.5f, 0.0f, 1.0f); // purple for undefined + glBegin(GL_LINE_STRIP); + for (unsigned p = 0; p < E[e].size(); p++) { + glVertex3f(E[e][p][0], E[e][p][1], E[e][p][2]); + } + glEnd(); + break; + } + } + glEndList(); //end the display list + } + glCallList(dlist); // render the display list + }*/ + + ///render the network cylinder as a series of tubes + ///colors are based on metric values + //void glCylinder(float sigma) { + // if (!glIsList(dlist)) { //if dlist isn't a display list, create it + // dlist = glGenLists(1); //generate a display list + // glNewList(dlist, GL_COMPILE); //start a new display list + // for (unsigned e = 0; e < E.size(); e++) { //for each edge in the network + // for (unsigned p = 1; p < E[e].size() - 1; p++) { // for each point on that edge + // stim::circle C1 = E[e].circ(p - 1); + // stim::circle C2 = E[e].circ(p); + // C1.set_R(2.5*sigma); // scale the circle to the same + // C2.set_R(2.5*sigma); + // std::vector< stim::vec3 >Cp1 = C1.glpoints(20); + // std::vector< stim::vec3 >Cp2 = C2.glpoints(20); + // glBegin(GL_QUAD_STRIP); + // for (unsigned i = 0; i < Cp1.size(); i++) { // for every point on the circle(+1 means closing the circle) + // glVertex3f(Cp1[i][0], Cp1[i][1], Cp1[i][2]); + // glVertex3f(Cp2[i][0], Cp2[i][1], Cp2[i][2]); + // glTexCoord1f(E[e].r(p)); + // } + // glEnd(); + // } //set the texture coordinate based on the specified magnitude index + // } + // for (unsigned n = 0; n < V.size(); n++) { + // size_t num = V[n].e[0].size(); //store the number of outgoing edge of that vertex + // if (num != 0) { //if it has outgoing edge + // unsigned idx = V[n].e[0][0]; //find the index of first outgoing edge of that vertex + // glTexCoord1f(E[idx].r(0)); //bind the texture as metric of first point on that edge + // } + // else { + // unsigned idx = V[n].e[1][0]; //find the index of first incoming edge of that vertex + // glTexCoord1f(E[idx].r(E[idx].size() - 1)); //bind the texture as metric of last point on that edge + // } + // renderBall(V[n][0], V[n][1], V[n][2], 2.5*sigma, 20); + // } + // glEndList(); //end the display list + // } + // glCallList(dlist); //render the display list + //} + +}; //end stim::gl_network class +}; //end stim namespace + + + +#endif \ No newline at end of file diff --git a/main.cu b/main.cu index 3d2292b..1ed1314 100644 --- a/main.cu +++ b/main.cu @@ -7,8 +7,8 @@ #include #include #include -#include -#include +#include "gl_network.h" +#include "network.h" #include // OpenGL includes @@ -929,7 +929,7 @@ int main(int argc, char* argv[]) args.add("help", "prints this help"); args.add("sigma", "force a sigma value to specify the tolerance of the network comparison", "3"); args.add("gui", "display the network or network comparison using OpenGL"); - args.add("device", "choose specific device to run", "0"); + args.add("device", "choose specific device to run (use -1 for CPU only)", "0"); args.add("features", "save features to a CSV file, specify file name"); args.add("mapping", "mapping input according to similarity"); args.add("stack", "load the image stacks"); diff --git a/network.h b/network.h new file mode 100644 index 0000000..ac84601 --- /dev/null +++ b/network.h @@ -0,0 +1,1520 @@ +#ifndef STIM_NETWORK_H +#define STIM_NETWORK_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "cylinder.h" +#include +#include +#include +//********************help function******************** +// gaussian_function +CUDA_CALLABLE float gaussianFunction(float x, float std = 25) { return exp(-x / (2 * std*std)); } // std default sigma value is 25 + +// compute metric in parallel +#ifdef __CUDACC__ +template +__global__ void find_metric_parallel(T* M, size_t n, T* D, float sigma){ + size_t x = blockDim.x * blockIdx.x + threadIdx.x; + if(x >= n) return; + M[x] = 1.0f - gaussianFunction(D[x], sigma); +} + +//find the corresponding edge index from array index +__global__ void find_edge_index_parallel(size_t* I, size_t n, unsigned* R, size_t* E, size_t ne){ + size_t x = blockDim.x * blockIdx.x + threadIdx.x; + if(x >= n) return; + unsigned i = 0; + size_t N = 0; + for(unsigned e = 0; e < ne; e++){ + N += E[e]; + if(I[x] < N){ + R[x] = i; + break; + } + i++; + } +} +#endif + +//hard-coded factor +int threshold_fac; + +namespace stim{ +/** This is the a class that interfaces with gl_spider in order to store the currently + * segmented network. The following data is stored and can be extracted: + * 1)Network geometry and centerline. + * 2)Network connectivity (a graph of nodes and edges), reconstructed using kdtree. +*/ + +template +class network{ + + ///Each edge is a fiber with two nodes. + ///Each node is an in index to the endpoint of the fiber in the nodes array. + class edge : public cylinder + { + public: + + unsigned int v[2]; //unique id's designating the starting and ending + // default constructor + edge() : cylinder() { + v[1] = (unsigned)(-1); v[0] = (unsigned)(-1); + } + /// Constructor - creates an edge from a list of points by calling the stim::fiber constructor +/* + ///@param v0, the starting index. + ///@param v1, the ending index. + ///@param sz, the number of point in the fiber. + edge(unsigned int v0, unsigned int v1, unsigned int sz) : cylinder( + { + + } +*/ + edge(std::vector > p, std::vector s) + : cylinder(p,s) + { + } + ///@param p is an array of positions in space + edge(stim::centerline p) : cylinder(p){} + + /// Copy constructor creates an edge from a cylinder + edge(stim::cylinder f) : cylinder(f) {} + + /// Resamples an edge by calling the fiber resampling function + edge resample(T spacing){ + edge e(cylinder::resample(spacing)); //call the fiber->edge constructor + e.v[0] = v[0]; //copy the vertex data + e.v[1] = v[1]; + + return e; //return the new edge + } + + /// Output the edge information as a string + std::string str(){ + std::stringstream ss; + ss<<"("<::size()<<")\tl = "<length()<<"\t"< split(unsigned int idx){ + + std::vector< stim::cylinder > C; + C.resize(2); + C = (*this).cylinder::split(idx); + std::vector E(C.size()); + + for(unsigned e = 0; e < E.size(); e++){ + E[e] = C[e]; + } + return E; + } + + /// operator for writing the edge information into a binary .nwt file. + friend std::ofstream& operator<<(std::ofstream& out, const edge& e) + { + out.write(reinterpret_cast(&e.v[0]), sizeof(unsigned int)); ///write the starting point. + out.write(reinterpret_cast(&e.v[1]), sizeof(unsigned int)); ///write the ending point. + unsigned int sz = e.size(); ///write the number of point in the edge. + out.write(reinterpret_cast(&sz), sizeof(unsigned int)); + for(int i = 0; i < sz; i++) ///write each point + { + stim::vec3 point = e[i]; + out.write(reinterpret_cast(&point[0]), 3*sizeof(T)); + // for(int j = 0; j < nmags(); j++) //future code for multiple mags + // { + out.write(reinterpret_cast(&e.R[i]), sizeof(T)); ///write the radius + //std::cout << point.str() << " " << e.R[i] << std::endl; + // } + } + return out; //return stream + } + + /// operator for reading an edge from a binary .nwt file. + friend std::ifstream& operator>>(std::ifstream& in, edge& e) + { + unsigned int v0, v1, sz; + in.read(reinterpret_cast(&v0), sizeof(unsigned int)); //read the staring point. + in.read(reinterpret_cast(&v1), sizeof(unsigned int)); //read the ending point + in.read(reinterpret_cast(&sz), sizeof(unsigned int)); //read the number of points in the edge +// stim::centerline temp = stim::centerline(sz); //allocate the new edge +// e = edge(temp); + std::vector > p(sz); + std::vector r(sz); + for(int i = 0; i < sz; i++) //set the points and radii to the newly read values + { + stim::vec3 point; + in.read(reinterpret_cast(&point[0]), 3*sizeof(T)); + p[i] = point; + T mag; + // for(int j = 0; j < nmags(); j++) ///future code for mags + // { + in.read(reinterpret_cast(&mag), sizeof(T)); + r[i] = mag; + //std::cout << point.str() << " " << mag << std::endl; + // } + } + e = edge(p,r); + e.v[0] = v0; e.v[1] = v1; + return in; + } + }; + + ///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. + class vertex : public stim::vec3 + { + public: + //std::vector edges; //indices of edges connected to this node. + std::vector e[2]; //indices of edges going out (e[0]) and coming in (e[1]) + //stim::vec3 p; //position of this node in physical space. + //default constructor + vertex() : stim::vec3() + { + } + //constructor takes a stim::vec + vertex(stim::vec3 p) : stim::vec3(p){} + + /// Output the vertex information as a string + std::string + str(){ + std::stringstream ss; + ss<<"\t(x, y, z) = "<::str(); + + if(e[0].size() > 0){ + ss<<"\t> "; + for(unsigned int o = 0; o < e[0].size(); o++) + ss< 0){ + ss<<"\t< "; + for(unsigned int i = 0; i < e[1].size(); i++) + ss<(&v.ptr[0]), 3*sizeof(T)); ///write physical vertex location + out.write(reinterpret_cast(&s0), sizeof(unsigned int)); ///write the number of "outgoing edges" + out.write(reinterpret_cast(&s1), sizeof(unsigned int)); ///write the number of "incoming edges" + if (s0 != 0) + out.write(reinterpret_cast(&v.e[0][0]), sizeof(unsigned int)*v.e[0].size()); ///write the "outgoing edges" + if (s1 != 0) + out.write(reinterpret_cast(&v.e[1][0]), sizeof(unsigned int)*v.e[1].size()); ///write the "incoming edges" + return out; + } + + ///operator for reading the vector out of the stream; + friend std::ifstream& operator>>(std::ifstream& in, vertex& v) + { + in.read(reinterpret_cast(&v[0]), 3*sizeof(T)); ///read the physical position + unsigned int s[2]; + in.read(reinterpret_cast(&s[0]), 2*sizeof(unsigned int)); ///read the sizes of incoming and outgoing edge arrays + + std::vector one(s[0]); + std::vector two(s[1]); + v.e[0] = one; + v.e[1] = two; + if (one.size() != 0) + in.read(reinterpret_cast(&v.e[0][0]), s[0] * sizeof(unsigned int)); ///read the arrays of "outgoing edges" + if (two.size() != 0) + in.read(reinterpret_cast(&v.e[1][0]), s[1] * sizeof(unsigned int)); ///read the arrays of "incoming edges" + return in; + } + + }; + +protected: + + std::vector E; //list of edges + std::vector V; //list of vertices. + +public: + + ///default constructor + network() + { + + } + + ///constructor with a file to load. + network(std::string fileLocation) + { + load_obj(fileLocation); + } + + ///Returns the number of edges in the network. + unsigned int edges(){ + return E.size(); + } + + ///Returns the number of nodes in the network. + unsigned int vertices(){ + return V.size(); + } + + ///Returns the radius at specific point in the edge + T get_r(unsigned e, unsigned i) { + return E[e].r(i); + } + + ///Returns the average radius of specific edge + T get_average_r(unsigned e) { + T result = 0.0; + unsigned n = E[e].size(); + for (unsigned p = 0; p < n; p++) + result += E[e].r(p); + + return (T)result / n; + } + + ///Returns the length of current edge + T get_l(unsigned e) { + return E[e].length(); + } + + ///Returns the start vertex of current edge + size_t get_start_vertex(unsigned e) { + return E[e].v[0]; + } + + ///Returns the end vertex of current edge + size_t get_end_vertex(unsigned e) { + return E[e].v[1]; + } + + ///Returns one vertex + stim::vec3 get_vertex(unsigned i) { + return V[i]; + } + + ///Returns the boundary vertices' indices + std::vector get_boundary_vertex() { + std::vector result; + + for (unsigned v = 0; v < V.size(); v++) { + if (V[v].e[0].size() + V[v].e[1].size() == 1) { // boundary vertex + result.push_back(v); + } + } + + return result; + } + + ///Set radius + void set_r(unsigned e, std::vector radius) { + E[e].cylinder::copy_r(radius); + } + + void set_r(unsigned e, T radius) { + for (size_t i = 0; i < E[e].size(); i++) + E[e].cylinder::set_r(i, radius); + } + //scale the network by some constant value + // I don't think these work?????? + /*std::vector operator*(T s){ + for (unsigned i=0; i< vertices; i ++ ){ + V[i] = V[i] * s; + } + return V; + } + + std::vector operator*(vec s){ + for (unsigned i=0; i< vertices; i ++ ){ + for (unsigned dim = 0 ; dim< 3; dim ++){ + V[i][dim] = V[i][dim] * s[dim]; + } + } + return V; + }*/ + + // Returns an average of branching index in the network + double BranchingIndex(){ + double B=0; + for(unsigned v=0; v < V.size(); v ++){ + B += ((V[v].e[0].size()) + (V[v].e[1].size())); + } + B = B / V.size(); + return B; + + } + + // Returns number of branch points in thenetwork + unsigned int BranchP(){ + unsigned int B=0; + unsigned int c; + for(unsigned v=0; v < V.size(); v ++){ + c = ((V[v].e[0].size()) + (V[v].e[1].size())); + if (c > 2){ + B += 1;} + } + return B; + + } + + // Returns number of end points (tips) in thenetwork + unsigned int EndP(){ + unsigned int B=0; + unsigned int c; + for(unsigned v=0; v < V.size(); v ++){ + c = ((V[v].e[0].size()) + (V[v].e[1].size())); + if (c == 1){ + B += 1;} + } + return B; + + } + + //// Returns a dictionary with the key as the vertex + //std::map,unsigned int> DegreeDict(){ + // std::map,unsigned int> dd; + // unsigned int c = 0; + // for(unsigned v=0; v < V.size(); v ++){ + // c = ((V[v].e[0].size()) + (V[v].e[1].size())); + // dd[V[v]] = c; + // } + // return dd; + //} + + //// Return number of branching stems + //unsigned int Stems(){ + // unsigned int s = 0; + // std::map,unsigned int> dd; + // dd = DegreeDict(); + // //for(unsigned v=0; v < V.size(); v ++){ + // // V[v].e[0]. + // return s; + //} + + //Calculate Metrics--------------------------------------------------- + // Returns an average of fiber/edge lengths in the network + double Lengths(){ + stim::vec L; + double sumLength = 0; + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network + L.push_back(E[e].length()); //append the edge length + sumLength = sumLength + E[e].length(); + } + double avg = sumLength / E.size(); + return avg; + } + + + // Returns an average of tortuosities in the network + double Tortuosities(){ + stim::vec t; + stim::vec id1, id2; // starting and ending vertices of the edge + double distance;double tortuosity;double sumTortuosity = 0; + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network + id1 = E[e][0]; //get the edge starting point + id2 = E[e][E[e].size() - 1]; //get the edge ending point + distance = (id1 - id2).len(); //displacement between the starting and ending points + if(distance > 0){ + tortuosity = E[e].length()/ distance ; // tortuoisty = edge length / edge displacement + } + else{ + tortuosity = 0;} + t.push_back(tortuosity); + sumTortuosity += tortuosity; + } + double avg = sumTortuosity / E.size(); + return avg; + } + + // Returns average contraction of the network + double Contractions(){ + stim::vec t; + stim::vec id1, id2; // starting and ending vertices of the edge + double distance;double contraction;double sumContraction = 0; + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network + id1 = E[e][0]; //get the edge starting point + id2 = E[e][E[e].size() - 1]; //get the edge ending point + distance = (id1 - id2).len(); //displacement between the starting and ending points + contraction = distance / E[e].length(); // tortuoisty = edge length / edge displacement + t.push_back(contraction); + sumContraction += contraction; + } + double avg = sumContraction / E.size(); + return avg; + } + + // returns average fractal dimension of the branches of the network + double FractalDimensions(){ + stim::vec t; + stim::vec id1, id2; // starting and ending vertices of the edge + double distance;double fract;double sumFractDim = 0; + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network + id1 = E[e][0]; //get the edge starting point + id2 = E[e][E[e].size() - 1]; //get the edge ending point + distance = (id1 - id2).len(); //displacement between the starting and ending points + fract = std::log(distance) / std::log(E[e].length()); // tortuoisty = edge length / edge displacement + t.push_back(sumFractDim); + sumFractDim += fract; + } + double avg = sumFractDim / E.size(); + return avg; + } + + //returns a cylinder represented a given fiber (based on edge index) + stim::cylinder get_cylinder(unsigned e){ + return E[e]; //return the specified edge (casting it to a fiber) + } + + /// subdivide current network + void subdivision() { + + std::vector ori_index; // original index + std::vector new_index; // new index + std::vector nE; // new edge + std::vector nV; // new vector + unsigned id = 0; + unsigned num_edge = (*this).E.size(); + + for (unsigned i = 0; i < num_edge; i++) { + if (E[i].size() == 2) { // if current edge can't be subdivided + stim::centerline line(2); + for (unsigned k = 0; k < 2; k++) + line[k] = E[i][k]; + line.update(); + + edge new_edge(line); + + vertex new_vertex = new_edge[0]; + id = E[i].v[0]; + auto position = std::find(ori_index.begin(), ori_index.end(), id); + if (position == ori_index.end()) { // new vertex + ori_index.push_back(id); + new_index.push_back(nV.size()); + + new_vertex.e[0].push_back(nE.size()); + new_edge.v[0] = nV.size(); + nV.push_back(new_vertex); // push back vertex as a new vertex + } + else { // existing vertex + int k = std::distance(ori_index.begin(), position); + new_edge.v[0] = new_index[k]; + nV[new_index[k]].e[0].push_back(nE.size()); + } + + new_vertex = new_edge[1]; + id = E[i].v[1]; + position = std::find(ori_index.begin(), ori_index.end(), id); + if (position == ori_index.end()) { // new vertex + ori_index.push_back(id); + new_index.push_back(nV.size()); + + new_vertex.e[1].push_back(nE.size()); + new_edge.v[1] = nV.size(); + nV.push_back(new_vertex); // push back vertex as a new vertex + } + else { // existing vertex + int k = std::distance(ori_index.begin(), position); + new_edge.v[1] = new_index[k]; + nV[new_index[k]].e[1].push_back(nE.size()); + } + + nE.push_back(new_edge); + + nE[nE.size() - 1].cylinder::set_r(0, E[i].cylinder::r(0)); + nE[nE.size() - 1].cylinder::set_r(1, E[i].cylinder::r(1)); + } + else { // subdivide current edge + for (unsigned j = 0; j < E[i].size() - 1; j++) { + stim::centerline line(2); + for (unsigned k = 0; k < 2; k++) + line[k] = E[i][j + k]; + line.update(); + + edge new_edge(line); + + if (j == 0) { // edge contains original starting point + vertex new_vertex = new_edge[0]; + id = E[i].v[0]; + auto position = std::find(ori_index.begin(), ori_index.end(), id); + if (position == ori_index.end()) { // new vertex + ori_index.push_back(id); + new_index.push_back(nV.size()); + + new_vertex.e[0].push_back(nE.size()); + new_edge.v[0] = nV.size(); + nV.push_back(new_vertex); // push back vertex as a new vertex + } + else { // existing vertex + int k = std::distance(ori_index.begin(), position); + new_edge.v[0] = new_index[k]; + nV[new_index[k]].e[0].push_back(nE.size()); + } + + new_vertex = new_edge[1]; + new_vertex.e[1].push_back(nE.size()); + new_edge.v[1] = nV.size(); + nV.push_back(new_vertex); // push back internal point as a new vertex + + nE.push_back(new_edge); + } + + else if (j == E[i].size() - 2) { // edge contains original ending point + + vertex new_vertex = new_edge[1]; + nV[nV.size() - 1].e[0].push_back(nE.size()); + new_edge.v[0] = nV.size() - 1; + + id = E[i].v[1]; + auto position = std::find(ori_index.begin(), ori_index.end(), id); + if (position == ori_index.end()) { // new vertex + ori_index.push_back(id); + new_index.push_back(nV.size()); + + new_vertex.e[1].push_back(nE.size()); + new_edge.v[1] = nV.size(); + nV.push_back(new_vertex); // push back vertex as a new vertex + } + else { // existing vertex + int k = std::distance(ori_index.begin(), position); + new_edge.v[1] = new_index[k]; + nV[new_index[k]].e[1].push_back(nE.size()); + } + + nE.push_back(new_edge); + } + + else { + vertex new_vertex = new_edge[1]; + + nV[nV.size() - 1].e[0].push_back(nE.size()); + new_vertex.e[1].push_back(nE.size()); + new_edge.v[0] = nV.size() - 1; + new_edge.v[1] = nV.size(); + nV.push_back(new_vertex); + + nE.push_back(new_edge); + } + + // get radii + nE[nE.size() - 1].cylinder::set_r(0, E[i].cylinder::r(j)); + nE[nE.size() - 1].cylinder::set_r(1, E[i].cylinder::r(j + 1)); + } + } + } + + (*this).E = nE; + (*this).V = nV; + } + + //load a network from an OBJ file + void load_obj(std::string filename){ + + stim::obj O; //create an OBJ object + O.load(filename); //load the OBJ file as an object + + std::vector id2vert; //this list stores the OBJ vertex ID associated with each network vertex + + unsigned i[2]; //temporary, IDs associated with the first and last points in an OBJ line + + //for each line in the OBJ object + for(unsigned int l = 1; l <= O.numL(); l++){ + + std::vector< stim::vec > c; //allocate an array of points for the vessel centerline + O.getLine(l, c); //get the fiber centerline + + stim::centerline c3(c.size()); + for(size_t j = 0; j < c.size(); j++) + c3[j] = c[j]; + c3.update(); + + // edge new_edge = c3; ///This is dangerous. + edge new_edge(c3); + + //create an edge from the given centerline + unsigned int I = new_edge.size(); //calculate the number of points on the centerline + + //get the first and last vertex IDs for the line + std::vector< unsigned > id; //create an array to store the centerline point IDs + O.getLinei(l, id); //get the list of point IDs for the line + i[0] = id.front(); //get the OBJ ID for the first element of the line + i[1] = id.back(); //get the OBJ ID for the last element of the line + + std::vector::iterator it; //create an iterator for searching the id2vert array + unsigned it_idx; //create an integer for the id2vert entry index + + //find out if the nodes for this fiber have already been created + it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node + if(it == id2vert.end()){ //if i[0] hasn't already been used + vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position + bool flag = false; + unsigned j = 0; + for (; j < V.size(); j++) { // check whether current vertex is already exist + if (new_vertex == V[j]) { + flag = true; + break; + } + } + if (!flag) { // unique one + new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing + new_edge.v[0] = V.size(); //add the new edge to the edge + V.push_back(new_vertex); //add the new vertex to the vertex list + id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list + } + else { + V[j].e[0].push_back(E.size()); + new_edge.v[0] = j; + } + } + else{ //if the vertex already exists + it_idx = std::distance(id2vert.begin(), it); + V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing + new_edge.v[0] = it_idx; + } + + it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID + if(it == id2vert.end()){ //if i[1] hasn't already been used + vertex new_vertex = new_edge[I-1]; //create a new vertex, assign it a position + bool flag = false; + unsigned j = 0; + for (; j < V.size(); j++) { // check whether current vertex is already exist + if (new_vertex == V[j]) { + flag = true; + break; + } + } + if (!flag) { + new_vertex.e[1].push_back(E.size()); //add the current edge as incoming + new_edge.v[1] = V.size(); //add the new vertex to the edge + V.push_back(new_vertex); //add the new vertex to the vertex list + id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list + } + else { + V[j].e[1].push_back(E.size()); + new_edge.v[1] = j; + } + } + else{ //if the vertex already exists + it_idx = std::distance(id2vert.begin(), it); + V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming + new_edge.v[1] = it_idx; + } + + E.push_back(new_edge); //push the edge to the list + + } + + // copy the radii information from OBJ + /*if (O.numVT()) { + unsigned k = 0; + for (unsigned i = 0; i < E.size(); i++) { + for (unsigned j = 0; j < E[i].size(); j++) { + E[i].cylinder::set_r(j, O.getVT(k)[0] / 2); + k++; + } + } + }*/ + // OBJ class assumes that in L the two values are equal + if (O.numVT()) { + std::vector< unsigned > id; //create an array to store the centerline point IDs + for (unsigned i = 0; i < O.numL(); i++) { + id.clear(); + O.getLinei(i + 1, id); //get the list of point IDs for the line + for (unsigned j = 0; j < id.size(); j++) + E[i].cylinder::set_r(j, O.getVT(id[j] - 1)[0] / 2); + } + } + } + + ///loads a .nwt file. Reads the header and loads the data into the network according to the header. + void + loadNwt(std::string filename) + { + int dims[2]; ///number of vertex, number of edges + readHeader(filename, &dims[0]); //read header + std::ifstream file; + file.open(filename.c_str(), std::ios::in | std::ios::binary); ///skip header information. + file.seekg(14+58+4+4, file.beg); + vertex v; + for(int i = 0; i < dims[0]; i++) ///for every vertex, read vertex, add to network. + { + file >> v; + V.push_back(v); +// std::cout << i << " " << v.str() << std::endl; + } + + std::cout << std::endl; + for(int i = 0; i < dims[1]; i++) ///for every edge, read edge, add to network. + { + edge e; + file >> e; + E.push_back(e); + //std::cout << i << " " << E[i].str() << std::endl; // not necessary? + } + file.close(); + } + + ///saves a .nwt file. Writes the header in raw text format, then saves the network as a binary file. + void + saveNwt(std::string filename) + { + writeHeader(filename); + std::ofstream file; + file.open(filename.c_str(), std::ios::out | std::ios::binary | std::ios::app); ///since we have written the header we are not appending. + for(int i = 0; i < V.size(); i++) ///look through the Vertices and write each one. + { +// std::cout << i << " " << V[i].str() << std::endl; + file << V[i]; + } + for(int i = 0; i < E.size(); i++) ///loop through the Edges and write each one. + { + //std::cout << i << " " << E[i].str() << std::endl; // not necesarry? + file << E[i]; + } + file.close(); + } + + + ///Writes the header information to a .nwt file. + void + writeHeader(std::string filename) + { + std::string magicString = "nwtFileFormat "; ///identifier for the file. + std::string desc = "fileid(14B), desc(58B), #vertices(4B), #edges(4B): bindata"; + int hNumVertices = V.size(); ///int byte header storing the number of vertices in the file + int hNumEdges = E.size(); ///int byte header storing the number of edges. + std::ofstream file; + file.open(filename.c_str(), std::ios::out | std::ios::binary); + std::cout << hNumVertices << " " << hNumEdges << std::endl; + file.write(reinterpret_cast(&magicString.c_str()[0]), 14); //write the file id + file.write(reinterpret_cast(&desc.c_str()[0]), 58); //write the description + file.write(reinterpret_cast(&hNumVertices), sizeof(int)); //write #vert. + file.write(reinterpret_cast(&hNumEdges), sizeof(int)); //write #edges +// file << magicString.c_str() << desc.c_str() << hNumVertices << hNumEdges; + file.close(); + + } + + ///Reads the header information from a .nwt file. + void + readHeader(std::string filename, int *dims) + { + char magicString[14]; ///id + char desc[58]; ///description + int hNumVertices; ///#vert + int hNumEdges; ///#edges + std::ifstream file; ////create stream + file.open(filename.c_str(), std::ios::in | std::ios::binary); + file.read(reinterpret_cast(&magicString[0]), 14); ///read the file id. + file.read(reinterpret_cast(&desc[0]), 58); ///read the description + file.read(reinterpret_cast(&hNumVertices), sizeof(int)); ///read the number of vertices + file.read(reinterpret_cast(&hNumEdges), sizeof(int)); ///read the number of edges +// std::cout << magicString << desc << hNumVertices << " " << hNumEdges << std::endl; + file.close(); ///close the file. + dims[0] = hNumVertices; ///fill the returned reference. + dims[1] = hNumEdges; + } + + //load a network from an SWC file + void load_swc(std::string filename) { + stim::swc S; // create swc variable + S.load(filename); // load the node information + S.create_tree(); // link those node according to their linking relationships as a tree + S.resample(); + + //NT.push_back(S.node[0].type); // set the neuronal_type value to the first vertex in the network + std::vector id2vert; // this list stores the SWC vertex ID associated with each network vertex + unsigned i[2]; // temporary, IDs associated with the first and last points + + for (unsigned int l = 0; l < S.numE(); l++) { // for every edge + //NT.push_back(S.node[l].type); + + std::vector< stim::vec3 > c; + S.get_points(l, c); + + stim::centerline c3(c.size()); // new fiber + + for (unsigned j = 0; j < c.size(); j++) + c3[j] = c[j]; // copy the points + + c3.update(); // update the L information + + stim::cylinder C3(c3); // create a new cylinder in order to copy the origin radius information + // upadate the R information + std::vector radius; + S.get_radius(l, radius); + + C3.copy_r(radius); + + edge new_edge(C3); // new edge + + //create an edge from the given centerline + unsigned int I = (unsigned)new_edge.size(); //calculate the number of points on the centerline + + //get the first and last vertex IDs for the line + i[0] = S.E[l].front(); + i[1] = S.E[l].back(); + + std::vector::iterator it; //create an iterator for searching the id2vert array + unsigned it_idx; //create an integer for the id2vert entry index + + //find out if the nodes for this fiber have already been created + it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node + if (it == id2vert.end()) { //if i[0] hasn't already been used + vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position + new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing + new_edge.v[0] = V.size(); //add the new edge to the edge + V.push_back(new_vertex); //add the new vertex to the vertex list + id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list + } + else { //if the vertex already exists + it_idx = std::distance(id2vert.begin(), it); + V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing + new_edge.v[0] = it_idx; + } + + it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID + if (it == id2vert.end()) { //if i[1] hasn't already been used + vertex new_vertex = new_edge[I - 1]; //create a new vertex, assign it a position + new_vertex.e[1].push_back(E.size()); //add the current edge as incoming + new_edge.v[1] = V.size(); //add the new vertex to the edge + V.push_back(new_vertex); //add the new vertex to the vertex list + id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list + } + else { //if the vertex already exists + it_idx = std::distance(id2vert.begin(), it); + V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming + new_edge.v[1] = it_idx; + } + + E.push_back(new_edge); //push the edge to the list + } + } + + /// Get adjacency matrix of the network + std::vector< typename std::vector > get_adj_mat() { + + unsigned n = V.size(); // get the number of vertices in the networks + + std::vector< typename std::vector > result(n, std::vector(n, 0)); // initialize every entry in the matrix to be 0 + result.resize(n); // resize rows + for (unsigned i = 0; i < n; i++) + result[i].resize(n); // resize columns + + for (unsigned i = 0; i < n; i++) { // for every vertex + unsigned num_out = V[i].e[0].size(); // number of outgoing edges of current vertex + if (num_out != 0) { + for (unsigned j = 0; j < num_out; j++) { + int edge_idx = V[i].e[0][j]; // get the jth out-going edge index of current vertex + int vertex_idx = E[edge_idx].v[1]; // get the ending vertex of specific out-going edge + result[i][vertex_idx] = 1; // can simply set to 1 if it is simple-graph + result[vertex_idx][i] = 1; // symmetric + } + } + } + + return result; + } + + /// Output the network as a string + std::string str(){ + + std::stringstream ss; + ss<<"Nodes ("< resample(T spacing){ + stim::network n; //create a new network that will be an exact copy, with resampled fibers + n.V = V; //copy all vertices + //n.NT = NT; //copy all the neuronal type information + n.E.resize(edges()); //allocate space for the edge list + + //copy all fibers, resampling them in the process + for(unsigned e = 0; e < edges(); e++){ //for each edge in the edge list + n.E[e] = E[e].resample(spacing); //resample the edge and copy it to the new network + } + + return n; //return the resampled network + } + + /// Calculate the total number of points on all edges. + unsigned total_points(){ + unsigned n = 0; + for(unsigned e = 0; e < E.size(); e++) + n += E[e].size(); + return n; + } + + //Copy the point cloud representing the centerline for the network into an array + void centerline_cloud(T* dst) { + size_t p; //stores the current edge point + size_t P; //stores the number of points in an edge + size_t i = 0; //index into the output array of points + for (size_t e = 0; e < E.size(); e++) { //for each edge in the network + P = E[e].size(); //get the number of points in this edge + for (p = 0; p < P; p++) { + dst[i * 3 + 0] = E[e][p][0]; + dst[i * 3 + 1] = E[e][p][1]; + dst[i * 3 + 2] = E[e][p][2]; + i++; + } + } + } + + // convert vec3 to array + void stim2array(float *a, stim::vec3 b){ + a[0] = b[0]; + a[1] = b[1]; + a[2] = b[2]; + } + + // convert vec3 to array in bunch + void edge2array(T* a, edge b){ + size_t n = b.size(); + for(size_t i = 0; i < n; i++){ + a[i * 3 + 0] = b[i][0]; + a[i * 3 + 1] = b[i][1]; + a[i * 3 + 2] = b[i][2]; + } + } + + // get list of metric + std::vector metric() { + std::vector result; + T m; + for (size_t e = 0; e < E.size(); e++) { + for (size_t p = 0; p < E[e].size(); p++) { + m = E[e].r(p); + result.push_back(m); + } + } + return result; + } + + /// Calculate the average magnitude across the entire network. + /// @param m is the magnitude value to use. The default is 0 (usually radius). + T average(unsigned m = 0){ + + T M, L; //allocate space for the total magnitude and length + M = L = 0; //initialize both the initial magnitude and length to zero + for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network + M += E[e].integrate(); //get the integrated magnitude + L += E[e].length(); //get the edge length + } + + return M / L; //return the average magnitude + } + + /// This function compares two networks and returns the percentage of the current network that is missing from A. + /// @param A is the network to compare to - the field is generated for A + /// @param sigma is the user-defined tolerance value - smaller values provide a stricter comparison + stim::network compare(stim::network A, float sigma, int device = -1){ + + stim::network R; //generate a network storing the result of the comparison + R = (*this); //initialize the result with the current network + + T *c; // centerline (array of double pointers) - points on kdtree must be double + size_t n_data = A.total_points(); // set the number of points + c = (T*) malloc(sizeof(T) * n_data * 3); // allocate an array to store all points in the data set + + unsigned t = 0; + for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in the network + for(unsigned p = 0; p < A.E[e].size(); p++){ //for each point in the edge + for(unsigned d = 0; d < 3; d++){ //for each coordinate + + c[t * 3 + d] = A.E[e][p][d]; //copy the point into the array c + } + t++; + } + } + + //generate a KD-tree for network A + size_t MaxTreeLevels = 3; // max tree level + +#ifdef __CUDACC__ + cudaSetDevice(device); + int current_device; + if (cudaGetDevice(¤t_device) == device) { + std::cout << "Using CUDA device " << device << " for calculations..." << std::endl; + } + stim::kdtree kdt; // initialize a pointer to a kd tree + + kdt.create(c, n_data, MaxTreeLevels); // build a KD tree + + for(unsigned e = 0; e < R.E.size(); e++){ //for each edge in A + //R.E[e].add_mag(0); //add a new magnitude for the metric + //size_t errormag_id = R.E[e].nmags() - 1; //get the id for the new magnitude + + size_t n = R.E[e].size(); // the number of points in current edge + T* queryPt = new T[3 * n]; + T* m1 = new T[n]; + T* dists = new T[n]; + size_t* nnIdx = new size_t[n]; + + T* d_dists; + T* d_m1; + cudaMalloc((void**)&d_dists, n * sizeof(T)); + cudaMalloc((void**)&d_m1, n * sizeof(T)); + + edge2array(queryPt, R.E[e]); + kdt.search(queryPt, n, nnIdx, dists); + + cudaMemcpy(d_dists, dists, n * sizeof(T), cudaMemcpyHostToDevice); // copy dists from host to device + + // configuration parameters + size_t threads = (1024>n)?n:1024; + size_t blocks = n/threads + (n%threads)?1:0; + + find_metric_parallel<<>>(d_m1, n, d_dists, sigma); //calculate the metric value based on the distance + + cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost); + + for(unsigned p = 0; p < n; p++){ + R.E[e].set_r(p, m1[p]); + } + + //d_set_mag<<>>(R.E[e].M, errormag_id, n, m1); + } + +#else + stim::kdtree kdt; + kdt.create(c, n_data, MaxTreeLevels); + + for(unsigned e = 0; e < R.E.size(); e++){ //for each edge in A + + size_t n = R.E[e].size(); // the number of points in current edge + T* query = new T[3 * n]; + T* m1 = new T[n]; + T* dists = new T[n]; + size_t* nnIdx = new size_t[n]; + + edge2array(query, R.E[e]); + + kdt.cpu_search(query, n, nnIdx, dists); //find the distance between A and the current network + + for(unsigned p = 0; p < R.E[e].size(); p++){ + m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma); //calculate the metric value based on the distance + R.E[e].set_r(p, m1[p]); //set the error for the second point in the segment + } + } +#endif + return R; //return the resulting network + } + + /// This function compares two networks and split the current one according to the nearest neighbor of each point in each edge + /// @param A is the network to split + /// @param B is the corresponding mapping network + /// @param sigma is the user-defined tolerance value - smaller values provide a stricter comparison + /// @param device is the device that user want to use + void split(stim::network A, stim::network B, float sigma, int device, float threshold){ + + T *c; + size_t n_data = B.total_points(); + c = (T*) malloc(sizeof(T) * n_data * 3); + + size_t NB = B.E.size(); // the number of edges in B + unsigned t = 0; + for(unsigned e = 0; e < NB; e++){ // for every edge in B + for(unsigned p = 0; p < B.E[e].size(); p++){ // for every points in B.E[e] + for(unsigned d = 0; d < 3; d++){ + + c[t * 3 + d] = B.E[e][p][d]; // convert to array + } + t++; + } + } + size_t MaxTreeLevels = 3; // max tree level + +#ifdef __CUDACC__ + cudaSetDevice(device); + stim::kdtree kdt; // initialize a pointer to a kd tree + + //compare each point in the current network to the field produced by A + kdt.create(c, n_data, MaxTreeLevels); // build a KD tree + + std::vector > relation; // the relationship between GT and T corresponding to NN + relation.resize(A.E.size()); + + for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in A + //A.E[e].add_mag(0); //add a new magnitude for the metric + //size_t errormag_id = A.E[e].nmags() - 1; + + size_t n = A.E[e].size(); // the number of edges in A + + T* queryPt = new T[3 * n]; // set of all the points in current edge + T* m1 = new T[n]; // array of metrics for every point in current edge + T* dists = new T[n]; // store the distances for every point in current edge + size_t* nnIdx = new size_t[n]; // store the indices for every point in current edge + + // define pointers in device + T* d_dists; + T* d_m1; + size_t* d_nnIdx; + + // allocate memory for defined pointers + cudaMalloc((void**)&d_dists, n * sizeof(T)); + cudaMalloc((void**)&d_m1, n * sizeof(T)); + cudaMalloc((void**)&d_nnIdx, n * sizeof(size_t)); + + edge2array(queryPt, A.E[e]); // convert edge to array + kdt.search(queryPt, n, nnIdx, dists); // search the tree to find the NN for every point in current edge + + cudaMemcpy(d_dists, dists, n * sizeof(T), cudaMemcpyHostToDevice); // copy dists from host to device + cudaMemcpy(d_nnIdx, nnIdx, n * sizeof(size_t), cudaMemcpyHostToDevice); // copy Idx from host to device + + // configuration parameters + size_t threads = (1024>n)?n:1024; // test to see whether the number of point in current edge is more than 1024 + size_t blocks = n/threads + (n%threads)?1:0; + + find_metric_parallel<<>>(d_m1, n, d_dists, sigma); // calculate the metrics in parallel + + cudaMemcpy(m1, d_m1, n * sizeof(T), cudaMemcpyDeviceToHost); + + for(unsigned p = 0; p < n; p++){ + A.E[e].set_r(p, m1[p]); // set the error(radius) value to every point in current edge + } + + relation[e].resize(n); // resize every edge relation size + + unsigned* d_relation; + cudaMalloc((void**)&d_relation, n * sizeof(unsigned)); // allocate memory + + std::vector edge_point_num(NB); // %th element is the number of points that %th edge has + for(unsigned ee = 0; ee < NB; ee++) + edge_point_num[ee] = B.E[ee].size(); + + size_t* d_edge_point_num; + cudaMalloc((void**)&d_edge_point_num, NB * sizeof(size_t)); + cudaMemcpy(d_edge_point_num, &edge_point_num[0], NB * sizeof(size_t), cudaMemcpyHostToDevice); + + find_edge_index_parallel<<>>(d_nnIdx, n, d_relation, d_edge_point_num, NB); // find the edge corresponding to the array index in parallel + + cudaMemcpy(&relation[e][0], d_relation, n * sizeof(unsigned), cudaMemcpyDeviceToHost); //copy relationship from device to host + } +#else + stim::kdtree kdt; + kdt.create(c, n_data, MaxTreeLevels); + + std::vector> relation; // the mapping relationship between two networks + relation.resize(A.E.size()); + for(unsigned i = 0; i < A.E.size(); i++) + relation[i].resize(A.E[i].size()); + + std::vector edge_point_num(NB); //%th element is the number of points that %th edge has + for(unsigned ee = 0; ee < NB; ee++) + edge_point_num[ee] = B.E[ee].size(); + + for(unsigned e = 0; e < A.E.size(); e++){ //for each edge in A + + size_t n = A.E[e].size(); //the number of edges in A + + T* queryPt = new T[3 * n]; + T* m1 = new T[n]; + T* dists = new T[n]; //store the distances + size_t* nnIdx = new size_t[n]; //store the indices + + edge2array(queryPt, A.E[e]); + kdt.search(queryPt, n, nnIdx, dists); + + for(unsigned p = 0; p < A.E[e].size(); p++){ + m1[p] = 1.0f - gaussianFunction((T)dists[p], sigma); //calculate the metric value based on the distance + A.E[e].set_r(p, m1[p]); //set the error for the second point in the segment + + unsigned id = 0; //mapping edge's idx + size_t num = 0; //total number of points before #th edge + for(unsigned i = 0; i < NB; i++){ + num += B.E[i].size(); + if(nnIdx[p] < num){ //find the edge it belongs to + relation[e][p] = id; + break; + } + id++; //current edge won't be the one, move to next edge + } + } + } +#endif + E = A.E; + V = A.V; + + unsigned int id = 0; // split value + for(unsigned e = 0; e < E.size(); e++){ // for every edge + for(unsigned p = 0; p < E[e].size() - 1; p++){ // for every point in each edge + int t = (int)(E[e].length() / sigma * 2); + if (t <= 20) + threshold_fac = E[e].size(); + else + threshold_fac = (E[e].length() / sigma * 2)/10; + if(relation[e][p] != relation[e][p + 1]){ // find the nearest edge changing point + id = p + 1; // if there is no change in NN + if(id < threshold_fac || (E[e].size() - id) < threshold_fac) + id = E[e].size() - 1; // extreme situation is not acceptable + else + break; + } + if(p == E[e].size() - 2) // if there is no splitting index, set the id to the last point index of current edge + id = E[e].size() - 1; + } + //unsigned errormag_id = E[e].nmags() - 1; + T G = 0; // test to see whether it has its nearest neighbor + for(unsigned i = 0; i < E[e].size(); i++) + G += E[e].r(i); // won't split special edges + if(G / E[e].size() > threshold) // should based on the color map + id = E[e].size() - 1; // set split idx to outgoing direction vertex + + std::vector tmpe; + tmpe.resize(2); + tmpe = E[e].split(id); + vertex tmpv = stim::vec3(-1, -1, 0); // store the split point as vertex + if(tmpe.size() == 2){ + relation.resize(relation.size() + 1); + for(unsigned d = id; d < E[e].size(); d++) + relation[relation.size() - 1].push_back(relation[e][d]); + tmpe[0].v[0] = E[e].v[0]; // begining vertex of first half edge -> original begining vertex + tmpe[1].v[1] = E[e].v[1]; // ending vertex of second half edge -> original ending vertex + tmpv = E[e][id]; + V.push_back(tmpv); + tmpe[0].v[1] = (unsigned)V.size() - 1; // ending vertex of first half edge -> new vertex + tmpe[1].v[0] = (unsigned)V.size() - 1; // begining vertex of second half edge -> new vertex + edge tmp(E[e]); + E[e] = tmpe[0]; // replace original edge by first half edge + E.push_back(tmpe[1]); // push second half edge to the last + V[V.size() - 1].e[1].push_back(e); // push first half edge to the incoming of new vertex + V[V.size() - 1].e[0].push_back((unsigned)E.size() - 1); // push second half edge to the outgoing of new vertex + for(unsigned i = 0; i < V[tmp.v[1]].e[1].size(); i++) // find the incoming edge of original ending vertex + if(V[tmp.v[1]].e[1][i] == e) + V[tmp.v[1]].e[1][i] = (unsigned)E.size() - 1; // set to new edge + } + } + } + + /// This function compares two splitted networks and yields a mapping relationship between them according to NN + /// @param B is the network that the current network is going to map to + /// @param C is the mapping relationship: C[e1] = _e1 means e1 edge in current network is mapping the _e1 edge in B + /// @param device is the device that user want to use + void mapping(stim::network B, std::vector &C, int device, float threshold){ + stim::network A; //generate a network storing the result of the comparison + A = (*this); + + size_t n = A.E.size(); // the number of edges in A + size_t NB = B.E.size(); // the number of edges in B + + C.resize(A.E.size()); + + T *c; // centerline (array of double pointers) - points on kdtree must be double + size_t n_data = B.total_points(); // set the number of points + c = (T*) malloc(sizeof(T) * n_data * 3); + + unsigned t = 0; + for(unsigned e = 0; e < NB; e++){ // for each edge in the network + for(unsigned p = 0; p < B.E[e].size(); p++){ // for each point in the edge + for(unsigned d = 0; d < 3; d++){ // for each coordinate + + c[t * 3 + d] = B.E[e][p][d]; + } + t++; + } + } + + //generate a KD-tree for network A + //float metric = 0.0; // initialize metric to be returned after comparing the network + size_t MaxTreeLevels = 3; // max tree level + +#ifdef __CUDACC__ + cudaSetDevice(device); + stim::kdtree kdt; // initialize a pointer to a kd tree + + kdt.create(c, n_data, MaxTreeLevels); // build a KD tree + + for(unsigned e = 0; e < n; e++){ //for each edge in A + //size_t errormag_id = A.E[e].nmags() - 1; //get the id for the new magnitude + + //pre-judge to get rid of impossibly mapping edges + T M = 0; + for(unsigned p = 0; p < A.E[e].size(); p++) + M += A.E[e].r(p); + M = M / A.E[e].size(); + if(M > threshold) + C[e] = (unsigned)-1; //set the nearest edge of impossibly mapping edges to maximum of unsigned + else{ + T* queryPt = new T[3]; + T* dists = new T[1]; + size_t* nnIdx = new size_t[1]; + + stim2array(queryPt, A.E[e][A.E[e].size()/2]); + kdt.search(queryPt, 1, nnIdx, dists); + + unsigned id = 0; //mapping edge's idx + size_t num = 0; //total number of points before #th edge + for(unsigned i = 0; i < NB; i++){ + num += B.E[i].size(); + if(nnIdx[0] < num){ + C[e] = id; + break; + } + id++; + } + } + } +#else + stim::kdtree kdt; + kdt.create(c, n_data, MaxTreeLevels); + T *dists = new T[1]; // near neighbor distances + size_t *nnIdx = new size_t[1]; // near neighbor indices // allocate near neigh indices + + stim::vec3 p0, p1; + T* queryPt = new T[3]; + + for(unsigned int e = 0; e < R.E.size(); e++){ // for each edge in A + T M; // the sum of metrics of current edge + for(unsigned p = 0; p < R.E[e].size(); p++) + M += A.E[e].r(p); + M = M / A.E[e].size(); + if(M > threshold) + C[e] = (unsigned)-1; + else{ // if it should have corresponding edge in B, then... + p1 = R.E[e][R.E[e].size()/2]; + stim2array(queryPt, p1); + kdt.cpu_search(queryPt, 1, nnIdx, dists); // search the tree + + unsigned id = 0; //mapping edge's idx + size_t num = 0; //total number of points before #th edge + for(unsigned i = 0; i < NB; i++){ + num += B.E[i].size(); + if(nnIdx[0] < num){ + C[e] = id; + break; + } + id++; + } + } + } +#endif + } + + /// Returns the number of magnitude values stored in each edge. This should be uniform across the network. + //unsigned nmags(){ + // return E[0].nmags(); + //} + // split a string in text by the character sep + stim::vec split(std::string &text, char sep) + { + stim::vec tokens; + std::size_t start = 0, end = 0; + while ((end = text.find(sep, start)) != std::string::npos) { + tokens.push_back(atof(text.substr(start, end - start).c_str())); + start = end + 1; + } + tokens.push_back(atof(text.substr(start).c_str())); + return tokens; + } + // load a network in text file to a network class + void load_txt(std::string filename) + { + std::vector file_contents; + std::ifstream file(filename.c_str()); + std::string line; + std::vector id2vert; //this list stores the vertex ID associated with each network vertex + //for each line in the text file, store them as strings in file_contents + while (std::getline(file, line)) + { + std::stringstream ss(line); + file_contents.push_back(ss.str()); + } + unsigned int numEdges = atoi(file_contents[0].c_str()); //number of edges in the network + unsigned int I = atoi(file_contents[1].c_str()) ; //calculate the number of points3d on the first edge + unsigned int count = 1; unsigned int k = 2; // count is global counter through the file contents, k is for the vertices on the edges + // for each edge in the network. + for (unsigned int i = 0; i < numEdges; i ++ ) + { + // pre allocate a position vector p with number of points3d on the edge p + std::vector< stim::vec > p(0, I); + // for each point on the nth edge + for (unsigned int j = k; j < I + k; j++) + { + // split the points3d of floats with separator space and form a float3 position vector out of them + p.push_back(split(file_contents[j], ' ')); + } + count += p.size() + 1; // increment count to point at the next edge in the network + I = atoi(file_contents[count].c_str()); // read in the points3d at the next edge and convert it to an integer + k = count + 1; + edge new_edge = p; // create an edge with a vector of points3d on the edge + E.push_back(new_edge); // push the edge into the network + } + unsigned int numVertices = atoi(file_contents[count].c_str()); // this line in the text file gives the number of distinct vertices + count = count + 1; // this line of text file gives the first verrtex + // push each vertex into V + for (unsigned int i = 0; i < numVertices; i ++) + { + vertex new_vertex = split(file_contents[count], ' '); + V.push_back(new_vertex); + count += atoi(file_contents[count + 1].c_str()) + 2; // Skip number of edge ids + 2 to point to the next vertex + } + } // end load_txt function + + // strTxt returns a string of edges + std::string + strTxt(std::vector< stim::vec > p) + { + std::stringstream ss; + std::stringstream oss; + for(unsigned int i = 0; i < p.size(); i++){ + ss.str(std::string()); + for(unsigned int d = 0; d < 3; d++){ + ss<