WebSVN – gelsvn – Diff – /branches/imesh-and-hmesh/src/CGLA/ArithMatFloat.h


#ifndef __CGLA_ARITHMATFLOAT_H__
#define __CGLA_ARITHMATFLOAT_H__
 
#include <vector>
#include <iostream>
 
#include "CGLA.h"
 
 
namespace CGLA 
{
 
  /** Basic class template for matrices.
		
  In this template a matrix is defined as an array of vectors. This may
  not in all cases be the most efficient but it has the advantage that 
  it is possible to use the double subscripting notation:
		
  T x = m[i][j]
 
  This template should be used through inheritance just like the 
  vector template */
  template <class VVT, class HVT, class MT, int ROWS>
    class ArithMatFloat
    { 
#define for_all_i(expr) for(int i=0;i<ROWS;i++) {expr;}
 
    public:
 
      /// Horizontal vector type
      typedef HVT HVectorType;
 
      /// Vertical vector type
      typedef VVT VVectorType;
 
      /// The type of a matrix element
      typedef typename HVT::ScalarType ScalarType;
		
    protected:
 
      /// The actual contents of the matrix.
      HVT data[ROWS];
 
    protected:
 
      /// Construct 0 matrix
      ArithMatFloat() 
	{
#ifndef NDEBUG
	  for_all_i(data[i]=HVT(CGLA_INIT_VALUE));
#endif
	}
 
      /// Construct a matrix where all entries are the same.
      explicit ArithMatFloat(ScalarType x)
	{
	  for_all_i(data[i] = HVT(x));
	}
 
      /// Construct a matrix where all rows are the same.
      explicit ArithMatFloat(HVT _a)
	{
	  for_all_i(data[i] = _a);
	}
 
 
      /// Construct a matrix with two rows.
      ArithMatFloat(HVT _a, HVT _b)
	{
	  assert(ROWS==2);
	  data[0] = _a;
	  data[1] = _b;
	}
 
      /// Construct a matrix with three rows.
      ArithMatFloat(HVT _a, HVT _b, HVT _c)
	{
	  assert(ROWS==3);
	  data[0] = _a;
	  data[1] = _b;
	  data[2] = _c;
	}
 
      /// Construct a matrix with four rows.
      ArithMatFloat(HVT _a, HVT _b, HVT _c, HVT _d)
	{
	  assert(ROWS==4);
	  data[0] = _a;
	  data[1] = _b;
	  data[2] = _c;
	  data[3] = _d;
	}
		
    public:
 
      /// Get vertical dimension of matrix 
      static int get_v_dim() {return VVT::get_dim();}
 
      /// Get horizontal dimension of matrix
      static int get_h_dim() {return HVT::get_dim();}
 
 
      /** Get const pointer to data array.
	  This function may be useful when interfacing with some other API 
	  such as OpenGL (TM). */
      const ScalarType* get() const 
	{
	  return data[0].get();
	}
 
      /** Get pointer to data array.
	  This function may be useful when interfacing with some other API 
	  such as OpenGL (TM). */
      ScalarType* get()
	{
	  return data[0].get();
	}
 
      //----------------------------------------------------------------------
      // index operators
      //----------------------------------------------------------------------
 
      /// Const index operator. Returns i'th row of matrix.
      const HVT& operator [] ( int i ) const
	{
	  assert(i<ROWS);
	  return data[i];
	}
 
      /// Non-const index operator. Returns i'th row of matrix.
      HVT& operator [] ( int i ) 
	{
	  assert(i<ROWS);
	  return data[i];
	}
 
      //----------------------------------------------------------------------
 
      /// Equality operator. 
      bool operator==(const MT& v) const 
	{
	  for_all_i(if (data[i] != v[i]) return false)
	    return true;
	}
 
      /// Inequality operator.
      bool operator!=(const MT& v) const 
	{
	  return !(*this==v);
	}
 
      //----------------------------------------------------------------------
 
      /// Multiply scalar onto matrix. All entries are multiplied by scalar.
      const MT operator * (ScalarType k) const
	{
	  MT v_new;
	  for_all_i(v_new[i] = data[i] * k);
	  return v_new;
	}
 
      /// Divide all entries in matrix by scalar.
      const MT operator / (ScalarType k) const
	{
	  MT v_new;
	  for_all_i(v_new[i] = data[i] / k);
	  return v_new;      
	}
 
      /// Assignment multiplication of matrix by scalar.
      const MT& operator *=(ScalarType k) 
	{
	  for_all_i(data[i] *= k); 
	  return static_cast<const MT&>(*this);
	}
 
      /// Assignment division of matrix by scalar.
      const MT& operator /=(ScalarType k) 
	{ 
	  for_all_i(data[i] /= k); 
	  return static_cast<const MT&>(*this);
	}
 
      //----------------------------------------------------------------------
 
      /// Add two matrices. 
      const MT operator + (const MT& m1) const
	{
	  MT v_new;
	  for_all_i(v_new[i] = data[i] + m1[i]);
	  return v_new;
	}
 
