Merge branch 'master' of https://github.com/libigl/libigl

Former-commit-id: 14414567054b102302f3ffa48656b7eac447b36c

include/igl/active_set.cpp

@@ -224,6 +224,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
     // Gather active constraints and resp. rhss
     PlainObjectBase<DerivedBeq> Beq_i;
     Beq_i.resize(Beq.rows()+as_ieq_count,1);
+    Beq_i.head(Beq.rows()) = Beq;
     {
       int k =0;
       for(int a=0;a<as_ieq.size();a++)
@@ -282,6 +283,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
       ret = SOLVER_STATUS_ERROR;
       break;
     }
+    //cout<<matlab_format((Aeq*Z-Beq).eval(),"cr")<<endl;
     //cout<<matlab_format(Z,"Z")<<endl;
 #ifdef ACTIVE_SET_CPP_DEBUG
     cout<<"  post"<<endl;

include/igl/colon.cpp

@@ -47,7 +47,7 @@ IGL_INLINE void igl::colon(
   I.resize(n);
   int i = 0;
   T v = (T)low;
-  while((low<hi && (H)v<=hi) || (low>hi && (H)v>=hi))
+  while((low==hi && (H)v==hi) || (low<hi && (H)v<=hi) || (low>hi && (H)v>=hi))
   {
     I(i) = v;
     v = v + (T)step;

include/igl/dot_row.cpp

@@ -7,23 +7,16 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "igl/dot_row.h"
 
-using namespace Eigen;
+template <typename DerivedV>
-
+IGL_INLINE Eigen::PlainObjectBase<DerivedV> igl::dot_row(
-namespace igl
+  const Eigen::PlainObjectBase<DerivedV>& A,
+  const Eigen::PlainObjectBase<DerivedV>& B
+  )
 {
+  assert(A.rows() == B.rows());
+  assert(A.cols() == B.cols());
 
-  template <typename DerivedV>
+  return (A.array() * B.array()).rowwise().sum();
-  IGL_INLINE Eigen::PlainObjectBase<DerivedV> dot_row(
-    const Eigen::PlainObjectBase<DerivedV>& A,
-    const Eigen::PlainObjectBase<DerivedV>& B
-    )
-  {
-    assert(A.rows() == B.rows());
-    assert(A.cols() == B.cols());
-
-    return (A.array() * B.array()).rowwise().sum();
-  }
-
 }
 
 #ifndef IGL_HEADER_ONLY

include/igl/dot_row.h

@@ -16,19 +16,16 @@ namespace igl
   // Compute the dot product between each row of A and B
   // Templates:
   //   DerivedV derived from vertex positions matrix type: i.e. MatrixXd
-  //   DerivedF derived from face indices matrix type: i.e. MatrixXi
   // Inputs:
   //   A  eigen matrix r by c
   //   B  eigen matrix r by c
-  // Outputs:
+  // Returns:
   //   d a column vector with r entries that contains the dot product of each corresponding row of A and B
   //
-  // See also: adjacency_matrix
   template <typename DerivedV>
   IGL_INLINE Eigen::PlainObjectBase<DerivedV> dot_row(
     const Eigen::PlainObjectBase<DerivedV>& A,
-    const Eigen::PlainObjectBase<DerivedV>& B
+    const Eigen::PlainObjectBase<DerivedV>& B);
-    );
 
 }
 

include/igl/find.cpp

