bone heat

Former-commit-id: c2d904ac4e45633fde179cc60c3d6a6b6d5981df

RELEASE_HISTORY.txt

@@ -1,3 +1,4 @@
+1.0.3  Bone heat method
 1.0.2  Bug fix in winding number code
 1.0.1  Bug fixes and more CGAL support
 1.0.0  Major beta release: many renames, tutorial, triangle wrapper, org. build

VERSION.txt

@@ -3,4 +3,4 @@
 # Anyone may increment Minor to indicate a small change.
 # Major indicates a large change or large number of changes (upload to website)
 # World indicates a substantial change or release
-1.0.2
+1.0.3

include/igl/embree/bone_heat.cpp

@@ -0,0 +1,116 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2013 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 "bone_heat.h"
+#include "EmbreeIntersector.h"
+#include "bone_visible.h"
+#include "../project_to_line_segment.h"
+#include "../cotmatrix.h"
+#include "../massmatrix.h"
+#include "../mat_min.h"
+#include <Eigen/Sparse>
+
+bool igl::bone_heat(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & C,
+  const Eigen::VectorXi & P,
+  const Eigen::MatrixXi & BE,
+  const Eigen::MatrixXi & CE,
+  Eigen::MatrixXd & W)
+{
+  using namespace std;
+  using namespace Eigen;
+  assert(CE.rows() == 0 && "Cage edges not supported.");
+  assert(C.cols() == V.cols() && "V and C should have same #cols");
+  assert(BE.cols() == 2 && "BE should have #cols=2");
+  assert(F.cols() == 3 && "F should contain triangles.");
+  assert(V.cols() == 3 && "V should contain 3D positions.");
+
+  const int n = V.rows();
+  const int np = P.rows();
+  const int nb = BE.rows();
+  const int m = np + nb;
+
+  // "double sided lighting"
+  MatrixXi FF;
+  FF.resize(F.rows()*2,F.cols());
+  FF << F, F.rowwise().reverse();
+  // Initialize intersector
+  EmbreeIntersector ei;
+  ei.init(V.cast<float>(),F.cast<int>());
+
+  typedef Matrix<bool,Dynamic,1> VectorXb;
+  typedef Matrix<bool,Dynamic,Dynamic> MatrixXb;
+  MatrixXb vis_mask(n,m);
+  // Distances
+  MatrixXd D(n,m);
+  // loop over points
+  for(int j = 0;j<np;j++)
+  {
+    const Vector3d p = C.row(P(j));
+    D.col(j) = (V.rowwise()-p.transpose()).rowwise().norm();
+    VectorXb vj;
+    bone_visible(V,F,ei,p,p,vj);
+    vis_mask.col(j) = vj;
+  }
+
+  // loop over bones
+  for(int j = 0;j<nb;j++)
+  {
+    const Vector3d s = C.row(BE(j,0));
+    const Vector3d d = C.row(BE(j,1));
+    VectorXd t,sqrD;
+    project_to_line_segment(V,s,d,t,sqrD);
+    D.col(np+j) = sqrD.array().sqrt();
+    VectorXb vj;
+    bone_visible(V,F,ei,s,d,vj);
+    vis_mask.col(np+j) = vj;
+  }
+
+  if(CE.rows() > 0)
+  {
+    cerr<<"Error: Cage edges are not supported. Ignored."<<endl;
+  }
+
+  MatrixXd PP = MatrixXd::Zero(n,m);
+  VectorXd min_D;
+  VectorXd Hdiag = VectorXd::Zero(n);
+  VectorXi J;
+  mat_min(D,2,min_D,J);
+  for(int i = 0;i<n;i++)
+  {
+    PP(i,J(i)) = 1;
+    if(vis_mask(i,J(i)))
+    {
+      double hii = pow(min_D(i),-2.); 
+      Hdiag(i) = (hii>1e10?1e10:hii);
+    }
+  }
+  SparseMatrix<double> Q,L,M;
+  cotmatrix(V,F,L);
+  massmatrix(V,F,MASSMATRIX_TYPE_DEFAULT,M);
+  const auto & H = Hdiag.asDiagonal();
+  Q = (-L+M*H);
+  SimplicialLLT <SparseMatrix<double > > llt;
+  llt.compute(Q);
+  switch(llt.info())
+  {
+    case Eigen::Success:
+      break;
+    case Eigen::NumericalIssue:
+      cerr<<"Error: Numerical issue."<<endl;
+      return false;
+    default:
+      cerr<<"Error: Other."<<endl;
+      return false;
+  }
+
+  const auto & rhs = M*H*PP;
+  W = llt.solve(rhs);
+  return true;
+}

