codebase / stimlib

Commit 7731cf246c65fa7c936fe89618a0701405f7382c

Authored by David Mayerich
1 parent 57729e5b

fixed precision problems

Showing 12 changed files with 46 additions and 1273 deletions   Show diff stats
Inline Side-by-side
math/complex.h
  Show/Hide comments   View file @7731cf2
... ... @@ -5,7 +5,7 @@ and can therefore be used in CUDA code and on CUDA devices.
5 5 #ifndef RTS_COMPLEX
6 6 #define RTS_COMPLEX
7 7  
8   -#include "cuda_callable.h"
  8 +#include "../cuda/callable.h"
9 9 #include <cmath>
10 10 #include <string>
11 11 #include <sstream>
... ... @@ -22,8 +22,8 @@ struct complex
22 22 //default constructor
23 23 CUDA_CALLABLE complex()
24 24 {
25   - r = 0.0;
26   - i = 0.0;
  25 + r = 0;
  26 + i = 0;
27 27 }
28 28  
29 29 //access methods
... ... @@ -235,7 +235,7 @@ struct complex
235 235  
236 236 //find the square root
237 237 T a_p = std::sqrt(a);
238   - T theta_p = theta/2.0;
  238 + T theta_p = theta/2.0f;
239 239  
240 240 //convert back to cartesian coordinates
241 241 result.r = a_p * std::cos(theta_p);
... ... @@ -262,7 +262,7 @@ struct complex
262 262  
263 263 CUDA_CALLABLE bool operator==(T rhs)
264 264 {
265   - if(r == rhs && i == (T)0.0)
  265 + if(r == rhs && i == 0)
266 266 return true;
267 267 return false;
268 268 }
... ...
math/complex.h~ deleted
View file @57729e5
1   -/*RTS Complex number class. This class is CUDA compatible,
2   -and can therefore be used in CUDA code and on CUDA devices.
3   -*/
4   -
5   -#ifndef RTS_COMPLEX
6   -#define RTS_COMPLEX
7   -
8   -#include "cuda_callable.h"
9   -#include <cmath>
10   -#include <string>
11   -#include <sstream>
12   -#include <iostream>
13   -
14   -namespace rts
15   -{
16   -
17   -template <class T>
18   -struct complex
19   -{
20   - T r, i;
21   -
22   - //default constructor
23   - CUDA_CALLABLE complex()
24   - {
25   - r = 0.0;
26   - i = 0.0;
27   - }
28   -
29   - //access methods
30   - CUDA_CALLABLE T real()
31   - {
32   - return r;
33   - }
34   -
35   - CUDA_CALLABLE T real(T r_val)
36   - {
37   - r = r_val;
38   - return r_val;
39   - }
40   -
41   - CUDA_CALLABLE T imag()
42   - {
43   - return i;
44   - }
45   - CUDA_CALLABLE T imag(T i_val)
46   - {
47   - i = i_val;
48   - return i_val;
49   - }
50   -
51   - //constructor when given real and imaginary values
52   - CUDA_CALLABLE complex(T r, T i)
53   - {
54   - this->r = r;
55   - this->i = i;
56   - }
57   -
58   - //return the current value multiplied by i
59   - CUDA_CALLABLE complex<T> imul()
60   - {
61   - complex<T> result;
62   - result.r = -i;
63   - result.i = r;
64   -
65   - return result;
66   - }
67   -
68   - //ARITHMETIC OPERATORS--------------------
69   -
70   - //binary + operator (returns the result of adding two complex values)
71   - CUDA_CALLABLE complex<T> operator+ (const complex<T> rhs)
72   - {
73   - complex<T> result;
74   - result.r = r + rhs.r;
75   - result.i = i + rhs.i;
76   - return result;
77   - }
78   -
79   - CUDA_CALLABLE complex<T> operator+ (const T rhs)
80   - {
81   - complex<T> result;
82   - result.r = r + rhs;
83   - result.i = i;
84   - return result;
85   - }
86   -
87   - //binary - operator (returns the result of adding two complex values)
88   - CUDA_CALLABLE complex<T> operator- (const complex<T> rhs)
89   - {
90   - complex<T> result;
91   - result.r = r - rhs.r;
92   - result.i = i - rhs.i;
93   - return result;
94   - }
95   -
96   - //binary - operator (returns the result of adding two complex values)
97   - CUDA_CALLABLE complex<T> operator- (const T rhs)
98   - {
99   - complex<T> result;
100   - result.r = r - rhs;
101   - result.i = i;
102   - return result;
103   - }
104   -
105   - //binary MULTIPLICATION operators (returns the result of multiplying complex values)
106   - CUDA_CALLABLE complex<T> operator* (const complex<T> rhs)
107   - {
108   - complex<T> result;
109   - result.r = r * rhs.r - i * rhs.i;
110   - result.i = r * rhs.i + i * rhs.r;
111   - return result;
112   - }
113   - CUDA_CALLABLE complex<T> operator* (const T rhs)
114   - {
115   - return complex<T>(r * rhs, i * rhs);
116   - }
117   -
118   - //binary DIVISION operators (returns the result of dividing complex values)
119   - CUDA_CALLABLE complex<T> operator/ (const complex<T> rhs)
120   - {
121   - complex<T> result;
122   - T denom = rhs.r * rhs.r + rhs.i * rhs.i;
123   - result.r = (r * rhs.r + i * rhs.i) / denom;
124   - result.i = (- r * rhs.i + i * rhs.r) / denom;
125   -
126   - return result;
127   - }
128   - CUDA_CALLABLE complex<T> operator/ (const T rhs)
129   - {
130   - return complex<T>(r / rhs, i / rhs);
131   - }
132   -
133   - //ASSIGNMENT operators-----------------------------------
134   - CUDA_CALLABLE complex<T> & operator=(const complex<T> &rhs)
135   - {
136   - //check for self-assignment
137   - if(this != &rhs)
138   - {
139   - this->r = rhs.r;
140   - this->i = rhs.i;
141   - }
142   - return *this;
143   - }
144   - CUDA_CALLABLE complex<T> & operator=(const T &rhs)
145   - {
146   - this->r = rhs;
147   - this->i = 0;
148   -
149   - return *this;
150   - }
151   -
152   - //arithmetic assignment operators
153   - CUDA_CALLABLE complex<T> operator+=(const complex<T> &rhs)
154   - {
155   - *this = *this + rhs;
156   - return *this;
157   - }
158   - CUDA_CALLABLE complex<T> operator+=(const T &rhs)
159   - {
160   - *this = *this + rhs;
161   - return *this;
162   - }
163   -
164   - CUDA_CALLABLE complex<T> operator*=(const complex<T> &rhs)
165   - {
166   - *this = *this * rhs;
167   - return *this;
168   - }
169   - CUDA_CALLABLE complex<T> operator*=(const T &rhs)
170   - {
171   - *this = *this * rhs;
172   - return *this;
173   - }
174   - //divide and assign
175   - CUDA_CALLABLE complex<T> operator/=(const complex<T> &rhs)
176   - {
177   - *this = *this / rhs;
178   - return *this;
179   - }
180   - CUDA_CALLABLE complex<T> operator/=(const T &rhs)
181   - {
182   - *this = *this / rhs;
183   - return *this;
184   - }
185   -
186   - //absolute value operator (returns the absolute value of the complex number)
187   - CUDA_CALLABLE T abs()
188   - {
189   - return std::sqrt(r * r + i * i);
190   - }
191   -
192   - CUDA_CALLABLE complex<T> log()
193   - {
194   - complex<T> result;
195   - result.r = std::log(std::sqrt(r * r + i * i));
196   - result.i = std::atan2(i, r);
197   -
198   -
199   - return result;
200   - }
201   -
202   - CUDA_CALLABLE complex<T> exp()
203   - {
204   - complex<T> result;
205   -
206   - T e_r = std::exp(r);
207   - result.r = e_r * std::cos(i);
208   - result.i = e_r * std::sin(i);
209   -
210   - return result;
211   - }
212   -
213   - /*CUDA_CALLABLE complex<T> pow(int y)
214   - {
215   -
216   - return pow((double)y);
217   - }*/
218   -
219   - CUDA_CALLABLE complex<T> pow(T y)
220   - {
221   - complex<T> result;
222   -
223   - result = log() * y;
224   -
225   - return result.exp();
226   - }
227   -
228   - CUDA_CALLABLE complex<T> sqrt()
229   - {
230   - complex<T> result;
231   -
232   - //convert to polar coordinates
233   - T a = std::sqrt(r*r + i*i);
234   - T theta = std::atan2(i, r);
235   -
236   - //find the square root
237   - T a_p = std::sqrt(a);
238   - T theta_p = theta/2.0;
239   -
240   - //convert back to cartesian coordinates
241   - result.r = a_p * std::cos(theta_p);
242   - result.i = a_p * std::sin(theta_p);
243   -
244   - return result;
245   - }
246   -
247   - std::string toStr()
248   - {
249   - std::stringstream ss;
250   - ss<<"("<<r<<","<<i<<")";
251   -
252   - return ss.str();
253   - }
254   -
255   - //COMPARISON operators
256   - CUDA_CALLABLE bool operator==(complex<T> rhs)
257   - {
258   - if(r == rhs.r && i == rhs.i)
259   - return true;
260   - return false;
261   - }
262   -
263   - CUDA_CALLABLE bool operator==(T rhs)
264   - {
265   - if(r == rhs && i == (T)0.0)
266   - return true;
267   - return false;
268   - }
269   -
270   -};
271   -
272   -} //end RTS namespace
273   -
274   -//addition
275   -template<typename T>
276   -CUDA_CALLABLE static rts::complex<T> operator+(const double a, const rts::complex<T> b)
277   -{
278   - return rts::complex<T>(a + b.r, b.i);
279   -}
280   -
281   -//subtraction with a real value
282   -template<typename T>
283   -CUDA_CALLABLE static rts::complex<T> operator-(const double a, const rts::complex<T> b)
284   -{
285   - return rts::complex<T>(a - b.r, -b.i);
286   -}
287   -
288   -//minus sign
289   -template<typename T>
290   -CUDA_CALLABLE static rts::complex<T> operator-(const rts::complex<T> &rhs)
291   -{
292   - return rts::complex<T>(-rhs.r, -rhs.i);
293   -}
294   -
295   -//multiply a T value by a complex value
296   -template<typename T>
297   -CUDA_CALLABLE static rts::complex<T> operator*(const double a, const rts::complex<T> b)
298   -{
299   - return rts::complex<T>((T)a * b.r, (T)a * b.i);
300   -}
301   -
302   -//divide a T value by a complex value
303   -template<typename T>
304   -CUDA_CALLABLE static rts::complex<T> operator/(const double a, const rts::complex<T> b)
305   -{
306   - //return complex<T>(a * b.r, a * b.i);
307   - rts::complex<T> result;
308   -
309   - T denom = b.r * b.r + b.i * b.i;
310   -
311   - result.r = (a * b.r) / denom;
312   - result.i = -(a * b.i) / denom;
313   -
314   - return result;
315   -}
316   -
317   -//POW function
318   -/*template<typename T>
319   -CUDA_CALLABLE static complex<T> pow(complex<T> x, int y)
320   -{
321   - return x.pow(y);
322   -}*/
323   -
324   -template<typename T>
325   -CUDA_CALLABLE static rts::complex<T> pow(rts::complex<T> x, T y)
326   -{
327   - return x.pow(y);
328   -}
329   -
330   -//log function
331   -template<typename T>
332   -CUDA_CALLABLE static rts::complex<T> log(rts::complex<T> x)
333   -{
334   - return x.log();
335   -}
336   -
337   -//exp function
338   -template<typename T>
339   -CUDA_CALLABLE static rts::complex<T> exp(rts::complex<T> x)
340   -{
341   - return x.exp();
342   -}
343   -
344   -//sqrt function
345   -template<typename T>
346   -CUDA_CALLABLE static rts::complex<T> sqrt(rts::complex<T> x)
347   -{
348   - return x.sqrt();
349   -}
350   -
351   -
352   -template <typename T>
353   -CUDA_CALLABLE static T abs(rts::complex<T> a)
354   -{
355   - return a.abs();
356   -}
357   -
358   -template <typename T>
359   -CUDA_CALLABLE static T real(rts::complex<T> a)
360   -{
361   - return a.r;
362   -}
363   -
364   -//template <typename T>
365   -CUDA_CALLABLE static float real(float a)
366   -{
367   - return a;
368   -}
369   -
370   -template <typename T>
371   -CUDA_CALLABLE static T imag(rts::complex<T> a)
372   -{
373   - return a.i;
374   -}
375   -
376   -//trigonometric functions
377   -template<class A>
378   -CUDA_CALLABLE rts::complex<A> sin(const rts::complex<A> x)
379   -{
380   - rts::complex<A> result;
381   - result.r = std::sin(x.r) * std::cosh(x.i);
382   - result.i = std::cos(x.r) * std::sinh(x.i);
383   -
384   - return result;
385   -}
386   -
387   -template<class A>
388   -CUDA_CALLABLE rts::complex<A> cos(const rts::complex<A> x)
389   -{
390   - rts::complex<A> result;
391   - result.r = std::cos(x.r) * std::cosh(x.i);
392   - result.i = -(std::sin(x.r) * std::sinh(x.i));
393   -
394   - return result;
395   -}
396   -
397   -
398   -template<class A>
399   -std::ostream& operator<<(std::ostream& os, rts::complex<A> x)
400   -{
401   - os<<x.toStr();
402   - return os;
403   -}
404   -
405   -#if __GNUC__ > 3 && __GNUC_MINOR__ > 7
406   -template<class T> using rtsComplex = rts::complex<T>;
407   -#endif
408   -
409   -
410   -
411   -#endif
math/matrix.h~ deleted
View file @57729e5
1   -#ifndef RTS_MATRIX_H
2   -#define RTS_MATRIX_H
3   -
4   -//#include "rts/vector.h"
5   -#include <string.h>
6   -#include <iostream>
7   -
8   -namespace rts
9   -{
10   -
11   -template <class T, int N>
12   -struct matrix
13   -{
14   - //the matrix will be stored in column-major order (compatible with OpenGL)
15   - T M[N*N];
16   -
17   - matrix()
18   - {
19   - for(int r=0; r<N; r++)
20   - for(int c=0; c<N; c++)
21   - if(r == c)
22   - (*this)(r, c) = 1;
23   - else
24   - (*this)(r, c) = 0;
25   - }
26   -
27   - T& operator()(int row, int col)
28   - {
29   - return M[col * N + row];
30   - }
31   -
32   - matrix<T, N> operator=(T rhs)
33   - {
34   - int Nsq = N*N;
35   - for(int i=0; i<Nsq; i++)
36   - M[i] = rhs;
37   -
38   - return *this;
39   - }
40   -
41   - /*matrix<T, N> operator=(matrix<T, N> rhs)
42   - {
43   - for(int i=0; i<N; i++)
44   - M[i] = rhs.M[i];
45   -
46   - return *this;
47   - }*/
48   -
49   - vector<T, N> operator*(vector<T, N> rhs)
50   - {
51   - vector<T, N> result;
52   -
53   - for(int r=0; r<N; r++)
54   - for(int c=0; c<N; c++)
55   - result[r] += (*this)(r, c) * rhs[c];
56   -
57   - return result;
58   - }
59   -
60   - std::string toStr()
61   - {
62   - std::stringstream ss;
63   -
64   - for(int r = 0; r < N; r++)
65   - {
66   - ss<<"| ";
67   - for(int c=0; c<N; c++)
68   - {
69   - ss<<(*this)(r, c)<<" ";
70   - }
71   - ss<<"|"<<std::endl;
72   - }
73   -
74   - return ss.str();
75   - }
76   -
77   -
78   -
79   -
80   -};
81   -
82   -} //end namespace rts
83   -
84   -template <typename T, int N>
85   -std::ostream& operator<<(std::ostream& os, rts::matrix<T, N> M)
86   -{
87   - os<<M.toStr();
88   - return os;
89   -}
90   -
91   -#if __GNUC__ > 3 && __GNUC_MINOR__ > 7
92   -template<class T, int N> using rtsMatrix = rts::matrix<T, N>;
93   -#endif
94   -
95   -#endif
math/point.h
  Show/Hide comments   View file @7731cf2