@@ -66,4 +66,5 @@ IGL_INLINE void igl::find(
 template void igl::find<double, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::find<double, 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::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh
+template void igl::find<double, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/frame_field_deformer.cpp

@@ -284,7 +284,7 @@ IGL_INLINE void Frame_field_deformer::compute_optimal_positions()
   for (int i=0;i<nfree;i++)
   {
     b.row(i) << 0.0, 0.0, 0.0;
-    for (int k=0;k<VT[i].size();k++)					// for all incident triangles
+    for (int k=0;k<(int)VT[i].size();k++)					// for all incident triangles
     {
       t = VT[i][k];												// incident tri
 			vi = (i==F(t,0))?0:(i==F(t,1))?1:(i==F(t,2))?2:3;	// index of i in t
@@ -345,7 +345,7 @@ IGL_INLINE void Frame_field_deformer::compute_optimal_positions()
   using namespace Eigen;
 
   WW.resize(F.rows());
-	for (int i=0;i<FF.size();i++)
+	for (int i=0;i<(int)FF.size();i++)
 	{
 		Vector3d v0,v1,v2;
 		v0 = FF[i].col(0);
@@ -360,7 +360,7 @@ IGL_INLINE void Frame_field_deformer::compute_optimal_positions()
 
 		// polar decomposition to discard rotational component (unnecessary but makes it easier)
 		Eigen::JacobiSVD<Matrix<double,3,3> > svd(AI, Eigen::ComputeFullU | Eigen::ComputeFullV );
-		Matrix<double,3,3>  au = svd.matrixU();
+		//Matrix<double,3,3>  au = svd.matrixU();
 		Matrix<double,3,3>  av = svd.matrixV();
 		DiagonalMatrix<double,3>	as(svd.singularValues());
 		WW[i] = av*as*av.transpose();

include/igl/intersect.cpp

@@ -10,6 +10,7 @@ template <class M>
 IGL_INLINE void igl::intersect(const M & A, const M & B, M & C)
 {
   // Stupid O(size(A) * size(B)) to do it
+  // Alec: This should be implemented by using unique and sort like `setdiff`
   M dyn_C(A.size() > B.size() ? A.size() : B.size(),1);
   // count of intersects
   int c = 0;

include/igl/is_border_vertex.h

@@ -20,7 +20,7 @@ namespace igl
   // Inputs:
   //   V  #V by dim list of vertex positions 
   //   F  #F by 3 list of triangle indices
-  // Returns vector of indices of vertices on open boundary of F
+  // Returns #V vector of bools revealing whether vertices are on boundary
   //
   // Known Bugs: does not depend on V
   // 

include/igl/is_symmetric.cpp

@@ -50,6 +50,7 @@ IGL_INLINE bool igl::is_symmetric(
 {
   using namespace Eigen;
   using namespace igl;
+  using namespace std;
   if(A.rows() != A.cols())
   {
     return false;
@@ -60,6 +61,10 @@ IGL_INLINE bool igl::is_symmetric(
   VectorXi AmATI,AmATJ;
   Matrix<AType,Dynamic,1> AmATV;
   find(AmAT,AmATI,AmATJ,AmATV);
+  if(AmATI.size() == 0)
+  {
+    return true;
+  }
   
   return AmATV.maxCoeff() < epsilon && AmATV.minCoeff() > -epsilon;
 }

include/igl/jet.cpp

@@ -94,7 +94,7 @@ IGL_INLINE void igl::jet(const T x_in, T & r, T & g, T & b)
     b = 0;
   }else
   {
-    r = (bone-(x-7./8.)/(1.-7./8.)*0.5);
+    r = (rone-(x-7./8.)/(1.-7./8.)*0.5);
     g = 0;
     b = 0;
   }
@@ -106,9 +106,19 @@ IGL_INLINE void igl::jet(
   const bool normalize,
   Eigen::PlainObjectBase<DerivedC> & C)
 {
-  C.resize(Z.rows(),3);
   const double min_z = (normalize?Z.minCoeff():0);
   const double max_z = (normalize?Z.maxCoeff():-1);
+  return jet(Z,min_z,max_z,C);
+}
+
+template <typename DerivedZ, typename DerivedC>
+IGL_INLINE void igl::jet(
+  const Eigen::PlainObjectBase<DerivedZ> & Z,
+  const double min_z,
+  const double max_z,
+  Eigen::PlainObjectBase<DerivedC> & C)
+{
+  C.resize(Z.rows(),3);
   for(int r = 0;r<Z.rows();r++)
   {
     jet((-min_z+Z(r,0))/(max_z-min_z),C(r,0),C(r,1),C(r,2));

include/igl/jet.h

@@ -46,6 +46,15 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedZ> & Z,
     const bool normalize,
     Eigen::PlainObjectBase<DerivedC> & C);
+  // Inputs:
+  //   min_z  value at blue
+  //   max_z  value at red
+  template <typename DerivedZ, typename DerivedC>
+  IGL_INLINE void jet(
+    const Eigen::PlainObjectBase<DerivedZ> & Z,
+    const double min_Z,
+    const double max_Z,
+    Eigen::PlainObjectBase<DerivedC> & C);
 };
 
 #ifdef IGL_HEADER_ONLY

include/igl/list_to_matrix.cpp

@@ -121,4 +121,5 @@ template bool igl::list_to_matrix<int, Eigen::PlainObjectBase<Eigen::Matrix<int,
 template bool igl::list_to_matrix<double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template bool igl::list_to_matrix<int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
 template bool igl::list_to_matrix<double, Eigen::Matrix<double, 1, -1, 1, 1, -1> >(std::vector<double, std::allocator<double> > const&, Eigen::Matrix<double, 1, -1, 1, 1, -1>&);
+template bool igl::list_to_matrix<unsigned long, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > >(std::vector<unsigned long, std::allocator<unsigned long> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/matrix_to_list.cpp

@@ -35,9 +35,12 @@ IGL_INLINE void igl::matrix_to_list(
   using namespace std;
   V.resize(M.size());
   // loop over rows
-  for(int i = 0;i<M.size();i++)
+  for(int j = 0;j<M.cols();j++)
   {
-    V[i] = M[i];
+    for(int i = 0;i<M.rows();i++)
+    {
+      V[i+j*M.rows()] = M(i,j);
+    }
   }
 }
 
@@ -61,4 +64,6 @@ template std::vector<Eigen::Matrix<int, -1, 1, 0, -1, 1>::Scalar, std::allocator
 template void igl::matrix_to_list<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, std::allocator<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar> >, std::allocator<std::vector<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, std::allocator<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar> > > >&);
 template void igl::matrix_to_list<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, std::allocator<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar> >, std::allocator<std::vector<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, std::allocator<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar> > > >&);
 template void igl::matrix_to_list<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, std::vector<Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar, std::allocator<Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar> >&);
+template void igl::matrix_to_list<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, std::allocator<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar> >&);
+template void igl::matrix_to_list<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, std::vector<Eigen::Matrix<int, -1, 1, 0, -1, 1>::Scalar, std::allocator<Eigen::Matrix<int, -1, 1, 0, -1, 1>::Scalar> >&);
 #endif

include/igl/min_quad_with_fixed.cpp

@@ -12,7 +12,7 @@
 #include "find.h"
 #include "sparse.h"
 #include "repmat.h"
-#include "lu_lagrange.h"
+//#include "lu_lagrange.h"
 #include "full.h"
 #include "matlab_format.h"
 #include "EPS.h"
@@ -109,7 +109,10 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
   }else
   {
     // determine if A(unknown,unknown) is symmetric and/or positive definite
-    data.Auu_sym = is_symmetric(Auu,EPS<double>());
+    VectorXi AuuI,AuuJ;
+    MatrixXd AuuV;
+    find(Auu,AuuI,AuuJ,AuuV);
+    data.Auu_sym = is_symmetric(Auu,EPS<double>()*AuuV.maxCoeff());
   }
 
   // Determine number of linearly independent constraints
@@ -144,6 +147,7 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
     nc = data.AeqTQR.rank();
     assert(nc<=neq);
     data.Aeq_li = nc == neq;
+    //cout<<"data.Aeq_li: "<<data.Aeq_li<<endl;
   }else
   {
     data.Aeq_li = true;
@@ -378,6 +382,8 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
   using namespace std;
   using namespace Eigen;
   using namespace igl;
+  typedef Matrix<T,Dynamic,1> VectorXT;
+  typedef Matrix<T,Dynamic,Dynamic> MatrixXT;
   // number of known rows
   int kr = data.known.size();
   if(kr!=0)
@@ -402,56 +408,55 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
 
   if(data.Aeq_li)
   {
+    // number of lagrange multipliers aka linear equality constraints
+    int neq = data.lagrange.size();
+    // append lagrange multiplier rhs's
+    VectorXT BBeq(B.size() + Beq.size());
+    BBeq << B, (Beq*-2.0);
+    // Build right hand side
+    VectorXT BBequl;
+    igl::slice(BBeq,data.unknown_lagrange,BBequl);
+    MatrixXT BBequlcols;
+    igl::repmat(BBequl,1,cols,BBequlcols);
+    MatrixXT NB;
+    if(kr == 0)
+    {
+      NB = BBequlcols;
+    }else
+    {
+      NB = data.preY * Y + BBequlcols;
+    }
 
-  // number of lagrange multipliers aka linear equality constraints
+    //std::cout<<"NB=["<<std::endl<<NB<<std::endl<<"];"<<std::endl;
-  int neq = data.lagrange.size();
+    //cout<<matlab_format(NB,"NB")<<endl;
-  // append lagrange multiplier rhs's
+    switch(data.solver_type)
-  Eigen::Matrix<T,Eigen::Dynamic,1> BBeq(B.size() + Beq.size());
-  BBeq << B, (Beq*-2.0);
-  // Build right hand side
-  Eigen::Matrix<T,Eigen::Dynamic,1> BBequl;
-  igl::slice(BBeq,data.unknown_lagrange,BBequl);
-  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> BBequlcols;
-  igl::repmat(BBequl,1,cols,BBequlcols);
-  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> NB;
-  if(kr == 0)
-  {
-    NB = BBequlcols;
-  }else
-  {
-    NB = data.preY * Y + BBequlcols;
-  }
-
-  //std::cout<<"NB=["<<std::endl<<NB<<std::endl<<"];"<<std::endl;
-  //cout<<matlab_format(NB,"NB")<<endl;
-  switch(data.solver_type)
-  {
-    case igl::min_quad_with_fixed_data<T>::LLT:
-      sol = data.llt.solve(NB);
-      break;
-    case igl::min_quad_with_fixed_data<T>::LDLT:
-      sol = data.ldlt.solve(NB);
-      break;
-    case igl::min_quad_with_fixed_data<T>::LU:
-      // Not a bottleneck
-      sol = data.lu.solve(NB);
-      break;
-    default:
-      cerr<<"Error: invalid solver type"<<endl;
-      return false;
-  }
-  //std::cout<<"sol=["<<std::endl<<sol<<std::endl<<"];"<<std::endl;
-  // Now sol contains sol/-0.5
-  sol *= -0.5;
-  // Now sol contains solution
-  // Place solution in Z
-  for(int i = 0;i<(sol.rows()-neq);i++)
-  {
-    for(int j = 0;j<sol.cols();j++)
     {
-      Z(data.unknown_lagrange(i),j) = sol(i,j);
+      case igl::min_quad_with_fixed_data<T>::LLT:
+        sol = data.llt.solve(NB);
+        break;
+      case igl::min_quad_with_fixed_data<T>::LDLT:
+        sol = data.ldlt.solve(NB);
+        break;
+      case igl::min_quad_with_fixed_data<T>::LU:
+        // Not a bottleneck
+        sol = data.lu.solve(NB);
+        break;
+      default:
+        cerr<<"Error: invalid solver type"<<endl;
+        return false;
+    }
+    //std::cout<<"sol=["<<std::endl<<sol<<std::endl<<"];"<<std::endl;
+    // Now sol contains sol/-0.5
+    sol *= -0.5;
+    // Now sol contains solution
+    // Place solution in Z
+    for(int i = 0;i<(sol.rows()-neq);i++)
+    {
+      for(int j = 0;j<sol.cols();j++)
+      {
+        Z(data.unknown_lagrange(i),j) = sol(i,j);
+      }
     }
-  }
   }else
   {
     assert(data.solver_type == min_quad_with_fixed_data<T>::QR_LLT);

include/igl/setdiff.cpp

@@ -0,0 +1,83 @@
+// 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 "setdiff.h"
+#include "colon.h"
+#include "list_to_matrix.h"
+#include "sort.h"
+#include "unique.h"
+
+template <
+  typename DerivedA,
+  typename DerivedB,
+  typename DerivedC,
+  typename DerivedIA>
+IGL_INLINE void igl::setdiff(
+  const Eigen::PlainObjectBase<DerivedA> & A,
+  const Eigen::PlainObjectBase<DerivedB> & B,
+  Eigen::PlainObjectBase<DerivedC> & C,
+  Eigen::PlainObjectBase<DerivedIA> & IA)
+{
+  using namespace Eigen;
+  using namespace std;
+  // boring base cases
+  if(A.size() == 0)
+  {
+    C.resize(0,1);
+    IA.resize(0,1);
+  }
+  if(B.size() == 0)
+  {
+    C.resize(A.size(),1);
+    copy(A.data(),A.data()+A.size(),C.data());
+    IA = igl::colon<typename DerivedIA::Scalar>(0,C.size()-1);
+  }
+
+  // Get rid of any duplicates
+  typedef Matrix<typename DerivedA::Scalar,Dynamic,1> VectorA;
+  typedef Matrix<typename DerivedB::Scalar,Dynamic,1> VectorB;
+  VectorA uA;
+  VectorB uB;
+  typedef PlainObjectBase<DerivedIA> IAType;
+  IAType uIA,uIuA,uIB,uIuB;
+  unique(A,uA,uIA,uIuA);
+  unique(B,uB,uIB,uIuB);
+
+  // Sort both
+  VectorA sA;
+  VectorB sB;
+  IAType sIA,sIB;
+  sort(uA,1,true,sA,sIA);
+  sort(uB,1,true,sB,sIB);
+
+  vector<typename DerivedB::Scalar> vC;
+  vector<typename DerivedIA::Scalar> vIA;
+  int bi = 0;
+  // loop over sA
+  bool past = false;
+  for(int a = 0;a<sA.size();a++)
+  {
+    while(!past && sA(a)>sB(bi))
+    {
+      bi++;
+      past = bi>=sB.size();
+    }
+    if(past || sA(a)<sB(bi))
+    {
+      // then sA(a) did not appear in sB
+      vC.push_back(sA(a));
+      vIA.push_back(uIA(sIA(a)));
+    }
+  }
+  list_to_matrix(vC,C);
+  list_to_matrix(vIA,IA);
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template void igl::setdiff<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -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<int, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/setdiff.h

@@ -0,0 +1,38 @@
+// 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_SETDIFF_H
+#define IGL_SETDIFF_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Set difference of elements of matrices
+  //
+  // Inputs:
+  //   A  m-long vector of indices
+  //   B  n-long vector of indices
+  // Outputs:
+  //   C  (k<=m)-long vector of unique elements appearing in A but not in B
+  //   IA  (k<=m)-long list of indices into A so that C = A(IA)
+  //
+  template <
+    typename DerivedA,
+    typename DerivedB,
+    typename DerivedC,
+    typename DerivedIA>
+  IGL_INLINE void setdiff(
+    const Eigen::PlainObjectBase<DerivedA> & A,
+    const Eigen::PlainObjectBase<DerivedB> & B,
+    Eigen::PlainObjectBase<DerivedC> & C,
+    Eigen::PlainObjectBase<DerivedIA> & IA);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "setdiff.cpp"
+#endif
+#endif

include/igl/sort.cpp

@@ -232,4 +232,5 @@ template void igl::sort<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<in
 template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
 template void igl::sort<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::sort<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&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/unique.cpp

@@ -11,6 +11,7 @@
 #include "SortableRow.h"
 #include "sortrows.h"
 #include "list_to_matrix.h"
+#include "matrix_to_list.h"
 
 #include <algorithm>
 #include <iostream>
@@ -61,6 +62,29 @@ IGL_INLINE void igl::unique(
 
 }
 
+template <
+  typename DerivedA,
+  typename DerivedC,
+  typename DerivedIA,
+  typename DerivedIC>
+IGL_INLINE void igl::unique(
+    const Eigen::PlainObjectBase<DerivedA> & A,
+    Eigen::PlainObjectBase<DerivedC> & C,
+    Eigen::PlainObjectBase<DerivedIA> & IA,
+    Eigen::PlainObjectBase<DerivedIC> & IC)
+{
+  using namespace std;
+  using namespace Eigen;
+  vector<typename DerivedA::Scalar > vA;
+  vector<typename DerivedC::Scalar > vC;
+  vector<size_t> vIA,vIC;
+  matrix_to_list(A,vA);
+  unique(vA,vC,vIA,vIC);
+  list_to_matrix(vC,C);
+  list_to_matrix(vIA,IA);
+  list_to_matrix(vIC,IC);
+}
+
 // Obsolete slow version converting to vectors
 // template <typename DerivedA, typename DerivedIA, typename DerivedIC>
 // IGL_INLINE void igl::unique_rows(
@@ -196,10 +220,13 @@ IGL_INLINE void igl::unique_rows(
 }
 
 #ifndef IGL_HEADER_ONLY
+// Explicit template specialization
 template void igl::unique<int>(std::vector<int, std::allocator<int> > const&, std::vector<int, std::allocator<int> >&, std::vector<size_t, std::allocator<size_t> >&, std::vector<size_t, std::allocator<size_t> >&);
 template void igl::unique_rows<Eigen::Matrix<int, -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<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::unique_rows<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<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::unique_rows<Eigen::Matrix<int, -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<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::unique_rows<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<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::unique_rows<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::unique<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::unique<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -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<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/unique.h

@@ -29,7 +29,16 @@ namespace igl
     std::vector<T> & C,
     std::vector<size_t> & IA,
     std::vector<size_t> & IC);
-
+  template <
+    typename DerivedA,
+    typename DerivedC,
+    typename DerivedIA,
+    typename DerivedIC>
+  IGL_INLINE void unique(
+      const Eigen::PlainObjectBase<DerivedA> & A,
+      Eigen::PlainObjectBase<DerivedC> & C,
+      Eigen::PlainObjectBase<DerivedIA> & IA,
+      Eigen::PlainObjectBase<DerivedIC> & IC);
   // Act like matlab's [C,IA,IC] = unique(X,'rows')
   //
   // Templates:

tutorial/303_LaplaceEquation/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(303_LaplaceEquation)
+
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${CMAKE_PROJECT_NAME} ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${CMAKE_PROJECT_NAME} ${SHARED_LIBRARIES})

tutorial/303_LaplaceEquation/main.cpp

@@ -0,0 +1,71 @@
+#include <igl/boundary_faces.h>
+#include <igl/colon.h>
+#include <igl/cotmatrix.h>
+#include <igl/jet.h>
+#include <igl/min_quad_with_fixed.h>
+#include <igl/readOFF.h>
+#include <igl/setdiff.h>
+#include <igl/slice.h>
+#include <igl/slice_into.h>
+#include <igl/unique.h>
+#include <igl/viewer/Viewer.h>
+#include <Eigen/Sparse>
+#include <iostream>
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+  MatrixXd V;
+  MatrixXi F;
+  igl::readOFF("../shared/camelhead.off",V,F);
+  // Find boundary edges
+  MatrixXi E;
+  igl::boundary_faces(F,E);
+  // Find boundary vertices
+  VectorXi b,IA,IC;
+  igl::unique(E,b,IA,IC);
+  // List of all vertex indices
+  VectorXi all,in;
+  igl::colon<int>(0,V.rows()-1,all);
+  // List of interior indices
+  igl::setdiff(all,b,in,IA);
+
+  // Construct and slice up Laplacian
+  SparseMatrix<double> L,L_in_in,L_in_b;
+  igl::cotmatrix(V,F,L);
+  igl::slice(L,in,in,L_in_in);
+  igl::slice(L,in,b,L_in_b);
+
+  // Dirichlet boundary conditions from z-coordinate
+  VectorXd bc;
+  VectorXd Z = V.col(2);
+  igl::slice(Z,b,bc);
+
+  // Solve PDE
+  SimplicialLLT<SparseMatrix<double > > solver(-L_in_in);
+  VectorXd Z_in = solver.solve(L_in_b*bc);
+  // slice into solution
+  igl::slice_into(Z_in,in,Z);
+
+  // Alternative, short hand
+  igl::min_quad_with_fixed_data<double> mqwf;
+  // Linear term is 0
+  VectorXd B = VectorXd::Zero(V.rows(),1);
+  // Empty constraints
+  VectorXd Beq;
+  SparseMatrix<double> Aeq;
+  // Our cotmatrix is _negative_ definite, so flip sign
+  igl::min_quad_with_fixed_precompute((-L).eval(),b,Aeq,true,mqwf);
+  igl::min_quad_with_fixed_solve(mqwf,B,bc,Beq,Z);
+
+  // Pseudo-color based on solution
+  MatrixXd C;
+  igl::jet(Z,true,C);
+
+  // Plot the mesh with pseudocolors
+  igl::Viewer viewer;
+  viewer.set_mesh(V, F);
+  viewer.options.show_lines = false;
+  viewer.set_colors(C);
+  viewer.launch();
+}

tutorial/304_LinearEqualityConstraints/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(304_LinearEqualityConstraints)
+
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${CMAKE_PROJECT_NAME} ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${CMAKE_PROJECT_NAME} ${SHARED_LIBRARIES})

