6b83807dc6
+ some triangulation stuff git-svn-id: https://free-cad.svn.sourceforge.net/svnroot/free-cad/trunk@5417 e8eeb9e2-ec13-0410-a4a9-efa5cf37419d
724 lines
25 KiB
C++
724 lines
25 KiB
C++
/***************************************************************************
|
|
* 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 <queue>
|
|
#endif
|
|
|
|
#include <Base/Console.h>
|
|
#include <Base/Exception.h>
|
|
#include "Triangulation.h"
|
|
#include "Approximation.h"
|
|
#include "Algorithm.h"
|
|
#include "MeshKernel.h"
|
|
|
|
#include <Mod/Mesh/App/WildMagic4/Wm4Delaunay2.h>
|
|
|
|
|
|
using namespace MeshCore;
|
|
|
|
|
|
AbstractPolygonTriangulator::AbstractPolygonTriangulator()
|
|
{
|
|
_discard = false;
|
|
}
|
|
|
|
AbstractPolygonTriangulator::~AbstractPolygonTriangulator()
|
|
{
|
|
}
|
|
|
|
void AbstractPolygonTriangulator::SetPolygon(const std::vector<Base::Vector3f>& raclPoints)
|
|
{
|
|
this->_points = raclPoints;
|
|
if (this->_points.size() > 0) {
|
|
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 (std::vector<Base::Vector3f>::const_iterator 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 (std::vector<Base::Vector3f>::const_iterator pt = _newpoints.begin(); pt != _newpoints.end(); ++pt)
|
|
added.push_back(_inverse * *pt);
|
|
return added;
|
|
}
|
|
|
|
Base::Matrix4D AbstractPolygonTriangulator::GetTransformToFitPlane() const
|
|
{
|
|
PlaneFit planeFit;
|
|
for (std::vector<Base::Vector3f>::const_iterator it = _points.begin(); it!=_points.end(); ++it)
|
|
planeFit.AddPoint(*it);
|
|
|
|
if (planeFit.Fit() == FLOAT_MAX)
|
|
throw Base::Exception("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] = ex.x; rInverse[0][1] = ey.x; rInverse[0][2] = ez.x; rInverse[0][3] = bs.x;
|
|
rInverse[1][0] = ex.y; rInverse[1][1] = ey.y; rInverse[1][2] = ez.y; rInverse[1][3] = bs.y;
|
|
rInverse[2][0] = ex.z; rInverse[2][1] = ey.z; rInverse[2][2] = ez.z; rInverse[2][3] = bs.z;
|
|
|
|
return rInverse;
|
|
}
|
|
|
|
std::vector<Base::Vector3f> AbstractPolygonTriangulator::ProjectToFitPlane()
|
|
{
|
|
std::vector<Base::Vector3f> proj = _points;
|
|
_inverse = GetTransformToFitPlane();
|
|
Base::Vector3f bs((float)_inverse[0][3], (float)_inverse[1][3], (float)_inverse[2][3]);
|
|
Base::Vector3f ex((float)_inverse[0][0], (float)_inverse[1][0], (float)_inverse[2][0]);
|
|
Base::Vector3f ey((float)_inverse[0][1], (float)_inverse[1][1], (float)_inverse[2][1]);
|
|
for (std::vector<Base::Vector3f>::iterator jt = proj.begin(); jt!=proj.end(); ++jt)
|
|
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((float)_inverse[0][3], (float)_inverse[1][3], (float)_inverse[2][3]);
|
|
Base::Vector3f ex((float)_inverse[0][0], (float)_inverse[1][0], (float)_inverse[2][0]);
|
|
Base::Vector3f ey((float)_inverse[0][1], (float)_inverse[1][1], (float)_inverse[2][1]);
|
|
|
|
for (std::vector<Base::Vector3f>::const_iterator it = points.begin(); it != points.end(); ++it) {
|
|
Base::Vector3f pt = *it;
|
|
pt.TransformToCoordinateSystem(bs, ex, ey);
|
|
polyFit.AddPoint(pt);
|
|
}
|
|
|
|
if (polyFit.CountPoints() >= uMinPts && polyFit.Fit() < FLOAT_MAX) {
|
|
for (std::vector<Base::Vector3f>::iterator pt = _newpoints.begin(); pt != _newpoints.end(); ++pt)
|
|
pt->z = (float)polyFit.Value(pt->x, pt->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->_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<unsigned long> 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()
|
|
{
|
|
}
|
|
|
|
EarClippingTriangulator::~EarClippingTriangulator()
|
|
{
|
|
}
|
|
|
|
bool EarClippingTriangulator::Triangulate()
|
|
{
|
|
_facets.clear();
|
|
_triangles.clear();
|
|
|
|
std::vector<Base::Vector3f> pts = ProjectToFitPlane();
|
|
std::vector<unsigned long> result;
|
|
|
|
// Invoke the triangulator to triangulate this polygon.
|
|
Triangulate::Process(pts,result);
|
|
|
|
// print out the results.
|
|
unsigned long tcount = result.size()/3;
|
|
|
|
bool ok = tcount+2 == _points.size();
|
|
if (tcount > _points.size())
|
|
return false; // no valid triangulation
|
|
|
|
MeshGeomFacet clFacet;
|
|
MeshFacet clTopFacet;
|
|
for (unsigned long 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 >= FLOAT_EPS) && (bCROSScp >= FLOAT_EPS) && (cCROSSap >= FLOAT_EPS));
|
|
}
|
|
|
|
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;
|
|
|
|
if (FLOAT_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<unsigned long> &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 m=0, v=nv-1; nv>2; ) {
|
|
/* if we loop, it is probably a non-simple polygon */
|
|
if (0 >= (count--)) {
|
|
//** Triangulate: ERROR - probable bad polygon!
|
|
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 );
|
|
|
|
m++;
|
|
|
|
/* remove v from remaining polygon */
|
|
for(s=v,t=v+1;t<nv;s++,t++) V[s] = V[t]; nv--;
|
|
|
|
/* resest error detection counter */
|
|
count = 2*nv;
|
|
}
|
|
}
|
|
|
|
delete [] V;
|
|
|
|
return true;
|
|
}
|
|
|
|
// -------------------------------------------------------------
|
|
|
|
QuasiDelaunayTriangulator::QuasiDelaunayTriangulator()
|
|
{
|
|
}
|
|
|
|
QuasiDelaunayTriangulator::~QuasiDelaunayTriangulator()
|
|
{
|
|
}
|
|
|
|
bool QuasiDelaunayTriangulator::Triangulate()
|
|
{
|
|
if (EarClippingTriangulator::Triangulate() == false)
|
|
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<unsigned long, unsigned long>, std::vector<unsigned long> > aEdge2Face;
|
|
for (std::vector<MeshFacet>::iterator pI = _facets.begin(); pI != _facets.end(); pI++) {
|
|
for (int i = 0; i < 3; i++) {
|
|
unsigned long ulPt0 = std::min<unsigned long>(pI->_aulPoints[i], pI->_aulPoints[(i+1)%3]);
|
|
unsigned long ulPt1 = std::max<unsigned long>(pI->_aulPoints[i], pI->_aulPoints[(i+1)%3]);
|
|
// ignore borderlines of the polygon
|
|
if ((ulPt1-ulPt0)%(_points.size()-1) > 1)
|
|
aEdge2Face[std::pair<unsigned long, unsigned long>(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<unsigned long, unsigned long> > aEdgeList;
|
|
std::map<std::pair<unsigned long, unsigned long>, std::vector<unsigned long> >::iterator pE;
|
|
for (pE = aEdge2Face.begin(); pE != aEdge2Face.end(); ++pE)
|
|
aEdgeList.push_back(pE->first);
|
|
|
|
// to be sure to avoid an endless loop
|
|
unsigned long 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<unsigned long, unsigned long> 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<unsigned long, unsigned long>, std::vector<unsigned long> >::iterator it;
|
|
// first facet
|
|
unsigned long ulPt0 = std::min<unsigned long>(rF1._aulPoints[i], rF1._aulPoints[(i+1)%3]);
|
|
unsigned long ulPt1 = std::max<unsigned long>(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<unsigned long>(rF2._aulPoints[i], rF2._aulPoints[(i+1)%3]);
|
|
ulPt1 = std::max<unsigned long>(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
|
|
unsigned long ulPt0 = std::min<unsigned long>(rF1._aulPoints[(side1+1)%3], rF2._aulPoints[(side2+1)%3]);
|
|
unsigned long ulPt1 = std::max<unsigned long>(rF1._aulPoints[(side1+1)%3], rF2._aulPoints[(side2+1)%3]);
|
|
std::pair<unsigned long, unsigned long> aNewEdge = std::make_pair(ulPt0, ulPt1);
|
|
aEdge2Face[aNewEdge] = pE->second;
|
|
aEdge2Face.erase(pE);
|
|
}
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
// -------------------------------------------------------------
|
|
|
|
namespace MeshCore {
|
|
namespace Triangulation {
|
|
struct Vertex2d_Less : public std::binary_function<const Base::Vector3f&, const Base::Vector3f&, bool>
|
|
{
|
|
bool operator()(const Base::Vector3f& p, const Base::Vector3f& q) const
|
|
{
|
|
if (fabs(p.x - q.x) < MeshDefinitions::_fMinPointDistanceD1) {
|
|
if (fabs(p.y - q.y) < MeshDefinitions::_fMinPointDistanceD1) {
|
|
return false; } else return p.y < q.y;
|
|
} else return p.x < q.x; return true;
|
|
}
|
|
};
|
|
struct Vertex2d_EqualTo : public std::binary_function<const Base::Vector3f&, const Base::Vector3f&, bool>
|
|
{
|
|
bool operator()(const Base::Vector3f& p, const Base::Vector3f& q) const
|
|
{
|
|
if (fabs(p.x - q.x) < MeshDefinitions::_fMinPointDistanceD1
|
|
&& fabs(p.y - q.y) < MeshDefinitions::_fMinPointDistanceD1)
|
|
return true; return false;
|
|
}
|
|
};
|
|
}
|
|
}
|
|
|
|
DelaunayTriangulator::DelaunayTriangulator()
|
|
{
|
|
}
|
|
|
|
DelaunayTriangulator::~DelaunayTriangulator()
|
|
{
|
|
}
|
|
|
|
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 (std::vector<Base::Vector3f>::iterator it = _points.begin(); it != _points.end(); it++) {
|
|
akVertex.push_back(Wm4::Vector2d(it->x, it->y));
|
|
}
|
|
|
|
Wm4::Delaunay2d del(akVertex.size(), &(akVertex[0]), 0.001, false, Wm4::Query::QT_INT64);
|
|
int iTQuantity = del.GetSimplexQuantity();
|
|
std::vector<int> aiTVertex(3*iTQuantity);
|
|
size_t uiSize = 3*iTQuantity*sizeof(int);
|
|
Wm4::System::Memcpy(&(aiTVertex[0]),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 = 0;
|
|
del.GetHull(iEQuantity,aiIndex);
|
|
int iUniqueVQuantity = del.GetUniqueVertexQuantity();
|
|
int iTVerify = 2*iUniqueVQuantity - 2 - iEQuantity;
|
|
(void)iTVerify; // avoid warning in release build
|
|
bool succeeded = (iTVerify == iTQuantity);
|
|
int iEVerify = 3*iUniqueVQuantity - 3 - iEQuantity;
|
|
(void)iEVerify; // avoid warning about unused variable
|
|
delete[] aiIndex;
|
|
|
|
MeshGeomFacet triangle;
|
|
MeshFacet facet;
|
|
for (int i = 0; i < iTQuantity; i++) {
|
|
for (int j=0; j<3; j++) {
|
|
facet._aulPoints[j] = aiTVertex[3*i+j];
|
|
triangle._aclPoints[j].x = (float)akVertex[aiTVertex[3*i+j]].X();
|
|
triangle._aclPoints[j].y = (float)akVertex[aiTVertex[3*i+j]].Y();
|
|
}
|
|
|
|
_triangles.push_back(triangle);
|
|
_facets.push_back(facet);
|
|
}
|
|
|
|
return succeeded;
|
|
}
|
|
|
|
// -------------------------------------------------------------
|
|
|
|
FlatTriangulator::FlatTriangulator()
|
|
{
|
|
}
|
|
|
|
FlatTriangulator::~FlatTriangulator()
|
|
{
|
|
}
|
|
|
|
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)
|
|
{
|
|
}
|
|
|
|
ConstraintDelaunayTriangulator::~ConstraintDelaunayTriangulator()
|
|
{
|
|
}
|
|
|
|
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
|