Stellarium 0.13.1
VecMath.hpp
1 /*
2  *
3  * Copyright (C) 2003 Fabien Chereau
4  *
5  * This program is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU General Public License
7  * as published by the Free Software Foundation; either version 2
8  * of the License, or (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA.
18  */
19 
20 // Template vector and matrix library.
21 // Use OpenGL compatible ordering ie. you can pass a matrix or vector to
22 // openGL functions without changes in the ordering
23 
24 #ifndef _VECMATH_H_
25 #define _VECMATH_H_
26 
27 #include "config.h"
28 
29 #include <cmath>
30 #include <QString>
31 
32 template<class T> class Vector2;
33 template<class T> class Vector3;
34 template<class T> class Vector4;
35 template<class T> class Matrix4;
36 
37 typedef Vector2<double> Vec2d;
38 typedef Vector2<float> Vec2f;
39 typedef Vector2<int> Vec2i;
40 
43 typedef Vector3<double> Vec3d;
44 
47 typedef Vector3<float> Vec3f;
48 
51 typedef Vector4<double> Vec4d;
52 
55 typedef Vector4<float> Vec4f;
56 
59 typedef Vector4<int> Vec4i;
60 
63 typedef Matrix4<double> Mat4d;
64 
67 typedef Matrix4<float> Mat4f;
68 
69 
73 template<class T> class Vector2
74 {
75 public:
76  inline Vector2();
77  inline Vector2(T, T);
78 
79  inline Vector2& operator=(const T*);
80  inline void set(T, T);
81 
82  inline bool operator==(const Vector2<T>&) const;
83  inline bool operator!=(const Vector2<T>&) const;
84 
85  inline const T& operator[](int x) const;
86  inline T& operator[](int);
87  inline operator const T*() const;
88  inline operator T*();
89 
90  inline void operator+=(const Vector2<T>&);
91  inline void operator-=(const Vector2<T>&);
92  inline void operator*=(T);
93  inline void operator/=(T);
94 
95  inline Vector2 operator-(const Vector2<T>&) const;
96  inline Vector2 operator+(const Vector2<T>&) const;
97 
98  inline Vector2 operator-() const;
99  inline Vector2 operator+() const;
100 
101  inline Vector2 operator*(T) const;
102  inline Vector2 operator/(T) const;
103 
104 
105  inline T dot(const Vector2<T>&) const;
106 
107  inline T length() const;
108  inline T lengthSquared() const;
109  inline void normalize();
110 
111  T v[2];
112 };
113 
114 
118 template<class T> class Vector3
119 {
120 public:
121  inline Vector3();
122  //inline Vector3(const Vector3&);
123  //template <class T2> inline Vector3(const Vector3<T2>&);
124  inline Vector3(T, T, T);
125  inline Vector3(T);
126 
127  //inline Vector3& operator=(const Vector3&);
128  inline Vector3& operator=(const T*);
129  //template <class T2> inline Vector3& operator=(const Vector3<T2>&);
130  inline void set(T, T, T);
131 
132  inline bool operator==(const Vector3<T>&) const;
133  inline bool operator!=(const Vector3<T>&) const;
134 
135  inline T& operator[](int);
136  inline const T& operator[](int) const;
137  inline operator const T*() const;
138  inline operator T*();
139  inline const T* data() const {return v;}
140  inline T* data() {return v;}
141 
142  inline void operator+=(const Vector3<T>&);
143  inline void operator-=(const Vector3<T>&);
144  inline void operator*=(T);
145  inline void operator/=(T);
146 
147  inline Vector3 operator-(const Vector3<T>&) const;
148  inline Vector3 operator+(const Vector3<T>&) const;
149 
150  inline Vector3 operator-() const;
151  inline Vector3 operator+() const;
152 
153  inline Vector3 operator*(T) const;
154  inline Vector3 operator/(T) const;
155 
156 
157  inline T dot(const Vector3<T>&) const;
158  inline Vector3 operator^(const Vector3<T>&) const;
159 
160  // Return latitude in rad
161  inline T latitude() const;
162  // Return longitude in rad
163  inline T longitude() const;
164 
165  // Distance in radian between two
166  inline T angle(const Vector3<T>&) const;
167  inline T angleNormalized(const Vector3<T>&) const;
168 
169  inline T length() const;
170  inline T lengthSquared() const;
171  inline void normalize();
172 
173  inline void transfo4d(const Mat4d&);
174  inline void transfo4d(const Mat4f&);
175  T v[3]; // The 3 values
176 
177  QString toString() const {return QString("[%1, %2, %3]").arg(v[0]).arg(v[1]).arg(v[2]);}
178  QString toStringLonLat() const {return QString("[") + QString::number(longitude()*180./M_PI, 'g', 12) + "," + QString::number(latitude()*180./M_PI, 'g', 12)+"]";}
179 };
180 
181 
185 template<class T> class Vector4
186 {
187 public:
188  inline Vector4();
189  inline Vector4(const Vector3<T>&);
190  inline Vector4(T, T, T);
191  inline Vector4(T, T, T, T);
192 
193  inline Vector4& operator=(const Vector3<T>&);
194  inline Vector4& operator=(const T*);
195  inline void set(T, T, T, T);
196 
197  inline bool operator==(const Vector4<T>&) const;
198  inline bool operator!=(const Vector4<T>&) const;
199 
200  inline T& operator[](int);
201  inline const T& operator[](int) const;
202  inline operator T*();
203  inline operator const T*() const;
204 
205  inline void operator+=(const Vector4<T>&);
206  inline void operator-=(const Vector4<T>&);
207  inline void operator*=(T);
208  inline void operator/=(T);
209 
210  inline Vector4 operator-(const Vector4<T>&) const;
211  inline Vector4 operator+(const Vector4<T>&) const;
212 
213  inline Vector4 operator-() const;
214  inline Vector4 operator+() const;
215 
216  inline Vector4 operator*(T) const;
217  inline Vector4 operator/(T) const;
218 
219 
220  inline T dot(const Vector4<T>&) const;
221 
222  inline T length() const;
223  inline T lengthSquared() const;
224  inline void normalize();
225 
226  inline void transfo4d(const Mat4d&);
227 
228  T v[4]; // The 4 values
229 };
230 
234 template<class T> class Matrix4
235 {
236  public:
237  Matrix4();
238  Matrix4(T,T,T,T,T,T,T,T,T,T,T,T,T,T,T,T);
239  Matrix4(const T*);
240  Matrix4(const Vector4<T>&, const Vector4<T>&,
241  const Vector4<T>&, const Vector4<T>&);
242 
243  inline Matrix4& operator=(const T*);
244  inline void set(T,T,T,T,T,T,T,T,T,T,T,T,T,T,T,T);
245 
246  inline T& operator[](int);
247  inline operator T*();
248  inline operator const T*() const;
249 
250  inline Matrix4 operator-(const Matrix4<T>&) const;
251  inline Matrix4 operator+(const Matrix4<T>&) const;
252  inline Matrix4 operator*(const Matrix4<T>&) const;
253 
254  inline Vector3<T> operator*(const Vector3<T>&) const;
255  inline Vector3<T> multiplyWithoutTranslation(const Vector3<T>& a) const;
256  inline Vector4<T> operator*(const Vector4<T>&) const;
257 
258  inline void transfo(Vector3<T>&) const;
259 
260  static Matrix4<T> identity();
261  static Matrix4<T> translation(const Vector3<T>&);
262 
263  // static Matrix4<T> rotation(const Vector3<T>&);
264  static Matrix4<T> rotation(const Vector3<T>&, T);
265  static Matrix4<T> xrotation(T);
266  static Matrix4<T> yrotation(T);
267  static Matrix4<T> zrotation(T);
268  static Matrix4<T> scaling(const Vector3<T>&);
269  static Matrix4<T> scaling(T);
270 
271  Matrix4<T> transpose() const;
272  Matrix4<T> inverse() const;
273 
274  inline Vector4<T> getRow(const int row) const;
275  inline Vector4<T> getColumn(const int column) const;
276 
277  inline void print(void) const;
278 
279  T r[16];
280 };
281 
283 template<class T> QDataStream& operator<<(QDataStream& out, const Vector2<T>& v) {out << v[0] << v[1]; return out;}
284 template<class T> QDataStream& operator<<(QDataStream& out, const Vector3<T>& v) {out << v[0] << v[1] << v[2]; return out;}
285 template<class T> QDataStream& operator<<(QDataStream& out, const Vector4<T>& v) {out << v[0] << v[1] << v[2] << v[3]; return out;}
286 template<class T> QDataStream& operator<<(QDataStream& out, const Matrix4<T>& m) {out << m[0] << m[1] << m[2] << m[3] << m[4] << m[5] << m[6] << m[7] << m[8] << m[9] << m[10] << m[11] << m[12] << m[13] << m[14] << m[15]; return out;}
287 
288 template<class T> QDataStream& operator>>(QDataStream& in, Vector2<T>& v) {in >> v[0] >> v[1]; return in;}
289 template<class T> QDataStream& operator>>(QDataStream& in, Vector3<T>& v) {in >> v[0] >> v[1] >> v[2]; return in;}
290 template<class T> QDataStream& operator>>(QDataStream& in, Vector4<T>& v) {in >> v[0] >> v[1] >> v[2] >> v[3]; return in;}
291 template<class T> QDataStream& operator>>(QDataStream& in, Matrix4<T>& m) {in >> m[0] >> m[1] >> m[2] >> m[3] >> m[4] >> m[5] >> m[6] >> m[7] >> m[8] >> m[9] >> m[10] >> m[11] >> m[12] >> m[13] >> m[14] >> m[15]; return in;}
292 
294 
295 template<class T> Vector2<T>::Vector2() {}
296 
297 template<class T> Vector2<T>::Vector2(T x, T y)
298 {
299  v[0]=x; v[1]=y;
300 }
301 
302 
303 template<class T> Vector2<T>& Vector2<T>::operator=(const T* a)
304 {
305  v[0]=a[0]; v[1]=a[1];
306  return *this;
307 }
308 
309 template<class T> void Vector2<T>::set(T x, T y)
310 {
311  v[0]=x; v[1]=y;
312 }
313 
314 
315 template<class T> bool Vector2<T>::operator==(const Vector2<T>& a) const
316 {
317  return (v[0] == a.v[0] && v[1] == a.v[1]);
318 }
319 
320 template<class T> bool Vector2<T>::operator!=(const Vector2<T>& a) const
321 {
322  return (v[0] != a.v[0] || v[1] != a.v[1]);
323 }
324 
325 template<class T> const T& Vector2<T>::operator[](int x) const
326 {
327  return v[x];
328 }
329 
330 template<class T> T& Vector2<T>::operator[](int x)
331 {
332  return v[x];
333 }
334 
335 template<class T> Vector2<T>::operator const T*() const
336 {
337  return v;
338 }
339 
340 template<class T> Vector2<T>::operator T*()
341 {
342  return v;
343 }
344 
345 
346 template<class T> void Vector2<T>::operator+=(const Vector2<T>& a)
347 {
348  v[0] += a.v[0]; v[1] += a.v[1];
349 }
350 
351 template<class T> void Vector2<T>::operator-=(const Vector2<T>& a)
352 {
353  v[0] -= a.v[0]; v[1] -= a.v[1];
354 }
355 
356 template<class T> void Vector2<T>::operator*=(T s)
357 {
358  v[0] *= s; v[1] *= s;
359 }
360 
361 template<class T> void Vector2<T>::operator/=(T s)
362 {
363  v[0] /= s; v[1] /= s;
364 }
365 
366 template<class T> Vector2<T> Vector2<T>::operator-() const
367 {
368  return Vector2<T>(-v[0], -v[1]);
369 }
370 
371 template<class T> Vector2<T> Vector2<T>::operator+() const
372 {
373  return *this;
374 }
375 
376 template<class T> Vector2<T> Vector2<T>::operator+(const Vector2<T>& b) const
377 {
378  return Vector2<T>(v[0] + b.v[0], v[1] + b.v[1]);
379 }
380 
381 template<class T> Vector2<T> Vector2<T>::operator-(const Vector2<T>& b) const
382 {
383  return Vector2<T>(v[0] - b.v[0], v[1] - b.v[1]);
384 }
385 
386 template<class T> Vector2<T> Vector2<T>::operator*(T s) const
387 {
388  return Vector2<T>(s * v[0], s * v[1]);
389 }
390 
391 template<class T> Vector2<T> Vector2<T>::operator/(T s) const
392 {
393  return Vector2<T>(v[0]/s, v[1]/s);
394 }
395 
396 
397 template<class T> T Vector2<T>::dot(const Vector2<T>& b) const
398 {
399  return v[0] * b.v[0] + v[1] * b.v[1];
400 }
401 
402 
403 template<class T> T Vector2<T>::length() const
404 {
405  return (T) std::sqrt(v[0] * v[0] + v[1] * v[1]);
406 }
407 
408 template<class T> T Vector2<T>::lengthSquared() const
409 {
410  return v[0] * v[0] + v[1] * v[1];
411 }
412 
413 template<class T> void Vector2<T>::normalize()
414 {
415  T s = (T) 1 / std::sqrt(v[0] * v[0] + v[1] * v[1]);
416  v[0] *= s;
417  v[1] *= s;
418 }
419 
420 // template<class T>
421 // std::ostream& operator<<(std::ostream &o,const Vector2<T> &v) {
422 // return o << '[' << v[0] << ',' << v[1] << ']';
423 // }
424 
426 
427 template<class T> Vector3<T>::Vector3() {}
428 
429 //template<class T> Vector3<T>::Vector3(const Vector3& a)
430 //{
431 // v[0]=a.v[0]; v[1]=a.v[1]; v[2]=a.v[2];
432 //}
433 
434 //template<class T> template<class T2> Vector3<T>::Vector3(const Vector3<T2>& a)
435 //{
436 // v[0]=a.v[0]; v[1]=a.v[1]; v[2]=a.v[2];
437 //}
438 
439 template<class T> Vector3<T>::Vector3(T x)
440 {
441  v[0]=x; v[1]=x; v[2]=x;
442 }
443 
444 template<class T> Vector3<T>::Vector3(T x, T y, T z)
445 {
446  v[0]=x; v[1]=y; v[2]=z;
447 }
448 
449 //template<class T> Vector3<T>& Vector3<T>::operator=(const Vector3& a)
450 //{
451 // v[0]=a.v[0]; v[1]=a.v[1]; v[2]=a.v[2];
452 // return *this;
453 //}
454 
455 //template<class T> template <class T2> Vector3<T>& Vector3<T>::operator=(const Vector3<T2>& a)
456 //{
457 // v[0]=a.v[0]; v[1]=a.v[1]; v[2]=a.v[2];
458 // return *this;
459 //}
460 
461 template<class T> Vector3<T>& Vector3<T>::operator=(const T* a)
462 {
463  v[0]=a[0]; v[1]=a[1]; v[2]=a[2];
464  return *this;
465 }
466 
467 template<class T> void Vector3<T>::set(T x, T y, T z)
468 {
469  v[0]=x; v[1]=y; v[2]=z;
470 }
471 
472 
473 template<class T> bool Vector3<T>::operator==(const Vector3<T>& a) const
474 {
475  return (v[0] == a.v[0] && v[1] == a.v[1] && v[2] == a.v[2]);
476 }
477 
478 template<class T> bool Vector3<T>::operator!=(const Vector3<T>& a) const
479 {
480  return (v[0] != a.v[0] || v[1] != a.v[1] || v[2] != a.v[2]);
481 }
482 
483 
484 template<class T> T& Vector3<T>::operator[](int x)
485 {
486  return v[x];
487 }
488 
489 template<class T> const T& Vector3<T>::operator[](int x) const
490 {
491  return v[x];
492 }
493 
494 template<class T> Vector3<T>::operator const T*() const
495 {
496  return v;
497 }
498 
499 template<class T> Vector3<T>::operator T*()
500 {
501  return v;
502 }
503 
504 template<class T> void Vector3<T>::operator+=(const Vector3<T>& a)
505 {
506  v[0] += a.v[0]; v[1] += a.v[1]; v[2] += a.v[2];
507 }
508 
509 template<class T> void Vector3<T>::operator-=(const Vector3<T>& a)
510 {
511  v[0] -= a.v[0]; v[1] -= a.v[1]; v[2] -= a.v[2];
512 }
513 
514 template<class T> void Vector3<T>::operator*=(T s)
515 {
516  v[0] *= s; v[1] *= s; v[2] *= s;
517 }
518 
519 template<class T> void Vector3<T>::operator/=(T s)
520 {
521  v[0] /= s; v[1] /= s; v[2] /= s;
522 }
523 
524 template<class T> Vector3<T> Vector3<T>::operator-() const
525 {
526  return Vector3<T>(-v[0], -v[1], -v[2]);
527 }
528 
529 template<class T> Vector3<T> Vector3<T>::operator+() const
530 {
531  return *this;
532 }
533 
534 template<class T> Vector3<T> Vector3<T>::operator+(const Vector3<T>& b) const
535 {
536  return Vector3<T>(v[0] + b.v[0], v[1] + b.v[1], v[2] + b.v[2]);
537 }
538 
539 template<class T> Vector3<T> Vector3<T>::operator-(const Vector3<T>& b) const
540 {
541  return Vector3<T>(v[0] - b.v[0], v[1] - b.v[1], v[2] - b.v[2]);
542 }
543 
544 template<class T> Vector3<T> Vector3<T>::operator*(T s) const
545 {
546  return Vector3<T>(s * v[0], s * v[1], s * v[2]);
547 }
548 
549 template<class T> Vector3<T> Vector3<T>::operator/(T s) const
550 {
551  return Vector3<T>(v[0]/s, v[1]/s, v[2]/s);
552 }
553 
554 
555 template<class T> T Vector3<T>::dot(const Vector3<T>& b) const
556 {
557  return v[0] * b.v[0] + v[1] * b.v[1] + v[2] * b.v[2];
558 }
559 
560 
561 // cross product
562 template<class T> Vector3<T> Vector3<T>::operator^(const Vector3<T>& b) const
563 {
564  return Vector3<T>(v[1] * b.v[2] - v[2] * b.v[1],
565  v[2] * b.v[0] - v[0] * b.v[2],
566  v[0] * b.v[1] - v[1] * b.v[0]);
567 }
568 
569 // Angle in radian between two vectors
570 template<class T> T Vector3<T>::angle(const Vector3<T>& b) const
571 {
572  const T cosAngle = dot(b)/sqrt(lengthSquared()*b.lengthSquared());
573  return cosAngle>=1 ? 0 : (cosAngle<=-1 ? M_PI : std::acos(cosAngle));
574 }
575 
576 // Angle in radian between two normalized vectors
577 template<class T> T Vector3<T>::angleNormalized(const Vector3<T>& b) const
578 {
579  const T cosAngle = dot(b);
580  return cosAngle>=1 ? 0 : (cosAngle<=-1 ? M_PI : std::acos(cosAngle));
581 }
582 
583 template<class T> T Vector3<T>::length() const
584 {
585  return (T) sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
586 }
587 
588 template<class T> T Vector3<T>::lengthSquared() const
589 {
590  return v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
591 }
592 
593 template<class T> void Vector3<T>::normalize()
594 {
595  T s = (T) (1. / std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]));
596  v[0] *= s;
597  v[1] *= s;
598  v[2] *= s;
599 }
600 
601 template<class T> void Vector3<T>::transfo4d(const Mat4d& m)
602 {
603  const T v0 = v[0];
604  const T v1 = v[1];
605  v[0]=m.r[0]*v0 + m.r[4]*v1 + m.r[8]*v[2] + m.r[12];
606  v[1]=m.r[1]*v0 + m.r[5]*v1 + m.r[9]*v[2] + m.r[13];
607  v[2]=m.r[2]*v0 + m.r[6]*v1 + m.r[10]*v[2] + m.r[14];
608 }
609 
610 template<class T> void Vector3<T>::transfo4d(const Mat4f& m)
611 {
612  const T v0 = v[0];
613  const T v1 = v[1];
614  v[0]=m.r[0]*v0 + m.r[4]*v1 + m.r[8]*v[2] + m.r[12];
615  v[1]=m.r[1]*v0 + m.r[5]*v1 + m.r[9]*v[2] + m.r[13];
616  v[2]=m.r[2]*v0 + m.r[6]*v1 + m.r[10]*v[2] + m.r[14];
617 }
618 
619 // Return latitude in rad
620 template<class T> T Vector3<T>::latitude() const
621 {
622  return std::asin(v[2]/length());
623 }
624 
625 // Return longitude in rad
626 template<class T> T Vector3<T>::longitude() const
627 {
628  return std::atan2(v[1],v[0]);
629 }
630 
631 
633 
634 template<class T> Vector4<T>::Vector4() {}
635 
636 template<class T> Vector4<T>::Vector4(const Vector3<T>& a)
637 {
638  v[0]=a.v[0]; v[1]=a.v[1]; v[2]=a.v[2]; v[3]=1;
639 }
640 
641 template<class T> Vector4<T>::Vector4(T x, T y, T z)
642 {
643  v[0]=x; v[1]=y; v[2]=z; v[3]=1;
644 }
645 
646 template<class T> Vector4<T>::Vector4(T x, T y, T z, T a)
647 {
648  v[0]=x; v[1]=y; v[2]=z; v[3]=a;
649 }
650 
651 template<class T> Vector4<T>& Vector4<T>::operator=(const Vector3<T>& a)
652 {
653  v[0]=a.v[0]; v[1]=a.v[1]; v[2]=a.v[2]; v[3]=1;
654  return *this;
655 }
656 
657 template<class T> Vector4<T>& Vector4<T>::operator=(const T* a)
658 {
659  v[0]=a[0]; v[1]=a[1]; v[2]=a[2]; v[3]=a[3];
660  return *this;
661 }
662 
663 template<class T> void Vector4<T>::set(T x, T y, T z, T a)
664 {
665  v[0]=x; v[1]=y; v[2]=z; v[3]=a;
666 }
667 
668 template<class T> bool Vector4<T>::operator==(const Vector4<T>& a) const
669 {
670  return (v[0] == a.v[0] && v[1] == a.v[1] && v[2] == a.v[2] && v[3] == a.v[3]);
671 }
672 
673 template<class T> bool Vector4<T>::operator!=(const Vector4<T>& a) const
674 {
675  return (v[0] != a.v[0] || v[1] != a.v[1] || v[2] != a.v[2] || v[3] != a.v[3]);
676 }
677 
678 template<class T> T& Vector4<T>::operator[](int x)
679 {
680  return v[x];
681 }
682 
683 template<class T> const T& Vector4<T>::operator[](int x) const
684 {
685  return v[x];
686 }
687 
688 template<class T> Vector4<T>::operator T*()
689 {
690  return v;
691 }
692 
693 template<class T> Vector4<T>::operator const T*() const
694 {
695  return v;
696 }
697 
698 template<class T> void Vector4<T>::operator+=(const Vector4<T>& a)
699 {
700  v[0] += a.v[0]; v[1] += a.v[1]; v[2] += a.v[2]; v[3] += a.v[3];
701 }
702 
703 template<class T> void Vector4<T>::operator-=(const Vector4<T>& a)
704 {
705  v[0] -= a.v[0]; v[1] -= a.v[1]; v[2] -= a.v[2]; v[3] -= a.v[3];
706 }
707 
708 template<class T> void Vector4<T>::operator*=(T s)
709 {
710  v[0] *= s; v[1] *= s; v[2] *= s; v[3] *= s;
711 }
712 
713 template<class T> void Vector4<T>::operator/=(T s)
714 {
715  v[0] /= s; v[1] /= s; v[2] /= s; v[3] /= s;
716 }
717 
718 template<class T> Vector4<T> Vector4<T>::operator-() const
719 {
720  return Vector4<T>(-v[0], -v[1], -v[2], -v[3]);
721 }
722 
723 template<class T> Vector4<T> Vector4<T>::operator+() const
724 {
725  return *this;
726 }
727 
728 template<class T> Vector4<T> Vector4<T>::operator+(const Vector4<T>& b) const
729 {
730  return Vector4<T>(v[0] + b.v[0], v[1] + b.v[1], v[2] + b.v[2], v[3] + b.v[3]);
731 }
732 
733 template<class T> Vector4<T> Vector4<T>::operator-(const Vector4<T>& b) const
734 {
735  return Vector4<T>(v[0] - b.v[0], v[1] - b.v[1], v[2] - b.v[2], v[3] - b.v[3]);
736 }
737 
738 template<class T> Vector4<T> Vector4<T>::operator*(T s) const
739 {
740  return Vector4<T>(s * v[0], s * v[1], s * v[2], s * v[3]);
741 }
742 
743 template<class T> Vector4<T> Vector4<T>::operator/(T s) const
744 {
745  return Vector4<T>(v[0]/s, v[1]/s, v[2]/s, v[3]/s);
746 }
747 
748 template<class T> T Vector4<T>::dot(const Vector4<T>& b) const
749 {
750  return v[0] * b.v[0] + v[1] * b.v[1] + v[2] * b.v[2] + v[3] * b.v[3];
751 }
752 
753 template<class T> T Vector4<T>::length() const
754 {
755  return (T) sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2] + v[3] * v[3]);
756 }
757 
758 template<class T> T Vector4<T>::lengthSquared() const
759 {
760  return v[0] * v[0] + v[1] * v[1] + v[2] * v[2] + v[3] * v[3];
761 }
762 
763 template<class T> void Vector4<T>::normalize()
764 {
765  T s = (T) (1. / sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2] + v[3] * v[3]));
766  v[0] *= s;
767  v[1] *= s;
768  v[2] *= s;
769  v[3] *= s;
770 }
771 
772 template<class T> void Vector4<T>::transfo4d(const Mat4d& m)
773 {
774  (*this)=m*(*this);
775 }
776 /*
777 template<class T>
778 std::ostream& operator<<(std::ostream &o,const Vector4<T> &v) {
779  return o << '[' << v[0] << ',' << v[1] << ',' << v[2] << ',' << v[3] << ']';
780 }*/
781 
782 
784 
785 template<class T> Matrix4<T>::Matrix4() {}
786 
787 template<class T> Matrix4<T>::Matrix4(const T* m)
788 {
789  memcpy(r,m,sizeof(T)*16);
790 }
791 
792 template<class T> Matrix4<T>::Matrix4(const Vector4<T>& v0,
793  const Vector4<T>& v1,
794  const Vector4<T>& v2,
795  const Vector4<T>& v3)
796 {
797  r[0] = v0.v[0];
798  r[1] = v0.v[1];
799  r[2] = v0.v[2];
800  r[3] = v0.v[3];
801  r[4] = v1.v[0];
802  r[5] = v1.v[1];
803  r[6] = v1.v[2];
804  r[7] = v1.v[3];
805  r[8] = v2.v[0];
806  r[9] = v2.v[1];
807  r[10] = v2.v[2];
808  r[11] = v2.v[3];
809  r[12] = v3.v[0];
810  r[13] = v3.v[1];
811  r[14] = v3.v[2];
812  r[15] = v3.v[3];
813 }
814 
815 
816 template<class T> Matrix4<T>::Matrix4(T a, T b, T c, T d, T e, T f, T g, T h, T i, T j, T k, T l, T m, T n, T o, T p)
817 {
818  r[0]=a; r[1]=b; r[2]=c; r[3]=d; r[4]=e; r[5]=f; r[6]=g; r[7]=h;
819  r[8]=i; r[9]=j; r[10]=k; r[11]=l; r[12]=m; r[13]=n; r[14]=o; r[15]=p;
820 }
821 
822 template<class T> void Matrix4<T>::set(T a, T b, T c, T d, T e, T f, T g, T h, T i, T j, T k, T l, T m, T n, T o, T p)
823 {
824  r[0]=a; r[1]=b; r[2]=c; r[3]=d; r[4]=e; r[5]=f; r[6]=g; r[7]=h;
825  r[8]=i; r[9]=j; r[10]=k; r[11]=l; r[12]=m; r[13]=n; r[14]=o; r[15]=p;
826 }
827 
828 template<class T> T& Matrix4<T>::operator[](int n)
829 {
830  return r[n];
831 }
832 
833 template<class T> Matrix4<T>::operator T*()
834 {
835  return r;
836 }
837 
838 template<class T> Matrix4<T>::operator const T*() const
839 {
840  return r;
841 }
842 
843 template<class T> Matrix4<T> Matrix4<T>::identity()
844 {
845  return Matrix4<T>( 1, 0, 0, 0,
846  0, 1, 0, 0,
847  0, 0, 1, 0,
848  0, 0, 0, 1 );
849 }
850 
851 
852 template<class T> Matrix4<T> Matrix4<T>::translation(const Vector3<T>& a)
853 {
854  return Matrix4<T>( 1, 0, 0, 0,
855  0, 1, 0, 0,
856  0, 0, 1, 0,
857  a.v[0], a.v[1], a.v[2], 1);
858 }
859 
860 
861 template<class T> Matrix4<T> Matrix4<T>::rotation(const Vector3<T>& axis, T angle)
862 {
863  Vector3<T> a(axis);
864  a.normalize();
865  const T c = (T) cos(angle);
866  const T s = (T) sin(angle);
867  const T d = 1-c;
868  return Matrix4<T>( a[0]*a[0]*d+c , a[1]*a[0]*d+a[2]*s, a[0]*a[2]*d-a[1]*s, 0,
869  a[0]*a[1]*d-a[2]*s, a[1]*a[1]*d+c , a[1]*a[2]*d+a[0]*s, 0,
870  a[0]*a[2]*d+a[1]*s, a[1]*a[2]*d-a[0]*s, a[2]*a[2]*d+c , 0,
871  0,0,0,1 );
872 }
873 
874 template<class T> Matrix4<T> Matrix4<T>::xrotation(T angle)
875 {
876  T c = (T) cos(angle);
877  T s = (T) sin(angle);
878 
879  return Matrix4<T>(1, 0, 0, 0,
880  0, c, s, 0,
881  0,-s, c, 0,
882  0, 0, 0, 1 );
883 }
884 
885 
886 template<class T> Matrix4<T> Matrix4<T>::yrotation(T angle)
887 {
888  T c = (T) cos(angle);
889  T s = (T) sin(angle);
890 
891  return Matrix4<T>( c, 0,-s, 0,
892  0, 1, 0, 0,
893  s, 0, c, 0,
894  0, 0, 0, 1 );
895 }
896 
897 
898 template<class T> Matrix4<T> Matrix4<T>::zrotation(T angle)
899 {
900  T c = (T) cos(angle);
901  T s = (T) sin(angle);
902 
903  return Matrix4<T>(c, s, 0, 0,
904  -s, c, 0, 0,
905  0, 0, 1, 0,
906  0, 0, 0, 1 );
907 }
908 
909 
910 template<class T> Matrix4<T> Matrix4<T>::scaling(const Vector3<T>& s)
911 {
912  return Matrix4<T>(s[0], 0 , 0 , 0,
913  0 ,s[1], 0 , 0,
914  0 , 0 ,s[2], 0,
915  0 , 0 , 0 , 1);
916 }
917 
918 
919 template<class T> Matrix4<T> Matrix4<T>::scaling(T scale)
920 {
921  return scaling(Vector3<T>(scale, scale, scale));
922 }
923 
924 // multiply column vector by a 4x4 matrix in homogeneous coordinate (use a[3]=1)
925 template<class T> Vector3<T> Matrix4<T>::operator*(const Vector3<T>& a) const
926 {
927  return Vector3<T>( r[0]*a.v[0] + r[4]*a.v[1] + r[8]*a.v[2] + r[12],
928  r[1]*a.v[0] + r[5]*a.v[1] + r[9]*a.v[2] + r[13],
929  r[2]*a.v[0] + r[6]*a.v[1] + r[10]*a.v[2] + r[14] );
930 }
931 
932 template<class T> Vector3<T> Matrix4<T>::multiplyWithoutTranslation(const Vector3<T>& a) const
933 {
934  return Vector3<T>( r[0]*a.v[0] + r[4]*a.v[1] + r[8]*a.v[2],
935  r[1]*a.v[0] + r[5]*a.v[1] + r[9]*a.v[2],
936  r[2]*a.v[0] + r[6]*a.v[1] + r[10]*a.v[2] );
937 }
938 
939 // multiply column vector by a 4x4 matrix in homogeneous coordinate (considere a[3]=1)
940 template<class T> Vector4<T> Matrix4<T>::operator*(const Vector4<T>& a) const
941 {
942  return Vector4<T>( r[0]*a.v[0] + r[4]*a.v[1] + r[8]*a.v[2] + r[12]*a.v[3],
943  r[1]*a.v[0] + r[5]*a.v[1] + r[9]*a.v[2] + r[13]*a.v[3],
944  r[2]*a.v[0] + r[6]*a.v[1] + r[10]*a.v[2] + r[14]*a.v[3] );
945 }
946 
947 template<class T> void Matrix4<T>::transfo(Vector3<T>& a) const
948 {
949  a.set( r[0]*a.v[0] + r[4]*a.v[1] + r[8]*a.v[2] + r[12],
950  r[1]*a.v[0] + r[5]*a.v[1] + r[9]*a.v[2] + r[13],
951  r[2]*a.v[0] + r[6]*a.v[1] + r[10]*a.v[2] + r[14]);
952 }
953 
954 template<class T> Matrix4<T> Matrix4<T>::transpose() const
955 {
956  return Matrix4<T>( r[0], r[4], r[8], r[12],
957  r[1], r[5], r[9], r[13],
958  r[2], r[6], r[10], r[14],
959  r[3], r[7], r[11], r[15]);
960 }
961 
962 template<class T> Matrix4<T> Matrix4<T>::operator*(const Matrix4<T>& a) const
963 {
964 #define MATMUL(R, C) (r[R] * a.r[C] + r[R+4] * a.r[C+1] + r[R+8] * a.r[C+2] + r[R+12] * a.r[C+3])
965  return Matrix4<T>( MATMUL(0,0), MATMUL(1,0), MATMUL(2,0), MATMUL(3,0),
966  MATMUL(0,4), MATMUL(1,4), MATMUL(2,4), MATMUL(3,4),
967  MATMUL(0,8), MATMUL(1,8), MATMUL(2,8), MATMUL(3,8),
968  MATMUL(0,12), MATMUL(1,12), MATMUL(2,12), MATMUL(3,12) );
969 #undef MATMUL
970 }
971 
972 
973 template<class T> Matrix4<T> Matrix4<T>::operator+(const Matrix4<T>& a) const
974 {
975  return Matrix4<T>( r[0]+a.r[0], r[1]+a.r[1], r[2]+a.r[2], r[3]+a.r[3],
976  r[4]+a.r[4], r[5]+a.r[5], r[6]+a.r[6], r[7]+a.r[7],
977  r[8]+a.r[8], r[9]+a.r[9], r[10]+a.r[10], r[11]+a.r[11],
978  r[12]+a.r[12], r[13]+a.r[13], r[14]+a.r[14], r[15]+a.r[15] );
979 }
980 
981 template<class T> Matrix4<T> Matrix4<T>::operator-(const Matrix4<T>& a) const
982 {
983  return Matrix4<T>( r[0]-a.r[0], r[1]-a.r[1], r[2]-a.r[2], r[3]-a.r[3],
984  r[4]-a.r[4], r[5]-a.r[5], r[6]-a.r[6], r[7]-a.r[7],
985  r[8]-a.r[8], r[9]-a.r[9], r[10]-a.r[10], r[11]-a.r[11],
986  r[12]-a.r[12], r[13]-a.r[13], r[14]-a.r[14], r[15]-a.r[15] );
987 }
988 
989 /*
990  * Code ripped from the GLU library
991  * Compute inverse of 4x4 transformation matrix.
992  * Code contributed by Jacques Leroy jle@star.be
993  * Return zero matrix on failure (singular matrix)
994  */
995 template<class T> Matrix4<T> Matrix4<T>::inverse() const
996 {
997  const T * m = r;
998  T out[16];
999 
1000 /* NB. OpenGL Matrices are COLUMN major. */
1001 #define SWAP_ROWS(a, b) { T *_tmp = a; (a)=(b); (b)=_tmp; }
1002 #define MAT(m,r,c) (m)[(c)*4+(r)]
1003 
1004  T wtmp[4][8];
1005  T m0, m1, m2, m3, s;
1006  T *r0, *r1, *r2, *r3;
1007 
1008  r0 = wtmp[0], r1 = wtmp[1], r2 = wtmp[2], r3 = wtmp[3];
1009 
1010  r0[0] = MAT(m, 0, 0), r0[1] = MAT(m, 0, 1),
1011  r0[2] = MAT(m, 0, 2), r0[3] = MAT(m, 0, 3),
1012  r0[4] = 1.0, r0[5] = r0[6] = r0[7] = 0.0,
1013  r1[0] = MAT(m, 1, 0), r1[1] = MAT(m, 1, 1),
1014  r1[2] = MAT(m, 1, 2), r1[3] = MAT(m, 1, 3),
1015  r1[5] = 1.0, r1[4] = r1[6] = r1[7] = 0.0,
1016  r2[0] = MAT(m, 2, 0), r2[1] = MAT(m, 2, 1),
1017  r2[2] = MAT(m, 2, 2), r2[3] = MAT(m, 2, 3),
1018  r2[6] = 1.0, r2[4] = r2[5] = r2[7] = 0.0,
1019  r3[0] = MAT(m, 3, 0), r3[1] = MAT(m, 3, 1),
1020  r3[2] = MAT(m, 3, 2), r3[3] = MAT(m, 3, 3),
1021  r3[7] = 1.0, r3[4] = r3[5] = r3[6] = 0.0;
1022 
1023  /* choose pivot - or die */
1024  if (fabs(r3[0]) > fabs(r2[0]))
1025  SWAP_ROWS(r3, r2);
1026  if (fabs(r2[0]) > fabs(r1[0]))
1027  SWAP_ROWS(r2, r1);
1028  if (fabs(r1[0]) > fabs(r0[0]))
1029  SWAP_ROWS(r1, r0);
1030  if (0.0 == r0[0])
1031  return Matrix4<T>();
1032 
1033  /* eliminate first variable */
1034  m1 = r1[0] / r0[0];
1035  m2 = r2[0] / r0[0];
1036  m3 = r3[0] / r0[0];
1037  s = r0[1];
1038  r1[1] -= m1 * s;
1039  r2[1] -= m2 * s;
1040  r3[1] -= m3 * s;
1041  s = r0[2];
1042  r1[2] -= m1 * s;
1043  r2[2] -= m2 * s;
1044  r3[2] -= m3 * s;
1045  s = r0[3];
1046  r1[3] -= m1 * s;
1047  r2[3] -= m2 * s;
1048  r3[3] -= m3 * s;
1049  s = r0[4];
1050  if (s != 0.0) {
1051  r1[4] -= m1 * s;
1052  r2[4] -= m2 * s;
1053  r3[4] -= m3 * s;
1054  }
1055  s = r0[5];
1056  if (s != 0.0) {
1057  r1[5] -= m1 * s;
1058  r2[5] -= m2 * s;
1059  r3[5] -= m3 * s;
1060  }
1061  s = r0[6];
1062  if (s != 0.0) {
1063  r1[6] -= m1 * s;
1064  r2[6] -= m2 * s;
1065  r3[6] -= m3 * s;
1066  }
1067  s = r0[7];
1068  if (s != 0.0) {
1069  r1[7] -= m1 * s;
1070  r2[7] -= m2 * s;
1071  r3[7] -= m3 * s;
1072  }
1073 
1074  /* choose pivot - or die */
1075  if (fabs(r3[1]) > fabs(r2[1]))
1076  SWAP_ROWS(r3, r2);
1077  if (fabs(r2[1]) > fabs(r1[1]))
1078  SWAP_ROWS(r2, r1);
1079  if (0.0 == r1[1])
1080  return Matrix4<T>();
1081 
1082  /* eliminate second variable */
1083  m2 = r2[1] / r1[1];
1084  m3 = r3[1] / r1[1];
1085  r2[2] -= m2 * r1[2];
1086  r3[2] -= m3 * r1[2];
1087  r2[3] -= m2 * r1[3];
1088  r3[3] -= m3 * r1[3];
1089  s = r1[4];
1090  if (0.0 != s) {
1091  r2[4] -= m2 * s;
1092  r3[4] -= m3 * s;
1093  }
1094  s = r1[5];
1095  if (0.0 != s) {
1096  r2[5] -= m2 * s;
1097  r3[5] -= m3 * s;
1098  }
1099  s = r1[6];
1100  if (0.0 != s) {
1101  r2[6] -= m2 * s;
1102  r3[6] -= m3 * s;
1103  }
1104  s = r1[7];
1105  if (0.0 != s) {
1106  r2[7] -= m2 * s;
1107  r3[7] -= m3 * s;
1108  }
1109 
1110  /* choose pivot - or die */
1111  if (fabs(r3[2]) > fabs(r2[2]))
1112  SWAP_ROWS(r3, r2);
1113  if (0.0 == r2[2])
1114  return Matrix4<T>();
1115 
1116  /* eliminate third variable */
1117  m3 = r3[2] / r2[2];
1118  r3[3] -= m3 * r2[3], r3[4] -= m3 * r2[4],
1119  r3[5] -= m3 * r2[5], r3[6] -= m3 * r2[6], r3[7] -= m3 * r2[7];
1120 
1121  /* last check */
1122  if (0.0 == r3[3])
1123  return Matrix4<T>();
1124 
1125  s = 1.0 / r3[3]; /* now back substitute row 3 */
1126  r3[4] *= s;
1127  r3[5] *= s;
1128  r3[6] *= s;
1129  r3[7] *= s;
1130 
1131  m2 = r2[3]; /* now back substitute row 2 */
1132  s = 1.0 / r2[2];
1133  r2[4] = s * (r2[4] - r3[4] * m2), r2[5] = s * (r2[5] - r3[5] * m2),
1134  r2[6] = s * (r2[6] - r3[6] * m2), r2[7] = s * (r2[7] - r3[7] * m2);
1135  m1 = r1[3];
1136  r1[4] -= r3[4] * m1, r1[5] -= r3[5] * m1,
1137  r1[6] -= r3[6] * m1, r1[7] -= r3[7] * m1;
1138  m0 = r0[3];
1139  r0[4] -= r3[4] * m0, r0[5] -= r3[5] * m0,
1140  r0[6] -= r3[6] * m0, r0[7] -= r3[7] * m0;
1141 
1142  m1 = r1[2]; /* now back substitute row 1 */
1143  s = 1.0 / r1[1];
1144  r1[4] = s * (r1[4] - r2[4] * m1), r1[5] = s * (r1[5] - r2[5] * m1),
1145  r1[6] = s * (r1[6] - r2[6] * m1), r1[7] = s * (r1[7] - r2[7] * m1);
1146  m0 = r0[2];
1147  r0[4] -= r2[4] * m0, r0[5] -= r2[5] * m0,
1148  r0[6] -= r2[6] * m0, r0[7] -= r2[7] * m0;
1149 
1150  m0 = r0[1]; /* now back substitute row 0 */
1151  s = 1.0 / r0[0];
1152  r0[4] = s * (r0[4] - r1[4] * m0), r0[5] = s * (r0[5] - r1[5] * m0),
1153  r0[6] = s * (r0[6] - r1[6] * m0), r0[7] = s * (r0[7] - r1[7] * m0);
1154 
1155  MAT(out, 0, 0) = r0[4];
1156  MAT(out, 0, 1) = r0[5], MAT(out, 0, 2) = r0[6];
1157  MAT(out, 0, 3) = r0[7], MAT(out, 1, 0) = r1[4];
1158  MAT(out, 1, 1) = r1[5], MAT(out, 1, 2) = r1[6];
1159  MAT(out, 1, 3) = r1[7], MAT(out, 2, 0) = r2[4];
1160  MAT(out, 2, 1) = r2[5], MAT(out, 2, 2) = r2[6];
1161  MAT(out, 2, 3) = r2[7], MAT(out, 3, 0) = r3[4];
1162  MAT(out, 3, 1) = r3[5], MAT(out, 3, 2) = r3[6];
1163  MAT(out, 3, 3) = r3[7];
1164 
1165  return Matrix4<T>(out);
1166 
1167 #undef MAT
1168 #undef SWAP_ROWS
1169 }
1170 
1171 template<class T> Vector4<T> Matrix4<T>::getRow(const int row) const
1172 {
1173  return Vector4<T>(r[0 + row], r[4 + row], r[8 + row], r[12 + row]);
1174 }
1175 
1176 template<class T> Vector4<T> Matrix4<T>::getColumn(const int column) const
1177 {
1178  return Vector4<T>(r[0 + column * 4], r[1 + column * 4], r[2 + column * 4], r[3 + column * 4]);
1179 }
1180 
1181 template<class T> void Matrix4<T>::print(void) const
1182 {
1183  printf("[%5.2lf %5.2lf %5.2lf %17.12le]\n"
1184  "[%5.2lf %5.2lf %5.2lf %17.12le]\n"
1185  "[%5.2lf %5.2lf %5.2lf %17.12le]\n"
1186  "[%5.2lf %5.2lf %5.2lf %17.12le]\n\n",
1187  r[0],r[4],r[8],r[12],
1188  r[1],r[5],r[9],r[13],
1189  r[2],r[6],r[10],r[14],
1190  r[3],r[7],r[11],r[15]);
1191 }
1192 
1193 
1194 template<class T> inline
1195 T operator*(const Vector2<T>&a,const Vector2<T>&b) {
1196  return a.v[0] * b.v[0] + a.v[1] * b.v[1];
1197 }
1198 
1199 template<class T> inline
1200 T operator*(const Vector3<T>&a,const Vector3<T>&b) {
1201  return a.v[0] * b.v[0] + a.v[1] * b.v[1] + a.v[2] * b.v[2];
1202 }
1203 
1204 template<class T> inline
1205 T operator*(const Vector4<T>&a,const Vector4<T>&b) {
1206  return a.v[0]*b.v[0] + a.v[1]*b.v[1] + a.v[2]*b.v[2] + a.v[3]*b.v[3];
1207 }
1208 
1209 template<class T> inline
1210 Vector2<T> operator*(T s,const Vector2<T>&v) {
1211  return Vector2<T>(s*v[0],s*v[1]);
1212 }
1213 
1214 template<class T> inline
1215 Vector3<T> operator*(T s,const Vector3<T>&v) {
1216  return Vector3<T>(s*v[0],s*v[1],s*v[2]);
1217 }
1218 
1219 template<class T> inline
1220 Vector4<T> operator*(T s,const Vector4<T>&v) {
1221  return Vector4<T>(s*v[0],s*v[1],s*v[2],s*v[3]);
1222 }
1223 
1224 Q_DECLARE_TYPEINFO(Vec2d, Q_PRIMITIVE_TYPE);
1225 Q_DECLARE_TYPEINFO(Vec2f, Q_PRIMITIVE_TYPE);
1226 Q_DECLARE_TYPEINFO(Vec2i, Q_PRIMITIVE_TYPE);
1227 Q_DECLARE_TYPEINFO(Vec3d, Q_PRIMITIVE_TYPE);
1228 Q_DECLARE_TYPEINFO(Vec3f, Q_PRIMITIVE_TYPE);
1229 Q_DECLARE_TYPEINFO(Vec4d, Q_PRIMITIVE_TYPE);
1230 Q_DECLARE_TYPEINFO(Vec4f, Q_PRIMITIVE_TYPE);
1231 Q_DECLARE_TYPEINFO(Vec4i, Q_PRIMITIVE_TYPE);
1232 Q_DECLARE_TYPEINFO(Mat4d, Q_PRIMITIVE_TYPE);
1233 Q_DECLARE_TYPEINFO(Mat4f, Q_PRIMITIVE_TYPE);
1234 
1235 #endif // _VECMATH_H_
A templatized 2d vector compatible with openGL.
Definition: VecMath.hpp:32
A templatized column-major 4x4 matrix compatible with openGL.
Definition: VecMath.hpp:35
A templatized 4d vector compatible with openGL.
Definition: VecMath.hpp:34
QDataStream & operator>>(QDataStream &in, SphericalRegionP &region)
Load the SphericalRegionP from a binary blob.
A templatized 3d vector compatible with openGL.
Definition: VecMath.hpp:33