Commit 0311ab81ccc840968f323354738811bd1a1ee33e

Authored by Pavel Govyadinov
1 parent eb119a4f

fixed some issues with std98 in filename, fixes to glObj, debugging code added, …

…Cylinder now uses centerline
stim/gl/gl_spider.h
@@ -24,7 +24,6 @@ @@ -24,7 +24,6 @@
24 #include <stim/visualization/glObj.h> 24 #include <stim/visualization/glObj.h>
25 #include <vector> 25 #include <vector>
26 #include <stim/cuda/branch_detection.cuh> 26 #include <stim/cuda/branch_detection.cuh>
27 -#include "../../../volume-spider/fiber.h"  
28 #include "../../../volume-spider/glnetwork.h" 27 #include "../../../volume-spider/glnetwork.h"
29 #include <stim/visualization/cylinder.h> 28 #include <stim/visualization/cylinder.h>
30 #include <stim/cuda/testKernel.cuh> 29 #include <stim/cuda/testKernel.cuh>
@@ -115,6 +114,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -115,6 +114,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
115 std::vector< stim::vec<float> > cM; //Magnitude of line currently being traced. 114 std::vector< stim::vec<float> > cM; //Magnitude of line currently being traced.
116 115
117 stim::glnetwork<float> nt; //object for storing the network. 116 stim::glnetwork<float> nt; //object for storing the network.
  117 + stim::glObj<float> sk;
118 118
119 stim::vec<float> rev; //reverse vector; 119 stim::vec<float> rev; //reverse vector;
120 stim::camera camSel; 120 stim::camera camSel;
@@ -228,51 +228,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -228,51 +228,6 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
228 } 228 }
229 229
230 230
231 - ///subject to change.  
232 - ///finds branches.  
233 - ///depreciated  
234 - void  
235 - branchDetection()  
236 - {  
237 - setMatrix();  
238 - glCallList(dList+3);  
239 - std::vector< stim::vec<float> > result = find_branch(  
240 - btexbufferID, GL_TEXTURE_2D, 16, 216);  
241 - stim::vec3<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);  
242 - if(!result.empty())  
243 - {  
244 - for(int i = 1; i < result.size(); i++)  
245 - {  
246 - stim::vec<float> cylp(  
247 - 0.5 * cos(2*M_PI*(result[i][1])),  
248 - 0.5 * sin(2*M_PI*(result[i][1])),  
249 - result[i][0]-0.5,  
250 - 1.0);  
251 - cylp = cT*cylp;  
252 -  
253 - stim::vec3<float> vec(  
254 - cylp[0]*S[0]*R[0],  
255 - cylp[1]*S[1]*R[1],  
256 - cylp[2]*S[2]*R[2]);  
257 - stim::vec3<float> seeddir(-p[0] + cylp[0]*S[0]*R[0],  
258 - -p[1] + cylp[1]*S[1]*R[1],  
259 - -p[2] + cylp[2]*S[2]*R[2]);  
260 - seeddir = seeddir.norm();  
261 - float seedm = m[0];  
262 -// Uncomment for global run  
263 - if(  
264 - !(vec[0] > size[0] || vec[1] > size[1]  
265 - || vec[2] > size[2] || vec[0] < 0  
266 - || vec[1] < 0 || vec[2] < 0))  
267 - {  
268 - setSeed(vec);  
269 - setSeedVec(seeddir);  
270 - setSeedMag(seedm);  
271 - }  
272 - }  
273 - }  
274 -  
275 - }  
276 231
277 232
278 ///finds all the branches in the a given fiber. 233 ///finds all the branches in the a given fiber.
@@ -314,6 +269,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -314,6 +269,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
314 stim::vec3<float> v = cyl.surf(pval, result[i][0]); 269 stim::vec3<float> v = cyl.surf(pval, result[i][0]);
315 stim::vec3<float> di = cyl.p(pval); 270 stim::vec3<float> di = cyl.p(pval);
316 float rad = cyl.r(pval); 271 float rad = cyl.r(pval);
  272 +
  273 + std::cout << v << " " << di << " " << rad << std::endl;
