manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Cleaned up harmonic.h (added standard header, needed meshassert)

This commit is contained in:
Paolo Cignoni 2014-11-03 15:00:06 +00:00
parent e1a327e556
commit db53a1ff06
1 changed files with 219 additions and 199 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

418
vcg/complex/algorithms/harmonic.h
View File

@ -1,11 +1,29 @@
/****************************************************************************
* VCGLib o o *
* Visual and Computer Graphics Library o o *
* _ O _ *
* Copyright(C) 2014 \/)\/ *
* Visual Computing Lab /\/| *
* ISTI - Italian National Research Council | *
* \ *
* All rights reserved. *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
* for more details. *
* *
****************************************************************************/
#ifndef __VCGLIB_HARMONIC_FIELD #ifndef __VCGLIB_HARMONIC_FIELD
#define __VCGLIB_HARMONIC_FIELD #define __VCGLIB_HARMONIC_FIELD
#include <vcg/complex/complex.h> #include <vcg/complex/complex.h>
#include <utility>
#include <vector>
#include <map>
#include <eigenlib/Eigen/Sparse> #include <eigenlib/Eigen/Sparse>
namespace vcg { namespace vcg {
@ -15,237 +33,239 @@ template <class MeshType, typename Scalar = double>
class Harmonic class Harmonic
{ {
public: public:
typedef typename MeshType::VertexType VertexType; typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::FaceType FaceType; typedef typename MeshType::FaceType FaceType;
typedef typename MeshType::CoordType CoordType; typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::ScalarType ScalarType; typedef typename MeshType::ScalarType ScalarType;
typedef double CoeffScalar; typedef double CoeffScalar;
typedef typename std::pair<VertexType *, Scalar> Constraint; typedef typename std::pair<VertexType *, Scalar> Constraint;
typedef typename std::vector<Constraint> ConstraintVec; typedef typename std::vector<Constraint> ConstraintVec;
typedef typename ConstraintVec::const_iterator ConstraintIt; typedef typename ConstraintVec::const_iterator ConstraintIt;
/** /**
* @brief ComputeScalarField * @brief ComputeScalarField
* Generates a scalar harmonic field over the mesh. * Generates a scalar harmonic field over the mesh.
* For more details see:\n Kai Xua, Hao Zhang, Daniel Cohen-Or, Yueshan Xionga,'Dynamic Harmonic Fields for Surface Processing'.\nin Computers & Graphics, 2009 * For more details see:\n Kai Xua, Hao Zhang, Daniel Cohen-Or, Yueshan Xionga,'Dynamic Harmonic Fields for Surface Processing'.\nin Computers & Graphics, 2009
* @param m the mesh * @param m the mesh
* @param constraints the Dirichlet boundary conditions in the form of vector of pairs <vertex pointer, value>. * @param constraints the Dirichlet boundary conditions in the form of vector of pairs <vertex pointer, value>.
* @param field the accessor to use to write the computed per-vertex values (must have the [ ] operator). * @param field the accessor to use to write the computed per-vertex values (must have the [ ] operator).
* @return true if the algorithm succeeds, false otherwise. * @return true if the algorithm succeeds, false otherwise.
* @note the algorithm has unexpected behavior if the mesh contains unreferenced vertices. * @note the algorithm has unexpected behavior if the mesh contains unreferenced vertices.
*/ */
template <typename ACCESSOR> template <typename ACCESSOR>
static bool ComputeScalarField(MeshType & m, const ConstraintVec & constraints, ACCESSOR field) static bool ComputeScalarField(MeshType & m, const ConstraintVec & constraints, ACCESSOR field)
{ {
typedef Eigen::SparseMatrix<CoeffScalar> SpMat; // sparse matrix type typedef Eigen::SparseMatrix<CoeffScalar> SpMat; // sparse matrix type
typedef Eigen::Triplet<CoeffScalar> Triple; // triplet type to fill the matrix typedef Eigen::Triplet<CoeffScalar> Triple; // triplet type to fill the matrix
RequirePerVertexFlags(m); RequirePerVertexFlags(m);
RequireCompactness(m); RequireCompactness(m);
RequireFFAdjacency(m); RequireFFAdjacency(m);
MeshAssert<MeshType>::FFAdjacencyIsInitialized(m);
MeshAssert<MeshType>::NoUnreferencedVertex(m);
if (constraints.empty()) if (constraints.empty())
return false; return false;
int n = m.VN(); int n = m.VN();
// Generate coefficients // Generate coefficients
std::vector<Triple> coeffs; // coefficients of the system std::vector<Triple> coeffs; // coefficients of the system
std::map<size_t,CoeffScalar> sums; // row sum of the coefficient matrix std::map<size_t,CoeffScalar> sums; // row sum of the coefficient matrix
vcg::tri::UpdateFlags<MeshType>::FaceClearV(m); vcg::tri::UpdateFlags<MeshType>::FaceClearV(m);
for (size_t i = 0; i < m.face.size(); ++i) for (size_t i = 0; i < m.face.size(); ++i)
{ {
FaceType & f = m.face[i]; FaceType & f = m.face[i];
assert(!f.IsV());
assert(!f.IsD()); f.SetV();
assert(!f.IsV());
f.SetV(); // Generate coefficients for each edge
for (int edge = 0; edge < 3; ++edge)
{
CoeffScalar weight;
WeightInfo res = CotangentWeightIfNotVisited(f, edge, weight);
// Generate coefficients for each edge if (res == EdgeAlreadyVisited) continue;
for (int edge = 0; edge < 3; ++edge) assert(res == Success);
{
CoeffScalar weight;
WeightInfo res = CotangentWeightIfNotVisited(f, edge, weight);
if (res == EdgeAlreadyVisited) continue; // Add the weight to the coefficients vector for both the vertices of the considered edge
assert(res == Success); size_t v0_idx = vcg::tri::Index(m, f.V0(edge));
size_t v1_idx = vcg::tri::Index(m, f.V1(edge));
// Add the weight to the coefficients vector for both the vertices of the considered edge coeffs.push_back(Triple(v0_idx, v1_idx, -weight));
size_t v0_idx = vcg::tri::Index(m, f.V0(edge)); coeffs.push_back(Triple(v1_idx, v0_idx, -weight));
size_t v1_idx = vcg::tri::Index(m, f.V1(edge));
coeffs.push_back(Triple(v0_idx, v1_idx, -weight)); // Add the weight to the row sum
coeffs.push_back(Triple(v1_idx, v0_idx, -weight)); sums[v0_idx] += weight;
sums[v1_idx] += weight;
}
}
// Add the weight to the row sum // Setup the system matrix
sums[v0_idx] += weight; SpMat laplaceMat; // the system to be solved
sums[v1_idx] += weight; laplaceMat.resize(n, n); // eigen initializes it to zero
} laplaceMat.reserve(coeffs.size());
} for (std::map<size_t,CoeffScalar>::const_iterator it = sums.begin(); it != sums.end(); ++it)
{
// Setup the system matrix coeffs.push_back(Triple(it->first, it->first, it->second));
SpMat laplaceMat; // the system to be solved }
laplaceMat.resize(n, n); // eigen initializes it to zero laplaceMat.setFromTriplets(coeffs.begin(), coeffs.end());
laplaceMat.reserve(coeffs.size());
for (std::map<size_t,CoeffScalar>::const_iterator it = sums.begin(); it != sums.end(); ++it)
{
coeffs.push_back(Triple(it->first, it->first, it->second));
}
laplaceMat.setFromTriplets(coeffs.begin(), coeffs.end());
// Setting the constraints // Setting the constraints
const CoeffScalar alpha = pow(10, 8); // penalty factor alpha const CoeffScalar alpha = pow(10.0, 8.0); // penalty factor alpha
Eigen::Matrix<CoeffScalar, Eigen::Dynamic, 1> b, x; // Unknown and known terms vectors // const CoeffScalar alpha = CoeffScalar(1e5); // penalty factor alpha
b.setZero(n);
for (ConstraintIt it=constraints.begin(); it!=constraints.end(); it++) Eigen::Matrix<CoeffScalar, Eigen::Dynamic, 1> b, x; // Unknown and known terms vectors
{ b.setZero(n);
size_t v_idx = vcg::tri::Index(m, it->first);
b(v_idx) = alpha * it->second;
laplaceMat.coeffRef(v_idx, v_idx) += alpha;
}
// Perform matrix decomposition for (ConstraintIt it=constraints.begin(); it!=constraints.end(); it++)
Eigen::SimplicialLDLT<SpMat> solver; {
solver.compute(laplaceMat); size_t v_idx = vcg::tri::Index(m, it->first);
// TODO eventually use another solver (e.g. CHOLMOD for dynamic setups) b(v_idx) = alpha * it->second;
if(solver.info() != Eigen::Success) laplaceMat.coeffRef(v_idx, v_idx) += alpha;
{ }
// decomposition failed
switch (solver.info())
{
// possible errors
case Eigen::NumericalIssue :
case Eigen::NoConvergence :
case Eigen::InvalidInput :
default : return false;
}
}
// Solve the system: laplacianMat x = b // Perform matrix decomposition
x = solver.solve(b); Eigen::SimplicialLDLT<SpMat> solver;
if(solver.info() != Eigen::Success) solver.compute(laplaceMat);
{ // TODO eventually use another solver (e.g. CHOLMOD for dynamic setups)
return false; if(solver.info() != Eigen::Success)
} {
// decomposition failed
switch (solver.info())
{
// possible errors
case Eigen::NumericalIssue :
case Eigen::NoConvergence :
case Eigen::InvalidInput :
default : return false;
}
}
// Set field value using the provided handle // Solve the system: laplacianMat x = b
for (size_t i = 0; int(i) < n; ++i) x = solver.solve(b);
{ if(solver.info() != Eigen::Success)
field[i] = Scalar(x[i]); {
} return false;
}
return true; // Set field value using the provided handle
} for (size_t i = 0; int(i) < n; ++i)
{
field[i] = Scalar(x[i]);
}
enum WeightInfo return true;
{ }
Success = 0,
EdgeAlreadyVisited enum WeightInfo
}; {
Success = 0,
EdgeAlreadyVisited
};
/** /**
* @brief CotangentWeightIfNotVisited computes the cotangent weighting for an edge * @brief CotangentWeightIfNotVisited computes the cotangent weighting for an edge
* (if it has not be done yet). * (if it has not be done yet).
* This must be ensured setting the visited flag on the face once all edge weights have been computed. * This must be ensured setting the visited flag on the face once all edge weights have been computed.
* @param f the face * @param f the face
* @param edge the edge of the provided face for which we compute the weight * @param edge the edge of the provided face for which we compute the weight
* @param weight the computed weight (output) * @param weight the computed weight (output)
* @return Success if everything is fine, EdgeAlreadyVisited if the weight * @return Success if everything is fine, EdgeAlreadyVisited if the weight
* for the considered edge has been already computed. * for the considered edge has been already computed.
* @note the mesh must have the face-face topology updated * @note the mesh must have the face-face topology updated
*/ */
template <typename ScalarT> template <typename ScalarT>
static WeightInfo CotangentWeightIfNotVisited(const FaceType &f, int edge, ScalarT & weight) static WeightInfo CotangentWeightIfNotVisited(const FaceType &f, int edge, ScalarT & weight)
{ {
const FaceType * fp = f.cFFp(edge); const FaceType * fp = f.cFFp(edge);
if (fp != NULL && fp != &f) if (fp != NULL && fp != &f)
{ {
// not a border edge // not a border edge
if (fp->IsV()) return EdgeAlreadyVisited; if (fp->IsV()) return EdgeAlreadyVisited;
} }
weight = CotangentWeight<ScalarT>(f, edge); weight = CotangentWeight<ScalarT>(f, edge);
return Success; return Success;
} }
/** /**
* @brief ComputeWeight computes the cotangent weighting for an edge * @brief ComputeWeight computes the cotangent weighting for an edge
* @param f the face * @param f the face
* @param edge the edge of the provided face for which we compute the weight * @param edge the edge of the provided face for which we compute the weight
* @return the computed weight * @return the computed weight
* @note the mesh must have the face-face topology updated * @note the mesh must have the face-face topology updated
*/ */
template <typename ScalarT> template <typename ScalarT>
static ScalarT CotangentWeight(const FaceType &f, int edge) static ScalarT CotangentWeight(const FaceType &f, int edge)
{ {
assert(edge >= 0 && edge < 3); assert(edge >= 0 && edge < 3);
// get the adjacent face // get the adjacent face
const FaceType * fp = f.cFFp(edge); const FaceType * fp = f.cFFp(edge);
// v0 // v0
// /|\ // /|\
// / | \ // / | \
// / | \ // / | \
// / | \ // / | \
// va\ | /vb // va\ | /vb
// \ | / // \ | /
// \ | / // \ | /
// \|/ // \|/
// v1 // v1
ScalarT cotA = 0; ScalarT cotA = 0;
ScalarT cotB = 0; ScalarT cotB = 0;
// Get the edge (a pair of vertices) // Get the edge (a pair of vertices)
VertexType * v0 = f.cV0(edge); VertexType * v0 = f.cV0(edge);
VertexType * v1 = f.cV1(edge); VertexType * v1 = f.cV1(edge);
if (fp != NULL && if (fp != NULL &&
fp != &f) fp != &f)
{ {
// not a border edge // not a border edge
VertexType * vb = fp->cV2(f.cFFi(edge)); VertexType * vb = fp->cV2(f.cFFi(edge));
ScalarT angleB = ComputeAngle<ScalarT>(v0, vb, v1); ScalarT angleB = ComputeAngle<ScalarT>(v0, vb, v1);
cotB = vcg::math::Cos(angleB) / vcg::math::Sin(angleB); cotB = vcg::math::Cos(angleB) / vcg::math::Sin(angleB);
} }
VertexType * va = f.cV2(edge); VertexType * va = f.cV2(edge);
ScalarT angleA = ComputeAngle<ScalarT>(v0, va, v1); ScalarT angleA = ComputeAngle<ScalarT>(v0, va, v1);
cotA = vcg::math::Cos(angleA) / vcg::math::Sin(angleA); cotA = vcg::math::Cos(angleA) / vcg::math::Sin(angleA);
return (cotA + cotB) / 2; return (cotA + cotB) / 2;
} }
template <typename ScalarT> template <typename ScalarT>
static ScalarT ComputeAngle(const VertexType * a, const VertexType * b, const VertexType * c) static ScalarT ComputeAngle(const VertexType * a, const VertexType * b, const VertexType * c)
{ {
// a // a
// / // /
// / // /
// / // /
// / ___ compute the angle in b // / ___ compute the angle in b
// b \ // b \
// \ // \
// \ // \
// \ // \
// c // c
assert(a != NULL && b != NULL && c != NULL); assert(a != NULL && b != NULL && c != NULL);
Point3<ScalarT> A,B,C; Point3<ScalarT> A,B,C;
A.Import(a->P()); A.Import(a->P());
B.Import(b->P()); B.Import(b->P());
C.Import(c->P()); C.Import(c->P());
ScalarT angle = vcg::Angle(A - B, C - B); ScalarT angle = vcg::Angle(A - B, C - B);
return angle; return angle;
} }
}; };
} }