... ... @@ -19,16 +19,27 @@ struct point
19 19 }
20 20  
21 21 //efficiency constructor, makes construction easier for 1D-4D vectors
22   - CUDA_CALLABLE point(T x, T y = (T)0.0, T z = (T)0.0, T w = (T)0.0)
  22 + CUDA_CALLABLE point(T x)
23 23 {
24   - if(N >= 1)
25   - p[0] = x;
26   - if(N >= 2)
27   - p[1] = y;
28   - if(N >= 3)
29   - p[2] = z;
30   - if(N >= 4)
31   - p[3] = w;
  24 + p[0] = x;
  25 + }
  26 + CUDA_CALLABLE point(T x, T y)
  27 + {
  28 + p[0] = x;
  29 + p[1] = y;
  30 + }
  31 + CUDA_CALLABLE point(T x, T y, T z)
  32 + {
  33 + p[0] = x;
  34 + p[1] = y;
  35 + p[2] = z;
  36 + }
  37 + CUDA_CALLABLE point(T x, T y, T z, T w)
  38 + {
  39 + p[0] = x;
  40 + p[1] = y;
  41 + p[2] = z;
  42 + p[3] = w;
32 43 }
33 44  
34 45 //arithmetic operators
... ...
math/point.h~ deleted
View file @57729e5
1   -#ifndef RTS_rtsPoint_H
2   -#define RTS_rtsPoint_H
3   -
4   -#include "rts/math/vector.h"
5   -#include <string.h>
6   -#include "rts/cuda/callable.h"
7   -
8   -namespace rts
9   -{
10   -
11   -template <class T, int N>
12   -struct point
13   -{
14   - T p[N];
15   -
16   - CUDA_CALLABLE point()
17   - {
18   -
19   - }
20   -
21   - //efficiency constructor, makes construction easier for 1D-4D vectors
22   - CUDA_CALLABLE point(T x, T y = (T)0.0, T z = (T)0.0, T w = (T)0.0)
23   - {
24   - if(N >= 1)
25   - p[0] = x;
26   - if(N >= 2)
27   - p[1] = y;
28   - if(N >= 3)
29   - p[2] = z;
30   - if(N >= 4)
31   - p[3] = w;
32   - }
33   -
34   - //arithmetic operators
35   - CUDA_CALLABLE rts::point<T, N> operator+(vector<T, N> v)
36   - {
37   - rts::point<T, N> r;
38   -
39   - //calculate the position of the resulting point
40   - for(int i=0; i<N; i++)
41   - r.p[i] = p[i] + v.v[i];
42   -
43   - return r;
44   - }
45   - CUDA_CALLABLE rts::point<T, N> operator-(vector<T, N> v)
46   - {
47   - rts::point<T, N> r;
48   -
49   - //calculate the position of the resulting point
50   - for(int i=0; i<N; i++)
51   - r.p[i] = p[i] - v.v[i];
52   -
53   - return r;
54   - }
55   - CUDA_CALLABLE vector<T, N> operator-(point<T, N> rhs)
56   - {
57   - vector<T, N> r;
58   -
59   - //calculate the position of the resulting point
60   - for(int i=0; i<N; i++)
61   - r.v[i] = p[i] - rhs.p[i];
62   -
63   - return r;
64   - }
65   - CUDA_CALLABLE rts::point<T, N> operator*(T rhs)
66   - {
67   - rts::point<T, N> r;
68   -
69   - //calculate the position of the resulting point
70   - for(int i=0; i<N; i++)
71   - r.p[i] = p[i] * rhs;
72   -
73   - return r;
74   - }
75   -
76   - CUDA_CALLABLE point(const T(&data)[N])
77   - {
78   - memcpy(p, data, sizeof(T) * N);
79   - }
80   -
81   - std::string toStr()
82   - {
83   - std::stringstream ss;
84   -
85   - ss<<"(";
86   - for(int i=0; i<N; i++)
87   - {
88   - ss<<p[i];
89   - if(i != N-1)
90   - ss<<", ";
91   - }
92   - ss<<")";
93   -
94   - return ss.str();
95   - }
96   -
97   - //bracket operator
98   - CUDA_CALLABLE T& operator[](int i)
99   - {
100   - return p[i];
101   - }
102   -
103   -};
104   -
105   -} //end namespace rts
106   -
107   -template <typename T, int N>
108   -std::ostream& operator<<(std::ostream& os, rts::point<T, N> p)
109   -{
110   - os<<p.toStr();
111   - return os;
112   -}
113   -
114   -//arithmetic
115   -template <typename T, int N>
116   -CUDA_CALLABLE rts::point<T, N> operator*(T lhs, rts::point<T, N> rhs)
117   -{
118   - rts::point<T, N> r;
119   -
120   - return rhs * lhs;
121   -}
122   -
123   -#if __GNUC__ > 3 && __GNUC_MINOR__ > 7
124   -template<class T, int N> using rtsPoint = rts::point<T, N>;
125   -#endif
126   -
127   -#endif
math/quad.h
  Show/Hide comments   View file @7731cf2