317 if( 274 if(
318 !(v[0] > size[0] || v[1] > size[1] 275 !(v[0] > size[0] || v[1] > size[1]
319 || v[2] > size[2] || v[0] < 0 276 || v[2] > size[2] || v[0] < 0
@@ -768,7 +725,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -768,7 +725,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
768 findOptimalScale(); 725 findOptimalScale();
769 Unbind(); 726 Unbind();
770 Bind(btexbufferID, bfboID, 27); 727 Bind(btexbufferID, bfboID, 27);
771 - branchDetection(); 728 + branchDetection2();
772 Unbind(); 729 Unbind();
773 return current_cost; 730 return current_cost;
774 } 731 }
@@ -1297,7 +1254,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1297,7 +1254,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1297 void 1254 void
1298 saveNetwork(std::string name) 1255 saveNetwork(std::string name)
1299 { 1256 {
1300 - stim::glObj<float> sk; 1257 +/* stim::glObj<float> sk;
1301 for(int i = 0; i < nt.sizeE(); i++) 1258 for(int i = 0; i < nt.sizeE(); i++)
1302 { 1259 {
1303 std::vector<stim::vec< float > > cm = nt.getEdgeCenterLineMag(i); 1260 std::vector<stim::vec< float > > cm = nt.getEdgeCenterLineMag(i);
@@ -1310,7 +1267,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1310,7 +1267,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1310 } 1267 }
1311 sk.End(); 1268 sk.End();
1312 } 1269 }
1313 - sk.save(name); 1270 +*/ sk.save(name);
1314 } 1271 }
1315 1272
1316 ///Depreciated, but might be reused later() 1273 ///Depreciated, but might be reused later()
@@ -1318,7 +1275,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1318,7 +1275,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1318 stim::glObj<float> 1275 stim::glObj<float>
1319 getNetwork() 1276 getNetwork()
1320 { 1277 {
1321 - 1278 + return sk;
1322 } 1279 }
1323 1280
1324 ///returns a COPY of the entire stim::glnetwork object. 1281 ///returns a COPY of the entire stim::glnetwork object.
@@ -1387,14 +1344,14 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1387,14 +1344,14 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1387 void 1344 void
1388 MonteCarloDirectionVectors(int nSamples, float solidAngle = 2*M_PI) 1345 MonteCarloDirectionVectors(int nSamples, float solidAngle = 2*M_PI)
1389 { 1346 {
1390 - float PHI[2], Z[2], range;  
1391 - PHI[0] = asin(solidAngle/2);  
1392 - PHI[1] = asin(0); 1347 +// float PHI[2];//, Z[2];//, range;
  1348 +// PHI[0] = asin(solidAngle/2);
  1349 +// PHI[1] = asin(0);
1393 1350
1394 - Z[0] = cos(PHI[0]);  
1395 - Z[1] = cos(PHI[1]); 1351 +// Z[0] = cos(PHI[0]);
  1352 +// Z[1] = cos(PHI[1]);
1396 1353
1397 - range = Z[0] - Z[1]; 1354 +// range = Z[0] - Z[1];
1398 1355
1399 std::vector<stim::vec3<float> > vecsUni; 1356 std::vector<stim::vec3<float> > vecsUni;
1400 for(int i = 0; i < numSamplesPos; i++) 1357 for(int i = 0; i < numSamplesPos; i++)
@@ -1506,10 +1463,13 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1506,10 +1463,13 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1506 setPosition(curSeed); 1463 setPosition(curSeed);
1507 setDirection(curSeedVec); 1464 setDirection(curSeedVec);
1508 setMagnitude(curSeedMag); 1465 setMagnitude(curSeedMag);
  1466 + std::cout << "This was done" << std::endl;
  1467 + findOptimalScale();
  1468 +
1509 // cL.push_back(curSeed); 1469 // cL.push_back(curSeed);
1510 // cM.push_back(curSeedMag); 1470 // cM.push_back(curSeedMag);
1511 // cD.push_back(curSeedMag); 1471 // cD.push_back(curSeedMag);
1512 - pair<stim::fiber<float>, int> a = traceLine(p, m, min_cost); 1472 + traceLine(p, m, min_cost);
1513 } 1473 }
1514 } 1474 }
1515 1475
@@ -1556,17 +1516,17 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1556,17 +1516,17 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1556 ds[0], ds[1], ds[2], 1516 ds[0], ds[1], ds[2],
1557 ups[0], ups[1], ups[2]); 1517 ups[0], ups[1], ups[2]);
1558 1518
1559 -// sk.Render();  
1560 - nt.Render(); 1519 + sk.Render();
  1520 +// nt.Render();
1561 1521
1562 CHECK_OPENGL_ERROR 1522 CHECK_OPENGL_ERROR
1563 1523
1564 1524
1565 -// glLoadName((int) sk.numL());  
1566 - glLoadName(nt.sizeE()); 1525 + glLoadName((int) sk.numL());
  1526 +// glLoadName(nt.sizeE());
1567 1527
1568 -// sk.RenderLine(cL);  
1569 - nt.RenderLine(cL); 1528 + sk.RenderLine(cL);
  1529 +// nt.RenderLine(cL);