include/igl/embree/bone_heat.h

@@ -0,0 +1,44 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2013 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_BONE_HEAT_H
+#define IGL_BONE_HEAT_H
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  // BONE_HEAT  Compute skinning weights W given a surface mesh (V,F) and an
+  // internal skeleton (C,BE) according to "Automatic Rigging" [Baran and
+  // Popovic 2007].
+  //
+  // Inputs:
+  //   V  #V by 3 list of mesh vertex positions
+  //   F  #F by 3 list of mesh corner indices into V
+  //   C  #C by 3 list of joint locations
+  //   P  #P list of point handle indices into C
+  //   BE  #BE by 2 list of bone edge indices into C
+  //   CE  #CE by 2 list of cage edge indices into **P**
+  // Outputs:
+  //   W  #V by #P+#BE matrix of weights.
+  // Returns true only on success.
+  //
+  IGL_INLINE bool bone_heat(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXd & C,
+    const Eigen::VectorXi & P,
+    const Eigen::MatrixXi & BE,
+    const Eigen::MatrixXi & CE,
+    Eigen::MatrixXd & W);
+};
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "bone_heat.cpp"
+#endif
+
+#endif

include/igl/embree/bone_visible.cpp

@@ -6,7 +6,6 @@
 // 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 "bone_visible.h"
-#include "EmbreeIntersector.h"
 #include <igl/project_to_line.h>
 #include <igl/EPS.h>
 #include <igl/Timer.h>