... ... @@ -69,11 +69,11 @@ struct quad
69 69  
70 70 //calculate point B
71 71 rts::point<T, 3> B;
72   - B = A + v0 * 0.5 + v1 * 0.5;
  72 + B = A + v0 * 0.5f + v1 * 0.5f;
73 73  
74 74 //calculate rtsPoint C
75 75 rts::point<T, 3> C;
76   - C = A + v0 * 0.5 - v1 * 0.5;
  76 + C = A + v0 * 0.5f - v1 * 0.5f;
77 77  
78 78 //calculate X and Y
79 79 X = B - A;
... ... @@ -110,7 +110,7 @@ struct quad
110 110 Y = Y * height;
111 111  
112 112 //set the corner of the plane
113   - A = c - X * 0.5 - Y * 0.5;
  113 + A = c - X * 0.5f - Y * 0.5f;
114 114  
115 115 std::cout<<X<<std::endl;
116 116 }
... ... @@ -155,13 +155,13 @@ struct quad
155 155 //scales the plane by a scalar value
156 156  
157 157 //compute the center point
158   - rts::point<T, N> c = A + X*0.5 + Y*0.5;
  158 + rts::point<T, N> c = A + X*0.5f + Y*0.5f;
159 159  
160 160 //create the new quadangle
161 161 quad<T, N> result;
162 162 result.X = X * rhs;
163 163 result.Y = Y * rhs;
164   - result.A = c - result.X*0.5 - result.Y*0.5;
  164 + result.A = c - result.X*0.5f - result.Y*0.5f;
