Program Listing for File SurfaceMesh.h¶
↰ Return to documentation for file (include/gamer/SurfaceMesh.h)
/*
* ***************************************************************************
* This file is part of the GAMer software.
* Copyright (C) 2016-2018
* by Christopher Lee, John Moody, Rommie Amaro, J. Andrew McCammon,
* and Michael Holst
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library 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
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
* ***************************************************************************
*/
#pragma once
#include <iostream>
#include <stdexcept>
#include <string>
#include <memory>
#include <unordered_set>
#include <utility>
#include <casc/casc>
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include "gamer/Vertex.h"
namespace gamer
{
struct SMGlobal
{
int marker;
float volumeConstraint;
bool useVolumeConstraint;
bool ishole;
SMGlobal(int marker = -1, float volumeConstraint = -1, bool useVolumeConstraint = false, bool ishole = false) :
marker(marker), volumeConstraint(volumeConstraint), useVolumeConstraint(useVolumeConstraint), ishole(ishole) {
}
};
struct SMVertex : Vertex {
Vector normal;
using Vertex::Vertex;
};
struct SMEdge {
bool selected;
SMEdge() : SMEdge(0) {
}
SMEdge(bool select) : selected(select) {
}
};
struct SMFaceProperties
{
int marker;
bool selected;
Vector normal;
SMFaceProperties(int marker, bool selected) :
marker(marker), selected(selected) {
}
};
struct SMFace : casc::Orientable, SMFaceProperties
{
SMFace() : SMFace(Orientable{0}, SMFaceProperties{-1, false}) {
}
SMFace(int marker, bool selected) : SMFace(Orientable{0}, SMFaceProperties{marker, selected}) {
}
SMFace(int orient, int marker, bool selected) : SMFace(Orientable{orient}, SMFaceProperties{marker, selected}) {
}
SMFace(Orientable orient, SMFaceProperties prop)
: Orientable(orient), SMFaceProperties(prop)
{
}
friend std::ostream& operator<<(std::ostream& output, const SMFace& f)
{
output << "SMFace("
<< "m:" << f.marker
<< ";sel:" << std::boolalpha << f.selected
<< ";o:" << f.orientation << ")";
return output;
}
std::string to_string() const
{
std::ostringstream output;
output << *this;
return output.str();
}
};
namespace surfmesh_detail
{
struct surfmesh_traits
{
using KeyType = int;
using NodeTypes = util::type_holder<SMGlobal, SMVertex, SMEdge, SMFace>;
using EdgeTypes = util::type_holder<casc::Orientable, casc::Orientable, casc::Orientable>;
};
} // end namespace surfmesh_detail
using SurfaceMesh = casc::simplicial_complex<surfmesh_detail::surfmesh_traits>;
tensor<double, 3, 2> getTangent(const SurfaceMesh& mesh, SurfaceMesh::SimplexID<1> vertexID);
tensor<double, 3, 2> getTangent(const SurfaceMesh& mesh, SurfaceMesh::SimplexID<3> faceID);
Vector getNormalFromTangent(const tensor<double, 3, 2> tangent);
Vector getNormal(const SurfaceMesh& mesh, SurfaceMesh::SimplexID<1> vertexID);
Vector getNormal(const SurfaceMesh& mesh, SurfaceMesh::SimplexID<3> faceID);
namespace surfacemesh_detail
{
void decimateVertex(SurfaceMesh& mesh, SurfaceMesh::SimplexID<1> vertexID, std::size_t rings = 2);
tensor<double, 3, 2> computeLocalStructureTensor(
const SurfaceMesh& mesh,
const SurfaceMesh::SimplexID<1> vertexID,
const int rings);
tensor<double, 3, 2> computeLSTFromCache(
const SurfaceMesh& mesh,
const SurfaceMesh::SimplexID<1> vertexID,
const int rings);
template <std::size_t dimension>
auto getTangentH(const SurfaceMesh& mesh,
const tensor<double, dimension, 1>& origin,
SurfaceMesh::SimplexID<SurfaceMesh::topLevel> curr)
{
return (*curr).orientation;
}
template <std::size_t level, std::size_t dimension>
auto getTangentH(const SurfaceMesh& mesh,
const tensor<double, dimension, 1>& origin,
SurfaceMesh::SimplexID<level> curr)
{
tensor<double, dimension, SurfaceMesh::topLevel - level> rval;
auto cover = mesh.get_cover(curr);
for (auto alpha : cover)
{
auto edge = *mesh.get_edge_up(curr, alpha);
const auto& v = (*mesh.get_simplex_up({alpha})).position;
auto next = mesh.get_simplex_up(curr, alpha);
rval += edge.orientation * (v-origin) * getTangentH(mesh, origin, next);
}
return rval/cover.size();
}
template <std::size_t dimension>
auto getTangentF(const SurfaceMesh& mesh,
const tensor<double, dimension, 1>& origin,
SurfaceMesh::SimplexID<SurfaceMesh::topLevel> curr,
std::set<SurfaceMesh::KeyType>& cover)
{
return (*curr).orientation;
}
template <std::size_t level, std::size_t dimension>
auto getTangentF(const SurfaceMesh& mesh,
const tensor<double, dimension, 1>& origin,
SurfaceMesh::SimplexID<level> curr,
std::set<SurfaceMesh::KeyType>& cover)
{
tensor<double, dimension, SurfaceMesh::topLevel - level> rval;
for (auto alpha : cover)
{
auto edge = *mesh.get_edge_up(curr, alpha);
const auto& v = (*mesh.get_simplex_up({alpha})).position;
auto next = mesh.get_simplex_up(curr, alpha);
auto coverup = cover;
coverup.erase(alpha);
rval += edge.orientation * (v-origin) * getTangentF(mesh, origin, next, coverup);
}
return rval/cover.size();
}
template <class K>
struct initLocalOrientation {};
template <std::size_t k>
struct initLocalOrientation<std::integral_constant<std::size_t, k> >
{
template <typename Iterator>
static void apply(SurfaceMesh& mesh,
const std::set<int>&& names,
Iterator begin,
Iterator end)
{
std::vector<SurfaceMesh::SimplexID<k+1> > next;
for (auto curr = begin; curr != end; ++curr)
{
auto currSimplexID = *curr;
for (auto a : mesh.get_cover(currSimplexID))
{
// Look for key a in names
auto find = names.find(a);
if (find != names.end())
{
next.push_back(mesh.get_simplex_up(currSimplexID, a));
int orient = 1;
for (auto b : mesh.get_name(currSimplexID))
{
if (a > b)
{
if (a > b)
{
orient *= -1;
}
else
{
break;
}
}
}
(*mesh.get_edge_up(currSimplexID, a)).orientation = orient;
}
}
}
initLocalOrientation<std::integral_constant<std::size_t, k+1> >::apply(mesh, std::move(names), next.begin(), next.end());
}
};
template <>
struct initLocalOrientation<std::integral_constant<std::size_t, SurfaceMesh::topLevel> > {
template <typename Iterator>
static void apply(SurfaceMesh& mesh, const std::set<int>&& names, Iterator begin, Iterator end){
}
};
bool computeLocalOrientation(SurfaceMesh& mesh, const std::vector<SurfaceMesh::SimplexID<2> >& edgeList);
void weightedVertexSmooth(SurfaceMesh& mesh, SurfaceMesh::SimplexID<1> vertexID, int rings);
Vector weightedVertexSmoothCache(SurfaceMesh& mesh,
SurfaceMesh::SimplexID<1> vertexID,
std::size_t rings);
void barycenterVertexSmooth(SurfaceMesh& mesh, SurfaceMesh::SimplexID<1> vertexID);
void edgeFlip(SurfaceMesh& mesh, SurfaceMesh::SimplexID<2> edgeID);
void edgeFlipCache(SurfaceMesh& mesh, SurfaceMesh::SimplexID<2> edgeID);
template <class Inserter>
void selectFlipEdges(const SurfaceMesh& mesh,
bool preserveRidges,
std::function<bool(const SurfaceMesh&, const SurfaceMesh::SimplexID<2> &)>&& checkFlip,
Inserter iter)
{
casc::NodeSet<SurfaceMesh::SimplexID<2> > ignoredEdges;
for (auto edgeID : mesh.get_level_id<2>())
{
if ((*edgeID).selected == true)
{
if (!ignoredEdges.count(edgeID))
{
auto up = mesh.get_cover(edgeID);
// The mesh is not a surface mesh...
if (up.size() > 2)
{
// std::cerr << "This edge participates in more than 2
// faces. "
// << "Returning..." << std::endl;
throw std::runtime_error("SurfaceMesh is not pseudomanifold. Found an edge connected to more than 2 faces.");
}
else if (up.size() < 2) // Edge is a boundary
{
// std::cerr << "This edge participates in fewer than 2
// faces. "
// << "Returning..." << std::endl;
continue;
}
// Check if the edge is a part of a tetrahedron.
if (mesh.exists<2>({up[0], up[1]}))
{
// std::cerr << "Found a tetrahedron cannot edge flip."
// << std::endl;
continue;
}
// Check if we're on a ridge. This prevents folding also.
if (preserveRidges)
{
auto a = getNormal(mesh, mesh.get_simplex_up(edgeID, up[0]));
auto b = getNormal(mesh, mesh.get_simplex_up(edgeID, up[1]));
auto val = angle(a, b);
if (val > 60)
{
continue;
}
}
// Check the flip using user function
if (checkFlip(mesh, edgeID))
{
*iter++ = edgeID; // Insert into edges to flip
// The local topology will be changed by edge flip.
// Don't flip edges which share a common face.
std::set<SurfaceMesh::SimplexID<2> > tmpIgnored;
kneighbors(mesh, edgeID, 3, tmpIgnored);
ignoredEdges.insert(tmpIgnored.begin(), tmpIgnored.end());
// Local neighborhood append. Larger neighborhood selected
// above appears to work better...
// neighbors(mesh, edgeID, std::inserter(ignoredEdges,
// ignoredEdges.end()));
}
}
}
}
}
bool checkEdgeFlip(const SurfaceMesh& mesh,
bool preserveRidges,
SurfaceMesh::SimplexID<2> edgeID,
std::function<bool(const SurfaceMesh&, const SurfaceMesh::SimplexID<2> &)>&& checkFlip
);
bool checkFlipAngle(const SurfaceMesh& mesh, const SurfaceMesh::SimplexID<2>& edgeID);
int checkFlipValenceExcess(const SurfaceMesh& mesh, const SurfaceMesh::SimplexID<2>& edgeID);
void normalSmoothH(SurfaceMesh& mesh, const SurfaceMesh::SimplexID<1> vertexID, const double k);
void findHoles(const SurfaceMesh& mesh,
std::vector<std::vector<SurfaceMesh::SimplexID<2> > >& holeList);
bool orderBoundaryEdgeRing(const SurfaceMesh& mesh,
std::set<SurfaceMesh::SimplexID<2> >& unvisitedBdryEdges,
std::vector<SurfaceMesh::SimplexID<1> >& visitedVerts,
std::vector<SurfaceMesh::SimplexID<2> >& bdryRing);
void edgeRingToVertices(const SurfaceMesh& mesh,
std::vector<SurfaceMesh::SimplexID<2> >& edgeRing,
std::back_insert_iterator<std::vector<SurfaceMesh::SimplexID<1> > > iter);
void triangulateHoleHelper(SurfaceMesh& mesh,
std::vector<SurfaceMesh::SimplexID<1> >& boundary,
const SMFace& fdata,
std::back_insert_iterator<std::vector<SurfaceMesh::SimplexID<2> > > iter);
void triangulateHole(SurfaceMesh& mesh,
std::vector<SurfaceMesh::SimplexID<1> >& sortedVerts,
const SMFace& fdata,
std::vector<SurfaceMesh::SimplexID<2> >& edgeList);
template <typename Complex>
struct CopyHelper
{
using SimplexSet = typename casc::SimplexSet<Complex>;
using KeyType = typename Complex::KeyType;
template <std::size_t k>
static void apply(Complex& before,
Complex& after,
const SimplexSet& S)
{
for (auto sID : casc::get<k>(S))
{
auto name = before.get_name(sID);
auto data = *sID;
after.insert(name, data);
}
}
};
template <typename Iterator>
void vertexGrabber(const SurfaceMesh& F,
int need,
std::vector<SurfaceMesh::SimplexID<1>>& nbors,
Iterator begin,
Iterator end)
{
if (need <= 0) return;
std::set<SurfaceMesh::SimplexID<1>> next;
for (; begin != end; ++begin)
{
for (auto a : F.get_cover(*begin))
{
auto id = F.get_simplex_up(*begin, a);
for (auto b : F.get_name(id))
{
if (b != a){
auto nbor = F.get_simplex_down(id, b);
if (std::find(nbors.begin(), nbors.end(), nbor) == nbors.end())
{
// Haven't visited so push into next ring and
nbors.push_back(nbor);
next.insert(nbor);
if (--need == 0) return;
}
}
}
}
}
vertexGrabber(F, need, nbors, next.begin(), next.end());
}
template <class Complex>
void vertexGrabber(const Complex& F,
int need,
std::vector<typename Complex::template SimplexID<1>>& nbors,
typename Complex::template SimplexID<1> vid)
{
if (need <= 0) return;
std::set<SurfaceMesh::SimplexID<1>> next;
for (auto a : F.get_cover(vid))
{
auto id = F.get_simplex_up(vid, a);
for (auto b : F.get_name(id))
{
if (b != a){
auto nbor = F.get_simplex_down(id, b);
if (std::find(nbors.begin(), nbors.end(), nbor) == nbors.end())
{
// Haven't visited so push into next ring and
nbors.push_back(nbor);
next.insert(nbor);
if (--need == 0) return;
}
}
}
}
vertexGrabber(F, need, nbors, next.begin(), next.end());
}
} // end namespace surfacemesh_detail
std::unique_ptr<SurfaceMesh> readOFF(const std::string& filename);
void writeOFF(const std::string& filename, const SurfaceMesh& mesh);
std::unique_ptr<SurfaceMesh> readOBJ(const std::string& filename);
void writeOBJ(const std::string& filename, const SurfaceMesh& mesh);
void print(const SurfaceMesh& mesh);
void printQualityInfo(const std::string& filename, const SurfaceMesh& mesh);
void generateHistogram(const SurfaceMesh& mesh);
std::tuple<double, double, int, int> getMinMaxAngles(const SurfaceMesh& mesh,
double maxMinAngle, double minMaxAngle);
double getArea(const SurfaceMesh& mesh);
double getArea(const SurfaceMesh& mesh, SurfaceMesh::SimplexID<3> faceID);
double getArea(Vertex a, Vertex b, Vertex c);
REAL getArea(std::array<Vertex, 3> t);
double getVolume(const SurfaceMesh& mesh);
bool hasHole(const SurfaceMesh& mesh);
inline std::size_t getValence(const SurfaceMesh& mesh, const SurfaceMesh::SimplexID<1> vertexID){
return mesh.get_cover(vertexID).size();
}
std::tuple<REAL*, REAL*, REAL*, REAL*, std::map<typename SurfaceMesh::KeyType, typename SurfaceMesh::KeyType> >
curvatureViaMDSB(const SurfaceMesh& mesh);
std::tuple<REAL*, REAL*, REAL*, REAL*, std::map<typename SurfaceMesh::KeyType, typename SurfaceMesh::KeyType> >
curvatureViaJets(const SurfaceMesh& mesh, std::size_t dJet = 2, std::size_t dPrime = 2);
// void osculatingJets(const SurfaceMesh&mesh, std::size_t dJet = 2, std::size_t
// dPrime = 2);
//
// @param mesh The mesh
// @param[in] v Displacement vector
//
void translate(SurfaceMesh& mesh, Vector v);
void translate(SurfaceMesh& mesh, double dx, double dy, double dz);
void scale(SurfaceMesh& mesh, Vector v);
void scale(SurfaceMesh& mesh, double sx, double sy, double sz);
void scale(SurfaceMesh& mesh, double s);
std::pair<Vector, double> getCenterRadius(SurfaceMesh& mesh);
void centeralize(SurfaceMesh& mesh);
void smoothMesh(SurfaceMesh& mesh, int maxIter, bool preserveRidges, std::size_t rings = 2, bool verbose = false);
void coarse(SurfaceMesh& mesh, double coarseRate, double flatRate, double denseWeight, std::size_t rings = 2, bool verbose = false);
void coarse_dense(SurfaceMesh& mesh, REAL threshold, REAL weight, std::size_t rings = 2, bool verbose = false);
void coarse_flat(SurfaceMesh& mesh, REAL threshold, REAL weight, std::size_t rings = 2, bool verbose = false);
void normalSmooth(SurfaceMesh& mesh);
void fillHoles(SurfaceMesh& mesh);
void flipNormals(SurfaceMesh& mesh);
std::unique_ptr<SurfaceMesh> refineMesh(const SurfaceMesh& mesh);
std::unique_ptr<SurfaceMesh> sphere(int order);
std::unique_ptr<SurfaceMesh> cube(int order);
std::vector<std::unique_ptr<SurfaceMesh> > splitSurfaces(SurfaceMesh& mesh);
void cacheNormals(SurfaceMesh& mesh);
std::tuple<bool, int, int, int> getBettiNumbers(SurfaceMesh& mesh);
} // end namespace gamer