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

Former-commit-id: cc67ceaa330dca53aea9faaefce074ffd9db6dc8

README.md

@@ -107,13 +107,13 @@ For more information see our [tutorial](tutorial/tutorial.html).
 Libigl compartmentalizes its **optional** dependences via its directory
 organization in the `include/` folder. All header files located _directly_ in
 the `include/igl/` folder have only stl and Eigen as dependencies. For example,
-all of the headers that depend on CGAL are located in `include/igl/cgal`. For a
-full list of _optional_ dependencies check `optional/CMakeLists.txt`.
+all of the headers that depend on CGAL are located in `include/igl/copyleft/cgal`.
+For a full list of _optional_ dependencies check `optional/CMakeLists.txt`.
 
 ### GCC and the Optional CGAL Dependency
-The `include/igl/cgal/*.h` headers depend on CGAL. It has come to our attention
-that CGAL does not work properly with GCC 4.8. To the best of our knowledge,
-GCC 4.7 and clang will work correctly.
+The `include/igl/copyleft/cgal/*.h` headers depend on CGAL. It has come to
+our attention that CGAL does not work properly with GCC 4.8. To the best of
+our knowledge, GCC 4.7 and clang will work correctly.
 
 ### OpenMP and Windows
 Some of our functions will take advantage of OpenMP if available. However, it

include/igl/hessian.cpp