165 165  
166 166 return result;
167 167  
... ... @@ -192,7 +192,7 @@ struct quad
192 192 T dc = (A+Y - p).len();
193 193 T dd = (A+X+Y - p).len();
194 194  
195   - return std::max( da, std::max(db, std::max(dc, dd) ) );
  195 + return fmax( da, fmax(db, fmax(dc, dd) ) );
196 196 }
197 197 };
198 198  
... ...
math/quad.h~ deleted
View file @57729e5
1   -#ifndef RTS_RECT_H
2   -#define RTS_RECT_H
3   -
4   -//enable CUDA_CALLABLE macro
5   -#include "rts/cuda/callable.h"
6   -#include "rts/math/vector.h"
7   -#include "rts/math/point.h"
8   -#include "rts/math/triangle.h"
9   -#include "rts/math/quaternion.h"
10   -#include <iostream>
11   -#include <algorithm>
12   -
13   -namespace rts{
14   -
15   -//template for a quadangle class in ND space
16   -template <class T, int N>
17   -struct quad
18   -{
19   - /*
20   - C------------------>O
21   - ^ ^
22   - | |
23   - Y |
24   - | |
25   - | |
26   - A---------X-------->B
27   - */
28   -
29   - /*T A[N];
30   - T B[N];
31   - T C[N];*/
32   -
33   - rts::point<T, N> A;
34   - rts::vector<T, N> X;
35   - rts::vector<T, N> Y;
36   -
37   -
38   - CUDA_CALLABLE quad()
39   - {
40   -
41   - }
42   -
43   - CUDA_CALLABLE quad(point<T, N> a, point<T, N> b, point<T, N> c)
44   - {
45   -
46   - A = a;
47   - X = b - a;
48   - Y = c - a;
49   -
50   - }
51   -
52   - /****************************************************************
53   - Constructor - create a quad from two points and a normal
54   - ****************************************************************/
55   - CUDA_CALLABLE quad(rts::point<T, N> pMin, rts::point<T, N> pMax, rts::vector<T, N> normal)
56   - {
57   -
58   - //assign the corner point
59   - A = pMin;
60   -
61   - //compute the vector from pMin to pMax
62   - rts::vector<T, 3> v0;
63   - v0 = pMax - pMin;
64   -
65   - //compute the cross product of A and the plane normal
66   - rts::vector<T, 3> v1;
67   - v1 = v0.cross(normal);
68   -
69   -
70   - //calculate point B
71   - rts::point<T, 3> B;
72   - B = A + v0 * 0.5 + v1 * 0.5;
73   -
74   - //calculate rtsPoint C
75   - rts::point<T, 3> C;
76   - C = A + v0 * 0.5 - v1 * 0.5;
77   -
78   - //calculate X and Y
79   - X = B - A;
80   - Y = C - A;
81   - }
82   -
83   - /*******************************************************************
84   - Constructor - create a quad from a position, normal, and rotation
85   - *******************************************************************/
86   - CUDA_CALLABLE quad(rts::point<T, N> c, rts::vector<T, N> normal, T width, T height, T theta)
87   - {
88   -
89   - //compute the X direction - start along world-space X
90   - Y = rts::vector<T, N>(0, 1, 0);
91   - if(Y == normal)
92   - Y = rts::vector<T, N>(0, 0, 1);
93   -
94   - X = Y.cross(normal).norm();
95   -
96   - std::cout<<X<<std::endl;
97   -
98   - //rotate the X axis by theta radians
99   - rts::quaternion<T> q;
100   - q.CreateRotation(theta, normal);
101   - X = q.toMatrix3() * X;
102   - Y = normal.cross(X);
103   -
104   - //normalize everything
105   - X = X.norm();
106   - Y = Y.norm();
107   -
108   - //scale to match the quad width and height
109   - X = X * width;
110   - Y = Y * height;
111   -
112   - //set the corner of the plane
113   - A = c - X * 0.5 - Y * 0.5;
114   -
115   - std::cout<<X<<std::endl;
116   - }
117   -
118   - /*******************************************
119   - Return the normal for the quad
120   - *******************************************/
121   - CUDA_CALLABLE rts::vector<T, N> n()
122   - {
123   - return (X.cross(Y)).norm();
124   - }
125   -
126   - CUDA_CALLABLE rts::point<T, N> p(T a, T b)
127   - {
128   - rts::point<T, N> result;
129   - //given the two parameters a, b = [0 1], returns the position in world space
130   - result = A + X * a + Y * b;
131   -
132   - return result;
133   - }
134   -
135   - CUDA_CALLABLE rts::point<T, N> operator()(T a, T b)
136   - {
137   - return p(a, b);
138   - }
139   -
140   - std::string toStr()
141   - {
142   - std::stringstream ss;
143   -
144   - ss<<"A = "<<A<<std::endl;
145   - ss<<"B = "<<A + X<<std::endl;
146   - ss<<"C = "<<A + X + Y<<std::endl;
147   - ss<<"D = "<<A + Y<<std::endl;
148   -
149   - return ss.str();
150   -
151   - }
152   -
153   - CUDA_CALLABLE quad<T, N> operator*(T rhs)
154   - {
155   - //scales the plane by a scalar value
156   -
157   - //compute the center point
158   - rts::point<T, N> c = A + X*0.5 + Y*0.5;
159   -
160   - //create the new quadangle
161   - quad<T, N> result;
162   - result.X = X * rhs;
163   - result.Y = Y * rhs;
164   - result.A = c - result.X*0.5 - result.Y*0.5;
165   -
166   - return result;
167   -
168   - }
169   -
170   - CUDA_CALLABLE T dist(point<T, N> p)
171   - {
172   - //compute the distance between a point and this quad
173   -
174   - //first break the quad up into two triangles
175   - triangle<T, N> T0(A, A+X, A+Y);
176   - triangle<T, N> T1(A+X+Y, A+X, A+Y);
177   -
178   -
179   - ptype d0 = T0.dist(p);
180   - ptype d1 = T1.dist(p);
181   -
182   - if(d0 < d1)
183   - return d0;
184   - else
185   - return d1;
186   - }
187   -
188   - CUDA_CALLABLE T dist_max(point<T, N> p)
189   - {
190   - T da = (A - p).len();
191   - T db = (A+X - p).len();
192   - T dc = (A+Y - p).len();
193   - T dd = (A+X+Y - p).len();
194   -
195   - return std::max( da, max(db, max(dc, dd) ) );
196   - }
197   -};
198   -
199   -} //end namespace rts
200   -
201   -template <typename T, int N>
202   -std::ostream& operator<<(std::ostream& os, rts::quad<T, N> R)
203   -{
204   - os<<R.toStr();
205   - return os;
206   -}
207   -
208   -
209   -#endif
math/quaternion.h
  Show/Hide comments   View file @7731cf2
