Adding biharmonic energy function

CMakeLists.txt

@@ -98,6 +98,7 @@ pybind11_add_module(gpytoolbox_bindings
 	"${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_upper_envelope.cpp"
 	"${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_read_ply.cpp"
 	"${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_write_ply.cpp"
+	"${CMAKE_CURRENT_SOURCE_DIR}/src/cpp/binding_curved_hessian_intrinsic.cpp"
 )
 
 pybind11_add_module(gpytoolbox_bindings_copyleft

src/cpp/binding_curved_hessian_intrinsic.cpp

@@ -0,0 +1,26 @@
+#include <igl/curved_hessian_energy.h>
+#include <igl/orient_halfedges.h>
+#include <pybind11/stl.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/eigen.h>
+#include <pybind11/functional.h>
+#include <string>
+
+using namespace Eigen;
+namespace py = pybind11;
+using EigenDStride = Stride<Eigen::Dynamic, Eigen::Dynamic>;
+template <typename MatrixType>
+using EigenDRef = Ref<MatrixType, 0, EigenDStride>; //allows passing column/row order matrices easily
+
+void binding_curved_hessian_intrinsic(py::module& m) {
+    m.def("_curved_hessian_intrinsic_cpp_impl",[](EigenDRef<MatrixXd> l_sq,
+                         EigenDRef<MatrixXi> f)
+        {
+            Eigen::MatrixXi E, oE;
+            igl::orient_halfedges(f, E, oE);
+            Eigen::SparseMatrix<double> Q;
+            igl::curved_hessian_energy_intrinsic(f, l_sq, E, oE, Q);
+            return Q;
+        });
+
+}

src/cpp/gpytoolbox_bindings_core.cpp

@@ -26,6 +26,7 @@ void binding_remesh_botsch(py::module& m);
 void binding_upper_envelope(py::module& m);
 void binding_read_ply(py::module& m);
 void binding_write_ply(py::module& m);
+void binding_curved_hessian_intrinsic(py::module& m);
 
 PYBIND11_MODULE(gpytoolbox_bindings, m) {
 
@@ -46,6 +47,7 @@ PYBIND11_MODULE(gpytoolbox_bindings, m) {
     binding_upper_envelope(m);
     binding_read_ply(m);
     binding_write_ply(m);
+    binding_curved_hessian_intrinsic(m);
 
     m.def("help", [&]() {printf("hi"); });
 }

src/gpytoolbox/__init__.py

@@ -29,6 +29,7 @@
 from .quadtree_laplacian import quadtree_laplacian
 from .quadtree_boundary import quadtree_boundary
 from .quadtree_children import quadtree_children
+from .grad_intrinsic import grad_intrinsic
 from .grad import grad
 from .doublearea import doublearea
 from .doublearea_intrinsic import doublearea_intrinsic
@@ -107,5 +108,7 @@
 from .read_dmat import read_dmat
 from .linear_blend_skinning import linear_blend_skinning
 from .barycenters import barycenters
+from .biharmonic_energy import biharmonic_energy
+from .biharmonic_energy_intrinsic import biharmonic_energy_intrinsic
 from .adjacency_matrix import adjacency_matrix
 from .connected_components import connected_components

src/gpytoolbox/biharmonic_energy.py

@@ -0,0 +1,86 @@
+import numpy as np
+import scipy as sp
+from .halfedge_lengths_squared import halfedge_lengths_squared
+from .biharmonic_energy_intrinsic import biharmonic_energy_intrinsic
+from .doublearea import doublearea
+from .boundary_vertices import boundary_vertices
+from .massmatrix import massmatrix
+from .grad import grad
+
+
+def biharmonic_energy(V,F,
+    bc='mixedfem_zero_neumann'):
+    """Constructs the biharmonic energy matrix Q such that for a per-vertex function u, the discrete biharmonic energy is u'Qu.
+
+    Parameters
+    ----------
+    V : (n,d) numpy array
+        vertex list of a triangle mesh
+    F : (m,3) numpy int array
+        face index list of a triangle mesh
+    bc : string, optional (default 'mixedfem_zero_neumann')
+         Which type of discretization and boundary condition to use.
+         Options are: {'mixedfem_zero_neumann', 'hessian', 'curved_hessian'}
+         - 'mixedfem_zero_neumann' implements the mixed finite element
+         discretization from Jacobson et al. 2010.
+         "Mixed Finite Elements for Variational Surface Modeling".
+         - 'hessian' implements the Hessian energy from Stein et al. 2018.
+         "Natural Boundary Conditions for Smoothing in Geometry Processing".
+         - 'curved_hessian' implements the Hessian energy from Stein et al.
+         2020 "A Smoothness Energy without Boundary Distortion for Curved Surfaces"
+         via libigl C++ binding.
+
+    Returns
+    -------
+    Q : (n,n) scipy csr_matrix
+        biharmonic energy matrix
+
+    Examples
+    --------
+    ```python
+    # Mesh in V,F
+    import gpytoolbox as gpy
+    Q = gpy.biharmonic_energy(V,F)
+    ```
+    """
+
+    if bc=='hessian':
+        return _hessian_energy(V, F)
+    else:
+        l_sq = halfedge_lengths_squared(V,F)
+        return biharmonic_energy_intrinsic(l_sq,F,n=V.shape[0],bc=bc)
+
+
+def _hessian_energy(V, F):
+    assert F.shape[1]==3, "Only works on triangle meshes."
+
+    # Q = G' * A * D * Mtilde * D' * A * G
+    n = V.shape[0]
+    m = F.shape[0]
+    dim = V.shape[1]
+
+    b = boundary_vertices(F)
+    i = np.setdiff1d(np.arange(n),b)
+    ni = len(i)
+
+    a = doublearea(V, F) / 2.
+    A = sp.sparse.spdiags([np.tile(a, dim)], 0,
+        m=dim*m, n=dim*m, format='csr')
+
+    M_d = massmatrix(V, F, type='voronoi').diagonal()[i]
+    M_d_inv = 1. / M_d
+    M_tilde_inv = sp.sparse.spdiags([np.tile(M_d_inv, dim**2)], 0,
+        m=ni*dim**2, n=ni*dim**2, format='csr')
+
+    G = grad(V,F)
+    Gi = G[:,i]
+    Dlist = []
+    for d in range(dim):
+        Dlist.append(Gi[d*m:(d+1)*m, :])
+    Dblock = sp.sparse.bmat([Dlist], format='csr')
+    D = sp.sparse.block_diag(dim*[Dblock], format='csr')
+
+    Q = G.transpose() * A * D * M_tilde_inv * D.transpose() * A * G
+
+    return Q
+

src/gpytoolbox/biharmonic_energy_intrinsic.py

@@ -0,0 +1,92 @@
+import numpy as np
+import scipy as sp
+from .cotangent_laplacian_intrinsic import cotangent_laplacian_intrinsic
+from .massmatrix_intrinsic import massmatrix_intrinsic
+
+def biharmonic_energy_intrinsic(l_sq,F,
+    n=None,
+    bc='mixedfem_zero_neumann'):
+    """Constructs the biharmonic energy matrix Q such that for a per-vertex function u, the discrete biharmonic energy is u'Qu, using only intrinsic information from the mesh.
+
+    Parameters
+    ----------
+    l_sq : (m,3) numpy array
+        squared halfedge lengths as computed by halfedge_lengths_squared
+    F : (m,3) numpy int array
+        face index list of a triangle mesh
+    n : int, optional (default None)
+        number of vertices in the mesh.
+        If absent, will try to infer from F.
+    bc : string, optional (default 'mixedfem_zero_neumann')
+         Which type of discretization and boundary condition to use.
+         Options are: {'mixedfem_zero_neumann', 'hessian', 'curved_hessian'}
+         - 'mixedfem_zero_neumann' implements the mixed finite element
+         discretization from Jacobson et al. 2010.
+         "Mixed Finite Elements for Variational Surface Modeling".
+         - 'curved_hessian' implements the Hessian energy from Stein et al.
+         2020 "A Smoothness Energy without Boundary Distortion for Curved Surfaces"
+         via libigl C++ binding.
+
+    Returns
+    -------
+    Q : (n,n) scipy csr_matrix
+        biharmonic energy matrix
+
+    Examples
+    --------
+    ```python
+    # Mesh in V,F
+    import gpytoolbox as gpy
+    l_sq = gpy.halfedge_lengths_squared(V,F)
+    Q = biharmonic_energy_intrinsic(l_sq,F)
+    ```
+
+    """
+
+    assert F.shape[1] == 3, "Only works on triangle meshes."
+    assert l_sq.shape == F.shape
+    assert np.all(l_sq>=0)
+
+    if n==None:
+        n = np.max(F)+1
+
+    if bc=='mixedfem_zero_neumann':
+        Q = _mixedfem_neumann_laplacian_energy(l_sq, F, n)
+    elif bc=='curved_hessian':
+        Q = _curved_hessian_energy(l_sq, F, n)
+    else:
+        assert False, "Invalid bc"
+
+    return Q
+
+
+def _mixedfem_neumann_laplacian_energy(l_sq, F, n):
+    # Q = L' * M^{-1} * L
+    L = cotangent_laplacian_intrinsic(l_sq, F, n=n)
+    M = massmatrix_intrinsic(l_sq, F, n=n, type='voronoi')
+    M_inv = M.power(-1) #Sparse matrix inverts componentwise
+    Q = L.transpose() * M_inv * L
+
+    return Q
+
+
+try:
+    # Import C++ bindings for curved Hessian
+    from gpytoolbox_bindings import _curved_hessian_intrinsic_cpp_impl
+    _CPP_CURVED_HESSIAN_INTRINSIC_AVAILABLE = True
+except Exception as e:
+    _CPP_CURVED_HESSIAN_INTRINSIC_AVAILABLE = False
+
+def _curved_hessian_energy(l_sq, F, n):
+    assert _CPP_CURVED_HESSIAN_INTRINSIC_AVAILABLE, \
+    "C++ bindings for curved Hessian not available."
+
+    Q = _curved_hessian_intrinsic_cpp_impl(l_sq.astype(np.float64),
+        F.astype(np.int32))
+    if Q.shape[0] != n:
+        Q = sp.sparse.block_diag([
+            Q, sp.sparse.csr_matrix((n-Q.shape[0], n-Q.shape[0]))],
+            format='csr')
+
+    return Q
+

src/gpytoolbox/grad.py

@@ -30,6 +30,7 @@ def grad(V,F=None):
 
     Notes
     -----
+    Adapted from https://github.com/alecjacobson/gptoolbox/blob/master/mesh/grad.m
 
     Examples
     --------

src/gpytoolbox/grad_intrinsic.py

@@ -0,0 +1,68 @@
+import numpy as np
+import scipy as sp
+
+from .grad import grad
+
+# Lifted from https://github.com/alecjacobson/gptoolbox/blob/master/mesh/grad.m
+
+def grad_intrinsic(l_sq,F,
+    n=None,):
+    """Intrinsic finite element gradient matrix
+
+    Given a triangle mesh, computes the finite element gradient matrix assuming
+    piecewise linear hat function basis.
+
+    Parameters
+    ----------
+    l_sq : (m,3) numpy array
+        squared halfedge lengths as computed by halfedge_lengths_squared
+    F : (m,3) numpy int array
+        face index list of a triangle mesh
+    n : int, optional (default None)
+        number of vertices in the mesh.
+        If absent, will try to infer from F.
+
+    Returns
+    -------
+    G : (2*m,n) scipy sparse.csr_matrix
+        Sparse FEM gradient matrix.
+        The first m rows ar ethe gradient with respect to the edge (1,2).
+        The m rows after that run in a pi/2 rotation counter-clockwise.
+
+    See Also
+    --------
+    cotangent_laplacian.
+
+    Notes
+    -----
+    Adapted from https://github.com/alecjacobson/gptoolbox/blob/master/mesh/grad_intrinsic.m
+
+    Examples
+    --------
+    ```python
+    import gpytoolbox as gpy
+    l = gpy.halfedge_lengths_squared(V,F)
+    G = gpy.grad_intrinsic(l_sq,F)
+    ```
+    """
+
+
+    assert F.shape[1] == 3, "Only works on triangle meshes."
+    assert l_sq.shape == F.shape
+    assert np.all(l_sq>=0)
+
+    if n==None:
+        n = np.max(F)+1
+    m = F.shape[0]
+
+    l0 = np.sqrt(l_sq[:,0])[:,None]
+    gx = (l_sq[:,1][:,None]-l_sq[:,0][:,None]-l_sq[:,2][:,None]) / (-2.*l0)
+    gy = np.sqrt(l_sq[:,2][:,None] - gx**2)
+    V2 = np.block([[gx,gy], [np.zeros((m,2))], [l0, np.zeros((m,1))]])
+    F2 = np.reshape(np.arange(3*m), (m,3), order='F')
+    G2 = grad(V2,F2)
+    P = sp.sparse.csr_matrix((np.ones(3*m),(F2.ravel(),F.ravel())),
+        shape=(3*m,n))
+    G = G2*P
+
+    return G