linear equality constraints

Former-commit-id: 9d193ada6f8afa42e8f7e7130115ae947df49f57

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/jet.cpp

@@ -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

tutorial/303_LaplaceEquation/main.cpp

@@ -36,7 +36,7 @@ int main(int argc, char *argv[])
   igl::slice(L,in,in,L_in_in);
   igl::slice(L,in,b,L_in_b);
 
-  // Dirichelet boundary conditions from z-coordinate
+  // Dirichlet boundary conditions from z-coordinate
   VectorXd bc;
   VectorXd Z = V.col(2);
   igl::slice(Z,b,bc);

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/images/cheburashka-biharmonic-leq.jpg.REMOVED.git-id

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

tutorial/readme.md

@@ -40,10 +40,11 @@ of these lecture notes links to a cross-platform example application.
           Laplacian](#alternativeconstuctionoflaplacian)
 * [Chapter 3: Matrices and Linear Algebra](#chapter3:matricesandlinearalgebra)
     * [301 Slice](#slice)
-    * [301 Sort](#sort)
+    * [302 Sort](#sort)
         * [Other Matlab-style functions](#othermatlab-stylefunctions) 
-    * [301 Laplace Equation](#laplaceequation)
+    * [303 Laplace Equation](#laplaceequation)
         * [Quadratic energy minimization](#quadraticenergyminimization)
+    * [304 Linear Equality Constraints](#linearequalityconstraints)
 
 # Compilation Instructions
 
@@ -568,11 +569,11 @@ boundary conditions.](images/camelhead-laplace-equation.jpg)
 The same Laplace equation may be equivalently derived by minimizing Dirichlet
 energy subject to the same boundary conditions:
 
- $\mathop{\text{minimize }}_z \int\limits_S \|\nabla z\|^2 dA$
+ $\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}  \mathbf{z}^T \mathbf{G}^T \mathbf{D}
+ $\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
@@ -582,7 +583,7 @@ solving such systems.
 Let's consider a general quadratic minimization problem subject to different
 common constraints:
 
- $$\mathop{\text{minimize }}_\mathbf{z}  \mathbf{z}^T \mathbf{Q} \mathbf{z} +
+ $$\mathop{\text{minimize }}_\mathbf{z}  \frac{1}{2}\mathbf{z}^T \mathbf{Q} \mathbf{z} +
  \mathbf{z}^T \mathbf{B} + \text{constant},$$
 
  subject to
@@ -629,6 +630,86 @@ 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
+
 
 [#meyer_2003]: Mark Meyer and Mathieu Desbrun and Peter Schröder and Alan H.  Barr,
 "Discrete Differential-Geometry Operators for Triangulated