... ... @@ -41,10 +41,10 @@ template&lt;typename T&gt;
41 41 void quaternion<T>::CreateRotation(T theta, T axis_x, T axis_y, T axis_z)
42 42 {
43 43 //assign the given Euler rotation to this quaternion
44   - w = (T)cos(theta/2.0);
45   - x = axis_x*(T)sin(theta/2.0);
46   - y = axis_y*(T)sin(theta/2.0);
47   - z = axis_z*(T)sin(theta/2.0);
  44 + w = (T)cos(theta/2);
  45 + x = axis_x*(T)sin(theta/2);
  46 + y = axis_y*(T)sin(theta/2);
  47 + z = axis_z*(T)sin(theta/2);
48 48 }
49 49  
50 50 template<typename T>
... ... @@ -93,20 +93,20 @@ matrix&lt;T, 3&gt; quaternion&lt;T&gt;::toMatrix3()
93 93 yy = y * y2; yz = y * z2; zz = z * z2;
94 94 wx = w * x2; wy = w * y2; wz = w * z2;
95 95  
96   - result(0, 0) = (T)1.0 - (yy + zz);
  96 + result(0, 0) = 1 - (yy + zz);
97 97 result(0, 1) = xy - wz;
98 98  
99 99 result(0, 2) = xz + wy;
100 100  
101 101 result(1, 0) = xy + wz;
102   - result(1, 1) = (T)1.0 - (xx + zz);
  102 + result(1, 1) = 1 - (xx + zz);
