generalized boundary_faces and fixed mat4_to_quat

Former-commit-id: 00ee86081b11a89935ff537778528c6844286ee5

Makefile.conf

@@ -27,9 +27,11 @@ ifeq ($(IGL_USERNAME),ajx)
 	IGL_WITH_MOSEK=1
 	IGL_WITH_PNG=1
 	# I don't use llvm
-	AFLAGS = -m64 -march="corei7-avx"
+#AFLAGS = -m64 -march="corei7-avx"
+	AFLAGS = -m64
 	OPENMP=-fopenmp
-	SSE=-mavx
+#SSE=-mavx
+	SSE=-msse
 endif
 
 ifeq ($(IGL_USERNAME),alecjaco) 

include/igl/boundary_faces.cpp

@@ -16,33 +16,57 @@ IGL_INLINE void igl::boundary_faces(
   using namespace std;
   using namespace igl;
 
+  if(T.size() == 0)
+  {
+    F.clear();
+    return;
+  }
+
+  int simplex_size = T[0].size();
   // Get a list of all faces
-  vector<vector<IntegerF> > allF(T.size()*4,vector<IntegerF>(3));
+  vector<vector<IntegerF> > allF(
+    T.size()*simplex_size,
+    vector<IntegerF>(simplex_size-1));
+
   // Gather faces, loop over tets
   for(int i = 0; i< (int)T.size();i++)
   {
-    assert(T[i].size() == 4);
-    // get face in correct order
-    allF[i*4+0][0] = T[i][1];
-    allF[i*4+0][1] = T[i][3];
-    allF[i*4+0][2] = T[i][2];
-    // get face in correct order
-    allF[i*4+1][0] = T[i][0];
-    allF[i*4+1][1] = T[i][2];
-    allF[i*4+1][2] = T[i][3];
-    // get face in correct order
-    allF[i*4+2][0] = T[i][0];
-    allF[i*4+2][1] = T[i][3];
-    allF[i*4+2][2] = T[i][1];
-    // get face in correct order
-    allF[i*4+3][0] = T[i][0];
-    allF[i*4+3][1] = T[i][1];
-    allF[i*4+3][2] = T[i][2];
+    assert(T[i].size() == simplex_size);
+    switch(simplex_size)
+    {
+      case 4:
+        // get face in correct order
+        allF[i*simplex_size+0][0] = T[i][1];
+        allF[i*simplex_size+0][1] = T[i][3];
+        allF[i*simplex_size+0][2] = T[i][2];
+        // get face in correct order
+        allF[i*simplex_size+1][0] = T[i][0];
+        allF[i*simplex_size+1][1] = T[i][2];
+        allF[i*simplex_size+1][2] = T[i][3];
+        // get face in correct order
+        allF[i*simplex_size+2][0] = T[i][0];
+        allF[i*simplex_size+2][1] = T[i][3];
+        allF[i*simplex_size+2][2] = T[i][1];
+        // get face in correct order
+        allF[i*simplex_size+3][0] = T[i][0];
+        allF[i*simplex_size+3][1] = T[i][1];
+        allF[i*simplex_size+3][2] = T[i][2];
+        break;
+      case 3:
+        allF[i*simplex_size+0][0] = T[i][0];
+        allF[i*simplex_size+0][1] = T[i][1];
+        allF[i*simplex_size+1][0] = T[i][1];
+        allF[i*simplex_size+1][1] = T[i][2];
+        allF[i*simplex_size+2][0] = T[i][2];
+        allF[i*simplex_size+2][1] = T[i][0];
+        break;
+    }
   }
   // Counts
   vector<int> C;
   face_occurences(allF,C);
 
+  // Q: Why not just count the number of ones?
   int twos = (int) count(C.begin(),C.end(),2);
   // Resize output to fit number of ones
   F.resize(allF.size() - twos);
@@ -68,7 +92,7 @@ IGL_INLINE void igl::boundary_faces(
   const Eigen::PlainObjectBase<DerivedT>& T,
   Eigen::PlainObjectBase<DerivedF>& F)
 {
-  assert(T.cols() == 0 || T.cols() == 4);
+  assert(T.cols() == 0 || T.cols() == 4 || T.cols() == 3);
   using namespace std;
   using namespace Eigen;
   using namespace igl;
@@ -79,6 +103,15 @@ IGL_INLINE void igl::boundary_faces(
   boundary_faces(vT,vF);
   list_to_matrix(vF,F);
 }
+
+template <typename DerivedT, typename DerivedF>
+IGL_INLINE Eigen::PlainObjectBase<DerivedF> igl::boundary_faces(
+  const Eigen::PlainObjectBase<DerivedT>& T)
+{
+  Eigen::Matrix<typename DerivedF::Scalar, Eigen::Dynamic, Eigen::Dynamic> F;
+  igl::boundary_faces(T,F);
+  return F;
+}
 #endif
 
 
@@ -86,5 +119,6 @@ IGL_INLINE void igl::boundary_faces(
 // Explicit template specialization
 template void igl::boundary_faces<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::boundary_faces<int, int>(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
+template Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > igl::boundary_faces(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
 #endif
 

include/igl/boundary_faces.h

@@ -10,15 +10,16 @@
 
 namespace igl
 {
-  // BOUNDARY_FACES Determine boundary faces of tetrahedra stored in T
+  // BOUNDARY_FACES Determine boundary faces (edges) of tetrahedra (triangles)
+  // stored in T
   //
   // Templates:
   //   IntegerT  integer-value: e.g. int
   //   IntegerF  integer-value: e.g. int
   // Input:
-  //  T  tetrahedron index list, m by 4, where m is the number of tetrahedra
+  //  T  tetrahedron (triangle) index list, m by 4 (3), where m is the number of tetrahedra
   // Output:
-  //  F  list of boundary faces, n by 3, where n is the number of boundary faces
+  //  F  list of boundary faces, n by 3 (2), where n is the number of boundary faces
   //
   //
   template <typename IntegerT, typename IntegerF>
@@ -34,6 +35,10 @@ namespace igl
   IGL_INLINE void boundary_faces(
     const Eigen::PlainObjectBase<DerivedT>& T,
     Eigen::PlainObjectBase<DerivedF>& F);
+  // Same as above but returns F
+  template <typename DerivedT, typename DerivedF>
+  IGL_INLINE Eigen::PlainObjectBase<DerivedF> boundary_faces(
+    const Eigen::PlainObjectBase<DerivedT>& T);
 #endif
 }
 

include/igl/mat_to_quat.cpp

@@ -8,20 +8,73 @@ static inline Q_type ReciprocalSqrt( const Q_type x )
   return 1.0/sqrt(x);
 }
 
-// Converts row major order matrix to quat
-// http://software.intel.com/sites/default/files/m/d/4/1/d/8/293748.pdf
+//// Converts row major order matrix to quat
+//// http://software.intel.com/sites/default/files/m/d/4/1/d/8/293748.pdf
+//template <typename Q_type>
+//IGL_INLINE void igl::mat4_to_quat(const Q_type * m, Q_type * q)
+//{
+//  Q_type t = + m[0 * 4 + 0] + m[1 * 4 + 1] + m[2 * 4 + 2] + 1.0f; 
+//  Q_type s = ReciprocalSqrt( t ) * 0.5f;
+//  q[3] = s * t;
+//  q[2] = ( m[0 * 4 + 1] - m[1 * 4 + 0] ) * s; 
+//  q[1] = ( m[2 * 4 + 0] - m[0 * 4 + 2] ) * s; 
+//  q[0] = ( m[1 * 4 + 2] - m[2 * 4 + 1] ) * s;
+//}
+
+// https://bmgame.googlecode.com/svn/idlib/math/Simd_AltiVec.cpp
 template <typename Q_type>
-IGL_INLINE void igl::mat4_to_quat(const Q_type * m, Q_type * q)
+IGL_INLINE void igl::mat4_to_quat(const Q_type * mat, Q_type * q)
 {
-  Q_type t = + m[0 * 4 + 0] + m[1 * 4 + 1] + m[2 * 4 + 2] + 1.0f; 
-  Q_type s = ReciprocalSqrt( t ) * 0.5f;
-  q[3] = s * t;
-  q[2] = ( m[0 * 4 + 1] - m[1 * 4 + 0] ) * s; 
-  q[1] = ( m[2 * 4 + 0] - m[0 * 4 + 2] ) * s; 
-  q[0] = ( m[1 * 4 + 2] - m[2 * 4 + 1] ) * s;
+  Q_type trace;
+  Q_type s;
+  Q_type t;
+  int i;
+  int j;
+  int k;
+  
+  static int next[3] = { 1, 2, 0 };
+
+  trace = mat[0 * 4 + 0] + mat[1 * 4 + 1] + mat[2 * 4 + 2];
+
+  if ( trace > 0.0f ) {
+
+    t = trace + 1.0f;
+    s = ReciprocalSqrt( t ) * 0.5f;
+
+    q[3] = s * t;
+    q[0] = ( mat[1 * 4 + 2] - mat[2 * 4 + 1] ) * s;
+    q[1] = ( mat[2 * 4 + 0] - mat[0 * 4 + 2] ) * s;
+    q[2] = ( mat[0 * 4 + 1] - mat[1 * 4 + 0] ) * s;
+
+  } else {
+
+    i = 0;
+    if ( mat[1 * 4 + 1] > mat[0 * 4 + 0] ) {
+      i = 1;
+    }
+    if ( mat[2 * 4 + 2] > mat[i * 4 + i] ) {
+      i = 2;
+    }
+    j = next[i];
+    k = next[j];
+
+    t = ( mat[i * 4 + i] - ( mat[j * 4 + j] + mat[k * 4 + k] ) ) + 1.0f;
+    s = ReciprocalSqrt( t ) * 0.5f;
+
+    q[i] = s * t;
+    q[3] = ( mat[j * 4 + k] - mat[k * 4 + j] ) * s;
+    q[j] = ( mat[i * 4 + j] + mat[j * 4 + i] ) * s;
+    q[k] = ( mat[i * 4 + k] + mat[k * 4 + i] ) * s;
+  }
+
+  //// Unused translation
+  //jq.t[0] = mat[0 * 4 + 3];
+  //jq.t[1] = mat[1 * 4 + 3];
+  //jq.t[2] = mat[2 * 4 + 3];
 }
 
 
+
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
 template void igl::mat4_to_quat<double>(double const*, double*);