kindred/create
1
1
Fork 0
Code Issues 37 Pull Requests Actions Packages Projects 9 Releases 2 Wiki Activity
Files
ba2c2ca5ad17fff53fee48826b0f2f669e6524fc
create/src/Mod/Mesh/App/Core/Triangulation.cpp

863 lines
29 KiB
C++
Raw Blame History

/***************************************************************************
* Copyright (c) 2005 Werner Mayer <wmayer[at]users.sourceforge.net> *
* *
* This file is part of the FreeCAD CAx development system. *
* *
* This library is free software; you can redistribute it and/or *
* modify it under the terms of the GNU Library General Public *
* License as published by the Free Software Foundation; either *
* version 2 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 Library General Public License for more details. *
* *
* You should have received a copy of the GNU Library General Public *
* License along with this library; see the file COPYING.LIB. If not, *
* write to the Free Software Foundation, Inc., 59 Temple Place, *
* Suite 330, Boston, MA 02111-1307, USA *
* *
***************************************************************************/
#include "PreCompiled.h"
#ifndef _PreComp_
#include <cmath>
#include <limits>
#include <queue>
#endif
#include <Base/Console.h>
#include <Base/Exception.h>
#include <Mod/Mesh/App/WildMagic4/Wm4Delaunay2.h>
#include "Approximation.h"
#include "MeshKernel.h"
#include "Triangulation.h"
using namespace MeshCore;
bool TriangulationVerifier::Accept(const Base::Vector3f& n,
const Base::Vector3f& p1,
const Base::Vector3f& p2,
const Base::Vector3f& p3) const
{
float ref_dist = (p2 - p1) * n;
float tri_dist = (p3 - p1) * n;
return (ref_dist * tri_dist <= 0.0F);
}
bool TriangulationVerifier::MustFlip(const Base::Vector3f& n1, const Base::Vector3f& n2) const
{
return n1.Dot(n2) <= 0.0F;
}
bool TriangulationVerifierV2::Accept(const Base::Vector3f& n,
const Base::Vector3f& p1,
const Base::Vector3f& p2,
const Base::Vector3f& p3) const
{
float ref_dist = (p2 - p1) * n;
float tri_dist = (p3 - p1) * n;
float prod = ref_dist * tri_dist;
(void)prod;
return true;
}
bool TriangulationVerifierV2::MustFlip(const Base::Vector3f& n1, const Base::Vector3f& n2) const
{
float dot = n1.Dot(n2);
(void)dot;
return false;
}
// ----------------------------------------------------------------------------
AbstractPolygonTriangulator::AbstractPolygonTriangulator()
: _discard {false}
, _verifier {new TriangulationVerifier()}
{}
AbstractPolygonTriangulator::~AbstractPolygonTriangulator()
{
delete _verifier;
}
TriangulationVerifier* AbstractPolygonTriangulator::GetVerifier() const
{
return _verifier;
}
void AbstractPolygonTriangulator::SetVerifier(TriangulationVerifier* v)
{
delete _verifier;
_verifier = v;
}
void AbstractPolygonTriangulator::SetPolygon(const std::vector<Base::Vector3f>& raclPoints)
{
this->_points = raclPoints;
if (!this->_points.empty()) {
if (this->_points.front() == this->_points.back()) {
this->_points.pop_back();
}
}
}
std::vector<Base::Vector3f> AbstractPolygonTriangulator::GetPolygon() const
{
return _points;
}
float AbstractPolygonTriangulator::GetLength() const
{
float len = 0.0F;
if (_points.size() > 2) {
for (auto it = _points.begin(); it != _points.end(); ++it) {
std::vector<Base::Vector3f>::const_iterator jt = it + 1;
if (jt == _points.end()) {
jt = _points.begin();
}
len += Base::Distance(*it, *jt);
}
}
return len;
}
std::vector<Base::Vector3f> AbstractPolygonTriangulator::AddedPoints() const
{
// Apply the inverse transformation to project back to world coordinates
std::vector<Base::Vector3f> added;
added.reserve(_newpoints.size());
for (auto point : _newpoints) {
added.push_back(_inverse * point);
}
return added;
}
Base::Matrix4D AbstractPolygonTriangulator::GetTransformToFitPlane() const
{
PlaneFit planeFit;
for (auto point : _points) {
planeFit.AddPoint(point);
}
if (planeFit.Fit() >= std::numeric_limits<float>::max()) {
throw Base::RuntimeError("Plane fit failed");
}
Base::Vector3f bs = planeFit.GetBase();
Base::Vector3f ex = planeFit.GetDirU();
Base::Vector3f ey = planeFit.GetDirV();
Base::Vector3f ez = planeFit.GetNormal();
// build the matrix for the inverse transformation
Base::Matrix4D rInverse;
rInverse.setToUnity();
rInverse[0][0] = static_cast<double>(ex.x);
rInverse[0][1] = static_cast<double>(ey.x);
rInverse[0][2] = static_cast<double>(ez.x);
rInverse[0][3] = static_cast<double>(bs.x);
rInverse[1][0] = static_cast<double>(ex.y);
rInverse[1][1] = static_cast<double>(ey.y);
rInverse[1][2] = static_cast<double>(ez.y);
rInverse[1][3] = static_cast<double>(bs.y);
rInverse[2][0] = static_cast<double>(ex.z);
rInverse[2][1] = static_cast<double>(ey.z);
rInverse[2][2] = static_cast<double>(ez.z);
rInverse[2][3] = static_cast<double>(bs.z);
return rInverse;
}
std::vector<Base::Vector3f> AbstractPolygonTriangulator::ProjectToFitPlane()
{
std::vector<Base::Vector3f> proj = _points;
_inverse = GetTransformToFitPlane();
Base::Vector3f bs(static_cast<float>(_inverse[0][3]),
static_cast<float>(_inverse[1][3]),
static_cast<float>(_inverse[2][3]));
Base::Vector3f ex(static_cast<float>(_inverse[0][0]),
static_cast<float>(_inverse[1][0]),
static_cast<float>(_inverse[2][0]));
Base::Vector3f ey(static_cast<float>(_inverse[0][1]),
static_cast<float>(_inverse[1][1]),
static_cast<float>(_inverse[2][1]));
for (auto& jt : proj) {
jt.TransformToCoordinateSystem(bs, ex, ey);
}
return proj;
}
void AbstractPolygonTriangulator::PostProcessing(const std::vector<Base::Vector3f>& points)
{
// For a good approximation we should have enough points, i.e. for 9 parameters
// for the fit function we should have at least 50 points.
unsigned int uMinPts = 50;
PolynomialFit polyFit;
Base::Vector3f bs(static_cast<float>(_inverse[0][3]),
static_cast<float>(_inverse[1][3]),
static_cast<float>(_inverse[2][3]));
Base::Vector3f ex(static_cast<float>(_inverse[0][0]),
static_cast<float>(_inverse[1][0]),
static_cast<float>(_inverse[2][0]));
Base::Vector3f ey(static_cast<float>(_inverse[0][1]),
static_cast<float>(_inverse[1][1]),
static_cast<float>(_inverse[2][1]));
for (auto pt : points) {
pt.TransformToCoordinateSystem(bs, ex, ey);
polyFit.AddPoint(pt);
}
if (polyFit.CountPoints() >= uMinPts && polyFit.Fit() < std::numeric_limits<float>::max()) {
for (auto& newpoint : _newpoints) {
newpoint.z = static_cast<float>(polyFit.Value(newpoint.x, newpoint.y));
}
}
}
MeshGeomFacet AbstractPolygonTriangulator::GetTriangle(const MeshPointArray& points,
const MeshFacet& facet) const
{
MeshGeomFacet triangle;
triangle._aclPoints[0] = points[facet._aulPoints[0]];
triangle._aclPoints[1] = points[facet._aulPoints[1]];
triangle._aclPoints[2] = points[facet._aulPoints[2]];
return triangle;
}
bool AbstractPolygonTriangulator::TriangulatePolygon()
{
try {
if (!this->_indices.empty() && this->_points.size() != this->_indices.size()) {
Base::Console().log("Triangulation: %d points <> %d indices\n",
_points.size(),
_indices.size());
return false;
}
bool ok = Triangulate();
if (ok) {
Done();
}
return ok;
}
catch (const Base::Exception& e) {
Base::Console().log("Triangulation: %s\n", e.what());
return false;
}
catch (const std::exception& e) {
Base::Console().log("Triangulation: %s\n", e.what());
return false;
}
catch (...) {
return false;
}
}
std::vector<PointIndex> AbstractPolygonTriangulator::GetInfo() const
{
return _info;
}
void AbstractPolygonTriangulator::Discard()
{
if (!_discard) {
_discard = true;
_info.pop_back();
}
}
void AbstractPolygonTriangulator::Reset()
{}
void AbstractPolygonTriangulator::Done()
{
_info.push_back(_points.size());
_discard = false;
}
// -------------------------------------------------------------
EarClippingTriangulator::EarClippingTriangulator() = default;
bool EarClippingTriangulator::Triangulate()
{
_facets.clear();
_triangles.clear();
std::vector<Base::Vector3f> pts = ProjectToFitPlane();
std::vector<PointIndex> result;
// Invoke the triangulator to triangulate this polygon.
Triangulate::Process(pts, result);
// print out the results.
size_t tcount = result.size() / 3;
bool ok = tcount + 2 == _points.size();
if (tcount > _points.size()) {
return false; // no valid triangulation
}
MeshGeomFacet clFacet;
MeshFacet clTopFacet;
for (size_t i = 0; i < tcount; i++) {
if (Triangulate::_invert) {
clFacet._aclPoints[0] = _points[result[i * 3 + 0]];
clFacet._aclPoints[2] = _points[result[i * 3 + 1]];
clFacet._aclPoints[1] = _points[result[i * 3 + 2]];
clTopFacet._aulPoints[0] = result[i * 3 + 0];
clTopFacet._aulPoints[2] = result[i * 3 + 1];
clTopFacet._aulPoints[1] = result[i * 3 + 2];
}
else {
clFacet._aclPoints[0] = _points[result[i * 3 + 0]];
clFacet._aclPoints[1] = _points[result[i * 3 + 1]];
clFacet._aclPoints[2] = _points[result[i * 3 + 2]];
clTopFacet._aulPoints[0] = result[i * 3 + 0];
clTopFacet._aulPoints[1] = result[i * 3 + 1];
clTopFacet._aulPoints[2] = result[i * 3 + 2];
}
_triangles.push_back(clFacet);
_facets.push_back(clTopFacet);
}
return ok;
}
float EarClippingTriangulator::Triangulate::Area(const std::vector<Base::Vector3f>& contour)
{
int n = contour.size();
float A = 0.0F;
for (int p = n - 1, q = 0; q < n; p = q++) {
A += contour[p].x * contour[q].y - contour[q].x * contour[p].y;
}
return A * 0.5F;
}
/*
InsideTriangle decides if a point P is Inside of the triangle
defined by A, B, C.
*/
bool EarClippingTriangulator::Triangulate::InsideTriangle(float Ax,
float Ay,
float Bx,
float By,
float Cx,
float Cy,
float Px,
float Py)
{
float ax {}, ay {}, bx {}, by {}, cx {}, cy {}, apx {}, apy {}, bpx {}, bpy {}, cpx {}, cpy {};
float cCROSSap {}, bCROSScp {}, aCROSSbp {};
ax = Cx - Bx;
ay = Cy - By;
bx = Ax - Cx;
by = Ay - Cy;
cx = Bx - Ax;
cy = By - Ay;
apx = Px - Ax;
apy = Py - Ay;
bpx = Px - Bx;
bpy = Py - By;
cpx = Px - Cx;
cpy = Py - Cy;
aCROSSbp = ax * bpy - ay * bpx;
cCROSSap = cx * apy - cy * apx;
bCROSScp = bx * cpy - by * cpx;
return ((aCROSSbp >= std::numeric_limits<float>::epsilon())
&& (bCROSScp >= std::numeric_limits<float>::epsilon())
&& (cCROSSap >= std::numeric_limits<float>::epsilon()));
}
bool EarClippingTriangulator::Triangulate::Snip(const std::vector<Base::Vector3f>& contour,
int u,
int v,
int w,
int n,
int* V)
{
int p {};
float Ax {}, Ay {}, Bx {}, By {}, Cx {}, Cy {}, Px {}, Py {};
Ax = contour[V[u]].x;
Ay = contour[V[u]].y;
Bx = contour[V[v]].x;
By = contour[V[v]].y;
Cx = contour[V[w]].x;
Cy = contour[V[w]].y;
constexpr float eps = std::numeric_limits<float>::epsilon();
if (eps > (((Bx - Ax) * (Cy - Ay)) - ((By - Ay) * (Cx - Ax)))) {
return false;
}
for (p = 0; p < n; p++) {
if ((p == u) || (p == v) || (p == w)) {
continue;
}
Px = contour[V[p]].x;
Py = contour[V[p]].y;
if (InsideTriangle(Ax, Ay, Bx, By, Cx, Cy, Px, Py)) {
return false;
}
}
return true;
}
bool EarClippingTriangulator::Triangulate::_invert = false;
bool EarClippingTriangulator::Triangulate::Process(const std::vector<Base::Vector3f>& contour,
std::vector<PointIndex>& result)
{
/* allocate and initialize list of Vertices in polygon */
int n = contour.size();
if (n < 3) {
return false;
}
int* V = new int[n];
/* we want a counter-clockwise polygon in V */
if (0.0F < Area(contour)) {
for (int v = 0; v < n; v++) {
V[v] = v;
}
_invert = true;
}
// for(int v=0; v<n; v++) V[v] = (n-1)-v;
else {
for (int v = 0; v < n; v++) {
V[v] = (n - 1) - v;
}
_invert = false;
}
int nv = n;
/* remove nv-2 Vertices, creating 1 triangle every time */
int count = 2 * nv; /* error detection */
for (int v = nv - 1; nv > 2;) {
/* if we loop, it is probably a non-simple polygon */
if (0 >= (count--)) {
//** Triangulate: ERROR - probable bad polygon!
delete[] V;
return false;
}
/* three consecutive vertices in current polygon, <u,v,w> */
int u = v;
if (nv <= u) {
u = 0; /* previous */
}
v = u + 1;
if (nv <= v) {
v = 0; /* new v */
}
int w = v + 1;
if (nv <= w) {
w = 0; /* next */
}
if (Snip(contour, u, v, w, nv, V)) {
int a {}, b {}, c {}, s {}, t {};
/* true names of the vertices */
a = V[u];
b = V[v];
c = V[w];
/* output Triangle */
result.push_back(a);
result.push_back(b);
result.push_back(c);
/* remove v from remaining polygon */
for (s = v, t = v + 1; t < nv; s++, t++) {
V[s] = V[t];
}
nv--;
/* reset error detection counter */
count = 2 * nv;
}
}
delete[] V;
return true;
}
// -------------------------------------------------------------
QuasiDelaunayTriangulator::QuasiDelaunayTriangulator() = default;
bool QuasiDelaunayTriangulator::Triangulate()
{
if (!EarClippingTriangulator::Triangulate()) {
return false; // no valid triangulation
}
// For each internal edge get the adjacent facets. When doing an edge swap we must update
// this structure.
std::map<std::pair<PointIndex, PointIndex>, std::vector<FacetIndex>> aEdge2Face;
for (auto pI = _facets.begin(); pI != _facets.end(); ++pI) {
for (int i = 0; i < 3; i++) {
PointIndex ulPt0 = std::min<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
PointIndex ulPt1 = std::max<PointIndex>(pI->_aulPoints[i], pI->_aulPoints[(i + 1) % 3]);
// ignore borderlines of the polygon
if ((ulPt1 - ulPt0) % (_points.size() - 1) > 1) {
aEdge2Face[std::pair<PointIndex, PointIndex>(ulPt0, ulPt1)].push_back(
pI - _facets.begin());
}
}
}
// fill up this list with all internal edges and perform swap edges until this list is empty
std::list<std::pair<PointIndex, PointIndex>> aEdgeList;
std::map<std::pair<PointIndex, PointIndex>, std::vector<FacetIndex>>::iterator pE;
for (pE = aEdge2Face.begin(); pE != aEdge2Face.end(); ++pE) {
aEdgeList.push_back(pE->first);
}
// to be sure to avoid an endless loop
size_t uMaxIter = 5 * aEdge2Face.size();
// Perform a swap edge where needed
while (!aEdgeList.empty() && uMaxIter > 0) {
// get the first edge and remove it from the list
std::pair<PointIndex, PointIndex> aEdge = aEdgeList.front();
aEdgeList.pop_front();
uMaxIter--;
// get the adjacent facets to this edge
pE = aEdge2Face.find(aEdge);
// this edge has been removed some iterations before
if (pE == aEdge2Face.end()) {
continue;
}
MeshFacet& rF1 = _facets[pE->second[0]];
MeshFacet& rF2 = _facets[pE->second[1]];
unsigned short side1 = rF1.Side(aEdge.first, aEdge.second);
Base::Vector3f cP1 = _points[rF1._aulPoints[side1]];
Base::Vector3f cP2 = _points[rF1._aulPoints[(side1 + 1) % 3]];
Base::Vector3f cP3 = _points[rF1._aulPoints[(side1 + 2) % 3]];
unsigned short side2 = rF2.Side(aEdge.first, aEdge.second);
Base::Vector3f cP4 = _points[rF2._aulPoints[(side2 + 2) % 3]];
MeshGeomFacet cT1(cP1, cP2, cP3);
float fMax1 = cT1.MaximumAngle();
MeshGeomFacet cT2(cP2, cP1, cP4);
float fMax2 = cT2.MaximumAngle();
MeshGeomFacet cT3(cP4, cP3, cP1);
float fMax3 = cT3.MaximumAngle();
MeshGeomFacet cT4(cP3, cP4, cP2);
float fMax4 = cT4.MaximumAngle();
float fMax12 = std::max<float>(fMax1, fMax2);
float fMax34 = std::max<float>(fMax3, fMax4);
// We must make sure that the two adjacent triangles builds a convex polygon, otherwise
// the swap edge operation is illegal
Base::Vector3f cU = cP2 - cP1;
Base::Vector3f cV = cP4 - cP3;
// build a helper plane through cP1 that must separate cP3 and cP4
Base::Vector3f cN1 = (cU % cV) % cU;
if (((cP3 - cP1) * cN1) * ((cP4 - cP1) * cN1) >= 0.0F) {
continue; // not convex
}
// build a helper plane through cP3 that must separate cP1 and cP2
Base::Vector3f cN2 = (cU % cV) % cV;
if (((cP1 - cP3) * cN2) * ((cP2 - cP3) * cN2) >= 0.0F) {
continue; // not convex
}
// ok, here we should perform a swap edge to minimize the maximum angle
if (fMax12 > fMax34) {
rF1._aulPoints[(side1 + 1) % 3] = rF2._aulPoints[(side2 + 2) % 3];
rF2._aulPoints[(side2 + 1) % 3] = rF1._aulPoints[(side1 + 2) % 3];
// adjust the edge list
for (int i = 0; i < 3; i++) {
std::map<std::pair<PointIndex, PointIndex>, std::vector<FacetIndex>>::iterator it;
// first facet
PointIndex ulPt0 =
std::min<PointIndex>(rF1._aulPoints[i], rF1._aulPoints[(i + 1) % 3]);
PointIndex ulPt1 =
std::max<PointIndex>(rF1._aulPoints[i], rF1._aulPoints[(i + 1) % 3]);
it = aEdge2Face.find(std::make_pair(ulPt0, ulPt1));
if (it != aEdge2Face.end()) {
if (it->second[0] == pE->second[1]) {
it->second[0] = pE->second[0];
}
else if (it->second[1] == pE->second[1]) {
it->second[1] = pE->second[0];
}
aEdgeList.push_back(it->first);
}
// second facet
ulPt0 = std::min<PointIndex>(rF2._aulPoints[i], rF2._aulPoints[(i + 1) % 3]);
ulPt1 = std::max<PointIndex>(rF2._aulPoints[i], rF2._aulPoints[(i + 1) % 3]);
it = aEdge2Face.find(std::make_pair(ulPt0, ulPt1));
if (it != aEdge2Face.end()) {
if (it->second[0] == pE->second[0]) {
it->second[0] = pE->second[1];
}
else if (it->second[1] == pE->second[0]) {
it->second[1] = pE->second[1];
}
aEdgeList.push_back(it->first);
}
}
// Now we must remove the edge and replace it through the new edge
PointIndex ulPt0 = std::min<PointIndex>(rF1._aulPoints[(side1 + 1) % 3],
rF2._aulPoints[(side2 + 1) % 3]);
PointIndex ulPt1 = std::max<PointIndex>(rF1._aulPoints[(side1 + 1) % 3],
rF2._aulPoints[(side2 + 1) % 3]);
std::pair<PointIndex, PointIndex> aNewEdge = std::make_pair(ulPt0, ulPt1);
aEdge2Face[aNewEdge] = pE->second;
aEdge2Face.erase(pE);
}
}
return true;
}
// -------------------------------------------------------------
namespace MeshCore
{
namespace Triangulation
{
struct Vertex2d_Less
{
bool operator()(const Base::Vector3f& p, const Base::Vector3f& q) const
{
if (std::fabs(p.x - q.x) < MeshDefinitions::_fMinPointDistanceD1) {
if (std::fabs(p.y - q.y) < MeshDefinitions::_fMinPointDistanceD1) {
return false;
}
return p.y < q.y;
}
return p.x < q.x;
}
};
struct Vertex2d_EqualTo
{
bool operator()(const Base::Vector3f& p, const Base::Vector3f& q) const
{
if (std::fabs(p.x - q.x) < MeshDefinitions::_fMinPointDistanceD1
&& std::fabs(p.y - q.y) < MeshDefinitions::_fMinPointDistanceD1) {
return true;
}
return false;
}
};
} // namespace Triangulation
} // namespace MeshCore
DelaunayTriangulator::DelaunayTriangulator() = default;
bool DelaunayTriangulator::Triangulate()
{
// before starting the triangulation we must make sure that all polygon
// points are different
std::vector<Base::Vector3f> aPoints = _points;
// sort the points ascending x,y coordinates
std::sort(aPoints.begin(), aPoints.end(), Triangulation::Vertex2d_Less());
// if there are two adjacent points whose distance is less then an epsilon
if (std::adjacent_find(aPoints.begin(), aPoints.end(), Triangulation::Vertex2d_EqualTo())
< aPoints.end()) {
return false;
}
_facets.clear();
_triangles.clear();
std::vector<Wm4::Vector2d> akVertex;
akVertex.reserve(_points.size());
for (const auto& point : _points) {
akVertex.emplace_back(static_cast<double>(point.x), static_cast<double>(point.y));
}
Wm4::Delaunay2d del(static_cast<int>(akVertex.size()),
akVertex.data(),
0.001,
false,
Wm4::Query::QT_INT64);
int iTQuantity = del.GetSimplexQuantity();
auto numFaces = static_cast<std::size_t>(iTQuantity);
std::vector<int> aiTVertex(3 * numFaces);
bool succeeded = false;
if (numFaces > 0) {
size_t uiSize = 3 * numFaces * sizeof(int);
Wm4::System::Memcpy(aiTVertex.data(), uiSize, del.GetIndices(), uiSize);
// If H is the number of hull edges and N is the number of vertices,
// then the triangulation must have 2*N-2-H triangles and 3*N-3-H
// edges.
int iEQuantity = 0;
int* aiIndex = nullptr;
del.GetHull(iEQuantity, aiIndex);
int iUniqueVQuantity = del.GetUniqueVertexQuantity();
int iTVerify = 2 * iUniqueVQuantity - 2 - iEQuantity;
(void)iTVerify; // avoid warning in release build
succeeded = (iTVerify == iTQuantity);
int iEVerify = 3 * iUniqueVQuantity - 3 - iEQuantity;
(void)iEVerify; // avoid warning about unused variable
delete[] aiIndex;
}
MeshGeomFacet triangle;
MeshFacet facet;
for (std::size_t i = 0; i < numFaces; i++) {
for (std::size_t j = 0; j < 3; j++) {
auto index = static_cast<size_t>(aiTVertex[3 * i + j]);
facet._aulPoints[j] = static_cast<PointIndex>(index);
triangle._aclPoints[j].x = static_cast<float>(akVertex[index].X());
triangle._aclPoints[j].y = static_cast<float>(akVertex[index].Y());
}
_triangles.push_back(triangle);
_facets.push_back(facet);
}
return succeeded;
}
// -------------------------------------------------------------
FlatTriangulator::FlatTriangulator() = default;
bool FlatTriangulator::Triangulate()
{
_newpoints.clear();
// before starting the triangulation we must make sure that all polygon
// points are different
std::vector<Base::Vector3f> aPoints = ProjectToFitPlane();
std::vector<Base::Vector3f> tmp = aPoints;
// sort the points ascending x,y coordinates
std::sort(tmp.begin(), tmp.end(), Triangulation::Vertex2d_Less());
// if there are two adjacent points whose distance is less then an epsilon
if (std::adjacent_find(tmp.begin(), tmp.end(), Triangulation::Vertex2d_EqualTo()) < tmp.end()) {
return false;
}
_facets.clear();
_triangles.clear();
// Todo: Implement algorithm for constraint delaunay triangulation
QuasiDelaunayTriangulator tria;
tria.SetPolygon(this->GetPolygon());
bool succeeded = tria.TriangulatePolygon();
this->_facets = tria.GetFacets();
this->_triangles = tria.GetTriangles();
return succeeded;
}
void FlatTriangulator::PostProcessing(const std::vector<Base::Vector3f>&)
{}
// -------------------------------------------------------------
ConstraintDelaunayTriangulator::ConstraintDelaunayTriangulator(float area)
: fMaxArea(area)
{
// silent warning: -Wunused-private-field
(void)fMaxArea;
}
bool ConstraintDelaunayTriangulator::Triangulate()
{
_newpoints.clear();
// before starting the triangulation we must make sure that all polygon
// points are different
std::vector<Base::Vector3f> aPoints = ProjectToFitPlane();
std::vector<Base::Vector3f> tmp = aPoints;
// sort the points ascending x,y coordinates
std::sort(tmp.begin(), tmp.end(), Triangulation::Vertex2d_Less());
// if there are two adjacent points whose distance is less then an epsilon
if (std::adjacent_find(tmp.begin(), tmp.end(), Triangulation::Vertex2d_EqualTo()) < tmp.end()) {
return false;
}
_facets.clear();
_triangles.clear();
// Todo: Implement algorithm for constraint delaunay triangulation
QuasiDelaunayTriangulator tria;
tria.SetPolygon(this->GetPolygon());
bool succeeded = tria.TriangulatePolygon();
this->_facets = tria.GetFacets();
this->_triangles = tria.GetTriangles();
return succeeded;
}
// -------------------------------------------------------------
#if 0
Triangulator::Triangulator(const MeshKernel& k, bool flat) : _kernel(k)
{
}
Triangulator::~Triangulator()
{
}
bool Triangulator::Triangulate()
{
return false;
}
MeshGeomFacet Triangulator::GetTriangle(const MeshPointArray&,
const MeshFacet& facet) const
{
return MeshGeomFacet();
}
void Triangulator::PostProcessing(const std::vector<Base::Vector3f>&)
{
}
void Triangulator::Discard()
{
AbstractPolygonTriangulator::Discard();
}
void Triangulator::Reset()
{
}
#endif
Reference in New Issue View Git Blame Copy Permalink