phynexis-deps/igl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added vector Dirichlet and curved Hessian functionality (#1560)

Browse Source 
* Added vector Dirichlet and curved Hessian functionality, as well as tutorials.

* Applied jdumas changes

* added github actions for this branch

* more jdumas changes

* Changed int to appropriate scalar

* moved around Int definition

* forgot braces somewhere

* forgot braces

* Applied changes suggested by jdumas

* forgot braces

* remove custom github workflow
This commit is contained in:
Oded Stein
2020-09-26 20:06:26 -07:00
committed by GitHub
parent 28008528b4
commit 8ea0b8c6dc
32 changed files with 2377 additions and 75 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/igl/average_from_edges_onto_vertices.cpp
+68
View File
@@ -0,0 +1,68 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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 "average_from_edges_onto_vertices.h"
template<typename DerivedF,typename DerivedE,typename DerivedoE,
typename DeriveduE,typename DeriveduV>
IGL_INLINE void
igl::average_from_edges_onto_vertices(
const Eigen::MatrixBase<DerivedF> &F,
const Eigen::MatrixBase<DerivedE> &E,
const Eigen::MatrixBase<DerivedoE> &oE,
const Eigen::MatrixBase<DeriveduE> &uE,
Eigen::PlainObjectBase<DeriveduV> &uV)
{
using Scalar = typename DeriveduE::Scalar;
using VecX = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
using Int = typename DerivedF::Scalar;
assert(E.rows()==F.rows() && "E does not match dimensions of F.");
assert(oE.rows()==F.rows() && "oE does not match dimensions of F.");
assert(E.cols()==3 && F.cols()==3 && oE.cols()==3 &&
"This method is for triangle meshes.");
const Int n = F.maxCoeff()+1;
VecX edgesPerVertex(n);
edgesPerVertex.setZero();
uV.resize(n,1);
uV.setZero();
for(Eigen::Index i=0; i<F.rows(); ++i) {
for(int j=0; j<3; ++j) {
if(oE(i,j)<0) {
continue;
}
const Int e = E(i,j);
const Int vi=F(i,(j+1)%3), vj=F(i,(j+2)%3);
//Count vertex valence
++edgesPerVertex(vi);
++edgesPerVertex(vj);
//Average uE value onto vertices
uV(vi) += uE(e);
uV(vj) += uE(e);
}
}
//Divide by valence
for(Int i=0; i<n; ++i) {
const Scalar valence = edgesPerVertex(i);
if(valence>0) {
uV(i) /= valence;
}
}
}
#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
template void igl::average_from_edges_onto_vertices<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::PartialReduxExpr<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::internal::member_norm<double>, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::PartialReduxExpr<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::internal::member_norm<double>, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
template void igl::average_from_edges_onto_vertices<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<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
template void igl::average_from_edges_onto_vertices<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<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
#endif
include/igl/average_from_edges_onto_vertices.h
+39
View File
@@ -0,0 +1,39 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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_AVERAGE_FROM_EDGES_ONTO_VERTICES_H
#define IGL_AVERAGE_FROM_EDGES_ONTO_VERTICES_H
#include "igl_inline.h"
#include <Eigen/Dense>
namespace igl
{
// Move a scalar field defined on edges to vertices by averaging
//
// Input:
// F: triangle mesh connectivity
// E, oE: mapping from halfedges to edges and orientation as generated by
// orient_halfedges
// uE: scalar field defined on edges, one per edge
//
// Output:
// uV: scalar field defined on vertices
template<typename DerivedF,typename DerivedE,typename DerivedoE,
typename DeriveduE,typename DeriveduV>
IGL_INLINE void average_from_edges_onto_vertices(
const Eigen::MatrixBase<DerivedF> &F,
const Eigen::MatrixBase<DerivedE> &E,
const Eigen::MatrixBase<DerivedoE> &oE,
const Eigen::MatrixBase<DeriveduE> &uE,
Eigen::PlainObjectBase<DeriveduV> &uV);
}
#ifndef IGL_STATIC_LIBRARY
# include "average_from_edges_onto_vertices.cpp"
#endif
#endif
include/igl/cr_vector_curvature_correction.cpp
+189
View File
@@ -0,0 +1,189 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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 "cr_vector_curvature_correction.h"
#include "orient_halfedges.h"
#include "gaussian_curvature.h"
#include "squared_edge_lengths.h"
#include "doublearea.h"
#include "boundary_loop.h"
#include "internal_angles.h"
#include "PI.h"
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarK>
IGL_INLINE void
igl::cr_vector_curvature_correction(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarK>& K)
{
Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic>
l_sq;
squared_edge_lengths(V, F, l_sq);
cr_vector_curvature_correction_intrinsic(F, l_sq, E, oE, K);
}
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarK>
IGL_INLINE void
igl::cr_vector_curvature_correction(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
Eigen::PlainObjectBase<DerivedE>& E,
Eigen::PlainObjectBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarK>& K)
{
if(E.rows()!=F.rows() || E.cols()!=F.cols() || oE.rows()!=F.rows() ||
oE.cols()!=F.cols()) {
orient_halfedges(F, E, oE);
}
const Eigen::PlainObjectBase<DerivedE>& cE = E;
const Eigen::PlainObjectBase<DerivedOE>& coE = oE;
cr_vector_curvature_correction(V, F, cE, coE, K);
}
template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
typename DerivedOE, typename ScalarK>
IGL_INLINE void
igl::cr_vector_curvature_correction_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarK>& K)
{
Eigen::Matrix<typename DerivedL_sq::Scalar,Eigen::Dynamic,Eigen::Dynamic>
theta;
internal_angles_using_squared_edge_lengths(l_sq, theta);
cr_vector_curvature_correction_intrinsic(F, l_sq, theta, E, oE, K);
}
template <typename DerivedF, typename DerivedL_sq, typename Derivedtheta,
typename DerivedE, typename DerivedOE,
typename ScalarK>
IGL_INLINE void
igl::cr_vector_curvature_correction_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<Derivedtheta>& theta,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarK>& K)
{
// Compute the angle defect kappa, set it to 0 at the boundary
const typename DerivedF::Scalar n = F.maxCoeff() + 1;
Eigen::Matrix<typename DerivedL_sq::Scalar,Eigen::Dynamic,1> kappa(n);
kappa.setZero();
for(Eigen::Index i=0; i<F.rows(); ++i) {
for(int j=0; j<3; ++j) {
kappa(F(i,j)) -= theta(i,j);
}
}
kappa.array() += 2 * PI;
std::vector<std::vector<typename DerivedF::Scalar> > b;
boundary_loop(F, b);
for(const auto& loop : b) {
for(auto v : loop) {
kappa(v) = 0;
}
}
cr_vector_curvature_correction_intrinsic(F, l_sq, theta, kappa, E, oE, K);
}
template <typename DerivedF, typename DerivedL_sq, typename Derivedtheta,
typename Derivedkappa, typename DerivedE, typename DerivedOE,
typename ScalarK>
IGL_INLINE void
igl::cr_vector_curvature_correction_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<Derivedtheta>& theta,
const Eigen::MatrixBase<Derivedkappa>& kappa,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarK>& K)
{
assert(F.cols()==3 && "Faces have three vertices");
assert(E.rows()==F.rows() && E.cols()==F.cols() && oE.rows()==F.rows() &&
theta.rows()==F.rows() && theta.cols()==F.cols() &&
oE.cols()==F.cols() && "Wrong dimension in edge vectors");
assert(kappa.rows()==F.maxCoeff()+1 &&
"Wrong dimension in theta or kappa");
const Eigen::Index m = F.rows();
const typename DerivedE::Scalar nE = E.maxCoeff() + 1;
//Divide kappa by the actual angle sum to weigh consistently.
Derivedtheta angleSum = Derivedtheta::Zero(kappa.rows(), 1);
for(Eigen::Index i=0; i<F.rows(); ++i) {
for(int j=0; j<3; ++j) {
angleSum(F(i,j)) += theta(i,j);
}
}
const Eigen::Matrix<typename Derivedkappa::Scalar, Eigen::Dynamic, 1>
scaledKappa = kappa.array() / angleSum.array();
std::vector<Eigen::Triplet<ScalarK> > tripletList;
tripletList.reserve(10*3*m);
for(Eigen::Index f=0; f<m; ++f) {
for(int e=0; e<3; ++e) {
const ScalarK eij=l_sq(f,e), ejk=l_sq(f,(e+1)%3),
eki=l_sq(f,(e+2)%3); //These are squared quantities.
const ScalarK lens = sqrt(eij*eki);
const ScalarK o = oE(f,e)*oE(f,(e+2)%3);
const typename DerivedF::Scalar i=F(f,(e+1)%3), j=F(f,(e+2)%3), k=F(f,e);
const ScalarK ki=scaledKappa(i)*theta(f,(e+1)%3),
kj=scaledKappa(j)*theta(f,(e+2)%3), kk=scaledKappa(k)*theta(f,e);
const ScalarK costhetaidiv = (eij-ejk+eki)/(2.*lens);
const ScalarK sinthetaidiv = sqrt( (1.-pow(eij-ejk+eki,2)/
(4.*eij*eki)) );
const ScalarK Corrijij = (ki+kj+kk);
tripletList.emplace_back(E(f,e), E(f,e), Corrijij);
tripletList.emplace_back(E(f,e)+nE, E(f,e)+nE, Corrijij);
const ScalarK Corrijki = -o*(ki-kj-kk)*costhetaidiv;
tripletList.emplace_back(E(f,e), E(f,(e+2)%3), Corrijki);
tripletList.emplace_back(E(f,(e+2)%3), E(f,e), Corrijki);
tripletList.emplace_back(E(f,e)+nE, E(f,(e+2)%3)+nE, Corrijki);
tripletList.emplace_back(E(f,(e+2)%3)+nE, E(f,e)+nE, Corrijki);
const ScalarK Corrijkiperp = o*(ki-kj-kk)*sinthetaidiv;
tripletList.emplace_back(E(f,e), E(f,(e+2)%3)+nE, Corrijkiperp);
tripletList.emplace_back(E(f,(e+2)%3)+nE, E(f,e), Corrijkiperp);
tripletList.emplace_back(E(f,e)+nE, E(f,(e+2)%3), -Corrijkiperp);
tripletList.emplace_back(E(f,(e+2)%3), E(f,e)+nE, -Corrijkiperp);
}
}
K.resize(2*nE, 2*nE);
K.setFromTriplets(tripletList.begin(), tripletList.end());
}
#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
template void igl::cr_vector_curvature_correction_intrinsic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
template void igl::cr_vector_curvature_correction<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::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::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::SparseMatrix<double, 0, int>&);
#endif
include/igl/cr_vector_curvature_correction.h
+115
View File
@@ -0,0 +1,115 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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_CR_VECTOR_CURVATURE_CORRECTION_H
#define IGL_CR_VECTOR_CURVATURE_CORRECTION_H
#include "igl_inline.h"
#include <Eigen/Core>
#include <Eigen/Sparse>
namespace igl
{
// Computes the vector Crouzeix-Raviart curvature correction
// term of Oded Stein, Alec Jacobson, Max Wardetzky, Eitan
// Grinspun, 2020. "A Smoothness Energy without Boundary Distortion for
// Curved Surfaces", but using the basis functions by Oded Stein,
// Max Wardetzky, Alec Jacobson, Eitan Grinspun, 2020.
// "A Simple Discretization of the Vector Dirichlet Energy"
//
// Inputs:
// V, F: input mesh
// E: a mapping from each halfedge to each edge, as computed with
// orient_halfedges.
// will be computed if not provided.
// oE: the orientation of each halfedge compared to the orientation of the
// actual edge, as computed with orient_halfedges.
// will be computed if not provided.
//
// Outputs:
// K: computed curvature correction matrix
// E, oE: these are computed if they are not present, as described above
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarK>
IGL_INLINE void
cr_vector_curvature_correction(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarK>& K);
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarK>
IGL_INLINE void
cr_vector_curvature_correction(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
Eigen::PlainObjectBase<DerivedE>& E,
Eigen::PlainObjectBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarK>& K);
// Version that uses intrinsic quantities as input
//
// Inputs:
// F: input mesh connectivity
// l_sq: squared edge lengths of each halfedge
// theta: the tip angles at each halfedge
// kappa: the Gaussian curvature at each vertex
// E: a mapping from each halfedge to each edge.
// oE: the orientation of each halfedge compared to the orientation of the
// actual edge.
//
// Outputs:
// K: computed curvature correction matrix
template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
typename DerivedOE, typename ScalarK>
IGL_INLINE void
cr_vector_curvature_correction_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarK>& K);
template <typename DerivedF, typename DerivedL_sq, typename Derivedtheta,
typename DerivedE, typename DerivedOE,
typename ScalarK>
IGL_INLINE void
cr_vector_curvature_correction_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<Derivedtheta>& theta,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarK>& K);
template <typename DerivedF, typename DerivedL_sq, typename Derivedtheta,
typename Derivedkappa, typename DerivedE, typename DerivedOE,
typename ScalarK>
IGL_INLINE void
cr_vector_curvature_correction_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<Derivedtheta>& theta,
const Eigen::MatrixBase<Derivedkappa>& kappa,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarK>& K);
}
#ifndef IGL_STATIC_LIBRARY
# include "cr_vector_curvature_correction.cpp"
#endif
#endif
include/igl/cr_vector_laplacian.cpp
+128
View File
@@ -0,0 +1,128 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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 "cr_vector_laplacian.h"
#include <vector>
#include "orient_halfedges.h"
#include "doublearea.h"
#include "squared_edge_lengths.h"
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarL>
IGL_INLINE void
igl::cr_vector_laplacian(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarL>& L)
{
Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic>
l_sq;
squared_edge_lengths(V, F, l_sq);
cr_vector_laplacian_intrinsic(F, l_sq, E, oE, L);
}
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarL>
IGL_INLINE void
igl::cr_vector_laplacian(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
Eigen::PlainObjectBase<DerivedE>& E,
Eigen::PlainObjectBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarL>& L)
{
if(E.rows()!=F.rows() || E.cols()!=F.cols() || oE.rows()!=F.rows() ||
oE.cols()!=F.cols()) {
orient_halfedges(F, E, oE);
}
const Eigen::PlainObjectBase<DerivedE>& cE = E;
const Eigen::PlainObjectBase<DerivedOE>& coE = oE;
cr_vector_laplacian(V, F, cE, coE, L);
}
template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
typename DerivedOE, typename ScalarL>
IGL_INLINE void
igl::cr_vector_laplacian_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarL>& L)
{
Eigen::Matrix<typename DerivedL_sq::Scalar, Eigen::Dynamic, Eigen::Dynamic>
dA;
DerivedL_sq l_sqrt = l_sq.array().sqrt().matrix();
doublearea(l_sqrt, dA);
cr_vector_laplacian_intrinsic(F, l_sq, dA, E, oE, L);
}
template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
typename DerivedE, typename DerivedOE, typename ScalarL>
IGL_INLINE void
igl::cr_vector_laplacian_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DeriveddA>& dA,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarL>& L)
{
assert(F.cols()==3 && "Faces have three vertices");
assert(E.rows()==F.rows() && E.cols()==F.cols() && oE.rows()==F.rows() &&
oE.cols()==F.cols() && "Wrong dimension in edge vectors");
assert(l_sq.rows()==F.rows() && l_sq.cols()==3 && "l_sq dimensions wrong");
assert(dA.size()==F.rows() && "dA dimensions wrong");
const Eigen::Index m = F.rows();
const typename DerivedE::Scalar nE = E.maxCoeff() + 1;
std::vector<Eigen::Triplet<ScalarL> > tripletList;
tripletList.reserve(10*3*m);
for(Eigen::Index f=0; f<m; ++f) {
for(int e=0; e<3; ++e) {
const ScalarL eij=l_sq(f,e), ejk=l_sq(f,(e+1)%3),
eki=l_sq(f,(e+2)%3); //These are squared quantities.
const ScalarL lens = sqrt(eij*eki);
const ScalarL o = oE(f,e)*oE(f,(e+2)%3);
tripletList.emplace_back(E(f,e), E(f,e), 2./dA(f) * eij);
tripletList.emplace_back(E(f,e)+nE, E(f,e)+nE, 2./dA(f) * eij);
const ScalarL Dijki = o * pow(eij-ejk+eki,2)/(2.*lens*dA(f));
tripletList.emplace_back(E(f,e), E(f,(e+2)%3), Dijki);
tripletList.emplace_back(E(f,(e+2)%3), E(f,e), Dijki);
tripletList.emplace_back(E(f,e)+nE, E(f,(e+2)%3)+nE, Dijki);
tripletList.emplace_back(E(f,(e+2)%3)+nE, E(f,e)+nE, Dijki);
const ScalarL Dijkiperp = -o * (eij-ejk+eki)/lens;
tripletList.emplace_back(E(f,e), E(f,(e+2)%3)+nE, Dijkiperp);
tripletList.emplace_back(E(f,(e+2)%3)+nE, E(f,e), Dijkiperp);
tripletList.emplace_back(E(f,e)+nE, E(f,(e+2)%3), -Dijkiperp);
tripletList.emplace_back(E(f,(e+2)%3), E(f,e)+nE, -Dijkiperp);
}
}
L.resize(2*nE, 2*nE);
L.setFromTriplets(tripletList.begin(), tripletList.end());
}
#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
template void igl::cr_vector_laplacian<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::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::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::SparseMatrix<double, 0, int>&);
template void igl::cr_vector_laplacian_intrinsic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -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/cr_vector_laplacian.h
+99
View File
@@ -0,0 +1,99 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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_CR_VECTOR_LAPLACIAN_H
#define IGL_CR_VECTOR_LAPLACIAN_H
#include "igl_inline.h"
#include <Eigen/Core>
#include <Eigen/Sparse>
namespace igl
{
// Computes the CR vector Laplacian matrix.
// See Oded Stein, Max Wardetzky, Alec Jacobson, Eitan Grinspun, 2020.
// "A Simple Discretization of the Vector Dirichlet Energy"
//
// Inputs:
// V, F: input mesh
// E: a mapping from each halfedge to each edge, as computed with
// orient_halfedges.
// will be computed if not provided.
// oE: the orientation of each halfedge compared to the orientation of the
// actual edge, as computed with orient_halfedges.
// will be computed if not provided.
//
// Outputs:
// L: computed Laplacian matrix
// E, oE: these are computed if they are not present, as described above
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarL>
IGL_INLINE void
cr_vector_laplacian(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarL>& L);
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarL>
IGL_INLINE void
cr_vector_laplacian(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
Eigen::PlainObjectBase<DerivedE>& E,
Eigen::PlainObjectBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarL>& L);
// Version that uses intrinsic quantities as input
//
// Inputs:
// F: input mesh connectivity
// l_sq: squared edge lengths of each halfedge
// dA: double area of each face
// E: a mapping from each halfedge to each edge.
// oE: the orientation of each halfedge compared to the orientation of the
// actual edge.
//
// Outputs:
// L: computed Laplacian matrix
template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
typename DerivedOE, typename ScalarL>
IGL_INLINE void
cr_vector_laplacian_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarL>& L);
template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
typename DerivedE, typename DerivedOE, typename ScalarL>
IGL_INLINE void
cr_vector_laplacian_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DeriveddA>& dA,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarL>& L);
}
#ifndef IGL_STATIC_LIBRARY
# include "cr_vector_laplacian.cpp"
#endif
#endif
include/igl/cr_vector_mass.cpp
+112
View File
@@ -0,0 +1,112 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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 "cr_vector_mass.h"
#include <vector>
#include "orient_halfedges.h"
#include "doublearea.h"
#include "squared_edge_lengths.h"
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarM>
IGL_INLINE void
igl::cr_vector_mass(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarM>& M)
{
Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic>
l_sq;
squared_edge_lengths(V, F, l_sq);
cr_vector_mass_intrinsic(F, l_sq, E, oE, M);
}
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarM>
IGL_INLINE void
igl::cr_vector_mass(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
Eigen::PlainObjectBase<DerivedE>& E,
Eigen::PlainObjectBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarM>& M)
{
if(E.rows()!=F.rows() || E.cols()!=F.cols() || oE.rows()!=F.rows() ||
oE.cols()!=F.cols()) {
orient_halfedges(F, E, oE);
}
const Eigen::PlainObjectBase<DerivedE>& cE = E;
const Eigen::PlainObjectBase<DerivedOE>& coE = oE;
cr_vector_mass(V, F, cE, coE, M);
}
template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
typename DerivedOE, typename ScalarM>
IGL_INLINE void
igl::cr_vector_mass_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarM>& M)
{
Eigen::Matrix<typename DerivedL_sq::Scalar, Eigen::Dynamic, Eigen::Dynamic>
dA;
DerivedL_sq l_sqrt = l_sq.array().sqrt().matrix();
doublearea(l_sqrt, dA);
cr_vector_mass_intrinsic(F, l_sq, dA, E, oE, M);
}
template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
typename DerivedE, typename DerivedOE, typename ScalarM>
IGL_INLINE void
igl::cr_vector_mass_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DeriveddA>& dA,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarM>& M)
{
assert(F.cols()==3 && "Faces have three vertices");
assert(E.rows()==F.rows() && E.cols()==F.cols() && oE.rows()==F.rows() &&
oE.cols()==F.cols() && "Wrong dimension in edge vectors");
const Eigen::Index m = F.rows();
const typename DerivedE::Scalar nE = E.maxCoeff() + 1;
std::vector<Eigen::Triplet<ScalarM> > tripletList;
tripletList.reserve(2*3*m);
for(Eigen::Index f=0; f<m; ++f) {
for(int e=0; e<3; ++e) {
const typename DerivedF::Scalar v1=F(f,(e+1)%3), v2=F(f,(e+2)%3);
//Scaled
const ScalarM entry = dA(f) / 6;
tripletList.emplace_back(E(f,e), E(f,e), entry);
tripletList.emplace_back(E(f,e)+nE, E(f,e)+nE, entry);
}
}
M.resize(2*nE, 2*nE);
M.setFromTriplets(tripletList.begin(), tripletList.end());
}
#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
template void igl::cr_vector_mass<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::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::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::SparseMatrix<double, 0, int>&);
template void igl::cr_vector_mass_intrinsic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -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/cr_vector_mass.h
+100
View File
@@ -0,0 +1,100 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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_CR_VECTOR_MASS
#define IGL_CR_VECTOR_MASS
#include "igl_inline.h"
#include <Eigen/Core>
#include <Eigen/Sparse>
namespace igl
{
// Computes the CR vector mass matrix, using an arrangement of all parallel
// degrees of freedom first, and all perpendicular degrees of freedom next.
// See Oded Stein, Max Wardetzky, Alec Jacobson, Eitan Grinspun, 2020.
// "A Simple Discretization of the Vector Dirichlet Energy"
//
// Inputs:
// V, F: input mesh
// E: a mapping from each halfedge to each edge, as computed with
// orient_halfedges.
// will be computed if not provided.
// oE: the orientation of each halfedge compared to the orientation of the
// actual edge, as computed with orient_halfedges.
// will be computed if not provided.
//
// Outputs:
// M: computed mass matrix
// E, oE: these are computed if they are not present, as described above
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarM>
IGL_INLINE void
cr_vector_mass(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarM>& M);
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarM>
IGL_INLINE void
cr_vector_mass(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
Eigen::PlainObjectBase<DerivedE>& E,
Eigen::PlainObjectBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarM>& M);
// Version that uses intrinsic quantities as input
//
// Inputs:
// F: input mesh connectivity
// l_sq: squared edge lengths of each halfedge
// dA: double area of each face
// E: a mapping from each halfedge to each edge.
// oE: the orientation of each halfedge compared to the orientation of the
// actual edge.
//
// Outputs:
// M: computed mass matrix
template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
typename DerivedOE, typename ScalarM>
IGL_INLINE void
cr_vector_mass_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarM>& M);
template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
typename DerivedE, typename DerivedOE, typename ScalarM>
IGL_INLINE void
cr_vector_mass_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DeriveddA>& dA,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarM>& M);
}
#ifndef IGL_STATIC_LIBRARY
# include "cr_vector_mass.cpp"
#endif
#endif
include/igl/crouzeix_raviart_cotmatrix.cpp
+1
View File
@@ -97,4 +97,5 @@ void igl::crouzeix_raviart_cotmatrix(
// Explicit template instantiation
// generated by autoexplicit.sh
template void igl::crouzeix_raviart_cotmatrix<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::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::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int>&);
template void igl::crouzeix_raviart_cotmatrix<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::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::MatrixBase<Eigen::Matrix<int, -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/crouzeix_raviart_massmatrix.cpp
+1
View File
@@ -82,4 +82,5 @@ void igl::crouzeix_raviart_massmatrix(
template void igl::crouzeix_raviart_massmatrix<double, 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::Matrix<int, -1, 1, 0, -1, 1> >(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>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
template void igl::crouzeix_raviart_massmatrix<double, 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::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int>&);
template void igl::crouzeix_raviart_massmatrix<float, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<float, 0, int>&);
template void igl::crouzeix_raviart_massmatrix<double, 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::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -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/curved_hessian_energy.cpp
+127
View File
@@ -0,0 +1,127 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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 "curved_hessian_energy.h"
#include "orient_halfedges.h"
#include "doublearea.h"
#include "squared_edge_lengths.h"
#include "cr_vector_laplacian.h"
#include "cr_vector_mass.h"
#include "cr_vector_curvature_correction.h"
#include "scalar_to_cr_vector_gradient.h"
template <typename DerivedV, typename DerivedF, typename ScalarQ>
IGL_INLINE void
igl::curved_hessian_energy(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
Eigen::SparseMatrix<ScalarQ>& Q)
{
Eigen::MatrixXi E, oE;
curved_hessian_energy(V, F, E, oE, Q);
}
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarQ>
IGL_INLINE void
igl::curved_hessian_energy(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarQ>& Q)
{
Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic>
l_sq;
squared_edge_lengths(V, F, l_sq);
curved_hessian_energy_intrinsic(F, l_sq, E, oE, Q);
}
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarQ>
IGL_INLINE void
igl::curved_hessian_energy(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
Eigen::PlainObjectBase<DerivedE>& E,
Eigen::PlainObjectBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarQ>& Q)
{
if(E.rows()!=F.rows() || E.cols()!=F.cols() || oE.rows()!=F.rows() ||
oE.cols()!=F.cols()) {
orient_halfedges(F, E, oE);
}
const Eigen::PlainObjectBase<DerivedE>& cE = E;
const Eigen::PlainObjectBase<DerivedOE>& coE = oE;
curved_hessian_energy(V, F, cE, coE, Q);
}
template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
typename DerivedOE, typename ScalarQ>
IGL_INLINE void
igl::curved_hessian_energy_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarQ>& Q)
{
Eigen::Matrix<typename DerivedL_sq::Scalar, Eigen::Dynamic, Eigen::Dynamic>
dA;
Eigen::Matrix<typename DerivedL_sq::Scalar, Eigen::Dynamic, Eigen::Dynamic>
l_sqrt = l_sq.array().sqrt().matrix();
doublearea(l_sqrt, dA);
curved_hessian_energy_intrinsic(F, l_sq, dA, E, oE, Q);
}
template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
typename DerivedE, typename DerivedOE, typename ScalarQ>
IGL_INLINE void
igl::curved_hessian_energy_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DeriveddA>& dA,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarQ>& Q)
{
//Matrices that need to be combined
Eigen::SparseMatrix<ScalarQ> M, D, L, K;
cr_vector_mass_intrinsic(F, l_sq, dA, E, oE, M);
scalar_to_cr_vector_gradient_intrinsic(F, l_sq, dA, E, oE, D);
cr_vector_laplacian_intrinsic(F, l_sq, dA, E, oE, L);
cr_vector_curvature_correction_intrinsic(F, l_sq, E, oE, K);
//Invert M
std::vector<Eigen::Triplet<ScalarQ> > tripletListMi;
for(Eigen::Index k=0; k<M.outerSize(); ++k) {
for(typename Eigen::SparseMatrix<ScalarQ>::InnerIterator it(M,k);
it; ++it) {
if(it.value() > 0) {
tripletListMi.emplace_back(it.row(), it.col(), 1./it.value());
}
}
}
Eigen::SparseMatrix<ScalarQ> Mi(M.rows(), M.cols());
Mi.setFromTriplets(tripletListMi.begin(), tripletListMi.end());
//Hessian energy matrix
Q = D.transpose()*Mi*(L + K)*Mi*D;
}
#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
template void igl::curved_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/curved_hessian_energy.h
+114
View File
@@ -0,0 +1,114 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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_CURVED_HESSIAN_ENERGY_H
#define IGL_CURVED_HESSIAN_ENERGY_H
#include "igl_inline.h"
#include <Eigen/Core>
#include <Eigen/Sparse>
namespace igl
{
// Computes the curved Hessian energy using the Crouzeix-Raviart
// discretization.
// See Oded Stein, Alec Jacobson, Max Wardetzky, Eitan Grinspun, 2020.
// "A Smoothness Energy without Boundary Distortion for Curved Surfaces"
//
// Inputs:
// V, F: input mesh
//
// Outputs:
// Q: computed Hessian energy matrix
template <typename DerivedV, typename DerivedF, typename ScalarQ>
IGL_INLINE void
curved_hessian_energy(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
Eigen::SparseMatrix<ScalarQ>& Q);
// Version that exposes the edge orientation used.
//
// Inputs:
// V, F: input mesh
// E: a mapping from each halfedge to each edge, as computed with
// orient_halfedges.
// will be computed if not provided.
// oE: the orientation of each halfedge compared to the orientation of the
// actual edge, as computed with orient_halfedges.
// will be computed if not provided.
//
// Outputs:
// Q: computed Hessian energy matrix
// E, oE: these are computed if they are not present, as described above
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarQ>
IGL_INLINE void
curved_hessian_energy(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarQ>& Q);
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarQ>
IGL_INLINE void
curved_hessian_energy(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
Eigen::PlainObjectBase<DerivedE>& E,
Eigen::PlainObjectBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarQ>& Q);
// Version that uses intrinsic quantities as input
//
// Inputs:
// F: input mesh connectivity
// l_sq: squared edge lengths of each halfedge
// dA: double area of each face
// E: a mapping from each halfedge to each edge.
// oE: the orientation of each halfedge compared to the orientation of the
// actual edge.
//
// Outputs:
// Q: computed Hessian energy matrix
template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
typename DerivedOE, typename ScalarQ>
IGL_INLINE void
curved_hessian_energy_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarQ>& Q);
template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
typename DerivedE, typename DerivedOE, typename ScalarQ>
IGL_INLINE void
curved_hessian_energy_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DeriveddA>& dA,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarQ>& Q);
}
#ifndef IGL_STATIC_LIBRARY
# include "curved_hessian_energy.cpp"
#endif
#endif
include/igl/edge_midpoints.cpp
+48
View File
@@ -0,0 +1,48 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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 "edge_midpoints.h"
template<typename DerivedV,typename DerivedF,typename DerivedE,
typename DerivedoE, typename Derivedmps>
IGL_INLINE void
igl::edge_midpoints(
const Eigen::MatrixBase<DerivedV> &V,
const Eigen::MatrixBase<DerivedF> &F,
const Eigen::MatrixBase<DerivedE> &E,
const Eigen::MatrixBase<DerivedoE> &oE,
Eigen::PlainObjectBase<Derivedmps> &mps)
{
assert(E.rows()==F.rows() && "E does not match dimensions of F.");
assert(oE.rows()==F.rows() && "oE does not match dimensions of F.");
assert(E.cols()==3 && F.cols()==3 && oE.cols()==3 &&
"This method is for triangle meshes.");
assert(F.maxCoeff()<V.rows() && "V does not seem to belong to F.");
using ScalarE = typename DerivedE::Scalar;
using ScalarF = typename DerivedF::Scalar;
const ScalarE m = E.maxCoeff()+1;
mps.resize(m, V.cols());
for(Eigen::Index i=0; i<F.rows(); ++i) {
for(int j=0; j<3; ++j) {
if(oE(i,j)<0) {
continue;
}
const ScalarE e = E(i,j);
const ScalarF vi=F(i,(j+1)%3), vj=F(i,(j+2)%3);
mps.row(e) = 0.5*(V.row(vi) + V.row(vj));
}
}
}
#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
template void igl::edge_midpoints<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::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
#endif
include/igl/edge_midpoints.h
+38
View File
@@ -0,0 +1,38 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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_EDGE_MIDPOINTS_H
#define IGL_EDGE_MIDPOINTS_H
#include "igl_inline.h"
#include <Eigen/Dense>
namespace igl
{
// Computes the midpoints of edges in a triangle mesh.
//
// Input:
// V, F: triangle mesh
// E, oE: mapping from halfedges to edges and orientation as generated by
// orient_halfedges
//
// Output:
// mps: edge midpoints, one per edge in E
template<typename DerivedV,typename DerivedF,typename DerivedE,
typename DerivedoE, typename Derivedmps>
IGL_INLINE void edge_midpoints(
const Eigen::MatrixBase<DerivedV> &V,
const Eigen::MatrixBase<DerivedF> &F,
const Eigen::MatrixBase<DerivedE> &E,
const Eigen::MatrixBase<DerivedoE> &oE,
Eigen::PlainObjectBase<Derivedmps> &mps);
}
#ifndef IGL_STATIC_LIBRARY
# include "edge_midpoints.cpp"
#endif
#endif
include/igl/edge_vectors.cpp
+94
View File
@@ -0,0 +1,94 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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 "edge_vectors.h"
#include <Eigen/Geometry>
#include "per_face_normals.h"
#include "PI.h"
template<typename DerivedV,typename DerivedF,typename DerivedE,
typename DerivedoE, typename Derivedvec>
IGL_INLINE void
igl::edge_vectors(
const Eigen::MatrixBase<DerivedV> &V,
const Eigen::MatrixBase<DerivedF> &F,
const Eigen::MatrixBase<DerivedE> &E,
const Eigen::MatrixBase<DerivedoE> &oE,
Eigen::PlainObjectBase<Derivedvec> &vec)
{
Eigen::Matrix<typename Derivedvec::Scalar, Eigen::Dynamic, Eigen::Dynamic>
dummy;
edge_vectors<false>(V, F, E, oE, vec, dummy);
}
template<bool computePerpendicular,
typename DerivedV,typename DerivedF,typename DerivedE,
typename DerivedoE, typename DerivedvecParallel,
typename DerivedvecPerpendicular>
IGL_INLINE void
igl::edge_vectors(
const Eigen::MatrixBase<DerivedV> &V,
const Eigen::MatrixBase<DerivedF> &F,
const Eigen::MatrixBase<DerivedE> &E,
const Eigen::MatrixBase<DerivedoE> &oE,
Eigen::PlainObjectBase<DerivedvecParallel> &vecParallel,
Eigen::PlainObjectBase<DerivedvecPerpendicular> &vecPerpendicular)
{
using Scalar = typename DerivedvecParallel::Scalar;
using MatX = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
assert(E.rows()==F.rows() && "E does not match dimensions of F.");
assert(oE.rows()==F.rows() && "oE does not match dimensions of F.");
assert(E.cols()==3 && F.cols()==3 && oE.cols()==3 &&
"This method is for triangle meshes.");
assert(F.maxCoeff()<V.rows() && "V does not seem to belong to F.");
const typename DerivedE::Scalar m = E.maxCoeff()+1;
//Compute edge-based normal
MatX N, edgeN(m, 3);
edgeN.setZero();
per_face_normals(V, F, N);
for(Eigen::Index i=0; i<E.rows(); ++i) {
for(int j=0; j<3; ++j) {
edgeN.row(E(i,j)) += N.row(i);
}
}
edgeN.rowwise().normalize();
//Compute edge vectors
vecParallel.resize(m, 3);
if(computePerpendicular) { //This should ideally be an if constexpr
vecPerpendicular.resize(m, 3);
}
for(Eigen::Index i=0; i<E.rows(); ++i) {
for(int j=0; j<3; ++j) {
if(oE(i,j)<0) {
continue;
}
const typename DerivedE::Scalar e=E(i,j);
const typename DerivedF::Scalar vi=F(i,(j+1)%3), vj=F(i,(j+2)%3);
vecParallel.row(e) = (V.row(vj)-V.row(vi)).normalized();
if(computePerpendicular) { //This should ideally be an if constexpr
vecPerpendicular.row(e) =
Eigen::AngleAxis<Scalar>(0.5*PI, edgeN.row(e)) *
vecParallel.row(e).transpose();
}
}
}
}
#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
template void igl::edge_vectors<true, 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::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
#endif
include/igl/edge_vectors.h
+68
View File
@@ -0,0 +1,68 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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_EDGE_VECTORS_H
#define IGL_EDGE_VECTORS_H
#include "igl_inline.h"
#include <Eigen/Dense>
namespace igl
{
// Computes the normalized edge vectors for edges in a triangle mesh
//
// Input:
// V, F: triangle mesh
// E, oE: mapping from halfedges to edges and orientation as generated by
// orient_halfedges
//
// Output:
// vec: normalized edge vectors for each unique edge in E, according to order
// of E.
template<typename DerivedV,typename DerivedF,typename DerivedE,
typename DerivedoE, typename Derivedvec>
IGL_INLINE void edge_vectors(
const Eigen::MatrixBase<DerivedV> &V,
const Eigen::MatrixBase<DerivedF> &F,
const Eigen::MatrixBase<DerivedE> &E,
const Eigen::MatrixBase<DerivedoE> &oE,
Eigen::PlainObjectBase<Derivedvec> &vec);
// Computes the normalized edge vectors for edges in a triangle mesh
//
// Input:
// V, F: triangle mesh
// E, oE: mapping from halfedges to edges and orientation as generated by
// orient_halfedges
// template parameter computePerpendicular: whether to compute
// vecPerpendicular or not.
//
// Output:
// vecParallel: normalized edge vectors for each unique edge in E, according
// to order of E.
// vecPerpendicular: tangent unit perpendicular vector to each edge, according
// to orientation in oE, corresponds to each vector in
// vecParallel rotated by pi/2 around an edge-based normal.
//
template<bool computePerpendicular=true,
typename DerivedV,typename DerivedF,typename DerivedE,
typename DerivedoE, typename DerivedvecParallel,
typename DerivedvecPerpendicular>
IGL_INLINE void edge_vectors(
const Eigen::MatrixBase<DerivedV> &V,
const Eigen::MatrixBase<DerivedF> &F,
const Eigen::MatrixBase<DerivedE> &E,
const Eigen::MatrixBase<DerivedoE> &oE,
Eigen::PlainObjectBase<DerivedvecParallel> &vecParallel,
Eigen::PlainObjectBase<DerivedvecPerpendicular> &vecPerpendicular);
}
#ifndef IGL_STATIC_LIBRARY
# include "edge_vectors.cpp"
#endif
#endif
include/igl/internal_angles.cpp
+1
View File
@@ -101,4 +101,5 @@ template void igl::internal_angles<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eig
template void igl::internal_angles<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
template void igl::internal_angles<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
template void igl::internal_angles_using_squared_edge_lengths<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
template void igl::internal_angles_using_squared_edge_lengths<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
#endif
include/igl/min_quad_with_fixed.cpp
+1
View File
@@ -600,4 +600,5 @@ template bool igl::min_quad_with_fixed_solve<double, Eigen::Matrix<double, -1, 1
template bool igl::min_quad_with_fixed_solve<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::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
template bool igl::min_quad_with_fixed_solve<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::Matrix<double, -1, -1, 0, -1, -1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
template bool igl::min_quad_with_fixed<double, Eigen::Matrix<int, -1, -1, 0, -1, -1>, 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::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
template bool igl::min_quad_with_fixed<double, Eigen::Matrix<int, -1, 1, 0, -1, 1>, 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::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
#endif
include/igl/orient_halfedges.cpp
+48
View File
@@ -0,0 +1,48 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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 "orient_halfedges.h"
#include "oriented_facets.h"
#include "unique_simplices.h"
template <typename DerivedF, typename DerivedE, typename DerivedOE>
IGL_INLINE void
igl::orient_halfedges(
const Eigen::MatrixBase<DerivedF>& F,
Eigen::PlainObjectBase<DerivedE>& E,
Eigen::PlainObjectBase<DerivedOE>& oE)
{
assert(F.cols()==3 && "This only works for triangle meshes.");
using Int = typename DerivedF::Scalar;
const Eigen::Index m = F.rows();
DerivedE allE, EE;
oriented_facets(F, allE);
Eigen::Matrix<Int, Eigen::Dynamic, 1> IA, IC;
unique_simplices(allE, EE, IA, IC);
E.resize(m, 3);
oE.resize(m, 3);
for(Eigen::Index f=0; f<m; ++f) {
for(int e=0; e<3; ++e) {
const Int ind = f + m*e;
E(f,e) = IC(ind);
assert((EE(E(f,e),0)==allE(ind,0) || EE(E(f,e),0)==allE(ind,1)) &&
"Something is wrong in the edge matrix.");
oE(f,e) = EE(E(f,e),0)==allE(ind,0) ? 1 : -1;
}
}
}
#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
template void igl::orient_halfedges<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::MatrixBase<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/orient_halfedges.h
+41
View File
@@ -0,0 +1,41 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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_ORIENT_HALFEDGES_H
#define IGL_ORIENT_HALFEDGES_H
#include "igl_inline.h"
#include <Eigen/Core>
namespace igl
{
// Orients halfedges for a triangle mesh, assigning them to a unique edge.
//
// Inputs:
// F: input mesh connectivity
//
// Outputs:
// E: a mapping from each halfedge to each edge
// oE: the orientation of each halfedge compared to the orientation of the
// actual edge. Every edge appears positively oriented exactly once.
template <typename DerivedF, typename DerivedE, typename DerivedOE>
IGL_INLINE void
orient_halfedges(
const Eigen::MatrixBase<DerivedF>& F,
Eigen::PlainObjectBase<DerivedE>& E,
Eigen::PlainObjectBase<DerivedOE>& oE);
}
#ifndef IGL_STATIC_LIBRARY
# include "orient_halfedges.cpp"
#endif
#endif
include/igl/scalar_to_cr_vector_gradient.cpp
+116
View File
@@ -0,0 +1,116 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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 "scalar_to_cr_vector_gradient.h"
#include "orient_halfedges.h"
#include "doublearea.h"
#include "squared_edge_lengths.h"
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarG>
IGL_INLINE void
igl::scalar_to_cr_vector_gradient(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarG>& G)
{
Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic>
l_sq;
squared_edge_lengths(V, F, l_sq);
scalar_to_cr_vector_gradient_intrinsic(F, l_sq, E, oE, G);
}
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarG>
IGL_INLINE void
igl::scalar_to_cr_vector_gradient(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
Eigen::PlainObjectBase<DerivedE>& E,
Eigen::PlainObjectBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarG>& G)
{
if(E.rows()!=F.rows() || E.cols()!=F.cols() || oE.rows()!=F.rows() ||
oE.cols()!=F.cols()) {
orient_halfedges(F, E, oE);
}
const Eigen::PlainObjectBase<DerivedE>& cE = E;
const Eigen::PlainObjectBase<DerivedOE>& coE = oE;
scalar_to_cr_vector_gradient(V, F, cE, coE, G);
}
template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
typename DerivedOE, typename ScalarG>
IGL_INLINE void
igl::scalar_to_cr_vector_gradient_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarG>& G)
{
Eigen::Matrix<typename DerivedL_sq::Scalar, Eigen::Dynamic, Eigen::Dynamic>
dA;
DerivedL_sq l_sqrt = l_sq.array().sqrt().matrix();
doublearea(l_sqrt, dA);
scalar_to_cr_vector_gradient_intrinsic(F, l_sq, dA, E, oE, G);
}
template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
typename DerivedE, typename DerivedOE, typename ScalarG>
IGL_INLINE void
igl::scalar_to_cr_vector_gradient_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DeriveddA>& dA,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarG>& G)
{
assert(F.cols()==3 && "Faces have three vertices");
assert(E.rows()==F.rows() && E.cols()==F.cols() && oE.rows()==F.rows() &&
oE.cols()==F.cols() && "Wrong dimension in edge vectors");
const Eigen::Index m = F.rows();
const typename DerivedF::Scalar n = F.maxCoeff() + 1;
const typename DerivedE::Scalar nE = E.maxCoeff() + 1;
std::vector<Eigen::Triplet<ScalarG> > tripletList;
tripletList.reserve(5*3*m);
for(Eigen::Index f=0; f<m; ++f) {
for(int e=0; e<3; ++e) {
const typename DerivedF::Scalar i=F(f,(e+1)%3), j=F(f,(e+2)%3), k=F(f,e);
const ScalarG o=oE(f,e),
eij=l_sq(f,e), ejk=l_sq(f,(e+1)%3), eki=l_sq(f,(e+2)%3); //These are squared quantities.
const ScalarG s_eij = sqrt(eij);
tripletList.emplace_back(E(f,e), i, -o*dA(f)/(6.*s_eij));
tripletList.emplace_back(E(f,e)+nE, i, -o*(eij+ejk-eki)/(12.*s_eij));
tripletList.emplace_back(E(f,e), j, o*dA(f)/(6.*s_eij));
tripletList.emplace_back(E(f,e)+nE, j, -o*(eij-ejk+eki)/(12.*s_eij));
tripletList.emplace_back(E(f,e)+nE, k, o*s_eij/6.);
}
}
G.resize(2*nE, n);
G.setFromTriplets(tripletList.begin(), tripletList.end());
}
#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
template void igl::scalar_to_cr_vector_gradient_intrinsic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -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/scalar_to_cr_vector_gradient.h
+100
View File
@@ -0,0 +1,100 @@
// This file is part of libigl, a simple c++ geometry processing library.
//
// Copyright (C) 2020 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_SCALAR_TO_CR_VECTOT_GRADIENT_H
#define IGL_SCALAR_TO_CR_VECTOT_GRADIENT_H
#include "igl_inline.h"
#include <Eigen/Core>
#include <Eigen/Sparse>
namespace igl
{
// Computes the gradient matrix with hat functions on the right, and
// vector CR functions on the left.
// See Oded Stein, Max Wardetzky, Alec Jacobson, Eitan Grinspun, 2020.
// "A Simple Discretization of the Vector Dirichlet Energy"
//
// Inputs:
// V, F: input mesh
// E: a mapping from each halfedge to each edge, as computed with
// orient_halfedges.
// will be computed if not provided.
// oE: the orientation of each halfedge compared to the orientation of the
// actual edge, as computed with orient_halfedges.
// will be computed if not provided.
//
// Outputs:
// G: computed gradient matrix
// E, oE: these are computed if they are not present, as described above
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarG>
IGL_INLINE void
scalar_to_cr_vector_gradient(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarG>& G);
template <typename DerivedV, typename DerivedF, typename DerivedE,
typename DerivedOE, typename ScalarG>
IGL_INLINE void
scalar_to_cr_vector_gradient(
const Eigen::MatrixBase<DerivedV>& V,
const Eigen::MatrixBase<DerivedF>& F,
Eigen::PlainObjectBase<DerivedE>& E,
Eigen::PlainObjectBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarG>& G);
// Version that uses intrinsic quantities as input
//
// Inputs:
// F: input mesh connectivity
// l_sq: squared edge lengths of each halfedge
// dA: double area of each face
// E: a mapping from each halfedge to each edge.
// oE: the orientation of each halfedge compared to the orientation of the
// actual edge.
//
// Outputs:
// G: computed gradient matrix
template <typename DerivedF, typename DerivedL_sq, typename DerivedE,
typename DerivedOE, typename ScalarG>
IGL_INLINE void
scalar_to_cr_vector_gradient_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarG>& G);
template <typename DerivedF, typename DerivedL_sq, typename DeriveddA,
typename DerivedE, typename DerivedOE, typename ScalarG>
IGL_INLINE void
scalar_to_cr_vector_gradient_intrinsic(
const Eigen::MatrixBase<DerivedF>& F,
const Eigen::MatrixBase<DerivedL_sq>& l_sq,
const Eigen::MatrixBase<DeriveddA>& dA,
const Eigen::MatrixBase<DerivedE>& E,
const Eigen::MatrixBase<DerivedOE>& oE,
Eigen::SparseMatrix<ScalarG>& G);
}
#ifndef IGL_STATIC_LIBRARY
# include "scalar_to_cr_vector_gradient.cpp"
#endif
#endif
tests/include/igl/cr_vector_curvature_correction.cpp
+84
View File
File diff suppressed because one or more lines are too long
tests/include/igl/cr_vector_laplacian.cpp
+37
View File
File diff suppressed because one or more lines are too long
tests/include/igl/curved_hessian_energy.cpp
+46
View File
File diff suppressed because one or more lines are too long
tests/include/igl/orient_halfedges.cpp
+102
View File
@@ -0,0 +1,102 @@
#include <test_common.h>
#include <igl/orient_halfedges.h>
#include <igl/is_border_vertex.h>
#include <igl/edges.h>
#include <igl/remove_unreferenced.h>
#include <igl/triangle_triangle_adjacency.h>
#include <igl/EPS.h>
#include <vector>
TEST_CASE("orient_halfedges: sanity checks", "[igl]")
{
const auto meshes = test_common::manifold_meshes();
for(const auto& mesh : meshes) {
Eigen::MatrixXd V;
Eigen::MatrixXi F;
igl::read_triangle_mesh(test_common::data_path(mesh), V, F);
Eigen::MatrixXi I;
igl::remove_unreferenced(Eigen::MatrixXd(V), Eigen::MatrixXi(F), V, F,
I);
Eigen::MatrixXi TT, TTi;
igl::triangle_triangle_adjacency(F, TT, TTi);
// Fix mis-match convention
{
Eigen::PermutationMatrix<3,3> perm(3);
perm.indices() = Eigen::Vector3i(1,2,0);
TT = (TT*perm).eval();
TTi = (TTi*perm).eval();
for(int i=0;i<TTi.rows();i++) {
for(int j=0;j<TTi.cols();j++) {
TTi(i,j)=TTi(i,j)==-1?-1:(TTi(i,j)+3-1)%3;
}
}
}
Eigen::MatrixXi E, oE;
igl::orient_halfedges(F, E, oE);
const int m = E.maxCoeff()+1;
std::vector<bool> b(m);
for(int i=0; i<F.rows(); ++i) {
for(int j=0; j<3; ++j) {
b[E(i,j)] = TT(i,j)<0;
}
}
int nb = 0;
for(int i=0; i<b.size(); ++i) {
if(b[i]) {
++nb;
}
}
//Perform a variety of sanity checks.
//Correct number of edges
Eigen::MatrixXi sanityEdges;
igl::edges(F, sanityEdges);
REQUIRE(m == sanityEdges.rows());
//All border halfedges edges have orientation 1 only, all others
// orientations 1 and -1, so oE must sum to the number of border
// edges.
REQUIRE(nb == oE.array().sum());
//Every border halfedge has orientation 1. Every other edge appeats
// with orientation 1 and -1.
std::vector<int> appeared1(m,0), appearedm1(m,0);
for(int i=0; i<F.rows(); ++i) {
for(int j=0; j<3; ++j) {
if(oE(i,j)==1) {
++appeared1[E(i,j)];
} else if(oE(i,j)==-1) {
++appearedm1[E(i,j)];
} else {
REQUIRE(false); //Only 1 and -1 should occur.
}
}
}
for(int i=0; i<m; ++i) {
REQUIRE(appeared1[i]==1);
if(b[i]) {
REQUIRE(appearedm1[i]==0);
} else {
REQUIRE(appearedm1[i]==1);
}
}
//Two opposite halfedges always map to the same edge
for(int i=0; i<F.rows(); ++i) {
for(int j=0; j<3; ++j) {
if(TT(i,j)>=0) {
REQUIRE(E(i,j) == E(TT(i,j),TTi(i,j)));
}
}
}
}
}
tutorial/712_DataSmoothing/main.cpp
+151 -75
View File
@@ -1,5 +1,6 @@
#include <igl/read_triangle_mesh.h>
#include <igl/hessian_energy.h>
#include <igl/curved_hessian_energy.h>
#include <igl/massmatrix.h>
#include <igl/cotmatrix.h>
#include <igl/isolines_map.h>
@@ -7,6 +8,7 @@
#include <igl/vertex_components.h>
#include <igl/remove_unreferenced.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/heat_geodesics.h>
#include <Eigen/Core>
#include <Eigen/SparseCholesky>
@@ -22,86 +24,160 @@
int main(int argc, char * argv[])
{
typedef Eigen::SparseMatrix<double> SparseMat;
//Read our mesh
Eigen::MatrixXd V;
Eigen::MatrixXi F;
if(!igl::read_triangle_mesh(
argc>1?argv[1]: TUTORIAL_SHARED_PATH "/beetle.off",V,F)) {
std::cout << "Failed to load mesh." << std::endl;
typedef Eigen::SparseMatrix<double> SparseMat;
srand(57);
//Read our mesh
Eigen::MatrixXd V;
Eigen::MatrixXi F;
if(!igl::read_triangle_mesh(
argc>1?argv[1]: TUTORIAL_SHARED_PATH "/beetle.off",V,F)) {
std::cout << "Failed to load mesh." << std::endl;
}
//Constructing an exact function to smooth
igl::HeatGeodesicsData<double> hgData;
igl::heat_geodesics_precompute(V, F, hgData);
Eigen::VectorXd heatDist;
Eigen::VectorXi gamma(1); gamma << 1947; //1631;
igl::heat_geodesics_solve(hgData, gamma, heatDist);
Eigen::VectorXd zexact =
0.1*(heatDist.array() + (-heatDist.maxCoeff())).pow(2)
+ 3*V.block(0,1,V.rows(),1).array().cos();
//Make the exact function noisy
const double s = 0.1*(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);
SparseMat QcH;
igl::curved_hessian_energy(V, F, QcH);
//Solve to find Laplacian-smoothed Hessian-smoothed, and
// curved-Hessian-smoothed solutions
const double al = 3e-7;
Eigen::SimplicialLDLT<SparseMat> lapSolver(al*QL + (1.-al)*M);
Eigen::VectorXd zl = lapSolver.solve(al*M*znoisy);
const double ah = 2e-7;
Eigen::SimplicialLDLT<SparseMat> hessSolver(ah*QH + (1.-ah)*M);
Eigen::VectorXd zh = hessSolver.solve(ah*M*znoisy);
const double ach = 3e-7;
Eigen::SimplicialLDLT<SparseMat> curvedHessSolver(al*QcH + (1.-ach)*M);
Eigen::VectorXd zch = curvedHessSolver.solve(ach*M*znoisy);
//Viewer that shows all functions: zexact, znoisy, zl, zh
igl::opengl::glfw::Viewer viewer;
viewer.data().set_mesh(V,F);
viewer.data().show_lines = false;
viewer.callback_key_down =
[&](igl::opengl::glfw::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;
case '5':
z = &zch;
break;
default:
return false;
}
//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::opengl::glfw::Viewer viewer;
viewer.data().set_mesh(V,F);
viewer.data().show_lines = false;
viewer.callback_key_down =
[&](igl::opengl::glfw::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;
}
viewer.data().set_data(*z);
return true;
};
std::cout << R"(Usage:
viewer.data().set_data(*z);
return true;
};
std::cout << R"(Smoothing a noisy function.
Usage:
1 Show original function
2 Show noisy function
3 Biharmonic smoothing (zero Neumann boundary)
4 Biharmonic smoothing (natural Hessian boundary)
4 Biharmonic smoothing (natural planar Hessian boundary)
5 Biharmonic smoothing (natural curved Hessian boundary)
)";
Eigen::MatrixXd CM;
igl::parula(Eigen::VectorXd::LinSpaced(21,0,1).eval(),false,CM);
igl::isolines_map(Eigen::MatrixXd(CM),CM);
viewer.data().set_colormap(CM);
viewer.data().set_data(znoisy);
viewer.launch();
Eigen::MatrixXd CM;
igl::parula(Eigen::VectorXd::LinSpaced(21,0,1).eval(),false,CM);
igl::isolines_map(Eigen::MatrixXd(CM),CM);
viewer.data().set_colormap(CM);
viewer.data().set_data(znoisy);
viewer.launch();
//Constructing a step function to smooth
Eigen::VectorXd zstep = Eigen::VectorXd::Zero(V.rows());
for(int i=0; i<V.rows(); ++i) {
zstep(i) = V(i,2)<-0.25 ? 1. : (V(i,2)>0.31 ? 2. : 0);
}
//Smooth that function
const double sl = 2e-5;
Eigen::SimplicialLDLT<SparseMat> stepLapSolver(sl*QL + (1.-sl)*M);
Eigen::VectorXd stepzl = stepLapSolver.solve(al*M*zstep);
const double sh = 6e-6;
Eigen::SimplicialLDLT<SparseMat> stepHessSolver(sh*QH + (1.-sh)*M);
Eigen::VectorXd stepzh = stepHessSolver.solve(ah*M*zstep);
const double sch = 2e-5;
Eigen::SimplicialLDLT<SparseMat> stepCurvedHessSolver(sl*QcH + (1.-sch)*M);
Eigen::VectorXd stepzch = stepCurvedHessSolver.solve(ach*M*zstep);
//Display functions
igl::opengl::glfw::Viewer viewer2;
viewer2.data().set_mesh(V,F);
viewer2.data().show_lines = false;
viewer2.callback_key_down =
[&](igl::opengl::glfw::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 = &zstep;
break;
case '2':
z = &stepzl;
break;
case '3':
z = &stepzh;
break;
case '4':
z = &stepzch;
break;
default:
return false;
}
viewer.data().set_data(*z);
return true;
};
std::cout << R"(Smoothing a step function.
Usage:
1 Show step function
2 Biharmonic smoothing (zero Neumann boundary)
3 Biharmonic smoothing (natural planar Hessian boundary)
4 Biharmonic smoothing (natural curved Hessian boundary)
return 0;
)";
viewer2.data().set_colormap(CM);
viewer2.data().set_data(zstep);
viewer2.launch();
return 0;
}
tutorial/721_VectorFieldSmoothing/CMakeLists.txt
+5
View File
@@ -0,0 +1,5 @@
get_filename_component(PROJECT_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
project(${PROJECT_NAME})
add_executable(${PROJECT_NAME}_bin main.cpp)
target_link_libraries(${PROJECT_NAME}_bin igl::core igl::opengl igl::opengl_glfw tutorials)
tutorial/721_VectorFieldSmoothing/main.cpp
Executable
+147
View File
@@ -0,0 +1,147 @@
#include <igl/read_triangle_mesh.h>
#include <igl/parula.h>
#include <igl/remove_unreferenced.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/per_face_normals.h>
#include <igl/orient_halfedges.h>
#include <igl/cr_vector_laplacian.h>
#include <igl/cr_vector_mass.h>
#include <igl/edge_midpoints.h>
#include <igl/edge_vectors.h>
#include <igl/average_from_edges_onto_vertices.h>
#include <igl/PI.h>
#include <Eigen/Core>
#include <Eigen/SparseCholesky>
#include <Eigen/Geometry>
#include <iostream>
#include <set>
#include <limits>
#include <stdlib.h>
#include "tutorial_shared_path.h"
int main(int argc, char * argv[])
{
typedef Eigen::SparseMatrix<double> SparseMat;
//Constants used for smoothing
const double howMuchToSmoothBy = 1e-1;
const int howManySmoothingInterations = 50;
//Read our mesh
Eigen::MatrixXd V;
Eigen::MatrixXi F;
if(!igl::read_triangle_mesh
(argc>1?argv[1]: TUTORIAL_SHARED_PATH "/elephant.obj",V,F)) {
std::cout << "Failed to load mesh." << std::endl;
}
//Orient edges for plotting
Eigen::MatrixXi E, oE;
igl::orient_halfedges(F, E, oE);
//Compute edge midpoints & edge vectors
Eigen::MatrixXd edgeMps, parVec, perpVec;
igl::edge_midpoints(V, F, E, oE, edgeMps);
igl::edge_vectors(V, F, E, oE, parVec, perpVec);
//Constructing a function to add noise to
const auto zraw_function = [] (const Eigen::Vector3d& x) {
return Eigen::Vector3d(0.2*x(1) + cos(2*x(1)+0.2),
0.5*x(0) + 0.15,
0.3*cos(0.2+igl::PI*x(2)));
};
Eigen::VectorXd zraw(2*edgeMps.rows());
for(int i=0; i<edgeMps.rows(); ++i) {
const Eigen::Vector3d f = zraw_function(edgeMps.row(i));
zraw(i) = f.dot(parVec.row(i));
zraw(i+edgeMps.rows()) = f.dot(perpVec.row(i));
}
//Add noise
srand(71);
const double l = 15;
Eigen::VectorXd znoisy = zraw + l*Eigen::VectorXd::Random(zraw.size());
//Denoise function using the vector Dirichlet energy
Eigen::VectorXd zsmoothed = znoisy;
for(int i=0; i<howManySmoothingInterations; ++i) {
//Compute Laplacian and mass matrix
SparseMat L, M;
igl::cr_vector_mass(V, F, E, oE, M);
igl::cr_vector_laplacian(V, F, E, oE, L);
//Implicit step
Eigen::SimplicialLDLT<SparseMat> rhsSolver(M + howMuchToSmoothBy*L);
zsmoothed = rhsSolver.solve(M*zsmoothed);
}
//Convert vector fields for plotting
const auto cr_result_to_vecs_and_colors = [&]
(const Eigen::VectorXd& z, Eigen::MatrixXd& vecs, Eigen::MatrixXd& colors) {
vecs.resize(edgeMps.rows(), 3);
for(int i=0; i<edgeMps.rows(); ++i) {
vecs.row(i) = z(i)*parVec.row(i)
+ z(i+edgeMps.rows())*perpVec.row(i);
}
igl::average_from_edges_onto_vertices
(F, E, oE, vecs.rowwise().norm(), colors);
};
Eigen::MatrixXd noisyvecs, noisycolors, smoothedvecs, smoothedcolors,
rawvecs, rawcolors;
cr_result_to_vecs_and_colors(znoisy, noisyvecs, noisycolors);
cr_result_to_vecs_and_colors(zsmoothed, smoothedvecs, smoothedcolors);
cr_result_to_vecs_and_colors(zraw, rawvecs, rawcolors);
//Viewer that shows noisy and denoised functions
igl::opengl::glfw::Viewer viewer;
viewer.data().set_mesh(V,F);
viewer.data().show_lines = false;
viewer.callback_key_down =
[&](igl::opengl::glfw::Viewer & viewer, unsigned char key, int mod)->bool
{
const Eigen::MatrixXd *vecs, *colors;
switch(key) {
case '1':
vecs = &rawvecs;
colors = &rawcolors;
break;
case '2':
vecs = &noisyvecs;
colors = &noisycolors;
break;
case '3':
vecs = &smoothedvecs;
colors = &smoothedcolors;
break;
default:
return false;
}
viewer.data().set_data(*colors);
viewer.data().clear_edges();
const double s = 0.08; //How much to scale vectors during plotting
viewer.data().add_edges(edgeMps, edgeMps + s*(*vecs),
Eigen::RowVector3d(0.1, 0.1, 0.1));
return true;
};
std::cout << R"(Usage:
1 Show raw function
2 Show noisy function
3 Show smoothed function
)";
Eigen::MatrixXd CM;
igl::parula(Eigen::VectorXd::LinSpaced
(500,znoisy.minCoeff(), znoisy.maxCoeff()).eval(), true, CM);
viewer.data().set_colormap(CM);
viewer.callback_key_down(viewer, '1', 0);
viewer.launch();
return 0;
}
tutorial/722_VectorParallelTransport/CMakeLists.txt
+5
View File
@@ -0,0 +1,5 @@
get_filename_component(PROJECT_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
project(${PROJECT_NAME})
add_executable(${PROJECT_NAME}_bin main.cpp)
target_link_libraries(${PROJECT_NAME}_bin igl::core igl::opengl igl::opengl_glfw tutorials)
tutorial/722_VectorParallelTransport/main.cpp
Executable
+150
View File
@@ -0,0 +1,150 @@
#include <igl/read_triangle_mesh.h>
#include <igl/parula.h>
#include <igl/remove_unreferenced.h>
#include <igl/opengl/glfw/Viewer.h>
#include <igl/per_face_normals.h>
#include <igl/orient_halfedges.h>
#include <igl/cr_vector_laplacian.h>
#include <igl/cr_vector_mass.h>
#include <igl/crouzeix_raviart_cotmatrix.h>
#include <igl/crouzeix_raviart_massmatrix.h>
#include <igl/edge_midpoints.h>
#include <igl/edge_vectors.h>
#include <igl/average_from_edges_onto_vertices.h>
#include <igl/min_quad_with_fixed.h>
#include <igl/heat_geodesics.h>
#include <Eigen/Core>
#include <Eigen/SparseCholesky>
#include <Eigen/Geometry>
#include <iostream>
#include <set>
#include <limits>
#include <stdlib.h>
#include "tutorial_shared_path.h"
int main(int argc, char * argv[])
{
typedef Eigen::SparseMatrix<double> SparseMat;
typedef Eigen::Matrix<double, 1, 1> Vector1d;
typedef Eigen::Matrix<int, 1, 1> Vector1i;
//Constants used for smoothing
const double howMuchToSmoothBy = 1e-1;
const int howManySmoothingInterations = 50;
//Read our mesh
Eigen::MatrixXd V;
Eigen::MatrixXi F;
if(!igl::read_triangle_mesh
(argc>1?argv[1]: TUTORIAL_SHARED_PATH "/cheburashka.off",V,F)) {
std::cout << "Failed to load mesh." << std::endl;
}
//Compute vector Laplacian and mass matrix
Eigen::MatrixXi E, oE;//Compute Laplacian and mass matrix
SparseMat vecL, vecM;
igl::cr_vector_mass(V, F, E, oE, vecM);
igl::cr_vector_laplacian(V, F, E, oE, vecL);
const int m = vecL.rows()/2; //The number of edges in the mesh
//Convert the E / oE matrix format to list of edges / EMAP format required
// by the functions constructing scalar Crouzeix-Raviart functions
Eigen::MatrixXi Elist(m,2), EMAP(3*F.rows(),1);
for(int i=0; i<F.rows(); ++i) {
for(int j=0; j<3; ++j) {
const int e = E(i,j);
EMAP(i+j*F.rows()) = e;
if(oE(i,j)>0) {
Elist.row(e) << F(i, (j+1)%3), F(i, (j+2)%3);
}
}
}
SparseMat scalarL, scalarM;
igl::crouzeix_raviart_massmatrix(V, F, Elist, EMAP, scalarM);
igl::crouzeix_raviart_cotmatrix(V, F, Elist, EMAP, scalarL);
//Compute edge midpoints & edge vectors
Eigen::MatrixXd edgeMps, parVec, perpVec;
igl::edge_midpoints(V, F, E, oE, edgeMps);
igl::edge_vectors(V, F, E, oE, parVec, perpVec);
//Perform the vector heat method
const int initialIndex = 14319;
const double initialPara=0.95, initialPerp=0.08;
const double t = 0.01;
SparseMat Aeq;
Eigen::VectorXd Beq;
Eigen::VectorXi known = Eigen::Vector2i(initialIndex, initialIndex+m);
Eigen::VectorXd knownVals = Eigen::Vector2d(initialPara, initialPerp);
Eigen::VectorXd Y0 = Eigen::VectorXd::Zero(2*m), Yt;
Y0(initialIndex) = initialPara; Y0(initialIndex+m) = initialPerp;
igl::min_quad_with_fixed
(SparseMat(vecM+t*vecL), Eigen::VectorXd(-vecM*Y0), known, knownVals,
Aeq, Beq, false, Yt);
Eigen::VectorXd u0 = Eigen::VectorXd::Zero(m), ut;
u0(initialIndex) = sqrt(initialPara*initialPara + initialPerp*initialPerp);
Eigen::VectorXi knownScal = Vector1i(initialIndex);
Eigen::VectorXd knownScalVals = Vector1d(u0(initialIndex));
igl::min_quad_with_fixed
(SparseMat(scalarM+t*scalarL), Eigen::VectorXd(-scalarM*u0), knownScal,
knownScalVals, Aeq, Beq, false, ut);
Eigen::VectorXd phi0 = Eigen::VectorXd::Zero(m), phit;
phi0(initialIndex) = 1;
Eigen::VectorXd knownScalValsPhi = Vector1d(1);
igl::min_quad_with_fixed
(SparseMat(scalarM+t*scalarL), Eigen::VectorXd(-scalarM*phi0), knownScal,
knownScalValsPhi, Aeq, Beq, false, phit);
Eigen::ArrayXd Xtfactor = ut.array() /
(phit.array() * (Yt.array().segment(0,m)*Yt.array().segment(0,m)
+ Yt.array().segment(m,m)*Yt.array().segment(m,m)).sqrt());
Eigen::VectorXd Xt(2*m);
Xt.segment(0,m) = Xtfactor * Yt.segment(0,m).array();
Xt.segment(m,m) = Xtfactor * Yt.segment(m,m).array();
//Compute scalar heat colors
igl::HeatGeodesicsData<double> hgData;
igl::heat_geodesics_precompute(V, F, hgData);
Eigen::VectorXd heatColor;
Eigen::VectorXi gamma = Elist.row(initialIndex);
igl::heat_geodesics_solve(hgData, gamma, heatColor);
//Convert vector field for plotting
Eigen::MatrixXd vecs(m, 3);
for(int i=0; i<edgeMps.rows(); ++i) {
vecs.row(i) = Xt(i)*parVec.row(i) + Xt(i+edgeMps.rows())*perpVec.row(i);
}
//Viewer that shows parallel transported vector
igl::opengl::glfw::Viewer viewer;
viewer.data().set_mesh(V,F);
viewer.data().show_lines = false;
viewer.data().set_data(heatColor.maxCoeff()-heatColor.array(), //invert colormap
igl::COLOR_MAP_TYPE_VIRIDIS);
const double s = 0.012; //How much to scale vectors during plotting
Eigen::MatrixXd vecColors(m, 3);
for(int i=0; i<m; ++i) {
vecColors.row(i) << 0.1, 0.1, 0.1;
}
vecColors.row(initialIndex) << 0.9, 0.1, 0.1;
viewer.data().add_edges(edgeMps, edgeMps + s*vecs, vecColors);
std::cout << R"(The red vector is parallel transported to every point on the surface.
The surface is shaded by geodesic distance from the red vector.
)"
<< std::endl;
viewer.launch();
return 0;
}
tutorial/CMakeLists.txt
+2
View File
@@ -166,4 +166,6 @@ if(TUTORIALS_CHAPTER7)
add_subdirectory("718_IterativeClosestPoint")
add_subdirectory("719_ExplodedView")
add_subdirectory("720_BlueNoise")
add_subdirectory("721_VectorFieldSmoothing")
add_subdirectory("722_VectorParallelTransport")
endif()