diff --git a/stim/biomodels/centerline.h b/stim/biomodels/centerline.h index 63d8f38..e118a5e 100644 --- a/stim/biomodels/centerline.h +++ b/stim/biomodels/centerline.h @@ -190,7 +190,7 @@ public: 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 + //double sum=0; //distance summation size_t N = size(); diff --git a/stim/biomodels/network.h b/stim/biomodels/network.h index 0d17ae7..b54aeb4 100644 --- a/stim/biomodels/network.h +++ b/stim/biomodels/network.h @@ -33,16 +33,15 @@ class network{ public: unsigned v[2]; //unique id's designating the starting and ending // default constructor - edge() : cylinder() - { + 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 p is an array of positions in space - edge(std::vector< stim::vec3 > p) : cylinder(p){} + edge(centerline p) : cylinder(p){} - /// Copy constructor creates an edge from a fiber + /// Copy constructor creates an edge from a cylinder edge(stim::cylinder f) : cylinder(f) {} /// Resamples an edge by calling the fiber resampling function @@ -291,7 +290,7 @@ public: std::vector< stim::vec > c; //allocate an array of points for the vessel centerline O.getLine(l, c); //get the fiber centerline - std::vector< stim::vec3 > c3(c.size()); + stim::centerline c3(c.size()); for(size_t j = 0; j < c.size(); j++) c3[j] = c[j]; @@ -431,7 +430,7 @@ public: /// @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){ + 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 @@ -471,6 +470,7 @@ public: T* queryPt = new T[3]; 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 for(unsigned p = 0; p < R.E[e].size(); p++){ //for each point in the edge @@ -479,7 +479,7 @@ public: kdt.search(queryPt, 1, nnIdx, dists); //find the distance between A and the current network m1 = 1.0f - gaussianFunction((T)dists[0], sigma); //calculate the metric value based on the distance - R.E[e].set_mag(m1, p, 1); //set the error for the second point in the segment + R.E[e].set_mag(errormag_id, p, m1); //set the error for the second point in the segment } } diff --git a/stim/visualization/cylinder.h b/stim/visualization/cylinder.h index 1f1fb10..58a2b8b 100644 --- a/stim/visualization/cylinder.h +++ b/stim/visualization/cylinder.h @@ -11,19 +11,10 @@ namespace stim template class cylinder : public centerline { protected: - //stim::circle s; //an arbitrary circle - //std::vector > e; //an array of circles that store the centerline - //std::vector > norms; std::vector< stim::vec3 > U; //stores the array of U vectors defining the Frenet frame std::vector< std::vector > M; //stores a list of magnitudes for each point in the centerline (assuming mags[0] is the radius) - //std::vector< T > L; //length of the cylinder at each position (pre-integration) - - //using stim::centerline::c; - //using stim::centerline::N; - using stim::centerline::size; - using stim::centerline::at; using stim::centerline::findIdx; //calculates the U values for each point to initialize the frenet frame @@ -40,33 +31,18 @@ protected: U[size() - 1] = c.U; //for the last point, duplicate the final frenet frame vector } - ///adds a magnitude to each point in the cylinder - void add_mag(T 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 - M[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 - M[i].push_back(val[i]); //add this value to the magnitude vector at each point - } - public: - 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 = "<< M[i][0]<<" u = "<::size; + using stim::centerline::at; + + cylinder() : centerline(){} + + cylinder(centerline& c) : centerline(c) { + init(); } - cylinder(centerline& c, T r = 0) : centerline(c) { + cylinder(centerline& c, T r) : centerline(c) { init(); add_mag(r); } @@ -81,7 +57,7 @@ public: ///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 m(T l, int idx, size_t i = 0) { + T m(T l, int idx, size_t i) { T a = (l - L[idx]) / (L[idx + 1] - L[idx]); T v2 = (M[idx][i] + (M[idx + 1][i] - M[idx][i])*a); @@ -100,6 +76,37 @@ public: return m(l, idx, i); } + /// 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 m(unsigned i, unsigned m) { + return M[i][m]; + } + + ///adds a magnitude to each point in the cylinder + void add_mag(T 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 + M[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 + M[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_mag(size_t m, size_t i, T v) { + M[i][m] = v; + } + + size_t nmags() { + if (M.size() == 0) return 0; + else return M[0].size(); + } + ///Returns a circle representing the cylinder cross section at point i stim::circle circ(size_t i, size_t m = 0) { return stim::circle(at(i), M[i][m], d(i), U[i]); @@ -153,6 +160,61 @@ public: 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 = " << M[i][0] << " 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(size_t m = 0) { + + T sum = 0; //initialize the integral to zero + T m0, m1; //allocate space for both magnitudes in a single segment + m0 = M[0][m]; //initialize the first point and magnitude to the first point in the cylinder + T len = L[0]; //allocate space for the segment length + + + for (unsigned p = 1; p < size(); p++) { //for every consecutive point in the cylinder + m1 = M[p][m]; + if (p > 1) len = (L[p - 1] - L[p - 2]); //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.M[i][mag] = m(t, mag); //retrieve the interpolated magnitude from the current cylinder and store it in the new one + } + } + + + return c; + } + /* ///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and magnitudes (inM) diff --git a/stim/visualization/gl_network.h b/stim/visualization/gl_network.h index f1b357f..f4f25f9 100644 --- a/stim/visualization/gl_network.h +++ b/stim/visualization/gl_network.h @@ -55,7 +55,7 @@ public: 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].ri(p, m)); //set the texture coordinate based on the specified magnitude index + glTexCoord1f(E[e].m(p, m)); //set the texture coordinate based on the specified magnitude index } glEnd(); } -- libgit2 0.21.4