diff --git a/common/CMakeLists.txt b/common/CMakeLists.txt
index b4d676ef67..c701b467cc 100644
--- a/common/CMakeLists.txt
+++ b/common/CMakeLists.txt
@@ -1,4 +1,3 @@
-
include_directories(BEFORE ${INC_BEFORE})
include_directories(
./dialogs
@@ -37,6 +36,7 @@ set(GAL_SRCS
gal/stroke_font.cpp
gal/color4d.cpp
view/wx_view_controls.cpp
+ geometry/hetriang.cpp
# OpenGL GAL
gal/opengl/opengl_gal.cpp
diff --git a/include/ttl/halfedge/hedart.h b/include/ttl/halfedge/hedart.h
new file mode 100644
index 0000000000..f85678963a
--- /dev/null
+++ b/include/ttl/halfedge/hedart.h
@@ -0,0 +1,150 @@
+/*
+ * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
+ * Applied Mathematics, Norway.
+ *
+ * Contact information: E-mail: tor.dokken@sintef.no
+ * SINTEF ICT, Department of Applied Mathematics,
+ * P.O. Box 124 Blindern,
+ * 0314 Oslo, Norway.
+ *
+ * This file is part of TTL.
+ *
+ * TTL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * TTL 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public
+ * License along with TTL. If not, see
+ * .
+ *
+ * In accordance with Section 7(b) of the GNU Affero General Public
+ * License, a covered work must retain the producer line in every data
+ * file that is created or manipulated using TTL.
+ *
+ * Other Usage
+ * You can be released from the requirements of the license by purchasing
+ * a commercial license. Buying such a license is mandatory as soon as you
+ * develop commercial activities involving the TTL library without
+ * disclosing the source code of your own applications.
+ *
+ * This file may be used in accordance with the terms contained in a
+ * written agreement between you and SINTEF ICT.
+ */
+
+#ifndef _HALF_EDGE_DART_
+#define _HALF_EDGE_DART_
+
+
+#include
+
+
+namespace hed {
+
+
+ //------------------------------------------------------------------------------------------------
+ // Dart class for the half-edge data structure
+ //------------------------------------------------------------------------------------------------
+
+ /** \class Dart
+ * \brief \b %Dart class for the half-edge data structure.
+ *
+ * See \ref api for a detailed description of how the member functions
+ * should be implemented.
+ */
+
+ class Dart {
+
+ EdgePtr edge_;
+ bool dir_; // true if dart is counterclockwise in face
+
+ public:
+ /// Default constructor
+ Dart() { dir_ = true; }
+
+ /// Constructor
+ Dart(const EdgePtr& edge, bool dir = true) { edge_ = edge; dir_ = dir; }
+
+ /// Copy constructor
+ Dart(const Dart& dart) { edge_ = dart.edge_; dir_ = dart.dir_; }
+
+ /// Destructor
+ ~Dart() {}
+
+ /// Assignment operator
+ Dart& operator = (const Dart& dart) {
+ if (this == &dart)
+ return *this;
+ edge_ = dart.edge_;
+ dir_ = dart.dir_;
+ return *this;
+ }
+
+ /// Comparing dart objects
+ bool operator==(const Dart& dart) const {
+ if (dart.edge_ == edge_ && dart.dir_ == dir_)
+ return true;
+ return false;
+ }
+
+ /// Comparing dart objects
+ bool operator!=(const Dart& dart) const {
+ return !(dart==*this);
+ }
+
+ /// Maps the dart to a different node
+ Dart& alpha0() { dir_ = !dir_; return *this; }
+
+ /// Maps the dart to a different edge
+ Dart& alpha1() {
+ if (dir_) {
+ edge_ = edge_->getNextEdgeInFace()->getNextEdgeInFace();
+ dir_ = false;
+ }
+ else {
+ edge_ = edge_->getNextEdgeInFace();
+ dir_ = true;
+ }
+ return *this;
+ }
+
+ /// Maps the dart to a different triangle. \b Note: the dart is not changed if it is at the boundary!
+ Dart& alpha2() {
+ if (edge_->getTwinEdge()) {
+ edge_ = edge_->getTwinEdge();
+ dir_ = !dir_;
+ }
+ // else, the dart is at the boundary and should not be changed
+ return *this;
+ }
+
+
+ // Utilities not required by TTL
+ // -----------------------------
+
+ /** @name Utilities not required by TTL */
+ //@{
+
+ void init(const EdgePtr& edge, bool dir = true) { edge_ = edge; dir_ = dir; }
+
+ double x() const { return getNode()->GetX(); } // x-coordinate of source node
+ double y() const { return getNode()->GetY(); } // y-coordinate of source node
+
+ bool isCounterClockWise() const { return dir_; }
+
+ const NodePtr& getNode() const { return dir_ ? edge_->getSourceNode() : edge_->getTargetNode(); }
+ const NodePtr& getOppositeNode() const { return dir_ ? edge_->getTargetNode() : edge_->getSourceNode(); }
+ EdgePtr& getEdge() { return edge_; }
+
+ //@} // End of Utilities not required by TTL
+
+ };
+
+}; // End of hed namespace
+
+#endif
diff --git a/include/ttl/halfedge/hetraits.h b/include/ttl/halfedge/hetraits.h
new file mode 100644
index 0000000000..e25e993e0e
--- /dev/null
+++ b/include/ttl/halfedge/hetraits.h
@@ -0,0 +1,300 @@
+/*
+ * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
+ * Applied Mathematics, Norway.
+ *
+ * Contact information: E-mail: tor.dokken@sintef.no
+ * SINTEF ICT, Department of Applied Mathematics,
+ * P.O. Box 124 Blindern,
+ * 0314 Oslo, Norway.
+ *
+ * This file is part of TTL.
+ *
+ * TTL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * TTL 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public
+ * License along with TTL. If not, see
+ * .
+ *
+ * In accordance with Section 7(b) of the GNU Affero General Public
+ * License, a covered work must retain the producer line in every data
+ * file that is created or manipulated using TTL.
+ *
+ * Other Usage
+ * You can be released from the requirements of the license by purchasing
+ * a commercial license. Buying such a license is mandatory as soon as you
+ * develop commercial activities involving the TTL library without
+ * disclosing the source code of your own applications.
+ *
+ * This file may be used in accordance with the terms contained in a
+ * written agreement between you and SINTEF ICT.
+ */
+
+#ifndef _HALF_EDGE_TRAITS_
+#define _HALF_EDGE_TRAITS_
+
+
+#include
+#include
+
+
+namespace hed {
+
+
+ //------------------------------------------------------------------------------------------------
+ // Traits class for the half-edge data structure
+ //------------------------------------------------------------------------------------------------
+
+ /** \struct TTLtraits
+ * \brief \b Traits class (static struct) for the half-edge data structure.
+ *
+ * The member functions are those required by different function templates
+ * in the TTL. Documentation is given here to explain what actions
+ * should be carried out on the actual data structure as required by the functions
+ * in the \ref ttl namespace.
+ *
+ * The source code of \c %HeTraits.h shows how the traits class is implemented for the
+ * half-edge data structure.
+ *
+ * \see \ref api
+ *
+ */
+
+ struct TTLtraits {
+
+ // The actual triangulation object
+ static Triangulation* triang_;
+
+ /** The floating point type used in calculations
+ * involving scalar products and cross products.
+ */
+ typedef double real_type;
+
+
+ //----------------------------------------------------------------------------------------------
+ // ------------------------------- Geometric Predicates Group ---------------------------------
+ //----------------------------------------------------------------------------------------------
+
+ /** @name Geometric Predicates */
+ //@{
+
+ //----------------------------------------------------------------------------------------------
+ /** Scalar product between two 2D vectors represented as darts.\n
+ *
+ * ttl_util::scalarProduct2d can be used.
+ */
+ static real_type scalarProduct2d(const Dart& v1, const Dart& v2) {
+ Dart v10 = v1; v10.alpha0();
+ Dart v20 = v2; v20.alpha0();
+ return ttl_util::scalarProduct2d(v10.x()-v1.x(), v10.y()-v1.y(),
+ v20.x()-v2.x(), v20.y()-v2.y());
+ }
+
+
+ //----------------------------------------------------------------------------------------------
+ /** Scalar product between two 2D vectors.
+ * The first vector is represented by a dart \e v, and the second
+ * vector has direction from the source node of \e v to the point \e p.\n
+ *
+ * ttl_util::scalarProduct2d can be used.
+ */
+ static real_type scalarProduct2d(const Dart& v, const NodePtr& p) {
+ Dart d0 = v; d0.alpha0();
+ return ttl_util::scalarProduct2d(d0.x() - v.x(), d0.y() - v.y(),
+ p->GetX() - v.x(), p->GetY() - v.y());
+ }
+
+
+ //----------------------------------------------------------------------------------------------
+ /** Cross product between two vectors in the plane represented as darts.
+ * The z-component of the cross product is returned.\n
+ *
+ * ttl_util::crossProduct2d can be used.
+ */
+ static real_type crossProduct2d(const Dart& v1, const Dart& v2) {
+ Dart v10 = v1; v10.alpha0();
+ Dart v20 = v2; v20.alpha0();
+ return ttl_util::crossProduct2d(v10.x()-v1.x(), v10.y()-v1.y(),
+ v20.x()-v2.x(), v20.y()-v2.y());
+ }
+
+
+ //----------------------------------------------------------------------------------------------
+ /** Cross product between two vectors in the plane.
+ * The first vector is represented by a dart \e v, and the second
+ * vector has direction from the source node of \e v to the point \e p.
+ * The z-component of the cross product is returned.\n
+ *
+ * ttl_util::crossProduct2d can be used.
+ */
+ static real_type crossProduct2d(const Dart& v, const NodePtr& p) {
+ Dart d0 = v; d0.alpha0();
+ return ttl_util::crossProduct2d(d0.x() - v.x(), d0.y() - v.y(),
+ p->GetX() - v.x(), p->GetY() - v.y());
+ }
+
+
+ //----------------------------------------------------------------------------------------------
+ /** Let \e n1 and \e n2 be the nodes associated with two darts, and let \e p
+ * be a point in the plane. Return a positive value if \e n1, \e n2,
+ * and \e p occur in counterclockwise order; a negative value if they occur
+ * in clockwise order; and zero if they are collinear.
+ */
+ static real_type orient2d(const Dart& n1, const Dart& n2, const NodePtr& p) {
+ real_type pa[2]; real_type pb[2]; real_type pc[2];
+ pa[0] = n1.x(); pa[1] = n1.y();
+ pb[0] = n2.x(); pb[1] = n2.y();
+ pc[0] = p->GetX(); pc[1] = p->GetY();
+ return ttl_util::orient2dfast(pa, pb, pc);
+ }
+
+
+ //----------------------------------------------------------------------------------------------
+ /** This is the same predicate as represented with the function above,
+ * but with a slighty different interface:
+ * The last parameter is given as a dart where the source node of the dart
+ * represents a point in the plane.
+ * This function is required for constrained triangulation.
+ */
+ static real_type orient2d(const Dart& n1, const Dart& n2, const Dart& p) {
+ real_type pa[2]; real_type pb[2]; real_type pc[2];
+ pa[0] = n1.x(); pa[1] = n1.y();
+ pb[0] = n2.x(); pb[1] = n2.y();
+ pc[0] = p.x(); pc[1] = p.y();
+ return ttl_util::orient2dfast(pa, pb, pc);
+ }
+
+ //@} // End of Geometric Predicates Group
+
+
+ // A rationale for directing these functions to traits is:
+ // e.g., constraints
+
+ //----------------------------------------------------------------------------------------------
+ /* Checks if the edge associated with \e dart should be swapped
+ * according to the Delaunay criterion.
+ *
+ * \note
+ * This function is also present in the TTL as ttl::swapTestDelaunay.
+ * Thus, the function can be implemented simply as:
+ * \code
+ * { return ttl::swapTestDelaunay(dart); }
+ * \endcode
+ */
+ //static bool swapTestDelaunay(const Dart& dart) {
+ // return ttl::swapTestDelaunay(dart);
+ //}
+
+
+ //----------------------------------------------------------------------------------------------
+ /* Checks if the edge associated with \e dart can be swapped, i.e.,
+ * if the edge is a diagonal in a (strictly) convex quadrilateral.
+ * This function is also present as ttl::swappableEdge.
+ */
+ //static bool swappableEdge(const Dart& dart) {
+ // return ttl::swappableEdge(dart);
+ //}
+
+
+ //----------------------------------------------------------------------------------------------
+ /* Checks if the edge associated with \e dart should be \e fixed, meaning
+ * that it should never be swapped. ??? Use when constraints.
+ */
+ //static bool fixedEdge(const Dart& dart) {
+ // return dart.getEdge()->isConstrained();
+ //}
+
+
+ //----------------------------------------------------------------------------------------------
+ // ----------------------- Functions for Delaunay Triangulation Group -------------------------
+ //----------------------------------------------------------------------------------------------
+
+ /** @name Functions for Delaunay Triangulation */
+ //@{
+
+ //----------------------------------------------------------------------------------------------
+ /** Swaps the edge associated with \e dart in the actual data structure.
+ *
+ *
+ * \image html swapEdge.gif
+ *
+ *
+ * \param dart
+ * Some of the functions require a dart as output.
+ * If this is required by the actual function, the dart should be delivered
+ * back in a position as seen if it was glued to the edge when swapping (rotating)
+ * the edge CCW; see the figure.
+ *
+ * \note
+ * - If the edge is \e constrained, or if it should not be swapped for
+ * some other reason, this function need not do the actual swap of the edge.
+ * - Some functions in TTL require that \c swapEdge is implemented such that
+ * darts outside the quadrilateral are not affected by the swap.
+ */
+ static void swapEdge(Dart& dart) {
+ if (!dart.getEdge()->isConstrained()) triang_->swapEdge(dart.getEdge());
+ }
+
+
+ //----------------------------------------------------------------------------------------------
+ /** Splits the triangle associated with \e dart in the actual data structure into
+ * three new triangles joining at \e point.
+ *
+ *
+ * \image html splitTriangle.gif
+ *
+ *
+ * \param dart
+ * Output: A CCW dart incident with the new node; see the figure.
+ */
+ static void splitTriangle(Dart& dart, NodePtr point) {
+ EdgePtr edge = triang_->splitTriangle(dart.getEdge(), point);
+ dart.init(edge);
+ }
+
+ //@} // End of Functions for Delaunay Triangulation group
+
+
+ //----------------------------------------------------------------------------------------------
+ // --------------------------- Functions for removing nodes Group -----------------------------
+ //----------------------------------------------------------------------------------------------
+
+ /** @name Functions for removing nodes */
+ //@{
+
+ //----------------------------------------------------------------------------------------------
+ /** The reverse operation of TTLtraits::splitTriangle.
+ * This function is only required for functions that involve
+ * removal of interior nodes; see for example ttl::removeInteriorNode.
+ *
+ *
+ * \image html reverse_splitTriangle.gif
+ *
+ */
+ static void reverse_splitTriangle(Dart& dart) {
+ triang_->reverse_splitTriangle(dart.getEdge());
+ }
+
+
+ //----------------------------------------------------------------------------------------------
+ /** Removes a triangle with an edge at the boundary of the triangulation
+ * in the actual data structure
+ */
+ static void removeBoundaryTriangle(Dart& d) {
+ triang_->removeTriangle(d.getEdge());
+ }
+
+ //@} // End of Functions for removing nodes Group
+
+ };
+
+}; // End of hed namespace
+
+#endif
diff --git a/include/ttl/halfedge/hetriang.h b/include/ttl/halfedge/hetriang.h
new file mode 100644
index 0000000000..ccee31bc0c
--- /dev/null
+++ b/include/ttl/halfedge/hetriang.h
@@ -0,0 +1,334 @@
+/*
+ * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
+ * Applied Mathematics, Norway.
+ * Copyright (C) 2013 CERN
+ * @author Maciej Suminski
+ *
+ * Contact information: E-mail: tor.dokken@sintef.no
+ * SINTEF ICT, Department of Applied Mathematics,
+ * P.O. Box 124 Blindern,
+ * 0314 Oslo, Norway.
+ *
+ * This file is part of TTL.
+ *
+ * TTL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * TTL 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public
+ * License along with TTL. If not, see
+ * .
+ *
+ * In accordance with Section 7(b) of the GNU Affero General Public
+ * License, a covered work must retain the producer line in every data
+ * file that is created or manipulated using TTL.
+ *
+ * Other Usage
+ * You can be released from the requirements of the license by purchasing
+ * a commercial license. Buying such a license is mandatory as soon as you
+ * develop commercial activities involving the TTL library without
+ * disclosing the source code of your own applications.
+ *
+ * This file may be used in accordance with the terms contained in a
+ * written agreement between you and SINTEF ICT.
+ */
+
+#ifndef _HE_TRIANG_H_
+#define _HE_TRIANG_H_
+
+
+#define TTL_USE_NODE_ID // Each node gets it's own unique id
+#define TTL_USE_NODE_FLAG // Each node gets a flag (can be set to true or false)
+
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+//--------------------------------------------------------------------------------------------------
+// The half-edge data structure
+//--------------------------------------------------------------------------------------------------
+
+namespace hed {
+ // Helper typedefs
+ class Node;
+ class Edge;
+ typedef boost::shared_ptr NodePtr;
+ typedef boost::shared_ptr EdgePtr;
+ typedef std::vector NodesContainer;
+
+ //------------------------------------------------------------------------------------------------
+ // Node class for data structures
+ //------------------------------------------------------------------------------------------------
+
+ /** \class Node
+ * \brief \b Node class for data structures (Inherits from HandleId)
+ *
+ * \note
+ * - To enable node IDs, TTL_USE_NODE_ID must be defined.
+ * - To enable node flags, TTL_USE_NODE_FLAG must be defined.
+ * - TTL_USE_NODE_ID and TTL_USE_NODE_FLAG should only be enabled if this functionality is
+ * required by the application, because they increase the memory usage for each Node object.
+ */
+
+ class Node {
+
+protected:
+#ifdef TTL_USE_NODE_FLAG
+ /// TTL_USE_NODE_FLAG must be defined
+ bool flag_;
+#endif
+
+#ifdef TTL_USE_NODE_ID
+ /// TTL_USE_NODE_ID must be defined
+ static int id_count;
+
+ /// A unique id for each node (TTL_USE_NODE_ID must be defined)
+ int id_;
+#endif
+
+ int x_, y_;
+
+ unsigned int refCount_;
+
+public:
+ /// Constructor
+ Node( int x = 0, int y = 0 ) :
+#ifdef TTL_USE_NODE_FLAG
+ flag_( false ),
+#endif
+#ifdef TTL_USE_NODE_ID
+ id_( id_count++ ),
+#endif
+ x_( x ), y_( y ), refCount_( 0 ) {}
+
+ /// Destructor
+ ~Node() {}
+
+ /// Returns the x-coordinate
+ int GetX() const { return x_; }
+
+ /// Returns the y-coordinate
+ int GetY() const { return y_; }
+
+#ifdef TTL_USE_NODE_ID
+ /// Returns the id (TTL_USE_NODE_ID must be defined)
+ int Id() const { return id_; }
+#endif
+
+#ifdef TTL_USE_NODE_FLAG
+ /// Sets the flag (TTL_USE_NODE_FLAG must be defined)
+ void SetFlag(bool aFlag) { flag_ = aFlag; }
+
+ /// Returns the flag (TTL_USE_NODE_FLAG must be defined)
+ const bool& GetFlag() const { return flag_; }
+#endif
+
+ void IncRefCount() { refCount_++; }
+ void DecRefCount() { refCount_--; }
+ unsigned int GetRefCount() const { return refCount_; }
+ }; // End of class Node
+
+
+ //------------------------------------------------------------------------------------------------
+ // Edge class in the half-edge data structure
+ //------------------------------------------------------------------------------------------------
+
+ /** \class Edge
+ * \brief \b %Edge class in the in the half-edge data structure.
+ */
+
+ class Edge {
+ public:
+ /// Constructor
+ Edge() : weight_(0)
+ { flags_.isLeadingEdge_ = false; flags_.isConstrained_ = false; }
+
+ /// Destructor
+ virtual ~Edge() {}
+
+ /// Sets the source node
+ void setSourceNode(const NodePtr& node) { sourceNode_ = node; }
+
+ /// Sets the next edge in face
+ void setNextEdgeInFace(const EdgePtr& edge) { nextEdgeInFace_ = edge; }
+
+ /// Sets the twin edge
+ void setTwinEdge(const EdgePtr& edge) { twinEdge_ = edge; }
+
+ /// Sets the edge as a leading edge
+ void setAsLeadingEdge(bool val=true) { flags_.isLeadingEdge_ = val; }
+
+ /// Checks if an edge is a leading edge
+ bool isLeadingEdge() const { return flags_.isLeadingEdge_; }
+
+ /// Sets the edge as a constrained edge
+ void setConstrained(bool val=true) { flags_.isConstrained_ = val;
+ if (twinEdge_) twinEdge_->flags_.isConstrained_ = val; }
+
+ /// Checks if an edge is constrained
+ bool isConstrained() const { return flags_.isConstrained_; }
+
+ /// Returns the twin edge
+ const EdgePtr& getTwinEdge() const { return twinEdge_; };
+
+ /// Returns the next edge in face
+ const EdgePtr& getNextEdgeInFace() const { return nextEdgeInFace_; }
+
+ /// Retuns the source node
+ virtual const NodePtr& getSourceNode() const { return sourceNode_; }
+
+ /// Returns the target node
+ virtual const NodePtr& getTargetNode() const { return getNextEdgeInFace()->getSourceNode(); }
+
+ void setWeight( unsigned int weight ) { weight_ = weight; }
+
+ unsigned int getWeight() const { return weight_; }
+
+ protected:
+ NodePtr sourceNode_;
+ EdgePtr twinEdge_;
+ EdgePtr nextEdgeInFace_;
+ unsigned int weight_;
+
+ struct {
+ bool isLeadingEdge_;
+ bool isConstrained_;
+ } flags_;
+ }; // End of class Edge
+
+
+ /** \class EdgeMST
+ * \brief \b %Specialization of Edge class to be used for Minimum Spanning Tree algorithm.
+ */
+ class EdgeMST : public Edge
+ {
+ private:
+ NodePtr target_;
+
+ public:
+ EdgeMST( const NodePtr& source, const NodePtr& target, unsigned int weight = 0 ) :
+ target_(target)
+ { sourceNode_ = source; weight_ = weight; }
+
+ ~EdgeMST() {};
+
+ /// @copydoc Edge::setSourceNode()
+ const NodePtr& getTargetNode() const { return target_; }
+ };
+
+
+ //------------------------------------------------------------------------------------------------
+ class Dart; // Forward declaration (class in this namespace)
+
+ //------------------------------------------------------------------------------------------------
+ // Triangulation class in the half-edge data structure
+ //------------------------------------------------------------------------------------------------
+
+ /** \class Triangulation
+ * \brief \b %Triangulation class for the half-edge data structure with adaption to TTL.
+ */
+
+ class Triangulation {
+
+ protected:
+ list leadingEdges_; // one half-edge for each arc
+ void addLeadingEdge(EdgePtr& edge) {
+ edge->setAsLeadingEdge();
+ leadingEdges_.push_front( edge );
+ }
+ bool removeLeadingEdgeFromList(EdgePtr& leadingEdge);
+ void cleanAll();
+
+ public:
+ /// Default constructor
+ Triangulation() {}
+
+ /// Copy constructor
+ Triangulation(const Triangulation& tr) {
+ std::cout << "Triangulation: Copy constructor not present - EXIT.";
+ exit(-1);
+ }
+
+ /// Destructor
+ ~Triangulation() { cleanAll(); }
+
+ /// Creates a Delaunay triangulation from a set of points
+ void createDelaunay(NodesContainer::iterator first,
+ NodesContainer::iterator last);
+
+ /// Creates an initial Delaunay triangulation from two enclosing triangles
+ // When using rectangular boundary - loop through all points and expand.
+ // (Called from createDelaunay(...) when starting)
+ EdgePtr initTwoEnclosingTriangles(NodesContainer::iterator first,
+ NodesContainer::iterator last);
+
+
+ // These two functions are required by TTL for Delaunay triangulation
+
+ /// Swaps the edge associated with diagonal
+ void swapEdge(EdgePtr& diagonal);
+
+ /// Splits the triangle associated with edge into three new triangles joining at point
+ EdgePtr splitTriangle(EdgePtr& edge, NodePtr& point);
+
+
+ // Functions required by TTL for removing nodes in a Delaunay triangulation
+
+ /// Removes the boundary triangle associated with edge
+ void removeTriangle(EdgePtr& edge); // boundary triangle required
+
+ /// The reverse operation of removeTriangle
+ void reverse_splitTriangle(EdgePtr& edge);
+
+ /// Creates an arbitrary CCW dart
+ Dart createDart();
+
+ /// Returns a list of "triangles" (one leading half-edge for each triangle)
+ const list& getLeadingEdges() const { return leadingEdges_; }
+
+ /// Returns the number of triangles
+ int noTriangles() const { return (int)leadingEdges_.size(); }
+
+ /// Returns a list of half-edges (one half-edge for each arc)
+ list* getEdges(bool skip_boundary_edges = false) const;
+
+#ifdef TTL_USE_NODE_FLAG
+ /// Sets flag in all the nodes
+ void flagNodes(bool flag) const;
+
+ /// Returns a list of nodes. This function requires TTL_USE_NODE_FLAG to be defined. \see Node.
+ list* getNodes() const;
+#endif
+
+ /// Swaps edges until the triangulation is Delaunay (constrained edges are not swapped)
+ void optimizeDelaunay();
+
+ /// Checks if the triangulation is Delaunay
+ bool checkDelaunay() const;
+
+ /// Returns an arbitrary interior node (as the source node of the returned edge)
+ EdgePtr getInteriorNode() const;
+
+ /// Returns an arbitrary boundary edge
+ EdgePtr getBoundaryEdge() const;
+
+ /// Print edges for plotting with, e.g., gnuplot
+ void printEdges(std::ofstream& os) const;
+
+ }; // End of class Triangulation
+
+
+}; // End of hed namespace
+
+#endif
diff --git a/include/ttl/ttl.h b/include/ttl/ttl.h
new file mode 100644
index 0000000000..003e0d81c6
--- /dev/null
+++ b/include/ttl/ttl.h
@@ -0,0 +1,1907 @@
+/*
+ * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
+ * Applied Mathematics, Norway.
+ *
+ * Contact information: E-mail: tor.dokken@sintef.no
+ * SINTEF ICT, Department of Applied Mathematics,
+ * P.O. Box 124 Blindern,
+ * 0314 Oslo, Norway.
+ *
+ * This file is part of TTL.
+ *
+ * TTL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * TTL 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public
+ * License along with TTL. If not, see
+ * .
+ *
+ * In accordance with Section 7(b) of the GNU Affero General Public
+ * License, a covered work must retain the producer line in every data
+ * file that is created or manipulated using TTL.
+ *
+ * Other Usage
+ * You can be released from the requirements of the license by purchasing
+ * a commercial license. Buying such a license is mandatory as soon as you
+ * develop commercial activities involving the TTL library without
+ * disclosing the source code of your own applications.
+ *
+ * This file may be used in accordance with the terms contained in a
+ * written agreement between you and SINTEF ICT.
+ */
+
+#ifndef _TTL_H_
+#define _TTL_H_
+
+
+#include
+#include
+
+// Debugging
+#ifdef DEBUG_TTL
+ static void errorAndExit(char* message) {
+ cout << "\n!!! ERROR: " << message << " !!!\n" << endl;
+ exit(-1);
+ }
+#endif
+
+ using std::list;
+
+
+// Next on TOPOLOGY:
+// - get triangle strips
+// - weighted graph, algorithms using a weight (real) for each edge,
+// e.g. an "abstract length". Use for minimum spanning tree
+// or some arithmetics on weights?
+// - Circulators as defined in CGAL with more STL compliant code
+
+// - analyze in detail locateFace: e.g. detect 0-orbit in case of infinite loop
+// around a node etc.
+
+
+/** \brief Main interface to TTL
+*
+* This namespace contains the basic generic algorithms for the TTL,
+* the Triangulation Template Library.\n
+*
+* Examples of functionality are:
+* - Incremental Delaunay triangulation
+* - Constrained triangulation
+* - Insert/remove nodes and constrained edges
+* - Traversal operations
+* - Misc. queries for extracting information for visualisation systems etc.
+*
+* \par General requirements and assumptions:
+* - \e DartType and \e TraitsType should be implemented in accordance with the description
+* in \ref api.
+* - A \b "Requires:" section in the documentation of a function template
+* shows which functionality is required in \e TraitsType to
+* support that specific function.\n
+* Functionalty required in \e DartType is the same (almost) for all
+* function templates; see \ref api and the example referred to.
+* - When a reference to a \e dart object is passed to a function in TTL,
+* it is assumed that it is oriented \e counterclockwise (CCW) in a triangle
+* unless it is explicitly mentioned that it can also be \e clockwise (CW).
+* The same applies for a dart that is passed from a function in TTL to
+* the users TraitsType class (or struct).
+* - When an edge (represented with a dart) is swapped, it is assumed that darts
+* outside the quadrilateral where the edge is a diagonal are not affected by
+* the swap. Thus, \ref hed::TTLtraits::swapEdge "TraitsType::swapEdge"
+* must be implemented in accordance with this rule.
+*
+* \par Glossary:
+* - General terms are explained in \ref api.
+* - \e CCW - counterclockwise
+* - \e CW - clockwise
+* - \e 0_orbit, \e 1_orbit and \e 2_orbit: A sequence of darts around
+* a node, around an edge and in a triangle respectively;
+* see ttl::get_0_orbit_interior and ttl::get_0_orbit_boundary
+* - \e arc - In a triangulation an arc is equivalent with an edge
+*
+* \see
+* \ref ttl_util and \ref api
+*
+* \author
+* �yvind Hjelle, oyvindhj@ifi.uio.no
+*/
+
+
+namespace ttl {
+
+
+#ifndef DOXYGEN_SHOULD_SKIP_THIS
+ //------------------------------------------------------------------------------------------------
+ // ----------------------------------- Forward declarations -------------------------------------
+ //------------------------------------------------------------------------------------------------
+
+#if ((_MSC_VER > 0) && (_MSC_VER < 1300))
+#else
+
+ // Delaunay Triangulation
+ // ----------------------
+ template
+ bool insertNode(DartType& dart, PointType& point);
+
+ template
+ void removeRectangularBoundary(DartType& dart);
+
+ template
+ void removeNode(DartType& dart);
+
+ template
+ void removeBoundaryNode(DartType& dart);
+
+ template
+ void removeInteriorNode(DartType& dart);
+
+
+ // Topological and Geometric Queries
+ // ---------------------------------
+ template
+ bool locateFaceSimplest(const PointType& point, DartType& dart);
+
+ template
+ bool locateTriangle(const PointType& point, DartType& dart);
+
+ template
+ bool inTriangleSimplest(const PointType& point, const DartType& dart);
+
+ template
+ bool inTriangle(const PointType& point, const DartType& dart);
+
+ template
+ void getBoundary(const DartType& dart, DartListType& boundary);
+
+ template
+ bool isBoundaryEdge(const DartType& dart);
+
+ template
+ bool isBoundaryFace(const DartType& dart);
+
+ template
+ bool isBoundaryNode(const DartType& dart);
+
+ template
+ int getDegreeOfNode(const DartType& dart);
+
+ template
+ void get_0_orbit_interior(const DartType& dart, DartListType& orbit);
+
+ template
+ void get_0_orbit_boundary(const DartType& dart, DartListType& orbit);
+
+ template
+ bool same_0_orbit(const DartType& d1, const DartType& d2);
+
+ template
+ bool same_1_orbit(const DartType& d1, const DartType& d2);
+
+ template
+ bool same_2_orbit(const DartType& d1, const DartType& d2);
+
+ template
+ bool swappableEdge(const DartType& dart, bool allowDegeneracy = false);
+
+ template
+ void positionAtNextBoundaryEdge(DartType& dart);
+
+ template
+ bool convexBoundary(const DartType& dart);
+
+
+ // Utilities for Delaunay Triangulation
+ // ------------------------------------
+ template
+ void optimizeDelaunay(DartListType& elist);
+
+ template
+ void optimizeDelaunay(DartListType& elist, const typename DartListType::iterator end);
+
+ template
+ bool swapTestDelaunay(const DartType& dart, bool cycling_check = false);
+
+ template
+ void recSwapDelaunay(DartType& diagonal);
+
+ template
+ void swapEdgesAwayFromInteriorNode(DartType& dart, ListType& swapped_edges);
+
+ template
+ void swapEdgesAwayFromBoundaryNode(DartType& dart, ListType& swapped_edges);
+
+ template
+ void swapEdgeInList(const typename DartListType::iterator& it, DartListType& elist);
+
+
+ // Constrained Triangulation
+ // -------------------------
+ template
+ DartType insertConstraint(DartType& dstart, DartType& dend, bool optimize_delaunay);
+
+#endif
+
+#endif // DOXYGEN_SHOULD_SKIP_THIS
+
+
+ //------------------------------------------------------------------------------------------------
+ // ------------------------------- Delaunay Triangulation Group ---------------------------------
+ //------------------------------------------------------------------------------------------------
+
+ /** @name Delaunay Triangulation */
+ //@{
+
+ //------------------------------------------------------------------------------------------------
+ /** Inserts a new node in an existing Delaunay triangulation and
+ * swaps edges to obtain a new Delaunay triangulation.
+ * This is the basic function for incremental Delaunay triangulation.
+ * When starting from a set of points, an initial Delaunay triangulation
+ * can be created as two triangles forming a rectangle that contains
+ * all the points.
+ * After \c insertNode has been called repeatedly with all the points,
+ * ttl::removeRectangularBoundary can be called to remove triangles
+ * at the boundary of the triangulation so that the boundary
+ * form the convex hull of the points.
+ *
+ * Note that this incremetal scheme will run much faster if the points
+ * have been sorted lexicographically on \e x and \e y.
+ *
+ * \param dart
+ * An arbitrary CCW dart in the tringulation.\n
+ * Output: A CCW dart incident to the new node.
+ *
+ * \param point
+ * A point (node) to be inserted in the triangulation.
+ *
+ * \retval bool
+ * \c true if \e point was inserted; \c false if not.\n
+ * If \e point is outside the triangulation, or the input dart is not valid,
+ * \c false is returned.
+ *
+ * \require
+ * - \ref hed::TTLtraits::splitTriangle "TraitsType::splitTriangle" (DartType&, const PointType&)
+ *
+ * \using
+ * - ttl::locateTriangle
+ * - ttl::recSwapDelaunay
+ *
+ * \note
+ * - For efficiency reasons \e dart should be close to the insertion \e point.
+ *
+ * \see
+ * ttl::removeRectangularBoundary
+ */
+ template
+ bool insertNode(DartType& dart, PointType& point) {
+
+ bool found = ttl::locateTriangle(point, dart);
+ if (!found) {
+#ifdef DEBUG_TTL
+ cout << "ERROR: Triangulation::insertNode: NO triangle found. /n";
+#endif
+ return false;
+ }
+
+ // ??? can we hide the dart? this is not possible if one triangle only
+ TraitsType::splitTriangle(dart, point);
+
+ DartType d1 = dart;
+ d1.alpha2().alpha1().alpha2().alpha0().alpha1();
+
+ DartType d2 = dart;
+ d2.alpha1().alpha0().alpha1();
+
+ // Preserve a dart as output incident to the node and CCW
+ DartType d3 = dart;
+ d3.alpha2();
+ dart = d3; // and see below
+ //DartType dsav = d3;
+ d3.alpha0().alpha1();
+
+ //if (!TraitsType::fixedEdge(d1) && !ttl::isBoundaryEdge(d1)) {
+ if (!ttl::isBoundaryEdge(d1)) {
+ d1.alpha2();
+ recSwapDelaunay(d1);
+ }
+
+ //if (!TraitsType::fixedEdge(d2) && !ttl::isBoundaryEdge(d2)) {
+ if (!ttl::isBoundaryEdge(d2)) {
+ d2.alpha2();
+ recSwapDelaunay(d2);
+ }
+
+ // Preserve the incoming dart as output incident to the node and CCW
+ //d = dsav.alpha2();
+ dart.alpha2();
+ //if (!TraitsType::fixedEdge(d3) && !ttl::isBoundaryEdge(d3)) {
+ if (!ttl::isBoundaryEdge(d3)) {
+ d3.alpha2();
+ recSwapDelaunay(d3);
+ }
+
+ return true;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ // Private/Hidden function (might change later)
+ template
+ void insertNodes(ForwardIterator first, ForwardIterator last, DartType& dart) {
+
+ // Assumes that the dereferenced point objects are pointers.
+ // References to the point objects are then passed to TTL.
+
+ ForwardIterator it;
+ for (it = first; it != last; ++it) {
+ bool status = insertNode(dart, **it);
+ }
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Removes the rectangular boundary of a triangulation as a final step of an
+ * incremental Delaunay triangulation.
+ * The four nodes at the corners will be removed and the resulting triangulation
+ * will have a convex boundary and be Delaunay.
+ *
+ * \param dart
+ * A CCW dart at the boundary of the triangulation\n
+ * Output: A CCW dart at the new boundary
+ *
+ * \using
+ * - ttl::removeBoundaryNode
+ *
+ * \note
+ * - This function requires that the boundary of the triangulation is
+ * a rectangle with four nodes (one in each corner).
+ */
+ template
+ void removeRectangularBoundary(DartType& dart) {
+
+ DartType d_next = dart;
+ DartType d_iter;
+
+ for (int i = 0; i < 4; i++) {
+ d_iter = d_next;
+ d_next.alpha0();
+ ttl::positionAtNextBoundaryEdge(d_next);
+ ttl::removeBoundaryNode(d_iter);
+ }
+
+ dart = d_next; // Return a dart at the new boundary
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Removes the node associated with \e dart and
+ * updates the triangulation to be Delaunay.
+ *
+ * \using
+ * - ttl::removeBoundaryNode if \e dart represents a node at the boundary
+ * - ttl::removeInteriorNode if \e dart represents an interior node
+ *
+ * \note
+ * - The node cannot belong to a fixed (constrained) edge that is not
+ * swappable. (An endless loop is likely to occur in this case).
+ */
+ template
+ void removeNode(DartType& dart) {
+
+ if (ttl::isBoundaryNode(dart))
+ ttl::removeBoundaryNode(dart);
+ else
+ ttl::removeInteriorNode(dart);
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Removes the boundary node associated with \e dart and
+ * updates the triangulation to be Delaunay.
+ *
+ * \using
+ * - ttl::swapEdgesAwayFromBoundaryNode
+ * - ttl::optimizeDelaunay
+ *
+ * \require
+ * - \ref hed::TTLtraits::removeBoundaryTriangle "TraitsType::removeBoundaryTriangle" (Dart&)
+ */
+ template
+ void removeBoundaryNode(DartType& dart) {
+
+ // ... and update Delaunay
+ // - CCW dart must be given (for remove)
+ // - No dart is delivered back now (but this is possible if
+ // we assume that there is not only one triangle left in the triangulation.
+
+ // Position at boundary edge and CCW
+ if (!ttl::isBoundaryEdge(dart)) {
+ dart.alpha1(); // ensures that next function delivers back a CCW dart (if the given dart is CCW)
+ ttl::positionAtNextBoundaryEdge(dart);
+ }
+
+ list swapped_edges;
+ ttl::swapEdgesAwayFromBoundaryNode(dart, swapped_edges);
+
+ // Remove boundary triangles and remove the new boundary from the list
+ // of swapped edges, see below.
+ DartType d_iter = dart;
+ DartType dnext = dart;
+ bool bend = false;
+ while (bend == false) {
+ dnext.alpha1().alpha2();
+ if (ttl::isBoundaryEdge(dnext))
+ bend = true; // Stop when boundary
+
+ // Generic: Also remove the new boundary from the list of swapped edges
+ DartType n_bedge = d_iter;
+ n_bedge.alpha1().alpha0().alpha1().alpha2(); // new boundary edge
+
+ // ??? can we avoid find if we do this in swap away?
+ typename list::iterator it;
+ it = find(swapped_edges.begin(), swapped_edges.end(), n_bedge);
+
+ if (it != swapped_edges.end())
+ swapped_edges.erase(it);
+
+ // Remove the boundary triangle
+ TraitsType::removeBoundaryTriangle(d_iter);
+ d_iter = dnext;
+ }
+
+ // Optimize Delaunay
+ typedef list DartListType;
+ ttl::optimizeDelaunay(swapped_edges);
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Removes the interior node associated with \e dart and
+ * updates the triangulation to be Delaunay.
+ *
+ * \using
+ * - ttl::swapEdgesAwayFromInteriorNode
+ * - ttl::optimizeDelaunay
+ *
+ * \require
+ * - \ref hed::TTLtraits::reverse_splitTriangle "TraitsType::reverse_splitTriangle" (Dart&)
+ *
+ * \note
+ * - The node cannot belong to a fixed (constrained) edge that is not
+ * swappable. (An endless loop is likely to occur in this case).
+ */
+ template
+ void removeInteriorNode(DartType& dart) {
+
+ // ... and update to Delaunay.
+ // Must allow degeneracy temporarily, see comments in swap edges away
+ // Assumes:
+ // - revese_splitTriangle does not affect darts
+ // outside the resulting triangle.
+
+ // 1) Swaps edges away from the node until degree=3 (generic)
+ // 2) Removes the remaining 3 triangles and creates a new to fill the hole
+ // unsplitTriangle which is required
+ // 3) Runs LOP on the platelet to obtain a Delaunay triangulation
+ // (No dart is delivered as output)
+
+ // Assumes dart is counterclockwise
+
+ list swapped_edges;
+ ttl::swapEdgesAwayFromInteriorNode(dart, swapped_edges);
+
+ // The reverse operation of split triangle:
+ // Make one triangle of the three triangles at the node associated with dart
+ // TraitsType::
+ TraitsType::reverse_splitTriangle(dart);
+
+ // ???? Not generic yet if we are very strict:
+ // When calling unsplit triangle, darts at the three opposite sides may
+ // change!
+ // Should we hide them longer away??? This is possible since they cannot
+ // be boundary edges.
+ // ----> Or should we just require that they are not changed???
+
+ // Make the swapped-away edges Delaunay.
+ // Note the theoretical result: if there are no edges in the list,
+ // the triangulation is Delaunay already
+
+ ttl::optimizeDelaunay(swapped_edges);
+ }
+
+ //@} // End of Delaunay Triangulation Group
+
+
+ //------------------------------------------------------------------------------------------------
+ // -------------------------- Topological and Geometric Queries Group ---------------------------
+ //------------------------------------------------------------------------------------------------
+
+ /** @name Topological and Geometric Queries */
+ //@{
+
+ //------------------------------------------------------------------------------------------------
+ // Private/Hidden function (might change later)
+ template
+ bool isMemberOfFace(const TopologyElementType& topologyElement, const DartType& dart) {
+
+ // Check if the given topology element (node, edge or face) is a member of the face
+ // Assumes:
+ // - DartType::isMember(TopologyElementType)
+
+ DartType dart_iter = dart;
+ do {
+ if (dart_iter.isMember(topologyElement))
+ return true;
+ dart_iter.alpha0().alpha1();
+ } while (dart_iter != dart);
+
+ return false;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ // Private/Hidden function (might change later)
+ template
+ bool locateFaceWithNode(const NodeType& node, DartType& dart_iter) {
+ // Locate a face in the topology structure with the given node as a member
+ // Assumes:
+ // - TraitsType::orient2d(DartType, DartType, NodeType)
+ // - DartType::isMember(NodeType)
+ // - Note that if false is returned, the node might still be in the
+ // topology structure. Application programmer
+ // should check all if by hypothesis the node is in the topology structure;
+ // see doc. on locateTriangle.
+
+ bool status = locateFaceSimplest(node, dart_iter);
+ if (status == false)
+ return status;
+
+ // True was returned from locateFaceSimplest, but if the located triangle is
+ // degenerate and the node is on the extension of the edges,
+ // the node might still be inside. Check if node is a member and return false
+ // if not. (Still the node might be in the topology structure, see doc. above
+ // and in locateTriangle(const PointType& point, DartType& dart_iter)
+
+ return isMemberOfFace(node, dart_iter);
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Locates the face containing a given point.
+ * It is assumed that the tessellation (e.g. a triangulation) is \e regular in the sense that
+ * there are no holes, the boundary is convex and there are no degenerate faces.
+ *
+ * \param point
+ * A point to be located
+ *
+ * \param dart
+ * An arbitrary CCW dart in the triangulation\n
+ * Output: A CCW dart in the located face
+ *
+ * \retval bool
+ * \c true if a face is found; \c false if not.
+ *
+ * \require
+ * - \ref hed::TTLtraits::orient2d "TraitsType::orient2d" (DartType&, DartType&, PointType&)
+ *
+ * \note
+ * - If \c false is returned, \e point may still be inside a face if the tessellation is not
+ * \e regular as explained above.
+ *
+ * \see
+ * ttl::locateTriangle
+ */
+ template
+ bool locateFaceSimplest(const PointType& point, DartType& dart) {
+ // Not degenerate triangles if point is on the extension of the edges
+ // But inTriangle may be called in case of true (may update to inFace2)
+ // Convex boundary
+ // no holes
+ // convex faces (works for general convex faces)
+ // Not specialized for triangles, but ok?
+ //
+ // TraitsType::orint2d(PointType) is the half open half-plane defined
+ // by the dart:
+ // n1 = dart.node()
+ // n2 = dart.alpha0().node
+ // Only the following gives true:
+ // ((n2->x()-n1->x())*(point.y()-n1->y()) >= (point.x()-n1->x())*(n2->y()-n1->y()))
+
+ DartType dart_start;
+ dart_start = dart;
+ DartType dart_prev;
+
+ DartType d0;
+ for (;;) {
+ d0 = dart;
+ d0.alpha0();
+ if (TraitsType::orient2d(dart, d0, point) >= 0) {
+ dart.alpha0().alpha1();
+ if (dart == dart_start)
+ return true; // left to all edges in face
+ }
+ else {
+ dart_prev = dart;
+ dart.alpha2();
+ if (dart == dart_prev)
+ return false; // iteration to outside boundary
+
+ dart_start = dart;
+ dart_start.alpha0();
+
+ dart.alpha1(); // avoid twice on same edge and ccw in next
+ }
+ }
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Locates the triangle containing a given point.
+ * It is assumed that the triangulation is \e regular in the sense that there
+ * are no holes and the boundary is convex.
+ * This function deals with degeneracy to some extent, but round-off errors may still
+ * lead to a wrong result if triangles are degenerate.
+ *
+ * \param point
+ * A point to be located
+ *
+ * \param dart
+ * An arbitrary CCW dart in the triangulation\n
+ * Output: A CCW dart in the located triangle
+ *
+ * \retval bool
+ * \c true if a triangle is found; \c false if not.\n
+ * If \e point is outside the triangulation, in which case \c false is returned,
+ * then the edge associated with \e dart will be at the boundary of the triangulation.
+ *
+ * \using
+ * - ttl::locateFaceSimplest
+ * - ttl::inTriangle
+ */
+ template
+ bool locateTriangle(const PointType& point, DartType& dart) {
+ // The purpose is to have a fast and stable procedure that
+ // i) avoids concluding that a point is inside a triangle if it is not inside
+ // ii) avoids infinite loops
+
+ // Thus, if false is returned, the point might still be inside a triangle in
+ // the triangulation. But this will probably only occur in the following cases:
+ // i) There are holes in the triangulation which causes the procedure to stop.
+ // ii) The boundary of the triangulation is not convex.
+ // ii) There might be degenerate triangles interior to the triangulation, or on the
+ // the boundary, which in some cases might cause the procedure to stop there due
+ // to the logic of the algorithm.
+
+ // It is the application programmer's responsibility to check further if false is
+ // returned. For example, if by hypothesis the point is inside a triangle
+ // in the triangulation and and false is returned, then all triangles in the
+ // triangulation should be checked by the application. This can be done using
+ // the function:
+ // bool inTriangle(const PointType& point, const DartType& dart).
+
+
+ // Assumes:
+ // - crossProduct2d, scalarProduct2d etc., see functions called
+
+ bool status = locateFaceSimplest(point, dart);
+ if (status == false)
+ return status;
+
+ // There may be degeneracy, i.e., the point might be outside the triangle
+ // on the extension of the edges of a degenerate triangle.
+
+ // The next call returns true if inside a non-degenerate or a degenerate triangle,
+ // but false if the point coincides with the "supernode" in the case where all
+ // edges are degenerate.
+ return inTriangle(point, dart);
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Checks if \e point is inside the triangle associated with \e dart.
+ * A fast and simple function that does not deal with degeneracy.
+ *
+ * \param dart
+ * A CCW dart in the triangle
+ *
+ * \require
+ * - \ref hed::TTLtraits::orient2d "TraitsType::orient2d" (DartType&, DartType&, PointType&)
+ *
+ * \see
+ * ttl::inTriangle for a more robust function
+ */
+ template
+ bool inTriangleSimplest(const PointType& point, const DartType& dart) {
+
+ // Fast and simple: Do not deal with degenerate faces, i.e., if there is
+ // degeneracy, true will be returned if the point is on the extension of the
+ // edges of a degenerate triangle
+
+ DartType d_iter = dart;
+ DartType d0 = d_iter;
+ d0.alpha0();
+ if (!TraitsType::orient2d(d_iter, d0, point) >= 0)
+ return false;
+
+ d_iter.alpha0().alpha1();
+ d0 = d_iter;
+ d0.alpha0();
+ if (!TraitsType::orient2d(d_iter, d0, point) >= 0)
+ return false;
+
+ d_iter.alpha0().alpha1();
+ d0 = d_iter;
+ d0.alpha0();
+ if (!TraitsType::orient2d(d_iter, d0, point) >= 0)
+ return false;
+
+ return true;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Checks if \e point is inside the triangle associated with \e dart.
+ * This function deals with degeneracy to some extent, but round-off errors may still
+ * lead to wrong result if the triangle is degenerate.
+ *
+ * \param dart
+ * A CCW dart in the triangle
+ *
+ * \require
+ * - \ref hed::TTLtraits::crossProduct2d "TraitsType::crossProduct2d" (DartType&, PointType&)
+ * - \ref hed::TTLtraits::scalarProduct2d "TraitsType::scalarProduct2d" (DartType&, PointType&)
+ *
+ * \see
+ * ttl::inTriangleSimplest
+ */
+ template
+ bool inTriangle(const PointType& point, const DartType& dart) {
+
+ // SHOULD WE INCLUDE A STRATEGY WITH EDGE X e_1 ETC? TO GUARANTEE THAT
+ // ONLY ON ONE EDGE? BUT THIS DOES NOT SOLVE PROBLEMS WITH
+ // notInE1 && notInE1.neghbour ?
+
+
+ // Returns true if inside (but not necessarily strictly inside)
+ // Works for degenerate triangles, but not when all edges are degenerate,
+ // and the point coincides with all nodes;
+ // then false is always returned.
+
+ typedef typename TraitsType::real_type real_type;
+
+ DartType dart_iter = dart;
+
+ real_type cr1 = TraitsType::crossProduct2d(dart_iter, point);
+ if (cr1 < 0)
+ return false;
+
+ dart_iter.alpha0().alpha1();
+ real_type cr2 = TraitsType::crossProduct2d(dart_iter, point);
+
+ if (cr2 < 0)
+ return false;
+
+ dart_iter.alpha0().alpha1();
+ real_type cr3 = TraitsType::crossProduct2d(dart_iter, point);
+ if (cr3 < 0)
+ return false;
+
+ // All cross products are >= 0
+ // Check for degeneracy
+
+ if (cr1 != 0 || cr2 != 0 || cr3 != 0)
+ return true; // inside non-degenerate face
+
+ // All cross-products are zero, i.e. degenerate triangle, check if inside
+ // Strategy: d.scalarProduct2d >= 0 && alpha0(d).d.scalarProduct2d >= 0 for one of
+ // the edges. But if all edges are degenerate and the point is on (all) the nodes,
+ // then "false is returned".
+
+ DartType dart_tmp = dart_iter;
+ real_type sc1 = TraitsType::scalarProduct2d(dart_tmp,point);
+ real_type sc2 = TraitsType::scalarProduct2d(dart_tmp.alpha0(), point);
+
+ if (sc1 >= 0 && sc2 >= 0) {
+ // test for degenerate edge
+ if (sc1 != 0 || sc2 != 0)
+ return true; // interior to this edge or on a node (but see comment above)
+ }
+
+ dart_tmp = dart_iter.alpha0().alpha1();
+ sc1 = TraitsType::scalarProduct2d(dart_tmp,point);
+ sc2 = TraitsType::scalarProduct2d(dart_tmp.alpha0(),point);
+ if (sc1 >= 0 && sc2 >= 0) {
+ // test for degenerate edge
+ if (sc1 != 0 || sc2 != 0)
+ return true; // interior to this edge or on a node (but see comment above)
+ }
+
+ dart_tmp = dart_iter.alpha1();
+ sc1 = TraitsType::scalarProduct2d(dart_tmp,point);
+ sc2 = TraitsType::scalarProduct2d(dart_tmp.alpha0(),point);
+ if (sc1 >= 0 && sc2 >= 0) {
+ // test for degenerate edge
+ if (sc1 != 0 || sc2 != 0)
+ return true; // interior to this edge or on a node (but see comment above)
+ }
+
+ // Not on any of the edges of the degenerate triangle.
+ // The only possibility for the point to be "inside" is that all edges are degenerate
+ // and the point coincide with all nodes. So false is returned in this case.
+
+ return false;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ // Private/Hidden function (might change later)
+ template
+ void getAdjacentTriangles(const DartType& dart, DartType& t1, DartType& t2, DartType& t3) {
+
+ DartType dart_iter = dart;
+
+ // add first
+ if (dart_iter.alpha2() != dart) {
+ t1 = dart_iter;
+ dart_iter = dart;
+ }
+
+ // add second
+ dart_iter.alpha0();
+ dart_iter.alpha1();
+ DartType dart_prev = dart_iter;
+ if ((dart_iter.alpha2()) != dart_prev) {
+ t2 = dart_iter;
+ dart_iter = dart_prev;
+ }
+
+ // add third
+ dart_iter.alpha0();
+ dart_iter.alpha1();
+ dart_prev = dart_iter;
+ if ((dart_iter.alpha2()) != dart_prev)
+ t3 = dart_iter;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Gets the boundary as sequence of darts, where the edges associated with the darts are boundary
+ * edges, given a dart with an associating edge at the boundary of a topology structure.
+ * The first dart in the sequence will be the given one, and the others will have the same
+ * orientation (CCW or CW) as the first.
+ * Assumes that the given dart is at the boundary.
+ *
+ * \param dart
+ * A dart at the boundary (CCW or CW)
+ *
+ * \param boundary
+ * A sequence of darts, where the associated edges are the boundary edges
+ *
+ * \require
+ * - DartListType::push_back (DartType&)
+ */
+ template
+ void getBoundary(const DartType& dart, DartListType& boundary) {
+ // assumes the given dart is at the boundary (by edge)
+
+ DartType dart_iter(dart);
+ boundary.push_back(dart_iter); // Given dart as first element
+ dart_iter.alpha0();
+ positionAtNextBoundaryEdge(dart_iter);
+
+ while (dart_iter != dart) {
+ boundary.push_back(dart_iter);
+ dart_iter.alpha0();
+ positionAtNextBoundaryEdge(dart_iter);
+ }
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /*
+ // Asumes a fixed point (a boundary edge) is given
+ //
+ template
+ class boundary_1_Iterator { // i.e. "circulator"
+
+ DartType current_;
+ public:
+ boundaryEdgeIterator(const DartType& dart) {current_ = dart;}
+ DartType& operator * () const {return current_;}
+ void operator ++ () {current_.alpha0(); positionAtNextBoundaryEdge(current_);}
+ };
+ */
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Checks if the edge associated with \e dart is at
+ * the boundary of the triangulation.
+ *
+ * \par Implements:
+ * \code
+ * DartType dart_iter = dart;
+ * if (dart_iter.alpha2() == dart)
+ * return true;
+ * else
+ * return false;
+ * \endcode
+ */
+ template
+ bool isBoundaryEdge(const DartType& dart) {
+
+ DartType dart_iter = dart;
+ if (dart_iter.alpha2() == dart)
+ return true;
+ else
+ return false;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Checks if the face associated with \e dart is at
+ * the boundary of the triangulation.
+ */
+ template
+ bool isBoundaryFace(const DartType& dart) {
+
+ // Strategy: boundary if alpha2(d)=d
+
+ DartType dart_iter(dart);
+ DartType dart_prev;
+
+ do {
+ dart_prev = dart_iter;
+
+ if (dart_iter.alpha2() == dart_prev)
+ return true;
+ else
+ dart_iter = dart_prev; // back again
+
+ dart_iter.alpha0();
+ dart_iter.alpha1();
+
+ } while (dart_iter != dart);
+
+ return false;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Checks if the node associated with \e dart is at
+ * the boundary of the triangulation.
+ */
+ template
+ bool isBoundaryNode(const DartType& dart) {
+
+ // Strategy: boundary if alpha2(d)=d
+
+ DartType dart_iter(dart);
+ DartType dart_prev;
+
+ // If input dart is reached again, then internal node
+ // If alpha2(d)=d, then boundary
+
+ do {
+ dart_iter.alpha1();
+ dart_prev = dart_iter;
+ dart_iter.alpha2();
+
+ if (dart_iter == dart_prev)
+ return true;
+
+ } while (dart_iter != dart);
+
+ return false;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Returns the degree of the node associated with \e dart.
+ *
+ * \par Definition:
+ * The \e degree (or valency) of a node \e V in a triangulation,
+ * is defined as the number of edges incident with \e V, i.e.,
+ * the number of edges joining \e V with another node in the triangulation.
+ */
+ template
+ int getDegreeOfNode(const DartType& dart) {
+
+ DartType dart_iter(dart);
+ DartType dart_prev;
+
+ // If input dart is reached again, then interior node
+ // If alpha2(d)=d, then boundary
+
+ int degree = 0;
+ bool boundaryVisited = false;
+ do {
+ dart_iter.alpha1();
+ degree++;
+ dart_prev = dart_iter;
+
+ dart_iter.alpha2();
+
+ if (dart_iter == dart_prev) {
+ if (!boundaryVisited) {
+ boundaryVisited = true;
+ // boundary is reached first time, count in the reversed direction
+ degree++; // count the start since it is not done above
+ dart_iter = dart;
+ dart_iter.alpha2();
+ }
+ else
+ return degree;
+ }
+
+ } while (dart_iter != dart);
+
+ return degree;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ // Modification of getDegreeOfNode:
+ // Strategy, reverse the list and start in the other direction if the boundary
+ // is reached. NB. copying of darts but ok., or we could have collected pointers,
+ // but the memory management.
+
+ // NOTE: not symmetry if we choose to collect opposite edges
+ // now we collect darts with radiating edges
+
+ // Remember that we must also copy the node, but ok with push_back
+ // The size of the list will be the degree of the node
+
+ // No CW/CCW since topology only
+
+
+ // Each dart consists of an incident edge and an adjacent node.
+ // But note that this is only how we interpret the dart in this implementation.
+ // Given this list, how can we find the opposite edges:
+ // We can perform alpha1 on each, but for boundary nodes we will get one edge twice.
+ // But this is will always be the last dart!
+ // The darts in the list are in sequence and starts with the alpha0(dart)
+ // alpha0, alpha1 and alpha2
+
+ // Private/Hidden function
+ template
+ void getNeighborNodes(const DartType& dart, std::list& node_list, bool& boundary) {
+
+ DartType dart_iter(dart);
+
+ dart_iter.alpha0(); // position the dart at an opposite node
+
+ DartType dart_prev = dart_iter;
+
+ bool start_at_boundary = false;
+ dart_iter.alpha2();
+ if (dart_iter == dart_prev)
+ start_at_boundary = true;
+ else
+ dart_iter = dart_prev; // back again
+
+ DartType dart_start = dart_iter;
+
+ do {
+ node_list.push_back(dart_iter);
+ dart_iter.alpha1();
+ dart_iter.alpha0();
+ dart_iter.alpha1();
+ dart_prev = dart_iter;
+ dart_iter.alpha2();
+ if (dart_iter == dart_prev) {
+ // boundary reached
+ boundary = true;
+ if (start_at_boundary == true) {
+ // add the dart which now is positioned at the opposite boundary
+ node_list.push_back(dart_iter);
+ return;
+ }
+ else {
+ // call the function again such that we start at the boundary
+ // first clear the list and reposition to the initial node
+ dart_iter.alpha0();
+ node_list.clear();
+ getNeighborNodes(dart_iter, node_list, boundary);
+ return; // after one recursive step
+ }
+ }
+ } while (dart_iter != dart_start);
+
+ boundary = false;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Gets the 0-orbit around an interior node.
+ *
+ * \param dart
+ * A dart (CCW or CW) positioned at an \e interior node.
+ *
+ * \retval orbit
+ * Sequence of darts with one orbit for each arc. All the darts have the same
+ * orientation (CCW or CW) as \e dart, and \e dart is the first element
+ * in the sequence.
+ *
+ * \require
+ * - DartListType::push_back (DartType&)
+ *
+ * \see
+ * ttl::get_0_orbit_boundary
+ */
+ template
+ void get_0_orbit_interior(const DartType& dart, DartListType& orbit) {
+
+ DartType d_iter = dart;
+ orbit.push_back(d_iter);
+ d_iter.alpha1().alpha2();
+
+ while (d_iter != dart) {
+ orbit.push_back(d_iter);
+ d_iter.alpha1().alpha2();
+ }
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Gets the 0-orbit around a node at the boundary
+ *
+ * \param dart
+ * A dart (CCW or CW) positioned at a \e boundary \e node and at a \e boundary \e edge.
+ *
+ * \retval orbit
+ * Sequence of darts with one orbit for each arc. All the darts, \e exept \e the \e last one,
+ * have the same orientation (CCW or CW) as \e dart, and \e dart is the first element
+ * in the sequence.
+ *
+ * \require
+ * - DartListType::push_back (DartType&)
+ *
+ * \note
+ * - The last dart in the sequence have opposite orientation compared to the others!
+ *
+ * \see
+ * ttl::get_0_orbit_interior
+ */
+ template
+ void get_0_orbit_boundary(const DartType& dart, DartListType& orbit) {
+
+ DartType dart_prev;
+ DartType d_iter = dart;
+ do {
+ orbit.push_back(d_iter);
+ d_iter.alpha1();
+ dart_prev = d_iter;
+ d_iter.alpha2();
+ } while (d_iter != dart_prev);
+
+ orbit.push_back(d_iter); // the last one with opposite orientation
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Checks if the two darts belong to the same 0-orbit, i.e.,
+ * if they share a node.
+ * \e d1 and/or \e d2 can be CCW or CW.
+ *
+ * (This function also examines if the the node associated with
+ * \e d1 is at the boundary, which slows down the function (slightly).
+ * If it is known that the node associated with \e d1 is an interior
+ * node and a faster version is needed, the user should implement his/her
+ * own version.)
+ */
+ template
+ bool same_0_orbit(const DartType& d1, const DartType& d2) {
+
+ // Two copies of the same dart
+ DartType d_iter = d2;
+ DartType d_end = d2;
+
+ if (ttl::isBoundaryNode(d_iter)) {
+ // position at both boundary edges
+ ttl::positionAtNextBoundaryEdge(d_iter);
+ d_end.alpha1();
+ ttl::positionAtNextBoundaryEdge(d_end);
+ }
+
+ for (;;) {
+ if (d_iter == d1)
+ return true;
+ d_iter.alpha1();
+ if (d_iter == d1)
+ return true;
+ d_iter.alpha2();
+ if (d_iter == d_end)
+ break;
+ }
+
+ return false;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Checks if the two darts belong to the same 1-orbit, i.e.,
+ * if they share an edge.
+ * \e d1 and/or \e d2 can be CCW or CW.
+ */
+ template
+ bool same_1_orbit(const DartType& d1, const DartType& d2) {
+
+ DartType d_iter = d2;
+ // (Also works at the boundary)
+ if (d_iter == d1 || d_iter.alpha0() == d1 || d_iter.alpha2() == d1 || d_iter.alpha0() == d1)
+ return true;
+ return false;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Checks if the two darts belong to the same 2-orbit, i.e.,
+ * if they lie in the same triangle.
+ * \e d1 and/or \e d2 can be CCW or CW
+ */
+ template
+ bool same_2_orbit(const DartType& d1, const DartType& d2) {
+
+ DartType d_iter = d2;
+ if (d_iter == d1 || d_iter.alpha0() == d1 ||
+ d_iter.alpha1() == d1 || d_iter.alpha0() == d1 ||
+ d_iter.alpha1() == d1 || d_iter.alpha0() == d1)
+ return true;
+ return false;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ // Private/Hidden function
+ template
+ bool degenerateTriangle(const DartType& dart) {
+
+ // Check if triangle is degenerate
+ // Assumes CCW dart
+
+ DartType d1 = dart;
+ DartType d2 = d1;
+ d2.alpha1();
+ if (TraitsType::crossProduct2d(d1,d2) == 0)
+ return true;
+
+ return false;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Checks if the edge associated with \e dart is swappable, i.e., if the edge
+ * is a diagonal in a \e strictly convex (or convex) quadrilateral.
+ *
+ * \param allowDegeneracy
+ * If set to true, the function will also return true if the numerical calculations
+ * indicate that the quadrilateral is convex only, and not necessarily strictly
+ * convex.
+ *
+ * \require
+ * - \ref hed::TTLtraits::crossProduct2d "TraitsType::crossProduct2d" (Dart&, Dart&)
+ */
+ template
+ bool swappableEdge(const DartType& dart, bool allowDegeneracy) {
+
+ // How "safe" is it?
+
+ if (isBoundaryEdge(dart))
+ return false;
+
+ // "angles" are at the diagonal
+ DartType d1 = dart;
+ d1.alpha2().alpha1();
+ DartType d2 = dart;
+ d2.alpha1();
+ if (allowDegeneracy) {
+ if (TraitsType::crossProduct2d(d1,d2) < 0.0)
+ return false;
+ }
+ else {
+ if (TraitsType::crossProduct2d(d1,d2) <= 0.0)
+ return false;
+ }
+
+ // Opposite side (still angle at the diagonal)
+ d1 = dart;
+ d1.alpha0();
+ d2 = d1;
+ d1.alpha1();
+ d2.alpha2().alpha1();
+
+ if (allowDegeneracy) {
+ if (TraitsType::crossProduct2d(d1,d2) < 0.0)
+ return false;
+ }
+ else {
+ if (TraitsType::crossProduct2d(d1,d2) <= 0.0)
+ return false;
+ }
+ return true;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Given a \e dart, CCW or CW, positioned in a 0-orbit at the boundary of a tessellation.
+ * Position \e dart at a boundary edge in the same 0-orbit.\n
+ * If the given \e dart is CCW, \e dart is positioned at the left boundary edge
+ * and will be CW.\n
+ * If the given \e dart is CW, \e dart is positioned at the right boundary edge
+ * and will be CCW.
+ *
+ * \note
+ * - The given \e dart must have a source node at the boundary, otherwise an
+ * infinit loop occurs.
+ */
+ template
+ void positionAtNextBoundaryEdge(DartType& dart) {
+
+ DartType dart_prev;
+
+ // If alpha2(d)=d, then boundary
+
+ //old convention: dart.alpha0();
+ do {
+ dart.alpha1();
+ dart_prev = dart;
+ dart.alpha2();
+ } while (dart != dart_prev);
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Checks if the boundary of a triangulation is convex.
+ *
+ * \param dart
+ * A CCW dart at the boundary of the triangulation
+ *
+ * \require
+ * - \ref hed::TTLtraits::crossProduct2d "TraitsType::crossProduct2d" (const Dart&, const Dart&)
+ */
+ template
+ bool convexBoundary(const DartType& dart) {
+
+ list blist;
+ ttl::getBoundary(dart, blist);
+
+ int no;
+ no = (int)blist.size();
+ typename list::const_iterator bit = blist.begin();
+ DartType d1 = *bit;
+ ++bit;
+ DartType d2;
+ bool convex = true;
+ for (; bit != blist.end(); ++bit) {
+ d2 = *bit;
+ double crossProd = TraitsType::crossProduct2d(d1, d2);
+ if (crossProd < 0.0) {
+ //cout << "!!! Boundary is NOT convex: crossProd = " << crossProd << endl;
+ convex = false;
+ return convex;
+ }
+ d1 = d2;
+ }
+
+ // Check the last angle
+ d2 = *blist.begin();
+ double crossProd = TraitsType::crossProduct2d(d1, d2);
+ if (crossProd < 0.0) {
+ //cout << "!!! Boundary is NOT convex: crossProd = " << crossProd << endl;
+ convex = false;
+ }
+
+ //if (convex)
+ // cout << "\n---> Boundary is convex\n" << endl;
+ //cout << endl;
+ return convex;
+ }
+
+ //@} // End of Topological and Geometric Queries Group
+
+
+ //------------------------------------------------------------------------------------------------
+ // ------------------------ Utilities for Delaunay Triangulation Group --------------------------
+ //------------------------------------------------------------------------------------------------
+
+ /** @name Utilities for Delaunay Triangulation */
+ //@{
+
+ //------------------------------------------------------------------------------------------------
+ /** Optimizes the edges in the given sequence according to the
+ * \e Delaunay criterion, i.e., such that the edge will fullfill the
+ * \e circumcircle criterion (or equivalently the \e MaxMin
+ * angle criterion) with respect to the quadrilaterals where
+ * they are diagonals.
+ *
+ * \param elist
+ * The sequence of edges
+ *
+ * \require
+ * - \ref hed::TTLtraits::swapEdge "TraitsType::swapEdge" (DartType& \e dart)\n
+ * \b Note: Must be implemented such that \e dart is delivered back in a position as
+ * seen if it was glued to the edge when swapping (rotating) the edge CCW
+ *
+ * \using
+ * - ttl::swapTestDelaunay
+ */
+ template
+ void optimizeDelaunay(DartListType& elist) {
+ optimizeDelaunay(elist, elist.end());
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ template
+ void optimizeDelaunay(DartListType& elist, const typename DartListType::iterator end) {
+
+ // CCW darts
+ // Optimize here means Delaunay, but could be any criterion by
+ // requiring a "should swap" in the traits class, or give
+ // a function object?
+ // Assumes that elist has only one dart for each arc.
+ // Darts outside the quadrilateral are preserved
+
+ // For some data structures it is possible to preserve
+ // all darts when swapping. Thus a preserve_darts_when swapping
+ // ccould be given to indicate this and we would gain performance by avoiding
+ // find in list.
+
+ // Requires that swap retuns a dart in the "same position when rotated CCW"
+ // (A vector instead of a list may be better.)
+
+ // First check that elist is not empty
+ if (elist.empty())
+ return;
+
+ // Avoid cycling by more extensive circumcircle test
+ bool cycling_check = true;
+ bool optimal = false;
+ typename DartListType::iterator it;
+
+ typename DartListType::iterator end_opt = end;
+
+ // Hmm... The following code is trying to derefence an iterator that may
+ // be invalid. This may lead to debug error on Windows, so we comment out
+ // this code. Checking elist.empty() above will prevent some
+ // problems...
+ //
+ // last_opt is passed the end of the "active list"
+ //typename DartListType::iterator end_opt;
+ //if (*end != NULL)
+ // end_opt = end;
+ //else
+ // end_opt = elist.end();
+
+ while(!optimal) {
+ optimal = true;
+ for (it = elist.begin(); it != end_opt; ++it) {
+ if (ttl::swapTestDelaunay(*it, cycling_check)) {
+
+ // Preserve darts. Potential darts in the list are:
+ // - The current dart
+ // - the four CCW darts on the boundary of the quadrilateral
+ // (the current arc has only one dart)
+
+ ttl::swapEdgeInList(it, elist);
+
+ optimal = false;
+ } // end if should swap
+ } // end for
+ } // end pass
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Checks if the edge associated with \e dart should be swapped according
+ * to the \e Delaunay criterion, i.e., the \e circumcircle criterion (or
+ * equivalently the \e MaxMin angle criterion).
+ *
+ * \param cycling_check
+ * Must be set to \c true when used in connection with optimization algorithms,
+ * e.g., optimizeDelaunay. This will avoid cycling and infinite loops in nearly
+ * neutral cases.
+ *
+ * \require
+ * - \ref hed::TTLtraits::scalarProduct2d "TraitsType::scalarProduct2d" (DartType&, DartType&)
+ * - \ref hed::TTLtraits::crossProduct2d "TraitsType::crossProduct2d" (DartType&, DartType&)
+ */
+ template
+#if ((_MSC_VER > 0) && (_MSC_VER < 1300))//#ifdef _MSC_VER
+ bool swapTestDelaunay(const DartType& dart, bool cycling_check = false) {
+#else
+ bool swapTestDelaunay(const DartType& dart, bool cycling_check) {
+#endif
+
+ // The general strategy is taken from Cline & Renka. They claim that
+ // their algorithm insure numerical stability, but experiments show
+ // that this is not correct for neutral, or almost neutral cases.
+ // I have extended this strategy (without using tolerances) to avoid
+ // cycling and infinit loops when used in connection with LOP algorithms;
+ // see the comments below.
+
+ typedef typename TraitsType::real_type real_type;
+
+ if (isBoundaryEdge(dart))
+ return false;
+
+ DartType v11 = dart;
+ v11.alpha1().alpha0();
+ DartType v12 = v11;
+ v12.alpha1();
+
+ DartType v22 = dart;
+ v22.alpha2().alpha1().alpha0();
+ DartType v21 = v22;
+ v21.alpha1();
+
+ real_type cos1 = TraitsType::scalarProduct2d(v11,v12);
+ real_type cos2 = TraitsType::scalarProduct2d(v21,v22);
+
+ // "Angles" are opposite to the diagonal.
+ // The diagonals should be swapped iff (t1+t2) .gt. 180
+ // degrees. The following two tests insure numerical
+ // stability according to Cline & Renka. But experiments show
+ // that cycling may still happen; see the aditional test below.
+ if (cos1 >= 0 && cos2 >= 0) // both angles are grater or equual 90
+ return false;
+ if (cos1 < 0 && cos2 < 0) // both angles are less than 90
+ return true;
+
+ real_type sin1 = TraitsType::crossProduct2d(v11,v12);
+ real_type sin2 = TraitsType::crossProduct2d(v21,v22);
+ real_type sin12 = sin1*cos2 + cos1*sin2;
+ if (sin12 >= 0) // equality represents a neutral case
+ return false;
+
+ if (cycling_check) {
+ // situation so far is sin12 < 0. Test if this also
+ // happens for the swapped edge.
+
+ // The numerical calculations so far indicate that the edge is
+ // not Delaunay and should not be swapped. But experiments show that
+ // in neutral cases, or almost neutral cases, it may happen that
+ // the swapped edge may again be found to be not Delaunay and thus
+ // be swapped if we return true here. This may lead to cycling and
+ // an infinte loop when used, e.g., in connection with optimizeDelaunay.
+ //
+ // In an attempt to avoid this we test if the swapped edge will
+ // also be found to be not Delaunay by repeating the last test above
+ // for the swapped edge.
+ // We now rely on the general requirement for TraitsType::swapEdge which
+ // should deliver CCW dart back in "the same position"; see the general
+ // description. This will insure numerical stability as the next calculation
+ // is the same as if this function was called again with the swapped edge.
+ // Cycling is thus impossible provided that the initial tests above does
+ // not result in ambiguity (and they should probably not do so).
+
+ v11.alpha0();
+ v12.alpha0();
+ v21.alpha0();
+ v22.alpha0();
+ // as if the edge was swapped/rotated CCW
+ cos1 = TraitsType::scalarProduct2d(v22,v11);
+ cos2 = TraitsType::scalarProduct2d(v12,v21);
+ sin1 = TraitsType::crossProduct2d(v22,v11);
+ sin2 = TraitsType::crossProduct2d(v12,v21);
+ sin12 = sin1*cos2 + cos1*sin2;
+ if (sin12 < 0) {
+ // A neutral case, but the tests above lead to swapping
+ return false;
+ }
+ }
+
+ return true;
+ }
+
+
+ //-----------------------------------------------------------------------
+ //
+ // x
+ //" / \ "
+ // / | \ Darts:
+ //oe2 / | \ oe2 = oppEdge2
+ // x....|....x
+ // \ d| d/ d = diagonal (input and output)
+ // \ | /
+ // oe1 \ / oe1 = oppEdge1
+ // x
+ //
+ //-----------------------------------------------------------------------
+ /** Recursively swaps edges in the triangulation according to the \e Delaunay criterion.
+ *
+ * \param diagonal
+ * A CCW dart representing the edge where the recursion starts from.
+ *
+ * \require
+ * - \ref hed::TTLtraits::swapEdge "TraitsType::swapEdge" (DartType&)\n
+ * \b Note: Must be implemented such that the darts outside the quadrilateral
+ * are not affected by the swap.
+ *
+ * \using
+ * - Calls itself recursively
+ */
+ template
+ void recSwapDelaunay(DartType& diagonal) {
+
+ if (!ttl::swapTestDelaunay(diagonal))
+ // ??? ttl::swapTestDelaunay also checks if boundary, so this can be optimized
+ return;
+
+ // Get the other "edges" of the current triangle; see illustration above.
+ DartType oppEdge1 = diagonal;
+ oppEdge1.alpha1();
+ bool b1;
+ if (ttl::isBoundaryEdge(oppEdge1))
+ b1 = true;
+ else {
+ b1 = false;
+ oppEdge1.alpha2();
+ }
+
+
+ DartType oppEdge2 = diagonal;
+ oppEdge2.alpha0().alpha1().alpha0();
+ bool b2;
+ if (ttl::isBoundaryEdge(oppEdge2))
+ b2 = true;
+ else {
+ b2 = false;
+ oppEdge2.alpha2();
+ }
+
+ // Swap the given diagonal
+ TraitsType::swapEdge(diagonal);
+
+ if (!b1)
+ recSwapDelaunay(oppEdge1);
+ if (!b2)
+ recSwapDelaunay(oppEdge2);
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Swaps edges away from the (interior) node associated with
+ * \e dart such that that exactly three edges remain incident
+ * with the node.
+ * This function is used as a first step in ttl::removeInteriorNode
+ *
+ * \retval dart
+ * A CCW dart incident with the node
+ *
+ * \par Assumes:
+ * - The node associated with \e dart is interior to the
+ * triangulation.
+ *
+ * \require
+ * - \ref hed::TTLtraits::swapEdge "TraitsType::swapEdge" (DartType& \e dart)\n
+ * \b Note: Must be implemented such that \e dart is delivered back in a position as
+ * seen if it was glued to the edge when swapping (rotating) the edge CCW
+ *
+ * \note
+ * - A degenerate triangle may be left at the node.
+ * - The function is not unique as it depends on which dart
+ * at the node that is given as input.
+ *
+ * \see
+ * ttl::swapEdgesAwayFromBoundaryNode
+ */
+ template
+ void swapEdgesAwayFromInteriorNode(DartType& dart, ListType& swapped_edges) {
+
+ // Same iteration as in fixEdgesAtCorner, but not boundary
+ DartType dnext = dart;
+
+ // Allow degeneracy, otherwise we might end up with degree=4.
+ // For example, the reverse operation of inserting a point on an
+ // existing edge gives a situation where all edges are non-swappable.
+ // Ideally, degeneracy in this case should be along the actual node,
+ // but there is no strategy for this now.
+ // ??? An alternative here is to wait with degeneracy till we get an
+ // infinite loop with degree > 3.
+ bool allowDegeneracy = true;
+
+ int degree = ttl::getDegreeOfNode(dart);
+ DartType d_iter;
+ while (degree > 3) {
+ d_iter = dnext;
+ dnext.alpha1().alpha2();
+
+ if (ttl::swappableEdge(d_iter, allowDegeneracy)) {
+ TraitsType::swapEdge(d_iter); // swap the edge away
+ // Collect swapped edges in the list
+ // "Hide" the dart on the other side of the edge to avoid it being changed for
+ // other swaps
+ DartType swapped_edge = d_iter; // it was delivered back
+ swapped_edge.alpha2().alpha0(); // CCW (if not at boundary)
+ swapped_edges.push_back(swapped_edge);
+
+ degree--;
+ }
+ }
+ // Output, incident to the node
+ dart = dnext;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Swaps edges away from the (boundary) node associated with
+ * \e dart in such a way that when removing the edges that remain incident
+ * with the node, the boundary of the triangulation will be convex.
+ * This function is used as a first step in ttl::removeBoundaryNode
+ *
+ * \retval dart
+ * A CCW dart incident with the node
+ *
+ * \require
+ * - \ref hed::TTLtraits::swapEdge "TraitsType::swapEdge" (DartType& \e dart)\n
+ * \b Note: Must be implemented such that \e dart is delivered back in a position as
+ * seen if it was glued to the edge when swapping (rotating) the edge CCW
+ *
+ * \par Assumes:
+ * - The node associated with \e dart is at the boundary of the triangulation.
+ *
+ * \see
+ * ttl::swapEdgesAwayFromInteriorNode
+ */
+ template
+ void swapEdgesAwayFromBoundaryNode(DartType& dart, ListType& swapped_edges) {
+
+ // All darts that are swappable.
+ // To treat collinear nodes at an existing boundary, we must allow degeneracy
+ // when swapping to the boundary.
+ // dart is CCW and at the boundary.
+ // The 0-orbit runs CCW
+ // Deliver the dart back in the "same position".
+ // Assume for the swap in the traits class:
+ // - A dart on the swapped edge is delivered back in a position as
+ // seen if it was glued to the edge when swapping (rotating) the edge CCW
+
+ //int degree = ttl::getDegreeOfNode(dart);
+
+passes:
+
+ // Swap swappable edges that radiate from the node away
+ DartType d_iter = dart; // ???? can simply use dart
+ d_iter.alpha1().alpha2(); // first not at boundary
+ DartType d_next = d_iter;
+ bool bend = false;
+ bool swapped_next_to_boundary = false;
+ bool swapped_in_pass = false;
+
+ bool allowDegeneracy; // = true;
+ DartType tmp1, tmp2;
+
+ while (!bend) {
+
+ d_next.alpha1().alpha2();
+ if (ttl::isBoundaryEdge(d_next))
+ bend = true; // then it is CW since alpha2
+
+ // To allow removing among collinear nodes at the boundary,
+ // degenerate triangles must be allowed
+ // (they will be removed when used in connection with removeBoundaryNode)
+ tmp1 = d_iter; tmp1.alpha1();
+ tmp2 = d_iter; tmp2.alpha2().alpha1(); // don't bother with boundary (checked later)
+
+ if (ttl::isBoundaryEdge(tmp1) && ttl::isBoundaryEdge(tmp2))
+ allowDegeneracy = true;
+ else
+ allowDegeneracy = false;
+
+ if (ttl::swappableEdge(d_iter, allowDegeneracy)) {
+ TraitsType::swapEdge(d_iter);
+
+ // Collect swapped edges in the list
+ // "Hide" the dart on the other side of the edge to avoid it being changed for
+ // other swapps
+ DartType swapped_edge = d_iter; // it was delivered back
+ swapped_edge.alpha2().alpha0(); // CCW
+ swapped_edges.push_back(swapped_edge);
+
+ //degree--; // if degree is 2, or bend=true, we are done
+ swapped_in_pass = true;
+ if (bend)
+ swapped_next_to_boundary = true;
+ }
+ if (!bend)
+ d_iter = d_next;
+ }
+
+ // Deliver a dart as output in the same position as the incoming dart
+ if (swapped_next_to_boundary) {
+ // Assume that "swapping is CCW and dart is preserved in the same position
+ d_iter.alpha1().alpha0().alpha1(); // CW and see below
+ }
+ else {
+ d_iter.alpha1(); // CW and see below
+ }
+ ttl::positionAtNextBoundaryEdge(d_iter); // CCW
+
+ dart = d_iter; // for next pass or output
+
+ // If a dart was swapped in this iteration we must run it more
+ if (swapped_in_pass)
+ goto passes;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Swap the the edge associated with iterator \e it and update affected darts
+ * in \e elist accordingly.
+ * The darts affected by the swap are those in the same quadrilateral.
+ * Thus, if one want to preserve one or more of these darts on should
+ * keep them in \e elist.
+ */
+ template
+ void swapEdgeInList(const typename DartListType::iterator& it, DartListType& elist) {
+
+ typename DartListType::iterator it1, it2, it3, it4;
+ DartType dart(*it);
+
+ //typename TraitsType::DartType d1 = dart; d1.alpha2().alpha1();
+ //typename TraitsType::DartType d2 = d1; d2.alpha0().alpha1();
+ //typename TraitsType::DartType d3 = dart; d3.alpha0().alpha1();
+ //typename TraitsType::DartType d4 = d3; d4.alpha0().alpha1();
+ DartType d1 = dart; d1.alpha2().alpha1();
+ DartType d2 = d1; d2.alpha0().alpha1();
+ DartType d3 = dart; d3.alpha0().alpha1();
+ DartType d4 = d3; d4.alpha0().alpha1();
+
+ // Find pinters to the darts that may change.
+ // ??? Note, this is not very efficient since we must use find, which is O(N),
+ // four times.
+ // - Solution?: replace elist with a vector of pair (dart,number)
+ // and avoid find?
+ // - make a function for swapping generically?
+ // - sould we use another container type or,
+ // - erase them and reinsert?
+ // - or use two lists?
+ it1 = find(elist.begin(), elist.end(), d1);
+ it2 = find(elist.begin(), elist.end(), d2);
+ it3 = find(elist.begin(), elist.end(), d3);
+ it4 = find(elist.begin(), elist.end(), d4);
+
+ TraitsType::swapEdge(dart);
+ // Update the current dart which may have changed
+ *it = dart;
+
+ // Update darts that may have changed again (if they were present)
+ // Note that dart is delivered back after swapping
+ if (it1 != elist.end()) {
+ d1 = dart; d1.alpha1().alpha0();
+ *it1 = d1;
+ }
+ if (it2 != elist.end()) {
+ d2 = dart; d2.alpha2().alpha1();
+ *it2 = d2;
+ }
+ if (it3 != elist.end()) {
+ d3 = dart; d3.alpha2().alpha1().alpha0().alpha1();
+ *it3 = d3;
+ }
+ if (it4 != elist.end()) {
+ d4 = dart; d4.alpha0().alpha1();
+ *it4 = d4;
+ }
+ }
+
+ //@} // End of Utilities for Delaunay Triangulation Group
+
+}; // End of ttl namespace scope (but other files may also contain functions for ttl)
+
+
+ //------------------------------------------------------------------------------------------------
+ // ----------------------------- Constrained Triangulation Group --------------------------------
+ //------------------------------------------------------------------------------------------------
+
+ // Still namespace ttl
+
+#include
+
+#endif // _TTL_H_
diff --git a/include/ttl/ttl_constr.h b/include/ttl/ttl_constr.h
new file mode 100644
index 0000000000..67b46fcbf0
--- /dev/null
+++ b/include/ttl/ttl_constr.h
@@ -0,0 +1,632 @@
+/*
+ * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
+ * Applied Mathematics, Norway.
+ *
+ * Contact information: E-mail: tor.dokken@sintef.no
+ * SINTEF ICT, Department of Applied Mathematics,
+ * P.O. Box 124 Blindern,
+ * 0314 Oslo, Norway.
+ *
+ * This file is part of TTL.
+ *
+ * TTL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * TTL 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public
+ * License along with TTL. If not, see
+ * .
+ *
+ * In accordance with Section 7(b) of the GNU Affero General Public
+ * License, a covered work must retain the producer line in every data
+ * file that is created or manipulated using TTL.
+ *
+ * Other Usage
+ * You can be released from the requirements of the license by purchasing
+ * a commercial license. Buying such a license is mandatory as soon as you
+ * develop commercial activities involving the TTL library without
+ * disclosing the source code of your own applications.
+ *
+ * This file may be used in accordance with the terms contained in a
+ * written agreement between you and SINTEF ICT.
+ */
+
+#ifndef _TTL_CONSTR_H_
+#define _TTL_CONSTR_H_
+
+
+#include
+#include
+
+
+// Debugging
+#ifdef DEBUG_TTL_CONSTR_PLOT
+ #include
+ static ofstream ofile_constr("qweCons.dat");
+#endif
+
+
+//using namespace std;
+
+/** \brief Constrained Delaunay triangulation
+*
+* Basic generic algorithms in TTL for inserting a constrained edge between two existing nodes.\n
+*
+* See documentation for the namespace ttl for general requirements and assumptions.
+*
+* \author
+* Øyvind Hjelle, oyvindhj@ifi.uio.no
+*/
+
+namespace ttl_constr {
+
+ // ??? A constant used to evluate a numerical expression against a user spesified
+ // roundoff-zero number
+#ifdef DEBUG_TTL_CONSTR
+ static const double ROUNDOFFZERO = 0.0; // 0.1e-15;
+#endif
+
+
+ //------------------------------------------------------------------------------------------------
+ /* Checks if \e dart has start and end points in \e dstart and \e dend.
+ *
+ * \param dart
+ * The dart that should be controlled to see if it's the constraint
+ *
+ * \param dstart
+ * A CCW dart with the startnode of the constraint as the startnode
+ *
+ * \param dend
+ * A CCW dart with the endnode of the constraint as the startnode
+ *
+ * \retval bool
+ * A bool confirming that it's the constraint or not
+ *
+ * \using
+ * ttl::same_0_orbit
+ */
+ template
+ bool isTheConstraint(const DartType& dart, const DartType& dstart, const DartType& dend) {
+ DartType d0 = dart;
+ d0.alpha0(); // CW
+ if ((ttl::same_0_orbit(dstart, dart) && ttl::same_0_orbit(dend, d0)) ||
+ (ttl::same_0_orbit(dstart, d0) && ttl::same_0_orbit(dend, dart))) {
+ return true;
+ }
+ return false;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /* Checks if \e d1 and \e d2 are on the same side of the line between \e dstart and \e dend.
+ * (The start nodes of \e d1 and \e d2 represent an edge).
+ *
+ * \param dstart
+ * A CCW dart with the start node of the constraint as the source node of the dart.
+ *
+ * \param dend
+ * A CCW dart with the end node of the constraint as the source node of the dart.
+ *
+ * \param d1
+ * A CCW dart with the first node as the start node of the dart.
+ *
+ * \param d2
+ * A CCW dart with the other node as the start node of the dart.
+ *
+ * \using
+ * TraitsType::orient2d
+ */
+ template
+ bool crossesConstraint(DartType& dstart, DartType& dend, DartType& d1, DartType& d2) {
+
+ typename TraitsType::real_type orient_1 = TraitsType::orient2d(dstart,d1,dend);
+ typename TraitsType::real_type orient_2 = TraitsType::orient2d(dstart,d2,dend);
+ // ??? Should we refine this? e.g. find if (dstart,dend) (d1,d2) represent the same edge
+ if ((orient_1 <= 0 && orient_2 <= 0) || (orient_1 >= 0 && orient_2 >= 0))
+ return false;
+
+ return true;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /* Return the dart \e d making the smallest non-negative angle,
+ * as calculated with: orient2d(dstart, d.alpha0(), dend),
+ * at the 0-orbit of dstart.
+ * If (dstart,dend) is a CCW boundary edge \e d will be CW, otherwise CCW (since CCW in)
+ * at the 0-orbit of dstart.
+ *
+ * \par Assumes:
+ * - CCW dstart and dend, but returned dart can be CW at the boundary.
+ * - Boundary is convex?
+ *
+ * \param dstart
+ * A CCW dart dstart
+ *
+ * \param dend
+ * A CCW dart dend
+ *
+ * \retval DartType
+ * The dart \e d making the smallest positive (or == 0) angle
+ *
+ * \using
+ * ttl::isBoundaryNode
+ * ttl::positionAtNextBoundaryEdge
+ * TraitsType::orient2d
+ */
+ template
+ DartType getAtSmallestAngle(const DartType& dstart, const DartType& dend) {
+
+ // - Must boundary be convex???
+ // - Handle the case where the constraint is already present???
+ // - Handle dstart and/or dend at the boundary
+ // (dstart and dend may define a boundary edge)
+
+ DartType d_iter = dstart;
+ if (ttl::isBoundaryNode(d_iter)) {
+ d_iter.alpha1(); // CW
+ ttl::positionAtNextBoundaryEdge(d_iter); // CCW (was rotated CW to the boundary)
+ }
+
+ // assume convex boundary; see comments
+
+ DartType d0 = d_iter;
+ d0.alpha0();
+ bool ccw = true; // the rotation later
+ typename TraitsType::real_type o_iter = TraitsType::orient2d(d_iter, d0, dend);
+ if (o_iter == 0) { // collinear BUT can be on "back side"
+ d0.alpha1().alpha0(); // CW
+ if (TraitsType::orient2d(dstart, dend, d0) > 0)
+ return d_iter; //(=dstart) collinear
+ else {
+ // collinear on "back side"
+ d_iter.alpha1().alpha2(); // assume convex boundary
+ ccw = true;
+ }
+ }
+ else if (o_iter < 0) {
+ // Prepare for rotating CW and with d_iter CW
+ d_iter.alpha1();
+ ccw = false;
+ }
+
+ // Set first angle
+ d0 = d_iter; d0.alpha0();
+ o_iter = TraitsType::orient2d(dstart, d0, dend);
+
+ typename TraitsType::real_type o_next;
+
+ // Rotate towards the constraint CCW or CW.
+ // Here we assume that the boundary is convex.
+ DartType d_next = d_iter;
+ for (;;) {
+ d_next.alpha1(); // CW !!! (if ccw == true)
+ d0 = d_next; d0.alpha0();
+ o_next = TraitsType::orient2d(dstart, d0, dend);
+
+ if (ccw && o_next < 0) // and o_iter > 0
+ return d_iter;
+ else if (!ccw && o_next > 0)
+ return d_next; // CCW
+ else if (o_next == 0) {
+ if (ccw)
+ return d_next.alpha2(); // also ok if boundary
+ else
+ return d_next;
+ }
+
+ // prepare next
+ d_next.alpha2(); // CCW if ccw
+ d_iter = d_next; // also ok if boundary CCW if ccw == true
+ }
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /* This function finds all the edges in the triangulation crossing
+ * the spesified constraint and puts them in a list.
+ * In the case of collinearity, an attempt is made to detect this.
+ * The first collinear node between dstart and dend is then returned.
+ *
+ * Strategy:
+ * - Iterate such that \e d_iter is always strictly "below" the constraint
+ * as seen with \e dstart to the left and \e dend to the right.
+ * - Add CCW darts, whose edges intersect the constrait, to a list.
+ * These edges are found by the orient2d predicate:
+ * If two nodes of an edge are on opposite sides of the constraint,
+ * the edge between them intersect.
+ * - Must handle collinnear cases, i.e., if a node falls on the constraint,
+ * and possibly restarting collection of edges. Detecting collinearity
+ * heavily relies on the orient2d predicate which is provided by the
+ * traits class.
+ *
+ * Action:
+ * 1) Find cone/opening angle containing \e dstart and \e dend
+ * 2) Find first edge from the first 0-orbit that intersects
+ * 3) Check which of the two opposite that intersects
+ *
+ * 1)
+ * Rotate CCW and find the (only) case where \e d_iter and \e d_next satisfy:
+ * - orient2d(d_iter, d_iter.alpha0(), dend) > 0
+ * - orient2d(d_next, d_next.alpha0(), dend) < 0
+ *
+ * - check if we are done, i.e., if (d_next.alpha0() == my_dend)
+ * - Note also the situation if, e.g., the constraint is a boundary edge in which case
+ * \e my_dend wil be CW
+ *
+ * \param dstart
+ * A CCW dart with the startnode of the constraint as the startnode
+ *
+ * \param dend
+ * A CCW dart with the endnode of the constraint as the startnode
+ *
+ * \param elist
+ * A list where all the edges crossing the spesified constraint will be put
+ *
+ * \retval dartType
+ * Returns the next "collinear" starting node such that dend is returned when done.
+ */
+ template
+ DartType findCrossingEdges(const DartType& dstart, const DartType& dend, ListType& elist) {
+
+ const DartType my_start = getAtSmallestAngle(dstart, dend);
+ DartType my_end = getAtSmallestAngle(dend, dstart);
+
+ DartType d_iter = my_start;
+ if (d_iter.alpha0().alpha2() == my_end)
+ return d_iter; // The constraint is an existing edge and we are done
+
+ // Facts/status so far:
+ // - my_start and my_end are now both CCW and the constraint is not a boundary edge.
+ // - Further, the constraint is not one single existing edge, but it might be a collection
+ // of collinear edges in which case we return the current collinear edge
+ // and calling this function until all are collected.
+
+ my_end.alpha1(); // CW! // ??? this is probably ok for testing now?
+
+ d_iter = my_start;
+ d_iter.alpha0().alpha1(); // alpha0 is downwards or along the constraint
+
+ // Facts:
+ // - d_iter is guaranteed to intersect, but can be in start point.
+ // - d_iter.alpha0() is not at dend yet
+ typename TraitsType::real_type orient = TraitsType::orient2d(dstart, d_iter, dend);
+
+ // Use round-off error/tolerance or rely on the orient2d predicate ???
+ // Make a warning message if orient != exact 0
+ if (orient == 0)
+ return d_iter;
+
+#ifdef DEBUG_TTL_CONSTR
+ else if (fabs(orient) <= ROUNDOFFZERO) {
+ cout << "The darts are not exactly colinear, but |d1 x d2| <= " << ROUNDOFFZERO << endl;
+ return d_iter; // collinear, not done (and not collect in the list)
+ }
+#endif
+
+ // Collect intersecting edges
+ // --------------------------
+ elist.push_back(d_iter); // The first with interior intersection point
+
+ // Facts, status so far:
+ // - The first intersecting edge is now collected
+ // (- d_iter.alpha0() is still not at dend)
+
+ // d_iter should always be the edge that intersects and be below or on the constraint
+ // One of the two edges opposite to d_iter must intersect, or we have collinearity
+
+ // Note: Almost collinear cases can be handled on the
+ // application side with orient2d. We should probably
+ // return an int and the application will set it to zero
+ for(;;) {
+ // assume orient have been calc. and collinearity has been tested,
+ // above the first time and below later
+ d_iter.alpha2().alpha1(); // 2a same node
+
+ DartType d0 = d_iter;
+ d0.alpha0(); // CW
+ if (d0 == my_end)
+ return dend; // WE ARE DONE (but can we enter an endless loop???)
+
+ // d_iter or d_iter.alpha0().alpha1() must intersect
+ orient = TraitsType::orient2d(dstart, d0, dend);
+
+ if (orient == 0)
+ return d0.alpha1();
+
+#ifdef DEBUG_TTL_CONSTR
+ else if (fabs(orient) <= ROUNDOFFZERO) {
+ return d0.alpha1(); // CCW, collinear
+ }
+#endif
+
+ else if (orient > 0) { // orient > 0 and still below
+ // This one must intersect!
+ d_iter = d0.alpha1();
+ }
+ elist.push_back(d_iter);
+ }
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /* This function recives a constrained edge and a list of all the edges crossing a constraint.
+ * It then swaps the crossing edges away from the constraint. This is done according to a
+ * scheme suggested by Dyn, Goren & Rippa (slightly modified).
+ * The resulting triangulation is a constrained one, but not necessarily constrained Delaunay.
+ * In other to run optimization later to obtain a constrained Delaunay triangulation,
+ * the swapped edges are maintained in a list.
+ *
+ * Strategy :
+ * - Situation A: Run through the list and swap crossing edges away from the constraint.
+ * All the swapped edges are moved to the end of the list, and are "invisible" to this procedure.
+ * - Situation B: We may come in a situation where none of the crossing edges can be swapped away
+ * from the constraint.
+ * Then we follow the strategy of Dyn, Goren & Rippa and allow edges to be swapped,
+ * even if they are not swapped away from the constraint.
+ * These edges are NOT moved to the end of the list. They are later swapped to none-crossing
+ * edges when the locked situation is solved.
+ * - We keep on swapping edges in Situation B until we have iterated trough the list.
+ * We then resume Situation A.
+ * - This is done until the list is virtualy empty. The resulting \c elist has the constraint
+ * as the last element.
+ *
+ * \param dstart
+ * A CCW dart dstart
+ *
+ * \param dend
+ * A CCW dart dend
+ *
+ * \param elist
+ * A list containing all the edges crossing the spesified constraint
+ *
+ * \using
+ * ttl::swappableEdge
+ * ttl::swapEdgeInList
+ * ttl::crossesConstraint
+ * ttl::isTheConstraint
+ */
+ template
+ void transformToConstraint(DartType& dstart, DartType& dend, std::list& elist) {
+
+ typename list::iterator it, used;
+
+ // We may enter in a situation where dstart and dend are altered because of a swap.
+ // (The general rule is that darts inside the actual quadrilateral can be changed,
+ // but not those outside.)
+ // So, we need some look-ahead strategies for dstart and dend and change these
+ // after a swap if necessary.
+
+ int dartsInList = (int)elist.size();
+ if (dartsInList == 0)
+ return;
+
+ bool erase; // indicates if an edge is swapped away from the constraint such that it can be
+ // moved to the back of the list
+ bool locked = false;
+ do {
+ int noswap = 0;
+ it = elist.begin();
+
+ // counts how many edges that have been swapped per list-cycle
+ int counter = 1;
+ while(it != elist.end()) { // ??? change this test with counter > dartsInList
+ erase = false;
+ // Check if our virtual end of the list has been crossed. It breaks the
+ // while and starts all over again in the do-while loop
+ if (counter > dartsInList)
+ break;
+
+ if (ttl::swappableEdge(*it, true)) {
+ // Dyn & Goren & Rippa 's notation:
+ // The node assosiated with dart *it is denoted u_m. u_m has edges crossing the constraint
+ // named w_1, ... , w_r . The other node to the edge assosiated with dart *it is w_s.
+ // We want to swap from edge u_m<->w_s to edge w_{s-1}<->w_{s+1}.
+ DartType op1 = *it;
+ DartType op2 = op1;
+ op1.alpha1().alpha0(); //finds dart with node w_{s-1}
+ op2.alpha2().alpha1().alpha0(); // (CW) finds dart with node w_{s+1}
+ DartType tmp = *it; tmp.alpha0(); // Dart with assosiated node opposite to node of *it allong edge
+ // If there is a locked situation we swap, even if the result is crossing the constraint
+ // If there is a looked situation, but we do an ordinary swap, it should be treated as
+ // if we were not in a locked situation!!
+
+ // The flag swap_away indicates if the edge is swapped away from the constraint such that
+ // it does not cross the constraint.
+ bool swap_away = (crossesConstraint(dstart, dend, *it, tmp) &&
+ !crossesConstraint(dstart, dend, op1, op2));
+ if (swap_away || locked) {
+ // Do a look-ahead to see if dstart and/or dend are in the quadrilateral
+ // If so, we mark it with a flag to make sure we update them after the swap
+ // (they may have been changed during the swap according to the general rule!)
+ bool start = false;
+ bool end = false;
+
+ DartType d = *it;
+ if (d.alpha1().alpha0() == dstart)
+ start = true;
+ d = *it;
+ if (d.alpha2().alpha1().alpha0().alpha1() == dend)
+ end = true;
+
+ // This is the only place swapping is called when inserting a constraint
+ ttl::swapEdgeInList(it,elist);
+
+ // If we, during look-ahead, found that dstart and/or dend were in the quadrilateral,
+ // we update them.
+ if (end)
+ dend = *it;
+ if (start) {
+ dstart = *it;
+ dstart.alpha0().alpha2();
+ }
+
+ if (swap_away) { // !locked || //it should be sufficient with swap_away ???
+ noswap++;
+ erase = true;
+ }
+
+ if (isTheConstraint(*it, dstart, dend)) {
+ // Move the constraint to the end of the list
+ DartType the_constraint = *it;
+ elist.erase(it);
+ elist.push_back(the_constraint);
+ return;
+ } //endif
+ } //endif
+ } //endif "swappable edge"
+
+
+ // Move the edge to the end of the list if it was swapped away from the constraint
+ if (erase) {
+ used = it;
+ elist.push_back(*it);
+ ++it;
+ elist.erase(used);
+ --dartsInList;
+ }
+ else {
+ ++it;
+ ++counter;
+ }
+
+ } //end while
+
+ if (noswap == 0)
+ locked = true;
+
+ } while (dartsInList != 0);
+
+
+#ifdef DEBUG_TTL_CONSTR
+ // We will never enter here. (If elist is empty, we return above).
+ cout << "??????? ERROR 2, should never enter here ????????????????????????? SKIP ???? " << endl;
+ exit(-1);
+#endif
+
+ }
+
+}; // End of ttl_constr namespace scope
+
+
+namespace ttl { // (extension)
+
+ /** @name Constrained (Delaunay) Triangulation */
+ //@{
+
+ //------------------------------------------------------------------------------------------------
+ /** Inserts a constrained edge between two existing nodes in a triangulation.
+ * If the constraint falls on one or more existing nodes and this is detected by the
+ * predicate \c TraitsType::orient2d, which should return zero in this case, the
+ * constraint is split. Otherwise a degenerate triangle will be made along
+ * the constraint.
+ *
+ * \param dstart
+ * A CCW dart with the start node of the constraint as the source node
+ *
+ * \param dend
+ * A CCW dart with the end node of the constraint as the source node
+ *
+ * \param optimize_delaunay
+ * If set to \c true, the resulting triangulation will be
+ * a \e constrained \e Delaunay \e triangulation. If set to \c false, the resulting
+ * triangulation will not necessarily be of constrained Delaunay type.
+ *
+ * \retval DartType
+ * A dart representing the constrained edge.
+ *
+ * \require
+ * - \ref hed::TTLtraits::orient2d "TraitsType::orient2d" (DartType&, DartType&, PointType&)
+ * - \ref hed::TTLtraits::swapEdge "TraitsType::swapEdge" (DartType&)
+ *
+ * \using
+ * - ttl::optimizeDelaunay if \e optimize_delaunay is set to \c true
+ *
+ * \par Assumes:
+ * - The constrained edge must be inside the existing triangulation (and it cannot
+ * cross the boundary of the triangulation).
+ */
+ template
+ DartType insertConstraint(DartType& dstart, DartType& dend, bool optimize_delaunay) {
+
+ // Assumes:
+ // - It is the users responsibility to avoid crossing constraints
+ // - The constraint cannot cross the boundary, i.e., the boundary must be
+ // convex in the area of crossing edges.
+ // - dtart and dend are preserved (same node associated.)
+
+
+ // Find edges crossing the constraint and put them in elist.
+ // If findCrossingEdges reaches a Node lying on the constraint, this function
+ // calls itself recursively.
+
+ // RECURSION
+ list elist;
+ DartType next_start = ttl_constr::findCrossingEdges(dstart, dend, elist);
+
+ // If there are no crossing edges (elist is empty), we assume that the constraint
+ // is an existing edge.
+ // In this case, findCrossingEdges returns the constraint.
+ // Put the constraint in the list to fit with the procedures below
+ // (elist can also be empty in the case of invalid input data (the constraint is in
+ // a non-convex area) but this is the users responsibility.)
+
+ //by Thomas Sevaldrud if (elist.size() == 0)
+ //by Thomas Sevaldrud elist.push_back(next_start);
+
+ // findCrossingEdges stops if it finds a node lying on the constraint.
+ // A dart with this node as start node is returned
+ // We call insertConstraint recursivly until the received dart is dend
+ if (!ttl::same_0_orbit(next_start, dend)) {
+
+#ifdef DEBUG_TTL_CONSTR_PLOT
+ cout << "RECURSION due to collinearity along constraint" << endl;
+#endif
+
+ insertConstraint(next_start, dend, optimize_delaunay);
+ }
+
+ // Swap edges such that the constraint edge is present in the transformed triangulation.
+ if (elist.size() > 0) // by Thomas Sevaldrud
+ ttl_constr::transformToConstraint(dstart, next_start, elist);
+
+#ifdef DEBUG_TTL_CONSTR_PLOT
+ cout << "size of elist = " << elist.size() << endl;
+ if (elist.size() > 0) {
+ DartType the_constraint = elist.back();
+ ofile_constr << the_constraint.x() << " " << the_constraint.y() << " " << 0 << endl;
+ the_constraint.alpha0();
+ ofile_constr << the_constraint.x() << " " << the_constraint.y() << " " << 0 << endl << endl;
+ }
+#endif
+
+ // Optimize to constrained Delaunay triangulation if required.
+ typename list::iterator end_opt = elist.end();
+ if (optimize_delaunay) {
+
+ // Indicate that the constrained edge, which is the last element in the list,
+ // should not be swapped
+ --end_opt;
+ ttl::optimizeDelaunay(elist, end_opt);
+ }
+
+ if(elist.size() == 0) // by Thomas Sevaldrud
+ return next_start; // by Thomas Sevaldrud
+
+ // Return the constraint, which is still the last element in the list
+ end_opt = elist.end();
+ --end_opt;
+ return *end_opt;
+ }
+
+ //@} // End of Constrained Triangulation Group
+
+}; // End of ttl namespace scope (extension)
+
+#endif // _TTL_CONSTR_H_
diff --git a/include/ttl/ttl_util.h b/include/ttl/ttl_util.h
new file mode 100644
index 0000000000..cc311edb69
--- /dev/null
+++ b/include/ttl/ttl_util.h
@@ -0,0 +1,314 @@
+/*
+ * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
+ * Applied Mathematics, Norway.
+ *
+ * Contact information: E-mail: tor.dokken@sintef.no
+ * SINTEF ICT, Department of Applied Mathematics,
+ * P.O. Box 124 Blindern,
+ * 0314 Oslo, Norway.
+ *
+ * This file is part of TTL.
+ *
+ * TTL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * TTL 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 Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public
+ * License along with TTL. If not, see
+ * .
+ *
+ * In accordance with Section 7(b) of the GNU Affero General Public
+ * License, a covered work must retain the producer line in every data
+ * file that is created or manipulated using TTL.
+ *
+ * Other Usage
+ * You can be released from the requirements of the license by purchasing
+ * a commercial license. Buying such a license is mandatory as soon as you
+ * develop commercial activities involving the TTL library without
+ * disclosing the source code of your own applications.
+ *
+ * This file may be used in accordance with the terms contained in a
+ * written agreement between you and SINTEF ICT.
+ */
+
+#ifndef _TTL_UTIL_H_
+#define _TTL_UTIL_H_
+
+
+#include
+#include
+
+
+#ifdef _MSC_VER
+# if _MSC_VER < 1300
+# include
+# endif
+#endif
+
+
+//using namespace std;
+
+
+/** \brief Utilities
+*
+* This name space contains utility functions for TTL.\n
+*
+* Point and vector algebra such as scalar product and cross product
+* between vectors are implemented here.
+* These functions are required by functions in the \ref ttl namespace,
+* where they are assumed to be present in the \ref hed::TTLtraits "TTLtraits" class.
+* Thus, the user can call these functions from the traits class.
+* For efficiency reasons, the user may consider implementing these
+* functions in the the API directly on the actual data structure;
+* see \ref api.
+*
+* \note
+* - Cross product between vectors in the xy-plane delivers a scalar,
+* which is the z-component of the actual cross product
+* (the x and y components are both zero).
+*
+* \see
+* ttl and \ref api
+*
+* \author
+* Øyvind Hjelle, oyvindhj@ifi.uio.no
+*/
+
+
+namespace ttl_util {
+
+
+ //------------------------------------------------------------------------------------------------
+ // ------------------------------ Computational Geometry Group ----------------------------------
+ //------------------------------------------------------------------------------------------------
+
+ /** @name Computational geometry */
+ //@{
+
+ //------------------------------------------------------------------------------------------------
+ /** Scalar product between two 2D vectors.
+ *
+ * \par Returns:
+ * \code
+ * dx1*dx2 + dy1*dy2
+ * \endcode
+ */
+ template
+ real_type scalarProduct2d(real_type dx1, real_type dy1, real_type dx2, real_type dy2) {
+ return dx1*dx2 + dy1*dy2;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Cross product between two 2D vectors. (The z-component of the actual cross product.)
+ *
+ * \par Returns:
+ * \code
+ * dx1*dy2 - dy1*dx2
+ * \endcode
+ */
+ template
+ real_type crossProduct2d(real_type dx1, real_type dy1, real_type dx2, real_type dy2) {
+ return dx1*dy2 - dy1*dx2;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /** Returns a positive value if the 2D nodes/points \e pa, \e pb, and
+ * \e pc occur in counterclockwise order; a negative value if they occur
+ * in clockwise order; and zero if they are collinear.
+ *
+ * \note
+ * - This is a finite arithmetic fast version. It can be made more robust using
+ * exact arithmetic schemes by Jonathan Richard Shewchuk. See
+ * http://www-2.cs.cmu.edu/~quake/robust.html
+ */
+ template
+ real_type orient2dfast(real_type pa[2], real_type pb[2], real_type pc[2]) {
+ real_type acx = pa[0] - pc[0];
+ real_type bcx = pb[0] - pc[0];
+ real_type acy = pa[1] - pc[1];
+ real_type bcy = pb[1] - pc[1];
+ return acx * bcy - acy * bcx;
+ }
+
+
+ //------------------------------------------------------------------------------------------------
+ /* Scalar product between 2D vectors represented as darts.
+ *
+ * \par Requires:
+ * - real_type DartType::x()
+ * - real_type DartType::y()
+ */
+ /*
+ template
+ typename TTLtraits::real_type scalarProduct2d(const DartType& d1, const DartType& d2) {
+
+ DartType d10 = d1;
+ d10.alpha0();
+
+ DartType d20 = d2;
+ d20.alpha0();
+
+ return scalarProduct2d(d10.x() - d1.x(), d10.y() - d1.y(), d20.x() - d2.x(), d20.y() - d2.y());
+ }
+ */
+
+
+ //------------------------------------------------------------------------------------------------
+ /* Scalar product between 2D vectors.
+ * The first vector is represented by the given dart, and the second vector has
+ * direction from the node of the given dart - and to the given point.
+ *
+ * \par Requires:
+ * - real_type DartType::x(), real_type DartType::y()
+ * - real_type PointType2d::x(), real_type PointType2d::y()
+ */
+ /*
+ template
+ typename TTLtraits::real_type scalarProduct2d(const typename TTLtraits::DartType& d,
+ const typename TTLtraits::PointType2d& p) {
+ typename TTLtraits::DartType d0 = d;
+ d0.alpha0();
+
+ return scalarProduct2d(d0.x() - d.x(), d0.y() - d.y(), p.x() - d.x(), p.y() - d.y());
+ }
+ */
+
+
+ //------------------------------------------------------------------------------------------------
+ /* Cross product between 2D vectors represented as darts.
+ *
+ * \par Requires:
+ * - real_type DartType::x(), real_type DartType::y()
+ */
+ /*
+ template
+ typename TTLtraits::real_type crossProduct2d(const typename TTLtraits::DartType& d1,
+ const typename TTLtraits::DartType& d2) {
+
+ TTLtraits::DartType d10 = d1;
+ d10.alpha0();
+
+ TTLtraits::DartType d20 = d2;
+ d20.alpha0();
+
+ return crossProduct2d(d10.x() - d1.x(), d10.y() - d1.y(), d20.x() - d2.x(), d20.y() - d2.y());
+ }
+ */
+
+
+ //------------------------------------------------------------------------------------------------
+ /* Cross product between 2D vectors.
+ * The first vector is represented by the given dart, and the second vector has
+ * direction from the node associated with given dart - and to the given point.
+ *
+ * \par Requires:
+ * - real_type DartType::x()
+ * - real_type DartType::y()
+ */
+ /*
+ template
+ typename TTLtraits::real_type crossProduct2d(const typename TTLtraits::DartType& d,
+ const typename TTLtraits::PointType2d& p) {
+
+ TTLtraits::DartType d0 = d;
+ d0.alpha0();
+
+ return crossProduct2d(d0.x() - d.x(), d0.y() - d.y(), p.x() - d.x(), p.y() - d.y());
+ }
+ */
+ // Geometric predicates; see more robust schemes by Jonathan Richard Shewchuk at
+ // http://www.cs.cmu.edu/~quake/robust.html
+
+
+ //------------------------------------------------------------------------------------------------
+ /* Return a positive value if the 2d nodes/points \e d, \e d.alpha0(), and
+ * \e p occur in counterclockwise order; a negative value if they occur
+ * in clockwise order; and zero if they are collinear. The
+ * result is also a rough approximation of twice the signed
+ * area of the triangle defined by the three points.
+ *
+ * \par Requires:
+ * - DartType::x(), DartType::y(),
+ * - PointType2d::x(), PointType2d::y()
+ */
+ /*
+ template
+ typename TTLtraits::real_type orient2dfast(const DartType& n1, const DartType& n2,
+ const PointType2d& p) {
+ return ((n2.x()-n1.x())*(p.y()-n1.y()) - (p.x()-n1.x())*(n2.y()-n1.y()));
+ }
+ */
+
+ //@} // End of Computational geometry
+
+
+ //------------------------------------------------------------------------------------------------
+ // ---------------------------- Utilities Involving Points Group --------------------------------
+ //------------------------------------------------------------------------------------------------
+
+ /** @name Utilities involving points */
+ //@{
+
+ //------------------------------------------------------------------------------------------------
+ /** Creates random data on the unit square.
+ *
+ * \param noPoints
+ * Number of random points to be generated
+ *
+ * \param seed
+ * Initial value for pseudorandom number generator
+ *
+ * \require
+ * - Constructor \c PointType::PointType(double x, double y).\n
+ * For example, one can use \c pair.
+ *
+ * \note
+ * - To deduce template argument for PointType the function must be
+ * called with the syntax: \c createRandomData(...) where \c MyPoint
+ * is the actual point type.
+ */
+ template
+ std::vector* createRandomData(int noPoints, int seed=1) {
+
+#ifdef _MSC_VER
+ srand(seed);
+#else
+ srand48((long int)seed);
+#endif
+
+ double x, y;
+ std::vector* points = new std::vector(noPoints);
+ typename std::vector::iterator it;
+ for (it = points->begin(); it != points->end(); ++it) {
+
+#ifdef _MSC_VER
+ int random = rand();
+ x = ((double)random/(double)RAND_MAX);
+ random = rand();
+ y = ((double)random/(double)RAND_MAX);
+ *it = new PointType(x,y);
+#else
+ double random = drand48();
+ x = random;
+ random = drand48();
+ y = random;
+ *it = new PointType(x,y);
+#endif
+
+ }
+ return points;
+ }
+
+ //@} // End of Utilities involving points
+
+}; // End of ttl_util namespace scope
+
+#endif // _TTL_UTIL_H_
diff --git a/include/wxBasePcbFrame.h b/include/wxBasePcbFrame.h
index f89a4ed383..1eaae3598f 100644
--- a/include/wxBasePcbFrame.h
+++ b/include/wxBasePcbFrame.h
@@ -183,7 +183,7 @@ public:
* BOARD.
* @param aBoard The BOARD to put into the frame.
*/
- void SetBoard( BOARD* aBoard );
+ virtual void SetBoard( BOARD* aBoard );
BOARD* GetBoard() const
{
@@ -191,8 +191,6 @@ public:
return m_Pcb;
}
- void ViewReloadBoard( const BOARD* aBoard ) const;
-
/**
* Function SetFootprintLibTable
* set the footprint library table to \a aFootprintLibTable.
@@ -717,8 +715,6 @@ public:
void OnUpdateSelectGrid( wxUpdateUIEvent& aEvent );
void OnUpdateSelectZoom( wxUpdateUIEvent& aEvent );
- virtual void UseGalCanvas( bool aEnable );
-
DECLARE_EVENT_TABLE()
};
diff --git a/include/wxPcbStruct.h b/include/wxPcbStruct.h
index 17be2ccdfd..c6b67b8086 100644
--- a/include/wxPcbStruct.h
+++ b/include/wxPcbStruct.h
@@ -225,7 +225,7 @@ public:
const wxPoint& pos, const wxSize& size,
long style = KICAD_DEFAULT_DRAWFRAME_STYLE );
- ~PCB_EDIT_FRAME();
+ virtual ~PCB_EDIT_FRAME();
void OnQuit( wxCommandEvent& event );
@@ -603,6 +603,13 @@ public:
*/
void Show3D_Frame( wxCommandEvent& event );
+ /**
+ * Function UseGalCanvas
+ * Enables/disables GAL canvas.
+ * @param aEnable determines if GAL should be active or not.
+ */
+ void UseGalCanvas( bool aEnable );
+
/**
* Function ChangeCanvas
* switches currently used canvas (default / Cairo / OpenGL).
@@ -888,6 +895,15 @@ public:
*/
bool Clear_Pcb( bool aQuery );
+ /// @copydoc PCB_BASE_FRAME::SetBoard()
+ void SetBoard( BOARD* aBoard );
+
+ /**
+ * Function ViewReloadBoard
+ * adds all items from the current board to the VIEW, so they can be displayed by GAL.
+ */
+ void ViewReloadBoard( const BOARD* aBoard ) const;
+
// Drc control
/* function GetDrcController
diff --git a/pcbnew/CMakeLists.txt b/pcbnew/CMakeLists.txt
index f243ce065f..b5f5e621fe 100644
--- a/pcbnew/CMakeLists.txt
+++ b/pcbnew/CMakeLists.txt
@@ -211,6 +211,8 @@ set( PCBNEW_CLASS_SRCS
print_board_functions.cpp
printout_controler.cpp
ratsnest.cpp
+ ratsnest_data.cpp
+ ratsnest_viewitem.cpp
# specctra.cpp #moved in pcbcommon lib
# specctra_export.cpp
# specctra_keywords.cpp
@@ -370,8 +372,8 @@ if( KICAD_SCRIPTING_MODULES )
polygon
bitmaps
gal
- ${GLEW_LIBRARIES}
- ${CAIRO_LIBRARIES}
+ ${GLEW_LIBRARIES}
+ ${CAIRO_LIBRARIES}
${wxWidgets_LIBRARIES}
${OPENGL_LIBRARIES}
${GDI_PLUS_LIBRARIES}
diff --git a/pcbnew/basepcbframe.cpp b/pcbnew/basepcbframe.cpp
index 89c5a8866e..7b3f66df58 100644
--- a/pcbnew/basepcbframe.cpp
+++ b/pcbnew/basepcbframe.cpp
@@ -54,6 +54,8 @@
#include
#include
#include
+#include
+#include
#include
#include
@@ -79,6 +81,7 @@ const LAYER_NUM PCB_BASE_FRAME::GAL_LAYER_ORDER[] =
ITEM_GAL_LAYER( MOD_TEXT_FR_VISIBLE ),
ITEM_GAL_LAYER( MOD_REFERENCES_VISIBLE), ITEM_GAL_LAYER( MOD_VALUES_VISIBLE ),
+ ITEM_GAL_LAYER( RATSNEST_VISIBLE ),
ITEM_GAL_LAYER( VIAS_HOLES_VISIBLE ), ITEM_GAL_LAYER( PADS_HOLES_VISIBLE ),
ITEM_GAL_LAYER( VIAS_VISIBLE ), ITEM_GAL_LAYER( PADS_VISIBLE ),
@@ -171,102 +174,6 @@ void PCB_BASE_FRAME::SetBoard( BOARD* aBoard )
{
delete m_Pcb;
m_Pcb = aBoard;
-
- if( m_galCanvas )
- {
- KIGFX::VIEW* view = m_galCanvas->GetView();
-
- try
- {
- ViewReloadBoard( m_Pcb );
- }
- catch( const std::exception& ex )
- {
- DBG(printf( "ViewReloadBoard: exception: %s\n", ex.what() );)
- }
-
- // update the tool manager with the new board and its view.
- if( m_toolManager )
- m_toolManager->SetEnvironment( m_Pcb, view, m_galCanvas->GetViewControls(), this );
- }
-}
-
-
-void PCB_BASE_FRAME::ViewReloadBoard( const BOARD* aBoard ) const
-{
- KIGFX::VIEW* view = m_galCanvas->GetView();
- view->Clear();
-
- // All of PCB drawing elements should be added to the VIEW
- // in order to be displayed
-
- // Load zones
- for( int i = 0; i < aBoard->GetAreaCount(); ++i )
- {
- view->Add( (KIGFX::VIEW_ITEM*) ( aBoard->GetArea( i ) ) );
- }
-
- // Load drawings
- for( BOARD_ITEM* drawing = aBoard->m_Drawings; drawing; drawing = drawing->Next() )
- {
- view->Add( drawing );
- }
-
- // Load tracks
- for( TRACK* track = aBoard->m_Track; track; track = track->Next() )
- {
- view->Add( track );
- }
-
- // Load modules and its additional elements
- for( MODULE* module = aBoard->m_Modules; module; module = module->Next() )
- {
- // Load module's pads
- for( D_PAD* pad = module->Pads().GetFirst(); pad; pad = pad->Next() )
- {
- view->Add( pad );
- }
-
- // Load module's drawing (mostly silkscreen)
- for( BOARD_ITEM* drawing = module->GraphicalItems().GetFirst(); drawing;
- drawing = drawing->Next() )
- {
- view->Add( drawing );
- }
-
- // Load module's texts (name and value)
- view->Add( &module->Reference() );
- view->Add( &module->Value() );
-
- // Add the module itself
- view->Add( module );
- }
-
- // Segzones (equivalent of ZONE_CONTAINER for legacy boards)
- for( SEGZONE* zone = aBoard->m_Zone; zone; zone = zone->Next() )
- {
- view->Add( zone );
- }
-
- // Add an entry for the worksheet layout
- KIGFX::WORKSHEET_VIEWITEM* worksheet = new KIGFX::WORKSHEET_VIEWITEM(
- std::string( aBoard->GetFileName().mb_str() ),
- std::string( GetScreenDesc().mb_str() ),
- &GetPageSettings(), &GetTitleBlock() );
- BASE_SCREEN* screen = GetScreen();
- if( screen != NULL )
- {
- worksheet->SetSheetNumber( GetScreen()->m_ScreenNumber );
- worksheet->SetSheetCount( GetScreen()->m_NumberOfScreens );
- }
-
- view->Add( worksheet );
-
- view->SetPanBoundary( worksheet->ViewBBox() );
- view->RecacheAllItems( true );
-
- if( m_galCanvasActive )
- m_galCanvas->Refresh();
}
@@ -604,17 +511,6 @@ void PCB_BASE_FRAME::OnUpdateSelectZoom( wxUpdateUIEvent& aEvent )
}
-void PCB_BASE_FRAME::UseGalCanvas( bool aEnable )
-{
- EDA_DRAW_FRAME::UseGalCanvas( aEnable );
-
- m_toolManager->SetEnvironment( m_Pcb, m_galCanvas->GetView(),
- m_galCanvas->GetViewControls(), this );
-
- ViewReloadBoard( m_Pcb );
-}
-
-
void PCB_BASE_FRAME::ProcessItemSelection( wxCommandEvent& aEvent )
{
int id = aEvent.GetId();
@@ -908,6 +804,7 @@ void PCB_BASE_FRAME::LoadSettings()
view->SetRequired( SOLDERMASK_N_BACK, ITEM_GAL_LAYER( PAD_BK_VISIBLE ) );
view->SetLayerTarget( ITEM_GAL_LAYER( GP_OVERLAY ), KIGFX::TARGET_OVERLAY );
+ view->SetLayerTarget( ITEM_GAL_LAYER( RATSNEST_VISIBLE ), KIGFX::TARGET_OVERLAY );
// Apply layer coloring scheme & display options
if( view->GetPainter() )
diff --git a/pcbnew/class_board.cpp b/pcbnew/class_board.cpp
index ca97243fc2..22b082fc9d 100644
--- a/pcbnew/class_board.cpp
+++ b/pcbnew/class_board.cpp
@@ -42,6 +42,7 @@
#include
#include
#include
+#include
#include
#include
@@ -102,11 +103,15 @@ BOARD::BOARD() :
m_NetClasses.GetDefault()->SetParams();
SetCurrentNetClass( m_NetClasses.GetDefault()->GetName() );
+
+ m_ratsnest = new RN_DATA( this );
}
BOARD::~BOARD()
{
+ delete m_ratsnest;
+
while( m_ZoneDescriptorList.size() )
{
ZONE_CONTAINER* area_to_remove = m_ZoneDescriptorList[0];
diff --git a/pcbnew/class_board.h b/pcbnew/class_board.h
index 4cd30d1f6f..cd20413baa 100644
--- a/pcbnew/class_board.h
+++ b/pcbnew/class_board.h
@@ -56,6 +56,7 @@ class MARKER_PCB;
class MSG_PANEL_ITEM;
class NETLIST;
class REPORTER;
+class RN_DATA;
// non-owning container of item candidates when searching for items on the same track.
@@ -225,6 +226,7 @@ private:
EDA_RECT m_BoundingBox;
NETINFO_LIST m_NetInfo; ///< net info list (name, design constraints ..
+ RN_DATA* m_ratsnest;
BOARD_DESIGN_SETTINGS m_designSettings;
ZONE_SETTINGS m_zoneSettings;
@@ -355,6 +357,16 @@ public:
*/
BOARD_ITEM* Remove( BOARD_ITEM* aBoardItem );
+ /**
+ * Function GetRatsnest()
+ * returns list of missing connections between components/tracks.
+ * @return RATSNEST* is an object that contains informations about missing connections.
+ */
+ RN_DATA* GetRatsnest() const
+ {
+ return m_ratsnest;
+ }
+
/**
* Function DeleteMARKERs
* deletes ALL MARKERS from the board.
diff --git a/pcbnew/pcbframe.cpp b/pcbnew/pcbframe.cpp
index 7e415cd831..c639f88024 100644
--- a/pcbnew/pcbframe.cpp
+++ b/pcbnew/pcbframe.cpp
@@ -59,6 +59,13 @@
#include
#include
+#include
+#include
+#include
+#include
+#include
+#include
+
#include
#include
@@ -329,6 +336,16 @@ PCB_EDIT_FRAME::PCB_EDIT_FRAME( wxWindow* parent, const wxString& title,
SetBoard( new BOARD() );
+ if( m_galCanvas )
+ {
+ ViewReloadBoard( m_Pcb );
+
+ // update the tool manager with the new board and its view.
+ if( m_toolManager )
+ m_toolManager->SetEnvironment( m_Pcb, m_galCanvas->GetView(),
+ m_galCanvas->GetViewControls(), this );
+ }
+
// Create the PCB_LAYER_WIDGET *after* SetBoard():
wxFont font = wxSystemSettings::GetFont( wxSYS_DEFAULT_GUI_FONT );
@@ -522,6 +539,106 @@ PCB_EDIT_FRAME::~PCB_EDIT_FRAME()
}
+void PCB_EDIT_FRAME::SetBoard( BOARD* aBoard )
+{
+ PCB_BASE_FRAME::SetBoard( aBoard );
+
+ if( m_galCanvas )
+ {
+ ViewReloadBoard( aBoard );
+
+ // update the tool manager with the new board and its view.
+ if( m_toolManager )
+ m_toolManager->SetEnvironment( aBoard, m_galCanvas->GetView(),
+ m_galCanvas->GetViewControls(), this );
+ }
+}
+
+
+void PCB_EDIT_FRAME::ViewReloadBoard( const BOARD* aBoard ) const
+{
+ KIGFX::VIEW* view = m_galCanvas->GetView();
+ view->Clear();
+
+ // All of PCB drawing elements should be added to the VIEW
+ // in order to be displayed
+
+ // Load zones
+ for( int i = 0; i < aBoard->GetAreaCount(); ++i )
+ {
+ view->Add( (KIGFX::VIEW_ITEM*) ( aBoard->GetArea( i ) ) );
+ }
+
+ // Load drawings
+ for( BOARD_ITEM* drawing = aBoard->m_Drawings; drawing; drawing = drawing->Next() )
+ {
+ view->Add( drawing );
+ }
+
+ // Load tracks
+ for( TRACK* track = aBoard->m_Track; track; track = track->Next() )
+ {
+ view->Add( track );
+ }
+
+ // Load modules and its additional elements
+ for( MODULE* module = aBoard->m_Modules; module; module = module->Next() )
+ {
+ // Load module's pads
+ for( D_PAD* pad = module->Pads().GetFirst(); pad; pad = pad->Next() )
+ {
+ view->Add( pad );
+ }
+
+ // Load module's drawing (mostly silkscreen)
+ for( BOARD_ITEM* drawing = module->GraphicalItems().GetFirst(); drawing;
+ drawing = drawing->Next() )
+ {
+ view->Add( drawing );
+ }
+
+ // Load module's texts (name and value)
+ view->Add( &module->Reference() );
+ view->Add( &module->Value() );
+
+ // Add the module itself
+ view->Add( module );
+ }
+
+ // Segzones (equivalent of ZONE_CONTAINER for legacy boards)
+ for( SEGZONE* zone = aBoard->m_Zone; zone; zone = zone->Next() )
+ {
+ view->Add( zone );
+ }
+
+ // Add an entry for the worksheet layout
+ KIGFX::WORKSHEET_VIEWITEM* worksheet = new KIGFX::WORKSHEET_VIEWITEM(
+ std::string( aBoard->GetFileName().mb_str() ),
+ std::string( GetScreenDesc().mb_str() ),
+ &GetPageSettings(), &GetTitleBlock() );
+ BASE_SCREEN* screen = GetScreen();
+ if( screen != NULL )
+ {
+ worksheet->SetSheetNumber( GetScreen()->m_ScreenNumber );
+ worksheet->SetSheetCount( GetScreen()->m_NumberOfScreens );
+ }
+
+ view->Add( worksheet );
+
+ // Add an entry for the ratsnest
+ RN_DATA* ratsnest = aBoard->GetRatsnest();
+ ratsnest->ProcessBoard();
+ ratsnest->Recalculate();
+ view->Add( new KIGFX::RATSNEST_VIEWITEM( ratsnest ) );
+
+ view->SetPanBoundary( worksheet->ViewBBox() );
+ view->RecacheAllItems( true );
+
+ if( m_galCanvasActive )
+ m_galCanvas->Refresh();
+}
+
+
bool PCB_EDIT_FRAME::isAutoSaveRequired() const
{
return GetScreen()->IsSave();
@@ -633,6 +750,17 @@ void PCB_EDIT_FRAME::Show3D_Frame( wxCommandEvent& event )
}
+void PCB_EDIT_FRAME::UseGalCanvas( bool aEnable )
+{
+ EDA_DRAW_FRAME::UseGalCanvas( aEnable );
+
+ m_toolManager->SetEnvironment( m_Pcb, m_galCanvas->GetView(),
+ m_galCanvas->GetViewControls(), this );
+
+ ViewReloadBoard( m_Pcb );
+}
+
+
void PCB_EDIT_FRAME::SwitchCanvas( wxCommandEvent& aEvent )
{
int id = aEvent.GetId();
@@ -795,7 +923,7 @@ void PCB_EDIT_FRAME::setHighContrastLayer( LAYER_NUM aLayer )
GetNetnameLayer( aLayer ), ITEM_GAL_LAYER( VIAS_VISIBLE ),
ITEM_GAL_LAYER( VIAS_HOLES_VISIBLE ), ITEM_GAL_LAYER( PADS_VISIBLE ),
ITEM_GAL_LAYER( PADS_HOLES_VISIBLE ), ITEM_GAL_LAYER( PADS_NETNAMES_VISIBLE ),
- ITEM_GAL_LAYER( GP_OVERLAY )
+ ITEM_GAL_LAYER( GP_OVERLAY ), ITEM_GAL_LAYER( RATSNEST_VISIBLE )
};
for( unsigned int i = 0; i < sizeof( layers ) / sizeof( LAYER_NUM ); ++i )
@@ -835,7 +963,7 @@ void PCB_EDIT_FRAME::setTopLayer( LAYER_NUM aLayer )
GetNetnameLayer( aLayer ), ITEM_GAL_LAYER( VIAS_VISIBLE ),
ITEM_GAL_LAYER( VIAS_HOLES_VISIBLE ), ITEM_GAL_LAYER( PADS_VISIBLE ),
ITEM_GAL_LAYER( PADS_HOLES_VISIBLE ), ITEM_GAL_LAYER( PADS_NETNAMES_VISIBLE ),
- ITEM_GAL_LAYER( GP_OVERLAY ), DRAW_N
+ ITEM_GAL_LAYER( GP_OVERLAY ), ITEM_GAL_LAYER( RATSNEST_VISIBLE ), DRAW_N
};
for( unsigned int i = 0; i < sizeof( layers ) / sizeof( LAYER_NUM ); ++i )
diff --git a/pcbnew/ratsnest_data.cpp b/pcbnew/ratsnest_data.cpp
new file mode 100644
index 0000000000..85ca642808
--- /dev/null
+++ b/pcbnew/ratsnest_data.cpp
@@ -0,0 +1,856 @@
+/*
+ * This program source code file is part of KICAD, a free EDA CAD application.
+ *
+ * Copyright (C) 2013 CERN
+ * @author Maciej Suminski
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program 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 General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, you may find one here:
+ * http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
+ * or you may search the http://www.gnu.org website for the version 2 license,
+ * or you may write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
+ */
+
+/**
+ * @file ratsnest_data.cpp
+ * @brief Class that computes missing connections on a PCB.
+ */
+
+#include
+
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+uint64_t getDistance( const RN_NODE_PTR& aNode1, const RN_NODE_PTR& aNode2 )
+{
+ // Drop the least significant bits to avoid overflow
+ int64_t x = ( aNode1->GetX() - aNode2->GetX() ) >> 16;
+ int64_t y = ( aNode1->GetY() - aNode2->GetY() ) >> 16;
+
+ // We do not need sqrt() here, as the distance is computed only for comparison
+ return ( x * x + y * y );
+}
+
+
+bool sortDistance( const RN_NODE_PTR& aOrigin, const RN_NODE_PTR& aNode1,
+ const RN_NODE_PTR& aNode2 )
+{
+ return getDistance( aOrigin, aNode1 ) < getDistance( aOrigin, aNode2 );
+}
+
+
+bool sortWeight( const RN_EDGE_PTR& aEdge1, const RN_EDGE_PTR& aEdge2 )
+{
+ return aEdge1->getWeight() < aEdge2->getWeight();
+}
+
+
+bool sortArea( const RN_POLY& aP1, const RN_POLY& aP2 )
+{
+ return aP1.m_bbox.GetArea() < aP2.m_bbox.GetArea();
+}
+
+
+bool isEdgeConnectingNode( const RN_EDGE_PTR& aEdge, const RN_NODE_PTR& aNode )
+{
+ return ( aEdge->getSourceNode().get() == aNode.get() ) ||
+ ( aEdge->getTargetNode().get() == aNode.get() );
+}
+
+
+std::vector* kruskalMST( RN_LINKS::RN_EDGE_LIST& aEdges,
+ const std::vector& aNodes )
+{
+ unsigned int nodeNumber = aNodes.size();
+ unsigned int mstExpectedSize = nodeNumber - 1;
+ unsigned int mstSize = 0;
+
+ // The output
+ std::vector* mst = new std::vector;
+ mst->reserve( mstExpectedSize );
+
+ // Set tags for marking cycles
+ boost::unordered_map tags;
+ unsigned int tag = 0;
+ BOOST_FOREACH( const RN_NODE_PTR& node, aNodes )
+ tags[node] = tag++;
+
+ // Lists of nodes connected together (subtrees) to detect cycles in the graph
+ std::vector > cycles( nodeNumber );
+ for( unsigned int i = 0; i < nodeNumber; ++i )
+ cycles[i].push_back( i );
+
+ // Kruskal algorithm requires edges to be sorted by their weight
+ aEdges.sort( sortWeight );
+
+ while( mstSize < mstExpectedSize && !aEdges.empty() )
+ {
+ RN_EDGE_PTR& dt = *aEdges.begin();
+
+ int srcTag = tags[dt->getSourceNode()];
+ int trgTag = tags[dt->getTargetNode()];
+
+ // Check if by adding this edge we are going to join two different forests
+ if( srcTag != trgTag )
+ {
+ // Update tags
+ std::list::iterator it, itEnd;
+ for( it = cycles[trgTag].begin(), itEnd = cycles[trgTag].end(); it != itEnd; ++it )
+ tags[aNodes[*it]] = srcTag;
+
+ // Move nodes that were marked with old tag to the list marked with the new tag
+ cycles[srcTag].splice( cycles[srcTag].end(), cycles[trgTag] );
+
+ if( dt->getWeight() == 0 ) // Skip already existing connections (weight == 0)
+ mstExpectedSize--;
+ else
+ {
+ mst->push_back( dt );
+ ++mstSize;
+ }
+ }
+
+ // Remove the edge that was just processed
+ aEdges.erase( aEdges.begin() );
+ }
+
+ // Probably we have discarded some of edges, so reduce the size
+ mst->resize( mstSize );
+
+ return mst;
+}
+
+
+void RN_NET::validateEdge( RN_EDGE_PTR& aEdge )
+{
+ RN_NODE_PTR source = aEdge->getSourceNode();
+ RN_NODE_PTR target = aEdge->getTargetNode();
+ bool valid = true;
+
+ // If any of nodes belonging to the edge has the flag set,
+ // change it to the closest node that has flag cleared
+ if( source->GetFlag() )
+ {
+ valid = false;
+
+ std::list closest = GetClosestNodes( source, WITHOUT_FLAG() );
+ BOOST_FOREACH( RN_NODE_PTR& node, closest )
+ {
+ if( node && node != target )
+ {
+ source = node;
+ break;
+ }
+ }
+ }
+
+ if( target->GetFlag() )
+ {
+ valid = false;
+
+ WITHOUT_FLAG without_flag;
+ std::list closest = GetClosestNodes( target, WITHOUT_FLAG() );
+ BOOST_FOREACH( RN_NODE_PTR& node, closest )
+ {
+ if( node && node != source )
+ {
+ target = node;
+ break;
+ }
+ }
+ }
+
+ // Replace an invalid edge with new, valid one
+ if( !valid )
+ aEdge.reset( new RN_EDGE_MST( source, target ) );
+}
+
+
+const RN_NODE_PTR& RN_LINKS::AddNode( int aX, int aY )
+{
+ RN_NODE_SET::iterator node;
+ bool wasNewElement;
+
+ boost::tie( node, wasNewElement ) = m_nodes.emplace( boost::make_shared( aX, aY ) );
+ (*node)->IncRefCount(); // TODO use the shared_ptr use_count
+
+ return *node;
+}
+
+
+void RN_LINKS::RemoveNode( const RN_NODE_PTR& aNode )
+{
+ aNode->DecRefCount(); // TODO use the shared_ptr use_count
+
+ if( aNode->GetRefCount() == 0 )
+ m_nodes.erase( aNode );
+}
+
+
+const RN_EDGE_PTR& RN_LINKS::AddConnection( const RN_NODE_PTR& aNode1, const RN_NODE_PTR& aNode2,
+ unsigned int aDistance )
+{
+ m_edges.push_back( boost::make_shared( aNode1, aNode2, aDistance ) );
+
+ return m_edges.back();
+}
+
+
+void RN_LINKS::RemoveConnection( const RN_EDGE_PTR& aEdge )
+{
+ m_edges.remove( aEdge );
+}
+
+
+void RN_NET::compute()
+{
+ const RN_LINKS::RN_NODE_SET& boardNodes = m_links.GetNodes();
+ const RN_LINKS::RN_EDGE_LIST& boardEdges = m_links.GetConnections();
+
+ // Special case that does need so complicated algorithm
+ if( boardNodes.size() == 2 )
+ {
+ m_rnEdges.reset( new std::vector( 0 ) );
+
+ // Check if the only possible connection exists
+ if( boardEdges.size() == 0 )
+ {
+ RN_LINKS::RN_NODE_SET::iterator last = ++boardNodes.begin();
+
+ // There can be only one possible connection, but it is missing
+ m_rnEdges->push_back( boost::make_shared( *boardNodes.begin(), *last ) );
+ }
+
+ return;
+ }
+ else if( boardNodes.size() == 1 ) // This case is even simpler
+ {
+ m_rnEdges.reset( new std::vector( 0 ) );
+
+ return;
+ }
+
+ // Move and sort (sorting speeds up) all nodes to a vector for the Delaunay triangulation
+ std::vector nodes( boardNodes.size() );
+ std::partial_sort_copy( boardNodes.begin(), boardNodes.end(), nodes.begin(), nodes.end() );
+
+ TRIANGULATOR triangulator;
+ triangulator.createDelaunay( nodes.begin(), nodes.end() );
+ boost::scoped_ptr triangEdges( triangulator.getEdges() );
+
+ // Compute weight/distance for edges resulting from triangulation
+ RN_LINKS::RN_EDGE_LIST::iterator eit, eitEnd;
+ for( eit = (*triangEdges).begin(), eitEnd = (*triangEdges).end(); eit != eitEnd; ++eit )
+ (*eit)->setWeight( getDistance( (*eit)->getSourceNode(), (*eit)->getTargetNode() ) );
+
+ // Add the currently existing connections list to the results of triangulation
+ std::copy( boardEdges.begin(), boardEdges.end(), std::front_inserter( *triangEdges ) );
+
+ // Get the minimal spanning tree
+ m_rnEdges.reset( kruskalMST( *triangEdges, nodes ) );
+}
+
+
+void RN_NET::clearNode( const RN_NODE_PTR& aNode )
+{
+ std::vector::iterator newEnd;
+
+ // Remove all ratsnest edges for associated with the node
+ newEnd = std::remove_if( m_rnEdges->begin(), m_rnEdges->end(),
+ boost::bind( isEdgeConnectingNode, _1, aNode ) );
+
+ m_rnEdges->resize( std::distance( m_rnEdges->begin(), newEnd ) );
+}
+
+
+RN_POLY::RN_POLY( const CPolyPt* aBegin, const CPolyPt* aEnd, const ZONE_CONTAINER* aParent,
+ RN_LINKS& aConnections, const BOX2I& aBBox ) :
+ m_parent( aParent), m_begin( aBegin ), m_end( aEnd ), m_bbox( aBBox )
+{
+ m_node = aConnections.AddNode( m_begin->x, m_begin->y );
+
+ // Mark it as not feasible as a destination of ratsnest edges
+ // (edges coming out from a polygon vertex look weird)
+ m_node->SetFlag( true );
+}
+
+
+bool RN_POLY::HitTest( const RN_NODE_PTR& aNode ) const
+{
+ long xt = aNode->GetX();
+ long yt = aNode->GetY();
+
+ // If the point lies outside the bounding box, there is no point to check it further
+ if( !m_bbox.Contains( xt, yt ) )
+ return false;
+
+ long xNew, yNew, xOld, yOld, x1, y1, x2, y2;
+ bool inside = false;
+
+ // For the first loop we have to use the last point as the previous point
+ xOld = m_end->x;
+ yOld = m_end->y;
+
+ for( const CPolyPt* point = m_begin; point <= m_end; ++point )
+ {
+ xNew = point->x;
+ yNew = point->y;
+
+ // Swap points if needed, so always x2 >= x1
+ if( xNew > xOld )
+ {
+ x1 = xOld; y1 = yOld;
+ x2 = xNew; y2 = yNew;
+ }
+ else
+ {
+ x1 = xNew; y1 = yNew;
+ x2 = xOld; y2 = yOld;
+ }
+
+ if( ( xNew < xt ) == ( xt <= xOld ) /* edge "open" at left end */
+ && ( yt - y1 ) * ( x2 - x1 ) < ( y2 - y1 ) * ( xt - x1 ) )
+ {
+ inside = !inside;
+ }
+
+ xOld = xNew;
+ yOld = yNew;
+ }
+
+ return inside;
+}
+
+
+void RN_NET::Update()
+{
+ // Add edges resulting from nodes being connected by zones
+ processZones();
+
+ compute();
+
+ BOOST_FOREACH( RN_EDGE_PTR& edge, *m_rnEdges )
+ validateEdge( edge );
+
+ m_dirty = false;
+}
+
+
+void RN_NET::AddItem( const D_PAD* aPad )
+{
+ RN_NODE_PTR nodePtr = m_links.AddNode( aPad->GetPosition().x, aPad->GetPosition().y );
+ m_pads[aPad] = nodePtr;
+
+ m_dirty = true;
+}
+
+
+void RN_NET::AddItem( const SEGVIA* aVia )
+{
+ m_vias[aVia] = m_links.AddNode( aVia->GetPosition().x, aVia->GetPosition().y );
+
+ m_dirty = true;
+}
+
+
+void RN_NET::AddItem( const TRACK* aTrack )
+{
+ RN_NODE_PTR start = m_links.AddNode( aTrack->GetStart().x, aTrack->GetStart().y );
+ RN_NODE_PTR end = m_links.AddNode( aTrack->GetEnd().x, aTrack->GetEnd().y );
+
+ m_tracks[aTrack] = m_links.AddConnection( start, end );
+
+ m_dirty = true;
+}
+
+
+void RN_NET::AddItem( const ZONE_CONTAINER* aZone )
+{
+ // Prepare a list of polygons (every zone can contain one or more polygons)
+ const std::vector& polyPoints = aZone->GetFilledPolysList().GetList();
+ if( polyPoints.size() == 0 )
+ return;
+
+ // Origin and end of bounding box for a polygon
+ VECTOR2I origin( polyPoints[0].x, polyPoints[0].y );
+ VECTOR2I end( polyPoints[0].x, polyPoints[0].y );
+ int idxStart = 0;
+
+ // Extract polygons from zones
+ for( unsigned int i = 0; i < polyPoints.size(); ++i )
+ {
+ const CPolyPt& point = polyPoints[i];
+
+ // Determine bounding box
+ if( point.x < origin.x )
+ origin.x = point.x;
+ else if( point.x > end.x )
+ end.x = point.x;
+
+ if( point.y < origin.y )
+ origin.y = point.y;
+ else if( point.y > end.y )
+ end.y = point.y;
+
+ if( point.end_contour )
+ {
+ // The last vertex is enclosing the polygon (it repeats at the beginning and
+ // at the end), so we skip it
+ m_zonePolygons[aZone].push_back( RN_POLY( &polyPoints[idxStart], &point, aZone,
+ m_links, BOX2I( origin, end - origin ) ) );
+
+ idxStart = i + 1;
+ origin.x = polyPoints[idxStart].x;
+ origin.y = polyPoints[idxStart].y;
+ end.x = polyPoints[idxStart].x;
+ end.y = polyPoints[idxStart].y;
+ }
+ }
+
+ m_dirty = true;
+}
+
+
+void RN_NET::RemoveItem( const D_PAD* aPad )
+{
+ RN_NODE_PTR& node = m_pads[aPad];
+ if( !node )
+ return;
+
+ // Remove edges associated with the node
+ clearNode( node );
+ m_links.RemoveNode( node );
+
+ m_pads.erase( aPad );
+
+ m_dirty = true;
+}
+
+
+void RN_NET::RemoveItem( const SEGVIA* aVia )
+{
+ RN_NODE_PTR& node = m_vias[aVia];
+ if( !node )
+ return;
+
+ // Remove edges associated with the node
+ clearNode( node );
+ m_links.RemoveNode( node );
+
+ m_vias.erase( aVia );
+
+ m_dirty = true;
+}
+
+
+void RN_NET::RemoveItem( const TRACK* aTrack )
+{
+ RN_EDGE_PTR& edge = m_tracks[aTrack];
+ if( !edge )
+ return;
+
+ // Save nodes, so they can be cleared later
+ const RN_NODE_PTR& aBegin = edge->getSourceNode();
+ const RN_NODE_PTR& aEnd = edge->getTargetNode();
+ m_links.RemoveConnection( edge );
+
+ // Remove nodes associated with the edge. It is done in a safe way, there is a check
+ // if nodes are not used by other edges.
+ clearNode( aBegin );
+ clearNode( aEnd );
+ m_links.RemoveNode( aBegin );
+ m_links.RemoveNode( aEnd );
+
+ m_tracks.erase( aTrack );
+
+ m_dirty = true;
+}
+
+
+void RN_NET::RemoveItem( const ZONE_CONTAINER* aZone )
+{
+ // Remove all subpolygons that make the zone
+ std::deque& polygons = m_zonePolygons[aZone];
+ BOOST_FOREACH( RN_POLY& polygon, polygons )
+ m_links.RemoveNode( polygon.GetNode() );
+ polygons.clear();
+
+ // Remove all connections added by the zone
+ std::deque& edges = m_zoneConnections[aZone];
+ BOOST_FOREACH( RN_EDGE_PTR& edge, edges )
+ m_links.RemoveConnection( edge );
+ edges.clear();
+
+ m_dirty = true;
+}
+
+
+const RN_NODE_PTR RN_NET::GetClosestNode( const RN_NODE_PTR& aNode ) const
+{
+ const RN_LINKS::RN_NODE_SET& nodes = m_links.GetNodes();
+ RN_LINKS::RN_NODE_SET::const_iterator it, itEnd;
+
+ unsigned int minDistance = std::numeric_limits::max();
+ RN_NODE_PTR closest;
+
+ for( it = nodes.begin(), itEnd = nodes.end(); it != itEnd; ++it )
+ {
+ // Obviously the distance between node and itself is the shortest,
+ // that's why we have to skip it
+ if( *it != aNode )
+ {
+ unsigned int distance = getDistance( *it, aNode );
+ if( distance < minDistance )
+ {
+ minDistance = distance;
+ closest = *it;
+ }
+ }
+ }
+
+ return closest;
+}
+
+
+const RN_NODE_PTR RN_NET::GetClosestNode( const RN_NODE_PTR& aNode, RN_NODE_FILTER aFilter ) const
+{
+ const RN_LINKS::RN_NODE_SET& nodes = m_links.GetNodes();
+ RN_LINKS::RN_NODE_SET::const_iterator it, itEnd;
+
+ unsigned int minDistance = std::numeric_limits::max();
+ RN_NODE_PTR closest;
+
+ for( it = nodes.begin(), itEnd = nodes.end(); it != itEnd; ++it )
+ {
+ RN_NODE_PTR baseNode = *it;
+
+ // Obviously the distance between node and itself is the shortest,
+ // that's why we have to skip it
+ if( *it != aNode && aFilter( baseNode ) )
+ {
+ unsigned int distance = getDistance( *it, aNode );
+ if( distance < minDistance )
+ {
+ minDistance = distance;
+ closest = *it;
+ }
+ }
+ }
+
+ return closest;
+}
+
+
+std::list RN_NET::GetClosestNodes( const RN_NODE_PTR& aNode, int aNumber ) const
+{
+ std::list closest;
+ const RN_LINKS::RN_NODE_SET& nodes = m_links.GetNodes();
+
+ // Copy nodes
+ BOOST_FOREACH( const RN_NODE_PTR& node, nodes )
+ closest.push_back( node );
+
+ // Sort by the distance from aNode
+ closest.sort( boost::bind( sortDistance, aNode, _1, _2 ) );
+
+ // Remove the first node (==aNode), as it is surely located within the smallest distance
+ closest.pop_front();
+
+ // Trim the result to the asked size
+ if( aNumber > 0 )
+ closest.resize( std::min( (size_t)aNumber, nodes.size() ) );
+
+ return closest;
+}
+
+
+std::list RN_NET::GetClosestNodes( const RN_NODE_PTR& aNode,
+ RN_NODE_FILTER aFilter, int aNumber ) const
+{
+ std::list closest;
+ const RN_LINKS::RN_NODE_SET& nodes = m_links.GetNodes();
+
+ // Copy nodes
+ BOOST_FOREACH( const RN_NODE_PTR& node, nodes )
+ closest.push_back( node );
+
+ // Sort by the distance from aNode
+ closest.sort( boost::bind( sortDistance, aNode, _1, _2 ) );
+
+ // Remove the first node (==aNode), as it is surely located within the smallest distance
+ closest.pop_front();
+
+ // Filter out by condition
+ std::remove_if( closest.begin(), closest.end(), aFilter );
+
+ // Trim the result to the asked size
+ if( aNumber > 0 )
+ closest.resize( std::min( static_cast( aNumber ), nodes.size() ) );
+
+ return closest;
+}
+
+
+std::list RN_NET::GetNodes( const BOARD_CONNECTED_ITEM* aItem ) const
+{
+ std::list nodes;
+
+ switch( aItem->Type() )
+ {
+ case PCB_PAD_T:
+ {
+ const D_PAD* pad = static_cast( aItem );
+ nodes.push_back( m_pads.at( pad ) );
+ }
+ break;
+
+ case PCB_VIA_T:
+ {
+ const SEGVIA* via = static_cast( aItem );
+ nodes.push_back( m_vias.at( via ) );
+ }
+ break;
+
+ case PCB_TRACE_T:
+ {
+ const TRACK* track = static_cast( aItem );
+ RN_EDGE_PTR edge = m_tracks.at( track );
+
+ nodes.push_back( edge->getSourceNode() );
+ nodes.push_back( edge->getTargetNode() );
+ }
+ break;
+
+ default:
+ break;
+ }
+
+ return nodes;
+}
+
+
+void RN_NET::ClearSimple()
+{
+ BOOST_FOREACH( const RN_NODE_PTR& node, m_simpleNodes )
+ node->SetFlag( false );
+
+ m_simpleNodes.clear();
+}
+
+
+void RN_DATA::AddSimple( const BOARD_CONNECTED_ITEM* aItem )
+{
+ int net = aItem->GetNet();
+ if( net < 1 ) // do not process unconnected items
+ return;
+
+ // Get list of nodes responding to the item
+ std::list nodes = m_nets[net].GetNodes( aItem );
+ std::list::iterator it, itEnd;
+
+ for( it = nodes.begin(), itEnd = nodes.end(); it != itEnd; ++it )
+ m_nets[net].AddSimpleNode( *it );
+}
+
+
+void RN_DATA::AddSimple( const MODULE* aModule )
+{
+ for( const D_PAD* pad = aModule->Pads().GetFirst(); pad; pad = pad->Next() )
+ AddSimple( pad );
+}
+
+
+void RN_NET::processZones()
+{
+ BOOST_FOREACH( std::deque& edges, m_zoneConnections | boost::adaptors::map_values )
+ {
+ BOOST_FOREACH( RN_EDGE_PTR& edge, edges )
+ m_links.RemoveConnection( edge );
+
+ edges.clear();
+ }
+
+ RN_LINKS::RN_NODE_SET candidates = m_links.GetNodes();
+
+ BOOST_FOREACH( std::deque& polygons, m_zonePolygons | boost::adaptors::map_values )
+ {
+ RN_LINKS::RN_NODE_SET::iterator point, pointEnd;
+ std::deque::iterator poly, polyEnd;
+
+ // Sorting by area should speed up the processing, as smaller polygons are computed
+ // faster and may reduce the number of points for further checks
+ std::sort( polygons.begin(), polygons.end(), sortArea );
+
+ for( poly = polygons.begin(), polyEnd = polygons.end(); poly != polyEnd; ++poly )
+ {
+ point = candidates.begin();
+ pointEnd = candidates.end();
+
+ while( point != pointEnd )
+ {
+ if( poly->HitTest( *point ) )
+ {
+ const RN_EDGE_PTR& connection = m_links.AddConnection( poly->GetNode(), *point );
+ m_zoneConnections[poly->GetParent()].push_back( connection );
+
+ // This point already belongs to a polygon, we do not need to check it anymore
+ point = candidates.erase( point );
+ pointEnd = candidates.end();
+ }
+ else
+ {
+ ++point;
+ }
+ }
+ }
+ }
+}
+
+
+void RN_DATA::updateNet( int aNetCode )
+{
+ assert( aNetCode < (int) m_nets.size() );
+ if( aNetCode < 1 )
+ return;
+
+ m_nets[aNetCode].ClearSimple();
+ m_nets[aNetCode].Update();
+}
+
+
+void RN_DATA::Update( const BOARD_CONNECTED_ITEM* aItem )
+{
+ int net = aItem->GetNet();
+ if( net < 1 ) // do not process unconnected items
+ return;
+
+ switch( aItem->Type() )
+ {
+ case PCB_PAD_T:
+ {
+ const D_PAD* pad = static_cast( aItem );
+ m_nets[net].RemoveItem( pad );
+ m_nets[net].AddItem( pad );
+ }
+ break;
+
+ case PCB_TRACE_T:
+ {
+ const TRACK* track = static_cast( aItem );
+ m_nets[net].RemoveItem( track );
+ m_nets[net].AddItem( track );
+ }
+ break;
+
+ case PCB_VIA_T:
+ {
+ const SEGVIA* via = static_cast( aItem );
+ m_nets[net].RemoveItem( via );
+ m_nets[net].AddItem( via );
+ }
+ break;
+
+ case PCB_ZONE_AREA_T:
+ {
+ const ZONE_CONTAINER* zone = static_cast( aItem );
+ m_nets[net].RemoveItem( zone);
+ m_nets[net].AddItem( zone );
+ }
+ break;
+
+ default:
+ break;
+ }
+}
+
+
+void RN_DATA::Update( const MODULE* aModule )
+{
+ for( const D_PAD* pad = aModule->Pads().GetFirst(); pad; pad = pad->Next() )
+ {
+ int net = pad->GetNet();
+
+ if( net > 0 ) // do not process unconnected items
+ {
+ m_nets[net].RemoveItem( pad );
+ m_nets[net].AddItem( pad );
+ }
+ }
+}
+
+
+void RN_DATA::ProcessBoard()
+{
+ m_nets.clear();
+ m_nets.resize( m_board->GetNetCount() );
+
+ // Iterate over all items that may need to be connected
+ for( MODULE* module = m_board->m_Modules; module; module = module->Next() )
+ {
+ for( D_PAD* pad = module->Pads().GetFirst(); pad; pad = pad->Next() )
+ m_nets[pad->GetNet()].AddItem( pad );
+ }
+
+ for( TRACK* track = m_board->m_Track; track; track = track->Next() )
+ {
+ if( track->Type() == PCB_VIA_T )
+ m_nets[track->GetNet()].AddItem( static_cast( track ) );
+ else if( track->Type() == PCB_TRACE_T )
+ m_nets[track->GetNet()].AddItem( track );
+ }
+
+ for( int i = 0; i < m_board->GetAreaCount(); ++i )
+ {
+ ZONE_CONTAINER* zone = m_board->GetArea( i );
+ m_nets[zone->GetNet()].AddItem( zone );
+ }
+}
+
+
+void RN_DATA::Recalculate( int aNet )
+{
+ if( aNet < 0 ) // Recompute everything
+ {
+ // Start with net number 1, as 0 stand for not connected
+ for( unsigned int i = 1; i < m_board->GetNetCount(); ++i )
+ {
+ if( m_nets[i].IsDirty() )
+ updateNet( i );
+ }
+ }
+ else // Recompute only specific net
+ {
+ updateNet( aNet );
+ }
+}
+
+
+void RN_DATA::ClearSimple()
+{
+ BOOST_FOREACH( RN_NET& net, m_nets )
+ net.ClearSimple();
+}
diff --git a/pcbnew/ratsnest_data.h b/pcbnew/ratsnest_data.h
new file mode 100644
index 0000000000..dc5e4016b1
--- /dev/null
+++ b/pcbnew/ratsnest_data.h
@@ -0,0 +1,595 @@
+/*
+ * This program source code file is part of KICAD, a free EDA CAD application.
+ *
+ * Copyright (C) 2013 CERN
+ * @author Maciej Suminski
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program 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 General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, you may find one here:
+ * http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
+ * or you may search the http://www.gnu.org website for the version 2 license,
+ * or you may write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
+ */
+
+/**
+ * @file ratsnest_data.h
+ * @brief Class that computes missing connections on a PCB.
+ */
+
+#ifndef RATSNEST_DATA_H
+#define RATSNEST_DATA_H
+
+#include
+#include
+
+#include