Theoretica
A C++ numerical and automatic mathematical library
transform.h
Go to the documentation of this file.
1 
5 
6 #ifndef THEORETICA_TRANSFORM_H
7 #define THEORETICA_TRANSFORM_H
8 
9 #include "./algebra.h"
10 #include "../core/error.h"
11 
12 
13 namespace theoretica {
14 
15  namespace algebra {
16 
20  template<typename Matrix>
21  inline Matrix identity(unsigned int rows = 0, unsigned int cols = 0) {
22 
23  Matrix m;
24  if(rows && cols)
25  m.resize(rows, cols);
26 
27  make_identity(m);
28  return m;
29  }
30 
31 
38  template<typename Vector, typename Matrix>
39  inline Matrix& diagonal(Matrix& res, const Vector& v) {
40 
41  if(v.size() != res.cols()) {
42  TH_MATH_ERROR("algebra::mat_diagonal", v.size(), INVALID_ARGUMENT);
43  mat_error(res);
44  return res;
45  }
46 
47  if(v.size() != res.rows()) {
48  TH_MATH_ERROR("algebra::mat_diagonal", v.size(), INVALID_ARGUMENT);
49  mat_error(res);
50  return res;
51  }
52 
53  for (unsigned int i = 0; i < res.rows(); ++i)
54  for (unsigned int j = 0; j < res.cols(); ++j)
55  res(i, j) = (i == j) ? v[i] : 0;
56 
57  return res;
58  }
59 
60 
69  template<typename Matrix, typename Vector>
70  inline Matrix diagonal(
71  const Vector& v, unsigned int rows = 0, unsigned int cols = 0) {
72 
73  Matrix m;
74  if(rows && cols)
75  m.resize(rows, cols);
76 
77  diagonal(m, v);
78  return m;
79  }
80 
81 
87  template<typename Matrix, typename Vector>
88  inline Matrix translation(
89  const Vector& v, unsigned int rows = 0, unsigned int cols = 0) {
90 
91  Matrix m;
92  if(rows && cols)
93  m.resize(rows, cols);
94 
95  if(v.size() != (m.rows() - 1)) {
96  TH_MATH_ERROR("algebra::translation", v.size(), INVALID_ARGUMENT);
97  mat_error(m);
98  return m;
99  }
100 
101  make_identity(m);
102 
103  // The translation matrix in projective
104  // geometry is constructed by setting
105  // the last column's values to the vector,
106  // while other elements form the identity
107  for (unsigned int i = 0; i < m.rows() - 1; ++i)
108  m(i, m.cols() - 1) = v[i];
109 
110  return m;
111  }
112 
113 
117  template<typename Matrix>
118  inline Matrix rotation_2d(
119  real theta, unsigned int rows = 0, unsigned int cols = 0) {
120 
121  Matrix m;
122  if(rows && cols)
123  m.resize(rows, cols);
124 
125  if(m.rows() < 2) {
126  TH_MATH_ERROR("algebra::rotation_2d", m.rows(), INVALID_ARGUMENT);
127  mat_error(m);
128  return m;
129  }
130 
131  if(m.cols() < 2) {
132  TH_MATH_ERROR("algebra::rotation_2d", m.cols(), INVALID_ARGUMENT);
133  mat_error(m);
134  return m;
135  }
136 
137  const real s = sin(theta);
138  const real c = cos(theta);
139 
140  if(m.rows() > 2 || m.cols() > 2)
141  make_identity(m);
142 
143  m(0, 0) = c;
144  m(0, 1) = -s;
145  m(1, 0) = s;
146  m(1, 1) = c;
147 
148  return m;
149  }
150 
151 
158  template<typename Matrix, typename Vector>
159  inline Matrix rotation_3d(
160  real theta, const Vector& axis, unsigned int rows = 0, unsigned int cols = 0) {
161 
162  Matrix m;
163  if(rows && cols)
164  m.resize(rows, cols);
165 
166  if(axis.size() < 3) {
167  TH_MATH_ERROR("algebra::rotation_3d", axis.size(), INVALID_ARGUMENT);
168  mat_error(m);
169  return m;
170  }
171 
172  if(m.rows() < 3) {
173  TH_MATH_ERROR("algebra::rotation_3d", m.rows(), INVALID_ARGUMENT);
174  mat_error(m);
175  return m;
176  }
177 
178  if(m.cols() < 3) {
179  TH_MATH_ERROR("algebra::rotation_3d", m.cols(), INVALID_ARGUMENT);
180  mat_error(m);
181  return m;
182  }
183 
184  if(m.rows() > 3 || m.cols() > 3)
185  make_identity(m);
186 
187  const real s = sin(theta);
188  const real c = cos(theta);
189 
190  const real Rx = (real) axis[0];
191  const real Ry = (real) axis[1];
192  const real Rz = (real) axis[2];
193 
194  const real cm1 = (1 - c);
195 
196  m(0, 0) = c + Rx * Rx * cm1;
197  m(0, 1) = Rx * Ry * cm1 - Rz * s;
198  m(0, 2) = Rx * Rz * cm1 - Ry * s;
199 
200  m(1, 0) = Ry * Rx * cm1 + Rz * s;
201  m(1, 1) = c + Ry * Ry * cm1;
202  m(1, 2) = Ry * Rz * cm1 - Rx * s;
203 
204  m(2, 0) = Rz * Rx * cm1 - Ry * s;
205  m(2, 1) = Rz * Ry * cm1 + Rx * s;
206  m(2, 2) = c + Rz * Rz * cm1;
207 
208  return m;
209  }
210 
211 
217  template<typename Matrix>
218  inline Matrix rotation_3d_xaxis(
219  real theta, unsigned int rows = 0, unsigned int cols = 0) {
220 
221  Matrix m;
222  if(rows && cols)
223  m.resize(rows, cols);
224 
225  if(m.rows() < 3) {
226  TH_MATH_ERROR("algebra::rotation_3d_xaxis", m.rows(), INVALID_ARGUMENT);
227  mat_error(m);
228  return m;
229  }
230 
231  if(m.cols() < 3) {
232  TH_MATH_ERROR("algebra::rotation_3d_xaxis", m.cols(), INVALID_ARGUMENT);
233  mat_error(m);
234  return m;
235  }
236 
237  if(m.rows() > 3 || m.cols() > 3)
238  make_identity(m);
239 
240  const real s = sin(theta);
241  const real c = cos(theta);
242 
243  m(0, 0) = 1;
244  m(1, 1) = c;
245  m(2, 2) = c;
246 
247  m(1, 2) = -s;
248  m(2, 1) = s;
249 
250  return m;
251  }
252 
253 
259  template<typename Matrix>
260  inline Matrix rotation_3d_yaxis(
261  real theta, unsigned int rows = 0, unsigned int cols = 0) {
262 
263  Matrix m;
264  if(rows && cols)
265  m.resize(rows, cols);
266 
267  if(m.rows() < 3) {
268  TH_MATH_ERROR("algebra::rotation_3d_yaxis", m.rows(), INVALID_ARGUMENT);
269  mat_error(m);
270  return m;
271  }
272 
273  if(m.cols() < 3) {
274  TH_MATH_ERROR("algebra::rotation_3d_yaxis", m.cols(), INVALID_ARGUMENT);
275  mat_error(m);
276  return m;
277  }
278 
279  if(m.rows() > 3 || m.cols() > 3)
280  make_identity(m);
281 
282  const real s = sin(theta);
283  const real c = cos(theta);
284 
285  m(0, 0) = c;
286  m(1, 1) = 1;
287  m(2, 2) = c;
288  m(3, 3) = 1;
289 
290  m(0, 2) = s;
291  m(2, 0) = -s;
292 
293  return m;
294  }
295 
296 
302  template<typename Matrix>
303  inline Matrix rotation_3d_zaxis(
304  real theta, unsigned int rows = 0, unsigned int cols = 0) {
305 
306  Matrix m;
307  if(rows && cols)
308  m.resize(rows, cols);
309 
310  if(m.rows() < 3) {
311  TH_MATH_ERROR("algebra::rotation_3d_zaxis", m.rows(), INVALID_ARGUMENT);
312  mat_error(m);
313  return m;
314  }
315 
316  if(m.cols() < 3) {
317  TH_MATH_ERROR("algebra::rotation_3d_zaxis", m.cols(), INVALID_ARGUMENT);
318  mat_error(m);
319  return m;
320  }
321 
322  if(m.rows() > 3 || m.cols() > 3)
323  make_identity(m);
324 
325  const real s = sin(theta);
326  const real c = cos(theta);
327 
328  m(0, 0) = c;
329  m(1, 1) = c;
330  m(2, 2) = 1;
331  m(3, 3) = 1;
332 
333  m(0, 1) = -s;
334  m(1, 0) = s;
335 
336  return m;
337  }
338 
339 
341  template<typename Matrix>
342  inline Matrix perspective(
343  real left, real right, real bottom,
344  real top, real near, real far,
345  unsigned int rows = 0, unsigned int cols = 0) {
346 
347  Matrix m;
348  if(rows && cols)
349  m.resize(rows, cols);
350 
351  if(m.rows() < 4) {
352  TH_MATH_ERROR("algebra::perspective", m.rows(), INVALID_ARGUMENT);
353  mat_error(m);
354  return m;
355  }
356 
357  if(m.cols() < 4) {
358  TH_MATH_ERROR("algebra::perspective", m.cols(), INVALID_ARGUMENT);
359  mat_error(m);
360  return m;
361  }
362 
363  mat_zeroes(m);
364 
365  m(0, 0) = 2 * near / (right - left);
366  m(2, 0) = (right + left) / (right - left);
367  m(1, 1) = 2 * near / (top - bottom);
368  m(2, 1) = (top + bottom) / (top - bottom);
369  m(2, 2) = -(far + near) / (far - near);
370  m(3, 2) = -(2 * far * near) / (far - near);
371  m(2, 3) = -1;
372  m(3, 3) = 0;
373 
374  return m;
375  }
376 
377 
380  template<typename Matrix>
381  inline Matrix perspective_fov(
382  real fov, real aspect, real near, real far,
383  unsigned int rows = 0, unsigned int cols = 0) {
384 
385  Matrix m;
386  if(rows && cols)
387  m.resize(rows, cols);
388 
389  if(m.rows() < 4) {
390  TH_MATH_ERROR("algebra::perspective_fov", m.rows(), INVALID_ARGUMENT);
391  mat_error(m);
392  return m;
393  }
394 
395  if(m.cols() < 4) {
396  TH_MATH_ERROR("algebra::perspective_fov", m.cols(), INVALID_ARGUMENT);
397  mat_error(m);
398  return m;
399  }
400 
401  const real height = near * tan(radians(fov / 2.f));
402  const real width = height * aspect;
403 
404  return perspective<Matrix>(-width, width, -height, height, near, far);
405  }
406 
407 
409  template<typename Matrix>
410  inline Matrix ortho(
411  real left, real right, real bottom,
412  real top, real near, real far,
413  unsigned int rows = 0, unsigned int cols = 0) {
414 
415  Matrix m;
416  if(rows && cols)
417  m.resize(rows, cols);
418 
419  if(m.rows() < 4) {
420  TH_MATH_ERROR("algebra::ortho", m.rows(), INVALID_ARGUMENT);
421  mat_error(m);
422  return m;
423  }
424 
425  if(m.cols() < 4) {
426  TH_MATH_ERROR("algebra::ortho", m.cols(), INVALID_ARGUMENT);
427  mat_error(m);
428  return m;
429  }
430 
431  mat_zeroes(m);
432 
433  m(0, 0) = 2 / (right - left);
434  m(3, 0) = -(right + left) / (right - left);
435  m(1, 1) = 2 / (top - bottom);
436  m(3, 1) = -(top + bottom) / (top - bottom);
437  m(2, 2) = -2 / (far - near);
438  m(3, 2) = -(far + near) / (far - near);
439 
440  return m;
441  }
442 
443 
446  template<
447  typename Matrix, typename Vector1,
448  typename Vector2, typename Vector3>
449  inline Matrix look_at(Vector1 camera, Vector2 target, Vector3 up) {
450 
451  // Construct an orthonormal basis
452 
453  Vector1 x_axis, y_axis, z_axis;
454  x_axis.resize(3);
455  y_axis.resize(3);
456  z_axis.resize(3);
457 
458  // z = target - camera
459  vec_diff(z_axis, target, camera);
460  make_normalized(z_axis);
461 
462  // x = z X up
463  vec_copy(x_axis, cross(z_axis, up));
464  make_normalized(x_axis);
465 
466  // y = x X z
467  vec_copy(y_axis, cross(x_axis, z_axis));
468 
469  // Negate z_axis to have a right-handed system
470  vec_scalmul(-1.0, z_axis);
471 
472  // Construct the rotation and translation matrix
473  Matrix m;
474  m.resize(4, 4);
475 
476  m(0, 0) = x_axis[1];
477  m(0, 1) = x_axis[2];
478  m(0, 2) = x_axis[3];
479  m(0, 3) = dot(camera, vec_scalmul(-1.0, x_axis));
480 
481  m(1, 0) = y_axis[1];
482  m(1, 1) = y_axis[2];
483  m(1, 2) = y_axis[3];
484  m(1, 3) = dot(camera, vec_scalmul(-1.0, y_axis));
485 
486  m(2, 0) = z_axis[1];
487  m(2, 1) = z_axis[2];
488  m(2, 2) = z_axis[3];
489  m(2, 3) = dot(camera, vec_scalmul(-1.0, z_axis));
490 
491  m(3, 0) = 0;
492  m(3, 1) = 0;
493  m(3, 2) = 0;
494  m(3, 3) = 1;
495 
496  return m;
497  }
498 
499 
501  template<typename Matrix>
502  inline Matrix symplectic(unsigned int rows = 0, unsigned int cols = 0) {
503 
504  Matrix m;
505  if(rows && cols)
506  m.resize(rows, cols);
507 
508  if(rows != cols || (rows % 2 != 0)) {
509  TH_MATH_ERROR("algebra::symplectic", rows, INVALID_ARGUMENT);
510  mat_error(m);
511  return m;
512  }
513 
514  const unsigned int half = m.rows() / 2;
515  mat_zeroes(m);
516 
517  for (unsigned int i = 0; i < half; ++i) {
518  m(i, i + half) = 1;
519  m(i + half, i) = -1;
520  }
521 
522  return m;
523  }
524 
525 
534  template<typename Matrix>
535  inline Matrix hilbert(unsigned int rows = 0) {
536 
537  Matrix H;
538  if (rows)
539  H.resize(rows, rows);
540 
541  for (unsigned int i = 0; i < H.rows(); ++i)
542  for (unsigned int j = 0; j < H.cols(); ++j)
543  H(i, j) = 1.0 / (i + j + 1);
544 
545  return H;
546  }
547 
548 
554  template<typename Vector1, typename Vector2>
555  inline Vector1 sphere_inversion(
556  const Vector1& p, const Vector2& c = Vector2(0), real r = 1) {
557 
558  Vector1 q = p - c;
559  return c + q * square(r / q.norm());
560  }
561  }
562 }
563 
564 #endif
algebra.h
Linear algebra routines.
TH_MATH_ERROR
#define TH_MATH_ERROR(F_NAME, VALUE, EXCEPTION)
TH_MATH_ERROR is a macro which throws exceptions or modifies errno (depending on which compiling opti...
Definition: error.h:219
theoretica::algebra::mat_zeroes
Matrix & mat_zeroes(Matrix &m)
Overwrite a matrix with all zeroes.
Definition: algebra.h:93
theoretica::algebra::symplectic
Matrix symplectic(unsigned int rows=0, unsigned int cols=0)
Generate a NxN symplectic matrix where N is even.
Definition: transform.h:502
theoretica::algebra::cross
Vector1 cross(const Vector1 &v1, const Vector2 &v2)
Compute the cross product between two tridimensional vectors.
Definition: algebra.h:373
theoretica::algebra::make_identity
Matrix & make_identity(Matrix &m)
Overwrite a matrix with the identity matrix.
Definition: algebra.h:77
theoretica::algebra::make_normalized
Vector & make_normalized(Vector &v)
Normalize a given vector overwriting it.
Definition: algebra.h:327
theoretica::algebra::look_at
Matrix look_at(Vector1 camera, Vector2 target, Vector3 up)
Return a transformation matrix that points the field of view towards a given point from the <camera> ...
Definition: transform.h:449
theoretica::algebra::rotation_2d
Matrix rotation_2d(real theta, unsigned int rows=0, unsigned int cols=0)
Returns a matrix representing a 2D rotation.
Definition: transform.h:118
theoretica::algebra::sphere_inversion
Vector1 sphere_inversion(const Vector1 &p, const Vector2 &c=Vector2(0), real r=1)
Sphere inversion of a point with respect to a sphere of radius r centered at a point c.
Definition: transform.h:555
theoretica::algebra::hilbert
Matrix hilbert(unsigned int rows=0)
Construct the Hilbert matrix of arbitrary dimension.
Definition: transform.h:535
theoretica::algebra::dot
auto dot(const Vector1 &v, const Vector2 &w)
Computes the dot product between two vectors, using the Hermitian form if needed.
Definition: algebra.h:351
theoretica::algebra::rotation_3d_xaxis
Matrix rotation_3d_xaxis(real theta, unsigned int rows=0, unsigned int cols=0)
Returns a matrix representing a 3D rotation around the x axis.
Definition: transform.h:218
theoretica::algebra::mat_error
Matrix & mat_error(Matrix &m)
Overwrite the given matrix with the error matrix with NaN values on the diagonal and zeroes everywher...
Definition: algebra.h:44
theoretica::algebra::diagonal
Matrix & diagonal(Matrix &res, const Vector &v)
Construct a matrix with the given vector as diagonal and zeroes everywhere else, overwriting the give...
Definition: transform.h:39
theoretica::algebra::vec_scalmul
Vector & vec_scalmul(Field a, Vector &v)
Multiply a vector by a scalar of any compatible type.
Definition: algebra.h:634
theoretica::algebra::perspective_fov
Matrix perspective_fov(real fov, real aspect, real near, real far, unsigned int rows=0, unsigned int cols=0)
Returns a perspective projection matrix, using the Field of View as parameter.
Definition: transform.h:381
theoretica::algebra::rotation_3d_yaxis
Matrix rotation_3d_yaxis(real theta, unsigned int rows=0, unsigned int cols=0)
Returns a matrix representing a 3D rotation around the y axis.
Definition: transform.h:260
theoretica::algebra::rotation_3d
Matrix rotation_3d(real theta, const Vector &axis, unsigned int rows=0, unsigned int cols=0)
Returns a matrix representing a 3D rotation around a given axis.
Definition: transform.h:159
theoretica::algebra::vec_copy
Vector1 & vec_copy(Vector1 &dest, const Vector2 &src)
Copy a vector by overwriting another.
Definition: algebra.h:143
theoretica::algebra::vec_diff
Vector2 & vec_diff(Vector1 &v1, const Vector2 &v2)
Subtract two vectors and store the result in the first vector.
Definition: algebra.h:1183
theoretica::algebra::ortho
Matrix ortho(real left, real right, real bottom, real top, real near, real far, unsigned int rows=0, unsigned int cols=0)
Returns an orthogonal projection matrix.
Definition: transform.h:410
theoretica::algebra::identity
Matrix identity(unsigned int rows=0, unsigned int cols=0)
Returns the identity matrix.
Definition: transform.h:21
theoretica::algebra::rotation_3d_zaxis
Matrix rotation_3d_zaxis(real theta, unsigned int rows=0, unsigned int cols=0)
Returns a matrix representing a 3D rotation around the z axis.
Definition: transform.h:303
theoretica::algebra::perspective
Matrix perspective(real left, real right, real bottom, real top, real near, real far, unsigned int rows=0, unsigned int cols=0)
Returns a perspective projection matrix.
Definition: transform.h:342
theoretica::algebra::translation
Matrix translation(const Vector &v, unsigned int rows=0, unsigned int cols=0)
Returns a translation matrix, that is, an NxN matrix which describes a translation in N-1 dimensions.
Definition: transform.h:88
theoretica
Main namespace of the library which contains all functions and objects.
Definition: algebra.h:27
theoretica::real
double real
A real number, defined as a floating point type.
Definition: constants.h:188
theoretica::radians
TH_CONSTEXPR real radians(real degrees)
Convert degrees to radians.
Definition: real_analysis.h:1253
theoretica::cos
dual2 cos(dual2 x)
Compute the cosine of a second order dual number.
Definition: dual2_functions.h:82
theoretica::tan
dual2 tan(dual2 x)
Compute the tangent of a second order dual number.
Definition: dual2_functions.h:94
theoretica::sin
dual2 sin(dual2 x)
Compute the sine of a second order dual number.
Definition: dual2_functions.h:70
theoretica::square
dual2 square(dual2 x)
Return the square of a second order dual number.
Definition: dual2_functions.h:23