103 103  
104 104 result(1, 2) = yz - wx;
105 105  
106 106 result(2, 0) = xz - wy;
107 107 result(2, 1) = yz + wx;
108 108  
109   - result(2, 2) = (T)1.0 - (xx + yy);
  109 + result(2, 2) = 1 - (xx + yy);
110 110  
111 111 return result;
112 112 }
... ... @@ -127,22 +127,22 @@ matrix&lt;T, 4&gt; quaternion&lt;T&gt;::toMatrix4()
127 127 yy = y * y2; yz = y * z2; zz = z * z2;
128 128 wx = w * x2; wy = w * y2; wz = w * z2;
129 129  
130   - result(0, 0) = (T)1.0 - (yy + zz);
  130 + result(0, 0) = 1 - (yy + zz);
131 131 result(0, 1) = xy - wz;
132 132  
133 133 result(0, 2) = xz + wy;
134 134  
135 135 result(1, 0) = xy + wz;
136   - result(1, 1) = (T)1.0 - (xx + zz);
  136 + result(1, 1) = 1 - (xx + zz);
137 137  
138 138 result(1, 2) = yz - wx;
139 139  
140 140 result(2, 0) = xz - wy;
141 141 result(2, 1) = yz + wx;
142 142  
143   - result(2, 2) = (T)1.0 - (xx + yy);
  143 + result(2, 2) = 1 - (xx + yy);
144 144  
145   - result(3, 3) = (T)1.0;
  145 + result(3, 3) = 1;
