better cross, volume and new barycentric coordinates (tets only)

Former-commit-id: 9d14fb6aecab2eac0d4d6266e5172a6e9e818b43

include/igl/barycentric_coordinates.cpp

@@ -0,0 +1,52 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "barycentric_coordinates.h"
+#include "volume.h"
+
+template <
+  typename DerivedP,
+  typename DerivedA,
+  typename DerivedB,
+  typename DerivedC,
+  typename DerivedD,
+  typename DerivedL>
+IGL_INLINE void igl::barycentric_coordinates(
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedA> & A,
+  const Eigen::PlainObjectBase<DerivedB> & B,
+  const Eigen::PlainObjectBase<DerivedC> & C,
+  const Eigen::PlainObjectBase<DerivedD> & D,
+  Eigen::PlainObjectBase<DerivedL> & L)
+{
+  using namespace Eigen;
+  assert(P.cols() == 3 && "query must be in 3d");
+  assert(A.cols() == 3 && "corners must be in 3d");
+  assert(B.cols() == 3 && "corners must be in 3d");
+  assert(C.cols() == 3 && "corners must be in 3d");
+  assert(D.cols() == 3 && "corners must be in 3d");
+  assert(P.rows() == A.rows() && "Must have same number of queries as corners");
+  assert(A.rows() == B.rows() && "Corners must be same size");
+  assert(A.rows() == C.rows() && "Corners must be same size");
+  assert(A.rows() == D.rows() && "Corners must be same size");
+  typedef Matrix<typename DerivedL::Scalar,DerivedL::RowsAtCompileTime,1> 
+    VectorXS;
+  // Total volume
+  VectorXS vol,LA,LB,LC,LD;
+  volume(B,D,C,P,LA);
+  volume(A,C,D,P,LB);
+  volume(A,D,B,P,LC);
+  volume(A,B,C,P,LD);
+  volume(A,B,C,D,vol);
+  L.resize(P.rows(),4);
+  L<<LA,LB,LC,LD;
+  L.array().colwise() /= vol.array();
+}
+
+#ifndef IGL_HEADER_ONLY
+template void igl::barycentric_coordinates<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/barycentric_coordinates.h

@@ -0,0 +1,45 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_BARYCENTRIC_COORDINATES_H
+#define IGL_BARYCENTRIC_COORDINATES_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Compute barycentric coordinates in a tet
+  //
+  // Inputs:
+  //   P  #P by 3 Query points in 3d
+  //   A  #P by 3 Tet corners in 3d
+  //   B  #P by 3 Tet corners in 3d
+  //   C  #P by 3 Tet corners in 3d
+  //   D  #P by 3 Tet corners in 3d
+  // Outputs:
+  //   L  #P by 4 list of barycentric coordinates
+  //   
+  template <
+    typename DerivedP,
+    typename DerivedA,
+    typename DerivedB,
+    typename DerivedC,
+    typename DerivedD,
+    typename DerivedL>
+  IGL_INLINE void barycentric_coordinates(
+    const Eigen::PlainObjectBase<DerivedP> & P,
+    const Eigen::PlainObjectBase<DerivedA> & A,
+    const Eigen::PlainObjectBase<DerivedB> & B,
+    const Eigen::PlainObjectBase<DerivedC> & C,
+    const Eigen::PlainObjectBase<DerivedD> & D,
+    Eigen::PlainObjectBase<DerivedL> & L);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "barycentric_coordinates.cpp"
+#endif
+
+#endif

include/igl/cross.cpp

@@ -17,3 +17,29 @@ IGL_INLINE void igl::cross(
   out[1] = a[2]*b[0]-a[0]*b[2];
   out[2] = a[0]*b[1]-a[1]*b[0];
 }