tutorial/304_LinearEqualityConstraints/main.cpp

@@ -0,0 +1,95 @@
+#include <igl/boundary_faces.h>
+#include <igl/cotmatrix.h>
+#include <igl/invert_diag.h>
+#include <igl/jet.h>
+#include <igl/massmatrix.h>
+#include <igl/min_quad_with_fixed.h>
+#include <igl/readOFF.h>
+#include <igl/viewer/Viewer.h>
+#include <Eigen/Sparse>
+#include <iostream>
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+  MatrixXd V;
+  MatrixXi F;
+  igl::readOFF("../shared/cheburashka.off",V,F);
+
+  // Two fixed points
+  VectorXi b(2,1);
+  // Left hand, left foot
+  b<<4331,5957;
+  VectorXd bc(2,1);
+  bc<<1,-1;
+
+  // Construct Laplacian and mass matrix
+  SparseMatrix<double> L,M,Minv,Q;
+  igl::cotmatrix(V,F,L);
+  igl::massmatrix(V,F,igl::MASSMATRIX_VORONOI,M);
+  igl::invert_diag(M,Minv);
+  // Bi-Laplacian
+  Q = L * (Minv * L);
+  // Zero linear term
+  VectorXd B = VectorXd::Zero(V.rows(),1);
+
+  VectorXd Z,Z_const;
+  {
+    // Alternative, short hand
+    igl::min_quad_with_fixed_data<double> mqwf;
+    // Empty constraints
+    VectorXd Beq;
+    SparseMatrix<double> Aeq;
+    igl::min_quad_with_fixed_precompute(Q,b,Aeq,true,mqwf);
+    igl::min_quad_with_fixed_solve(mqwf,B,bc,Beq,Z);
+  }
+
+  {
+    igl::min_quad_with_fixed_data<double> mqwf;
+    // Constraint forcing difference of two points to be 0
+    SparseMatrix<double> Aeq(1,V.rows());
+    // Right hand, right foot
+    Aeq.insert(0,6074) = 1;
+    Aeq.insert(0,6523) = -1;
+    Aeq.makeCompressed();
+    VectorXd Beq(1,1);
+    Beq(0) = 0;
+    igl::min_quad_with_fixed_precompute(Q,b,Aeq,true,mqwf);
+    igl::min_quad_with_fixed_solve(mqwf,B,bc,Beq,Z_const);
+  }
+
+
+  // Pseudo-color based on solution
+  struct Data{
+    MatrixXd C,C_const;
+  } data;
+  // Use same color axes
+  const double min_z = std::min(Z.minCoeff(),Z_const.minCoeff());
+  const double max_z = std::max(Z.maxCoeff(),Z_const.maxCoeff());
+  igl::jet(      Z,min_z,max_z,data.C);
+  igl::jet(Z_const,min_z,max_z,data.C_const);
+
+  // Plot the mesh with pseudocolors
+  igl::Viewer viewer;
+  viewer.set_mesh(V, F);
+  viewer.options.show_lines = false;
+  viewer.set_colors(data.C);
+
+  viewer.callback_key_down = 
+    [](igl::Viewer& viewer,unsigned char key,int mod)->bool
+    {
+      if(key == ' ')
+      {
+        Data & data = *static_cast<Data*>(viewer.callback_key_down_data);
+        static bool toggle = true;
+        viewer.set_colors(toggle?data.C_const:data.C);
+        toggle = !toggle;
+        return true;
+      }else
+      {
+        return false;
+      }
+    };
+  viewer.callback_key_down_data = &data;
+  viewer.launch();
+}

