codebase / stimlib

Commit ff3296427a3c1b3f9a1725a9a83dba900cb100da

Authored by Pavel Govyadinov
1 parent 029c021a

Stable fixes to the base classes. All working. No error here, I think

Showing 7 changed files with 155 additions and 39 deletions   Show diff stats
Inline Side-by-side
stim/cuda/branch_detection.cuh
  Show/Hide comments   View file @ff32964
... ... @@ -19,7 +19,7 @@ find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y, i
19 19 float* cpuCenters = (float*) malloc(x*y*sizeof(float));
20 20 int sizek = 7;
21 21  
22   - stringstream name;
  22 + std::stringstream name;
23 23  
24 24  
25 25 cpuCenters = stim::cuda::get_centers(texbufferID, texType, x, y, sizek, sigma, conn, threshold, iter);
... ...
stim/cuda/filter.cuh
  Show/Hide comments   View file @ff32964
... ... @@ -11,7 +11,7 @@
11 11 #include <stim/cuda/cudatools/devices.h>
12 12 #include <stim/cuda/cudatools/threads.h>
13 13 #include <stim/cuda/cuda_texture.cuh>
14   -#include <stim/cuda/ivote.cuh>
  14 +#include <stim/iVote/ivote2/local_max.cuh>
15 15 #include <stim/cuda/arraymath.cuh>
16 16  
17 17 #define IMAD(a,b,c) ( __mul24((a), (b)) + (c) )
... ... @@ -180,7 +180,7 @@ namespace stim
180 180 #endif
181 181  
182 182  
183   - stim::cuda::gpu_local_max<float>(centers, res, threshold, conn, DIM_X, DIM_Y);
  183 + stim::cuda::gpu_local_max<float>(centers, res, conn, DIM_X, DIM_Y);
184 184  
185 185 cudaDeviceSynchronize();
186 186  
... ...
stim/gl/gl_spider.h
  Show/Hide comments   View file @ff32964
... ... @@ -20,12 +20,11 @@
20 20 #include <stim/cuda/arraymath.cuh>
21 21 #include <stim/cuda/cuda_texture.cuh>
22 22 #include <stim/cuda/cudatools.h>
23   -#include <stim/cuda/ivote.cuh>
24 23 #include <stim/visualization/glObj.h>
25 24 #include <vector>
26 25 #include <stack>
27 26 #include <stim/cuda/branch_detection.cuh>
28   -#include "../../../volume-spider/glnetwork.h"
  27 +#include <stim/visualization/gl_network.h>
29 28 #include <stim/visualization/cylinder.h>
30 29 #include <iostream>
31 30 #include <fstream>
... ... @@ -72,7 +71,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
72 71 std::vector<float> mV; //A list of all the size vectors.
73 72 std::vector<float> lV; //A list of all the size vectors.
74 73  
75   - stim::matrix<float, 4> cT; //current Transformation matrix (tissue)->(texture)
  74 + stim::matrix<float> cT; //new matrix class //current Transformation matrix (tissue)->(texture)
76 75 GLuint texID; //OpenGL ID for the texture to be traced
77 76 stim::vec3<float> S; //Size of a voxel in the volume.
78 77 stim::vec3<float> R; //Dimensions of the volume.
... ... @@ -707,7 +706,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
707 706 m);
708 707 //get and store the current transformation matrix for later use.
709 708 glGetFloatv(GL_TEXTURE_MATRIX, curTrans);
710   - cT.set(curTrans);
  709 + cT = stim::matrix<float>(4,4, &curTrans[0]); //new matrix
711 710  
712 711 CHECK_OPENGL_ERROR
713 712 //revert back to default gl mode.
... ... @@ -960,6 +959,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
960 959 iter_dir = 0;
961 960 iter_siz = 0;
962 961 #endif
  962 + cT = stim::matrix<float>(4,4); //new matrix class.
