Subversion Repositories gelsvn

Rev

Rev 10 | Rev 43 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

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