tutorial/305_QuadraticProgramming/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(305_QuadraticProgramming)
+
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${CMAKE_PROJECT_NAME} ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${CMAKE_PROJECT_NAME} ${SHARED_LIBRARIES})

tutorial/305_QuadraticProgramming/main.cpp

@@ -0,0 +1,94 @@
+#include <igl/active_set.h>
+#include <igl/boundary_faces.h>
+#include <igl/cotmatrix.h>
+#include <igl/invert_diag.h>
+#include <igl/jet.h>
+#include <igl/massmatrix.h>
+#include <igl/readOFF.h>
+#include <igl/viewer/Viewer.h>
+#include <Eigen/Sparse>
+#include <iostream>
+
+  
+Eigen::VectorXi b;
+Eigen::VectorXd B,bc,lx,ux,Beq,Bieq,Z;
+Eigen::SparseMatrix<double> Q,Aeq,Aieq;
+
+void solve(igl::Viewer &viewer)
+{
+  using namespace std;
+  igl::active_set_params as;
+  as.max_iter = 8;
+  igl::active_set(Q,B,b,bc,Aeq,Beq,Aieq,Bieq,lx,ux,as,Z);
+  // Pseudo-color based on solution
+  Eigen::MatrixXd C;
+  igl::jet(Z,0,1,C);
+  viewer.set_colors(C);
+}
+
+bool key_down(igl::Viewer &viewer, unsigned char key, int mod)
+{
+  switch(key)
+  {
+    case '.':
+      Beq(0) *= 2.0;
+      solve(viewer);
+      return true;
+    case ',':
+      Beq(0) /= 2.0;
+      solve(viewer);
+      return true;
+    case ' ':
+      solve(viewer);
+      return true;
+    default:
+      return false;
+  }
+}
+
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+  MatrixXd V;
+  MatrixXi F;
+  igl::readOFF("../shared/cheburashka.off",V,F);
+
+  // Plot the mesh
+  igl::Viewer viewer;
+  viewer.set_mesh(V, F);
+  viewer.options.show_lines = false;
+  viewer.callback_key_down = &key_down;
+
+  // One fixed point
+  b.resize(1,1);
+  // point on belly.
+  b<<2556;
+  bc.resize(1,1);
+  bc<<1;
+
+  // Construct Laplacian and mass matrix
+  SparseMatrix<double> L,M,Minv;
+  igl::cotmatrix(V,F,L);
+  igl::massmatrix(V,F,igl::MASSMATRIX_VORONOI,M);
+  //M = (M/M.diagonal().maxCoeff()).eval();
+  igl::invert_diag(M,Minv);
+  // Bi-Laplacian
+  Q = L.transpose() * (Minv * L);
+  // Zero linear term
+  B = VectorXd::Zero(V.rows(),1);
+
+  // Lower and upper bound
+  lx = VectorXd::Zero(V.rows(),1);
+  ux = VectorXd::Ones(V.rows(),1);
+
+  // Equality constraint constrain solution to sum to 1
+  Beq.resize(1,1);
+  Beq(0) = 0.08;
+  Aeq = M.diagonal().transpose().sparseView();
+  // (Empty inequality constraints)
+  solve(viewer);
+
+  viewer.launch();
+}