1570 1530
1571 // glPopName(); 1531 // glPopName();
1572 glFlush(); 1532 glFlush();
@@ -1615,8 +1575,9 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1615,8 +1575,9 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1615 { 1575 {
1616 cL.clear(); 1576 cL.clear();
1617 cM.clear(); 1577 cM.clear();
1618 - } 1578 + }
1619 1579
  1580 +/*
1620 void 1581 void
1621 addToNetwork(pair<stim::fiber<float>, int> in, stim::vec3<float> spos, 1582 addToNetwork(pair<stim::fiber<float>, int> in, stim::vec3<float> spos,
1622 stim::vec<float> smag, stim::vec3<float> sdir) 1583 stim::vec<float> smag, stim::vec3<float> sdir)
@@ -1664,13 +1625,23 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1664,13 +1625,23 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1664 network_time += nt * 1000.0; 1625 network_time += nt * 1000.0;
1665 #endif 1626 #endif
1666 } 1627 }
  1628 +*/
  1629 + void
  1630 + addToNetwork(std::vector<stim::vec3<float> > L, std::vector<stim::vec<float> > M)
  1631 + {
  1632 + if(L.size() > 3)
  1633 + {
  1634 + sk.Begin(stim::OBJ_LINE);
  1635 + for(int i = 0; i < L.size(); i++)
  1636 + {
  1637 + sk.TexCoord(M[i][0]);
  1638 + sk.Vertex(L[i][0], L[i][1], L[i][2]);
  1639 + }
  1640 + sk.End();
1667 1641
1668 -// void  
1669 -// addToNetwork(pair<stim::fiber<float>, int> in, stim::vec3<float> spos,  
1670 -// stim::vec<float> smag, stim::vec3<float> sdir)  
1671 -// {  
1672 -//  
1673 -// } 1642 + nt.addEdge(L,M);
  1643 + }
  1644 + }
1674 1645
1675 1646
1676 void 1647 void
@@ -1681,11 +1652,13 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1681,11 +1652,13 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1681 1652
1682 } 1653 }
1683 1654
1684 - std::pair<stim::fiber<float>, int > 1655 + void
1685 traceLine(stim::vec3<float> pos, stim::vec<float> mag, int min_cost) 1656 traceLine(stim::vec3<float> pos, stim::vec<float> mag, int min_cost)
1686 { 1657 {
1687 //starting (seed) position and magnitude. 1658 //starting (seed) position and magnitude.
1688 stim::vec3<float> spos = getPosition(); 1659 stim::vec3<float> spos = getPosition();
  1660 +// std::cout << "I did this" << std::endl;
  1661 +// findOptimalScale();
1689 stim::vec<float> smag = getMagnitude(); 1662 stim::vec<float> smag = getMagnitude();
1690 stim::vec3<float> sdir = getDirection(); 1663 stim::vec3<float> sdir = getDirection();
1691 1664
@@ -1693,8 +1666,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1693,8 +1666,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1693 // sk.Begin(stim::OBJ_LINE); 1666 // sk.Begin(stim::OBJ_LINE);
1694 1667
1695 1668
1696 -// sk.createFromSelf(GL_SELECT);  
1697 - nt.createFromSelf(GL_SELECT); 1669 + sk.createFromSelf(GL_SELECT);
  1670 +// nt.createFromSelf(GL_SELECT);
1698 1671
1699 cL.push_back(pos); 1672 cL.push_back(pos);
1700 cM.push_back(mag); 1673 cM.push_back(mag);
@@ -1712,10 +1685,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1712,10 +1685,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1712 running = false; 1685 running = false;
1713 // sk.End(); 1686 // sk.End();
1714 branchDetection2(); 1687 branchDetection2();
1715 - pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);  
1716 - addToNetwork(a, spos, smag, sdir);  
1717 -// std::cout << "the cost of " << cost << " > " << min_cost << std::endl;  
1718 - return a; 1688 + addToNetwork(cL, cM);
  1689 + std::cout << "the cost of " << cost << " > " << min_cost << std::endl;
1719 break; 1690 break;
1720 } else { 1691 } else {
1721 //Have we found the edge of the map? 1692 //Have we found the edge of the map?
@@ -1726,10 +1697,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1726,10 +1697,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1726 { 1697 {
1727 running = false; 1698 running = false;
1728 branchDetection2(); 1699 branchDetection2();
1729 - pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);  
1730 - addToNetwork(a, spos, smag, sdir);  
1731 -// std::cout << "I hit and edge" << std::endl;  
1732 - return a; 1700 + addToNetwork(cL, cM);
  1701 + std::cout << "I hit and edge" << std::endl;
1733 break; 1702 break;
1734 } 1703 }
1735 //If this is the first step in the trace, 1704 //If this is the first step in the trace,
@@ -1744,10 +1713,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1744,10 +1713,8 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1744 if(mag[0] > 75 || mag[0] < 1){ 1713 if(mag[0] > 75 || mag[0] < 1){
1745 running = false; 1714 running = false;
1746 branchDetection2(); 1715 branchDetection2();
1747 - pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), -1);  
1748 - addToNetwork(a, spos, smag, sdir);  
1749 -// std::cout << "The templates are too big" << std::endl;  
1750 - return a; 1716 + addToNetwork(cL, cM);
  1717 + std::cout << "The templates are too big" << std::endl;
