edge collapsing and linear programs

Former-commit-id: e27876824d90a57356f90bae045956fac6fa0225

include/igl/circulation.cpp

@@ -0,0 +1,71 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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 "circulation.h"
+#include "list_to_matrix.h"
+
+IGL_INLINE std::vector<int> igl::circulation(
+  const int e,
+  const bool ccw,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXi & E,
+  const Eigen::VectorXi & EMAP,
+  const Eigen::MatrixXi & EF,
+  const Eigen::MatrixXi & EI)
+{
+  // prepare output
+  std::vector<int> N;
+  N.reserve(6);
+  const int m = F.rows();
+  const auto & step = [&](
+    const int e, 
+    const int ff,
+    int & ne, 
+    int & nf)
+  {
+    assert((EF(e,1) == ff || EF(e,0) == ff) && "e should touch ff");
+    //const int fside = EF(e,1)==ff?1:0;
+    const int nside = EF(e,0)==ff?1:0;
+    const int nv = EI(e,nside);
+    // get next face
+    nf = EF(e,nside);
+    // get next edge 
+    const int dir = ccw?-1:1;
+    ne = EMAP(nf+m*((nv+dir+3)%3));
+  };
+  // Always start with first face (ccw in step will be sure to turn right
+  // direction)
+  const int f0 = EF(e,0);
+  int fi = f0;
+  int ei = e;
+  while(true)
+  {
+    N.push_back(fi);
+    step(ei,fi,ei,fi);
+    // back to start?
+    if(fi == f0)
+    {
+      assert(ei == e);
+      break;
+    }
+  }
+  return N;
+}
+
+IGL_INLINE void igl::circulation(
+  const int e,
+  const bool ccw,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXi & E,
+  const Eigen::VectorXi & EMAP,
+  const Eigen::MatrixXi & EF,
+  const Eigen::MatrixXi & EI,
+  Eigen::VectorXi & vN)
+{
+  std::vector<int> N = circulation(e,ccw,F,E,EMAP,EF,EI);
+  igl::list_to_matrix(N,vN);
+}

include/igl/circulation.h

