← Back to branch summary

~siretart/ubuntu/utopic/blender/libav10

« back to all changes in this revision

Viewing changes to extern/libmv/libmv/numeric/numeric.h

  • Committer: Package Import Robot
  • Author(s): Matteo F. Vescovi
  • Date: 2012-07-23 08:54:18 UTC
  • mfrom: (14.2.16 sid)
  • mto: (14.2.19 sid)
  • mto: This revision was merged to the branch mainline in revision 42.
  • Revision ID: package-import@ubuntu.com-20120723085418-9foz30v6afaf5ffs
Tags: 2.63a-2
http://bugs.debian.org/658075
* debian/: Cycles support added (Closes: #658075)
  For now, this top feature has been enabled only
  on [any-amd64 any-i386] architectures because
  of OpenImageIO failing on all others
* debian/: scripts installation path changed
  from /usr/lib to /usr/share:
  + debian/patches/: patchset re-worked for path changing
  + debian/control: "Breaks" field added on yafaray-exporter

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/libmv/libmv/numeric/numeric.h
 
1
// Copyright (c) 2007, 2008 libmv authors.
 
2
//
 
3
// Permission is hereby granted, free of charge, to any person obtaining a copy
 
4
// of this software and associated documentation files (the "Software"), to
 
5
// deal in the Software without restriction, including without limitation the
 
6
// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
 
7
// sell copies of the Software, and to permit persons to whom the Software is
 
8
// furnished to do so, subject to the following conditions:
 
9
//
 
10
// The above copyright notice and this permission notice shall be included in
 
11
// all copies or substantial portions of the Software.
 
12
//
 
13
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 
14
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 
15
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 
16
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 
17
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 
18
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
 
19
// IN THE SOFTWARE.
 
20
//
 
21
// Matrix and vector classes, based on Eigen2.
 
22
//
 
23
// Avoid using Eigen2 classes directly; instead typedef them here.
 
24
 
 
25
#ifndef LIBMV_NUMERIC_NUMERIC_H
 
26
#define LIBMV_NUMERIC_NUMERIC_H
 
27
 
 
28
#include <Eigen/Cholesky>
 
29
#include <Eigen/Core>
 
30
#include <Eigen/Eigenvalues>
 
31
#include <Eigen/Geometry>
 
32
#include <Eigen/LU>
 
33
#include <Eigen/QR>
 
34
#include <Eigen/SVD>
 
35
 
 
36
#if (defined(_WIN32) || defined(__APPLE__) || defined(__FreeBSD__)) && !defined(__MINGW64__)
 
37
  void static sincos (double x, double *sinx, double *cosx) {
 
38
    *sinx = sin(x);
 
39
    *cosx = cos(x);
 
40
  }
 
41
#endif //_WIN32 || __APPLE__
 
42
 
 
43
#if (defined(WIN32) || defined(WIN64)) && !defined(__MINGW32__)
 
44
  inline long lround(double d) {
 
45
    return (long)(d>0 ? d+0.5 : ceil(d-0.5));
 
46
  }
 
47
  inline int round(double d) {
 
48
    return (d>0) ? int(d+0.5) : int(d-0.5);
 
49
  }
 
50
  typedef unsigned int uint;
 
51
#endif //_WIN32
 
52
 
 
53
namespace libmv {
 
54
 
 
55
typedef Eigen::MatrixXd Mat;
 
56
typedef Eigen::VectorXd Vec;
 
57
 
 
58
typedef Eigen::MatrixXf Matf;
 
59
typedef Eigen::VectorXf Vecf;
 
60
 
 
61
typedef Eigen::Matrix<unsigned int, Eigen::Dynamic, Eigen::Dynamic> Matu;
 
62
typedef Eigen::Matrix<unsigned int, Eigen::Dynamic, 1> Vecu;
 
63
typedef Eigen::Matrix<unsigned int, 2, 1> Vec2u;
 
64
 
 
65
typedef Eigen::Matrix<double, 2, 2> Mat2;
 
66
typedef Eigen::Matrix<double, 2, 3> Mat23;
 
67
typedef Eigen::Matrix<double, 3, 3> Mat3;
 
68
typedef Eigen::Matrix<double, 3, 4> Mat34;
 
69
typedef Eigen::Matrix<double, 3, 5> Mat35;
 
70
typedef Eigen::Matrix<double, 4, 1> Mat41;
 
71
typedef Eigen::Matrix<double, 4, 3> Mat43;
 
72
typedef Eigen::Matrix<double, 4, 4> Mat4;
 
73
typedef Eigen::Matrix<double, 4, 6> Mat46;
 
74
typedef Eigen::Matrix<float, 2, 2> Mat2f;
 
75
typedef Eigen::Matrix<float, 2, 3> Mat23f;
 
76
typedef Eigen::Matrix<float, 3, 3> Mat3f;
 
77
typedef Eigen::Matrix<float, 3, 4> Mat34f;
 
78
typedef Eigen::Matrix<float, 3, 5> Mat35f;
 
79
typedef Eigen::Matrix<float, 4, 3> Mat43f;
 
80
typedef Eigen::Matrix<float, 4, 4> Mat4f;
 
81
typedef Eigen::Matrix<float, 4, 6> Mat46f;
 
82
 
 
83
typedef Eigen::Matrix<double, 3, 3, Eigen::RowMajor> RMat3;
 
84
typedef Eigen::Matrix<double, 4, 4, Eigen::RowMajor> RMat4;
 
85
 
 
86
typedef Eigen::Matrix<double, 2, Eigen::Dynamic> Mat2X;
 
87
typedef Eigen::Matrix<double, 3, Eigen::Dynamic> Mat3X;
 
88
typedef Eigen::Matrix<double, 4, Eigen::Dynamic> Mat4X;
 
89
typedef Eigen::Matrix<double, Eigen::Dynamic, 2> MatX2;
 
90
typedef Eigen::Matrix<double, Eigen::Dynamic, 3> MatX3;
 
91
typedef Eigen::Matrix<double, Eigen::Dynamic, 4> MatX4;
 
92
typedef Eigen::Matrix<double, Eigen::Dynamic, 5> MatX5;
 
93
typedef Eigen::Matrix<double, Eigen::Dynamic, 6> MatX6;
 
94
typedef Eigen::Matrix<double, Eigen::Dynamic, 7> MatX7;
 
95
typedef Eigen::Matrix<double, Eigen::Dynamic, 8> MatX8;
 
96
typedef Eigen::Matrix<double, Eigen::Dynamic, 9> MatX9;
 
97
typedef Eigen::Matrix<double, Eigen::Dynamic,15> MatX15;
 
98
typedef Eigen::Matrix<double, Eigen::Dynamic,16> MatX16;
 
99
 
 
100
typedef Eigen::Vector2d Vec2;
 
101
typedef Eigen::Vector3d Vec3;
 
102
typedef Eigen::Vector4d Vec4;
 
103
typedef Eigen::Matrix<double, 5, 1>  Vec5;
 
104
typedef Eigen::Matrix<double, 6, 1>  Vec6;
 
105
typedef Eigen::Matrix<double, 7, 1>  Vec7;
 
106
typedef Eigen::Matrix<double, 8, 1>  Vec8;
 
107
typedef Eigen::Matrix<double, 9, 1>  Vec9;
 
108
typedef Eigen::Matrix<double, 10, 1> Vec10;
 
109
typedef Eigen::Matrix<double, 11, 1> Vec11;
 
110
typedef Eigen::Matrix<double, 12, 1> Vec12;
 
111
typedef Eigen::Matrix<double, 13, 1> Vec13;
 
112
typedef Eigen::Matrix<double, 14, 1> Vec14;
 
113
typedef Eigen::Matrix<double, 15, 1> Vec15;
 
114
typedef Eigen::Matrix<double, 16, 1> Vec16;
 
115
typedef Eigen::Matrix<double, 17, 1> Vec17;
 
116
typedef Eigen::Matrix<double, 18, 1> Vec18;
 
117
typedef Eigen::Matrix<double, 19, 1> Vec19;
 
118
typedef Eigen::Matrix<double, 20, 1> Vec20;
 
119
 
 
120
typedef Eigen::Vector2f Vec2f;
 
121
typedef Eigen::Vector3f Vec3f;
 
122
typedef Eigen::Vector4f Vec4f;
 
123
 
 
124
typedef Eigen::VectorXi VecXi;
 
125
 
 
126
typedef Eigen::Vector2i Vec2i;
 
127
typedef Eigen::Vector3i Vec3i;
 
128
typedef Eigen::Vector4i Vec4i;
 
129
 
 
130
typedef Eigen::Matrix<float,
 
131
                      Eigen::Dynamic,
 
132
                      Eigen::Dynamic,
 
133
                      Eigen::RowMajor> RMatf;
 
134
 
 
135
typedef Eigen::NumTraits<double> EigenDouble;
 
136
 
 
137
using Eigen::Map;
 
138
using Eigen::Dynamic;
 
139
using Eigen::Matrix;
 
140
 
 
141
// Find U, s, and VT such that
 
142
//
 
143
//   A = U * diag(s) * VT
 
144
//
 
145
template <typename TMat, typename TVec>
 
146
inline void SVD(TMat *A, Vec *s, Mat *U, Mat *VT) {
 
147
  assert(0);
 
148
}
 
149
 
 
150
// Solve the linear system Ax = 0 via SVD. Store the solution in x, such that
 
151
// ||x|| = 1.0. Return the singluar value corresponding to the solution.
 
152
// Destroys A and resizes x if necessary.
 
153
// TODO(maclean): Take the SVD of the transpose instead of this zero padding.
 
154
template <typename TMat, typename TVec>
 
155
double Nullspace(TMat *A, TVec *nullspace) {
 
156
  Eigen::JacobiSVD<TMat> svd(*A, Eigen::ComputeFullV);
 
157
  (*nullspace) = svd.matrixV().col(A->cols()-1);
 
158
  if (A->rows() >= A->cols())
 
159
    return svd.singularValues()(A->cols()-1);
 
160
  else
 
161
    return 0.0;
 
162
}
 
163
 
 
164
// Solve the linear system Ax = 0 via SVD. Finds two solutions, x1 and x2, such
 
165
// that x1 is the best solution and x2 is the next best solution (in the L2
 
166
// norm sense). Store the solution in x1 and x2, such that ||x|| = 1.0. Return
 
167
// the singluar value corresponding to the solution x1.  Destroys A and resizes
 
168
// x if necessary.
 
169
template <typename TMat, typename TVec1, typename TVec2>
 
170
double Nullspace2(TMat *A, TVec1 *x1, TVec2 *x2) {
 
171
  Eigen::JacobiSVD<TMat> svd(*A, Eigen::ComputeFullV);
 
172
  *x1 = svd.matrixV().col(A->cols() - 1);
 
173
  *x2 = svd.matrixV().col(A->cols() - 2);
 
174
  if (A->rows() >= A->cols())
 
175
    return svd.singularValues()(A->cols()-1);
 
176
  else
 
177
    return 0.0;
 
178
}
 
179
 
 
180
// In place transpose for square matrices.
 
181
template<class TA>
 
182
inline void TransposeInPlace(TA *A) {
 
183
  *A = A->transpose().eval();
 
184
}
 
185
 
 
186
template<typename TVec>
 
187
inline double NormL1(const TVec &x) {
 
188
  return x.array().abs().sum();
 
189
}
 
190
 
 
191
template<typename TVec>
 
192
inline double NormL2(const TVec &x) {
 
193
  return x.norm();
 
194
}
 
195
 
 
196
template<typename TVec>
 
197
inline double NormLInfinity(const TVec &x) {
 
198
  return x.array().abs().maxCoeff();
 
199
}
 
200
 
 
201
template<typename TVec>
 
202
inline double DistanceL1(const TVec &x, const TVec &y) {
 
203
  return (x - y).array().abs().sum();
 
204
}
 
205
 
 
206
template<typename TVec>
 
207
inline double DistanceL2(const TVec &x, const TVec &y) {
 
208
  return (x - y).norm();
 
209
}
 
210
template<typename TVec>
 
211
inline double DistanceLInfinity(const TVec &x, const TVec &y) {
 
212
  return (x - y).array().abs().maxCoeff();
 
213
}
 
214
 
 
215
// Normalize a vector with the L1 norm, and return the norm before it was
 
216
// normalized.
 
217
template<typename TVec>
 
218
inline double NormalizeL1(TVec *x) {
 
219
  double norm = NormL1(*x);
 
220
  *x /= norm;
 
221
  return norm;
 
222
}
 
223
 
 
224
// Normalize a vector with the L2 norm, and return the norm before it was
 
225
// normalized.
 
226
template<typename TVec>
 
227
inline double NormalizeL2(TVec *x) {
 
228
  double norm = NormL2(*x);
 
229
  *x /= norm;
 
230
  return norm;
 
231
}
 
232
 
 
233
// Normalize a vector with the L^Infinity norm, and return the norm before it
 
234
// was normalized.
 
235
template<typename TVec>
 
236
inline double NormalizeLInfinity(TVec *x) {
 
237
  double norm = NormLInfinity(*x);
 
238
  *x /= norm;
 
239
  return norm;
 
240
}
 
241
 
 
242
// Return the square of a number.
 
243
template<typename T>
 
244
inline T Square(T x) {
 
245
  return x * x;
 
246
}
 
247
 
 
248
Mat3 RotationAroundX(double angle);
 
249
Mat3 RotationAroundY(double angle);
 
250
Mat3 RotationAroundZ(double angle);
 
251
 
 
252
// Returns the rotation matrix of a rotation of angle |axis| around axis.
 
253
// This is computed using the Rodrigues formula, see:
 
254
//  http://mathworld.wolfram.com/RodriguesRotationFormula.html
 
255
Mat3 RotationRodrigues(const Vec3 &axis);
 
256
 
 
257
// Make a rotation matrix such that center becomes the direction of the
 
258
// positive z-axis, and y is oriented close to up.
 
259
Mat3 LookAt(Vec3 center);
 
260
 
 
261
// Return a diagonal matrix from a vector containg the diagonal values.
 
262
template <typename TVec>
 
263
inline Mat Diag(const TVec &x) {
 
264
  return x.asDiagonal();
 
265
}
 
266
 
 
267
template<typename TMat>
 
268
inline double FrobeniusNorm(const TMat &A) {
 
269
  return sqrt(A.array().abs2().sum());
 
270
}
 
271
 
 
272
template<typename TMat>
 
273
inline double FrobeniusDistance(const TMat &A, const TMat &B) {
 
274
  return FrobeniusNorm(A - B);
 
275
}
 
276
 
 
277
inline Vec3 CrossProduct(const Vec3 &x, const Vec3 &y) {
 
278
  return x.cross(y);
 
279
}
 
280
 
 
281
Mat3 CrossProductMatrix(const Vec3 &x);
 
282
 
 
283
void MeanAndVarianceAlongRows(const Mat &A,
 
284
                              Vec *mean_pointer,
 
285
                              Vec *variance_pointer);
 
286
 
 
287
#if _WIN32
 
288
  // TODO(bomboze): un-#if this for both platforms once tested under Windows
 
289
  /* This solution was extensively discussed here http://forum.kde.org/viewtopic.php?f=74&t=61940 */
 
290
  #define SUM_OR_DYNAMIC(x,y) (x==Eigen::Dynamic||y==Eigen::Dynamic)?Eigen::Dynamic:(x+y)
 
291
 
 
292
  template<typename Derived1, typename Derived2>
 
293
  struct hstack_return {
 
294
    typedef typename Derived1::Scalar Scalar;
 
295
    enum {
 
296
         RowsAtCompileTime = Derived1::RowsAtCompileTime,
 
297
         ColsAtCompileTime = SUM_OR_DYNAMIC(Derived1::ColsAtCompileTime, Derived2::ColsAtCompileTime),
 
298
         Options = Derived1::Flags&Eigen::RowMajorBit ? Eigen::RowMajor : 0,
 
299
         MaxRowsAtCompileTime = Derived1::MaxRowsAtCompileTime,
 
300
         MaxColsAtCompileTime = SUM_OR_DYNAMIC(Derived1::MaxColsAtCompileTime, Derived2::MaxColsAtCompileTime)
 
301
    };
 
302
    typedef Eigen::Matrix<Scalar,
 
303
                RowsAtCompileTime,
 
304
                ColsAtCompileTime,
 
305
                Options,
 
306
                MaxRowsAtCompileTime,
 
307
                MaxColsAtCompileTime> type;
 
308
  };
 
309
 
 
310
  template<typename Derived1, typename Derived2>
 
311
  typename hstack_return<Derived1,Derived2>::type
 
312
  HStack (const Eigen::MatrixBase<Derived1>& lhs, const Eigen::MatrixBase<Derived2>& rhs) {
 
313
    typename hstack_return<Derived1,Derived2>::type res;
 
314
    res.resize(lhs.rows(), lhs.cols()+rhs.cols());
 
315
    res << lhs, rhs;
 
316
    return res;
 
317
  };
 
318
 
 
319
 
 
320
  template<typename Derived1, typename Derived2>
 
321
  struct vstack_return {
 
322
    typedef typename Derived1::Scalar Scalar;
 
323
    enum {
 
324
         RowsAtCompileTime = SUM_OR_DYNAMIC(Derived1::RowsAtCompileTime, Derived2::RowsAtCompileTime),
 
325
         ColsAtCompileTime = Derived1::ColsAtCompileTime,
 
326
         Options = Derived1::Flags&Eigen::RowMajorBit ? Eigen::RowMajor : 0,
 
327
         MaxRowsAtCompileTime = SUM_OR_DYNAMIC(Derived1::MaxRowsAtCompileTime, Derived2::MaxRowsAtCompileTime),
 
328
         MaxColsAtCompileTime = Derived1::MaxColsAtCompileTime
 
329
    };
 
330
    typedef Eigen::Matrix<Scalar,
 
331
                RowsAtCompileTime,
 
332
                ColsAtCompileTime,
 
333
                Options,
 
334
                MaxRowsAtCompileTime,
 
335
                MaxColsAtCompileTime> type;
 
336
  };
 
337
 
 
338
  template<typename Derived1, typename Derived2>
 
339
  typename vstack_return<Derived1,Derived2>::type
 
340
  VStack (const Eigen::MatrixBase<Derived1>& lhs, const Eigen::MatrixBase<Derived2>& rhs) {
 
341
    typename vstack_return<Derived1,Derived2>::type res;
 
342
    res.resize(lhs.rows()+rhs.rows(), lhs.cols());
 
343
    res << lhs, rhs;
 
344
    return res;
 
345
  };
 
346
 
 
347
 
 
348
#else //_WIN32
 
349
 
 
350
  // Since it is not possible to typedef privately here, use a macro.
 
351
  // Always take dynamic columns if either side is dynamic.
 
352
  #define COLS \
 
353
    ((ColsLeft == Eigen::Dynamic || ColsRight == Eigen::Dynamic) \
 
354
     ? Eigen::Dynamic : (ColsLeft + ColsRight))
 
355
 
 
356
  // Same as above, except that prefer fixed size if either is fixed.
 
357
  #define ROWS \
 
358
    ((RowsLeft == Eigen::Dynamic && RowsRight == Eigen::Dynamic) \
 
359
     ? Eigen::Dynamic \
 
360
     : ((RowsLeft == Eigen::Dynamic) \
 
361
        ? RowsRight \
 
362
        : RowsLeft \
 
363
       ) \
 
364
    )
 
365
 
 
366
  // TODO(keir): Add a static assert if both rows are at compiletime.
 
367
  template<typename T, int RowsLeft, int RowsRight, int ColsLeft, int ColsRight>
 
368
  Eigen::Matrix<T, ROWS, COLS>
 
369
  HStack(const Eigen::Matrix<T, RowsLeft,  ColsLeft>  &left,
 
370
         const Eigen::Matrix<T, RowsRight, ColsRight> &right) {
 
371
    assert(left.rows() == right.rows());
 
372
    int n = left.rows();
 
373
    int m1 = left.cols();
 
374
    int m2 = right.cols();
 
375
 
 
376
    Eigen::Matrix<T, ROWS, COLS> stacked(n, m1 + m2);
 
377
    stacked.block(0, 0,  n, m1) = left;
 
378
    stacked.block(0, m1, n, m2) = right;
 
379
    return stacked;
 
380
  }
 
381
 
 
382
  // Reuse the above macros by swapping the order of Rows and Cols. Nasty, but
 
383
  // the duplication is worse.
 
384
  // TODO(keir): Add a static assert if both rows are at compiletime.
 
385
  // TODO(keir): Mail eigen list about making this work for general expressions
 
386
  // rather than only matrix types.
 
387
  template<typename T, int RowsLeft, int RowsRight, int ColsLeft, int ColsRight>
 
388
  Eigen::Matrix<T, COLS, ROWS>
 
389
  VStack(const Eigen::Matrix<T, ColsLeft,  RowsLeft>  &top,
 
390
         const Eigen::Matrix<T, ColsRight, RowsRight> &bottom) {
 
391
    assert(top.cols() == bottom.cols());
 
392
     int n1 = top.rows();
 
393
    int n2 = bottom.rows();
 
394
    int m = top.cols();
 
395
 
 
396
    Eigen::Matrix<T, COLS, ROWS> stacked(n1 + n2, m);
 
397
    stacked.block(0,  0, n1, m) = top;
 
398
    stacked.block(n1, 0, n2, m) = bottom;
 
399
    return stacked;
 
400
  }
 
401
  #undef COLS
 
402
  #undef ROWS
 
403
#endif //_WIN32
 
404
 
 
405
 
 
406
 
 
407
void HorizontalStack(const Mat &left, const Mat &right, Mat *stacked);
 
408
 
 
409
template<typename TTop, typename TBot, typename TStacked>
 
410
void VerticalStack(const TTop &top, const TBot &bottom, TStacked *stacked) {
 
411
  assert(top.cols() == bottom.cols());
 
412
  int n1 = top.rows();
 
413
  int n2 = bottom.rows();
 
414
  int m = top.cols();
 
415
 
 
416
  stacked->resize(n1 + n2, m);
 
417
  stacked->block(0,  0, n1, m) = top;
 
418
  stacked->block(n1, 0, n2, m) = bottom;
 
419
}
 
420
 
 
421
void MatrixColumn(const Mat &A, int i, Vec2 *v);
 
422
void MatrixColumn(const Mat &A, int i, Vec3 *v);
 
423
void MatrixColumn(const Mat &A, int i, Vec4 *v);
 
424
 
 
425
template <typename TMat, typename TCols>
 
426
TMat ExtractColumns(const TMat &A, const TCols &columns) {
 
427
  TMat compressed(A.rows(), columns.size());
 
428
  for (int i = 0; i < columns.size(); ++i) {
 
429
    compressed.col(i) = A.col(columns[i]);
 
430
  }
 
431
  return compressed;
 
432
}
 
433
 
 
434
template <typename TMat, typename TDest>
 
435
void reshape(const TMat &a, int rows, int cols, TDest *b) {
 
436
  assert(a.rows()*a.cols() == rows*cols);
 
437
  b->resize(rows, cols);
 
438
  for (int i = 0; i < rows; i++) {
 
439
    for (int j = 0; j < cols; j++) {
 
440
      (*b)(i, j) = a[cols*i + j];
 
441
    }
 
442
  }
 
443
}
 
444
 
 
445
inline bool isnan(double i) {
 
446
#ifdef WIN32
 
447
  return _isnan(i) > 0;
 
448
#else
 
449
  return std::isnan(i);
 
450
#endif
 
451
}
 
452
 
 
453
/// Ceil function that has the same behaviour for positive
 
454
/// and negative values
 
455
template <typename FloatType>
 
456
FloatType ceil0(const FloatType& value) {
 
457
    FloatType result = std::ceil( std::fabs( value ) );
 
458
    return (value < 0.0) ? -result : result;
 
459
}
 
460
 
 
461
/// Returns the skew anti-symmetric matrix of a vector
 
462
inline Mat3 SkewMat(const Vec3 &x) {
 
463
  Mat3 skew;
 
464
  skew <<   0 , -x(2),  x(1),
 
465
          x(2),    0 , -x(0),
 
466
         -x(1),  x(0),    0;
 
467
  return skew;
 
468
}
 
469
/// Returns the skew anti-symmetric matrix of a vector with only
 
470
/// the first two (independent) lines
 
471
inline Mat23 SkewMatMinimal(const Vec2 &x) {
 
472
  Mat23 skew;
 
473
  skew << 0,-1,  x(1),
 
474
          1, 0, -x(0);
 
475
  return skew;
 
476
}
 
477
 
 
478
/// Returns the rotaiton matrix built from given vector of euler angles
 
479
inline Mat3 RotationFromEulerVector(Vec3 euler_vector) {
 
480
  double theta = euler_vector.norm();
 
481
  if (theta == 0.0) {
 
482
    return Mat3::Identity();
 
483
  }
 
484
  Vec3 w = euler_vector / theta;
 
485
  Mat3 w_hat = CrossProductMatrix(w);
 
486
  return Mat3::Identity() + w_hat*sin(theta) + w_hat*w_hat*(1 - cos(theta));
 
487
}
 
488
} // namespace libmv
 
489
 
 
490
#endif  // LIBMV_NUMERIC_NUMERIC_H