1751 break; 1718 break;
1752 } 1719 }
1753 else 1720 else
@@ -1755,12 +1722,10 @@ class gl_spider : public virtual gl_texture&lt;T&gt; @@ -1755,12 +1722,10 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
1755 h = selectObject(p, getDirection(), m[0]); 1722 h = selectObject(p, getDirection(), m[0]);
1756 //Have we hit something previously traced? 1723 //Have we hit something previously traced?
1757 if(h != -1){ 1724 if(h != -1){
1758 -// std::cout << "I hit the fiber " << h << std::endl; 1725 + std::cout << "I hit the fiber " << h << std::endl;
1759 running = false; 1726 running = false;
1760 branchDetection2(); 1727 branchDetection2();
1761 - pair<stim::fiber<float>, int> a(stim::fiber<float> (cL, cM), h);  
1762 - addToNetwork(a, spos, smag, sdir);  
1763 - return a; 1728 + addToNetwork(cL, cM);
1764 break; 1729 break;
1765 } 1730 }
1766 else { 1731 else {
stim/math/circle.h
@@ -48,9 +48,9 @@ public: @@ -48,9 +48,9 @@ public:
48 CUDA_CALLABLE 48 CUDA_CALLABLE
49 circle(T size, T z_pos = (T)0) : plane<T>() 49 circle(T size, T z_pos = (T)0) : plane<T>()
50 { 50 {
51 - init();  
52 center(stim::vec3<T>(0,0,z_pos)); 51 center(stim::vec3<T>(0,0,z_pos));
53 scale(size); 52 scale(size);
  53 + init();
54 } 54 }
55 55
56 ///create a rectangle from a center point, normal 56 ///create a rectangle from a center point, normal
@@ -88,6 +88,7 @@ public: @@ -88,6 +88,7 @@ public:
88 { 88 {
89 init(); 89 init();
90 setU(u); 90 setU(u);
  91 +// U = u;
91 center(c); 92 center(c);
92 normal(n); 93 normal(n);
93 scale(s); 94 scale(s);
@@ -174,6 +175,13 @@ public: @@ -174,6 +175,13 @@ public:
174 return p(x,y); 175 return p(x,y);
175 } 176 }
176 177
  178 + std::string str() const
  179 + {
  180 + std::stringstream ss;
  181 + ss << "(P=" << P.str() << ", N=" << N.str() << ", U=" << U.str() << ", Y=" << Y.str();
  182 + return ss.str();
  183 + }
  184 +
177 }; 185 };
178 } 186 }
179 #endif 187 #endif
stim/parser/filename.h
@@ -110,11 +110,11 @@ protected: @@ -110,11 +110,11 @@ protected:
110 unix_dir = unix_dir.substr(2, unix_dir.length()-2); //extract the directory structure 110 unix_dir = unix_dir.substr(2, unix_dir.length()-2); //extract the directory structure
111 } 111 }
112 112
113 - if(unix_dir.front() == '/'){ //if there is a leading slash 113 + if(unix_dir.at(0) == '/'){ //if there is a leading slash
114 relative = false; //the path is not relative 114 relative = false; //the path is not relative
115 unix_dir = unix_dir.substr(1, unix_dir.length() - 1); //remove the slash 115 unix_dir = unix_dir.substr(1, unix_dir.length() - 1); //remove the slash
116 } 116 }
117 - if(unix_dir.back() == '/') 117 + if(unix_dir.at(unix_dir.size()-1) == '/')
118 unix_dir = unix_dir.substr(0, unix_dir.length() - 1); 118 unix_dir = unix_dir.substr(0, unix_dir.length() - 1);
119 119
120 path = stim::parser::split(unix_dir, '/'); //split up the directory structure 120 path = stim::parser::split(unix_dir, '/'); //split up the directory structure
stim/visualization/cylinder.h
@@ -2,19 +2,27 @@ @@ -2,19 +2,27 @@
2 #define STIM_CYLINDER_H 2 #define STIM_CYLINDER_H
3 #include <iostream> 3 #include <iostream>
4 #include <stim/math/circle.h> 4 #include <stim/math/circle.h>
5 -#include <stim/math/vec3.h> 5 +#include <stim/biomodels/centerline.h>
6 6
7 7
8 namespace stim 8 namespace stim
9 { 9 {
10 template<typename T> 10 template<typename T>
11 class cylinder 11 class cylinder
  12 + : public centerline<T>
12 { 13 {
13 private: 14 private:
14 stim::circle<T> s; //an arbitrary circle 15 stim::circle<T> s; //an arbitrary circle
15 - std::vector<stim::circle<T> > e; 16 + std::vector<stim::circle<T> > e; //an array of circles that store the centerline
  17 +
  18 + std::vector<stim::vec3<T> > norms;
  19 + std::vector<stim::vec<T> > Us;
16 std::vector<stim::vec<T> > mags; 20 std::vector<stim::vec<T> > mags;
17 std::vector< T > L; //length of the cylinder at each position. 21 std::vector< T > L; //length of the cylinder at each position.
  22 +
  23 +
  24 + using stim::centerline<T>::c;
  25 + using stim::centerline<T>::N;
18 26
19 ///default init 27 ///default init
20 void 28 void
@@ -23,6 +31,20 @@ class cylinder @@ -23,6 +31,20 @@ class cylinder
23 31
24 } 32 }
25 33
  34 + ///Calculate the physical length of the fiber at each point in the fiber.
  35 + void
  36 + calculateL()
  37 + {
  38 +/* L.resize(inP.size());
  39 + T temp = (T)0;
  40 + L[0] = 0;
  41 + for(size_t i = 1; i < L.size(); i++)
  42 + {
  43 + temp += (inP[i-1] - inP[i]).len();
  44 + L[i] = temp;
  45 + }
  46 +*/ }
  47 +
26 ///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and radii (inM) 48 ///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and radii (inM)
27 void 49 void
28 init(std::vector<stim::vec3<T> > inP, std::vector<stim::vec<T> > inM) 50 init(std::vector<stim::vec3<T> > inP, std::vector<stim::vec<T> > inM)
@@ -31,6 +53,10 @@ class cylinder @@ -31,6 +53,10 @@ class cylinder
31 stim::vec3<float> v1; 53 stim::vec3<float> v1;
32 stim::vec3<float> v2; 54 stim::vec3<float> v2;
33 e.resize(inP.size()); 55 e.resize(inP.size());
  56 +
  57 + norms.resize(inP.size());
  58 + Us.resize(inP.size());
  59 +
34 if(inP.size() < 2) 60 if(inP.size() < 2)
35 return; 61 return;
36 62
@@ -46,7 +72,9 @@ class cylinder @@ -46,7 +72,9 @@ class cylinder
46 72
47 stim::vec3<T> dr = (inP[1] - inP[0]).norm(); 73 stim::vec3<T> dr = (inP[1] - inP[0]).norm();
48 s = stim::circle<T>(inP[0], inM[0][0], dr, stim::vec3<T>(1,0,0)); 74 s = stim::circle<T>(inP[0], inM[0][0], dr, stim::vec3<T>(1,0,0));
  75 + norms[0] = s.N;
49 e[0] = s; 76 e[0] = s;
  77 + Us[0] = e[0].U;
50 for(size_t i = 1; i < inP.size()-1; i++) 78 for(size_t i = 1; i < inP.size()-1; i++)
51 { 79 {
52 s.center(inP[i]); 80 s.center(inP[i]);
@@ -54,16 +82,24 @@ class cylinder @@ -54,16 +82,24 @@ class cylinder
54 v2 = (inP[i+1] - inP[i]).norm(); 82 v2 = (inP[i+1] - inP[i]).norm();
55 dr = (v1+v2).norm(); 83 dr = (v1+v2).norm();
56 s.normal(dr); 84 s.normal(dr);
  85 +
  86 + norms[i] = s.N;
  87 +
57 s.scale(inM[i][0]/inM[i-1][0]); 88 s.scale(inM[i][0]/inM[i-1][0]);
58 e[i] = s; 89 e[i] = s;
  90 + Us[i] = e[i].U;
59 } 91 }
60 92
61 int j = inP.size()-1; 93 int j = inP.size()-1;
62 s.center(inP[j]); 94 s.center(inP[j]);
63 dr = (inP[j] - inP[j-1]).norm(); 95 dr = (inP[j] - inP[j-1]).norm();
64 s.normal(dr); 96 s.normal(dr);
  97 +
  98 + norms[j] = s.N;
  99 +
65 s.scale(inM[j][0]/inM[j-1][0]); 100 s.scale(inM[j][0]/inM[j-1][0]);
66 e[j] = s; 101 e[j] = s;
  102 + Us[j] = e[j].U;
67 } 103 }
68 104
69 ///returns the direction vector at point idx. 105 ///returns the direction vector at point idx.
@@ -72,18 +108,42 @@ class cylinder @@ -72,18 +108,42 @@ class cylinder
72 { 108 {
73 if(idx == 0) 109 if(idx == 0)
74 { 110 {
75 - return (e[idx+1].P - e[idx].P).norm(); 111 + stim::vec3<T> temp(
  112 + c[idx+1][0]-c[idx][0],
  113 + c[idx+1][1]-c[idx][1],
  114 + c[idx+1][2]-c[idx][2]
  115 + );
  116 +// return (e[idx+1].P - e[idx].P).norm();
  117 + return (temp.norm());
76 } 118 }
77 - else if(idx == e.size()-1) 119 + else if(idx == N-1)
78 { 120 {
79 - return (e[idx].P - e[idx-1].P).norm(); 121 + stim::vec3<T> temp(
  122 + c[idx][0]-c[idx+1][0],
  123 + c[idx][1]-c[idx+1][1],
  124 + c[idx][2]-c[idx+1][2]
  125 + );
  126 + // return (e[idx].P - e[idx-1].P).norm();
  127 + return (temp.norm());
80 } 128 }
81 else 129 else
82 { 130 {
83 // return (e[idx+1].P - e[idx].P).norm(); 131 // return (e[idx+1].P - e[idx].P).norm();
84 - stim::vec3<float> v1 = (e[idx].P-e[idx-1].P).norm();  
85 - stim::vec3<float> v2 = (e[idx+1].P-e[idx].P).norm();  
86 - return (v1+v2).norm(); 132 +// stim::vec3<float> v1 = (e[idx].P-e[idx-1].P).norm();
  133 + stim::vec3<T> v1(
  134 + c[idx][0]-c[idx-1][0],
  135 + c[idx][1]-c[idx-1][1],
  136 + c[idx][2]-c[idx-1][2]
  137 + );
  138 +
  139 +// stim::vec3<float> v2 = (e[idx+1].P-e[idx].P).norm();
  140 + stim::vec3<T> v2(
  141 + c[idx+1][0]-c[idx][0],
  142 + c[idx+1][1]-c[idx][1],
  143 + c[idx+1][2]-c[idx][2]
  144 + );
  145 +
  146 + return (v1.norm()+v2.norm()).norm();
87 } 147 }
88 // return e[idx].N; 148 // return e[idx].N;
89 149
@@ -92,14 +152,17 @@ class cylinder @@ -92,14 +152,17 @@ class cylinder
92 stim::vec3<T> 152 stim::vec3<T>
93 d(T l, int idx) 153 d(T l, int idx)
94 { 154 {
95 - if(idx == 0 || idx == e.size()-1) 155 + if(idx == 0 || idx == N-1)
96 { 156 {
97 - return e[idx].N; 157 + return norms[idx];
  158 +// return e[idx].N;
98 } 159 }
99 else 160 else
100 { 161 {
  162 +
101 T rat = (l-L[idx])/(L[idx+1]-L[idx]); 163 T rat = (l-L[idx])/(L[idx+1]-L[idx]);
102 - return( e[idx].N + (e[idx+1].N - e[idx].N)*rat); 164 + return( norms[idx] + (norms[idx+1] - norms[idx])*rat);
  165 +// return( e[idx].N + (e[idx+1].N - e[idx].N)*rat);
103 } 166 }
104 } 167 }
105 168
@@ -137,6 +200,7 @@ class cylinder @@ -137,6 +200,7 @@ class cylinder
137 public: 200 public:
138 ///default constructor 201 ///default constructor
139 cylinder() 202 cylinder()
  203 + // : centerline<T>()
140 { 204 {
141 205
142 } 206 }
@@ -144,14 +208,18 @@ class cylinder @@ -144,14 +208,18 @@ class cylinder
144 ///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder. 208 ///constructor to create a cylinder from a set of points, radii, and the number of sides for the cylinder.
145 ///@param inP: Vector of stim vec3 composing the points of the centerline. 209 ///@param inP: Vector of stim vec3 composing the points of the centerline.
146 ///@param inM: Vector of stim vecs composing the radii of the centerline. 210 ///@param inM: Vector of stim vecs composing the radii of the centerline.
147 - cylinder(std::vector<stim::vec3<T> > inP, std::vector<stim::vec<T> > inM){ 211 + cylinder(std::vector<stim::vec3<T> > inP, std::vector<stim::vec<T> > inM)
  212 + : centerline<T>(inP)
  213 + {
148 init(inP, inM); 214 init(inP, inM);
149 } 215 }
150 216
151 217
152 ///Constructor defines a cylinder with centerline inP and magnitudes of zero 218 ///Constructor defines a cylinder with centerline inP and magnitudes of zero
153 ///@param inP: Vector of stim vec3 composing the points of the centerline 219 ///@param inP: Vector of stim vec3 composing the points of the centerline
154 - cylinder(std::vector< stim::vec3<T> > inP){ 220 + cylinder(std::vector< stim::vec3<T> > inP)
  221 + : centerline<T>(inP)
  222 + {
155 std::vector< stim::vec<T> > inM; //create an array of arbitrary magnitudes 223 std::vector< stim::vec<T> > inM; //create an array of arbitrary magnitudes
156 224
157 stim::vec<T> zero; 225 stim::vec<T> zero;
@@ -164,7 +232,7 @@ class cylinder @@ -164,7 +232,7 @@ class cylinder
164 ///Returns the number of points on the cylinder centerline 232 ///Returns the number of points on the cylinder centerline
165 233
166 unsigned int size(){ 234 unsigned int size(){
167 - return e.size(); 235 + return N;
168 } 236 }
169 237
170 238
@@ -180,9 +248,23 @@ class cylinder @@ -180,9 +248,23 @@ class cylinder
180 } 248 }
181 T l = pvalue*L[L.size()-1]; 249 T l = pvalue*L[L.size()-1];
182 int idx = findIdx(l); 250 int idx = findIdx(l);
183 - T rat = (l-L[idx])/(L[idx+1]-L[idx]);  
184 - return( e[idx].P + (e[idx+1].P-e[idx].P)*rat);  
185 - } 251 + return (p(l,idx));
  252 +/* T rat = (l-L[idx])/(L[idx+1]-L[idx]);
  253 + stim::vec3<T> v1(
  254 + c[idx][0],
  255 + c[idx][1],
  256 + c[idx][2]
  257 + );
  258 +
  259 + stim::vec3<T> v2(
  260 + c[idx+1][0],
  261 + c[idx+1][1],
  262 + c[idx+1][2]
  263 + );
  264 +
  265 +// return( e[idx].P + (e[idx+1].P-e[idx].P)*rat);
  266 + return( v1 + (v2 - v1)*rat);
  267 +*/ }