@@ -0,0 +1,56 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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_CIRCULATION_H
+#define IGL_CIRCULATION_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+namespace igl
+{
+  // Return list of faces around the end point of an edge. Assumes
+  // data-structures are built from an edge-manifold **closed** mesh.
+  //
+  // Inputs:
+  //   e  index into E of edge to circulate
+  //   ccw  whether to _continue_ in ccw direction of edge (circulate around
+  //     E(e,1))
+  //   F  #F by 3 list of face indices
+  //   E  #E by 2 list of edge indices
+  //   EMAP #F*3 list of indices into E, mapping each directed edge to unique
+  //     unique edge in E
+  //   EF  #E by 2 list of edge flaps, EF(e,0)=f means e=(i-->j) is the edge of
+  //     F(f,:) opposite the vth corner, where EI(e,0)=v. Similarly EF(e,1) "
+  //     e=(j->i)
+  //   EI  #E by 2 list of edge flap corners (see above).
+  // Returns list of faces touched by circulation.
+  //   
+  IGL_INLINE std::vector<int> circulation(
+    const int e,
+    const bool ccw,
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXi & E,
+    const Eigen::VectorXi & EMAP,
+    const Eigen::MatrixXi & EF,
+    const Eigen::MatrixXi & EI);
+  // Wrapper with VectorXi output.
+  IGL_INLINE void circulation(
+    const int e,
+    const bool ccw,
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXi & E,
+    const Eigen::VectorXi & EMAP,
+    const Eigen::MatrixXi & EF,
+    const Eigen::MatrixXi & EI,
+    Eigen::VectorXi & vN);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "circulation.cpp"
+#endif
+#endif

include/igl/collapse_edge.cpp

@@ -0,0 +1,156 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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 "collapse_edge.h"
+#include "circulation.h"
+#include "edge_collapse_is_valid.h"
+#include <vector>
+
+IGL_INLINE bool igl::collapse_edge(
+  const int e,
+  const Eigen::RowVectorXd & p,
+  Eigen::MatrixXd & V,
+  Eigen::MatrixXi & F,
+  Eigen::MatrixXi & E,
+  Eigen::VectorXi & EMAP,
+  Eigen::MatrixXi & EF,
+  Eigen::MatrixXi & EI,
+  int & a_e1,
+  int & a_e2,
+  int & a_f1,
+  int & a_f2)
+{
+  // Assign this to 0 rather than, say, -1 so that deleted elements will get
+  // draw as degenerate elements at vertex 0 (which should always exist and
+  // never get collapsed to anything else since it is the smallest index)
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  const int eflip = E(e,0)>E(e,1);
+  // source and destination
+  const int s = eflip?E(e,1):E(e,0);
+  const int d = eflip?E(e,0):E(e,1);
+
+  if(!edge_collapse_is_valid(e,F,E,EMAP,EF,EI))
+  {
+    return false;
+  }
+
+  // Important to grab neighbors of d before monkeying with edges
+  const std::vector<int> nV2Fd = circulation(e,!eflip,F,E,EMAP,EF,EI);
+
+  // The following implementation strongly relies on s<d
+  assert(s<d && "s should be less than d");
+  // move source and destination to midpoint
+  V.row(s) = p;
+  V.row(d) = p;
+
+  // Helper function to replace edge and associate information with NULL
+  const auto & kill_edge = [&E,&EI,&EF](const int e)
+  {
+    E(e,0) = IGL_COLLAPSE_EDGE_NULL;
+    E(e,1) = IGL_COLLAPSE_EDGE_NULL;
+    EF(e,0) = IGL_COLLAPSE_EDGE_NULL;
+    EF(e,1) = IGL_COLLAPSE_EDGE_NULL;
+    EI(e,0) = IGL_COLLAPSE_EDGE_NULL;
+    EI(e,1) = IGL_COLLAPSE_EDGE_NULL;
+  };
+
+  // update edge info
+  // for each flap
+  const int m = F.rows();
+  for(int side = 0;side<2;side++)
+  {
+    const int f = EF(e,side);
+    const int v = EI(e,side);
+    const int sign = (eflip==0?1:-1)*(1-2*side);
+    // next edge emanating from d
+    const int e1 = EMAP(f+m*((v+sign*1+3)%3));
+    // prev edge pointing to s
+    const int e2 = EMAP(f+m*((v+sign*2+3)%3));
+    assert(E(e1,0) == d || E(e1,1) == d);
+    assert(E(e2,0) == s || E(e2,1) == s);
+    // face adjacent to f on e1, also incident on d
+    const bool flip1 = EF(e1,1)==f;
+    const int f1 = flip1 ? EF(e1,0) : EF(e1,1);
+    assert(f1!=f);
+    assert(F(f1,0)==d || F(f1,1)==d || F(f1,2) == d);
+    // across from which vertex of f1 does e1 appear?
+    const int v1 = flip1 ? EI(e1,0) : EI(e1,1);
+    // Kill e1
+    kill_edge(e1);
+    // Kill f
+    F(f,0) = IGL_COLLAPSE_EDGE_NULL;
+    F(f,1) = IGL_COLLAPSE_EDGE_NULL;
+    F(f,2) = IGL_COLLAPSE_EDGE_NULL;
+    // map f1's edge on e1 to e2
+    assert(EMAP(f1+m*v1) == e1);
+    EMAP(f1+m*v1) = e2;
+    // side opposite f2, the face adjacent to f on e2, also incident on s
+    const int opp2 = (EF(e2,0)==f?0:1);
+    assert(EF(e2,opp2) == f);
+    EF(e2,opp2) = f1;
+    EI(e2,opp2) = v1;
+    // remap e2 from d to s
+    E(e2,0) = E(e2,0)==d ? s : E(e2,0);
+    E(e2,1) = E(e2,1)==d ? s : E(e2,1);
+    if(side==0)
+    {
+      a_e1 = e1;
+      a_f1 = f;
+    }else
+    {
+      a_e2 = e1;
+      a_f2 = f;
+    }
+  }
+
+  // finally, reindex faces and edges incident on d. Do this last so asserts
+  // make sense.
+  //
+  // Could actually skip first two, since those are always the two collpased
+  // faces.
+  for(auto f : nV2Fd)
+  {
+    for(int v = 0;v<3;v++)
+    {
+      if(F(f,v) == d)
+      {
+        const int flip1 = (EF(EMAP(f+m*((v+1)%3)),0)==f)?1:0;
+        const int flip2 = (EF(EMAP(f+m*((v+2)%3)),0)==f)?0:1;
+        assert(
+          E(EMAP(f+m*((v+1)%3)),flip1) == d ||
+          E(EMAP(f+m*((v+1)%3)),flip1) == s);
+        E(EMAP(f+m*((v+1)%3)),flip1) = s;
+        assert(
+          E(EMAP(f+m*((v+2)%3)),flip2) == d ||
+          E(EMAP(f+m*((v+2)%3)),flip2) == s);
+        E(EMAP(f+m*((v+2)%3)),flip2) = s;
+        F(f,v) = s;
+        break;
+      }
+    }
+  }
+  // Finally, "remove" this edge and its information
+  kill_edge(e);
+
+  return true;
+}
+
+IGL_INLINE bool igl::collapse_edge(
+  const int e,
+  const Eigen::RowVectorXd & p,
+  Eigen::MatrixXd & V,
+  Eigen::MatrixXi & F,
+  Eigen::MatrixXi & E,
+  Eigen::VectorXi & EMAP,
+  Eigen::MatrixXi & EF,
+  Eigen::MatrixXi & EI)
+{
+  int e1,e2,f1,f2;
+  return collapse_edge(e,p,V,F,E,EMAP,EF,EI,e1,e2,f1,f2);
+}

include/igl/collapse_edge.h

@@ -0,0 +1,70 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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_COLLAPSE_EDGE_H
+#define IGL_COLLAPSE_EDGE_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Assumes (V,F) is a closed manifold mesh (except for previouslly collapsed
+  // faces which should be set to: 
+  // [IGL_COLLAPSE_EDGE_NULL IGL_COLLAPSE_EDGE_NULL IGL_COLLAPSE_EDGE_NULL].
+  // Collapses exactly two faces and exactly 3 edges from E (e and one side of
+  // each face gets collapsed to the other). This is implemented in a way that
+  // it can be repeatedly called until satisfaction and then the garbage in F
+  // can be collected by removing NULL faces.
+  //
+  // Inputs:
+  //   e  index into E of edge to try to collapse. E(e,:) = [s d] or [d s] so
+  //     that s<d, then d is collapsed to s.
+  ///  p  dim list of vertex position where to place merged vertex
+  // Inputs/Outputs:
+  //   V  #V by dim list of vertex positions, lesser index of E(e,:) will be set
+  //     to midpoint of edge.
+  //   F  #F by 3 list of face indices into V.
+  //   E  #E by 2 list of edge indices into V.
+  //   EMAP #F*3 list of indices into E, mapping each directed edge to unique
+  //     unique edge in E
+  //   EF  #E by 2 list of edge flaps, EF(e,0)=f means e=(i-->j) is the edge of
+  //     F(f,:) opposite the vth corner, where EI(e,0)=v. Similarly EF(e,1) "
+  //     e=(j->i)
+  //   EI  #E by 2 list of edge flap corners (see above).
+  //   e1  index into E of edge collpased on left
+  //   e2  index into E of edge collpased on left
+  //   f1  index into E of edge collpased on left
+  //   f2  index into E of edge collpased on left
+  // Returns true if edge was collapsed
+  #define IGL_COLLAPSE_EDGE_NULL 0
+  IGL_INLINE bool collapse_edge(
+    const int e,
+    const Eigen::RowVectorXd & p,
+    Eigen::MatrixXd & V,
+    Eigen::MatrixXi & F,
+    Eigen::MatrixXi & E,
+    Eigen::VectorXi & EMAP,
+    Eigen::MatrixXi & EF,
+    Eigen::MatrixXi & EI,
+    int & e1,
+    int & e2,
+    int & f1,
+    int & f2);
+  IGL_INLINE bool collapse_edge(
+    const int e,
+    const Eigen::RowVectorXd & p,
+    Eigen::MatrixXd & V,
+    Eigen::MatrixXi & F,
+    Eigen::MatrixXi & E,
+    Eigen::VectorXi & EMAP,
+    Eigen::MatrixXi & EF,
+    Eigen::MatrixXi & EI);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "collapse_edge.cpp"
+#endif
+#endif

include/igl/edge_collapse_is_valid.cpp

@@ -0,0 +1,87 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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 "edge_collapse_is_valid.h"
+#include "collapse_edge.h"
+#include "circulation.h"
+#include "intersect.h"
+#include "unique.h"
+#include "list_to_matrix.h"
+#include <vector>
+
+IGL_INLINE bool igl::edge_collapse_is_valid(
+  const int e,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXi & E,
+  const Eigen::VectorXi & EMAP,
+  const Eigen::MatrixXi & EF,
+  const Eigen::MatrixXi & EI)
+{
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  // For consistency with collapse_edge.cpp, let's determine edge flipness
+  // (though not needed to check validity)
+  const int eflip = E(e,0)>E(e,1);
+  // source and destination
+  const int s = eflip?E(e,1):E(e,0);
+  const int d = eflip?E(e,0):E(e,1);
+
+  if(s == IGL_COLLAPSE_EDGE_NULL && d==IGL_COLLAPSE_EDGE_NULL)
+  {
+    return false;
+  }
+  // check if edge collapse is valid: intersection of vertex neighbors of s and
+  // d should be exactly 2+(s,d) = 4
+  // http://stackoverflow.com/a/27049418/148668
+  {
+    // all vertex neighbors around edge, including the two vertices of the edge
+    const auto neighbors = [](
+      const int e,
+      const bool ccw,
+      const Eigen::MatrixXi & F,
+      const Eigen::MatrixXi & E,
+      const Eigen::VectorXi & EMAP,
+      const Eigen::MatrixXi & EF,
+      const Eigen::MatrixXi & EI) 
+    {
+      vector<int> N,uN;
+      vector<int> V2Fe = circulation(e, ccw,F,E,EMAP,EF,EI);
+      for(auto f : V2Fe)
+      {
+        N.push_back(F(f,0));
+        N.push_back(F(f,1));
+        N.push_back(F(f,2));
+      }
+      vector<size_t> _1,_2;
+      igl::unique(N,uN,_1,_2);
+      VectorXi uNm;
+      list_to_matrix(uN,uNm);
+      return uNm;
+    };
+    VectorXi Ns = neighbors(e, eflip,F,E,EMAP,EF,EI);
+    VectorXi Nd = neighbors(e,!eflip,F,E,EMAP,EF,EI);
+    VectorXi Nint = igl::intersect(Ns,Nd);
+    if(Nint.size() != 4)
+    {
+      return false;
+    }
+    if(Ns.size() == 4 && Nd.size() == 4)
+    {
+      VectorXi NsNd(8);
+      NsNd<<Ns,Nd;
+      VectorXi Nun,_1,_2;
+      igl::unique(NsNd,Nun,_1,_2);
+      // single tet, don't collapse
+      if(Nun.size() == 4)
+      {
+        return false;
+      }
+    }
+  }
+  return true;
+}

include/igl/edge_collapse_is_valid.h

@@ -0,0 +1,44 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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_EDGE_COLLAPSE_IS_VALID_H
+#define IGL_EDGE_COLLAPSE_IS_VALID_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Assumes (V,F) is a closed manifold mesh (except for previouslly collapsed
+  // faces which should be set to: 
+  // [IGL_COLLAPSE_EDGE_NULL IGL_COLLAPSE_EDGE_NULL IGL_COLLAPSE_EDGE_NULL].
+  // Tests whether collapsing exactly two faces and exactly 3 edges from E (e
+  // and one side of each face gets collapsed to the other) will result in a
+  // mesh with the same topology.
+  //
+  // Inputs:
+  //   e  index into E of edge to try to collapse. E(e,:) = [s d] or [d s] so
+  //     that s<d, then d is collapsed to s.
+  //   F  #F by 3 list of face indices into V.
+  //   E  #E by 2 list of edge indices into V.
+  //   EMAP #F*3 list of indices into E, mapping each directed edge to unique
+  //     unique edge in E
+  //   EF  #E by 2 list of edge flaps, EF(e,0)=f means e=(i-->j) is the edge of
+  //     F(f,:) opposite the vth corner, where EI(e,0)=v. Similarly EF(e,1) "
+  //     e=(j->i)
+  //   EI  #E by 2 list of edge flap corners (see above).
+  // Returns true if edge collapse is valid
+  IGL_INLINE bool edge_collapse_is_valid(
+    const int e,
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXi & E,
+    const Eigen::VectorXi & EMAP,
+    const Eigen::MatrixXi & EF,
+    const Eigen::MatrixXi & EI);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "edge_collapse_is_valid.cpp"
+#endif
+#endif

include/igl/edge_flaps.cpp

@@ -0,0 +1,59 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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 "edge_flaps.h"
+#include "unique_edge_map.h"
+#include <vector>
+#include <cassert>
+
+IGL_INLINE void igl::edge_flaps(
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXi & E,
+  const Eigen::VectorXi & EMAP,
+  Eigen::MatrixXi & EF,
+  Eigen::MatrixXi & EI)
+{
+  EF.resize(E.rows(),2);
+  EI.resize(E.rows(),2);
+  // loop over all faces
+  for(int f = 0;f<F.rows();f++)
+  {
+    // loop over edges across from corners
+    for(int v = 0;v<3;v++)
+    {
+      // get edge id
+      const int e = EMAP(v*F.rows()+f);
+      // See if this is left or right flap w.r.t. edge orientation
+      if( F(f,(v+1)%3) == E(e,0) && F(f,(v+2)%3) == E(e,1))
+      {
+        EF(e,0) = f;
+        EI(e,0) = v;
+      }else
+      {
+        assert(F(f,(v+1)%3) == E(e,1) && F(f,(v+2)%3) == E(e,0));
+        EF(e,1) = f;
+        EI(e,1) = v;
+      }
+    }
+  }
+}
+
+IGL_INLINE void igl::edge_flaps(
+  const Eigen::MatrixXi & F,
+  Eigen::MatrixXi & E,
+  Eigen::VectorXi & EMAP,
+  Eigen::MatrixXi & EF,
+  Eigen::MatrixXi & EI)
+{
+  Eigen::MatrixXi allE;
+  std::vector<std::vector<int> > uE2E;
+  igl::unique_edge_map(F,allE,E,EMAP,uE2E);
+  // Const-ify to call overload
+  const auto & cE = E;
+  const auto & cEMAP = EMAP;
+  return edge_flaps(F,cE,cEMAP,EF,EI);
+}

include/igl/edge_flaps.h

@@ -0,0 +1,45 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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_EDGE_FLAPS_H
+#define IGL_EDGE_FLAPS_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Determine "edge flaps": two faces on either side of a unique edge (assumes
+  // edge-manifold mesh)
+  //
+  // Inputs:
+  //   F  #F by 3 list of face indices
+  //   E  #E by 2 list of edge indices into V.
+  //   EMAP #F*3 list of indices into E, mapping each directed edge to unique
+  //     unique edge in E
+  // Outputs:
+  //   EF  #E by 2 list of edge flaps, EF(e,0)=f means e=(i-->j) is the edge of
+  //     F(f,:) opposite the vth corner, where EI(e,0)=v. Similarly EF(e,1) "
+  //     e=(j->i)
+  //   EI  #E by 2 list of edge flap corners (see above).
+  IGL_INLINE void edge_flaps(
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXi & E,
+    const Eigen::VectorXi & EMAP,
+    Eigen::MatrixXi & EF,
+    Eigen::MatrixXi & EI);
+  // Only faces as input
+  IGL_INLINE void edge_flaps(
+    const Eigen::MatrixXi & F,
+    Eigen::MatrixXi & E,
+    Eigen::VectorXi & EMAP,
+    Eigen::MatrixXi & EF,
+    Eigen::MatrixXi & EI);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "edge_flaps.cpp"
+#endif
+
+#endif

include/igl/linprog.cpp

@@ -0,0 +1,305 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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 "linprog.h"
+#include "slice.h"
+#include "slice_into.h"
+#include "find.h"
+#include "matlab_format.h"
+#include "colon.h"
+#include <iostream>
+
+//#define IGL_LINPROG_VERBOSE
+IGL_INLINE bool igl::linprog(
+  const Eigen::VectorXd & c,
+  const Eigen::MatrixXd & _A,
+  const Eigen::VectorXd & b,
+  const int k,
+  Eigen::VectorXd & x)
+{
+  // This is a very literal translation of
+  // http://www.mathworks.com/matlabcentral/fileexchange/2166-introduction-to-linear-algebra/content/strang/linprog.m
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  bool success = true;
+  // number of constraints
+  const int m = _A.rows();
+  // number of original variables
+  const int n = _A.cols();
+  // number of iterations
+  int it = 0;
+  // maximum number of iterations
+  //const int MAXIT = 10*m;
+  const int MAXIT = 100*m;
+  // residual tolerance
+  const double tol = 1e-10;
+  const auto & sign = [](const Eigen::VectorXd & B) -> Eigen::VectorXd
+  {
+    Eigen::VectorXd Bsign(B.size());
+    for(int i = 0;i<B.size();i++)
+    {
+      Bsign(i) = B(i)>0?1:(B(i)<0?-1:0);
+    }
+    return Bsign;
+  };
+  // initial (inverse) basis matrix
+  VectorXd Dv = sign(sign(b).array()+0.5);
+  Dv.head(k).setConstant(1.);
+  MatrixXd D = Dv.asDiagonal();
+  // Incorporate slack variables
+  MatrixXd A(_A.rows(),_A.cols()+D.cols());
+  A<<_A,D;
+  // Initial basis
+  VectorXi B = igl::colon<int>(n,n+m-1);
+  // non-basis, may turn out that vector<> would be better here
+  VectorXi N = igl::colon<int>(0,n-1);
+  int j;
+  double bmin = b.minCoeff(&j);
+  int phase;
+  VectorXd xb;
+  VectorXd s;
+  VectorXi J;
+  if(k>0 && bmin<0)
+  {
+    phase = 1;
+    xb = VectorXd::Ones(m);
+    // super cost
+    s.resize(n+m+1);
+    s<<VectorXd::Zero(n+k),VectorXd::Ones(m-k+1);
+    N.resize(n+1);
+    N<<igl::colon<int>(0,n-1),B(j);
+    J.resize(B.size()-1);
+    // [0 1 2 3 4]
+    //      ^
+    // [0 1]
+    //      [3 4]
+    J.head(j) = B.head(j);
+    J.tail(B.size()-j-1) = B.tail(B.size()-j-1);
+    B(j) = n+m;
+    MatrixXd AJ;
+    igl::slice(A,J,2,AJ);
+    const VectorXd a = b - AJ.rowwise().sum();
+    {
+      MatrixXd old_A = A;
+      A.resize(A.rows(),A.cols()+a.cols());
+      A<<old_A,a;
+    }
+    D.col(j) = -a/a(j);
+    D(j,j) = 1./a(j);
+  }else if(k==m)
+  {
+    phase = 2;
+    xb = b;
+    s.resize(c.size()+m);
+    // cost function
+    s<<c,VectorXd::Zero(m);
+  }else //k = 0 or bmin >=0
+  {
+    phase = 1;
+    xb = b.array().abs();
+    s.resize(n+m);
+    // super cost
+    s<<VectorXd::Zero(n+k),VectorXd::Ones(m-k);
+  }
+  while(phase<3)
+  {
+    double df = -1;
+    int t = std::numeric_limits<int>::max();
+    // Lagrange mutipliers fro Ax=b
+    VectorXd yb = D.transpose() * igl::slice(s,B);
+    while(true)
+    {
+      if(MAXIT>0 && it>=MAXIT)
+      {
+#ifdef IGL_LINPROG_VERBOSE
+        cerr<<"linprog: warning! maximum iterations without convergence."<<endl;
+#endif
+        success = false;
+        break;
+      }
+      // no freedom for minimization
+      if(N.size() == 0)
+      {
+        break;
+      }
+      // reduced costs
+      VectorXd sN = igl::slice(s,N);
+      MatrixXd AN = igl::slice(A,N,2);
+      VectorXd r = sN - AN.transpose() * yb;
+      int q;
+      // determine new basic variable
+      double rmin = r.minCoeff(&q);
+      // optimal! infinity norm
+      if(rmin>=-tol*(sN.array().abs().maxCoeff()+1))
+      {
+        break;
+      }
+      // increment iteration count
+      it++;
+      // apply Bland's rule to avoid cycling
+      if(df>=0)
+      {
+        if(MAXIT == -1)
+        {
+#ifdef IGL_LINPROG_VERBOSE
+          cerr<<"linprog: warning! degenerate vertex"<<endl;
+#endif
+          success = false;
+        }
+        igl::find((r.array()<0).eval(),J);
+        double Nq = igl::slice(N,J).minCoeff();
+        // again seems like q is assumed to be a scalar though matlab code
+        // could produce a vector for multiple matches
+        (N.array()==Nq).cast<int>().maxCoeff(&q);
+      }
+      VectorXd d = D*A.col(N(q));
+      VectorXi I;
+      igl::find((d.array()>tol).eval(),I);
+      if(I.size() == 0)
+      {
+#ifdef IGL_LINPROG_VERBOSE
+        cerr<<"linprog: warning! solution is unbounded"<<endl;
+#endif
+        // This seems dubious:
+        it=-it;
+        success = false;
+        break;
+      }
+      VectorXd xbd = igl::slice(xb,I).array()/igl::slice(d,I).array();
+      // new use of r
+      int p;
+      {
+        double r;
+        r = xbd.minCoeff(&p);
+        p = I(p);
+        // apply Bland's rule to avoid cycling
+        if(df>=0)
+        {
+          igl::find((xbd.array()==r).eval(),J);
+          double Bp = igl::slice(B,igl::slice(I,J)).minCoeff();
+          // idiotic way of finding index in B of Bp
+          // code down the line seems to assume p is a scalar though the matlab
+          // code could find a vector of matches)
+          (B.array()==Bp).cast<int>().maxCoeff(&p);
+        }
+        // update x
+        xb -= r*d;
+        xb(p) = r;
+        // change in f
+        df = r*rmin;
+      }
+      // row vector
+      RowVectorXd v = D.row(p)/d(p);
+      yb += v.transpose() * (s(N(q)) - d.transpose()*igl::slice(s,B));
+      d(p)-=1;
+      // update inverse basis matrix
+      D = D - d*v;
+      t = B(p);
+      B(p) = N(q);
+      if(t>(n+k-1))
+      {
+        // remove qth entry from N
+        VectorXi old_N = N;
+        N.resize(N.size()-1);
+        N.head(q) = old_N.head(q);
+        N.head(q) = old_N.head(q);
+        N.tail(old_N.size()-q-1) = old_N.tail(old_N.size()-q-1);
+      }else
+      {
+        N(q) = t;
+      }
+    }
+    // iterative refinement
+    xb = (xb+D*(b-igl::slice(A,B,2)*xb)).eval();
+    // must be due to rounding
+    VectorXi I;
+    igl::find((xb.array()<0).eval(),I);
+    if(I.size()>0)
+    {
+      // so correct
+      VectorXd Z = VectorXd::Zero(I.size(),1);
+      igl::slice_into(Z,I,xb);
+    }
+    // B, xb,n,m,res=A(:,B)*xb-b
+    if(phase == 2 || it<0)
+    {
+      break;
+    }
+    if(xb.transpose()*igl::slice(s,B) > tol)
+    {
+      it = -it;
+#ifdef IGL_LINPROG_VERBOSE
+      cerr<<"linprog: warning, no feasible solution"<<endl;
+#endif
+      success = false;
+      break;
+    }
+    // re-initialize for Phase 2
+    phase = phase+1;
+    s*=1e6*c.array().abs().maxCoeff();
+    s.head(n) = c;
+  }
+  x.resize(std::max(B.maxCoeff()+1,n));
+  igl::slice_into(xb,B,x);
+  x = x.head(n).eval();
+  return success;
+}
+
+IGL_INLINE bool igl::linprog(
+  const Eigen::VectorXd & f,
+  const Eigen::MatrixXd & A,
+  const Eigen::VectorXd & b,
+  const Eigen::MatrixXd & B,
+  const Eigen::VectorXd & c,
+  Eigen::VectorXd & x)
+{
+  using namespace Eigen;
+  using namespace igl;
+  using namespace std;
+  const int m = A.rows();
+  const int n = A.cols();
+  const int p = B.rows();
+  MatrixXd Im = MatrixXd::Identity(m,m);
+  MatrixXd AS(m,n+m);
+  AS<<A,Im;
+  MatrixXd bS = b.array().abs();
+  for(int i = 0;i<m;i++)
+  {
+    const auto & sign = [](double x)->double
+    {
+      return (x<0?-1:(x>0?1:0));
+    };
+    AS.row(i) *= sign(b(i));
+  }
+  MatrixXd In = MatrixXd::Identity(n,n);
+  MatrixXd P(n+m,2*n+m);
+  P<<              In, -In, MatrixXd::Zero(n,m),
+     MatrixXd::Zero(m,2*n), Im;
+  MatrixXd ASP = AS*P;
+  MatrixXd BSP(0,2*n+m);
+  if(p>0)
+  {
+    MatrixXd BS(p,2*n);
+    BS<<B,MatrixXd::Zero(p,n);
+    BSP = BS*P;
+  }
+
+  VectorXd fSP = VectorXd::Ones(2*n+m);
+  fSP.head(2*n) = P.block(0,0,n,2*n).transpose()*f;
+  const VectorXd & cc = fSP;
+
+  MatrixXd AA(m+p,2*n+m);
+  AA<<ASP,BSP;
+  VectorXd bb(m+p);
+  bb<<bS,c;
+
+  VectorXd xxs;
+  bool ret = linprog(cc,AA,bb,0,xxs);
+  x = P.block(0,0,n,2*n+m)*xxs;
+  return ret;
+}

include/igl/linprog.h

@@ -0,0 +1,65 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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_LINPROG_H
+#define IGL_LINPROG_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Solve a linear program given in "standard form"
+  //
+  // min  f'x
+  // s.t. A(    1:k,:) x <= b(1:k)
+  //      A(k+1:end,:) x = b(k+1:end)
+  //   ** x >= 0 **
+  //
+  // In contrast to other APIs the entries in b may be negative.
+  //
+  // Inputs:
+  //   c  #x list of linear coefficients
+  //   A  #A by #x matrix of linear constraint coefficients
+  //   b  #A list of linear constraint right-hand sides
+  //   k  number of inequality constraints as first rows of A,b
+  // Outputs:
+  //   x  #x solution vector
+  //
+  IGL_INLINE bool linprog(
+    const Eigen::VectorXd & c,
+    const Eigen::MatrixXd & A,
+    const Eigen::VectorXd & b,
+    const int k,
+    Eigen::VectorXd & f);
+  
+  // Wrapper in friendlier general form (no implicit bounds on x)
+  //
+  // min  f'x
+  // s.t. A x <= b
+  //      B x = c
+  //
+  // Inputs:
+  //   f  #x list of linear coefficients
+  //   A  #A by #x matrix of linear inequality constraint coefficients
+  //   b  #A list of linear constraint right-hand sides
+  //   B  #B by #x matrix of linear equality constraint coefficients
+  //   c  #B list of linear constraint right-hand sides
+  // Outputs:
+  //   x  #x solution vector
+  //
+  IGL_INLINE bool linprog(
+    const Eigen::VectorXd & f,
+    const Eigen::MatrixXd & A,
+    const Eigen::VectorXd & b,
+    const Eigen::MatrixXd & B,
+    const Eigen::VectorXd & c,
+    Eigen::VectorXd & x);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "linprog.cpp"
+#endif
+#endif

include/igl/mosek/mosek_linprog.cpp

@@ -0,0 +1,162 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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 "mosek_linprog.h"
+#include "../mosek/mosek_guarded.h"
+#include "../harwell_boeing.h"
+#include <limits>
+#include <cmath>
+#include <vector>
+
+IGL_INLINE bool igl::mosek_linprog(
+  const Eigen::VectorXd & c,
+  const Eigen::SparseMatrix<double> & A,
+  const Eigen::VectorXd & lc,
+  const Eigen::VectorXd & uc,
+  const Eigen::VectorXd & lx,
+  const Eigen::VectorXd & ux,
+  Eigen::VectorXd & x)
+{
+  // variables for mosek task, env and result code
+  MSKenv_t env;
+  // Create the MOSEK environment
+  mosek_guarded(MSK_makeenv(&env,NULL));
+  // initialize mosek environment
+  mosek_guarded(MSK_initenv(env));
+  const bool ret = mosek_linprog(c,A,lc,uc,lx,ux,env,x);
+  MSK_deleteenv(&env);
+  return ret;
+}
+
+IGL_INLINE bool igl::mosek_linprog(
+  const Eigen::VectorXd & c,
+  const Eigen::SparseMatrix<double> & A,
+  const Eigen::VectorXd & lc,
+  const Eigen::VectorXd & uc,
+  const Eigen::VectorXd & lx,
+  const Eigen::VectorXd & ux,
+  const MSKenv_t & env,
+  Eigen::VectorXd & x)
+{
+  // following http://docs.mosek.com/7.1/capi/Linear_optimization.html
+  using namespace std;
+  // number of constraints
+  const int m = A.rows();
+  // number of variables
+  const int n = A.cols();
+
+
+  vector<double> vAv;
+  vector<int> vAri,vAcp;
+  int nr;
+  harwell_boeing(A,nr,vAv,vAri,vAcp);
+
+  MSKtask_t task;
+  // Create the optimization task
+  mosek_guarded(MSK_maketask(env,m,n,&task));
+  // no threads
+  mosek_guarded(MSK_putintparam(task,MSK_IPAR_NUM_THREADS,1));
+  if(m>0)
+  {
+    // Append 'm' empty constraints, the constrainst will initially have no
+    // bounds
+    mosek_guarded(MSK_appendcons(task,m));
+  }
+  mosek_guarded(MSK_appendvars(task,n));
+
+  
+  const auto & key = [](const double lxj, const double uxj) ->
+    MSKboundkeye
+  {
+    MSKboundkeye k = MSK_BK_FR;
+    if(isfinite(lxj) && isfinite(uxj))
+    {
+      if(lxj == uxj)
+      {
+        k = MSK_BK_RA;
+      }else{
+        k = MSK_BK_FX;
+      }
+    }else if(isfinite(lxj))
+    {
+      k = MSK_BK_LO;
+    }else if(isfinite(uxj))
+    {
+      k = MSK_BK_UP;
+    }
+    return k;
+  };
+
+  // loop over variables
+  for(int j = 0;j<n;j++)
+  {
+    if(c.size() > 0)
+    {
+      // Set linear term c_j in the objective
+      mosek_guarded(MSK_putcj(task,j,c(j)));
+    }
+
+    // Set constant bounds on variable j
+    const double lxj = lx.size()>0?lx[j]:-numeric_limits<double>::infinity();
+    const double uxj = ux.size()>0?ux[j]: numeric_limits<double>::infinity();
+    mosek_guarded(MSK_putvarbound(task,j,key(lxj,uxj),lxj,uxj));
+
+    if(m>0)
+    {
+      // Input column j of A
+      mosek_guarded(
+        MSK_putacol(
+          task,
+          j,
+          vAcp[j+1]-vAcp[j],
+          &vAri[vAcp[j]],
+          &vAv[vAcp[j]])
+        );
+    }
+  }
+  // loop over constraints
+  for(int i = 0;i<m;i++)
+  {
+    // Set constraint bounds for row i
+    const double lci = lc.size()>0?lc[i]:-numeric_limits<double>::infinity();
+    const double uci = uc.size()>0?uc[i]: numeric_limits<double>::infinity();
+    mosek_guarded(MSK_putconbound(task,i,key(lci,uci),lci,uci));
+  }
+
+  // Now the optimizer has been prepared
+  MSKrescodee trmcode;
+  // run the optimizer
+  mosek_guarded(MSK_optimizetrm(task,&trmcode));
+  // Get status
+  MSKsolstae solsta;
+  MSK_getsolsta (task,MSK_SOL_ITR,&solsta);
+  bool success = false;
+  switch(solsta)
+  {
+    case MSK_SOL_STA_OPTIMAL:   
+    case MSK_SOL_STA_NEAR_OPTIMAL:
+      x.resize(n);
+      /* Request the basic solution. */ 
+      MSK_getxx(task,MSK_SOL_BAS,x.data()); 
+      success = true;
+      break;
+    case MSK_SOL_STA_DUAL_INFEAS_CER:
+    case MSK_SOL_STA_PRIM_INFEAS_CER:
+    case MSK_SOL_STA_NEAR_DUAL_INFEAS_CER:
+    case MSK_SOL_STA_NEAR_PRIM_INFEAS_CER:  
+      //printf("Primal or dual infeasibility certificate found.\n");
+      break;
+    case MSK_SOL_STA_UNKNOWN:
+      //printf("The status of the solution could not be determined.\n");
+      break;
+    default:
+      //printf("Other solution status.");
+      break;
+  }
+  MSK_deletetask(&task);
+  return success;
+}

include/igl/mosek/mosek_linprog.h

@@ -0,0 +1,56 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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_MOSEK_LINPROG_H
+#define IGL_MOSEK_LINPROG_H
+#include "../igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+#include <mosek.h>
+namespace igl
+{
+  // Solve a linear program using mosek:
+  // 
+  // min c'x
+  // s.t. lc <= A x <= uc
+  //      lx <= x <= ux
+  //
+  // Inputs:
+  //   c  #x list of linear objective coefficients
+  //   A  #A by #x matrix of linear inequality constraint coefficients
+  //   lc  #A list of lower constraint bounds
+  //   uc  #A list of upper constraint bounds
+  //   lx  #x list of lower variable bounds
+  //   ux  #x list of upper variable bounds
+  // Outputs:
+  //   x  #x list of solution values
+  // Returns true iff success.
+  IGL_INLINE bool mosek_linprog(
+    const Eigen::VectorXd & c,
+    const Eigen::SparseMatrix<double> & A,
+    const Eigen::VectorXd & lc,
+    const Eigen::VectorXd & uc,
+    const Eigen::VectorXd & lx,
+    const Eigen::VectorXd & ux,
+    Eigen::VectorXd & x);
+  // Wrapper that keeps mosek environment alive (if licence checking is
+  // becoming a bottleneck)
+  IGL_INLINE bool mosek_linprog(
+    const Eigen::VectorXd & c,
+    const Eigen::SparseMatrix<double> & A,
+    const Eigen::VectorXd & lc,
+    const Eigen::VectorXd & uc,
+    const Eigen::VectorXd & lx,
+    const Eigen::VectorXd & ux,
+    const MSKenv_t & env,
+    Eigen::VectorXd & x);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "mosek_linprog.cpp"
+#endif
+#endif