vcglib: eigen.h Source File

00001 /****************************************************************************
00002 * VCGLib                                                            o o     *
00003 * Visual and Computer Graphics Library                            o     o   *
00004 *                                                                _   O  _   *
00005 * Copyright(C) 2004                                                \/)\/    *
00006 * Visual Computing Lab                                            /\/|      *
00007 * ISTI - Italian National Research Council                           |      *
00008 *                                                                    \      *
00009 * All rights reserved.                                                      *
00010 *                                                                           *
00011 * This program is free software; you can redistribute it and/or modify      *
00012 * it under the terms of the GNU General Public License as published by      *
00013 * the Free Software Foundation; either version 2 of the License, or         *
00014 * (at your option) any later version.                                       *
00015 *                                                                           *
00016 * This program is distributed in the hope that it will be useful,           *
00017 * but WITHOUT ANY WARRANTY; without even the implied warranty of            *
00018 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
00019 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt)          *
00020 * for more details.                                                         *
00021 *                                                                           *
00022 ****************************************************************************/
00023 
00024 #ifndef EIGEN_VCGLIB
00025 #define EIGEN_VCGLIB
00026 
00027 // TODO enable the vectorization
00028 // #define EIGEN_DONT_VECTORIZE
00029 #define EIGEN_MATRIXBASE_PLUGIN <vcg/math/eigen_matrixbase_addons.h>
00030 #define EIGEN_MATRIX_PLUGIN <vcg/math/eigen_matrix_addons.h>
00031 
00032 // forward declarations
00033 namespace Eigen {
00034 
00035 template<typename T> struct ei_traits;
00036 template<typename A,typename B> struct ei_is_same_type;
00037 
00038 template<typename Derived1, typename Derived2, int Size> struct ei_lexi_comparison;
00039 
00040 template<typename Derived1, typename Derived2,
00041         bool SameType = ei_is_same_type<Derived1,Derived2>::ret,
00042         bool SameSize = Derived1::SizeAtCompileTime==Derived2::SizeAtCompileTime>
00043 struct ei_import_selector;
00044 
00045 template<typename XprType,
00046         int Rows = ei_traits<XprType>::RowsAtCompileTime,
00047         int Cols = ei_traits<XprType>::ColsAtCompileTime,
00048         int StorageOrder = ei_traits<XprType>::Flags&1,
00049         int MRows = ei_traits<XprType>::MaxRowsAtCompileTime,
00050         int MCols = ei_traits<XprType>::MaxColsAtCompileTime>
00051 struct ei_to_vcgtype;
00052 
00053 }
00054 
00055 #include "base.h"
00056 #include "../Eigen/LU"
00057 #include "../Eigen/Geometry"
00058 #include "../Eigen/Array"
00059 #include "../Eigen/Core"
00060 
00061 // add support for unsigned char and short int
00062 namespace Eigen {
00063 template<> struct NumTraits<unsigned char>
00064 {
00065   typedef unsigned char Real;
00066   typedef float FloatingPoint;
00067   enum {
00068     IsComplex = 0,
00069     HasFloatingPoint = 0,
00070     ReadCost = 1,
00071     AddCost = 1,
00072     MulCost = 1
00073   };
00074 };
00075 
00076 template<> struct NumTraits<short int>
00077 {
00078   typedef short int Real;
00079   typedef float FloatingPoint;
00080   enum {
00081     IsComplex = 0,
00082     HasFloatingPoint = 0,
00083     ReadCost = 1,
00084     AddCost = 1,
00085     MulCost = 1
00086   };
00087 };
00088 
00089 // WARNING this is a default version provided so that Intersection() stuff can compile.
00090 // Indeed, the compiler try to instanciate all versions of Intersection() leading to
00091 // the instanciation of Eigen::Matrix<Face,...> !!!
00092 template<typename T> struct NumTraits
00093 {
00094         struct wrong_type
00095         {
00096                 wrong_type() { assert(0 && "Eigen: you are using a wrong scalar type" ); }
00097         };
00098 
00099         typedef wrong_type Real;
00100         typedef wrong_type FloatingPoint;
00101         enum {
00102                 IsComplex = 0,
00103                 HasFloatingPoint = 0,
00104                 ReadCost = 0,
00105                 AddCost = 0,
00106                 MulCost = 0
00107         };
00108 };
00109 
00110 // implementation of Lexicographic order comparison
00111 // TODO should use meta unrollers
00112 template<typename Derived1, typename Derived2> struct ei_lexi_comparison<Derived1,Derived2,2>
00113 {
00114         inline static bool less(const Derived1& a, const Derived2& b) {
00115                 return (a.coeff(1)!=b.coeff(1))?(a.coeff(1)< b.coeff(1)) : (a.coeff(0)<b.coeff(0));
00116         }
00117 
00118         inline static bool greater(const Derived1& a, const Derived2& b) {
00119                 return (a.coeff(1)!=b.coeff(1))?(a.coeff(1)> b.coeff(1)) : (a.coeff(0)>b.coeff(0));
00120         }
00121 
00122         inline static bool lessEqual(const Derived1& a, const Derived2& b) {
00123                 return (a.coeff(1)!=b.coeff(1))?(a.coeff(1)< b.coeff(1)) : (a.coeff(0)<=b.coeff(0));
00124         }
00125 
00126         inline static bool greaterEqual(const Derived1& a, const Derived2& b) {
00127                 return (a.coeff(1)!=b.coeff(1))?(a.coeff(1)> b.coeff(1)) : (a.coeff(0)>=b.coeff(0));
00128         }
00129 };
00130 
00131 template<typename Derived1, typename Derived2> struct ei_lexi_comparison<Derived1,Derived2,3>
00132 {
00133         inline static bool less(const Derived1& a, const Derived2& b) {
00134                 return  (a.coeff(2)!=b.coeff(2))?(a.coeff(2)< b.coeff(2)):
00135                                     (a.coeff(1)!=b.coeff(1))?(a.coeff(1)< b.coeff(1)) : (a.coeff(0)<b.coeff(0));
00136         }
00137 
00138         inline static bool greater(const Derived1& a, const Derived2& b) {
00139                 return  (a.coeff(2)!=b.coeff(2))?(a.coeff(2)> b.coeff(2)):
00140                                     (a.coeff(1)!=b.coeff(1))?(a.coeff(1)> b.coeff(1)) : (a.coeff(0)>b.coeff(0));
00141         }
00142 
00143         inline static bool lessEqual(const Derived1& a, const Derived2& b) {
00144                 return  (a.coeff(2)!=b.coeff(2))?(a.coeff(2)< b.coeff(2)):
00145                                     (a.coeff(1)!=b.coeff(1))?(a.coeff(1)< b.coeff(1)) : (a.coeff(0)<=b.coeff(0));
00146         }
00147 
00148         inline static bool greaterEqual(const Derived1& a, const Derived2& b) {
00149                 return  (a.coeff(2)!=b.coeff(2))?(a.coeff(2)> b.coeff(2)):
00150                                     (a.coeff(1)!=b.coeff(1))?(a.coeff(1)> b.coeff(1)) : (a.coeff(0)>=b.coeff(0));
00151         }
00152 };
00153 
00154 template<typename Derived1, typename Derived2> struct ei_lexi_comparison<Derived1,Derived2,4>
00155 {
00156         inline static bool less(const Derived1& a, const Derived2& b) {
00157                 return  (a.coeff(3)!=b.coeff(3))?(a.coeff(3)< b.coeff(3)) : (a.coeff(2)!=b.coeff(2))?(a.coeff(2)< b.coeff(2)):
00158                                     (a.coeff(1)!=b.coeff(1))?(a.coeff(1)< b.coeff(1)) : (a.coeff(0)<b.coeff(0));
00159         }
00160 
00161         inline static bool greater(const Derived1& a, const Derived2& b) {
00162                 return  (a.coeff(3)!=b.coeff(3))?(a.coeff(3)> b.coeff(3)) : (a.coeff(2)!=b.coeff(2))?(a.coeff(2)> b.coeff(2)):
00163                                     (a.coeff(1)!=b.coeff(1))?(a.coeff(1)> b.coeff(1)) : (a.coeff(0)>b.coeff(0));
00164         }
00165 
00166         inline static bool lessEqual(const Derived1& a, const Derived2& b) {
00167                 return  (a.coeff(3)!=b.coeff(3))?(a.coeff(3)< b.coeff(3)) : (a.coeff(2)!=b.coeff(2))?(a.coeff(2)< b.coeff(2)):
00168                                     (a.coeff(1)!=b.coeff(1))?(a.coeff(1)< b.coeff(1)) : (a.coeff(0)<=b.coeff(0));
00169         }
00170 
00171         inline static bool greaterEqual(const Derived1& a, const Derived2& b) {
00172                 return  (a.coeff(3)!=b.coeff(3))?(a.coeff(3)> b.coeff(3)) : (a.coeff(2)!=b.coeff(2))?(a.coeff(2)> b.coeff(2)):
00173                                     (a.coeff(1)!=b.coeff(1))?(a.coeff(1)> b.coeff(1)) : (a.coeff(0)>=b.coeff(0));
00174         }
00175 };
00176 
00177 // implementation of Import
00178 template<typename Derived1, typename Derived2>
00179 struct ei_import_selector<Derived1,Derived2,true,true>
00180 {
00181         static void run(Derived1& a, const Derived2& b) { a = b; }
00182 };
00183 
00184 template<typename Derived1, typename Derived2>
00185 struct ei_import_selector<Derived1,Derived2,false,true>
00186 {
00187         static void run(Derived1& a, const Derived2& b)
00188         { a = b.template cast<typename Derived1::Scalar>(); }
00189 };
00190 
00191 template<typename Derived1, typename Derived2>
00192 struct ei_import_selector<Derived1,Derived2,false,false>
00193 {
00194         static void run(Derived1& a, const Derived2& b)
00195         {
00196                 EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived1);
00197                 EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived1);
00198                 EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived2);
00199                 EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived2);
00200                 enum {
00201                         Size1 = Derived1::SizeAtCompileTime,
00202                         Size2 = Derived2::SizeAtCompileTime
00203                 };
00204                 assert(Size1<=4 && Size2<=4);
00205                 a.coeffRef(0) = Scalar(b.coeff(0));
00206                 if (Size1>1) { if (Size2>1) a.coeffRef(1) = Scalar(b.coeff(1)); else a.coeffRef(1) = 0; }
00207                 if (Size1>2) { if (Size2>2) a.coeffRef(2) = Scalar(b.coeff(2)); else a.coeffRef(2) = 0; }
00208                 if (Size1>3) { if (Size2>3) a.coeffRef(3) = Scalar(b.coeff(3)); else a.coeffRef(3) = 0; }
00209         }
00210 };
00211 
00212 // default implementation of ei_to_vcgtype
00213 // the specialization are with
00214 template<typename XprType,int Rows,int Cols,int StorageOrder,int MRows,int MCols>
00215 struct ei_to_vcgtype { typedef Matrix<typename XprType::Scalar,Rows,Cols,StorageOrder,MRows,MCols> type; };
00216 
00217 }
00218 
00219 #define VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, Op) \
00220         template<typename OtherDerived> \
00221         Derived& operator Op(const Eigen::MatrixBase<OtherDerived>& other) \
00222         { \
00223                 Base::operator Op(other.derived());  return *this;\
00224         } \
00225         Derived& operator Op(const Derived& other) \
00226         { \
00227                 Base::operator Op(other); return *this;\
00228         }
00229 
00230         #define VCG_EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, Op) \
00231         template<typename Other> \
00232         Derived& operator Op(const Other& scalar) \
00233         { \
00234                 Base::operator Op(scalar); return *this;\
00235         }
00236 
00237 #define VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Derived) \
00238         VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, =) \
00239         VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, +=) \
00240         VCG_EIGEN_INHERIT_ASSIGNMENT_OPERATOR(Derived, -=) \
00241         VCG_EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, *=) \
00242         VCG_EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, /=)
00243 
00244 
00245 namespace vcg {
00246 
00247 template<typename Derived1, typename Derived2>
00248 typename Eigen::ei_traits<Derived1>::Scalar
00249 Angle(const Eigen::MatrixBase<Derived1>& p1, const Eigen::MatrixBase<Derived2> & p2)
00250 {
00251         EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived1)
00252         EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived2)
00253         EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived1)
00254         EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived2)
00255         EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived1,Derived2)
00256         typedef typename Eigen::ei_traits<Derived1>::Scalar Scalar;
00257 
00258         Scalar w = p1.norm()*p2.norm();
00259         if(w==0) return Scalar(-1);
00260         Scalar t = (p1.dot(p2))/w;
00261         if(t>1) t = 1;
00262         else if(t<-1) t = -1;
00263         return vcg::math::Acos(t);
00264 }
00265 
00266 template<typename Derived1, typename Derived2>
00267 typename Eigen::ei_traits<Derived1>::Scalar
00268 AngleN(const Eigen::MatrixBase<Derived1>& p1, const Eigen::MatrixBase<Derived2> & p2)
00269 {
00270         EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived1)
00271         EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived2)
00272         EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived1)
00273         EIGEN_STATIC_ASSERT_FIXED_SIZE(Derived2)
00274         EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived1,Derived2)
00275         typedef typename Eigen::ei_traits<Derived1>::Scalar Scalar;
00276 
00277         Scalar t = (p1.dot(p2));
00278         if(t>1) t = 1;
00279         else if(t<-1) t = -1;
00280         return vcg::math::Acos(t);
00281 }
00282 
00283 template<typename Derived1>
00284 inline typename Eigen::ei_traits<Derived1>::Scalar Norm( const Eigen::MatrixBase<Derived1>& p)
00285 { return p.norm(); }
00286 
00287 template<typename Derived1>
00288 inline typename Eigen::ei_traits<Derived1>::Scalar SquaredNorm( const Eigen::MatrixBase<Derived1>& p)
00289 { return p.squaredNorm(); }
00290 
00291 template<typename Derived1, typename Derived2>
00292 inline typename Eigen::ei_traits<Derived1>::Scalar
00293 Distance(const Eigen::MatrixBase<Derived1>& p1, const Eigen::MatrixBase<Derived2> & p2)
00294 { return (p1-p2).norm(); }
00295 
00296 template<typename Derived1, typename Derived2>
00297 inline typename Eigen::ei_traits<Derived1>::Scalar
00298 SquaredDistance(const Eigen::MatrixBase<Derived1>& p1, const Eigen::MatrixBase<Derived2> & p2)
00299 { return (p1-p2).squaredNorm(); }
00300 
00301 template<typename Derived>
00302 inline const Eigen::CwiseUnaryOp<Eigen::ei_scalar_abs_op<typename Eigen::ei_traits<Derived>::Scalar>, Derived>
00303 Abs(const Eigen::MatrixBase<Derived>& p)
00304 { return p.cwise().abs(); }
00305 
00307 template<typename Scalar,int Size,int StorageOrder>
00308 EIGEN_DEPRECATED inline Eigen::Matrix<Scalar,Size,Size,StorageOrder>&
00309 Transpose(const Eigen::Matrix<Scalar,Size,Size,StorageOrder>& m)
00310 { return m.transposeInPlace(); return m; }
00311 
00312 template<typename Derived>
00313 inline const Eigen::CwiseBinaryOp<Eigen::ei_scalar_max_op<typename Eigen::ei_traits<Derived>::Scalar>,
00314                                                                                                                                         Derived,
00315                                                                                                                                         Eigen::NestByValue<typename Derived::ConstantReturnType> >
00316 LowClampToZero(const Eigen::MatrixBase<Derived>& p)
00317 { return p.cwise().max(Derived::Zero().nestByValue()); }
00318 
00319 }
00320 
00321 #endif