186 268
187 ///Returns a position vector at the given length into the fiber (based on the pvalue). 269 ///Returns a position vector at the given length into the fiber (based on the pvalue).
188 ///Interpolates the radius along the line. 270 ///Interpolates the radius along the line.
@@ -192,7 +274,19 @@ class cylinder @@ -192,7 +274,19 @@ class cylinder
192 p(T l, int idx) 274 p(T l, int idx)
193 { 275 {
194 T rat = (l-L[idx])/(L[idx+1]-L[idx]); 276 T rat = (l-L[idx])/(L[idx+1]-L[idx]);
195 - return( e[idx].P + (e[idx+1].P-e[idx].P)*rat); 277 + stim::vec3<T> v1(
  278 + c[idx][0],
  279 + c[idx][1],
  280 + c[idx][2]
  281 + );
  282 +
  283 + stim::vec3<T> v2(
  284 + c[idx+1][0],
  285 + c[idx+1][1],
  286 + c[idx+1][2]
  287 + );
  288 +// return( e[idx].P + (e[idx+1].P-e[idx].P)*rat);
  289 + return( v1 + (v2-v1)*rat);
196 // return( 290 // return(
197 // return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx]))); 291 // return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
198 } 292 }
@@ -209,8 +303,11 @@ class cylinder @@ -209,8 +303,11 @@ class cylinder
209 } 303 }
210 T l = pvalue*L[L.size()-1]; 304 T l = pvalue*L[L.size()-1];
211 int idx = findIdx(l); 305 int idx = findIdx(l);
212 - return (e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*((l-L[idx])/(L[idx+1]- L[idx])));  
213 - } 306 + return (r(l,idx));
  307 +/* T rat = (l-L[idx])/(L[idx+1]-L[idx]);
  308 +// return (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat);
  309 + return (e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*rat);
  310 +*/ }