      /// Subtract two matrices.
      const MT operator - (const MT& m1) const
	{
	  MT v_new;
	  for_all_i(v_new[i] = data[i] - m1[i]);
	  return v_new;
	}
 
      /// Assigment addition of matrices.
      const MT& operator +=(const MT& v) 
	{
	  for_all_i(data[i] += v[i]); 
	  return static_cast<const MT&>(*this);
	}
 
      /// Assigment subtraction of matrices.
      const MT& operator -=(const MT& v) 
	{
	  for_all_i(data[i] -= v[i]); 
	  return static_cast<const MT&>(*this);
	}
 
      //----------------------------------------------------------------------
 
      /// Negate matrix.
      const MT operator - () const
	{
	  MT v_new;
	  for_all_i(v_new[i] = - data[i]);
	  return v_new;
	}
 
#undef for_all_i  
 
    };
 
  /// Multiply scalar onto matrix
  template <class VVT, class HVT, class MT, int ROWS>
    inline const MT operator * (double k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v) 
    {
      return v * k;
    }
 
  /// Multiply scalar onto matrix
  template <class VVT, class HVT, class MT, int ROWS>
    inline const MT operator * (float k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v) 
    {
      return v * k;
    }
 
  /// Multiply scalar onto matrix
  template <class VVT, class HVT, class MT, int ROWS>
    inline const MT operator * (int k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v) 
    {
      return v * k;
    }
 
  /// Multiply vector onto matrix 
  template <class VVT, class HVT, class MT, int ROWS>
    inline VVT operator*(const ArithMatFloat<VVT,HVT,MT,ROWS>& m,const HVT& v) 
    {
      VVT v2;
      for(int i=0;i<ROWS;i++) v2[i] = dot(m[i], v);
      return v2;
    }
 
 
#ifndef WIN32
  /** Multiply two arbitrary matrices. 
      In principle, this function could return a matrix, but in general
      the new matrix will be of a type that is different from either of
      the two matrices that are multiplied together. We do not want to 
      return an ArithMatFloat - so it seems best to let the return value be
      a reference arg.
		
      This template can only be instantiated if the dimensions of the
      matrices match -- i.e. if the multiplication can actually be
      carried out. This is more type safe than the win32 version below.
  */
 
  template <class VVT, class HVT, 
    class HV1T, class VV2T,
    class MT1, class MT2, class MT,
    int ROWS1, int ROWS2>
    inline void mul(const ArithMatFloat<VVT,HV1T,MT1,ROWS1>& m1,
		    const ArithMatFloat<VV2T,HVT,MT2,ROWS2>& m2,
		    ArithMatFloat<VVT,HVT,MT,ROWS1>& m)
    {
      int cols = ArithMatFloat<VVT,HVT,MT,ROWS1>::get_h_dim();
      for(int i=0;i<ROWS1;i++)
	for(int j=0;j<cols;j++)
	  {
	    m[i][j] = 0;
	    for(int k=0;k<ROWS2;k++)
	      m[i][j] += m1[i][k] * m2[k][j]; 
	  }
    }
 
 
  /** Transpose. See the discussion on mul if you are curious as to why
      I don't simply return the transpose. */
  template <class VVT, class HVT, class M1T, class M2T, int ROWS, int COLS>
    inline void transpose(const ArithMatFloat<VVT,HVT,M1T,ROWS>& m,
			  ArithMatFloat<HVT,VVT,M2T,COLS>& m_new)
    {
      for(int i=0;i<M2T::get_v_dim();i++)
	for(int j=0;j<M2T::get_h_dim();j++)
	  m_new[i][j] = m[j][i];
    }
 
#else
 
  //----------------- win32 -------------------------------
  // Visual studio is not good at deducing the args. to these template functions.
  // This means that you can call the two functions below with 
  // matrices of wrong dimension.
 
  template <class M1, class M2, class M>
    inline void mul(const M1& m1, const M2& m2, M& m)
    {
      int cols = M::get_h_dim();
      int rows1 = M1::get_v_dim();
      int rows2 = M2::get_v_dim();
 
      for(int i=0;i<rows1;i++)
	for(int j=0;j<cols;j++)
	  {
	    m[i][j] = 0;
	    for(int k=0;k<rows2;k++)
	      m[i][j] += m1[i][k] * m2[k][j];
	  }
    }
 
 
  /** Transpose. See the discussion on mul if you are curious as to why
      I don't simply return the transpose. */
  template <class M1, class M2>
    inline void transpose(const M1& m1, M2& m2)
    {
      for(int i=0;i<M2::get_v_dim();i++)
	for(int j=0;j<M2::get_h_dim();j++)
	  m2[i][j] = m1[j][i];
    }
 
#endif
 
  /** Compute the outer product of a and b: a * transpose(b). This is 
      a matrix with a::rows and b::columns. */
  template <class VVT, class HVT, class MT, int ROWS>
    void outer_product(const VVT& a, const HVT& b, 
		       ArithMatFloat<VVT,HVT,MT,ROWS>& m)
    {
      int R = VVT::get_dim();
      int C = HVT::get_dim();
      for(int i=0;i<R;++i)
	for(int j=0;j<C;++j)
	  {
	    m[i][j] = a[i] * b[j];
	  }
    }
 