@@ -24,10 +23,6 @@ IGL_INLINE void igl::bone_visible(
   const Eigen::PlainObjectBase<DerivedSD> & d,
   Eigen::PlainObjectBase<Derivedflag>  & flag)
 {
-  using namespace igl;
-  using namespace std;
-  using namespace Eigen;
-  flag.resize(V.rows());
   // "double sided lighting"
   Eigen::PlainObjectBase<DerivedF> FF;
   FF.resize(F.rows()*2,F.cols());
@@ -35,13 +30,33 @@ IGL_INLINE void igl::bone_visible(
   // Initialize intersector
   EmbreeIntersector ei;
   ei.init(V.template cast<float>(),FF.template cast<int>());
+  return bone_visible(V,F,ei,s,d,flag);
+}
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedSD,
+  typename Derivedflag>
+IGL_INLINE void igl::bone_visible(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const igl::EmbreeIntersector & ei,
+  const Eigen::PlainObjectBase<DerivedSD> & s,
+  const Eigen::PlainObjectBase<DerivedSD> & d,
+  Eigen::PlainObjectBase<Derivedflag>  & flag)
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+  flag.resize(V.rows());
   const double sd_norm = (s-d).norm();
   // Embree seems to be parallel when constructing but not when tracing rays
 #pragma omp parallel for
   // loop over mesh vertices
   for(int v = 0;v<V.rows();v++)
   {
-    Vector3d Vv = V.row(v);
+    const Vector3d Vv = V.row(v);
     // Project vertex v onto line segment sd
     //embree.intersectSegment
     double t,sqrd;

include/igl/embree/bone_visible.h

@@ -9,29 +9,30 @@
 #define IGL_BONE_VISIBLE_H
 #include <igl/igl_inline.h>
 #include <Eigen/Core>
-//
-// BONE_VISIBLE  test whether vertices of mesh are "visible" to a given bone,
-// where "visible" is defined as in [Baran & Popovic 07]. Instead of checking
-// whether each point can see *any* of the bone, we just check if each point
-// can see its own projection onto the bone segment. In other words, we project
-// each vertex v onto the bone, projv. Then we check if there are any
-// intersections between the line segment (projv-->v) and the mesh.
-//
-// [flag] = bone_visible(V,F,s,d);
-//
-// Input:
-//    s  row vector of position of start end point of bone
-//    d  row vector of position of dest end point of bone
-//    V  #V by 3 list of vertex positions
-//    F  #F by 3 list of triangle indices
-// Output:
-//    flag  #V by 1 list of bools (true) visible, (false) obstructed
-//
-// Note: This checks for hits along the segment which are facing in *any*
-// direction from the ray.
-//
+#include "EmbreeIntersector.h"
 namespace igl
 {
+  //
+  // BONE_VISIBLE  test whether vertices of mesh are "visible" to a given bone,
+  // where "visible" is defined as in [Baran & Popovic 07]. Instead of checking
+  // whether each point can see *any* of the bone, we just check if each point
+  // can see its own projection onto the bone segment. In other words, we project
+  // each vertex v onto the bone, projv. Then we check if there are any
+  // intersections between the line segment (projv-->v) and the mesh.
+  //
+  // [flag] = bone_visible(V,F,s,d);
+  //
+  // Input:
+  //    V  #V by 3 list of vertex positions
+  //    F  #F by 3 list of triangle indices
+  //    s  row vector of position of start end point of bone
+  //    d  row vector of position of dest end point of bone
+  // Output:
+  //    flag  #V by 1 list of bools (true) visible, (false) obstructed
+  //
+  // Note: This checks for hits along the segment which are facing in *any*
+  // direction from the ray.
+  //
   template <
     typename DerivedV, 
     typename DerivedF, 
@@ -43,6 +44,20 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedSD> & s,
     const Eigen::PlainObjectBase<DerivedSD> & d,
     Eigen::PlainObjectBase<Derivedflag>  & flag);
+  // Inputs:
+  //  ei  EmbreeIntersector for mesh (V,F) should be double sided
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedSD,
+    typename Derivedflag>
+  IGL_INLINE void bone_visible(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const igl::EmbreeIntersector & ei,
+    const Eigen::PlainObjectBase<DerivedSD> & s,
+    const Eigen::PlainObjectBase<DerivedSD> & d,
+    Eigen::PlainObjectBase<Derivedflag>  & flag);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "bone_visible.cpp"

include/igl/min_quad_with_fixed.h

@@ -1,6 +1,6 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 //
-// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 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

include/igl/mosek/mosek_quadprog.cpp

@@ -16,11 +16,6 @@
 #include "../harwell_boeing.h"
 #include "../EPS.h"
 
-// Little function mosek needs in order to know how to print to std out
-static void MSKAPI printstr(void *handle, char str[])
-{
-  printf("%s",str);
-}
 
 igl::MosekData::MosekData()
 {
@@ -109,6 +104,11 @@ IGL_INLINE bool igl::mosek_quadprog(
   mosek_guarded(MSK_makeenv(&env,NULL,NULL,NULL,NULL));
 #endif
   ///* Directs the log stream to the 'printstr' function. */
+  //// Little function mosek needs in order to know how to print to std out
+  //const auto & printstr = [](void *handle, char str[])
+  //{
+  //  printf("%s",str);
+  //}
   //mosek_guarded(MSK_linkfunctoenvstream(env,MSK_STREAM_LOG,NULL,printstr));
   // initialize mosek environment
   mosek_guarded(MSK_initenv(env));

include/igl/project_to_line_segment.cpp

@@ -28,13 +28,14 @@ IGL_INLINE void igl::project_to_line_segment(
 #pragma omp parallel for if (np>10000)
   for(int p = 0;p<np;p++)
   {
+    const Eigen::PlainObjectBase<DerivedS> Pp = P.row(p);
     if(t(p)<0)
     {
-      sqrD(p) = (P.row(p)-S).squaredNorm();
+      sqrD(p) = (Pp-S).squaredNorm();
       t(p) = 0;
     }else if(t(p)>1)
     {
-      sqrD(p) = (P.row(p)-D).squaredNorm();
+      sqrD(p) = (Pp-D).squaredNorm();
       t(p) = 1;
     }
   }