+
+template <
+  typename DerivedA,
+  typename DerivedB,
+  typename DerivedC>
+IGL_INLINE void igl::cross(
+  const Eigen::PlainObjectBase<DerivedA> & A,
+  const Eigen::PlainObjectBase<DerivedB> & B,
+  Eigen::PlainObjectBase<DerivedC> & C)
+{
+  assert(A.cols() == 3 && "#cols should be 3");
+  assert(B.cols() == 3 && "#cols should be 3");
+  assert(A.rows() == B.rows() && "#rows in A and B should be equal");
+  C.resize(A.rows(),3);
+  for(int d = 0;d<3;d++)
+  {
+    C.col(d) = 
+      A.col((d+1)%3).array() * B.col((d+2)%3).array() -
+      A.col((d+2)%3).array() * B.col((d+1)%3).array();
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+template void igl::cross<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::cross<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/cross.h

@@ -8,6 +8,7 @@
 #ifndef IGL_CROSS_H
 #define IGL_CROSS_H
 #include "igl_inline.h"
+#include <Eigen/Core>
 namespace igl
 {
   // Computes out = cross(a,b)
@@ -16,10 +17,22 @@ namespace igl
   //   b  right 3d vector
   // Outputs:
   //   out  result 3d vector
+  IGL_INLINE void cross( const double *a, const double *b, double *out);
+  // Computes C = cross(A,B,2);
+  //
+  // Inputs:
+  //   A  #A by 3 list of row-vectors
+  //   B  #A by 3 list of row-vectors
+  // Outputs:
+  //   C  #A by 3 list of row-vectors
+  template <
+    typename DerivedA,
+    typename DerivedB,
+    typename DerivedC>
   IGL_INLINE void cross(
-    const double *a, 
-    const double *b,
-    double *out);
+    const Eigen::PlainObjectBase<DerivedA> & A,
+    const Eigen::PlainObjectBase<DerivedB> & B,
+    Eigen::PlainObjectBase<DerivedC> & C);
 }
 
 #ifdef IGL_HEADER_ONLY

include/igl/doublearea.cpp

@@ -60,6 +60,20 @@ IGL_INLINE void igl::doublearea(
   }
 }
 
+template <
+  typename DerivedA,
+  typename DerivedB,
+  typename DerivedC>
+IGL_INLINE typename DerivedA::Scalar igl::doublearea_single(
+  const Eigen::PlainObjectBase<DerivedA> & A, 
+  const Eigen::PlainObjectBase<DerivedB> & B, 
+  const Eigen::PlainObjectBase<DerivedC> & C)
+{
+  auto r = A-C;
+  auto s = B-C;
+  return r(0)*s(1) - r(1)*s(0);
+}
+
 template <typename Derivedl, typename DeriveddblA>
 IGL_INLINE void igl::doublearea( 
   const Eigen::PlainObjectBase<Derivedl> & ul,
@@ -108,4 +122,5 @@ template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::M
 template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::doublearea<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template Eigen::Matrix<double, 2, 1, 0, 2, 1>::Scalar igl::doublearea_single<Eigen::Matrix<double, 2, 1, 0, 2, 1>, Eigen::Matrix<double, 2, 1, 0, 2, 1>, Eigen::Matrix<double, 2, 1, 0, 2, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > const&);
 #endif

include/igl/doublearea.h

@@ -24,14 +24,24 @@ namespace igl
   //   V  #V by dim list of mesh vertex positions
   //   F  #F by simplex_size list of mesh faces (must be triangles)
   // Outputs:
-  //   dblA  #F list of triangle double areas
+  //   dblA  #F list of triangle double areas (SIGNED only for 2D input)
   //
-  // Note: THESE ARE *NOT* SIGNED. In matlab doublearea is signed in 2d
   template <typename DerivedV, typename DerivedF, typename DeriveddblA>
   IGL_INLINE void doublearea( 
     const Eigen::PlainObjectBase<DerivedV> & V, 
     const Eigen::PlainObjectBase<DerivedF> & F, 
     Eigen::PlainObjectBase<DeriveddblA> & dblA);
+  // Single triangle in 2D!
+  //
+  // This should handle streams of corners not just single corners
+  template <
+    typename DerivedA,
+    typename DerivedB,
+    typename DerivedC>
+  IGL_INLINE typename DerivedA::Scalar doublearea_single( 
+    const Eigen::PlainObjectBase<DerivedA> & A, 
+    const Eigen::PlainObjectBase<DerivedB> & B, 
+    const Eigen::PlainObjectBase<DerivedC> & C);
   // Same as above but use instrinsic edge lengths rather than (V,F) mesh
   // Templates:
   //   DerivedV  derived type of eigen matrix for V (e.g. derived from

include/igl/volume.cpp