tutorial/cmake/FindLIBIGL.cmake

@@ -48,6 +48,13 @@ if(LIBIGL_USE_STATIC_LIBRARY)
     set(LIBIGL_FOUND FALSE)
   endif(NOT LIBIGLMATLAB_LIBRARY)
   set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${LIBIGLMATLAB_LIBRARY})
+  find_package(Matlab REQUIRED)
+  if(MATLAB_FOUND)
+    set(LIBIGL_INCLUDE_DIR ${LIBIGL_INCLUDE_DIR}  ${MATLAB_INCLUDE_DIR})
+    set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES}  ${MATLAB_LIBRARIES})
+  else(MATLAB_FOUND)
+    set(LIBIGL_FOUND FALSE)
+  endif(MATLAB_FOUND)
   FIND_LIBRARY( LIBIGLSVD3X3_LIBRARY NAMES iglsvd3x3 PATHS ${LIBIGL_LIB_DIRS})
   if(NOT LIBIGLSVD3X3_LIBRARY)
     set(LIBIGL_FOUND FALSE)

tutorial/cmake/FindMATLAB.cmake

@@ -1,54 +1,210 @@
+# - this module looks for Matlab
+# Defines:
+#  MATLAB_INCLUDE_DIR: include path for mex.h, engine.h
+#  MATLAB_LIBRARIES:   required libraries: libmex, etc
+#  MATLAB_MEX_LIBRARY: path to libmex.lib
+#  MATLAB_MX_LIBRARY:  path to libmx.lib
+#  MATLAB_MAT_LIBRARY:  path to libmat.lib # added
+#  MATLAB_ENG_LIBRARY: path to libeng.lib
+#  MATLAB_ROOT: path to Matlab's root directory
+
+# This file is part of Gerardus
+#
+# This is a derivative work of file FindMatlab.cmake released with
+# CMake v2.8, because the original seems to be a bit outdated and
+# doesn't work with my Windows XP and Visual Studio 10 install
 #
-# Try to find a MATLAB installation
+# (Note that the original file does work for Ubuntu Natty)
-# Once done this will define
 #
-# MATLAB_FOUND
+# Author: Ramon Casero <rcasero@gmail.com>, Tom Doel
-# MATLAB_INCLUDE_DIR
+# Version: 0.2.3
-# MATLAB_LIBRARIES
+# $Rev$
+# $Date$
 #
+# The original file was copied from an Ubuntu Linux install
+# /usr/share/cmake-2.8/Modules/FindMatlab.cmake
 