214 311
215 ///Returns a radius at the given length into the fiber (based on the pvalue). 312 ///Returns a radius at the given length into the fiber (based on the pvalue).
216 ///Interpolates the position along the line. 313 ///Interpolates the position along the line.
@@ -220,7 +317,23 @@ class cylinder @@ -220,7 +317,23 @@ class cylinder
220 r(T l, int idx) 317 r(T l, int idx)
221 { 318 {
222 T rat = (l-L[idx])/(L[idx+1]-L[idx]); 319 T rat = (l-L[idx])/(L[idx+1]-L[idx]);
223 - return( e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*rat); 320 + T v1 = (e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*rat);
  321 + T v3 = (Us[idx].len() + (Us[idx+1].len() - Us[idx].len())*rat);
  322 + T v2 = (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat);
  323 +// std::cout << (float)v1 = (float) v2 << std::endl;
  324 + if(abs(v3 - v1) >= 10e-6)
  325 + {
  326 + std::cout << "-------------------------" << std::endl;
  327 + std::cout << e[idx].str() << std::endl << std::endl;
  328 + std::cout << Us[idx] << std::endl;
  329 + std::cout << (float)v1 - (float) v2 << std::endl;
  330 + std::cout << "failed" << std::endl;
  331 + }
  332 +// std::cout << e[idx].U.len() << " " << mags[idx][0] << std::endl;
  333 +// std::cout << v2 << std::endl;
  334 + return(v1);
  335 +// return (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat);
  336 + // (
224 } 337 }
225 338
226 /// Returns the magnitude at the given index 339 /// Returns the magnitude at the given index
@@ -233,7 +346,7 @@ class cylinder @@ -233,7 +346,7 @@ class cylinder
233 /// Adds a new magnitude value to all points 346 /// Adds a new magnitude value to all points
234 /// @param m is the starting value for the new magnitude 347 /// @param m is the starting value for the new magnitude
235 void add_mag(T m = 0){ 348 void add_mag(T m = 0){
236 - for(unsigned int p = 0; p < e.size(); p++) 349 + for(unsigned int p = 0; p < N; p++)
237 mags[p].push_back(m); 350 mags[p].push_back(m);
238 } 351 }
239 352
@@ -279,8 +392,8 @@ class cylinder @@ -279,8 +392,8 @@ class cylinder
279 getPoints(int sides) 392 getPoints(int sides)
280 { 393 {
281 std::vector<std::vector <vec3<T> > > points; 394 std::vector<std::vector <vec3<T> > > points;
282 - points.resize(e.size());  
283 - for(int i = 0; i < e.size(); i++) 395 + points.resize(N);
  396 + for(int i = 0; i < N; i++)
284 { 397 {
285 points[i] = e[i].getPoints(sides); 398 points[i] = e[i].getPoints(sides);
286 } 399 }
@@ -319,7 +432,7 @@ class cylinder @@ -319,7 +432,7 @@ class cylinder
319 T len = L[0]; //allocate space for the segment length 432 T len = L[0]; //allocate space for the segment length
320 433
321 //for every consecutive point in the cylinder 434 //for every consecutive point in the cylinder
322 - for(unsigned p = 1; p < e.size(); p++){ 435 + for(unsigned p = 1; p < N; p++){
323 436
324 //p1 = pos[p]; //get the position and magnitude for the next point 437 //p1 = pos[p]; //get the position and magnitude for the next point
325 m1 = mags[p][m]; 438 m1 = mags[p][m];
stim/visualization/glObj.h
@@ -35,9 +35,9 @@ private: @@ -35,9 +35,9 @@ private:
35 glListBase(dList); 35 glListBase(dList);
36 36
37 glMatrixMode(GL_PROJECTION); 37 glMatrixMode(GL_PROJECTION);
38 - glLoadIdentity; 38 + glLoadIdentity();
39 glMatrixMode(GL_MODELVIEW); 39 glMatrixMode(GL_MODELVIEW);
40 - glLoadIdentity; 40 + glLoadIdentity();
41 41
42 } 42 }
43 43
@@ -129,7 +129,7 @@ public: @@ -129,7 +129,7 @@ public:
129 } 129 }
130 130
131 void 131 void
132 - RenderLine(std::vector< stim::vec<T> > l) 132 + RenderLine(std::vector< stim::vec3<T> > l)
133 { 133 {
134 glColor3f(0.5, 1.0, 0.5); 134 glColor3f(0.5, 1.0, 0.5);
135 glLineWidth(3.0); 135 glLineWidth(3.0);
stim/visualization/obj.h
@@ -322,13 +322,6 @@ protected: @@ -322,13 +322,6 @@ protected:
322 return OBJ_INVALID; 322 return OBJ_INVALID;
323 } 323 }
324 324
325 - //private method returns the vertex indices for a line  
326 - std::vector<unsigned> get_l_v(unsigned l){  
327 -  
328 -  
329 -  
330 - }  
331 -  
332 public: 325 public:
333 /// Constructor loads a Wavefront OBJ file 326 /// Constructor loads a Wavefront OBJ file
334 obj(std::string filename){ 327 obj(std::string filename){