963 963 stepsize = 6.0;
964 964 n_pixels = 16.0;
965 965  
... ... @@ -1301,10 +1301,10 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
1301 1301 setSeeds(std::string file)
1302 1302 {
1303 1303 std::ifstream myfile(file.c_str()); ///open a stream
1304   - string line;
  1304 + std::string line;
1305 1305 if(myfile.is_open())
1306 1306 {
1307   - while (getline(myfile, line))
  1307 + while (std::getline(myfile, line))
1308 1308 {
1309 1309 float x, y, z, u, v, w, m; ///read the xyz uvw and m coordinates.
1310 1310 myfile >> x >> y >> z >> u >> v >> w >> m;
... ... @@ -1319,24 +1319,25 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
1319 1319 }
1320 1320  
1321 1321 ///Saves the network to a file.
1322   - void
1323   - saveNetwork(std::string name)
1324   - {
1325   - stim::glObj<float> sk1;
1326   - for(int i = 0; i < nt.sizeE(); i++)
1327   - {
1328   - std::vector<float> cm = nt.getEdgeCenterLineMag(i);
1329   - std::vector<stim::vec3< float > > ce = nt.getEdgeCenterLine(i);
1330   - sk1.Begin(stim::OBJ_LINE);
1331   - for(int j = 0; j < ce.size(); j++)
1332   - {
1333   - sk1.TexCoord(cm[j]);
1334   - sk1.Vertex(ce[j][0], ce[j][1], ce[j][2]);
1335   - }
1336   - sk1.End();
1337   - }
1338   - sk1.save(name);
1339   - }
  1322 + void
  1323 + saveNetwork(std::string name, int xoffset = 0, int yoffset = 0, int zoffset = 0)
  1324 + {
  1325 + stim::glObj<float> sk1;
  1326 + for(int i = 0; i < nt.sizeE(); i++)
  1327 + {
  1328 + std::vector<float> cm = nt.getEdgeCenterLineMag(i);
  1329 + std::vector<stim::vec3< float > > ce = nt.getEdgeCenterLine(i);
  1330 + sk1.Begin(stim::OBJ_LINE);
  1331 + for(int j = 0; j < ce.size(); j++)
  1332 + {
  1333 + sk1.TexCoord(cm[j]);
  1334 + sk1.Vertex(ce[j][0]+xoffset, ce[j][1]+yoffset,ce[j][2]+zoffset);
  1335 + }
  1336 + sk1.End();
  1337 + }
  1338 + sk1.save(name);
  1339 + }
  1340 +
1340 1341  
1341 1342 ///Depreciated, but might be reused later()
1342 1343 ///returns a COPY of the entire stim::glObj object.
... ... @@ -1415,7 +1416,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
1415 1416 void
1416 1417 printTransform()
1417 1418 {
1418   - std::cout << cT << std::endl;
  1419 + std::cout << cT.toStr() << std::endl; ///new matrix changes
1419 1420 }
1420 1421  
1421 1422  
... ...
stim/math/matrix.h
  Show/Hide comments   View file @ff32964
... ... @@ -40,7 +40,7 @@ namespace stim{
40 40 public:
41 41 /// Constructor opens a mat4 file for writing
42 42 mat4file(std::string filename) {
43   - matfile.open(filename, std::ios::binary);
  43 + matfile.open(filename.c_str(), std::ios::binary);
44 44 }
45 45  
46 46 bool is_open() {
... ... @@ -281,6 +281,26 @@ public:
281 281 return result;
282 282 }
283 283  
  284 + ///modified for gl_spider.
  285 + //vector multiplication, multiplies a stim::vec<T> by stim::matrix<T>
  286 + vec<T> operator*(const vec<T> rhs) const {
  287 + if(rhs.size() != C)
  288 + {
  289 + std::cerr << "size of vector != number of columns" << std::endl;
  290 + return NULL;
  291 + }
  292 + vec<T> result(R);
  293 + for(int row = 0; row < R; row++)
  294 + {
  295 + for(int col = 0; col < C; col++)
  296 + {
  297 + result[row] += get(row,col)*rhs[row];
  298 + }
  299 + }
  300 +
  301 + return result;
  302 + }
  303 +
284 304 //returns a pointer to the raw matrix data (in column major format)
285 305 T* data(){
286 306 return M;
... ...
stim/math/matrix_sq.h
  Show/Hide comments   View file @ff32964
... ... @@ -75,6 +75,15 @@ struct matrix_sq
75 75 return *this;
76 76 }
77 77  
  78 + // M/rhs
  79 + CUDA_CALLABLE matrix_sq<T, N> operator/(T rhs)
  80 + {
  81 + for(int i=0; i<N; i++)
  82 + for(int j=0 ; j<N; j++)
  83 + M[i*N+j] = M[i*N+j] / rhs;
  84 + return *this;
  85 + }
  86 +
78 87 template<typename Y>
79 88 vec<Y> operator*(vec<Y> rhs){
80 89 unsigned int M = rhs.size();
... ... @@ -116,6 +125,22 @@ struct matrix_sq
116 125 return ss.str();
117 126 }
118 127  
  128 + std::string toTensor()
  129 + {
  130 + std::stringstream ss;
  131 +
  132 + ss << (*this)(0, 0) << " ";
  133 + ss << (*this)(0, 1) << " ";
  134 + ss << (*this)(1, 1) << " ";
  135 + ss << (*this)(2, 0) << " ";
  136 + ss << (*this)(2, 1) << " ";
  137 + ss << (*this)(2, 2);
  138 +
  139 + return ss.str();
  140 + }
  141 +
  142 +
  143 +
119 144 static matrix_sq<T, N> identity() {
120 145 matrix_sq<T, N> I;
121 146 I = 0;
... ...
stim/math/spharmonics.h
  Show/Hide comments   View file @ff32964
... ... @@ -3,6 +3,7 @@
3 3  
4 4 #include <complex>
5 5 #include <boost/math/special_functions/spherical_harmonic.hpp>
  6 +#include <stim/visualization/colormap.h>
6 7 #include <stim/math/constants.h>
7 8 #include <stim/math/random.h>
8 9 #include <vector>
... ... @@ -181,23 +182,23 @@ namespace stim {
181 182 /// @param n is the number of points of the surface of the sphere used to create the PDF. DEFAULT 1000
182 183 /// @param norm, a boolean that sets where the output vectors will be normalized between 0 and 1.
183 184 /// @param
184   - /*void pdf(std::vector<stim::vec3<T> > sph_pts, unsigned int l, int m, stim::vec3<T> c = stim::vec3<T>(0, 0, 0), unsigned int n = 1000, bool norm = true, std::vector<T> w = std::vector<T>())
  185 + void pdf(std::vector<stim::vec3<T> > sph_pts, unsigned int l, int m, stim::vec3<T> c = stim::vec3<T>(0, 0, 0), unsigned int n = 1000, bool norm = true, std::vector<T> w = std::vector<T>(), int normfactor = 2)
185 186 {
186 187 std::vector<double> weights; ///the weight at each point on the surface of the sphere.
187 188 // weights.resize(n);
188   - unsigned int nP = sph_pts.size();
189   - std::vector<stim::vec3<T> > sphere = stim::Random<T>::sample_sphere(n, 1.0, stim::TAU);
  189 + unsigned int nP = sph_pts.size(); ///sph_pts is the vectors we want to fit a spehrical harmonic to.
  190 + std::vector<stim::vec3<T> > sphere = stim::Random<T>::sample_sphere(n, 1.0, stim::TAU); ///randomsample a sphere of radius 1 and return the sampled vectors.
190 191 if (w.size() < nP)
191   - w = std::vector<T>(nP, 1.0);
  192 + w = std::vector<T>(nP, 1.0); //1 weight per each m.
192 193  
193   - for (int i = 0; i < n; i++)
  194 + for (int i = 0; i < n; i++) //for each weight
194 195 {
195 196 T val = 0;
196   - for (int j = 0; j < nP; j++)
  197 + for (int j = 0; j < nP; j++) //for each vector
197 198 {
198 199 stim::vec3<T> temp = sph_pts[j] - c;
199 200 if (temp.dot(sphere[i]) > 0)
200   - val += pow(temp.dot(sphere[i]), 4)*w[j];
  201 + val += pow(temp.dot(sphere[i]), normfactor)*w[j];
201 202 }
202 203 weights.push_back(val);
203 204 }
... ... @@ -223,7 +224,7 @@ namespace stim {
223 224 }
224 225 }
225 226 mcEnd();
226   - }*/
  227 + }
227 228  
228 229 std::string str() {
229 230  
... ...
stim/visualization/gl_spharmonics.h
  Show/Hide comments   View file @ff32964
... ... @@ -51,8 +51,7 @@ namespace stim {
51 51 magnitude = true;
52 52 }
53 53  
54   - gl_spharmonics(stim::spharmonics<T> disp, stim::spharmonics<T> color, size_t slices) :
55   - gl_spharmonics<T>(slices)
  54 + gl_spharmonics(stim::spharmonics<T> disp, stim::spharmonics<T> color, size_t slices) : gl_spharmonics<T>(slices)
56 55 {
57 56 Sc = color;
58 57 Sd = disp;
... ... @@ -62,6 +61,76 @@ namespace stim {
62 61 if (dlist) glDeleteLists(dlist, 1); //delete the display list when the object is destroyed
63 62 }
64 63  
  64 + void save(std::string filename)
  65 + {
  66 +
  67 + if (!tex) {
  68 + init_tex();
  69 + }
  70 +
  71 + stim::obj<T> object;
  72 + size_t theta_i, phi_i;
  73 + T d_theta = (T)stim::TAU / (T)N;
  74 + T d_phi = (T)stim::PI / (T)(N-1);
  75 + object.matKd("sfunc.jpg");
  76 + for (phi_i = 1; phi_i < N; phi_i++) {
  77 + T phi = phi_i * d_phi;
  78 + object.Begin(OBJ_TRIANGLE_STRIP);
  79 + for (theta_i = 0; theta_i <= N; theta_i++) {
  80 + T theta = (N - theta_i) * d_theta;
  81 + float theta_t = 1 - (float)theta_i / (float)N;
  82 +
  83 + T r;
  84 + if (!displacement) r = 1; //if no displacement, set the r value to 1 (renders a unit sphere)
  85 + else r = Sd(theta, phi); //otherwise calculate the displacement value
  86 +
  87 + glColor3f(1.0f, 1.0f, 1.0f);
  88 + if (!colormap) { //if no colormap is being rendered
  89 + if (r < 0) glColor3f(1.0, 0.0, 0.0); //if r is negative, render it red
  90 + else glColor3f(0.0, 1.0, 0.0); //otherwise render in green
  91 + }
  92 + if (magnitude) { //if the magnitude is being displaced, calculate the magnitude of r
  93 + if (r < 0) r = -r;
  94 + }
  95 + stim::vec3<T> s(r, theta, phi);
  96 + stim::vec3<T> c = s.sph2cart();
  97 + stim::vec3<T> n; //allocate a value to store the normal
  98 + if (!displacement) n = c; //if there is no displacement, the normal is spherical
  99 + else n = Sd.dphi(theta, phi).cross(Sd.dtheta(theta, phi)); //otherwise calculate the normal as the cross product of derivatives
  100 +
  101 + object.TexCoord(theta_t, (float)phi_i / (float)N);
  102 + //std::cout << theta_t <<" "<<(float)phi_i / (float)N << "--- -------------";
  103 + object.Normal(n[0], n[1], n[2]);
  104 + object.Vertex(c[0], c[1], c[2]);
  105 +
  106 + T r1;
  107 + if (!displacement) r1 = 1;
  108 + else r1 = Sd(theta, phi - d_phi);
  109 + if (!colormap) { //if no colormap is being rendered
  110 + if (r1 < 0) glColor3f(1.0, 0.0, 0.0); //if r1 is negative, render it red
  111 + else glColor3f(0.0, 1.0, 0.0); //otherwise render in green
  112 + }
  113 + if (magnitude) { //if the magnitude is being rendered, calculate the magnitude of r
  114 + if (r1 < 0) r1 = -r1;
  115 + }
  116 + stim::vec3<T> s1(r1, theta, phi - d_phi);
  117 + stim::vec3<T> c1 = s1.sph2cart();
  118 + stim::vec3<T> n1;
  119 + if (!displacement) n1 = c1;
  120 + else n1 = Sd.dphi(theta, phi - d_phi).cross(Sd.dtheta(theta, phi - d_phi));
  121 +
  122 + //std::cout << theta_t << " " << (float)(phi_i - 1) / (float) N << std::endl;
  123 + object.TexCoord(theta_t, 1.0/(2*(N)) + (float)(phi_i-1) / (float)N);
  124 + object.Normal(n1[0], n1[1], n1[2]);
  125 + object.Vertex(c1[0], c1[1], c1[2]);
  126 + }
  127 + object.End();
  128 + }
  129 + object.matKd();
  130 + object.save(filename);
  131 + }
  132 +
  133 +
65 134 /// This function renders the spherical harmonic to the current OpenGL context
66 135 void render() {
67 136 //glShadeModel(GL_FLAT);
... ...