-FIND_PATH(MATLAB_INCLUDE_DIR engine.h
+#=============================================================================
-      PATHS
+# Copyright 2005-2009 Kitware, Inc.
-      /Applications/MATLAB_R2012b.app/extern/include
+#
-      /Applications/MATLAB_R2013a.app/extern/include
+# Distributed under the OSI-approved BSD License (the "License");
-      /Applications/MATLAB_R2013b.app/extern/include
+# see accompanying file Copyright.txt for details.
-      /Applications/MATLAB_R2014a.app/extern/include
+#
-      /Applications/MATLAB_R2014b.app/extern/include
+# This software is distributed WITHOUT ANY WARRANTY; without even the
-      NO_DEFAULT_PATH)
+# implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
-
+# See the License for more information.
-FIND_LIBRARY(MATLAB_LIBRARY1 eng
+#=============================================================================
-  PATHS
+# (To distribute this file outside of CMake, substitute the full
-    /Applications/MATLAB_R2012b.app/bin/maci64
+#  License text for the above reference.)
-    /Applications/MATLAB_R2013a.app/bin/maci64
+
-    /Applications/MATLAB_R2013b.app/bin/maci64
+SET(MATLAB_FOUND 0)
-    /Applications/MATLAB_R2014a.app/bin/maci64
+IF(WIN32)
-    /Applications/MATLAB_R2014b.app/bin/maci64
+  # Search for a version of Matlab available, starting from the most modern one to older versions
-  PATH_SUFFIXES`
+  FOREACH(MATVER "7.14" "7.11" "7.10" "7.9" "7.8" "7.7" "7.6" "7.5" "7.4")
-    a
+    IF((NOT DEFINED MATLAB_ROOT) 
-    lib64
+        OR ("${MATLAB_ROOT}" STREQUAL "")
-    lib
+        OR ("${MATLAB_ROOT}" STREQUAL "/registry"))
-    dylib
+      GET_FILENAME_COMPONENT(MATLAB_ROOT
-    NO_DEFAULT_PATH
+        "[HKEY_LOCAL_MACHINE\\SOFTWARE\\MathWorks\\MATLAB\\${MATVER};MATLABROOT]"
-)
+        ABSOLUTE)
+      SET(MATLAB_VERSION ${MATVER})
+    ENDIF((NOT DEFINED MATLAB_ROOT) 
+      OR ("${MATLAB_ROOT}" STREQUAL "")
+      OR ("${MATLAB_ROOT}" STREQUAL "/registry"))
+  ENDFOREACH(MATVER)
+  
+  # Directory name depending on whether the Windows architecture is 32
+  # bit or 64 bit
+  set(CMAKE_SIZEOF_VOID_P 8) # Note: For some weird reason this variable is undefined in my system...
+  IF(CMAKE_SIZEOF_VOID_P MATCHES "4")
+    SET(WINDIR "win32")
+  ELSEIF(CMAKE_SIZEOF_VOID_P MATCHES "8")
+    SET(WINDIR "win64")
+  ELSE(CMAKE_SIZEOF_VOID_P MATCHES "4")
+    MESSAGE(FATAL_ERROR 
+      "CMAKE_SIZEOF_VOID_P (${CMAKE_SIZEOF_VOID_P}) doesn't indicate a valid platform")
+  ENDIF(CMAKE_SIZEOF_VOID_P MATCHES "4")
+
+  # Folder where the MEX libraries are, depending of the Windows compiler
+  IF(${CMAKE_GENERATOR} MATCHES "Visual Studio 6")
+    SET(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/msvc60")
+  ELSEIF(${CMAKE_GENERATOR} MATCHES "Visual Studio 7")
+    # Assume people are generally using Visual Studio 7.1,
+    # if using 7.0 need to link to: ../extern/lib/${WINDIR}/microsoft/msvc70
+    SET(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/msvc71")
+    # SET(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/msvc70")
+  ELSEIF(${CMAKE_GENERATOR} MATCHES "Borland")
+    # Assume people are generally using Borland 5.4,
+    # if using 7.0 need to link to: ../extern/lib/${WINDIR}/microsoft/msvc70
+    SET(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/bcc54")
+    # SET(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/bcc50")
+    # SET(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/bcc51")
+  ELSEIF(${CMAKE_GENERATOR} MATCHES "Visual Studio*")
+    # If the compiler is Visual Studio, but not any of the specific
+    # versions above, we try our luck with the microsoft directory
+    SET(MATLAB_LIBRARIES_DIR "${MATLAB_ROOT}/extern/lib/${WINDIR}/microsoft/")
+  ELSE(${CMAKE_GENERATOR} MATCHES "Visual Studio 6")
+    MESSAGE(FATAL_ERROR "Generator not compatible: ${CMAKE_GENERATOR}")
+  ENDIF(${CMAKE_GENERATOR} MATCHES "Visual Studio 6")
+
+  # Get paths to the Matlab MEX libraries
+  FIND_LIBRARY(MATLAB_MEX_LIBRARY
+    libmex
+    ${MATLAB_LIBRARIES_DIR}
+    )
+  FIND_LIBRARY(MATLAB_MX_LIBRARY
+    libmx
+    ${MATLAB_LIBRARIES_DIR}
+    )
+  FIND_LIBRARY(MATLAB_MAT_LIBRARY
+    libmat
+    ${MATLAB_LIBRARIES_DIR}
+    )
+  FIND_LIBRARY(MATLAB_ENG_LIBRARY
+    libeng
+    ${MATLAB_LIBRARIES_DIR}
+    )
+
+  # Get path to the include directory
+  FIND_PATH(MATLAB_INCLUDE_DIR
+    "mex.h"
+    "${MATLAB_ROOT}/extern/include"
+    )
+
+ELSE(WIN32)
+
+  IF((NOT DEFINED MATLAB_ROOT) 
+      OR ("${MATLAB_ROOT}" STREQUAL ""))
+    # get path to the Matlab root directory
+    EXECUTE_PROCESS(
+      COMMAND which matlab
+      COMMAND xargs readlink
+      COMMAND xargs dirname
+      COMMAND xargs dirname
+      COMMAND xargs echo -n
+      OUTPUT_VARIABLE MATLAB_ROOT
+      )
+  ENDIF((NOT DEFINED MATLAB_ROOT) 
+    OR ("${MATLAB_ROOT}" STREQUAL ""))
+
+  # Check if this is a Mac
+  IF(${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
 
-FIND_LIBRARY(MATLAB_LIBRARY2 mx
+    SET(LIBRARY_EXTENSION .dylib)
-  PATHS
+
-    /Applications/MATLAB_R2012b.app/bin/maci64/
+    # If this is a Mac and the attempts to find MATLAB_ROOT have so far failed, 
-    /Applications/MATLAB_R2013a.app/bin/maci64/
+    # we look in the applications folder
-    /Applications/MATLAB_R2013b.app/bin/maci64/
+    IF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
-    /Applications/MATLAB_R2014a.app/bin/maci64/
+
-    /Applications/MATLAB_R2014b.app/bin/maci64/
+    # Search for a version of Matlab available, starting from the most modern one to older versions
-  PATH_SUFFIXES
+      FOREACH(MATVER "R2014a" "R2014a" "R2013b" "R2013a" "R2012b" "R2012a" "R2011b" "R2011a" "R2010b" "R2010a" "R2009b" "R2009a" "R2008b")
-    a
+        IF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
-    lib64
+          IF(EXISTS /Applications/MATLAB_${MATVER}.app)
-    lib
+            SET(MATLAB_ROOT /Applications/MATLAB_${MATVER}.app)
-    dylib
+    
-    NO_DEFAULT_PATH
+          ENDIF(EXISTS /Applications/MATLAB_${MATVER}.app)
+        ENDIF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
+      ENDFOREACH(MATVER)
+
+    ENDIF((NOT DEFINED MATLAB_ROOT) OR ("${MATLAB_ROOT}" STREQUAL ""))
+
+  ELSE(${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
+    SET(LIBRARY_EXTENSION .so)
+
+  ENDIF(${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
+
+  # Get path to the MEX libraries
+  EXECUTE_PROCESS(
+    #COMMAND find "${MATLAB_ROOT}/extern/lib" -name libmex${LIBRARY_EXTENSION} # Peter
+	COMMAND find "${MATLAB_ROOT}/bin" -name libmex${LIBRARY_EXTENSION} # Standard
+    COMMAND xargs echo -n
+    OUTPUT_VARIABLE MATLAB_MEX_LIBRARY
+    )
+  EXECUTE_PROCESS(
+    #COMMAND find "${MATLAB_ROOT}/extern/lib" -name libmx${LIBRARY_EXTENSION} # Peter
+	COMMAND find "${MATLAB_ROOT}/bin" -name libmx${LIBRARY_EXTENSION} # Standard
+    COMMAND xargs echo -n
+    OUTPUT_VARIABLE MATLAB_MX_LIBRARY
+    )
+  EXECUTE_PROCESS(
+    #COMMAND find "${MATLAB_ROOT}/extern/lib" -name libmat${LIBRARY_EXTENSION} # Peter
+	COMMAND find "${MATLAB_ROOT}/bin" -name libmat${LIBRARY_EXTENSION} # Standard
+    COMMAND xargs echo -n
+    OUTPUT_VARIABLE MATLAB_MAT_LIBRARY
+    )
+  EXECUTE_PROCESS(
+    #COMMAND find "${MATLAB_ROOT}/extern/lib" -name libeng${LIBRARY_EXTENSION} # Peter
+	COMMAND find "${MATLAB_ROOT}/bin" -name libeng${LIBRARY_EXTENSION} # Standard
+    COMMAND xargs echo -n
+    OUTPUT_VARIABLE MATLAB_ENG_LIBRARY
+    )
+
+  # Get path to the include directory
+  FIND_PATH(MATLAB_INCLUDE_DIR
+    "mex.h"
+    PATHS "${MATLAB_ROOT}/extern/include"
+    )
+
+ENDIF(WIN32)
+
+# This is common to UNIX and Win32:
+SET(MATLAB_LIBRARIES
+  ${MATLAB_MAT_LIBRARY}
+  ${MATLAB_MEX_LIBRARY}
+  ${MATLAB_MX_LIBRARY}
+  ${MATLAB_ENG_LIBRARY}
 )
 
-if(MATLAB_INCLUDE_DIR AND MATLAB_LIBRARY1 AND MATLAB_LIBRARY2)
+IF(MATLAB_INCLUDE_DIR AND MATLAB_LIBRARIES)
-	message(STATUS "Found MATLAB: ${MATLAB_INCLUDE_DIR}")
+  SET(MATLAB_FOUND 1)
-  set(MATLAB_LIBRARIES ${MATLAB_LIBRARY1} ${MATLAB_LIBRARY2})
+ENDIF(MATLAB_INCLUDE_DIR AND MATLAB_LIBRARIES)
-else(MATLAB_INCLUDE_DIR AND MATLAB_LIBRARY1 AND MATLAB_LIBRARY2)
+
-	message(FATAL_ERROR "could NOT find MATLAB")
+MARK_AS_ADVANCED(
-endif(MATLAB_INCLUDE_DIR AND MATLAB_LIBRARY1 AND MATLAB_LIBRARY2)
+  MATLAB_LIBRARIES
+  MATLAB_MEX_LIBRARY
+  MATLAB_MX_LIBRARY
+  MATLAB_ENG_LIBRARY
+  MATLAB_INCLUDE_DIR
+  MATLAB_FOUND
+  MATLAB_ROOT
+)

tutorial/images/camelhead-laplace-equation.jpg.REMOVED.git-id

@@ -0,0 +1 @@
+39ddde658cbb1e577caeb58083724a81f796956e

tutorial/images/cheburashka-biharmonic-leq.jpg.REMOVED.git-id

@@ -0,0 +1 @@
+3254fdfb5863a1fe149c14f0630a90c0b942d842

tutorial/images/cheburashka-multiscale-biharmonic-kernels.jpg.REMOVED.git-id

@@ -0,0 +1 @@
+e81844bdd47fc074d33e369495ea84569e66856d

tutorial/readme.md

@@ -37,11 +37,23 @@ of these lecture notes links to a cross-platform example application.
     * [204 Laplacian](#laplacian)
         * [Mass matrix](#massmatrix)
         * [Alternative construction of
-          Laplacian](#alternativeconstuctionoflaplacian)
+          Laplacian](#alternativeconstructionoflaplacian)
 * [Chapter 3: Matrices and Linear Algebra](#chapter3:matricesandlinearalgebra)
     * [301 Slice](#slice)
-    * [301 Sort](#sort)
+    * [302 Sort](#sort)
         * [Other Matlab-style functions](#othermatlab-stylefunctions) 
+    * [303 Laplace Equation](#laplaceequation)
+        * [Quadratic energy minimization](#quadraticenergyminimization)
+    * [304 Linear Equality Constraints](#linearequalityconstraints)
+    * [305 Quadratic Programming](#quadraticprogramming)
+* [Chapter 4: Shape Deformation](#chapter4:shapedeformation)
+    * [401 Biharmonic Deformation](#biharmonicdeformation)
+    * [402 Bounded Biharmonic Weights](#boundedbiharmonicweights)
+    * [403 Dual Quaternion Skinning](#dualquaternionskinning)
+    * [404 As-rigid-as-possible](#as-rigid-as-possible)
+    * [405 Fast automatic skinning
+      transformations](#fastautomaticskinningtransformations)
+
 
 # Compilation Instructions
 
@@ -367,7 +379,7 @@ An alternative construction of the discrete cotangent Laplacian is by
 "squaring" the discrete gradient operator. This may be derived by applying
 Green's identity (ignoring boundary conditions for the moment):
 
-  $\int_S \nabla f \nabla f dA = \int_S f \Delta f dA$
+  $\int_S \|\nabla f\|^2 dA = \int_S f \Delta f dA$
 
 Or in matrix form which is immediately translatable to code:
 
@@ -493,8 +505,263 @@ functionality as common matlab functions.
 - `igl::median` Compute the median per column
 - `igl::mode` Compute the mode per column
 - `igl::orth` Orthogonalization of a basis
+- `igl::setdiff` Set difference of matrix elements
 - `igl::speye` Identity as sparse matrix
 
+## Laplace Equation
+A common linear system in geometry processing is the Laplace equation:
+
+ $∆z = 0$
+
+subject to some boundary conditions, for example Dirichlet boundary conditions
+(fixed value):
+
+ $\left.z\right|_{\partial{S}} = z_{bc}$
+
+In the discrete setting, this begins with the linear system:
+
+ $\mathbf{L} \mathbf{z} = \mathbf{0}$
+
+where $\mathbf{L}$ is the $n \times n$ discrete Laplacian and $\mathbf{z}$ is a
+vector of per-vertex values. Most of $\mathbf{z}$ correspond to interior
+vertices and are unknown, but some of $\mathbf{z}$ represent values at boundary
+vertices. Their values are known so we may move their corresponding terms to
+the right-hand side.
+
+Conceptually, this is very easy if we have sorted $\mathbf{z}$ so that interior
+vertices come first and then boundary vertices:
+
+ $$\left(\begin{array}{cc}
+ \mathbf{L}_{in,in} & \mathbf{L}_{in,b}\\
+ \mathbf{L}_{b,in} & \mathbf{L}_{b,b}\end{array}\right) 
+ \left(\begin{array}{c}
+ \mathbf{z}_{in}\\
+ \mathbf{L}_{b}\end{array}\right) = 
+ \left(\begin{array}{c}
+ \mathbf{0}_{in}\\
+ \mathbf{*}_{b}\end{array}\right)$$ 
+
+The bottom block of equations is no longer meaningful so we'll only consider
+the top block:
+
+ $$\left(\begin{array}{cc}
+ \mathbf{L}_{in,in} & \mathbf{L}_{in,b}\end{array}\right) 
+ \left(\begin{array}{c}
+ \mathbf{z}_{in}\\
+ \mathbf{z}_{b}\end{array}\right) = 
+ \mathbf{0}_{in}$$
+
+Where now we can move known values to the right-hand side:
+
+ $$\mathbf{L}_{in,in} 
+ \mathbf{z}_{in} = -
+ \mathbf{L}_{in,b}
+ \mathbf{z}_{b}$$
+
+Finally we can solve this equation for the unknown values at interior vertices
+$\mathbf{z}_{in}$.
+
+However, probably our vertices are not sorted. One option would be to sort `V`,
+then proceed as above and then _unsort_ the solution `Z` to match `V`. However,
+this solution is not very general.
+
+With array slicing no explicit sort is needed. Instead we can _slice-out_
+submatrix blocks ($\mathbf{L}_{in,in}$, $\mathbf{L}_{in,b}$, etc.) and follow
+the linear algebra above directly. Then we can slice the solution _into_ the
+rows of `Z` corresponding to the interior vertices.
+
+![The `LaplaceEquation` example solves a Laplace equation with Dirichlet
+boundary conditions.](images/camelhead-laplace-equation.jpg)
+
+### Quadratic energy minimization
+
+The same Laplace equation may be equivalently derived by minimizing Dirichlet
+energy subject to the same boundary conditions:
+
+ $\mathop{\text{minimize }}_z \frac{1}{2}\int\limits_S \|\nabla z\|^2 dA$
+
+On our discrete mesh, recall that this becomes
+
+ $\mathop{\text{minimize }}_\mathbf{z}  \frac{1}{2}\mathbf{z}^T \mathbf{G}^T \mathbf{D}
+ \mathbf{G} \mathbf{z} \rightarrow \mathop{\text{minimize }}_\mathbf{z} \mathbf{z}^T \mathbf{L} \mathbf{z}$
+
+The general problem of minimizing some energy over a mesh subject to fixed
+value boundary conditions is so wide spread that libigl has a dedicated api for
+solving such systems. 
+
+Let's consider a general quadratic minimization problem subject to different
+common constraints:
+
+ $$\mathop{\text{minimize }}_\mathbf{z}  \frac{1}{2}\mathbf{z}^T \mathbf{Q} \mathbf{z} +
+ \mathbf{z}^T \mathbf{B} + \text{constant},$$
+
+ subject to
+ 
+ $$\mathbf{z}_b = \mathbf{z}_{bc} \text{ and } \mathbf{A}_{eq} \mathbf{z} =
+ \mathbf{B}_{eq},$$
+
+where 
+
+  - $\mathbf{Q}$ is a (usually sparse) $n \times n$ positive semi-definite
+    matrix of quadratic coefficients (Hessian), 
+  - $\mathbf{B}$ is a $n \times 1$ vector of linear coefficients, 
+  - $\mathbf{z}_b$ is a $|b| \times 1$ portion of
+$\mathbf{z}$ corresponding to boundary or _fixed_ vertices,
+  - $\mathbf{z}_{bc}$ is a $|b| \times 1$ vector of known values corresponding to
+    $\mathbf{z}_b$,
+  - $\mathbf{A}_{eq}$ is a (usually sparse) $m \times n$ matrix of linear
+    equality constraint coefficients (one row per constraint), and
+  - $\mathbf{B}_{eq}$ is a $m \times 1$ vector of linear equality constraint
+    right-hand side values.
+
+This specification is overly general as we could write $\mathbf{z}_b =
+\mathbf{z}_{bc}$ as rows of $\mathbf{A}_{eq} \mathbf{z} =
+\mathbf{B}_{eq}$, but these fixed value constraints appear so often that they
+merit a dedicated place in the API.
+
+In libigl, solving such quadratic optimization problems is split into two
+routines: precomputation and solve. Precomputation only depends on the
+quadratic coefficients, known value indices and linear constraint coefficients:
+
+```cpp
+igl::min_quad_with_fixed_data mqwf;
+igl::min_quad_with_fixed_precompute(Q,b,Aeq,true,mqwf);
+```
+
+The output is a struct `mqwf` which contains the system matrix factorization
+and is used during solving with arbitrary linear terms, known values, and
+constraint right-hand sides:
+
+```cpp
+igl::min_quad_with_fixed_solve(mqwf,B,bc,Beq,Z);
+```
+
+The output `Z` is a $n \times 1$ vector of solutions with fixed values
+correctly placed to match the mesh vertices `V`.
+
+## Linear Equality Constraints
+We saw above that `min_quad_with_fixed_*` in libigl provides a compact way to
+solve general quadratic programs. Let's consider another example, this time
+with active linear equality constraints. Specifically let's solve the
+`bi-Laplace equation` or equivalently minimize the Laplace energy:
+
+ $$\Delta^2 z = 0 \leftrightarrow \mathop{\text{minimize }}\limits_z \frac{1}{2}
+ \int\limits_S (\Delta z)^2 dA$$
+
+subject to fixed value constraints and a linear equality constraint:
+
+ $z_{a} = 1, z_{b} = -1$ and $z_{c} = z_{d}$.
+
+Notice that we can rewrite the last constraint in the familiar form from above:
+
+ $z_{c} - z_{d} = 0.$
+
+Now we can assembly `Aeq` as a $1 \times n$ sparse matrix with a coefficient
+$1$
+in the column corresponding to vertex $c$ and a $-1$ at $d$. The right-hand side
+`Beq` is simply zero.
+
+Internally, `min_quad_with_fixed_*` solves using the Lagrange Multiplier
+method. This method adds additional variables for each linear constraint (in
+general a $m \times 1$ vector of variables $\lambda$) and then solves the
+saddle problem:
+
+ $$\mathop{\text{find saddle }}_{\mathbf{z},\lambda}\, \frac{1}{2}\mathbf{z}^T \mathbf{Q} \mathbf{z} +
+  \mathbf{z}^T \mathbf{B} + \text{constant} + \lambda^T\left(\mathbf{A}_{eq}
+ \mathbf{z} - \mathbf{B}_{eq}\right)$$
+
+This can be rewritten in a more familiar form by stacking $\mathbf{z}$ and
+$\lambda$ into one $(m+n) \times 1$ vector of unknowns:
+
+ $$\mathop{\text{find saddle }}_{\mathbf{z},\lambda}\, 
+ \frac{1}{2}
+ \left(
+  \mathbf{z}^T 
+  \lambda^T
+ \right)
+ \left(
+  \begin{array}{cc}
+  \mathbf{Q}      & \mathbf{A}_{eq}^T\\
+  \mathbf{A}_{eq} & 0
+  \end{array}
+ \right)
+ \left(
+  \begin{array}{c}
+  \mathbf{z}\\
+  \lambda
+  \end{array}
+ \right) + 
+ \left(
+  \mathbf{z}^T 
+  \lambda^T
+ \right)
+ \left(
+  \begin{array}{c}
+  \mathbf{B}\\
+  -\mathbf{B}_{eq}
+  \end{array}
+  \right)
+  + \text{constant}$$
+
+Differentiating with respect to $\left( \mathbf{z}^T \lambda^T \right)$ reveals
+a linear system and we can solve for $\mathbf{z}$ and $\lambda$. The only
+difference from
+the straight quadratic
+_minimization_ system, is that 
+this saddle problem system will not be positive definite. Thus, we must use a
+different factorization technique (LDLT rather than LLT). Luckily, libigl's
+`min_quad_with_fixed_precompute` automatically chooses the correct solver in
+the presence of linear equality constraints.
+
+![The example `LinearEqualityConstraints` first solves with just fixed value
+constraints (left: 1 and -1 on the left hand and foot respectively), then
+solves with an additional linear equality constraint (right: points on right
+hand and foot constrained to be equal).](images/cheburashka-biharmonic-leq.jpg)
+
+## Quadratic Programming
+
+We can generalize the quadratic optimization in the previous section even more
+by allowing inequality constraints. Specifically box constraints (lower and
+upper bounds):
+
+ $\mathbf{l} \le \mathbf{z} \le \mathbf{u},$
+
+where $\mathbf{l},\mathbf{u}$ are $n \times 1$ vectors of lower and upper
+bounds
+and general linear inequality constraints:
+
+ $\mathbf{A}_{ieq} \mathbf{z} \le \mathbf{B}_{ieq},$
+
+where $\mathbf{A}_{ieq}$ is a $k \times n$ matrix of linear coefficients and
+$\mathbf{B}_{ieq}$ is a $k \times 1$ matrix of constraint right-hand sides.
+
+Again, we are overly general as the box constraints could be written as
+rows of the linear inequality constraints, but bounds appear frequently enough
+to merit a dedicated api.
+
+Libigl implements its own active set routine for solving _quadratric programs_
+(QPs). This algorithm works by iteratively "activating" violated inequality
+constraints by enforcing them as equalities and "deactivating" constraints
+which are no longer needed.
+
+After deciding which constraints are active each iteration simple reduces to a
+quadratic minimization subject to linear _equality_ constraints, and the method
+from the previous section is invoked. This is repeated until convergence.
+
+Currently the implementation is efficient for box constraints and sparse
+non-overlapping linear inequality constraints.
+
+Unlike alternative interior-point methods, the active set method benefits from
+a warm-start (initial guess for the solution vector $\mathbf{z}$).
+
+```cpp
+igl::active_set_params as;
+// Z is optional initial guess and output
+igl::active_set(Q,B,b,bc,Aeq,Beq,Aieq,Bieq,lx,ux,as,Z);
+```
+
+![The example `QuadraticProgramming` uses an active set solver to optimize
+discrete biharmonic kernels at multiple scales [#rustamov_2011][].](images/cheburashka-multiscale-biharmonic-kernels.jpg)
 
 [#meyer_2003]: Mark Meyer and Mathieu Desbrun and Peter Schröder and Alan H.  Barr,
 "Discrete Differential-Geometry Operators for Triangulated
@@ -506,3 +773,4 @@ _Algorithms and Interfaces for Real-Time Deformation of 2D and 3D Shapes_,
 2013.
 [#kazhdan_2012]: Michael Kazhdan, Jake Solomon, Mirela Ben-Chen,
 "Can Mean-Curvature Flow Be Made Non-Singular," 2012.
+[#rustamov_2011]: Raid M. Rustamov, "Multiscale Biharmonic Kernels", 2011.