146 146  
147 147 return result;
148 148 }
... ... @@ -150,7 +150,7 @@ matrix&lt;T, 4&gt; quaternion&lt;T&gt;::toMatrix4()
150 150 template<typename T>
151 151 quaternion<T>::quaternion()
152 152 {
153   - w=0.0; x=0.0; y=0.0; z=0.0;
  153 + w=0; x=0; y=0; z=0;
154 154 }
155 155  
156 156 template<typename T>
... ...
math/quaternion.h~ deleted
View file @57729e5
1   -#include "rts/math/matrix.h"
2   -
3   -#ifndef RTS_QUATERNION_H
4   -#define RTS_QUATERNION_H
5   -
6   -namespace rts{
7   -
8   -template<typename T>
9   -class quaternion
10   -{
11   -public:
12   - T w;
13   - T x;
14   - T y;
15   - T z;
16   -
17   - void normalize();
18   - void CreateRotation(T theta, T axis_x, T axis_y, T axis_z);
19   - void CreateRotation(T theta, vector<T, 3> axis);
20   - quaternion<T> operator*(quaternion<T> &rhs);
21   - matrix<T, 3> toMatrix3();
22   - matrix<T, 4> toMatrix4();
23   -
24   -
25   - quaternion();
26   - quaternion(T w, T x, T y, T z);
27   -
28   -};
29   -
30   -template<typename T>
31   -void quaternion<T>::normalize()
32   -{
33   - double length=sqrt(w*w + x*x + y*y + z*z);
34   - w=w/length;
35   - x=x/length;
36   - y=y/length;
37   - z=z/length;
38   -}
39   -
40   -template<typename T>
41   -void quaternion<T>::CreateRotation(T theta, T axis_x, T axis_y, T axis_z)
42   -{
43   - //assign the given Euler rotation to this quaternion
44   - w = (T)cos(theta/2.0);
45   - x = axis_x*(T)sin(theta/2.0);
46   - y = axis_y*(T)sin(theta/2.0);
47   - z = axis_z*(T)sin(theta/2.0);
48   -}
49   -
50   -template<typename T>
51   -void quaternion<T>::CreateRotation(T theta, vector<T, 3> axis)
52   -{
53   - CreateRotation(theta, axis[0], axis[1], axis[2]);
54   -}
55   -
56   -template<typename T>
57   -quaternion<T> quaternion<T>::operator *(quaternion<T> &param)
58   -{
59   - float A, B, C, D, E, F, G, H;
60   -
61   -
62   - A = (w + x)*(param.w + param.x);
63   - B = (z - y)*(param.y - param.z);
64   - C = (w - x)*(param.y + param.z);
65   - D = (y + z)*(param.w - param.x);
66   - E = (x + z)*(param.x + param.y);
67   - F = (x - z)*(param.x - param.y);
68   - G = (w + y)*(param.w - param.z);
69   - H = (w - y)*(param.w + param.z);
70   -
71   - quaternion<T> result;
72   - result.w = B + (-E - F + G + H) /2;
73   - result.x = A - (E + F + G + H)/2;
74   - result.y = C + (E - F + G - H)/2;
75   - result.z = D + (E - F - G + H)/2;
76   -
77   - return result;
78   -}
79   -
80   -template<typename T>
81   -matrix<T, 3> quaternion<T>::toMatrix3()
82   -{
83   - matrix<T, 3> result;
84   -
85   -
86   - T wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
87   -
88   -
89   - // calculate coefficients
90   - x2 = x + x; y2 = y + y;
91   - z2 = z + z;
92   - xx = x * x2; xy = x * y2; xz = x * z2;
93   - yy = y * y2; yz = y * z2; zz = z * z2;
94   - wx = w * x2; wy = w * y2; wz = w * z2;
95   -
96   - result(0, 0) = (T)1.0 - (yy + zz);
97   - result(0, 1) = xy - wz;
98   -
99   - result(0, 2) = xz + wy;
100   -
101   - result(1, 0) = xy + wz;
102   - result(1, 1) = (T)1.0 - (xx + zz);
103   -
104   - result(1, 2) = yz - wx;
105   -
106   - result(2, 0) = xz - wy;
107   - result(2, 1) = yz + wx;
108   -
109   - result(2, 2) = (T)1.0 - (xx + yy);
110   -
111   - return result;
112   -}
113   -
114   -template<typename T>
115   -matrix<T, 4> quaternion<T>::toMatrix4()
116   -{
117   - matrix<T, 4> result;
118   -
119   -
120   - T wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
121   -
122   -
123   - // calculate coefficients
124   - x2 = x + x; y2 = y + y;
125   - z2 = z + z;
126   - xx = x * x2; xy = x * y2; xz = x * z2;
127   - yy = y * y2; yz = y * z2; zz = z * z2;
128   - wx = w * x2; wy = w * y2; wz = w * z2;
129   -
130   - result(0, 0) = (T)1.0 - (yy + zz);
131   - result(0, 1) = xy - wz;
132   -
133   - result(0, 2) = xz + wy;
134   -
135   - result(1, 0) = xy + wz;
136   - result(1, 1) = (T)1.0 - (xx + zz);
137   -
138   - result(1, 2) = yz - wx;
139   -
140   - result(2, 0) = xz - wy;
141   - result(2, 1) = yz + wx;
142   -
143   - result(2, 2) = (T)1.0 - (xx + yy);
144   -
145   - result(3, 3) = (T)1.0;
146   -
147   - return result;
148   -}
149   -
150   -template<typename T>
151   -quaternion<T>::quaternion()
152   -{
153   - w=0.0; x=0.0; y=0.0; z=0.0;
154   -}
155   -
156   -template<typename T>
157   -quaternion<T>::quaternion(T c, T i, T j, T k)
158   -{
159   - w=c; x=i; y=j; z=k;
160   -}
161   -
162   -} //end rts namespace
163   -
164   -#if __GNUC__ > 3 && __GNUC_MINOR__ > 7
165   -template<class T> using rtsQuaternion = rts::quaternion<T>;
166   -#endif
167   -
168   -#endif
math/spherical_bessel.h
  Show/Hide comments   View file @7731cf2
... ... @@ -111,7 +111,7 @@ CUDA_CALLABLE void shift_sbesselj(int n, T x, T* b)//, T stability = 1.4)
111 111 //if(n > stability*x)
112 112 if(n > real(x))
113 113 if(real(bnew) < RTS_BESSEL_CONVERGENCE_MIN || real(bnew) > RTS_BESSEL_CONVERGENCE_MAX)
114   - bnew = 0.0;
  114 + bnew = 0;