  /** Copy a matrix to another matrix, cell by cell.
      This conversion that takes a const matrix as first argument
      (source) and a non-const matrix as second argument
      (destination). The contents of the first matrix is simply copied
      to the second matrix. 
			
      However, if the first matrix is	larger than the second,
      the cells outside the range of the destination are simply not
      copied. If the destination is larger, the cells outside the 
      range of the source matrix are not touched.
 
      An obvious use of this function is to copy a 3x3 rotation matrix
      into a 4x4 transformation matrix.
  */
 
  template <class M1, class M2>
    void copy_matrix(const M1& inmat, M2& outmat)
    {
      const int R = s_min(inmat.get_v_dim(), outmat.get_v_dim());
      const int C = s_min(inmat.get_h_dim(), outmat.get_h_dim());
      for(int i=0;i<R;++i)
	for(int j=0;j<C;++j)
	  outmat[i][j] = inmat[i][j];
    }
 
  /** Put to operator */
  template <class VVT, class HVT, class MT, int ROWS>
    inline std::ostream& 
    operator<<(std::ostream&os, const ArithMatFloat<VVT,HVT,MT,ROWS>& m)
    {
      os << "[\n";
      for(int i=0;i<ROWS;i++) os << "  " << m[i] << "\n";
      os << "]\n";
      return os;
    }
 
  /** Get from operator */
  template <class VVT, class HVT, class MT, int ROWS>
    inline std::istream& operator>>(std::istream&is, 
				    const ArithMatFloat<VVT,HVT,MT,ROWS>& m)
    {
      for(int i=0;i<ROWS;i++) is>>m[i];
      return is;
    }
}
#endif
 