@@ -0,0 +1,60 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Alec Jacobson <alecjacobson@gmail.com>
+//  and Oded Stein <oded.stein@columbia.edu>
+//
+// 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 "hessian.h"
+#include <vector>
+
+#include "grad.h"
+#include "igl/doublearea.h"
+#include "igl/repdiag.h"
+
+
+
+template <typename DerivedV, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::hessian(
+                             const Eigen::MatrixBase<DerivedV> & V,
+                             const Eigen::MatrixBase<DerivedF> & F,
+                             Eigen::SparseMatrix<Scalar>& H)
+{
+    typedef typename DerivedV::Scalar denseScalar;
+    typedef typename Eigen::Matrix<denseScalar, Eigen::Dynamic, 1> VecXd;
+    typedef typename Eigen::SparseMatrix<Scalar> SparseMat;
+    typedef typename Eigen::DiagonalMatrix
+                       <Scalar, Eigen::Dynamic, Eigen::Dynamic> DiagMat;
+    
+    int dim = V.cols();
+    assert((dim==2 || dim==3) &&
+           "The dimension of the vertices should be 2 or 3");
+    
+    //Construct the combined gradient matric
+    SparseMat G;
+    igl::grad(Eigen::PlainObjectBase<DerivedV>(V),
+              Eigen::PlainObjectBase<DerivedF>(F),
+              G, false);
+    SparseMat GG(F.rows(), dim*V.rows());
+    GG.reserve(G.nonZeros());
+    for(int i=0; i<dim; ++i)
+        GG.middleCols(i*G.cols(),G.cols()) = G.middleRows(i*F.rows(),F.rows());
+    SparseMat D;
+    igl::repdiag(GG,dim,D);
+    
+    //Compute area matrix
+    VecXd areas;
+    igl::doublearea(V, F, areas);
+    DiagMat A = (0.5*areas).replicate(dim,1).asDiagonal();
+    
+    //Compute FEM Hessian
+    H = D.transpose()*A*G;
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::hessian<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/hessian.h

@@ -0,0 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Alec Jacobson <alecjacobson@gmail.com>
+//  and Oded Stein <oded.stein@columbia.edu>
+//
+// 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_FEM_HESSIAN_H
+#define IGL_FEM_HESSIAN_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+
+namespace igl
+{
+    // Constructs the finite element Hessian matrix
+    // as described in https://arxiv.org/abs/1707.04348,
+    // Natural Boundary Conditions for Smoothing in Geometry Processing
+    // (Oded Stein, Eitan Grinspun, Max Wardetzky, Alec Jacobson)
+    // The interior vertices are NOT set to zero yet.
+    //
+    // Inputs:
+    //   V  #V by dim list of mesh vertex positions
+    //   F  #F by 3 list of mesh faces (must be triangles)
+    // Outputs:
+    //   H  #V by #V Hessian energy matrix, each column i
+    //     corresponding to V(i,:)
+    //
+    //
+    //
+    template <typename DerivedV, typename DerivedF, typename Scalar>
+    IGL_INLINE void hessian(
+                            const Eigen::MatrixBase<DerivedV> & V,
+                            const Eigen::MatrixBase<DerivedF> & F,
+                            Eigen::SparseMatrix<Scalar>& H);
+    
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "hessian.cpp"
+#endif
+
+#endif

include/igl/hessian_energy.cpp

@@ -0,0 +1,64 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Alec Jacobson <alecjacobson@gmail.com>
+//   and Oded Stein <oded.stein@columbia.edu>
+//
+// 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 "hessian_energy.h"
+#include <vector>
+
+#include "hessian.h"
+#include "massmatrix.h"
+#include "boundary_loop.h"
+
+
+template <typename DerivedV, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::hessian_energy(
+                                    const Eigen::MatrixBase<DerivedV> & V,
+                                    const Eigen::MatrixBase<DerivedF> & F,
+                                    Eigen::SparseMatrix<Scalar>& Q)
+{
+    typedef typename DerivedV::Scalar denseScalar;
+    typedef typename Eigen::Matrix<denseScalar, Eigen::Dynamic, 1> VecXd;
+    typedef typename Eigen::SparseMatrix<Scalar> SparseMat;
+    typedef typename Eigen::DiagonalMatrix
+                       <Scalar, Eigen::Dynamic, Eigen::Dynamic> DiagMat;
+    
+    int dim = V.cols();
+    assert((dim==2 || dim==3) &&
+           "The dimension of the vertices should be 2 or 3");
+    
+    SparseMat M;
+    igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_VORONOI,M);
+    
+    //Kill non-interior DOFs
+    VecXd Mint = M.diagonal();
+    std::vector<std::vector<int> > bdryLoop;
+    igl::boundary_loop(Eigen::PlainObjectBase<DerivedF>(F),bdryLoop);
+    for(const std::vector<int>& loop : bdryLoop)
+        for(const int& bdryVert : loop)
+            Mint(bdryVert) = 0.;
+    
+    //Invert Mint
+    for(int i=0; i<Mint.rows(); ++i)
+        if(Mint(i) > 0)
+            Mint(i) = 1./Mint(i);
+    
+    //Repeat Mint to form diaginal matrix
+    DiagMat stackedMinv = Mint.replicate(dim*dim,1).asDiagonal();
+    
+    //Compute squared Hessian
+    SparseMat H;
+    igl::hessian(V,F,H);
+    Q = H.transpose()*stackedMinv*H;
+    
+}
+
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::hessian_energy<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/hessian_energy.h

@@ -0,0 +1,45 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Alec Jacobson <alecjacobson@gmail.com>
+//  and Oded Stein <oded.stein@columbia.edu>
+//
+// 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_HESSIAN_ENERGY_H
+#define IGL_HESSIAN_ENERGY_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+
+namespace igl
+{
+    // Constructs the Hessian energy matrix using mixed FEM
+    // as described in https://arxiv.org/abs/1707.04348
+    // Natural Boundary Conditions for Smoothing in Geometry Processing
+    // (Oded Stein, Eitan Grinspun, Max Wardetzky, Alec Jacobson)
+    //
+    // Inputs:
+    //   V  #V by dim list of mesh vertex positions
+    //   F  #F by 3 list of mesh faces (must be triangles)
+    // Outputs:
+    //   Q  #V by #V Hessian energy matrix, each row/column i
+    //     corresponding to V(i,:)
+    //
+    //
+    //
+    template <typename DerivedV, typename DerivedF, typename Scalar>
+    IGL_INLINE void hessian_energy(
+                                   const Eigen::MatrixBase<DerivedV> & V,
+                                   const Eigen::MatrixBase<DerivedF> & F,
+                                   Eigen::SparseMatrix<Scalar>& Q);
+    
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "hessian_energy.cpp"
+#endif
+
+#endif

tutorial/712_DataSmoothing/CMakeLists.txt

@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 2.8.12)
+project(712_DataSmoothing)
+
+add_executable(${PROJECT_NAME}_bin
+  main.cpp)
+target_include_directories(${PROJECT_NAME}_bin PRIVATE ${LIBIGL_INCLUDE_DIRS})
+target_compile_definitions(${PROJECT_NAME}_bin PRIVATE ${LIBIGL_DEFINITIONS})
+target_link_libraries(${PROJECT_NAME}_bin ${LIBIGL_LIBRARIES} ${LIBIGL_VIEWER_EXTRA_LIBRARIES} ${LIBIGL_OPENGL_EXTRA_LIBRARIES} ${LIBIGL_OPENGL_GLFW_EXTRA_LIBRARIES})

tutorial/712_DataSmoothing/isolines.h

@@ -0,0 +1,90 @@
+
+#include <vector>
+#include <limits>
+#include <stdlib.h>
+
+#include <igl/remove_unreferenced.h>
+
+
+static void isolines(const Eigen::MatrixXd& V,
+    const Eigen::MatrixXi& F,
+    const Eigen::VectorXd& z,
+    const int grads,
+    Eigen::MatrixXd& isoV,
+    Eigen::MatrixXi& isoE)
+{
+    const double min = z.minCoeff(), max = z.maxCoeff();
+
+    //Following http://www.alecjacobson.com/weblog/?p=2529
+    Eigen::VectorXd iso(grads+1);
+    for(int i=0; i<iso.size(); ++i)
+        iso(i) = double(i)/double(grads)*(max-min) + min;
+
+    Eigen::MatrixXd t12(F.rows(), iso.size()), t23(F.rows(), iso.size()),
+        t31(F.rows(), iso.size());
+    for(int i=0; i<t12.rows(); ++i) {
+        const double z1=z(F(i,0)), z2=z(F(i,1)), z3=z(F(i,2));
+        const double s12 = z2-z1;
+        const double s23 = z3-z2;
+        const double s31 = z1-z3;
+        for(int j=0; j<t12.cols(); ++j) {
+            t12(i,j) = (iso(j)-z1) / s12;
+            t23(i,j) = (iso(j)-z2) / s23;
+            t31(i,j) = (iso(j)-z3) / s31;
+            if(t12(i,j)<0 || t12(i,j)>1)
+                t12(i,j) = std::numeric_limits<double>::quiet_NaN();
+            if(t23(i,j)<0 || t23(i,j)>1)
+                t23(i,j) = std::numeric_limits<double>::quiet_NaN();
+            if(t31(i,j)<0 || t31(i,j)>1)
+                t31(i,j) = std::numeric_limits<double>::quiet_NaN();
+        }
+    }
+
+    std::vector<int> F12, F23, F31, I12, I23, I31;
+    for(int i=0; i<t12.rows(); ++i) {
+        for(int j=0; j<t12.cols(); ++j) {
+            if(std::isfinite(t23(i,j)) && std::isfinite(t31(i,j))) {
+                F12.push_back(i);
+                I12.push_back(j);
+            }
+            if(std::isfinite(t31(i,j)) && std::isfinite(t12(i,j))) {
+                F23.push_back(i);
+                I23.push_back(j);
+            }
+            if(std::isfinite(t12(i,j)) && std::isfinite(t23(i,j))) {
+                F31.push_back(i);
+                I31.push_back(j);
+            }
+        }
+    }
+
+    const int K = F12.size()+F23.size()+F31.size();
+    isoV.resize(2*K, 3);
+    int b = 0;
+    for(int i=0; i<F12.size(); ++i) {
+        isoV.row(b+i) = (1.-t23(F12[i],I12[i]))*V.row(F(F12[i],1))
+            + t23(F12[i],I12[i])*V.row(F(F12[i],2));
+        isoV.row(K+b+i) = (1.-t31(F12[i],I12[i]))*V.row(F(F12[i],2))
+            + t31(F12[i],I12[i])*V.row(F(F12[i],0));
+    }
+    b += F12.size();
+    for(int i=0; i<F23.size(); ++i) {
+        isoV.row(b+i) = (1.-t31(F23[i],I23[i]))*V.row(F(F23[i],2))
+            + t31(F23[i],I23[i])*V.row(F(F23[i],0));
+        isoV.row(K+b+i) = (1.-t12(F23[i],I23[i]))*V.row(F(F23[i],0))
+            + t12(F23[i],I23[i])*V.row(F(F23[i],1));
+    }
+    b += F23.size();
+    for(int i=0; i<F31.size(); ++i) {
+        isoV.row(b+i) = (1.-t12(F31[i],I31[i]))*V.row(F(F31[i],0))
+            + t12(F31[i],I31[i])*V.row(F(F31[i],1));
+        isoV.row(K+b+i) = (1.-t23(F31[i],I31[i]))*V.row(F(F31[i],1))
+            + t23(F31[i],I31[i])*V.row(F(F31[i],2));
+    }
+
+    isoE.resize(K,2);
+    for(int i=0; i<K; ++i)
+        isoE.row(i) << i, K+i;
+
+}
+

tutorial/712_DataSmoothing/main.cpp

@@ -0,0 +1,109 @@
+#include <igl/read_triangle_mesh.h>
+#include <igl/hessian_energy.h>
+#include <igl/massmatrix.h>
+#include <igl/cotmatrix.h>
+#include <igl/jet.h>
+#include <igl/edges.h>
+#include <igl/components.h>
+#include <igl/remove_unreferenced.h>
+#include <igl/viewer/Viewer.h>
+
+#include <Eigen/Core>
+#include <Eigen/SparseCholesky>
+
+#include <iostream>
+#include <set>
+#include <limits>
+#include <stdlib.h>
+
+#include "tutorial_shared_path.h"
+
+#include "isolines.h"
+
+
+int main(int argc, char * argv[])
+{
+    typedef Eigen::SparseMatrix<double> SparseMat;
+
+    //Read our mesh
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F, E;
+    if(!igl::read_triangle_mesh(
+        argc>1?argv[1]: TUTORIAL_SHARED_PATH "/beetle.off",V,F)) {
+        std::cout << "Failed to load mesh." << std::endl;
+    }
+    igl::edges(F,E);
+
+    //Constructing an exact function to smooth
+    Eigen::VectorXd zexact = V.block(0,2,V.rows(),1).array()
+        + 0.5*V.block(0,1,V.rows(),1).array()
+        + V.block(0,1,V.rows(),1).array().pow(2)
+        + V.block(0,2,V.rows(),1).array().pow(3);
+
+    //Make the exact function noisy
+    srand(5);
+    const double s = 0.2*(zexact.maxCoeff() - zexact.minCoeff());
+    Eigen::VectorXd znoisy = zexact + s*Eigen::VectorXd::Random(zexact.size());
+
+    //Constructing the squared Laplacian and squared Hessian energy
+    SparseMat L, M;
+    igl::cotmatrix(V, F, L);
+    igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_BARYCENTRIC, M);
+    Eigen::SimplicialLDLT<SparseMat> solver(M);
+    SparseMat MinvL = solver.solve(L);
+    SparseMat QL = L.transpose()*MinvL;
+    SparseMat QH;
+    igl::hessian_energy(V, F, QH);
+
+    //Solve to find Laplacian-smoothed and Hessian-smoothed solutions
+    const double al = 8e-4;
+    Eigen::SimplicialLDLT<SparseMat> lapSolver(al*QL + (1.-al)*M);
+    Eigen::VectorXd zl = lapSolver.solve(al*M*znoisy);
+    const double ah = 5e-6;
+    Eigen::SimplicialLDLT<SparseMat> hessSolver(ah*QH + (1.-ah)*M);
+    Eigen::VectorXd zh = hessSolver.solve(ah*M*znoisy);
+
+    //Viewer that shows all functions: zexact, znoisy, zl, zh
+    igl::viewer::Viewer viewer;
+    viewer.data.set_mesh(V,F);
+    viewer.core.show_lines = false;
+    viewer.callback_key_down =
+      [&](igl::viewer::Viewer & viewer, unsigned char key, int mod)->bool {
+        //Graduate result to show isolines, then compute color matrix
+        const Eigen::VectorXd* z;
+        switch(key) {
+            case '1':
+                z = &zexact;
+                break;
+            case '2':
+                z = &znoisy;
+                break;
+            case '3':
+                z = &zl;
+                break;
+            case '4':
+                z = &zh;
+                break;
+            default:
+                return false;
+        }
+        Eigen::MatrixXd isoV;
+        Eigen::MatrixXi isoE;
+        if(key!='2')
+            isolines(V, F, *z, 30, isoV, isoE);
+        viewer.data.set_edges(isoV,isoE,Eigen::RowVector3d(0,0,0));
+        Eigen::MatrixXd colors;
+        igl::jet(*z, true, colors);
+        viewer.data.set_colors(colors);
+    };
+    std::cout << R"(Usage:
+1  Show original function
+2  Show noisy function
+3  Biharmonic smoothing (zero Neumann boundary)
+4  Biharmonic smoothing (natural Hessian boundary)
+
+)";
+    viewer.launch();
+
+    return 0;
+}

tutorial/CMakeLists.txt

@@ -186,4 +186,5 @@ if(TUTORIALS_CHAPTER7)
   endif()
   add_subdirectory("710_SLIM")
   add_subdirectory("711_Subdivision")
+  add_subdirectory("712_DataSmoothing")
 endif()

tutorial/images/712_beetles.jpg.REMOVED.git-id

@@ -0,0 +1 @@
+5513240871aa3e7a212ac8c2ae8cb466ffac3793

tutorial/tutorial.html.REMOVED.git-id

@@ -1 +1 @@
-5ae3f8bf55895f273b2e1417314d6ed544096ea6
+e0514ec101376a87462a079591da84d17750f2fa

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-4868a892c78947f78ff6985ff7ea6cc069dd1b4d
+58d89e2b385d7c8849f9e3d7d1ec66e970cad92d