[FileIO] Adding EarClipping triangulation ... again

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

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

GeoLib/EarClippingTriangulation.cpp

+/**
+ * \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

GeoLib/EarClippingTriangulation.h

+/**
+ * \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_ */

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

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:
     /**