115 115  
116 116 //shift and add the new value to the array
117 117 b[0] = b[1];
... ... @@ -130,8 +130,8 @@ CUDA_CALLABLE void shift_sbessely(int n, T x, T* b)//, T stability = 1.4)
130 130 if(bnew < RTS_BESSEL_MAXIMUM_FLOAT ||
131 131 (n > x && bnew > 0))
132 132 {
133   - bnew = (T)0;
134   - b[1] = (T)0;
  133 + bnew = 0;
  134 + b[1] = 0;
135 135 }
136 136  
137 137  
... ...
math/vector.h~ deleted
View file @57729e5
1   -#ifndef RTS_VECTOR_H
2   -#define RTS_VECTOR_H
3   -
4   -#include <iostream>
5   -#include <cmath>
6   -#include <sstream>
7   -//#include "rts/point.h"
8   -#include "rts/cuda/callable.h"
9   -
10   -namespace rts
11   -{
12   -
13   -
14   -
15   -template <class T, int N>
16   -struct vector
17   -{
18   - T v[N];
19   -
20   - CUDA_CALLABLE vector()
21   - {
22   - //memset(v, 0, sizeof(T) * N);
23   - for(int i=0; i<N; i++)
24   - v[i] = 0;
25   - }
26   -
27   - //efficiency constructor, makes construction easier for 1D-4D vectors
28   - CUDA_CALLABLE vector(T x, T y = (T)0.0, T z = (T)0.0, T w = (T)0.0)
29   - {
30   - if(N >= 1)
31   - v[0] = x;
32   - if(N >= 2)
33   - v[1] = y;
34   - if(N >= 3)
35   - v[2] = z;
36   - if(N >= 4)
37   - v[3] = w;
38   - }
39   -
40   - CUDA_CALLABLE vector(const T(&data)[N])
41   - {
42   - memcpy(v, data, sizeof(T) * N);
43   - }
44   -
45   - CUDA_CALLABLE T len()
46   - {
47   - //compute and return the vector length
48   - T sum_sq = (T)0;
49   - for(int i=0; i<N; i++)
50   - {
51   - sum_sq += v[i] * v[i];
52   - }
53   - return std::sqrt(sum_sq);
54   -
55   - }
56   -
57   - CUDA_CALLABLE vector<T, N> cart2sph()
58   - {
59   - //convert the vector from cartesian to spherical coordinates
60   - //x, y, z -> r, theta, phi (where theta = 0 to 2*pi)
61   -
62   - vector<T, N> sph;
63   - sph[0] = std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
64   - sph[1] = std::atan2(v[1], v[0]);
65   - sph[2] = std::acos(v[2] / sph[0]);
66   -
67   - return sph;
68   - }
69   -
70   - CUDA_CALLABLE vector<T, N> sph2cart()
71   - {
72   - //convert the vector from cartesian to spherical coordinates
73   - //r, theta, phi -> x, y, z (where theta = 0 to 2*pi)
74   -
75   - vector<T, N> cart;
76   - cart[0] = v[0] * std::cos(v[1]) * std::sin(v[2]);
77   - cart[1] = v[0] * std::sin(v[1]) * std::sin(v[2]);
78   - cart[2] = v[0] * std::cos(v[2]);
79   -
80   - return cart;
81   - }
82   -
83   - CUDA_CALLABLE vector<T, N> norm()
84   - {
85   - //compute and return the vector norm
86   - vector<T, N> result;
87   -
88   - //compute the vector length
89   - T l = len();
90   -
91   - //normalize
92   - for(int i=0; i<N; i++)
93   - {
94   - result.v[i] = v[i] / l;
95   - }
96   -
97   - return result;
98   - }
99   -
100   - CUDA_CALLABLE vector<T, 3> cross(vector<T, 3> rhs)
101   - {
102   - vector<T, 3> result;
103   -
104   - //compute the cross product (only valid for 3D vectors)
105   - result[0] = v[1] * rhs[2] - v[2] * rhs[1];
106   - result[1] = v[2] * rhs[0] - v[0] * rhs[2];
107   - result[2] = v[0] * rhs[1] - v[1] * rhs[0];
108   -
109   - return result;
110   - }
111   -
112   - CUDA_CALLABLE T dot(vector<T, N> rhs)
113   - {
114   - T result = (T)0;
115   -
116   - for(int i=0; i<N; i++)
117   - result += v[i] * rhs.v[i];
118   -
119   - return result;
120   -
121   - }
122   -
123   - //arithmetic
124   - CUDA_CALLABLE vector<T, N> operator+(vector<T, N> rhs)
125   - {
126   - vector<T, N> result;
127   -
128   - for(int i=0; i<N; i++)
129   - result.v[i] = v[i] + rhs.v[i];
130   -
131   - return result;
132   - }
133   - CUDA_CALLABLE vector<T, N> operator-(vector<T, N> rhs)
134   - {
135   - vector<T, N> result;
136   -
137   - for(int i=0; i<N; i++)
138   - result.v[i] = v[i] - rhs.v[i];
139   -
140   - return result;
141   - }
142   - CUDA_CALLABLE vector<T, N> operator*(T rhs)
143   - {
144   - vector<T, N> result;
145   -
146   - for(int i=0; i<N; i++)
147   - result.v[i] = v[i] * rhs;
148   -
149   - return result;
150   - }
151   - CUDA_CALLABLE vector<T, N> operator/(T rhs)
152   - {
153   - vector<T, N> result;
154   -
155   - for(int i=0; i<N; i++)
156   - result.v[i] = v[i] / rhs;
157   -
158   - return result;
159   - }
160   -
161   - CUDA_CALLABLE bool operator==(vector<T, N> rhs)
162   - {
163   - if ( (rhs.v[0] == v[0]) && (rhs.v[1] == v[1]) && (rhs.v[2] == v[2]) )
164   - return true;
165   -
166   - return false;
167   - }
168   -
169   - std::string toStr()
170   - {
171   - std::stringstream ss;
172   -
173   - ss<<"[";
174   - for(int i=0; i<N; i++)
175   - {
176   - ss<<v[i];
177   - if(i != N-1)
178   - ss<<", ";
179   - }
180   - ss<<"]";
181   -
182   - return ss.str();
183   - }
184   -
185   - //bracket operator
186   - CUDA_CALLABLE T& operator[](int i)
187   - {
188   - return v[i];
189   - }
190   -
191   -};
192   -
193   -
194   -} //end namespace rts
195   -
196   -template <typename T, int N>
197   -std::ostream& operator<<(std::ostream& os, rts::vector<T, N> v)
198   -{
199   - os<<v.toStr();
200   - return os;
201   -}
202   -
203   -//arithmetic operators
204   -template <typename T, int N>
205   -CUDA_CALLABLE rts::vector<T, N> operator-(rts::vector<T, N> v)
206   -{
207   - rts::vector<T, N> r;
208   -
209   - //negate the vector
210   - for(int i=0; i<N; i++)
211   - r.v[i] = -v.v[i];
212   -
213   - return r;
214   -}
215   -
216   -template <typename T, int N>
217   -CUDA_CALLABLE rts::vector<T, N> operator*(T lhs, rts::vector<T, N> rhs)
218   -{
219   - rts::vector<T, N> r;
220   -
221   - return rhs * lhs;
222   -}
223   -
224   -#if __GNUC__ > 3 && __GNUC_MINOR__ > 7
225   -template<class T, int N> using rtsVector = rts::vector<T, N>;
226   -#endif
227   -
228   -#endif
visualization/camera.h 100755 → 100644
  Show/Hide comments   View file @7731cf2