Rev 12	Rev 22
1	`#ifndef __CGLA_ARITHMATFLOAT_H__`	1	`#ifndef __CGLA_ARITHMATFLOAT_H__`
2	`#define __CGLA_ARITHMATFLOAT_H__`	2	`#define __CGLA_ARITHMATFLOAT_H__`
3		3
4	`#include <vector>`	4	`#include <vector>`
5	`#include <iostream>`	5	`#include <iostream>`
6		6
7	`#include "CGLA.h"`	7	`#include "CGLA.h"`
8		8
9		9
10	`namespace CGLA`	10	`namespace CGLA`
11	`{`	11	`{`
12		12
13	`/** Basic class template for matrices.`	13	`/** Basic class template for matrices.`
14		14
15	`In this template a matrix is defined as an array of vectors. This may`	15	`In this template a matrix is defined as an array of vectors. This may`
16	`not in all cases be the most efficient but it has the advantage that`	16	`not in all cases be the most efficient but it has the advantage that`
17	`it is possible to use the double subscripting notation:`	17	`it is possible to use the double subscripting notation:`
18		18
19	`T x = m[i][j]`	19	`T x = m[i][j]`
20		20
21	`This template should be used through inheritance just like the`	21	`This template should be used through inheritance just like the`
22	`vector template */`	22	`vector template */`
23	`template <class VVT, class HVT, class MT, int ROWS>`	23	`template <class VVT, class HVT, class MT, int ROWS>`
24	`class ArithMatFloat`	24	`class ArithMatFloat`
25	`{`	25	`{`
26	`#define for_all_i(expr) for(int i=0;i<ROWS;i++) {expr;}`	26	`#define for_all_i(expr) for(int i=0;i<ROWS;i++) {expr;}`
27		27
28	`public:`	28	`public:`
29		29
30	`/// Horizontal vector type`	30	`/// Horizontal vector type`
31	`typedef HVT HVectorType;`	31	`typedef HVT HVectorType;`
32		32
33	`/// Vertical vector type`	33	`/// Vertical vector type`
34	`typedef VVT VVectorType;`	34	`typedef VVT VVectorType;`
35		35
36	`/// The type of a matrix element`	36	`/// The type of a matrix element`
37	`typedef typename HVT::ScalarType ScalarType;`	37	`typedef typename HVT::ScalarType ScalarType;`
38		38
39	`protected:`	39	`protected:`
40		40
41	`/// The actual contents of the matrix.`	41	`/// The actual contents of the matrix.`
42	`HVT data[ROWS];`	42	`HVT data[ROWS];`
43		43
44	`protected:`	44	`protected:`
45		45
46	`/// Construct 0 matrix`	46	`/// Construct 0 matrix`
47	`ArithMatFloat()`	47	`ArithMatFloat()`
48	`{`	48	`{`
49	`#ifndef NDEBUG`	49	`#ifndef NDEBUG`
50	`for_all_i(data[i]=HVT(CGLA_INIT_VALUE));`	50	`for_all_i(data[i]=HVT(CGLA_INIT_VALUE));`
51	`#endif`	51	`#endif`
52	`}`	52	`}`
53		53
54	`/// Construct a matrix where all entries are the same.`	54	`/// Construct a matrix where all entries are the same.`
55	`explicit ArithMatFloat(ScalarType x)`	55	`explicit ArithMatFloat(ScalarType x)`
56	`{`	56	`{`
57	`for_all_i(data[i] = HVT(x));`	57	`for_all_i(data[i] = HVT(x));`
58	`}`	58	`}`
59		59
60	`/// Construct a matrix where all rows are the same.`	60	`/// Construct a matrix where all rows are the same.`
61	`explicit ArithMatFloat(HVT _a)`	61	`explicit ArithMatFloat(HVT _a)`
62	`{`	62	`{`
63	`for_all_i(data[i] = _a);`	63	`for_all_i(data[i] = _a);`
64	`}`	64	`}`
65		65
66		66
67	`/// Construct a matrix with two rows.`	67	`/// Construct a matrix with two rows.`
68	`ArithMatFloat(HVT _a, HVT _b)`	68	`ArithMatFloat(HVT _a, HVT _b)`
69	`{`	69	`{`
70	`assert(ROWS==2);`	70	`assert(ROWS==2);`
71	`data[0] = _a;`	71	`data[0] = _a;`
72	`data[1] = _b;`	72	`data[1] = _b;`
73	`}`	73	`}`
74		74
75	`/// Construct a matrix with three rows.`	75	`/// Construct a matrix with three rows.`
76	`ArithMatFloat(HVT _a, HVT _b, HVT _c)`	76	`ArithMatFloat(HVT _a, HVT _b, HVT _c)`
77	`{`	77	`{`
78	`assert(ROWS==3);`	78	`assert(ROWS==3);`
79	`data[0] = _a;`	79	`data[0] = _a;`
80	`data[1] = _b;`	80	`data[1] = _b;`
81	`data[2] = _c;`	81	`data[2] = _c;`
82	`}`	82	`}`
83		83
84	`/// Construct a matrix with four rows.`	84	`/// Construct a matrix with four rows.`
85	`ArithMatFloat(HVT _a, HVT _b, HVT _c, HVT _d)`	85	`ArithMatFloat(HVT _a, HVT _b, HVT _c, HVT _d)`
86	`{`	86	`{`
87	`assert(ROWS==4);`	87	`assert(ROWS==4);`
88	`data[0] = _a;`	88	`data[0] = _a;`
89	`data[1] = _b;`	89	`data[1] = _b;`
90	`data[2] = _c;`	90	`data[2] = _c;`
91	`data[3] = _d;`	91	`data[3] = _d;`
92	`}`	92	`}`
93		93
94	`public:`	94	`public:`
95		95
96	`/// Get vertical dimension of matrix`	96	`/// Get vertical dimension of matrix`
97	`static int get_v_dim() {return VVT::get_dim();}`	97	`static int get_v_dim() {return VVT::get_dim();}`
98		98
99	`/// Get horizontal dimension of matrix`	99	`/// Get horizontal dimension of matrix`
100	`static int get_h_dim() {return HVT::get_dim();}`	100	`static int get_h_dim() {return HVT::get_dim();}`
101		101
102		102
103	`/** Get const pointer to data array.`	103	`/** Get const pointer to data array.`
104	`This function may be useful when interfacing with some other API`	104	`This function may be useful when interfacing with some other API`
105	`such as OpenGL (TM). */`	105	`such as OpenGL (TM). */`
106	`const ScalarType* get() const`	106	`const ScalarType* get() const`
107	`{`	107	`{`
108	`return data[0].get();`	108	`return data[0].get();`
109	`}`	109	`}`
110		110
111	`/** Get pointer to data array.`	111	`/** Get pointer to data array.`
112	`This function may be useful when interfacing with some other API`	112	`This function may be useful when interfacing with some other API`
113	`such as OpenGL (TM). */`	113	`such as OpenGL (TM). */`
114	`ScalarType* get()`	114	`ScalarType* get()`
115	`{`	115	`{`
116	`return data[0].get();`	116	`return data[0].get();`
117	`}`	117	`}`
118		118
119	`//----------------------------------------------------------------------`	119	`//----------------------------------------------------------------------`
120	`// index operators`	120	`// index operators`
121	`//----------------------------------------------------------------------`	121	`//----------------------------------------------------------------------`
122		122
123	`/// Const index operator. Returns i'th row of matrix.`	123	`/// Const index operator. Returns i'th row of matrix.`
124	`const HVT& operator [] ( int i ) const`	124	`const HVT& operator [] ( int i ) const`
125	`{`	125	`{`
126	`assert(i<ROWS);`	126	`assert(i<ROWS);`
127	`return data[i];`	127	`return data[i];`
128	`}`	128	`}`
129		129
130	`/// Non-const index operator. Returns i'th row of matrix.`	130	`/// Non-const index operator. Returns i'th row of matrix.`
131	`HVT& operator [] ( int i )`	131	`HVT& operator [] ( int i )`
132	`{`	132	`{`
133	`assert(i<ROWS);`	133	`assert(i<ROWS);`
134	`return data[i];`	134	`return data[i];`
135	`}`	135	`}`
136		136
137	`//----------------------------------------------------------------------`	137	`//----------------------------------------------------------------------`
138		138
139	`/// Equality operator.`	139	`/// Equality operator.`
140	`bool operator==(const MT& v) const`	140	`bool operator==(const MT& v) const`
141	`{`	141	`{`
142	`for_all_i(if (data[i] != v[i]) return false)`	142	`for_all_i(if (data[i] != v[i]) return false)`
143	`return true;`	143	`return true;`
144	`}`	144	`}`
145		145
146	`/// Inequality operator.`	146	`/// Inequality operator.`
147	`bool operator!=(const MT& v) const`	147	`bool operator!=(const MT& v) const`
148	`{`	148	`{`
149	`return !(*this==v);`	149	`return !(*this==v);`
150	`}`	150	`}`
151		151
152	`//----------------------------------------------------------------------`	152	`//----------------------------------------------------------------------`
153		153
154	`/// Multiply scalar onto matrix. All entries are multiplied by scalar.`	154	`/// Multiply scalar onto matrix. All entries are multiplied by scalar.`
155	`const MT operator * (ScalarType k) const`	155	`const MT operator * (ScalarType k) const`
156	`{`	156	`{`
157	`MT v_new;`	157	`MT v_new;`
158	`for_all_i(v_new[i] = data[i] * k);`	158	`for_all_i(v_new[i] = data[i] * k);`
159	`return v_new;`	159	`return v_new;`
160	`}`	160	`}`
161		161
162	`/// Divide all entries in matrix by scalar.`	162	`/// Divide all entries in matrix by scalar.`
163	`const MT operator / (ScalarType k) const`	163	`const MT operator / (ScalarType k) const`
164	`{`	164	`{`
165	`MT v_new;`	165	`MT v_new;`
166	`for_all_i(v_new[i] = data[i] / k);`	166	`for_all_i(v_new[i] = data[i] / k);`
167	`return v_new;`	167	`return v_new;`
168	`}`	168	`}`
169		169
170	`/// Assignment multiplication of matrix by scalar.`	170	`/// Assignment multiplication of matrix by scalar.`
171	`const MT& operator *=(ScalarType k)`	171	`const MT& operator *=(ScalarType k)`
172	`{`	172	`{`
173	`for_all_i(data[i] *= k);`	173	`for_all_i(data[i] *= k);`
174	`return static_cast<const MT&>(*this);`	174	`return static_cast<const MT&>(*this);`
175	`}`	175	`}`
176		176
177	`/// Assignment division of matrix by scalar.`	177	`/// Assignment division of matrix by scalar.`
178	`const MT& operator /=(ScalarType k)`	178	`const MT& operator /=(ScalarType k)`
179	`{`	179	`{`
180	`for_all_i(data[i] /= k);`	180	`for_all_i(data[i] /= k);`
181	`return static_cast<const MT&>(*this);`	181	`return static_cast<const MT&>(*this);`
182	`}`	182	`}`
183		183
184	`//----------------------------------------------------------------------`	184	`//----------------------------------------------------------------------`
185		185
186	`/// Add two matrices.`	186	`/// Add two matrices.`
187	`const MT operator + (const MT& m1) const`	187	`const MT operator + (const MT& m1) const`
188	`{`	188	`{`
189	`MT v_new;`	189	`MT v_new;`
190	`for_all_i(v_new[i] = data[i] + m1[i]);`	190	`for_all_i(v_new[i] = data[i] + m1[i]);`
191	`return v_new;`	191	`return v_new;`
192	`}`	192	`}`
193		193
194	`/// Subtract two matrices.`	194	`/// Subtract two matrices.`
195	`const MT operator - (const MT& m1) const`	195	`const MT operator - (const MT& m1) const`
196	`{`	196	`{`
197	`MT v_new;`	197	`MT v_new;`
198	`for_all_i(v_new[i] = data[i] - m1[i]);`	198	`for_all_i(v_new[i] = data[i] - m1[i]);`
199	`return v_new;`	199	`return v_new;`
200	`}`	200	`}`
201		201
202	`/// Assigment addition of matrices.`	202	`/// Assigment addition of matrices.`
203	`const MT& operator +=(const MT& v)`	203	`const MT& operator +=(const MT& v)`
204	`{`	204	`{`
205	`for_all_i(data[i] += v[i]);`	205	`for_all_i(data[i] += v[i]);`
206	`return static_cast<const MT&>(*this);`	206	`return static_cast<const MT&>(*this);`
207	`}`	207	`}`
208		208
209	`/// Assigment subtraction of matrices.`	209	`/// Assigment subtraction of matrices.`
210	`const MT& operator -=(const MT& v)`	210	`const MT& operator -=(const MT& v)`
211	`{`	211	`{`
212	`for_all_i(data[i] -= v[i]);`	212	`for_all_i(data[i] -= v[i]);`
213	`return static_cast<const MT&>(*this);`	213	`return static_cast<const MT&>(*this);`
214	`}`	214	`}`
215		215
216	`//----------------------------------------------------------------------`	216	`//----------------------------------------------------------------------`
217		217
218	`/// Negate matrix.`	218	`/// Negate matrix.`
219	`const MT operator - () const`	219	`const MT operator - () const`
220	`{`	220	`{`
221	`MT v_new;`	221	`MT v_new;`
222	`for_all_i(v_new[i] = - data[i]);`	222	`for_all_i(v_new[i] = - data[i]);`
223	`return v_new;`	223	`return v_new;`
224	`}`	224	`}`
225		225
226	`#undef for_all_i`	226	`#undef for_all_i`
227		227
228	`};`	228	`};`
229		229
230	`/// Multiply scalar onto matrix`	230	`/// Multiply scalar onto matrix`
231	`template <class VVT, class HVT, class MT, int ROWS>`	231	`template <class VVT, class HVT, class MT, int ROWS>`
232	`inline const MT operator * (double k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v)`	232	`inline const MT operator * (double k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v)`
233	`{`	233	`{`
234	`return v * k;`	234	`return v * k;`
235	`}`	235	`}`
236		236
237	`/// Multiply scalar onto matrix`	237	`/// Multiply scalar onto matrix`
238	`template <class VVT, class HVT, class MT, int ROWS>`	238	`template <class VVT, class HVT, class MT, int ROWS>`
239	`inline const MT operator * (float k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v)`	239	`inline const MT operator * (float k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v)`
240	`{`	240	`{`
241	`return v * k;`	241	`return v * k;`
242	`}`	242	`}`
243		243
244	`/// Multiply scalar onto matrix`	244	`/// Multiply scalar onto matrix`
245	`template <class VVT, class HVT, class MT, int ROWS>`	245	`template <class VVT, class HVT, class MT, int ROWS>`
246	`inline const MT operator * (int k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v)`	246	`inline const MT operator * (int k, const ArithMatFloat<VVT,HVT,MT,ROWS>& v)`
247	`{`	247	`{`
248	`return v * k;`	248	`return v * k;`
249	`}`	249	`}`
250		250
251	`/// Multiply vector onto matrix`	251	`/// Multiply vector onto matrix`
252	`template <class VVT, class HVT, class MT, int ROWS>`	252	`template <class VVT, class HVT, class MT, int ROWS>`
253	`inline VVT operator*(const ArithMatFloat<VVT,HVT,MT,ROWS>& m,const HVT& v)`	253	`inline VVT operator*(const ArithMatFloat<VVT,HVT,MT,ROWS>& m,const HVT& v)`
254	`{`	254	`{`
255	`VVT v2;`	255	`VVT v2;`
256	`for(int i=0;i<ROWS;i++) v2[i] = dot(m[i], v);`	256	`for(int i=0;i<ROWS;i++) v2[i] = dot(m[i], v);`
257	`return v2;`	257	`return v2;`
258	`}`	258	`}`
259		259
260		260
261	`#ifndef WIN32`	261	`#ifndef WIN32`
262	`/** Multiply two arbitrary matrices.`	262	`/** Multiply two arbitrary matrices.`
263	`In principle, this function could return a matrix, but in general`	263	`In principle, this function could return a matrix, but in general`
264	`the new matrix will be of a type that is different from either of`	264	`the new matrix will be of a type that is different from either of`
265	`the two matrices that are multiplied together. We do not want to`	265	`the two matrices that are multiplied together. We do not want to`
266	`return an ArithMatFloat - so it seems best to let the return value be`	266	`return an ArithMatFloat - so it seems best to let the return value be`
267	`a reference arg.`	267	`a reference arg.`
268		268
269	`This template can only be instantiated if the dimensions of the`	269	`This template can only be instantiated if the dimensions of the`
270	`matrices match -- i.e. if the multiplication can actually be`	270	`matrices match -- i.e. if the multiplication can actually be`
271	`carried out. This is more type safe than the win32 version below.`	271	`carried out. This is more type safe than the win32 version below.`
272	`*/`	272	`*/`
273		273
274	`template <class VVT, class HVT,`	274	`template <class VVT, class HVT,`
275	`class HV1T, class VV2T,`	275	`class HV1T, class VV2T,`
276	`class MT1, class MT2, class MT,`	276	`class MT1, class MT2, class MT,`
277	`int ROWS1, int ROWS2>`	277	`int ROWS1, int ROWS2>`
278	`inline void mul(const ArithMatFloat<VVT,HV1T,MT1,ROWS1>& m1,`	278	`inline void mul(const ArithMatFloat<VVT,HV1T,MT1,ROWS1>& m1,`
279	`const ArithMatFloat<VV2T,HVT,MT2,ROWS2>& m2,`	279	`const ArithMatFloat<VV2T,HVT,MT2,ROWS2>& m2,`
280	`ArithMatFloat<VVT,HVT,MT,ROWS1>& m)`	280	`ArithMatFloat<VVT,HVT,MT,ROWS1>& m)`
281	`{`	281	`{`
282	`int cols = ArithMatFloat<VVT,HVT,MT,ROWS1>::get_h_dim();`	282	`int cols = ArithMatFloat<VVT,HVT,MT,ROWS1>::get_h_dim();`
283	`for(int i=0;i<ROWS1;i++)`	283	`for(int i=0;i<ROWS1;i++)`
284	`for(int j=0;j<cols;j++)`	284	`for(int j=0;j<cols;j++)`
285	`{`	285	`{`
286	`m[i][j] = 0;`	286	`m[i][j] = 0;`
287	`for(int k=0;k<ROWS2;k++)`	287	`for(int k=0;k<ROWS2;k++)`
288	`m[i][j] += m1[i][k] * m2[k][j];`	288	`m[i][j] += m1[i][k] * m2[k][j];`
289	`}`	289	`}`
290	`}`	290	`}`
291		291
292		292
293	`/** Transpose. See the discussion on mul if you are curious as to why`	293	`/** Transpose. See the discussion on mul if you are curious as to why`
294	`I don't simply return the transpose. */`	294	`I don't simply return the transpose. */`
295	`template <class VVT, class HVT, class M1T, class M2T, int ROWS, int COLS>`	295	`template <class VVT, class HVT, class M1T, class M2T, int ROWS, int COLS>`
296	`inline void transpose(const ArithMatFloat<VVT,HVT,M1T,ROWS>& m,`	296	`inline void transpose(const ArithMatFloat<VVT,HVT,M1T,ROWS>& m,`
297	`ArithMatFloat<HVT,VVT,M2T,COLS>& m_new)`	297	`ArithMatFloat<HVT,VVT,M2T,COLS>& m_new)`
298	`{`	298	`{`
299	`for(int i=0;i<M2T::get_v_dim();i++)`	299	`for(int i=0;i<M2T::get_v_dim();i++)`
300	`for(int j=0;j<M2T::get_h_dim();j++)`	300	`for(int j=0;j<M2T::get_h_dim();j++)`
301	`m_new[i][j] = m[j][i];`	301	`m_new[i][j] = m[j][i];`
302	`}`	302	`}`
303		303
304	`#else`	304	`#else`
305		305
306	`//----------------- win32 -------------------------------`	306	`//----------------- win32 -------------------------------`
307	`// Visual studio is not good at deducing the args. to these template functions.`	307	`// Visual studio is not good at deducing the args. to these template functions.`
308	`// This means that you can call the two functions below with`	308	`// This means that you can call the two functions below with`
309	`// matrices of wrong dimension.`	309	`// matrices of wrong dimension.`
310		310
311	`template <class M1, class M2, class M>`	311	`template <class M1, class M2, class M>`
312	`inline void mul(const M1& m1, const M2& m2, M& m)`	312	`inline void mul(const M1& m1, const M2& m2, M& m)`
313	`{`	313	`{`
314	`int cols = M::get_h_dim();`	314	`int cols = M::get_h_dim();`
315	`int rows1 = M1::get_v_dim();`	315	`int rows1 = M1::get_v_dim();`
316	`int rows2 = M2::get_v_dim();`	316	`int rows2 = M2::get_v_dim();`
317		317
318	`for(int i=0;i<rows1;i++)`	318	`for(int i=0;i<rows1;i++)`
319	`for(int j=0;j<cols;j++)`	319	`for(int j=0;j<cols;j++)`
320	`{`	320	`{`
321	`m[i][j] = 0;`	321	`m[i][j] = 0;`
322	`for(int k=0;k<rows2;k++)`	322	`for(int k=0;k<rows2;k++)`
323	`m[i][j] += m1[i][k] * m2[k][j];`	323	`m[i][j] += m1[i][k] * m2[k][j];`
324	`}`	324	`}`
325	`}`	325	`}`
326		326
327		327
328	`/** Transpose. See the discussion on mul if you are curious as to why`	328	`/** Transpose. See the discussion on mul if you are curious as to why`
329	`I don't simply return the transpose. */`	329	`I don't simply return the transpose. */`
330	`template <class M1, class M2>`	330	`template <class M1, class M2>`
331	`inline void transpose(const M1& m1, M2& m2)`	331	`inline void transpose(const M1& m1, M2& m2)`
332	`{`	332	`{`
333	`for(int i=0;i<M2::get_v_dim();i++)`	333	`for(int i=0;i<M2::get_v_dim();i++)`
334	`for(int j=0;j<M2::get_h_dim();j++)`	334	`for(int j=0;j<M2::get_h_dim();j++)`
335	`m2[i][j] = m1[j][i];`	335	`m2[i][j] = m1[j][i];`
336	`}`	336	`}`
337		337
338	`#endif`	338	`#endif`
339		339
340	`/** Compute the outer product of a and b: a * transpose(b). This is`	340	`/** Compute the outer product of a and b: a * transpose(b). This is`
341	`a matrix with a::rows and b::columns. */`	341	`a matrix with a::rows and b::columns. */`
342	`template <class VVT, class HVT, class MT, int ROWS>`	342	`template <class VVT, class HVT, class MT, int ROWS>`
343	`void outer_product(const VVT& a, const HVT& b,`	343	`void outer_product(const VVT& a, const HVT& b,`
344	`ArithMatFloat<VVT,HVT,MT,ROWS>& m)`	344	`ArithMatFloat<VVT,HVT,MT,ROWS>& m)`
345	`{`	345	`{`
346	`int R = VVT::get_dim();`	346	`int R = VVT::get_dim();`
347	`int C = HVT::get_dim();`	347	`int C = HVT::get_dim();`
348	`for(int i=0;i<R;++i)`	348	`for(int i=0;i<R;++i)`
349	`for(int j=0;j<C;++j)`	349	`for(int j=0;j<C;++j)`
350	`{`	350	`{`
351	`m[i][j] = a[i] * b[j];`	351	`m[i][j] = a[i] * b[j];`
352	`}`	352	`}`
353	`}`	353	`}`
354		354
355	`/** Copy a matrix to another matrix, cell by cell.`	355	`/** Copy a matrix to another matrix, cell by cell.`
356	`This conversion that takes a const matrix as first argument`	356	`This conversion that takes a const matrix as first argument`
357	`(source) and a non-const matrix as second argument`	357	`(source) and a non-const matrix as second argument`
358	`(destination). The contents of the first matrix is simply copied`	358	`(destination). The contents of the first matrix is simply copied`
359	`to the second matrix.`	359	`to the second matrix.`
360		360
361	`However, if the first matrix is larger than the second,`	361	`However, if the first matrix is larger than the second,`
362	`the cells outside the range of the destination are simply not`	362	`the cells outside the range of the destination are simply not`
363	`copied. If the destination is larger, the cells outside the`	363	`copied. If the destination is larger, the cells outside the`
364	`range of the source matrix are not touched.`	364	`range of the source matrix are not touched.`
365		365
366	`An obvious use of this function is to copy a 3x3 rotation matrix`	366	`An obvious use of this function is to copy a 3x3 rotation matrix`
367	`into a 4x4 transformation matrix.`	367	`into a 4x4 transformation matrix.`
368	`*/`	368	`*/`
369		369
370	`template <class M1, class M2>`	370	`template <class M1, class M2>`
371	`void copy_matrix(const M1& inmat, M2& outmat)`	371	`void copy_matrix(const M1& inmat, M2& outmat)`
372	`{`	372	`{`
373	`const int R = s_min(inmat.get_v_dim(), outmat.get_v_dim());`	373	`const int R = s_min(inmat.get_v_dim(), outmat.get_v_dim());`
374	`const int C = s_min(inmat.get_h_dim(), outmat.get_h_dim());`	374	`const int C = s_min(inmat.get_h_dim(), outmat.get_h_dim());`
375	`for(int i=0;i<R;++i)`	375	`for(int i=0;i<R;++i)`
376	`for(int j=0;j<C;++j)`	376	`for(int j=0;j<C;++j)`
377	`outmat[i][j] = inmat[i][j];`	377	`outmat[i][j] = inmat[i][j];`
378	`}`	378	`}`
379		379
380	`/** Put to operator */`	380	`/** Put to operator */`
381	`template <class VVT, class HVT, class MT, int ROWS>`	381	`template <class VVT, class HVT, class MT, int ROWS>`
382	`inline std::ostream&`	382	`inline std::ostream&`
383	`operator<<(std::ostream&os, const ArithMatFloat<VVT,HVT,MT,ROWS>& m)`	383	`operator<<(std::ostream&os, const ArithMatFloat<VVT,HVT,MT,ROWS>& m)`
384	`{`	384	`{`
385	`os << "[\n";`	385	`os << "[\n";`
386	`for(int i=0;i<ROWS;i++) os << " " << m[i] << "\n";`	386	`for(int i=0;i<ROWS;i++) os << " " << m[i] << "\n";`
387	`os << "]\n";`	387	`os << "]\n";`
388	`return os;`	388	`return os;`
389	`}`	389	`}`
390		390
391	`/** Get from operator */`	391	`/** Get from operator */`
392	`template <class VVT, class HVT, class MT, int ROWS>`	392	`template <class VVT, class HVT, class MT, int ROWS>`
393	`inline std::istream& operator>>(std::istream&is,`	393	`inline std::istream& operator>>(std::istream&is,`
394	`const ArithMatFloat<VVT,HVT,MT,ROWS>& m)`	394	`const ArithMatFloat<VVT,HVT,MT,ROWS>& m)`
395	`{`	395	`{`
396	`for(int i=0;i<ROWS;i++) is>>m[i];`	396	`for(int i=0;i<ROWS;i++) is>>m[i];`
397	`return is;`	397	`return is;`
398	`}`	398	`}`
399	`}`	399	`}`
400	`#endif`	400	`#endif`
401		401

Subversion Repositories gelsvn

(root)/branches/imesh-and-hmesh/src/CGLA/ArithMatFloat.h – Rev 12 → 22