@@ -6,13 +6,15 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "volume.h"
+#include "cross.h"
+#include <Eigen/Geometry>
 template <
   typename DerivedV, 
   typename DerivedT, 
   typename Derivedvol>
 IGL_INLINE void igl::volume(
-  Eigen::PlainObjectBase<DerivedV>& V,
-  Eigen::PlainObjectBase<DerivedT>& T,
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedT>& T,
   Eigen::PlainObjectBase<Derivedvol>& vol)
 {
   using namespace Eigen;
@@ -28,11 +30,48 @@ IGL_INLINE void igl::volume(
   }
 }
 
+template <
+  typename DerivedA,
+  typename DerivedB,
+  typename DerivedC,
+  typename DerivedD,
+  typename Derivedvol>
+IGL_INLINE void igl::volume(
+  const Eigen::PlainObjectBase<DerivedA> & A,
+  const Eigen::PlainObjectBase<DerivedB> & B,
+  const Eigen::PlainObjectBase<DerivedC> & C,
+  const Eigen::PlainObjectBase<DerivedD> & D,
+  Eigen::PlainObjectBase<Derivedvol> & vol)
+{
+  const auto & AmD = A-D;
+  const auto & BmD = B-D;
+  const auto & CmD = C-D;
+  Eigen::PlainObjectBase<DerivedA> BmDxCmD;
+  cross(BmD.eval(),CmD.eval(),BmDxCmD);
+  const auto & AmDdx = (AmD.array() * BmDxCmD.array()).rowwise().sum();
+  vol = -AmDdx/6.;
+}
+
+template <
+  typename VecA,
+  typename VecB,
+  typename VecC,
+  typename VecD>
+IGL_INLINE typename VecA::Scalar igl::volume_single(
+  const VecA & a,
+  const VecB & b,
+  const VecC & c,
+  const VecD & d)
+{
+  return -(a-d).dot((b-d).cross(c-d))/6.;
+}
+
+
 template <
   typename DerivedL, 
   typename Derivedvol>
 IGL_INLINE void igl::volume(
-  Eigen::PlainObjectBase<DerivedL>& L,
+  const Eigen::PlainObjectBase<DerivedL>& L,
   Eigen::PlainObjectBase<Derivedvol>& vol)
 {
   using namespace Eigen;
@@ -68,5 +107,8 @@ IGL_INLINE void igl::volume(
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
 // generated by autoexplicit.sh
-template void igl::volume<Eigen::Matrix<double, -1, 6, 0, -1, 6>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::volume<Eigen::Matrix<double, -1, 6, 0, -1, 6>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar igl::volume_single<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3> const&);
+template Eigen::Matrix<double, 3, 1, 0, 3, 1>::Scalar igl::volume_single<Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1> >(Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 1> const&);
+template void igl::volume<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/volume.h

@@ -27,9 +27,32 @@ namespace igl
     typename DerivedT, 
     typename Derivedvol>
   IGL_INLINE void volume(
-    Eigen::PlainObjectBase<DerivedV>& V,
-    Eigen::PlainObjectBase<DerivedT>& T,
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedT>& T,
     Eigen::PlainObjectBase<Derivedvol>& vol);
+  template <
+    typename DerivedA,
+    typename DerivedB,
+    typename DerivedC,
+    typename DerivedD,
+    typename Derivedvol>
+  IGL_INLINE void volume(
+    const Eigen::PlainObjectBase<DerivedA> & A,
+    const Eigen::PlainObjectBase<DerivedB> & B,
+    const Eigen::PlainObjectBase<DerivedC> & C,
+    const Eigen::PlainObjectBase<DerivedD> & D,
+    Eigen::PlainObjectBase<Derivedvol> & vol);
+  // Single tet
+  template <
+    typename VecA,
+    typename VecB,
+    typename VecC,
+    typename VecD>
+  IGL_INLINE typename VecA::Scalar volume_single(
+    const VecA & a,
+    const VecB & b,
+    const VecC & c,
+    const VecD & d);
   // Intrinsic version:
   //
   // Inputs:
@@ -38,7 +61,7 @@ namespace igl
     typename DerivedL, 
     typename Derivedvol>
   IGL_INLINE void volume(
-    Eigen::PlainObjectBase<DerivedL>& L,
+    const Eigen::PlainObjectBase<DerivedL>& L,
     Eigen::PlainObjectBase<Derivedvol>& vol);
 }