1 /*
2  * This program is free software; you can redistribute it and/or
3  * modify it under the terms of the GNU General Public License
4  * as published by the Free Software Foundation; either version 2
5  * of the License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program; if not, write to the Free Software Foundation,
14  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
15  */
16 
17 #pragma once
18 
19 /** \file
20  * \ingroup freestyle
21  * \brief Vectors and Matrices definition and manipulation
22  */
23 
24 #include <iostream>
25 #include <math.h>
26 #include <vector>
27 
28 #ifdef WITH_CXX_GUARDEDALLOC
29 #  include "MEM_guardedalloc.h"
30 #endif
31 
32 namespace Freestyle {
33 
34 namespace VecMat {
35 
36 namespace Internal {
37 template<bool B> struct is_false {
38 };
39 
40 template<> struct is_false<false> {
41   static inline void ensure()
42   {
43   }
44 };
45 }  // end of namespace Internal
46 
47 //
48 //  Vector class
49 //    - T: value type
50 //    - N: dimension
51 //
52 /////////////////////////////////////////////////////////////////////////////
53 
54 template<class T, unsigned N> class Vec {
55  public:
56   typedef T value_type;
57 
58   // constructors
59   inline Vec()
60   {
61     for (unsigned int i = 0; i < N; i++) {
62       this->_coord[i] = 0;
63     }
64   }
65 
66   ~Vec()
67   {
68     Internal::is_false<(N == 0)>::ensure();
69   }
70 
71   template<class U> explicit inline Vec(const U tab[N])
72   {
73     for (unsigned int i = 0; i < N; i++) {
74       this->_coord[i] = (T)tab[i];
75     }
76   }
77 
78   template<class U> explicit inline Vec(const std::vector<U> &tab)
79   {
80     for (unsigned int i = 0; i < N; i++) {
81       this->_coord[i] = (T)tab[i];
82     }
83   }
84 
85   template<class U> explicit inline Vec(const Vec<U, N> &v)
86   {
87     for (unsigned int i = 0; i < N; i++) {
88       this->_coord[i] = (T)v[i];
89     }
90   }
91 
92   // accessors
93   inline value_type operator[](const unsigned i) const
94   {
95     return this->_coord[i];
96   }
97 
98   inline value_type &operator[](const unsigned i)
99   {
100     return this->_coord[i];
101   }
102 
103   static inline unsigned dim()
104   {
105     return N;
106   }
107 
108   // various useful methods
109   inline value_type norm() const
110   {
111     return (T)sqrt((float)squareNorm());
112   }
113 
114   inline value_type squareNorm() const
115   {
116     return (*this) * (*this);
117   }
118 
119   inline Vec<T, N> &normalize()
120   {
121     value_type n = norm();
122     for (unsigned int i = 0; i < N; i++) {
123       this->_coord[i] /= n;
124     }
125     return *this;
126   }
127 
128   inline Vec<T, N> &normalizeSafe()
129   {
130     value_type n = norm();
131     if (n) {
132       for (unsigned int i = 0; i < N; i++) {
133         this->_coord[i] /= n;
134       }
135     }
136     return *this;
137   }
138 
139   // classical operators
140   inline Vec<T, N> operator+(const Vec<T, N> &v) const
141   {
142     Vec<T, N> res(v);
143     res += *this;
144     return res;
145   }
146 
147   inline Vec<T, N> operator-(const Vec<T, N> &v) const
148   {
149     Vec<T, N> res(*this);
150     res -= v;
151     return res;
152   }
153 
154   inline Vec<T, N> operator*(const typename Vec<T, N>::value_type r) const
155   {
156     Vec<T, N> res(*this);
157     res *= r;
158     return res;
159   }
160 
161   inline Vec<T, N> operator/(const typename Vec<T, N>::value_type r) const
162   {
163     Vec<T, N> res(*this);
164     if (r) {
165       res /= r;
166     }
167     return res;
168   }
169 
170   // dot product
171   inline value_type operator*(const Vec<T, N> &v) const
172   {
173     value_type sum = 0;
174     for (unsigned int i = 0; i < N; i++) {
175       sum += (*this)[i] * v[i];
176     }
177     return sum;
178   }
179 
180   template<class U> inline Vec<T, N> &operator=(const Vec<U, N> &v)
181   {
182     if (this != &v) {
183       for (unsigned int i = 0; i < N; i++) {
184         this->_coord[i] = (T)v[i];
185       }
186     }
187     return *this;
188   }
189 
190   template<class U> inline Vec<T, N> &operator+=(const Vec<U, N> &v)
191   {
192     for (unsigned int i = 0; i < N; i++) {
193       this->_coord[i] += (T)v[i];
194     }
195     return *this;
196   }
197 
198   template<class U> inline Vec<T, N> &operator-=(const Vec<U, N> &v)
199   {
200     for (unsigned int i = 0; i < N; i++) {
201       this->_coord[i] -= (T)v[i];
202     }
203     return *this;
204   }
205 
206   template<class U> inline Vec<T, N> &operator*=(const U r)
207   {
208     for (unsigned int i = 0; i < N; i++) {
209       this->_coord[i] *= r;
210     }
211     return *this;
212   }
213 
214   template<class U> inline Vec<T, N> &operator/=(const U r)
215   {
216     if (r) {
217       for (unsigned int i = 0; i < N; i++) {
218         this->_coord[i] /= r;
219       }
220     }
221     return *this;
222   }
223 
224   inline bool operator==(const Vec<T, N> &v) const
225   {
226     for (unsigned int i = 0; i < N; i++) {
227       if (this->_coord[i] != v[i]) {
228         return false;
229       }
230     }
231     return true;
232   }
233 
234   inline bool operator!=(const Vec<T, N> &v) const
235   {
236     for (unsigned int i = 0; i < N; i++) {
237       if (this->_coord[i] != v[i]) {
238         return true;
239       }
240     }
241     return false;
242   }
243 
244   inline bool operator<(const Vec<T, N> &v) const
245   {
246     for (unsigned int i = 0; i < N; i++) {
247       if (this->_coord[i] < v[i]) {
248         return true;
249       }
250       if (this->_coord[i] > v[i]) {
251         return false;
252       }
253       if (this->_coord[i] == v[i]) {
254         continue;
255       }
256     }
257     return false;
258   }
259 
260   inline bool operator>(const Vec<T, N> &v) const
261   {
262     for (unsigned int i = 0; i < N; i++) {
263       if (this->_coord[i] > v[i]) {
264         return true;
265       }
266       if (this->_coord[i] < v[i]) {
267         return false;
268       }
269       if (this->_coord[i] == v[i]) {
270         continue;
271       }
272     }
273     return false;
274   }
275 
276  protected:
277   value_type _coord[N];
278   enum {
279     _dim = N,
280   };
281 
282 #ifdef WITH_CXX_GUARDEDALLOC
283   MEM_CXX_CLASS_ALLOC_FUNCS("Freestyle:VecMat:Vec")
284 #endif
285 };
286 
287 //
288 //  Vec2 class (2D Vector)
289 //    - T: value type
290 //
291 /////////////////////////////////////////////////////////////////////////////
292 
293 template<class T> class Vec2 : public Vec<T, 2> {
294  public:
295   typedef typename Vec<T, 2>::value_type value_type;
296 
297   inline Vec2() : Vec<T, 2>()
298   {
299   }
300 
301   template<class U> explicit inline Vec2(const U tab[2]) : Vec<T, 2>(tab)
302   {
303   }
304 
305   template<class U> explicit inline Vec2(const std::vector<U> &tab) : Vec<T, 2>(tab)
306   {
307   }
308 
309   template<class U> inline Vec2(const Vec<U, 2> &v) : Vec<T, 2>(v)
310   {
311   }
312 
313   inline Vec2(const value_type x, const value_type y = 0) : Vec<T, 2>()
314   {
315     this->_coord[0] = (T)x;
316     this->_coord[1] = (T)y;
317   }
318 
319   inline value_type x() const
320   {
321     return this->_coord[0];
322   }
323 
324   inline value_type &x()
325   {
326     return this->_coord[0];
327   }
328 
329   inline value_type y() const
330   {
331     return this->_coord[1];
332   }
333 
334   inline value_type &y()
335   {
336     return this->_coord[1];
337   }
338 
339   inline void setX(const value_type v)
340   {
341     this->_coord[0] = v;
342   }
343 
344   inline void setY(const value_type v)
345   {
346     this->_coord[1] = v;
347   }
348 
349   // FIXME: hack swig -- no choice
350   inline Vec2<T> operator+(const Vec2<T> &v) const
351   {
352     Vec2<T> res(v);
353     res += *this;
354     return res;
355   }
356 
357   inline Vec2<T> operator-(const Vec2<T> &v) const
358   {
359     Vec2<T> res(*this);
360     res -= v;
361     return res;
362   }
363 
364   inline Vec2<T> operator*(const value_type r) const
365   {
366     Vec2<T> res(*this);
367     res *= r;
368     return res;
369   }
370 
371   inline Vec2<T> operator/(const value_type r) const
372   {
373     Vec2<T> res(*this);
374     if (r) {
375       res /= r;
376     }
377     return res;
378   }
379 
380   // dot product
381   inline value_type operator*(const Vec2<T> &v) const
382   {
383     value_type sum = 0;
384     for (unsigned int i = 0; i < 2; i++) {
385       sum += (*this)[i] * v[i];
386     }
387     return sum;
388   }
389 };
390 
391 //
392 //  HVec3 class (3D Vector in homogeneous coordinates)
393 //    - T: value type
394 //
395 /////////////////////////////////////////////////////////////////////////////
396 
397 template<class T> class HVec3 : public Vec<T, 4> {
398  public:
399   typedef typename Vec<T, 4>::value_type value_type;
400 
401   inline HVec3() : Vec<T, 4>()
402   {
403   }
404 
405   template<class U> explicit inline HVec3(const U tab[4]) : Vec<T, 4>(tab)
406   {
407   }
408 
409   template<class U> explicit inline HVec3(const std::vector<U> &tab) : Vec<T, 4>(tab)
410   {
411   }
412 
413   template<class U> inline HVec3(const Vec<U, 4> &v) : Vec<T, 4>(v)
414   {
415   }
416 
417   inline HVec3(const value_type sx,
418                const value_type sy = 0,
419                const value_type sz = 0,
420                const value_type s = 1)
421   {
422     this->_coord[0] = sx;
423     this->_coord[1] = sy;
424     this->_coord[2] = sz;
425     this->_coord[3] = s;
426   }
427 
428   template<class U> inline HVec3(const Vec<U, 3> &sv, const U s = 1)
429   {
430     this->_coord[0] = (T)sv[0];
431     this->_coord[1] = (T)sv[1];
432     this->_coord[2] = (T)sv[2];
433     this->_coord[3] = (T)s;
434   }
435 
436   inline value_type sx() const
437   {
438     return this->_coord[0];
439   }
440 
441   inline value_type &sx()
442   {
443     return this->_coord[0];
444   }
445 
446   inline value_type sy() const
447   {
448     return this->_coord[1];
449   }
450 
451   inline value_type &sy()
452   {
453     return this->_coord[1];
454   }
455 
456   inline value_type sz() const
457   {
458     return this->_coord[2];
459   }
460 
461   inline value_type &sz()
462   {
463     return this->_coord[2];
464   }
465 
466   inline value_type s() const
467   {
468     return this->_coord[3];
469   }
470 
471   inline value_type &s()
472   {
473     return this->_coord[3];
474   }
475 
476   // Access to non-homogeneous coordinates in 3D
477   inline value_type x() const
478   {
479     return this->_coord[0] / this->_coord[3];
480   }
481 
482   inline value_type y() const
483   {
484     return this->_coord[1] / this->_coord[3];
485   }
486 
487   inline value_type z() const
488   {
489     return this->_coord[2] / this->_coord[3];
490   }
491 };
492 
493 //
494 //  Vec3 class (3D Vec)
495 //    - T: value type
496 //
497 /////////////////////////////////////////////////////////////////////////////
498 template<class T> class Vec3 : public Vec<T, 3> {
499  public:
500   typedef typename Vec<T, 3>::value_type value_type;
501 
502   inline Vec3() : Vec<T, 3>()
503   {
504   }
505 
506   template<class U> explicit inline Vec3(const U tab[3]) : Vec<T, 3>(tab)
507   {
508   }
509 
510   template<class U> explicit inline Vec3(const std::vector<U> &tab) : Vec<T, 3>(tab)
511   {
512   }
513 
514   template<class U> inline Vec3(const Vec<U, 3> &v) : Vec<T, 3>(v)
515   {
516   }
517 
518   template<class U> inline Vec3(const HVec3<U> &v)
519   {
520     this->_coord[0] = (T)v.x();
521     this->_coord[1] = (T)v.y();
522     this->_coord[2] = (T)v.z();
523   }
524 
525   inline Vec3(const value_type x, const value_type y = 0, const value_type z = 0) : Vec<T, 3>()
526   {
527     this->_coord[0] = x;
528     this->_coord[1] = y;
529     this->_coord[2] = z;
530   }
531 
532   inline value_type x() const
533   {
534     return this->_coord[0];
535   }
536 
537   inline value_type &x()
538   {
539     return this->_coord[0];
540   }
541 
542   inline value_type y() const
543   {
544     return this->_coord[1];
545   }
546 
547   inline value_type &y()
548   {
549     return this->_coord[1];
550   }
551 
552   inline value_type z() const
553   {
554     return this->_coord[2];
555   }
556 
557   inline value_type &z()
558   {
559     return this->_coord[2];
560   }
561 
562   inline void setX(const value_type v)
563   {
564     this->_coord[0] = v;
565   }
566 
567   inline void setY(const value_type v)
568   {
569     this->_coord[1] = v;
570   }
571 
572   inline void setZ(const value_type v)
573   {
574     this->_coord[2] = v;
575   }
576 
577   // classical operators
578   // FIXME: hack swig -- no choice
579   inline Vec3<T> operator+(const Vec3<T> &v) const
580   {
581     Vec3<T> res(v);
582     res += *this;
583     return res;
584   }
585 
586   inline Vec3<T> operator-(const Vec3<T> &v) const
587   {
588     Vec3<T> res(*this);
589     res -= v;
590     return res;
591   }
592 
593   inline Vec3<T> operator*(const value_type r) const
594   {
595     Vec3<T> res(*this);
596     res *= r;
597     return res;
598   }
599 
600   inline Vec3<T> operator/(const value_type r) const
601   {
602     Vec3<T> res(*this);
603     if (r) {
604       res /= r;
605     }
606     return res;
607   }
608 
609   // dot product
610   inline value_type operator*(const Vec3<T> &v) const
611   {
612     value_type sum = 0;
613     for (unsigned int i = 0; i < 3; i++) {
614       sum += (*this)[i] * v[i];
615     }
616     return sum;
617   }
618 
619   // cross product for 3D Vectors
620   // FIXME: hack swig -- no choice
621   inline Vec3<T> operator^(const Vec3<T> &v) const
622   {
623     Vec3<T> res((*this)[1] * v[2] - (*this)[2] * v[1],
624                 (*this)[2] * v[0] - (*this)[0] * v[2],
625                 (*this)[0] * v[1] - (*this)[1] * v[0]);
626     return res;
627   }
628 
629   // cross product for 3D Vectors
630   template<typename U> inline Vec3<T> operator^(const Vec<U, 3> &v) const
631   {
632     Vec3<T> res((*this)[1] * v[2] - (*this)[2] * v[1],
633                 (*this)[2] * v[0] - (*this)[0] * v[2],
634                 (*this)[0] * v[1] - (*this)[1] * v[0]);
635     return res;
636   }
637 };
638 
639 //
640 //  Matrix class
641 //    - T: value type
642 //    - M: rows
643 //    - N: cols
644 //
645 /////////////////////////////////////////////////////////////////////////////
646 
647 // Dirty, but icc under Windows needs this
648 #define _SIZE (M * N)
649 
650 template<class T, unsigned M, unsigned N> class Matrix {
651  public:
652   typedef T value_type;
653 
654   inline Matrix()
655   {
656     for (unsigned int i = 0; i < _SIZE; i++) {
657       this->_coord[i] = 0;
658     }
659   }
660 
661   ~Matrix()
662   {
663     Internal::is_false<(M == 0)>::ensure();
664     Internal::is_false<(N == 0)>::ensure();
665   }
666 
667   template<class U> explicit inline Matrix(const U tab[_SIZE])
668   {
669     for (unsigned int i = 0; i < _SIZE; i++) {
670       this->_coord[i] = tab[i];
671     }
672   }
673 
674   template<class U> explicit inline Matrix(const std::vector<U> &tab)
675   {
676     for (unsigned int i = 0; i < _SIZE; i++) {
677       this->_coord[i] = tab[i];
678     }
679   }
680 
681   template<class U> inline Matrix(const Matrix<U, M, N> &m)
682   {
683     for (unsigned int i = 0; i < M; i++) {
684       for (unsigned int j = 0; j < N; j++) {
685         this->_coord[i * N + j] = (T)m(i, j);
686       }
687     }
688   }
689 
690   inline value_type operator()(const unsigned i, const unsigned j) const
691   {
692     return this->_coord[i * N + j];
693   }
694 
695   inline value_type &operator()(const unsigned i, const unsigned j)
696   {
697     return this->_coord[i * N + j];
698   }
699 
700   static inline unsigned rows()
701   {
702     return M;
703   }
704 
705   static inline unsigned cols()
706   {
707     return N;
708   }
709 
710   inline Matrix<T, M, N> &transpose() const
711   {
712     Matrix<T, N, M> res;
713     for (unsigned int i = 0; i < M; i++) {
714       for (unsigned int j = 0; j < N; j++) {
715         res(j, i) = this->_coord[i * N + j];
716       }
717     }
718     *this = res;
719     return *this;
720   }
721 
722   template<class U> inline Matrix<T, M, N> &operator=(const Matrix<U, M, N> &m)
723   {
724     if (this != &m) {
725       for (unsigned int i = 0; i < M; i++) {
726         for (unsigned int j = 0; j < N; j++) {
727           this->_coord[i * N + j] = (T)m(i, j);
728         }
729       }
730     }
731     return *this;
732   }
733 
734   template<class U> inline Matrix<T, M, N> &operator+=(const Matrix<U, M, N> &m)
735   {
736     for (unsigned int i = 0; i < M; i++) {
737       for (unsigned int j = 0; j < N; j++) {
738         this->_coord[i * N + j] += (T)m(i, j);
739       }
740     }
741     return *this;
742   }
743 
744   template<class U> inline Matrix<T, M, N> &operator-=(const Matrix<U, M, N> &m)
745   {
746     for (unsigned int i = 0; i < M; i++) {
747       for (unsigned int j = 0; j < N; j++) {
748         this->_coord[i * N + j] -= (T)m(i, j);
749       }
750     }
751     return *this;
752   }
753 
754   template<class U> inline Matrix<T, M, N> &operator*=(const U lambda)
755   {
756     for (unsigned int i = 0; i < M; i++) {
757       for (unsigned int j = 0; j < N; j++) {
758         this->_coord[i * N + j] *= lambda;
759       }
760     }
761     return *this;
762   }
763 
764   template<class U> inline Matrix<T, M, N> &operator/=(const U lambda)
765   {
766     if (lambda) {
767       for (unsigned int i = 0; i < M; i++) {
768         for (unsigned int j = 0; j < N; j++) {
769           this->_coord[i * N + j] /= lambda;
770         }
771       }
772     }
773     return *this;
774   }
775 
776  protected:
777   value_type _coord[_SIZE];
778 
779 #ifdef WITH_CXX_GUARDEDALLOC
780   MEM_CXX_CLASS_ALLOC_FUNCS("Freestyle:VecMat:Matrix")
781 #endif
782 };
783 
784 #undef _SIZE
785 
786 //
787 //  SquareMatrix class
788 //    - T: value type
789 //    - N: rows & cols
790 //
791 /////////////////////////////////////////////////////////////////////////////
792 
793 // Dirty, but icc under Windows needs this
794 #define _SIZE (N * N)
795 
796 template<class T, unsigned N> class SquareMatrix : public Matrix<T, N, N> {
797  public:
798   typedef T value_type;
799 
800   inline SquareMatrix() : Matrix<T, N, N>()
801   {
802   }
803 
804   template<class U> explicit inline SquareMatrix(const U tab[_SIZE]) : Matrix<T, N, N>(tab)
805   {
806   }
807 
808   template<class U> explicit inline SquareMatrix(const std::vector<U> &tab) : Matrix<T, N, N>(tab)
809   {
810   }
811 
812   template<class U> inline SquareMatrix(const Matrix<U, N, N> &m) : Matrix<T, N, N>(m)
813   {
814   }
815 
816   static inline SquareMatrix<T, N> identity()
817   {
818     SquareMatrix<T, N> res;
819     for (unsigned int i = 0; i < N; i++) {
820       res(i, i) = 1;
821     }
822     return res;
823   }
824 };
825 
826 #undef _SIZE
827 
828 //
829 // Vector external functions
830 //
831 /////////////////////////////////////////////////////////////////////////////
832 
833 #if 0
834 template<class T, unsigned N> inline Vec<T, N> operator+(const Vec<T, N> &v1, const Vec<T, N> &v2)
835 {
836   Vec<T, N> res(v1);
837   res += v2;
838   return res;
839 }
840 
841 template<class T, unsigned N> inline Vec<T, N> operator-(const Vec<T, N> &v1, const Vec<T, N> &v2)
842 {
843   Vec<T, N> res(v1);
844   res -= v2;
845   return res;
846 }
847 
848 template<class T, unsigned N>
849 inline Vec<T, N> operator*(const Vec<T, N> &v, const typename Vec<T, N>::value_type r)
850 {
851   Vec<T, N> res(v);
852   res *= r;
853   return res;
854 }
855 #endif
856 
857 template<class T, unsigned N>
858 inline Vec<T, N> operator*(const typename Vec<T, N>::value_type r, const Vec<T, N> &v)
859 {
860   Vec<T, N> res(v);
861   res *= r;
862   return res;
863 }
864 
865 #if 0
866 template<class T, unsigned N>
867 inline Vec<T, N> operator/(const Vec<T, N> &v, const typename Vec<T, N>::value_type r)
868 {
869   Vec<T, N> res(v);
870   if (r) {
871     res /= r;
872   }
873   return res;
874 }
875 
876 // dot product
877 template<class T, unsigned N>
878 inline typename Vec<T, N>::value_type operator*(const Vec<T, N> &v1, const Vec<T, N> &v2)
879 {
880   typename Vec<T, N>::value_type sum = 0;
881   for (unsigned int i = 0; i < N; i++) {
882     sum += v1[i] * v2[i];
883   }
884   return sum;
885 }
886 
887 // cross product for 3D Vectors
888 template<typename T> inline Vec3<T> operator^(const Vec<T, 3> &v1, const Vec<T, 3> &v2)
889 {
890   Vec3<T> res(
891       v1[1] * v2[2] - v1[2] * v2[1], v1[2] * v2[0] - v1[0] * v2[2], v1[0] * v2[1] - v1[1] * v2[0]);
892   return res;
893 }
894 #endif
895 
896 // stream operator
897 template<class T, unsigned N> inline std::ostream &operator<<(std::ostream &s, const Vec<T, N> &v)
898 {
899   unsigned int i;
900   s << "[";
901   for (i = 0; i < N - 1; i++) {
902     s << v[i] << ", ";
903   }
904   s << v[i] << "]";
905   return s;
906 }
907 
908 //
909 // Matrix external functions
910 //
911 /////////////////////////////////////////////////////////////////////////////
912 
913 template<class T, unsigned M, unsigned N>
914 inline Matrix<T, M, N> operator+(const Matrix<T, M, N> &m1, const Matrix<T, M, N> &m2)
915 {
916   Matrix<T, M, N> res(m1);
917   res += m2;
918   return res;
919 }
920 
921 template<class T, unsigned M, unsigned N>
922 inline Matrix<T, M, N> operator-(const Matrix<T, M, N> &m1, const Matrix<T, M, N> &m2)
923 {
924   Matrix<T, M, N> res(m1);
925   res -= m2;
926   return res;
927 }
928 
929 template<class T, unsigned M, unsigned N>
930 inline Matrix<T, M, N> operator*(const Matrix<T, M, N> &m1,
931                                  const typename Matrix<T, M, N>::value_type lambda)
932 {
933   Matrix<T, M, N> res(m1);
934   res *= lambda;
935   return res;
936 }
937 
938 template<class T, unsigned M, unsigned N>
939 inline Matrix<T, M, N> operator*(const typename Matrix<T, M, N>::value_type lambda,
940                                  const Matrix<T, M, N> &m1)
941 {
942   Matrix<T, M, N> res(m1);
943   res *= lambda;
944   return res;
945 }
946 
947 template<class T, unsigned M, unsigned N>
948 inline Matrix<T, M, N> operator/(const Matrix<T, M, N> &m1,
949                                  const typename Matrix<T, M, N>::value_type lambda)
950 {
951   Matrix<T, M, N> res(m1);
952   res /= lambda;
953   return res;
954 }
955 
956 template<class T, unsigned M, unsigned N, unsigned P>
957 inline Matrix<T, M, P> operator*(const Matrix<T, M, N> &m1, const Matrix<T, N, P> &m2)
958 {
959   unsigned int i, j, k;
960   Matrix<T, M, P> res;
961   typename Matrix<T, N, P>::value_type scale;
962 
963   for (j = 0; j < P; j++) {
964     for (k = 0; k < N; k++) {
965       scale = m2(k, j);
966       for (i = 0; i < N; i++) {
967         res(i, j) += m1(i, k) * scale;
968       }
969     }
970   }
971   return res;
972 }
973 
974 template<class T, unsigned M, unsigned N>
975 inline Vec<T, M> operator*(const Matrix<T, M, N> &m, const Vec<T, N> &v)
976 {
977   Vec<T, M> res;
978   typename Matrix<T, M, N>::value_type scale;
979 
980   for (unsigned int j = 0; j < M; j++) {
981     scale = v[j];
982     for (unsigned int i = 0; i < N; i++) {
983       res[i] += m(i, j) * scale;
984     }
985   }
986   return res;
987 }
988 
989 // stream operator
990 template<class T, unsigned M, unsigned N>
991 inline std::ostream &operator<<(std::ostream &s, const Matrix<T, M, N> &m)
992 {
993   unsigned int i, j;
994   for (i = 0; i < M; i++) {
995     s << "[";
996     for (j = 0; j < N - 1; j++) {
997       s << m(i, j) << ", ";
998     }
999     s << m(i, j) << "]" << std::endl;
1000   }
1001   return s;
1002 }
1003 
1004 }  // end of namespace VecMat
1005 
1006 } /* namespace Freestyle */
1007