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