From 949a5b25c2dbe83b2ba657b4d2bff284d9e16ed1 Mon Sep 17 00:00:00 2001
From: rinkk <karsten.rink@ufz.de>
Date: Fri, 28 May 2021 18:01:49 +0200
Subject: [PATCH] [FileIO] Adding EarClipping triangulation ... again
---
Applications/FileIO/Legacy/createSurface.cpp | 52 ++++
Applications/FileIO/Legacy/createSurface.h | 10 +-
GeoLib/EarClippingTriangulation.cpp | 306 +++++++++++++++++++
GeoLib/EarClippingTriangulation.h | 66 ++++
GeoLib/Polygon.cpp | 14 +
GeoLib/Polygon.h | 3 +
6 files changed, 449 insertions(+), 2 deletions(-)
create mode 100644 GeoLib/EarClippingTriangulation.cpp
create mode 100644 GeoLib/EarClippingTriangulation.h
diff --git a/Applications/FileIO/Legacy/createSurface.cpp b/Applications/FileIO/Legacy/createSurface.cpp
index bd956f30260..c0d6c6b62b5 100644
--- a/Applications/FileIO/Legacy/createSurface.cpp
+++ b/Applications/FileIO/Legacy/createSurface.cpp
@@ -17,11 +17,13 @@
#include "Applications/FileIO/Gmsh/GMSHInterface.h"
#include "BaseLib/Logging.h"
#include "BaseLib/StringTools.h"
+#include "GeoLib/EarClippingTriangulation.h"
#include "GeoLib/GEOObjects.h"
#include "GeoLib/Point.h"
#include "GeoLib/Polygon.h"
#include "GeoLib/Polyline.h"
#include "GeoLib/Surface.h"
+#include "GeoLib/Triangle.h"
#include "MeshLib/IO/readMeshFromFile.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/convertMeshToGeo.h"
@@ -122,4 +124,54 @@ bool createSurface(GeoLib::Polyline const& ply,
return true;
}
+GeoLib::Surface* createSurfaceWithEarClipping(GeoLib::Polyline const& line)
+{
+ if (!line.isClosed())
+ {
+ WARN("Error in Surface::createSurface() - Polyline is not closed.");
+ return nullptr;
+ }
+
+ if (line.getNumberOfPoints() <= 2)
+ {
+ WARN(
+ "Error in Surface::createSurface() - Polyline consists of less "
+ "than three points and therefore cannot be triangulated.");
+ return nullptr;
+ }
+
+ // create empty surface
+ GeoLib::Surface* sfc(new GeoLib::Surface(line.getPointsVec()));
+ GeoLib::Polygon* polygon(new GeoLib::Polygon(line));
+ std::list<GeoLib::Polygon*> const& list_of_simple_polygons =
+ polygon->computeListOfSimplePolygons();
+
+ for (std::list<GeoLib::Polygon*>::const_iterator simple_polygon_it(
+ list_of_simple_polygons.begin());
+ simple_polygon_it != list_of_simple_polygons.end();
+ ++simple_polygon_it)
+ {
+ std::list<GeoLib::Triangle> triangles;
+ GeoLib::EarClippingTriangulation(*simple_polygon_it, triangles);
+
+ // add Triangles to Surface
+ std::list<GeoLib::Triangle>::const_iterator it(triangles.begin());
+ while (it != triangles.end())
+ {
+ sfc->addTriangle((*it)[0], (*it)[1], (*it)[2]);
+ ++it;
+ }
+ }
+ // delete polygon;
+ if (sfc->getNumberOfTriangles() == 0)
+ {
+ WARN(
+ "Surface::createSurface(): Triangulation does not contain any "
+ "triangle.");
+ delete sfc;
+ return nullptr;
+ }
+ return sfc;
+}
+
} // namespace FileIO
diff --git a/Applications/FileIO/Legacy/createSurface.h b/Applications/FileIO/Legacy/createSurface.h
index 24a0347cd45..8b3fd65c8e6 100644
--- a/Applications/FileIO/Legacy/createSurface.h
+++ b/Applications/FileIO/Legacy/createSurface.h
@@ -15,12 +15,13 @@
namespace GeoLib
{
class Polyline;
+class Surface;
class GEOObjects;
-}
+} // namespace GeoLib
namespace FileIO
{
-/// Creates a plane surface from the given polyline. The polyline have to be
+/// Creates a plane surface from the given polyline. The polyline has to be
/// closed, i.e. the first and the last point have to be the identical. The
/// triangulation of the polyline is done by the finite element meshing tool
/// Gmsh. Finally, the resulting mesh is converted into a GeoLib::Surface which
@@ -30,4 +31,9 @@ bool createSurface(GeoLib::Polyline const& ply,
GeoLib::GEOObjects& geometries,
std::string const& geometry_name,
std::string const& gmsh_binary);
+
+/// Creates a plane surface from the given polyline. The polyline has to be
+/// closed, i.e. the first and the last point have to be the identical. The
+/// triangulation of the polyline is done via a simple ear clipping algorithm.
+GeoLib::Surface* createSurfaceWithEarClipping(GeoLib::Polyline const& line);
} // namespace FileIO
diff --git a/GeoLib/EarClippingTriangulation.cpp b/GeoLib/EarClippingTriangulation.cpp
new file mode 100644
index 00000000000..3b1bdcf1007
--- /dev/null
+++ b/GeoLib/EarClippingTriangulation.cpp
@@ -0,0 +1,306 @@
+/**
+ * \file
+ * \author Thomas Fischer
+ * \date 2011-02-23
+ * \brief Implementation of the EarClippingTriangulation class.
+ *
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "EarClippingTriangulation.h"
+
+//#include "BaseLib/uniqueInsert.h"
+
+#include "MathLib/GeometricBasics.h"
+
+#include "Polygon.h"
+#include "Triangle.h"
+#include "Point.h"
+
+template <typename Container>
+void uniquePushBack(Container& container,
+ typename Container::value_type const& element)
+{
+ if (std::find(container.begin(), container.end(), element) ==
+ container.end())
+ container.push_back(element);
+}
+
+namespace GeoLib
+{
+EarClippingTriangulation::EarClippingTriangulation(const GeoLib::Polygon* polygon,
+ std::list<GeoLib::Triangle> &triangles, bool rot)
+{
+ copyPolygonPoints (polygon);
+
+ if (rot) {
+ rotatePointsToXY (_pnts);
+ ensureCWOrientation ();
+ }
+
+ initVertexList ();
+ initLists ();
+ clipEars ();
+
+ std::vector<GeoLib::Point*> const& ref_pnts_vec (polygon->getPointsVec());
+ std::list<GeoLib::Triangle>::const_iterator it (_triangles.begin());
+ if (_original_orient == GeoLib::CW) {
+ while (it != _triangles.end()) {
+ const std::size_t i0 (polygon->getPointID ((*it)[0]));
+ const std::size_t i1 (polygon->getPointID ((*it)[1]));
+ const std::size_t i2 (polygon->getPointID ((*it)[2]));
+ triangles.push_back (GeoLib::Triangle (ref_pnts_vec, i0, i1, i2));
+ ++it;
+ }
+ } else {
+ std::size_t n_pnts (_pnts.size() - 1);
+ while (it != _triangles.end()) {
+ const std::size_t i0 (polygon->getPointID (n_pnts - (*it)[0]));
+ const std::size_t i1 (polygon->getPointID (n_pnts - (*it)[1]));
+ const std::size_t i2 (polygon->getPointID (n_pnts - (*it)[2]));
+ triangles.push_back (GeoLib::Triangle (ref_pnts_vec, i0, i1, i2));
+ ++it;
+ }
+ }
+}
+
+EarClippingTriangulation::~EarClippingTriangulation()
+{
+ const std::size_t n_pnts (_pnts.size());
+ for (std::size_t k(0); k<n_pnts; k++) {
+ delete _pnts[k];
+ }
+}
+
+void EarClippingTriangulation::copyPolygonPoints (const GeoLib::Polygon* polygon)
+{
+ // copy points - last point is identical to the first
+ std::size_t n_pnts (polygon->getNumberOfPoints() - 1);
+ for (std::size_t k(0); k < n_pnts; k++) {
+ _pnts.push_back (new GeoLib::Point (*(polygon->getPoint(k))));
+ }
+}
+
+void EarClippingTriangulation::ensureCWOrientation ()
+{
+ std::size_t n_pnts (_pnts.size());
+ // get the left most upper point
+ std::size_t min_x_max_y_idx (0); // for orientation check
+ for (std::size_t k(0); k<n_pnts; k++) {
+ if ((*(_pnts[k]))[0] <= (*(_pnts[min_x_max_y_idx]))[0]) {
+ if ((*(_pnts[k]))[0] < (*(_pnts[min_x_max_y_idx]))[0]) {
+ min_x_max_y_idx = k;
+ } else {
+ if ((*(_pnts[k]))[1] > (*(_pnts[min_x_max_y_idx]))[1]) {
+ min_x_max_y_idx = k;
+ }
+ }
+ }
+ }
+ // determine orientation
+ if (0 < min_x_max_y_idx && min_x_max_y_idx < n_pnts-1) {
+ _original_orient = GeoLib::getOrientation (
+ *_pnts[min_x_max_y_idx-1], *_pnts[min_x_max_y_idx], *_pnts[min_x_max_y_idx+1]);
+ } else {
+ if (0 == min_x_max_y_idx) {
+ _original_orient = GeoLib::getOrientation (*_pnts[n_pnts-1], *_pnts[0], *_pnts[1]);
+ } else {
+ _original_orient = GeoLib::getOrientation (*_pnts[n_pnts-2], *_pnts[n_pnts-1], *_pnts[0]);
+ }
+ }
+ if (_original_orient == GeoLib::CCW) {
+ // switch orientation
+ for (std::size_t k(0); k<n_pnts/2; k++) {
+ std::swap (_pnts[k], _pnts[n_pnts-1-k]);
+ }
+ }
+}
+
+bool EarClippingTriangulation::isEar(std::size_t v0, std::size_t v1, std::size_t v2) const
+{
+ for (std::list<std::size_t>::const_iterator it (_vertex_list.begin ());
+ it != _vertex_list.end(); ++it) {
+ if (*it != v0 && *it != v1 && *it != v2) {
+ if (MathLib::isPointInTriangle (*_pnts[*it], *_pnts[v0], *_pnts[v1], *_pnts[v2])) {
+ return false;
+ }
+ }
+ }
+ return true;
+}
+
+void EarClippingTriangulation::initVertexList ()
+{
+ std::size_t n_pnts (_pnts.size());
+ for (std::size_t k(0); k<n_pnts; k++)
+ _vertex_list.push_back (k);
+}
+
+void EarClippingTriangulation::initLists ()
+{
+ // go through points checking ccw, cw or collinear order and identifying ears
+ std::list<std::size_t>::iterator it (_vertex_list.begin()), prev(_vertex_list.end()), next;
+ --prev;
+ next = it;
+ ++next;
+ GeoLib::Orientation orientation;
+ bool first_run(true); // saves special handling of the last case with identical code
+ while (_vertex_list.size() >= 3 && first_run) {
+ if (next == _vertex_list.end()) {
+ first_run = false;
+ next = _vertex_list.begin();
+ }
+ orientation = getOrientation (*_pnts[*prev], *_pnts[*it], *_pnts[*next]);
+ if (orientation == GeoLib::COLLINEAR) {
+ WARN(
+ "EarClippingTriangulation::initLists(): collinear points "
+ "({:f}, {:f}, {:f}), ({:f}, {:f}, {:f}), ({:f}, {:f}, {:f})",
+ (*_pnts[*prev])[0], (*_pnts[*prev])[1], (*_pnts[*prev])[2],
+ (*_pnts[*it])[0], (*_pnts[*it])[1], (*_pnts[*it])[2],
+ (*_pnts[*next])[0], (*_pnts[*next])[1], (*_pnts[*next])[2]);
+ it = _vertex_list.erase (it);
+ ++next;
+ } else {
+ if (orientation == GeoLib::CW) {
+ _convex_vertex_list.push_back (*it);
+ if (isEar (*prev, *it, *next))
+ _ear_list.push_back (*it);
+ }
+ prev = it;
+ it = next;
+ ++next;
+ }
+ }
+}
+
+void EarClippingTriangulation::clipEars()
+{
+ std::list<std::size_t>::iterator it, prev, next;
+ // *** clip an ear
+ while (_vertex_list.size() > 3) {
+ // pop ear from list
+ std::size_t ear = _ear_list.front();
+ _ear_list.pop_front();
+ // remove ear tip from _convex_vertex_list
+ _convex_vertex_list.remove(ear);
+
+ // remove ear from vertex_list, apply changes to _ear_list, _convex_vertex_list
+ bool nfound(true);
+ it = _vertex_list.begin();
+ prev = _vertex_list.end();
+ --prev;
+ while (nfound && it != _vertex_list.end()) {
+ if (*it == ear) {
+ nfound = false;
+ it = _vertex_list.erase(it); // remove ear tip
+ next = it;
+ if (next == _vertex_list.end()) {
+ next = _vertex_list.begin();
+ prev = _vertex_list.end();
+ --prev;
+ }
+ // add triangle
+ _triangles.push_back(GeoLib::Triangle(_pnts, *prev, *next, ear));
+
+ // check the orientation of prevprev, prev, next
+ std::list<std::size_t>::iterator prevprev;
+ if (prev == _vertex_list.begin()) {
+ prevprev = _vertex_list.end();
+ } else {
+ prevprev = prev;
+ }
+ --prevprev;
+
+ // apply changes to _convex_vertex_list and _ear_list looking "backward"
+ GeoLib::Orientation orientation = GeoLib::getOrientation(*_pnts[*prevprev], *_pnts[*prev],
+ *_pnts[*next]);
+ if (orientation == GeoLib::CW) {
+ uniquePushBack(_convex_vertex_list, *prev);
+ // prev is convex
+ if (isEar(*prevprev, *prev, *next)) {
+ // prev is an ear tip
+ uniquePushBack(_ear_list, *prev);
+ } else {
+ // if necessary remove prev
+ _ear_list.remove(*prev);
+ }
+ } else {
+ // prev is not convex => reflex or collinear
+ _convex_vertex_list.remove(*prev);
+ _ear_list.remove(*prev);
+ if (orientation == GeoLib::COLLINEAR) {
+ prev = _vertex_list.erase(prev);
+ if (prev == _vertex_list.begin()) {
+ prev = _vertex_list.end();
+ --prev;
+ } else {
+ --prev;
+ }
+ }
+ }
+
+ // check the orientation of prev, next, nextnext
+ std::list<std::size_t>::iterator nextnext,
+ help_it(_vertex_list.end());
+ --help_it;
+ if (next == help_it) {
+ nextnext = _vertex_list.begin();
+ } else {
+ nextnext = next;
+ ++nextnext;
+ }
+
+ // apply changes to _convex_vertex_list and _ear_list looking "forward"
+ orientation = getOrientation(*_pnts[*prev], *_pnts[*next],
+ *_pnts[*nextnext]);
+ if (orientation == GeoLib::CW) {
+ uniquePushBack(_convex_vertex_list, *next);
+ // next is convex
+ if (isEar(*prev, *next, *nextnext)) {
+ // next is an ear tip
+ uniquePushBack(_ear_list, *next);
+ } else {
+ // if necessary remove *next
+ _ear_list.remove(*next);
+ }
+ } else {
+ // next is not convex => reflex or collinear
+ _convex_vertex_list.remove(*next);
+ _ear_list.remove(*next);
+ if (orientation == GeoLib::COLLINEAR) {
+ next = _vertex_list.erase(next);
+ if (next == _vertex_list.end())
+ next = _vertex_list.begin();
+ }
+ }
+ } else {
+ prev = it;
+ ++it;
+ }
+ }
+
+ }
+
+ // add last triangle
+ next = _vertex_list.begin();
+ prev = next;
+ ++next;
+ if (next == _vertex_list.end())
+ return;
+ it = next;
+ ++next;
+ if (next == _vertex_list.end())
+ return;
+
+ if (getOrientation(*_pnts[*prev], *_pnts[*it], *_pnts[*next]) == GeoLib::CCW)
+ _triangles.push_back(GeoLib::Triangle(_pnts, *prev, *it, *next));
+ else
+ _triangles.push_back(GeoLib::Triangle(_pnts, *prev, *next, *it));
+}
+
+} // end namespace GeoLib
diff --git a/GeoLib/EarClippingTriangulation.h b/GeoLib/EarClippingTriangulation.h
new file mode 100644
index 00000000000..4eec2a09ba3
--- /dev/null
+++ b/GeoLib/EarClippingTriangulation.h
@@ -0,0 +1,66 @@
+/**
+ * \file
+ * \author Thomas Fischer
+ * \date 2011-02-23
+ * \brief Definition of the EarClippingTriangulation class.
+ *
+ * \copyright
+ * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef EARCLIPPINGTRIANGULATION_H_
+#define EARCLIPPINGTRIANGULATION_H_
+
+#include <list>
+#include <vector>
+
+#include "AnalyticalGeometry.h"
+
+namespace GeoLib
+{
+
+class Polygon;
+class Triangle;
+
+class EarClippingTriangulation
+{
+public:
+ EarClippingTriangulation(const GeoLib::Polygon* ply,
+ std::list<GeoLib::Triangle> &triangles,
+ bool rot = true);
+ virtual ~EarClippingTriangulation();
+private:
+ /**
+ * copies the points of the polygon to the vector _pnts
+ */
+ inline void copyPolygonPoints (const GeoLib::Polygon* polygon);
+ inline void ensureCWOrientation ();
+
+ inline bool isEar(std::size_t v0, std::size_t v1, std::size_t v2) const;
+
+ inline void initVertexList ();
+ inline void initLists ();
+ inline void clipEars ();
+
+ /**
+ * a copy of the polygon points
+ */
+ std::vector<GeoLib::Point*> _pnts;
+ std::list<std::size_t> _vertex_list;
+ std::list<std::size_t> _convex_vertex_list;
+ std::list<std::size_t> _ear_list;
+
+ /**
+ * triangles of the triangulation (maybe in the wrong orientation)
+ */
+ std::list<GeoLib::Triangle> _triangles;
+
+ GeoLib::Orientation _original_orient;
+};
+} // end namespace GeoLib
+
+#endif /* EARCLIPPINGTRIANGULATION_H_ */
diff --git a/GeoLib/Polygon.cpp b/GeoLib/Polygon.cpp
index c74b550fc0f..dbcdaf8043a 100644
--- a/GeoLib/Polygon.cpp
+++ b/GeoLib/Polygon.cpp
@@ -615,4 +615,18 @@ bool operator==(Polygon const& lhs, Polygon const& rhs)
return false;
}
+std::list<Polygon*> const& Polygon::computeListOfSimplePolygons()
+{
+ splitPolygonAtPoint(_simple_polygon_list.begin());
+ splitPolygonAtIntersection(_simple_polygon_list.begin());
+
+ for (std::list<GeoLib::Polygon*>::iterator it(_simple_polygon_list.begin());
+ it != _simple_polygon_list.end();
+ ++it)
+ {
+ (*it)->initialise();
+ }
+ return _simple_polygon_list;
+}
+
} // end namespace GeoLib
diff --git a/GeoLib/Polygon.h b/GeoLib/Polygon.h
index f0a5d8b11e3..09f4ad3bc8e 100644
--- a/GeoLib/Polygon.h
+++ b/GeoLib/Polygon.h
@@ -100,6 +100,9 @@ public:
GeoLib::Point& intersection_pnt,
std::size_t& seg_num) const;
+ /// Subdivides a self-intersecting polygon into a list of non-intersecting shapes.
+ std::list<Polygon*> const& computeListOfSimplePolygons();
+
friend bool operator==(Polygon const& lhs, Polygon const& rhs);
private:
/**
--
GitLab