#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(); } 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(); } //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(); } //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 } } /// 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()); auto 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