codebase / stimlib

Commit 487a9b49549dba99d44b68b7fcfeb4cd521d11bb

Authored by David Mayerich
1 parent eaf2cb07

added the ability to load OBJ files into spherical harmonic functions

Showing 7 changed files with 558 additions and 458 deletions   Show diff stats
Inline Side-by-side
stim/math/vector.h renamed to stim/math/mathvec.h
  Show/Hide comments   View file @487a9b4
... ... @@ -4,7 +4,7 @@
4 4 #include <iostream>
5 5 #include <cmath>
6 6 #include <sstream>
7   -//#include "rts/math/point.h"
  7 +#include <vector>
8 8 #include "../cuda/callable.h"
9 9  
10 10  
... ... @@ -13,55 +13,53 @@ namespace stim
13 13  
14 14  
15 15  
16   -template <class T, int N=3>
17   -struct vec
  16 +template <class T>
  17 +struct vec : public std::vector<T>
18 18 {
19   - T v[N];
  19 + CUDA_CALLABLE vec(){
20 20  
21   - CUDA_CALLABLE vec()
22   - {
23   - //memset(v, 0, sizeof(T) * N);
24   - for(int i=0; i<N; i++)
25   - v[i] = 0;
26 21 }
27 22  
28   - //efficiency constructor, makes construction easier for 1D-4D vectors
  23 + //efficiency constructors, makes construction easier for 1D-4D vectors
29 24 CUDA_CALLABLE vec(T rhs)
30 25 {
31   - for(int i=0; i<N; i++)
32   - v[i] = rhs;
  26 + resize(1, rhs);
  27 + //for(int i=0; i<N; i++)
  28 + // v[i] = rhs;
33 29 }
34 30 CUDA_CALLABLE vec(T x, T y)
35 31 {
36   - v[0] = x;
37   - v[1] = y;
  32 + push_back(x);
  33 + push_back(y);
38 34 }
39 35 CUDA_CALLABLE vec(T x, T y, T z)
40 36 {
41   - v[0] = x;
42   - v[1] = y;
43   - v[2] = z;
  37 + push_back(x);
  38 + push_back(y);
  39 + push_back(z);
44 40 }
45 41 CUDA_CALLABLE vec(T x, T y, T z, T w)
46 42 {
47   - v[0] = x;
48   - v[1] = y;
49   - v[2] = z;
50   - v[3] = w;
  43 + push_back(x);
  44 + push_back(y);
  45 + push_back(z);
  46 + push_back(w);
51 47 }
52 48  
53 49 //copy constructor
54   - CUDA_CALLABLE vec( const vec<T, N>& other){
  50 + CUDA_CALLABLE vec( const vec<T>& other){
  51 + unsigned int N = other.size();
55 52 for(int i=0; i<N; i++)
56   - v[i] = other.v[i];
  53 + push_back(other[i]);
57 54 }
58 55  
59 56  
60 57 template< typename U >
61   - CUDA_CALLABLE operator vec<U, N>(){
62   - vec<U, N> result;
  58 + CUDA_CALLABLE operator vec<U>(){
  59 + unsigned int N = size();
  60 + vec<U> result;
63 61 for(int i=0; i<N; i++)
64   - result.v[i] = v[i];
  62 + result.push_back(at(i));
65 63  
66 64 return result;
67 65 }
... ... @@ -71,46 +69,54 @@ struct vec
71 69  
72 70 CUDA_CALLABLE T len() const
73 71 {
  72 + unsigned int N = size();
  73 +
74 74 //compute and return the vector length
75 75 T sum_sq = (T)0;
76 76 for(int i=0; i<N; i++)
77 77 {
78   - sum_sq += v[i] * v[i];
  78 + sum_sq += pow( at(i), 2 );
79 79 }
80 80 return sqrt(sum_sq);
81 81  
82 82 }
83 83  
84   - CUDA_CALLABLE vec<T, N> cart2sph() const
  84 + CUDA_CALLABLE vec<T> cart2sph() const
85 85 {
86 86 //convert the vector from cartesian to spherical coordinates
87 87 //x, y, z -> r, theta, phi (where theta = 0 to 2*pi)
88 88  
89   - vec<T, N> sph;
90   - sph[0] = std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
91   - sph[1] = std::atan2(v[1], v[0]);
92   - sph[2] = std::acos(v[2] / sph[0]);
  89 + vec<T> sph;
  90 + sph.push_back(std::sqrt(at(0)*at(0) + at(1)*at(1) + at(2)*at(2)));
  91 + sph.push_back(std::atan2(at(1), at(0)));
  92 +
  93 + if(sph[0] == 0)
  94 + sph.push_back(0);
  95 + else
  96 + sph.push_back(std::acos(at(2) / sph[0]));
93 97  
94 98 return sph;
95 99 }
96 100  
97   - CUDA_CALLABLE vec<T, N> sph2cart() const
  101 + CUDA_CALLABLE vec<T> sph2cart() const
98 102 {
99 103 //convert the vector from cartesian to spherical coordinates
100 104 //r, theta, phi -> x, y, z (where theta = 0 to 2*pi)
101 105  
102   - vec<T, N> cart;
103   - cart[0] = v[0] * std::cos(v[1]) * std::sin(v[2]);
104   - cart[1] = v[0] * std::sin(v[1]) * std::sin(v[2]);
105   - cart[2] = v[0] * std::cos(v[2]);
  106 + vec<T> cart;
  107 + cart.push_back(at(0) * std::cos(at(1)) * std::sin(at(2)));
  108 + cart.push_back(at(0) * std::sin(at(1)) * std::sin(at(2)));
  109 + cart.push_back(at(0) * std::cos(at(2)));
106 110  
107 111 return cart;
108 112 }
109 113  
110   - CUDA_CALLABLE vec<T, N> norm() const
  114 + CUDA_CALLABLE vec<T> norm() const
111 115 {
112   - //compute and return the vector norm
113   - vec<T, N> result;
  116 + unsigned int N = size();
  117 +
  118 + //compute and return the unit vector
  119 + vec<T> result;
114 120  
115 121 //compute the vector length
116 122 T l = len();
... ... @@ -118,125 +124,138 @@ struct vec
118 124 //normalize
119 125 for(int i=0; i<N; i++)
120 126 {
121   - result.v[i] = v[i] / l;
  127 + result.push_back(at(i) / l);
122 128 }
123 129  
124 130 return result;
125 131 }
126 132  
127   - CUDA_CALLABLE vec<T, 3> cross(const vec<T, 3> rhs) const
  133 + CUDA_CALLABLE vec<T> cross(const vec<T> rhs) const
128 134 {
129   - vec<T, 3> result;
  135 + vec<T> result;
130 136  
131 137 //compute the cross product (only valid for 3D vectors)
132   - result[0] = v[1] * rhs.v[2] - v[2] * rhs.v[1];
133   - result[1] = v[2] * rhs.v[0] - v[0] * rhs.v[2];
134   - result[2] = v[0] * rhs.v[1] - v[1] * rhs.v[0];
  138 + result.push_back(at(1) * rhs[2] - at(2) * rhs[1]);
  139 + result.push_back(at(2) * rhs[0] - at(0) * rhs[2]);
  140 + result.push_back(at(0) * rhs[1] - at(1) * rhs[0]);
135 141  
136 142 return result;
137 143 }
138 144  
139   - CUDA_CALLABLE T dot(vec<T, N> rhs) const
  145 + CUDA_CALLABLE T dot(vec<T> rhs) const
140 146 {
141 147 T result = (T)0;
142   -
  148 + unsigned int N = size();
143 149 for(int i=0; i<N; i++)
144   - result += v[i] * rhs.v[i];
  150 + result += at(i) * rhs[i];
145 151  
146 152 return result;
147 153  
148 154 }
149 155  
150 156 //arithmetic
151   - CUDA_CALLABLE vec<T, N> operator+(vec<T, N> rhs) const
  157 + CUDA_CALLABLE vec<T> operator+(vec<T> rhs) const
152 158 {
153   - vec<T, N> result;
  159 + vec<T> result;
154 160  
  161 + unsigned int N = size();
155 162 for(int i=0; i<N; i++)
156   - result.v[i] = v[i] + rhs.v[i];
  163 + result.push_back(at(i) + rhs[i]);
157 164  
158 165 return result;
159 166 }
160   - CUDA_CALLABLE vec<T, N> operator+(T rhs) const
  167 + CUDA_CALLABLE vec<T> operator+(T rhs) const
161 168 {
162   - vec<T, N> result;
  169 + vec<T> result;
163 170  
164 171 for(int i=0; i<N; i++)
165   - result.v[i] = v[i] + rhs;
  172 + result.push_back(at(i) + rhs);
166 173  
167 174 return result;
168 175 }
169   - CUDA_CALLABLE vec<T, N> operator-(vec<T, N> rhs) const
  176 + CUDA_CALLABLE vec<T> operator-(vec<T> rhs) const
170 177 {
171   - vec<T, N> result;
  178 + vec<T> result;
172 179  
  180 + unsigned int N = size();
173 181 for(int i=0; i<N; i++)
174   - result.v[i] = v[i] - rhs.v[i];
  182 + result.push_back(at(i) - rhs[i]);
175 183  
176 184 return result;
177 185 }
178   - CUDA_CALLABLE vec<T, N> operator*(T rhs) const
  186 + CUDA_CALLABLE vec<T> operator*(T rhs) const
179 187 {
180   - vec<T, N> result;
  188 + vec<T> result;
181 189  
  190 + unsigned int N = size();
182 191 for(int i=0; i<N; i++)
183   - result.v[i] = v[i] * rhs;
  192 + result.push_back(at(i) * rhs);
184 193  
185 194 return result;
186 195 }
187   - CUDA_CALLABLE vec<T, N> operator/(T rhs) const
  196 + CUDA_CALLABLE vec<T> operator/(T rhs) const
188 197 {
189   - vec<T, N> result;
  198 + vec<T> result;
190 199  
  200 + unsigned int N = size();
191 201 for(int i=0; i<N; i++)
192   - result.v[i] = v[i] / rhs;
  202 + result.push_back(at(i) / rhs);
193 203  
194 204 return result;
195 205 }
196   - CUDA_CALLABLE vec<T, N> operator*=(T rhs){
  206 + CUDA_CALLABLE vec<T> operator*=(T rhs){
  207 +
  208 + unsigned int N = size();
  209 + for(int i=0; i<N; i++)
  210 + at(i) = at(i) * rhs;
  211 + return *this;
  212 + }
  213 + CUDA_CALLABLE vec<T> operator+=(vec<T> rhs){
  214 + unsigned int N = size();
197 215 for(int i=0; i<N; i++)
198   - v[i] = v[i] * rhs;
  216 + at(i) += rhs[i];
199 217 return *this;
200 218 }
201   - CUDA_CALLABLE vec<T, N> & operator=(T rhs){
  219 + CUDA_CALLABLE vec<T> & operator=(T rhs){
  220 +
  221 + unsigned int N = size();
202 222 for(int i=0; i<N; i++)
203   - v[i] = rhs;
  223 + at(i) = rhs;
204 224 return *this;
205 225 }
206 226  
207 227 template<typename Y>
208   - CUDA_CALLABLE vec<T, N> & operator=(vec<Y, N> rhs){
  228 + CUDA_CALLABLE vec<T> & operator=(vec<Y> rhs){
  229 + unsigned int N = rhs.size();
  230 + resize(N);
  231 +
209 232 for(int i=0; i<N; i++)
210   - v[i] = rhs.v[i];
  233 + at(i) = rhs[i];
211 234 return *this;
212 235 }
213 236 //unary minus
214   - CUDA_CALLABLE vec<T, N> operator-() const{
215   - vec<T, N> r;
  237 + CUDA_CALLABLE vec<T> operator-() const{
  238 + vec<T> r;
216 239  
217 240 //negate the vector
  241 + unsigned int N = size();
218 242 for(int i=0; i<N; i++)
219   - r.v[i] = -v[i];
  243 + r.push_back(-at(i));
220 244  
221 245 return r;
222 246 }
223 247  
224   - CUDA_CALLABLE bool operator==(vec<T, N> rhs) const
225   - {
226   - if ( (rhs.v[0] == v[0]) && (rhs.v[1] == v[1]) && (rhs.v[2] == v[2]) )
227   - return true;
228   -
229   - return false;
230   - }
231 248  
232 249 std::string str() const
233 250 {
234 251 std::stringstream ss;
235 252  
  253 + unsigned int N = size();
  254 +
236 255 ss<<"[";
237 256 for(int i=0; i<N; i++)
238 257 {
239   - ss<<v[i];
  258 + ss<<at(i);
240 259 if(i != N-1)
241 260 ss<<", ";
242 261 }
... ... @@ -245,19 +264,13 @@ struct vec
245 264 return ss.str();
246 265 }
247 266  
248   - //bracket operator - allows assignment to the vector
249   - CUDA_CALLABLE T& operator[](const unsigned int i)
250   - {
251   - return v[i];
252   - }
253   -
254 267 };
255 268  
256 269  
257 270 } //end namespace rts
258 271  
259   -template <typename T, int N>
260   -std::ostream& operator<<(std::ostream& os, stim::vec<T, N> v)
  272 +template <typename T>
  273 +std::ostream& operator<<(std::ostream& os, stim::vec<T> v)
261 274 {
262 275 os<<v.str();
263 276 return os;
... ... @@ -265,10 +278,10 @@ std::ostream&amp; operator&lt;&lt;(std::ostream&amp; os, stim::vec&lt;T, N&gt; v)
265 278  
266 279  
267 280  
268   -template <typename T, int N>
269   -CUDA_CALLABLE stim::vec<T, N> operator*(T lhs, stim::vec<T, N> rhs)
  281 +template <typename T>
  282 +CUDA_CALLABLE stim::vec<T> operator*(T lhs, stim::vec<T> rhs)
270 283 {
271   - stim::vec<T, N> r;
  284 + stim::vec<T> r;
272 285  
273 286 return rhs * lhs;
274 287 }
... ...
stim/math/matrix.h
  Show/Hide comments   View file @487a9b4
... ... @@ -4,7 +4,7 @@
4 4 //#include "rts/vector.h"
5 5 #include <string.h>
6 6 #include <iostream>
7   -#include "vector.h"
  7 +#include "mathvec.h"
8 8 #include "../cuda/callable.h"
9 9  
10 10 namespace stim{
... ... @@ -40,9 +40,12 @@ struct matrix
40 40 }
41 41  
42 42 template<typename Y>
43   - CUDA_CALLABLE vec<Y, N> operator*(vec<Y, N> rhs)
  43 + CUDA_CALLABLE vec<Y> operator*(vec<Y> rhs)
44 44 {
45   - vec<Y, N> result;
  45 + unsigned int N = rhs.size();
  46 +
  47 + vec<Y> result;
  48 + result.resize(N);
46 49  
47 50 for(int r=0; r<N; r++)
48 51 for(int c=0; c<N; c++)
... ...
stim/math/quaternion.h
  Show/Hide comments   View file @487a9b4
... ... @@ -26,13 +26,13 @@ public:
26 26  
27 27 CUDA_CALLABLE void CreateRotation(T theta, T ux, T uy, T uz){
28 28  
29   - vec<T, 3> u(ux, uy, uz);
  29 + vec<T> u(ux, uy, uz);
30 30 CreateRotation(theta, u);
31 31 }
32 32  
33   - CUDA_CALLABLE void CreateRotation(T theta, vec<T, 3> u){
  33 + CUDA_CALLABLE void CreateRotation(T theta, vec<T> u){
34 34  
35   - vec<T, 3> u_hat = u.norm();
  35 + vec<T> u_hat = u.norm();
36 36  
37 37 //assign the given Euler rotation to this quaternion
38 38 w = (T)cos(theta/2);
... ... @@ -41,7 +41,7 @@ public:
41 41 z = u_hat[2]*(T)sin(theta/2);
42 42 }
43 43  
44   - CUDA_CALLABLE void CreateRotation(vec<T, 3> from, vec<T, 3> to){
  44 + CUDA_CALLABLE void CreateRotation(vec<T> from, vec<T> to){
45 45  
46 46 vec<T> r = from.cross(to); //compute the rotation vector
47 47 T theta = asin(r.len()); //compute the angle of the rotation about r
... ...
stim/math/spharmonics.h 0 → 100644
  Show/Hide comments   View file @487a9b4
  1 +#ifndef STIM_SPH_HARMONICS
  2 +#define STIM_SPH_HARMONICS
  3 +
  4 +
  5 +#include <boost/math/special_functions/spherical_harmonic.hpp>
  6 +#include <vector>
  7 +
  8 +#define PI 3.14159
  9 +#define WIRE_SCALE 1.001
  10 +namespace stim{
  11 +
  12 +template<class T>
  13 +class spharmonics{
  14 +
  15 +protected:
  16 +
  17 + std::vector<T> C; //list of SH coefficients
  18 +
  19 + unsigned int mcN; //number of Monte-Carlo samples
  20 +
  21 + //calculate the value of the SH basis function (l, m) at (theta, phi)
  22 + //here, theta = [0, PI], phi = [0, 2*PI]
  23 + double SH(int l, int m, double theta, double phi){
  24 + return boost::math::spherical_harmonic_r(l, m, phi, theta);
  25 + }
  26 +
  27 + unsigned int coeff_1d(unsigned int l, int m){
  28 + return pow(l + 1, 2) - (l - m) - 1;
  29 + }
  30 +
  31 +
  32 +
  33 +
  34 +public:
  35 +
  36 + void push(double c){
  37 + C.push_back(c);
  38 + }
  39 +
  40 + void resize(unsigned int n){
  41 + C.resize(n);
  42 + }
  43 +
  44 + void setc(unsigned int l, int m, T value){
  45 + unsigned int c = coeff_1d(l, m);
  46 + C[c] = value;
  47 + }
  48 +
  49 + void setc(unsigned int c, T value){
  50 + C[c] = value;
  51 + }
  52 +
  53 + /// Initialize Monte-Carlo sampling of a function using N spherical harmonics coefficients
  54 +
  55 + /// @param N is the number of spherical harmonics coefficients used to represent the user function
  56 + void mcBegin(unsigned int coefficients){
  57 + C.resize(coefficients, 0);
  58 + mcN = 0;
  59 + }
  60 +
  61 + void mcBegin(unsigned int l, int m){
  62 + unsigned int c = pow(l + 1, 2) - (l - m);
  63 + mcBegin(c);
  64 + }
  65 +
  66 + void mcSample(double theta, double phi, double val){
  67 +
  68 + int l, m;
  69 + double sh;
  70 +
  71 + l = m = 0;
  72 + for(unsigned int i = 0; i < C.size(); i++){
  73 +
  74 + sh = SH(l, m, theta, phi);
  75 + C[i] += sh * val;
  76 +
  77 + m++; //increment m
  78 +
  79 + //if we're in a new tier, increment l and set m = -l
  80 + if(m > l){
  81 + l++;
  82 + m = -l;
  83 + }
  84 + } //end for all coefficients
  85 +
  86 + //increment the number of samples
  87 + mcN++;
  88 +
  89 + } //end mcSample()
  90 +
  91 + void mcEnd(){
  92 +
  93 + //divide all coefficients by the number of samples
  94 + for(unsigned int i = 0; i < C.size(); i++)
  95 + C[i] /= mcN;
  96 + }
  97 +
  98 + /// Generates a PDF describing the probability distribution of points on a spherical surface
  99 +
  100 + /// @param sph_pts is a list of points in spherical coordinates (theta, phi) where theta = [0, 2pi] and phi = [0, pi]
  101 + /// @param l is the maximum degree of the spherical harmonic function
  102 + /// @param m is the maximum order
  103 + void pdf(std::vector<stim::vec<double>> sph_pts, unsigned int l, int m){
  104 +
  105 + mcBegin( l, m ); //begin spherical harmonic sampling
  106 +
  107 + unsigned int nP = sph_pts.size();
  108 +
  109 + for(unsigned int p = 0; p < nP; p++){
  110 + mcSample(sph_pts[p][1], sph_pts[p][2], 1.0);
  111 + }
  112 +
  113 + mcEnd();
  114 + }
  115 +
  116 + std::string str(){
  117 +
  118 + std::stringstream ss;
  119 +
  120 + int l, m;
  121 + l = m = 0;
  122 + for(unsigned int i = 0; i < C.size(); i++){
  123 +
  124 + ss<<C[i]<<'\t';
  125 +
  126 + m++; //increment m
  127 +
  128 + //if we're in a new tier, increment l and set m = -l
  129 + if(m > l){
  130 + l++;
  131 + m = -l;
  132 +
  133 + ss<<std::endl;
  134 +
  135 + }
  136 + }
  137 +
  138 + return ss.str();
  139 +
  140 +
  141 + }
  142 +
  143 + /// Returns the value of the function at the coordinate (theta, phi)
  144 +
  145 + /// @param theta = [0, 2pi]
  146 + /// @param phi = [0, pi]
  147 + double operator()(double theta, double phi){
  148 +
  149 + double fx = 0;
  150 +
  151 + int l = 0;
  152 + int m = 0;
  153 + for(unsigned int i = 0; i < C.size(); i++){
  154 + fx += C[i] * SH(l, m, theta, phi);
  155 + m++;
  156 + if(m > l){
  157 + l++;
  158 + m = -l;
  159 + }
  160 +
  161 + }
  162 +
  163 + return fx;
  164 + }
  165 +
  166 +}; //end class sph_harmonics
  167 +
  168 +
  169 +
  170 +
  171 +}
  172 +
  173 +
  174 +#endif
... ...
stim/visualization/camera.h
  Show/Hide comments   View file @487a9b4
1   -#include "../math/vector.h"
  1 +#include "../math/mathvec.h"
2 2 #include "../math/quaternion.h"
3 3 #include "../math/matrix.h"
4 4  
... ... @@ -11,32 +11,32 @@ namespace stim{
11 11  
12 12 class camera
13 13 {
14   - vec<float, 3> d; //direction that the camera is pointing
15   - vec<float, 3> p; //position of the camera
16   - vec<float, 3> up; //"up" direction
  14 + vec<float> d; //direction that the camera is pointing
  15 + vec<float> p; //position of the camera
  16 + vec<float> up; //"up" direction
17 17 float focus; //focal length of the camera
18 18 float fov;
19 19  
20 20 //private function makes sure that the up vector is orthogonal to the direction vector and both are normalized
21 21 void stabalize()
22 22 {
23   - vec<float, 3> side = up.cross(d);
  23 + vec<float> side = up.cross(d);
24 24 up = d.cross(side);
25 25 up = up.norm();
26 26 d = d.norm();
27 27 }
28 28  
29 29 public:
30   - void setPosition(vec<float, 3> pos)
  30 + void setPosition(vec<float> pos)
31 31 {
32 32 p = pos;
33 33 }
34   - void setPosition(float x, float y, float z){setPosition(vec<float, 3>(x, y, z));}
  34 + void setPosition(float x, float y, float z){setPosition(vec<float>(x, y, z));}
35 35  
36 36 void setFocalDistance(float distance){focus = distance;}
37 37 void setFOV(float field_of_view){fov = field_of_view;}
38 38  
39   - void LookAt(vec<float, 3> pos)
  39 + void LookAt(vec<float> pos)
40 40 {
41 41 //find the new direction
42 42 d = pos - p;
... ... @@ -47,22 +47,22 @@ public:
47 47 //stabalize the camera
48 48 stabalize();
49 49 }
50   - void LookAt(float px, float py, float pz){LookAt(vec<float, 3>(px, py, pz));}
51   - void LookAt(vec<float, 3> pos, vec<float, 3> new_up){up = new_up; LookAt(pos);}
52   - void LookAt(float px, float py, float pz, float ux, float uy, float uz){LookAt(vec<float, 3>(px, py, pz), vec<float, 3>(ux, uy, uz));}
  50 + void LookAt(float px, float py, float pz){LookAt(vec<float>(px, py, pz));}
  51 + void LookAt(vec<float> pos, vec<float> new_up){up = new_up; LookAt(pos);}
  52 + void LookAt(float px, float py, float pz, float ux, float uy, float uz){LookAt(vec<float>(px, py, pz), vec<float>(ux, uy, uz));}
53 53 void LookAtDolly(float lx, float ly, float lz)
54 54 {
55 55 //find the current focus point
56   - vec<float, 3> f = p + focus*d;
57   - vec<float, 3> T = vec<float, 3>(lx, ly, lz) - f;
  56 + vec<float> f = p + focus*d;
  57 + vec<float> T = vec<float>(lx, ly, lz) - f;
58 58 p = p + T;
59 59 }
60 60  
61   - void Dolly(vec<float, 3> direction)
  61 + void Dolly(vec<float> direction)
62 62 {
63 63 p = p+direction;
64 64 }
65   - void Dolly(float x, float y, float z){Dolly(vec<float, 3>(x, y, z));}
  65 + void Dolly(float x, float y, float z){Dolly(vec<float>(x, y, z));}
66 66 void Push(float delta)
67 67 {
68 68 if(delta > focus)
... ... @@ -80,7 +80,7 @@ public:
80 80 qx.CreateRotation(theta_x, up[0], up[1], up[2]);
81 81  
82 82 //y rotation is around the side axis
83   - vec<float, 3> side = up.cross(d);
  83 + vec<float> side = up.cross(d);
84 84 quaternion<float> qy;
85 85 qy.CreateRotation(theta_y, side[0], side[1], side[2]);
86 86  
... ... @@ -118,28 +118,28 @@ public:
118 118 void OrbitFocus(float theta_x, float theta_y)
119 119 {
120 120 //find the focal point
121   - vec<float, 3> focal_point = p + focus*d;
  121 + vec<float> focal_point = p + focus*d;
122 122  
123 123 //center the coordinate system on the focal point
124   - vec<float, 3> centered = p - (focal_point - vec<float, 3>(0, 0, 0));
  124 + vec<float> centered = p - (focal_point - vec<float>(0, 0, 0));
125 125  
126 126 //create the x rotation (around the up vector)
127 127 quaternion<float> qx;
128 128 qx.CreateRotation(theta_x, up[0], up[1], up[2]);
129   - centered = vec<float, 3>(0, 0, 0) + qx.toMatrix3()*(centered - vec<float, 3>(0, 0, 0));
  129 + centered = vec<float>(0, 0, 0) + qx.toMatrix3()*(centered - vec<float>(0, 0, 0));
130 130  
131 131 //get a side vector for theta_y rotation
132   - vec<float, 3> side = up.cross((vec<float, 3>(0, 0, 0) - centered).norm());
  132 + vec<float> side = up.cross((vec<float>(0, 0, 0) - centered).norm());
133 133  
134 134 quaternion<float> qy;
135 135 qy.CreateRotation(theta_y, side[0], side[1], side[2]);
136   - centered = vec<float, 3>(0, 0, 0) + qy.toMatrix3()*(centered - vec<float, 3>(0, 0, 0));
  136 + centered = vec<float>(0, 0, 0) + qy.toMatrix3()*(centered - vec<float>(0, 0, 0));
137 137  
138 138 //perform the rotation on the centered camera position
139 139 //centered = final.toMatrix()*centered;
140 140  
141 141 //re-position the camera
142   - p = centered + (focal_point - vec<float, 3>(0, 0, 0));
  142 + p = centered + (focal_point - vec<float>(0, 0, 0));
143 143  
144 144 //make sure we are looking at the focal point
145 145 LookAt(focal_point);
... ... @@ -151,17 +151,17 @@ public:
151 151  
152 152 void Slide(float u, float v)
153 153 {
154   - vec<float, 3> V = up.norm();
155   - vec<float, 3> U = up.cross(d).norm();
  154 + vec<float> V = up.norm();
  155 + vec<float> U = up.cross(d).norm();
156 156  
157 157 p = p + (V * v) + (U * u);
158 158 }
159 159  
160 160 //accessor methods
161   - vec<float, 3> getPosition(){return p;}
162   - vec<float, 3> getUp(){return up;}
163   - vec<float, 3> getDirection(){return d;}
164   - vec<float, 3> getLookAt(){return p + focus*d;}
  161 + vec<float> getPosition(){return p;}
  162 + vec<float> getUp(){return up;}
  163 + vec<float> getDirection(){return d;}
  164 + vec<float> getLookAt(){return p + focus*d;}
165 165 float getFOV(){return fov;}
166 166  
167 167 //output the camera settings
... ... @@ -182,9 +182,9 @@ public:
182 182 //constructor
183 183 camera()
184 184 {
185   - p = vec<float, 3>(0, 0, 0);
186   - d = vec<float, 3>(0, 0, 1);
187   - up = vec<float, 3>(0, 1, 0);
  185 + p = vec<float>(0, 0, 0);
  186 + d = vec<float>(0, 0, 1);
  187 + up = vec<float>(0, 1, 0);
188 188 focus = 1;
189 189  
190 190 }
... ...
stim/visualization/spharmonics.h renamed to stim/visualization/gl_spharmonics.h
  Show/Hide comments   View file @487a9b4
1   -#ifndef STIM_SPH_HARMONICS
2   -#define STIM_SPH_HARMONICS
3   -
4   -#include <GL/gl.h>
5   -
6   -#include <stim/gl/error.h>
7   -#include <stim/visualization/colormap.h>
8   -#include <vector>
9   -
10   -#define PI 3.14159
11   -#define WIRE_SCALE 1.001
12   -namespace stim{
13   -
14   - class sph_harmonics{
15   -
16   - private:
17   -
18   - double* func; //stores the raw function data (samples at each point)
19   -
20   - GLuint color_tex; //texture map that acts as a colormap for the spherical function
21   -
22   - unsigned int N; //resolution of the spherical grid
23   -
24   - std::vector<double> C; //list of SH coefficients
25   -
26   -
27   - //evaluates an associated Legendre polynomial (-l <= m <= l)
28   - double P(int l,int m,double x)
29   - {
30   - // evaluate an Associated Legendre Polynomial P(l,m,x) at x
31   - double pmm = 1.0;
32   - if(m>0) {
33   - double somx2 = sqrt((1.0-x)*(1.0+x));
34   - double fact = 1.0;
35   - for(int i=1; i<=m; i++) {
36   - pmm *= (-fact) * somx2;
37   - fact += 2.0;
38   - }
39   - }
40   - if(l==m) return pmm;
41   - double pmmp1 = x * (2.0*m+1.0) * pmm;
42   - if(l==m+1) return pmmp1;
43   - double pll = 0.0;
44   - for(int ll=m+2; ll<=l; ++ll) {
45   - pll = ( (2.0*ll-1.0)*x*pmmp1-(ll+m-1.0)*pmm ) / (ll-m);
46   - pmm = pmmp1;
47   - pmmp1 = pll;
48   - }
49   - return pll;
50   - }
51   -
52   - //recursively calculate a factorial given a positive integer n
53   - unsigned int factorial(unsigned int n) {
54   - if (n == 0)
55   - return 1;
56   - return n * factorial(n - 1);
57   - }
58   -
59   - //calculate the SH scaling constant
60   - double K(int l, int m){
61   -
62   - // renormalisation constant for SH function
63   - double temp = ((2.0*l+1.0)*factorial(l-m)) / (4.0*PI*factorial(l+m));
64   - return sqrt(temp);
65   - }
66   -
67   - //calculate the value of the SH basis function (l, m) at (theta, phi)
68   - //here, theta = [0, PI], phi = [0, 2*PI]
69   - double SH(int l, int m, double theta, double phi){
70   - // return a point sample of a Spherical Harmonic basis function
71   - // l is the band, range [0..N]
72   - // m in the range [-l..l]
73   - // theta in the range [0..Pi]
74   - // phi in the range [0..2*Pi]
75   - const double sqrt2 = sqrt(2.0);
76   - if(m==0) return K(l,0)*P(l,m,cos(theta));
77   - else if(m>0) return sqrt2*K(l,m)*cos(m*phi)*P(l,m,cos(theta));
78   - else return sqrt2*K(l,-m)*sin(-m*phi)*P(l,-m,cos(theta));
79   - }
80   -
81   - void gen_function(){
82   -
83   - //initialize the function to zero
84   - memset(func, 0, sizeof(double) * N * N);
85   -
86   - double theta, phi;
87   - double result;
88   - int l, m;
89   -
90   - l = m = 0;
91   - for(unsigned int c = 0; c < C.size(); c++){
92   -
93   -
94   - for(unsigned int xi = 0; xi < N; xi++)
95   - for(unsigned int yi = 0; yi < N; yi++){
96   -
97   - theta = (2 * PI) * ((double)xi / (N-1));
98   - phi = PI * ((double)yi / (N-1));
99   - result = C[c] * SH(l, m, phi, theta); //phi and theta are reversed here (damn physicists)
100   - func[yi * N + xi] += result;
101   - }
102   -
103   - m++; //increment m
104   -
105   - //if we're in a new tier, increment l and set m = -l
106   - if(m > l){
107   - l++;
108   - m = -l;
109   - }
110   - }
111   - }
112   -
113   - void gl_prep_draw(){
114   -
115   - //enable depth testing
116   - //this has to be used instead of culling because the sphere can have negative values
117   - glEnable(GL_DEPTH_TEST);
118   - glDepthMask(GL_TRUE);
119   - glEnable(GL_TEXTURE_2D); //enable 2D texture mapping
120   - }
121   -
122   - //draw a texture mapped sphere representing the function surface
123   - void gl_draw_sphere() {
124   -
125   - //PI is used to convert from spherical to cartesian coordinates
126   - //const double PI = 3.14159;
127   -
128   - //bind the 2D texture representing the color map
129   - glBindTexture(GL_TEXTURE_2D, color_tex);
130   -
131   - //Draw the Sphere
132   - int i, j;
133   -
134   - for(i = 1; i <= N-1; i++) {
135   - double phi0 = PI * ((double) (i - 1) / (N-1));
136   - double phi1 = PI * ((double) i / (N-1));
137   -
138   - glBegin(GL_QUAD_STRIP);
139   - for(j = 0; j <= N; j++) {
140   -
141   - //calculate the indices into the function array
142   - int phi0_i = i-1;
143   - int phi1_i = i;
144   - int theta_i = j;
145   - if(theta_i == N)
146   - theta_i = 0;
147   -
148   - double v0 = func[phi0_i * N + theta_i];
149   - double v1 = func[phi1_i * N + theta_i];
150   -
151   - v0 = fabs(v0);
152   - v1 = fabs(v1);
153   -
154   -
155   - double theta = 2 * PI * (double) (j - 1) / N;
156   - double x0 = v0 * cos(theta) * sin(phi0);
157   - double y0 = v0 * sin(theta) * sin(phi0);
158   - double z0 = v0 * cos(phi0);
159   -
160   - double x1 = v1 * cos(theta) * sin(phi1);
161   - double y1 = v1 * sin(theta) * sin(phi1);
162   - double z1 = v1 * cos(phi1);
163   -
164   - glTexCoord2f(theta / (2 * PI), phi0 / PI);
165   - glVertex3f(x0, y0, z0);
166   -
167   - glTexCoord2f(theta / (2 * PI), phi1 / PI);
168   - glVertex3f(x1, y1, z1);
169   - }
170   - glEnd();
171   - }
172   - }
173   -
174   - //draw a wire frame sphere representing the function surface
175   - void gl_draw_wireframe() {
176   -
177   - //PI is used to convert from spherical to cartesian coordinates
178   - //const double PI = 3.14159;
179   -
180   - //bind the 2D texture representing the color map
181   - glDisable(GL_TEXTURE_2D);
182   - glColor3f(0.0f, 0.0f, 0.0f);
183   -
184   - //Draw the Sphere
185   - int i, j;
186   -
187   - for(i = 1; i <= N-1; i++) {
188   - double phi0 = PI * ((double) (i - 1) / (N-1));
189   - double phi1 = PI * ((double) i / (N-1));
190   -
191   - glBegin(GL_LINE_STRIP);
192   - for(j = 0; j <= N; j++) {
193   -
194   - //calculate the indices into the function array
195   - int phi0_i = i-1;
196   - int phi1_i = i;
197   - int theta_i = j;
198   - if(theta_i == N)
199   - theta_i = 0;
200   -
201   - double v0 = func[phi0_i * N + theta_i];
202   - double v1 = func[phi1_i * N + theta_i];
203   -
204   - v0 = fabs(v0);
205   - v1 = fabs(v1);
206   -
207   -
208   - double theta = 2 * PI * (double) (j - 1) / N;
209   - double x0 = WIRE_SCALE * v0 * cos(theta) * sin(phi0);
210   - double y0 = WIRE_SCALE * v0 * sin(theta) * sin(phi0);
211   - double z0 = WIRE_SCALE * v0 * cos(phi0);
212   -
213   - double x1 = WIRE_SCALE * v1 * cos(theta) * sin(phi1);
214   - double y1 = WIRE_SCALE * v1 * sin(theta) * sin(phi1);
215   - double z1 = WIRE_SCALE * v1 * cos(phi1);
216   -
217   - glTexCoord2f(theta / (2 * PI), phi0 / PI);
218   - glVertex3f(x0, y0, z0);
219   -
220   - glTexCoord2f(theta / (2 * PI), phi1 / PI);
221   - glVertex3f(x1, y1, z1);
222   - }
223   - glEnd();
224   - }
225   - }
226   -
227   - void init(unsigned int n){
228   -
229   - //set the sphere resolution
230   - N = n;
231   -
232   - //allocate space for the color map
233   - unsigned int bytes = N * N * sizeof(unsigned char) * 3;
234   - unsigned char* color_image;
235   - color_image = (unsigned char*) malloc(bytes);
236   -
237   - //allocate space for the function
238   - func = (double*) malloc(N * N * sizeof(double));
239   -
240   - //generate a function (temporary)
241   - gen_function();
242   -
243   - //generate a colormap from the function
244   - stim::cpu2cpu<double>(func, color_image, N*N, stim::cmBrewer);
245   -
246   - //prep everything for drawing
247   - gl_prep_draw();
248   -
249   - //generate an OpenGL texture map in the current context
250   - glGenTextures(1, &color_tex);
251   - //bind the texture
252   - glBindTexture(GL_TEXTURE_2D, color_tex);
253   -
254   - //copy the color data from the buffer to the GPU
255   - glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, N, N, 0, GL_RGB, GL_UNSIGNED_BYTE, color_image);
256   -
257   - //initialize all of the texture parameters
258   - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
259   - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
260   - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
261   - glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
262   - glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
263   -
264   - //free the buffer
265   - free(color_image);
266   - }
267   -
268   -
269   - public:
270   -
271   - void glRender(){
272   - //set all OpenGL parameters required for drawing
273   - gl_prep_draw();
274   -
275   - //draw the sphere
276   - gl_draw_sphere();
277   - //gl_draw_wireframe();
278   -
279   - }
280   -
281   - void glInit(unsigned int n){
282   - init(n);
283   - }
284   -
285   - void push(double c){
286   - C.push_back(c);
287   - }
288   -
289   -
290   -
291   -
292   -
293   - }; //end class sph_harmonics
294   -
295   -
296   -
297   -
298   -}
299   -
300   -
301   -#endif
  1 +#ifndef STIM_GL_SPHARMONICS_H
  2 +#define STIM_GL_SPHARMONICS_H
  3 +
  4 +#include <GL/gl.h>
  5 +
  6 +#include <stim/gl/error.h>
  7 +#include <stim/visualization/colormap.h>
  8 +#include <stim/math/spharmonics.h>
  9 +
  10 +namespace stim{
  11 +
  12 +
  13 +template <typename T>
  14 +class gl_spharmonics : public spharmonics<T>{
  15 +
  16 +protected:
  17 +
  18 + T* func; //stores the raw function data (samples at each point)
  19 +
  20 + GLuint color_tex; //texture map that acts as a colormap for the spherical function
  21 +
  22 + unsigned int N; //resolution of the spherical grid
  23 +
  24 + void gen_function(){
  25 +
  26 + //initialize the function to zero
  27 + memset(func, 0, sizeof(double) * N * N);
  28 +
  29 + double theta, phi;
  30 + double result;
  31 + int l, m;
  32 +
  33 +
  34 + l = m = 0;
  35 + //for each coefficient
  36 + for(unsigned int c = 0; c < C.size(); c++){
  37 +
  38 + //iterate through the entire 2D grid representing the function
  39 + for(unsigned int xi = 0; xi < N; xi++){
  40 + for(unsigned int yi = 0; yi < N; yi++){
  41 +
  42 + //get the spherical coordinates for each grid point
  43 + theta = (2 * PI) * ((double)xi / (N-1));
  44 + phi = PI * ((double)yi / (N-1));
  45 +
  46 + //sum the contribution of the current spherical harmonic based on the coefficient
  47 + result = C[c] * SH(l, m, theta, phi);
  48 +
  49 + //store the result in a 2D array (which will later be used as a texture map)
  50 + func[yi * N + xi] += result;
  51 + }
  52 + }
  53 +
  54 + //keep track of m and l here
  55 + m++; //increment m
  56 + //if we're in a new tier, increment l and set m = -l
  57 + if(m > l){
  58 + l++;
  59 + m = -l;
  60 + }
  61 + }
  62 + }
  63 +
  64 + void gl_prep_draw(){
  65 +
  66 + //enable depth testing
  67 + //this has to be used instead of culling because the sphere can have negative values
  68 + glEnable(GL_DEPTH_TEST);
  69 + glDepthMask(GL_TRUE);
  70 + glEnable(GL_TEXTURE_2D); //enable 2D texture mapping
  71 + }
  72 +
  73 + //draw a texture mapped sphere representing the function surface
  74 + void gl_draw_sphere() {
  75 +
  76 + //bind the 2D texture representing the color map
  77 + glBindTexture(GL_TEXTURE_2D, color_tex);
  78 +
  79 + //Draw the Sphere
  80 + int i, j;
  81 +
  82 + for(i = 1; i <= N-1; i++) {
  83 + double phi0 = PI * ((double) (i - 1) / (N-1));
  84 + double phi1 = PI * ((double) i / (N-1));
  85 +
  86 + glBegin(GL_QUAD_STRIP);
  87 + for(j = 0; j <= N; j++) {
  88 +
  89 + //calculate the indices into the function array
  90 + int phi0_i = i-1;
  91 + int phi1_i = i;
  92 + int theta_i = j;
  93 + if(theta_i == N)
  94 + theta_i = 0;
  95 +
  96 + double v0 = func[phi0_i * N + theta_i];
  97 + double v1 = func[phi1_i * N + theta_i];
  98 +
  99 + v0 = fabs(v0);
  100 + v1 = fabs(v1);
  101 +
  102 +
  103 + double theta = 2 * PI * (double) (j - 1) / N;
  104 + double x0 = v0 * cos(theta) * sin(phi0);
  105 + double y0 = v0 * sin(theta) * sin(phi0);
  106 + double z0 = v0 * cos(phi0);
  107 +
  108 + double x1 = v1 * cos(theta) * sin(phi1);
  109 + double y1 = v1 * sin(theta) * sin(phi1);
  110 + double z1 = v1 * cos(phi1);
  111 +
  112 + glTexCoord2f(theta / (2 * PI), phi0 / PI);
  113 + glVertex3f(x0, y0, z0);
  114 +
  115 + glTexCoord2f(theta / (2 * PI), phi1 / PI);
  116 + glVertex3f(x1, y1, z1);
  117 + }
  118 + glEnd();
  119 + }
  120 + }
  121 +
  122 + void gl_init(unsigned int n){
  123 +
  124 + //set the sphere resolution
  125 + N = n;
  126 +
  127 + //allocate space for the color map
  128 + unsigned int bytes = N * N * sizeof(unsigned char) * 3;
  129 + unsigned char* color_image;
  130 + color_image = (unsigned char*) malloc(bytes);
  131 +
  132 + //allocate space for the function
  133 + func = (double*) malloc(N * N * sizeof(double));
  134 +
  135 + //generate a functional representation that will be used for the texture map and vertices
  136 + gen_function();
  137 +
  138 + //generate a colormap from the function
  139 + stim::cpu2cpu<double>(func, color_image, N*N, stim::cmBrewer);
  140 +
  141 + //prep everything for drawing
  142 + gl_prep_draw();
  143 +
  144 + //generate an OpenGL texture map in the current context
  145 + glGenTextures(1, &color_tex);
  146 + //bind the texture
  147 + glBindTexture(GL_TEXTURE_2D, color_tex);
  148 +
  149 + //copy the color data from the buffer to the GPU
  150 + glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, N, N, 0, GL_RGB, GL_UNSIGNED_BYTE, color_image);
  151 +
  152 + //initialize all of the texture parameters
  153 + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_REPEAT);
  154 + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
  155 + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
  156 + glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
  157 + glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE);
  158 +
  159 + //free the buffer
  160 + free(color_image);
  161 + }
  162 +
  163 +public:
  164 +
  165 + void glRender(){
  166 + //set all OpenGL parameters required for drawing
  167 + gl_prep_draw();
  168 +
  169 + //draw the sphere
  170 + gl_draw_sphere();
  171 + }
  172 +
  173 + void glInit(unsigned int n = 256){
  174 + gl_init(n);
  175 + gen_function();
  176 + }
  177 +
  178 +
  179 +}; //end gl_spharmonics
  180 +
  181 +
  182 +}; //end namespace stim
  183 +
  184 +
  185 +
  186 +
  187 +#endif
302 188 \ No newline at end of file
... ...
stim/visualization/obj.h
  Show/Hide comments   View file @487a9b4
... ... @@ -6,6 +6,7 @@
6 6 #include <fstream>
7 7 #include <stdlib.h>
8 8 #include <stim/parser/parser.h>
  9 +#include <stim/math/mathvec.h>
9 10  
10 11 namespace stim{
11 12  
... ... @@ -34,37 +35,14 @@ protected:
34 35  
35 36 enum token_type { OBJ_INVALID, OBJ_V, OBJ_VT, OBJ_VN, OBJ_P, OBJ_L, OBJ_F };
36 37  
37   - struct vertex : public std::vector<T>{
  38 + struct vertex : public stim::vec<T>{
38 39  
39   - using std::vector<T>::push_back;
40   - using std::vector<T>::size;
41   - using std::vector<T>::at;
  40 + using vec<T>::push_back;
  41 + using vec<T>::size;
  42 + using vec<T>::at;
42 43  
43 44 vertex(){}
44 45  
45   - //constructors for quickly building vertices with arbitrary numbers of components
46   - vertex(T x){
47   - push_back(x);
48   - }
49   -
50   - vertex(T x, T y){
51   - push_back(x);
52   - push_back(y);
53   - }
54   -
55   - vertex(T x, T y, T z){
56   - push_back(x);
57   - push_back(y);
58   - push_back(z);
59   - }
60   -
61   - vertex(T x, T y, T z, T w){
62   - push_back(x);
63   - push_back(y);
64   - push_back(z);
65   - push_back(w);
66   - }
67   -
68 46 //constructor creates a vertex from a line string
69 47 vertex(std::string line){
70 48  
... ... @@ -305,6 +283,11 @@ protected:
305 283 }
306 284  
307 285 public:
  286 + /// Constructor loads a Wavefront OBJ file
  287 + obj(std::string filename){
  288 + load(filename);
  289 + }
  290 +
308 291 //functions for setting the texture coordinate for the next vertex
309 292 void TexCoord(T x){ update_vt(vertex(x));}
310 293 void TexCoord(T x, T y){ update_vt(vertex(x, y));}
... ... @@ -504,6 +487,47 @@ public:
504 487  
505 488 }
506 489  
  490 + /// Return the number of position vertices in the OBJ model
  491 + unsigned int numV(){
  492 + return V.size();
  493 + }
  494 +
  495 + /// Retrieve the vertex stored in index i
  496 +
  497 + /// @param i is the desired vertex index
  498 + stim::vec<T> getV(unsigned int i){
  499 +
  500 + stim::vec<T> v = V[i];
  501 + return v;
  502 + }
  503 +
  504 + stim::vec<T> centroid(){
  505 +
  506 + //get the number of coordinates
  507 + unsigned int N = V[0].size();
  508 +
  509 + //allocate space for the minimum and maximum coordinate points (bounding box corners)
  510 + stim::vec<float> vmin, vmax;
  511 + vmin.resize(N);
  512 + vmax.resize(N);
  513 +
  514 + //find the minimum and maximum value for each coordinate
  515 + unsigned int NV = V.size();
  516 + for(unsigned int v = 0; v < NV; v++)
  517 + for(unsigned int i = 0; i < N; i++){
  518 +
  519 + if(V[v][i] < vmin[i])
  520 + vmin[i] = V[v][i];
  521 + if(V[v][i] > vmax[i])
  522 + vmax[i] = V[v][i];
  523 + }
  524 +
  525 + //find the centroid using the min and max points
  526 + stim::vec<T> c = vmin * 0.5 + vmax * 0.5;
  527 +
  528 + return c;
  529 + }
  530 +
507 531  
508 532 };
509 533  
... ...