Added ratsnest for GAL

common/CMakeLists.txt

-
 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

include/ttl/halfedge/hedart.h

+/*
+ * 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
+ * <http://www.gnu.org/licenses/>.
+ *
+ * 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 <ttl/halfedge/hetriang.h>
+
+
+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

include/ttl/halfedge/hetraits.h

+/*
+ * 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
+ * <http://www.gnu.org/licenses/>.
+ *
+ * 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 <ttl/halfedge/hetriang.h>
+#include <ttl/halfedge/hedart.h>
+
+
+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.<br>
+    *
+    *   \note
+    *   This function is also present in the TTL as ttl::swapTestDelaunay.<br>
+    *   Thus, the function can be implemented simply as:
+    *   \code
+    *   { return ttl::swapTestDelaunay<TTLtraits>(dart); }
+    *   \endcode
+    */
+    //static bool swapTestDelaunay(const Dart& dart) {
+    //  return ttl::swapTestDelaunay<TTLtraits>(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<TTLtraits>(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.
+    *
+    *   <center>
+    *   \image html swapEdge.gif
+    *   </center>
+    * 
+    *   \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.
+    *
+    *   <center>
+    *   \image html splitTriangle.gif
+    *   </center>
+    *
+    *   \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.
+    *
+    *   <center>
+    *   \image html reverse_splitTriangle.gif
+    *   </center>
+    */
+    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

include/ttl/halfedge/hetriang.h

+/*
+ * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
+ * Applied Mathematics, Norway.
+ * Copyright (C) 2013 CERN
+ * @author Maciej Suminski <maciej.suminski@cern.ch>
+ *
+ * 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
+ * <http://www.gnu.org/licenses/>.
+ *
+ * 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 <list>
+#include <vector>
+#include <iostream>
+#include <fstream>
+#include <ttl/ttl.h>
+#include <ttl/ttl_util.h>
+#include <boost/shared_ptr.hpp>
+
+//--------------------------------------------------------------------------------------------------
+// The half-edge data structure
+//--------------------------------------------------------------------------------------------------
+
+namespace hed {
+  // Helper typedefs
+  class Node;
+  class Edge;
+  typedef boost::shared_ptr<Node> NodePtr;
+  typedef boost::shared_ptr<Edge> EdgePtr;
+  typedef std::vector<NodePtr> 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<EdgePtr> 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<EdgePtr>& 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<EdgePtr>* 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<NodePtr>* 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

include/ttl/ttl.h

+/*
+ * 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
+ * <http://www.gnu.org/licenses/>.
+ *
+ * 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 <list>
+#include <iterator>
+
+// 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<class TraitsType, class DartType, class PointType>
+    bool insertNode(DartType& dart, PointType& point);
+
+  template<class TraitsType, class DartType>
+    void removeRectangularBoundary(DartType& dart);
+
+  template<class TraitsType, class DartType>
+    void removeNode(DartType& dart);
+
+  template<class TraitsType, class DartType>
+    void removeBoundaryNode(DartType& dart);
+
+  template<class TraitsType, class DartType>
+    void removeInteriorNode(DartType& dart);
+
+
+  // Topological and Geometric Queries
+  // ---------------------------------
+  template<class TraitsType, class PointType, class DartType>
+    bool locateFaceSimplest(const PointType& point, DartType& dart);
+
+  template<class TraitsType, class PointType, class DartType>
+    bool locateTriangle(const PointType& point, DartType& dart);
+
+  template<class TraitsType, class PointType, class DartType>
+    bool inTriangleSimplest(const PointType& point, const DartType& dart);
+
+  template<class TraitsType, class PointType, class DartType>
+    bool inTriangle(const PointType& point, const DartType& dart);
+
+  template<class DartType, class DartListType>
+    void getBoundary(const DartType& dart, DartListType& boundary);
+
+  template<class DartType>
+    bool isBoundaryEdge(const DartType& dart);
+
+  template<class DartType>
+    bool isBoundaryFace(const DartType& dart);
+
+  template<class DartType>
+    bool isBoundaryNode(const DartType& dart);
+
+  template<class DartType>
+    int getDegreeOfNode(const DartType& dart);
+
+  template<class DartType, class DartListType>
+    void get_0_orbit_interior(const DartType& dart, DartListType& orbit);
+
+  template<class DartType, class DartListType>
+    void get_0_orbit_boundary(const DartType& dart, DartListType& orbit);
+
+  template<class DartType>
+    bool same_0_orbit(const DartType& d1, const DartType& d2);
+
+  template<class DartType>
+    bool same_1_orbit(const DartType& d1, const DartType& d2);
+
+  template<class DartType>
+    bool same_2_orbit(const DartType& d1, const DartType& d2);
+
+  template <class TraitsType, class DartType>
+    bool swappableEdge(const DartType& dart, bool allowDegeneracy = false);
+
+  template<class DartType>
+    void positionAtNextBoundaryEdge(DartType& dart);
+
+  template<class TraitsType, class DartType>
+    bool convexBoundary(const DartType& dart);
+
+
+  // Utilities for Delaunay Triangulation
+  // ------------------------------------
+  template<class TraitsType, class DartType, class DartListType>
+    void optimizeDelaunay(DartListType& elist);
+
+  template <class TraitsType, class DartType, class DartListType>
+    void optimizeDelaunay(DartListType& elist, const typename DartListType::iterator end);
+
+  template<class TraitsType, class DartType>
+    bool swapTestDelaunay(const DartType& dart, bool cycling_check = false);
+
+  template<class TraitsType, class DartType>
+    void recSwapDelaunay(DartType& diagonal);
+
+  template<class TraitsType, class DartType, class ListType>
+    void swapEdgesAwayFromInteriorNode(DartType& dart, ListType& swapped_edges);
+
+  template<class TraitsType, class DartType, class ListType>
+    void swapEdgesAwayFromBoundaryNode(DartType& dart, ListType& swapped_edges);
+
+  template<class TraitsType, class DartType, class DartListType>
+    void swapEdgeInList(const typename DartListType::iterator& it, DartListType& elist);
+
+
+  // Constrained Triangulation
+  // -------------------------
+  template<class TraitsType, class DartType>
+    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 <class TraitsType, class DartType, class PointType>
+    bool insertNode(DartType& dart, PointType& point) {
+    
+    bool found = ttl::locateTriangle<TraitsType>(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<TraitsType>(d1);
+    }
+    
+    //if (!TraitsType::fixedEdge(d2) && !ttl::isBoundaryEdge(d2)) {
+    if (!ttl::isBoundaryEdge(d2)) {
+      d2.alpha2();
+      recSwapDelaunay<TraitsType>(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<TraitsType>(d3);
+    }
+    
+    return true;
+  }
+
+
+  //------------------------------------------------------------------------------------------------
+  // Private/Hidden function (might change later)
+  template <class TraitsType, class ForwardIterator, class DartType>
+    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<TraitsType>(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 <class TraitsType, class DartType>
+    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<TraitsType>(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 <class TraitsType, class DartType>
+    void removeNode(DartType& dart) {
+    
+    if (ttl::isBoundaryNode(dart))
+      ttl::removeBoundaryNode<TraitsType>(dart);
+    else
+      ttl::removeInteriorNode<TraitsType>(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 <class TraitsType, class DartType>
+    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<DartType> swapped_edges;
+    ttl::swapEdgesAwayFromBoundaryNode<TraitsType>(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<DartType>::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<DartType> DartListType;
+    ttl::optimizeDelaunay<TraitsType, DartType, DartListType>(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 <class TraitsType, class DartType>
+    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<DartType> swapped_edges;
+    ttl::swapEdgesAwayFromInteriorNode<TraitsType>(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<TraitsType, DartType>(swapped_edges);
+  }
+
+  //@} // End of Delaunay Triangulation Group
+
+
+  //------------------------------------------------------------------------------------------------
+  // -------------------------- Topological and Geometric Queries Group ---------------------------
+  //------------------------------------------------------------------------------------------------
+   
+  /** @name Topological and Geometric Queries */
+  //@{
+
+  //------------------------------------------------------------------------------------------------
+  // Private/Hidden function (might change later)
+  template <class TopologyElementType, class DartType>
+    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 <class TraitsType, class NodeType, class DartType>
+    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<TraitsType>(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 <class TraitsType, class PointType, class DartType>
+    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 <class TraitsType, class PointType, class DartType>
+    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<TraitsType>(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<TraitsType>(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 <class TraitsType, class PointType, class DartType>
+    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 <class TraitsType, class PointType, class DartType>
+    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 <class DartType>
+    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 <class DartType, class DartListType>
+    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 DartType>
+  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 <class DartType>
+    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 <class DartType>
+    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 <class DartType>
+    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 <class DartType>
+    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 <class DartType>
+    void getNeighborNodes(const DartType& dart, std::list<DartType>& 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 <class DartType, class DartListType>
+    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 <class DartType, class DartListType>
+    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 <class DartType>
+    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 <class DartType>
+    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 <class DartType>
+    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 <class TraitsType, class DartType>
+    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 <class TraitsType, class DartType>
+    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 <class DartType>
+    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 <class TraitsType, class DartType>
+    bool convexBoundary(const DartType& dart) {
+    
+    list<DartType> blist;
+    ttl::getBoundary(dart, blist);
+    
+    int no;
+    no = (int)blist.size();
+    typename list<DartType>::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 <class TraitsType, class DartType, class DartListType>
+    void optimizeDelaunay(DartListType& elist) {
+      optimizeDelaunay<TraitsType, DartType, DartListType>(elist, elist.end());
+  }
+
+
+  //------------------------------------------------------------------------------------------------
+  template <class TraitsType, class DartType, class DartListType>
+    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<TraitsType>(*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<TraitsType, DartType>(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 <class TraitsType, class DartType>
+#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 <class TraitsType, class DartType>
+    void recSwapDelaunay(DartType& diagonal) {
+
+    if (!ttl::swapTestDelaunay<TraitsType>(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<TraitsType>(oppEdge1);
+    if (!b2)
+      recSwapDelaunay<TraitsType>(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 <class TraitsType, class DartType, class ListType>
+    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<TraitsType>(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 <class TraitsType, class DartType, class ListType>
+    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<TraitsType>(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 <class TraitsType, class DartType, class DartListType>
+    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 <ttl/ttl_constr.h>
+
+#endif // _TTL_H_

include/ttl/ttl_constr.h

+/*
+ * 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
+ * <http://www.gnu.org/licenses/>.
+ *
+ * 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 <list>
+#include <cmath>
+
+
+// Debugging
+#ifdef DEBUG_TTL_CONSTR_PLOT
+  #include <fstream>
+  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 <class DartType>
+    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 <class TraitsType, class DartType>
+    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 <class TraitsType, class DartType>
+    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 <class TraitsType, class DartType, class ListType>
+    DartType findCrossingEdges(const DartType& dstart, const DartType& dend, ListType& elist) {
+    
+    const DartType my_start = getAtSmallestAngle<TraitsType>(dstart, dend);
+    DartType my_end   = getAtSmallestAngle<TraitsType>(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 <class TraitsType, class DartType>
+    void transformToConstraint(DartType& dstart, DartType& dend, std::list<DartType>& elist) {
+    
+    typename list<DartType>::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<TraitsType, DartType>(*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<TraitsType>(dstart, dend, *it, tmp) &&
+                            !crossesConstraint<TraitsType>(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<TraitsType, DartType>(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 <class TraitsType, class DartType>
+    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<DartType> elist;
+    DartType next_start = ttl_constr::findCrossingEdges<TraitsType>(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<TraitsType,DartType>(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<TraitsType>(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<DartType>::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<TraitsType, DartType>(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_

include/ttl/ttl_util.h

+/*
+ * 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
+ * <http://www.gnu.org/licenses/>.
+ *
+ * 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 <vector>
+#include <algorithm>
+
+
+#ifdef _MSC_VER
+#  if _MSC_VER < 1300
+#    include <minmax.h>
+#  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 <class real_type>
+    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 <class real_type>
+    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 <class real_type>
+    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 <class TTLtraits, class DartType>
+    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 <class TTLtraits>
+    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 <class TTLtraits>
+    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 <class TTLtraits>
+    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 <class TTLtraits, class DartType, class PointType2d>
+    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<double, double>.
+  *
+  *   \note 
+  *   - To deduce template argument for PointType the function must be
+  *    called with the syntax: \c createRandomData<MyPoint>(...) where \c MyPoint
+  *    is the actual point type.
+  */
+  template <class PointType>
+    std::vector<PointType*>* createRandomData(int noPoints, int seed=1) {
+    
+#ifdef _MSC_VER
+    srand(seed);
+#else
+    srand48((long int)seed);
+#endif
+    
+    double x, y;
+    std::vector<PointType*>* points = new std::vector<PointType*>(noPoints);
+    typename std::vector<PointType*>::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_

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()
 };
 

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

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

pcbnew/basepcbframe.cpp

@@ -54,6 +54,8 @@
 #include <trigo.h>
 #include <pcb_painter.h>
 #include <worksheet_viewitem.h>
+#include <ratsnest_data.h>
+#include <ratsnest_viewitem.h>
 
 #include <tool/tool_manager.h>
 #include <tool/tool_dispatcher.h>
@@ -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() )

pcbnew/class_board.cpp

@@ -42,6 +42,7 @@
 #include <pcb_netlist.h>
 #include <reporter.h>
 #include <base_units.h>
+#include <ratsnest_data.h>
 
 #include <pcbnew.h>
 #include <colors_selection.h>
@@ -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];

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.

pcbnew/pcbframe.cpp

@@ -59,6 +59,13 @@
 #include <view/view.h>
 #include <painter.h>
 
+#include <class_track.h>
+#include <class_board.h>
+#include <class_module.h>
+#include <worksheet_viewitem.h>
+#include <ratsnest_data.h>
+#include <ratsnest_viewitem.h>
+
 #include <tool/tool_manager.h>
 #include <tool/tool_dispatcher.h>
 
@@ -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 )

pcbnew/ratsnest_data.cpp

+/*
+ * This program source code file is part of KICAD, a free EDA CAD application.
+ *
+ * Copyright (C) 2013 CERN
+ * @author Maciej Suminski <maciej.suminski@cern.ch>
+ *
+ * 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 <ratsnest_data.h>
+
+#include <class_board.h>
+#include <class_module.h>
+#include <class_pad.h>
+#include <class_track.h>
+#include <class_zone.h>
+
+#include <boost/foreach.hpp>
+#include <boost/range/adaptor/map.hpp>
+#include <boost/scoped_ptr.hpp>
+#include <boost/make_shared.hpp>
+#include <boost/bind.hpp>
+
+#include <cassert>
+#include <algorithm>
+#include <limits>
+
+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<RN_EDGE_PTR>* kruskalMST( RN_LINKS::RN_EDGE_LIST& aEdges,
+                                      const std::vector<RN_NODE_PTR>& aNodes )
+{
+    unsigned int nodeNumber = aNodes.size();
+    unsigned int mstExpectedSize = nodeNumber - 1;
+    unsigned int mstSize = 0;
+
+    // The output
+    std::vector<RN_EDGE_PTR>* mst = new std::vector<RN_EDGE_PTR>;
+    mst->reserve( mstExpectedSize );
+
+    // Set tags for marking cycles
+    boost::unordered_map<RN_NODE_PTR, int> 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<std::list<int> > 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<int>::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<RN_NODE_PTR> 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<RN_NODE_PTR> 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<RN_NODE>( 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<RN_EDGE_MST>( 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<RN_EDGE_PTR>( 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<RN_EDGE_MST>( *boardNodes.begin(), *last ) );
+        }
+
+        return;
+    }
+    else if( boardNodes.size() == 1 )   // This case is even simpler
+    {
+        m_rnEdges.reset( new std::vector<RN_EDGE_PTR>( 0 ) );
+
+        return;
+    }
+
+    // Move and sort (sorting speeds up) all nodes to a vector for the Delaunay triangulation
+    std::vector<RN_NODE_PTR> 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<RN_LINKS::RN_EDGE_LIST> 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<RN_EDGE_PTR>::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<CPolyPt>& 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<RN_POLY>& 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<RN_EDGE_PTR>& 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<unsigned int>::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<unsigned int>::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_NODE_PTR> RN_NET::GetClosestNodes( const RN_NODE_PTR& aNode, int aNumber ) const
+{
+    std::list<RN_NODE_PTR> 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_NODE_PTR> RN_NET::GetClosestNodes( const RN_NODE_PTR& aNode,
+                                                RN_NODE_FILTER aFilter, int aNumber ) const
+{
+    std::list<RN_NODE_PTR> 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<size_t>( aNumber ), nodes.size() ) );
+
+    return closest;
+}
+
+
+std::list<RN_NODE_PTR> RN_NET::GetNodes( const BOARD_CONNECTED_ITEM* aItem ) const
+{
+    std::list<RN_NODE_PTR> nodes;
+
+    switch( aItem->Type() )
+    {
+    case PCB_PAD_T:
+    {
+        const D_PAD* pad = static_cast<const D_PAD*>( aItem );
+        nodes.push_back( m_pads.at( pad ) );
+    }
+    break;
+
+    case PCB_VIA_T:
+    {
+        const SEGVIA* via = static_cast<const SEGVIA*>( aItem );
+        nodes.push_back( m_vias.at( via ) );
+    }
+    break;
+
+    case PCB_TRACE_T:
+    {
+        const TRACK* track = static_cast<const TRACK*>( 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<RN_NODE_PTR> nodes = m_nets[net].GetNodes( aItem );
+    std::list<RN_NODE_PTR>::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<RN_EDGE_PTR>& 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<RN_POLY>& polygons, m_zonePolygons | boost::adaptors::map_values )
+    {
+        RN_LINKS::RN_NODE_SET::iterator point, pointEnd;
+        std::deque<RN_POLY>::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<const D_PAD*>( aItem );
+        m_nets[net].RemoveItem( pad );
+        m_nets[net].AddItem( pad );
+    }
+    break;
+
+    case PCB_TRACE_T:
+    {
+        const TRACK* track = static_cast<const TRACK*>( aItem );
+        m_nets[net].RemoveItem( track );
+        m_nets[net].AddItem( track );
+    }
+    break;
+
+    case PCB_VIA_T:
+    {
+        const SEGVIA* via = static_cast<const SEGVIA*>( aItem );
+        m_nets[net].RemoveItem( via );
+        m_nets[net].AddItem( via );
+    }
+    break;
+
+    case PCB_ZONE_AREA_T:
+    {
+        const ZONE_CONTAINER* zone = static_cast<const ZONE_CONTAINER*>( 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<SEGVIA*>( 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();
+}

pcbnew/ratsnest_data.h

+/*
+ * This program source code file is part of KICAD, a free EDA CAD application.
+ *
+ * Copyright (C) 2013 CERN
+ * @author Maciej Suminski <maciej.suminski@cern.ch>
+ *
+ * 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 <ttl/halfedge/hetriang.h>
+#include <ttl/halfedge/hetraits.h>
+
+#include <math/box2.h>
+
+#include <boost/unordered_set.hpp>
+#include <boost/unordered_map.hpp>
+
+class BOARD;
+class BOARD_ITEM;
+class BOARD_CONNECTED_ITEM;
+class MODULE;
+class D_PAD;
+class SEGVIA;
+class TRACK;
+class ZONE_CONTAINER;
+class CPolyPt;
+
+// Preserve KiCad coding style policy
+typedef hed::Node RN_NODE;
+typedef hed::NodePtr RN_NODE_PTR;
+typedef hed::Edge RN_EDGE;
+typedef hed::EdgePtr RN_EDGE_PTR;
+typedef hed::EdgeMST RN_EDGE_MST;
+typedef boost::shared_ptr<hed::EdgeMST> RN_EDGE_MST_PTR;
+typedef hed::Triangulation TRIANGULATOR;
+
+///> General interface for filtering out nodes in search functions.
+struct RN_NODE_FILTER : public std::unary_function<const RN_NODE_PTR&, bool>
+{
+    virtual ~RN_NODE_FILTER() {}
+    virtual bool operator()( const RN_NODE_PTR& aNode )
+    {
+        return true;        // By default everything passes
+    }
+};
+
+///> Filters out nodes that have the flag set.
+struct WITHOUT_FLAG : public RN_NODE_FILTER
+{
+    bool operator()( const RN_NODE_PTR& aNode )
+    {
+        return !aNode->GetFlag();
+    }
+};
+
+///> Functor comparing if two nodes are equal by their coordinates. It is required to make set of
+///> shared pointers work properly.
+struct RN_NODE_COMPARE : std::binary_function<RN_NODE_PTR, RN_NODE_PTR, bool>
+{
+    bool operator()( const RN_NODE_PTR& aNode1, const RN_NODE_PTR& aNode2 ) const
+    {
+        return ( aNode1->GetX() == aNode2->GetX() && aNode1->GetY() == aNode2->GetY() );
+    }
+};
+
+///> Functor calculating hash for a given node. It is required to make set of shared pointers
+///> work properly.
+struct RN_NODE_HASH : std::unary_function<RN_NODE_PTR, std::size_t>
+{
+    std::size_t operator()( const RN_NODE_PTR& aNode ) const
+    {
+        std::size_t hash = 2166136261u;
+
+        hash ^= aNode->GetX();
+        hash *= 16777619;
+        hash ^= aNode->GetY();
+
+        return hash;
+    }
+};
+
+
+/**
+ * Class RN_LINKS
+ * Manages data describing nodes and connections for a given net.
+ */
+class RN_LINKS
+{
+public:
+    // Helper typedefs
+    typedef boost::unordered_set<RN_NODE_PTR, RN_NODE_HASH, RN_NODE_COMPARE> RN_NODE_SET;
+    typedef std::list<RN_EDGE_PTR> RN_EDGE_LIST;
+
+    /**
+     * Function AddNode()
+     * Adds a node with given coordinates and returns pointer to the newly added node. If the node
+     * existed before, only appropriate pointer is returned.
+     * @param aX is the x coordinate of a node.
+     * @param aY is the y coordinate of a node.
+     * @return Pointer to the node with given coordinates.
+     */
+    const RN_NODE_PTR& AddNode( int aX, int aY );
+
+    /**
+     * Function RemoveNode()
+     * Removes a node described by a given node pointer.
+     * @param aNode is a pointer to node to be removed.
+     */
+    void RemoveNode( const RN_NODE_PTR& aNode );
+
+    /**
+     * Function GetNodes()
+     * Returns the set of currently used nodes.
+     * @return The set of currently used nodes.
+     */
+    const RN_NODE_SET& GetNodes() const
+    {
+        return m_nodes;
+    }
+
+    /**
+     * Function AddConnection()
+     * Adds a connection between two nodes and of given distance. Edges with distance equal 0 are
+     * considered to be existing connections. Distance different than 0 means that the connection
+     * is missing.
+     * @param aNode1 is the origin node of a new connection.
+     * @param aNode2 is the end node of a new connection.
+     * @param aDistance is the distance of the connection (0 means that nodes are actually
+     * connected, >0 means a missing connection).
+     */
+    const RN_EDGE_PTR& AddConnection( const RN_NODE_PTR& aNode1, const RN_NODE_PTR& aNode2,
+                                      unsigned int aDistance = 0 );
+
+    /**
+     * Function RemoveConnection()
+     * Removes a connection described by a given edge pointer.
+     * @param aEdge is a pointer to edge to be removed.
+     */
+    void RemoveConnection( const RN_EDGE_PTR& aEdge );
+
+    /**
+     * Function GetConnections()
+     * Returns the list of edges that currently connect nodes.
+     * @return the list of edges that currently connect nodes.
+     */
+    const RN_EDGE_LIST& GetConnections() const
+    {
+        return m_edges;
+    }
+
+protected:
+    ///> Set of nodes that are used are expected to be connected together.
+    RN_NODE_SET m_nodes;
+
+    ///> List of edges that currently connect nodes.
+    RN_EDGE_LIST m_edges;
+};
+
+
+/**
+ * Class RN_POLY
+ * Describes a single subpolygon (ZONE_CONTAINER is supposed to contain one or more of those) and
+ * performs fast point-inside-polygon test.
+ */
+class RN_POLY
+{
+public:
+    RN_POLY( const CPolyPt* aBegin, const CPolyPt* aEnd, const ZONE_CONTAINER* aParent,
+          RN_LINKS& aConnections, const BOX2I& aBBox );
+
+    /**
+     * Function GetNode()
+     * Returns node representing a polygon (it has the same coordinates as the first point of its
+     * bounding polyline.
+     */
+    const RN_NODE_PTR& GetNode() const
+    {
+        return m_node;
+    }
+
+
+    /**
+     * Function GetParent()
+     * Returns pointer to zone that is the owner of subpolygon.
+     * @return Pointer to zone that is the owner of subpolygon.
+     */
+    const ZONE_CONTAINER* GetParent() const
+    {
+        return m_parent;
+    }
+
+    /**
+     * Function HitTest()
+     * Tests if selected node is located within polygon boundaries.
+     * @param aNode is a node to be checked.
+     * @return True is the node is located within polygon boundaries.
+     */
+    bool HitTest( const RN_NODE_PTR& aNode ) const;
+
+private:
+    ///> Owner of this subpolygon.
+    const ZONE_CONTAINER* m_parent;
+
+    ///> Pointer to the first point of polyline bounding the polygon.
+    const CPolyPt* m_begin;
+
+    ///> Pointer to the last point of polyline bounding the polygon.
+    const CPolyPt* m_end;
+
+    ///> Bounding box of the polygon.
+    BOX2I m_bbox;
+
+    ///> Node representing a polygon (it has the same coordinates as the first point of its
+    ///> bounding polyline.
+    RN_NODE_PTR m_node;
+
+    friend bool sortArea( const RN_POLY& aP1, const RN_POLY& aP2 );
+};
+
+
+/**
+ * Class RN_NET
+ * Describes ratsnest for a single net.
+ */
+class RN_NET
+{
+public:
+    ///> Default constructor.
+    RN_NET() : m_dirty( true ), m_visible( true )
+    {}
+
+    /**
+     * Function SetVisible()
+     * Sets state of the visibility flag.
+     * @param aEnabled is new state. True if ratsnest for a given net is meant to be displayed,
+     * false otherwise.
+     */
+    void SetVisible( bool aEnabled )
+    {
+        m_visible = aEnabled;
+    }
+
+    /**
+     * Function IsVisible()
+     * Returns the visibility flag state.
+     * @return True if ratsnest for given net is set as visible, false otherwise,
+     */
+    bool IsVisible() const
+    {
+        return m_visible;
+    }
+
+    /**
+     * Function MarkDirty()
+     * Marks ratsnest for given net as 'dirty', ie. requiring recomputation.
+     */
+    void MarkDirty()
+    {
+        m_dirty = true;
+    }
+
+    /**
+     * Function IsDirty()
+     * Returns state of the 'dirty' flag, indicating that ratsnest for a given net is invalid
+     * and requires an update.
+     * @return True if ratsnest requires recomputation, false otherwise.
+     */
+    bool IsDirty() const
+    {
+        return m_dirty;
+    }
+
+    /**
+     * Function GetUnconnected()
+     * Returns pointer to a vector of edges that makes ratsnest for a given net.
+     * @return Pointer to a vector of edges that makes ratsnest for a given net.
+     */
+    const std::vector<RN_EDGE_PTR>* GetUnconnected() const
+    {
+        return m_rnEdges.get();
+    }
+
+    /**
+     * Function Update()
+     * Recomputes ratsnest for a net.
+     */
+    void Update();
+
+    /**
+     * Function AddItem()
+     * Adds an appropriate node associated with selected pad, so it is
+     * taken into account during ratsnest computations.
+     * @param aPad is a pad for which node is added.
+     */
+    void AddItem( const D_PAD* aPad );
+
+    /**
+     * Function AddItem()
+     * Adds an appropriate node associated with selected via, so it is
+     * taken into account during ratsnest computations.
+     * @param aVia is a via for which node is added.
+     */
+    void AddItem( const SEGVIA* aVia );
+
+    /**
+     * Function AddItem()
+     * Adds appropriate nodes and edges associated with selected track, so they are
+     * taken into account during ratsnest computations.
+     * @param aTrack is a track for which nodes and edges are added.
+     */
+    void AddItem( const TRACK* aTrack );
+
+    /**
+     * Function AddItem()
+     * Processes zone to split it into subpolygons and adds appropriate nodes for them, so they are
+     * taken into account during ratsnest computations.
+     * @param aZone is a zone to be processed.
+     */
+    void AddItem( const ZONE_CONTAINER* aZone );
+
+    /**
+     * Function RemoveItem()
+     * Removes all nodes and edges associated with selected pad, so they are not
+     * taken into account during ratsnest computations anymore.
+     * @param aPad is a pad for which nodes and edges are removed.
+     */
+    void RemoveItem( const D_PAD* aPad );
+
+    /**
+     * Function RemoveItem()
+     * Removes all nodes and edges associated with selected via, so they are not
+     * taken into account during ratsnest computations anymore.
+     * @param aVia is a via for which nodes and edges are removed.
+     */
+    void RemoveItem( const SEGVIA* aVia );
+
+    /**
+     * Function RemoveItem()
+     * Removes all nodes and edges associated with selected track, so they are not
+     * taken into account during ratsnest computations anymore.
+     * @param aTrack is a track for which nodes and edges are removed.
+     */
+    void RemoveItem( const TRACK* aTrack );
+
+    /**
+     * Function RemoveItem()
+     * Removes all nodes and edges associated with selected zone, so they are not
+     * taken into account during ratsnest computations anymore.
+     * @param aZone is a zone for which nodes and edges are removed.
+     */
+    void RemoveItem( const ZONE_CONTAINER* aZone );
+
+    /**
+     * Function GetNodes()
+     * Returns list of nodes that are associated with a given item.
+     * @param aItem is an item for which the list is generated.
+     * @return List of associated nodes.
+     */
+    std::list<RN_NODE_PTR> GetNodes( const BOARD_CONNECTED_ITEM* aItem ) const;
+
+    /**
+     * Function GetClosestNode()
+     * Returns a single node that lies in the shortest distance from a specific node.
+     * @param aNode is the node for which the closest node is searched.
+     */
+    const RN_NODE_PTR GetClosestNode( const RN_NODE_PTR& aNode ) const;
+
+    /**
+     * Function GetClosestNode()
+     * Returns a single node that lies in the shortest distance from a specific node and meets
+     * selected filter criterion..
+     * @param aNode is the node for which the closest node is searched.
+     * @param aFilter is a functor that filters nodes.
+     */
+    const RN_NODE_PTR GetClosestNode( const RN_NODE_PTR& aNode, RN_NODE_FILTER aFilter ) const;
+
+    /**
+     * Function GetClosestNodes()
+     * Returns list of nodes sorted by the distance from a specific node.
+     * @param aNode is the node for which the closest nodes are searched.
+     * @param aNumber is asked number of returned nodes. If it is negative then all nodes that
+     * belong to the same net are returned. If asked number is greater than number of possible
+     * nodes then the size of list is limited to number of possible nodes.
+     */
+    std::list<RN_NODE_PTR> GetClosestNodes( const RN_NODE_PTR& aNode, int aNumber = -1 ) const;
+
+    /**
+     * Function GetClosestNodes()
+     * Returns filtered list of nodes sorted by the distance from a specific node.
+     * @param aNode is the node for which the closest nodes are searched.
+     * @param aFilter is a functor that filters nodes.
+     * @param aNumber is asked number of returned nodes. If it is negative then all nodes that
+     * belong to the same net are returned. If asked number is greater than number of possible
+     * nodes then the size of list is limited to number of possible nodes.
+     */
+    std::list<RN_NODE_PTR> GetClosestNodes( const RN_NODE_PTR& aNode,
+                                            RN_NODE_FILTER aFilter, int aNumber = -1 ) const;
+
+    /**
+     * Function GetEdges()
+     * Returns pointer to the vector of edges that makes ratsnest for a given net.
+     * @return Pointer to the vector of edges that makes ratsnest for a given net
+     */
+    const std::vector<RN_EDGE_PTR>* GetEdges() const
+    {
+        return m_rnEdges.get();
+    }
+
+    /**
+     * Function AddSimpleNode()
+     * Changes drawing mode for a node to simple (ie. one ratsnest line per node).
+     * @param aNode is a node that changes its drawing mode..
+     */
+    void AddSimpleNode( RN_NODE_PTR& aNode )
+    {
+        m_simpleNodes.push_back( aNode );
+        aNode->SetFlag( true );
+    }
+
+    /**
+     * Function GetSimpleNodes()
+     * Returns list of nodes for which ratsnest is drawn in simple mode (ie. one
+     * ratsnest line per node).
+     * @return list of nodes for which ratsnest is drawn in simple mode.
+     */
+    const std::deque<RN_NODE_PTR>& GetSimpleNodes() const
+    {
+        return m_simpleNodes;
+    }
+
+    /**
+     * Function ClearSimple()
+     * Removes all nodes and edges that are used for displaying ratsnest in simple mode.
+     */
+    void ClearSimple();
+
+protected:
+    ///> Validates edge, ie. modifies source and target nodes for an edge
+    ///> to make sure that they are not ones with the flag set.
+    void validateEdge( RN_EDGE_PTR& aEdge );
+
+    ///> Removes all ratsnest edges for a given node.
+    void clearNode( const RN_NODE_PTR& aNode );
+
+    ///> Adds appropriate edges for nodes that are connected by zones.
+    void processZones();
+
+    ///> Recomputes ratsnset from scratch.
+    void compute();
+
+    ////> Stores information about connections for a given net.
+    RN_LINKS m_links;
+
+    ///> Vector of edges that makes ratsnest for a given net.
+    boost::shared_ptr< std::vector<RN_EDGE_PTR> > m_rnEdges;
+
+    ///> List of nodes for which ratsnest is drawn in simple mode.
+    std::deque<RN_NODE_PTR> m_simpleNodes;
+
+    ///> Flag indicating necessity of recalculation of ratsnest for a net.
+    bool m_dirty;
+
+    ///> Map that associates nodes in the ratsnest model to respective nodes.
+    boost::unordered_map<const D_PAD*, RN_NODE_PTR> m_pads;
+
+    ///> Map that associates nodes in the ratsnest model to respective vias.
+    boost::unordered_map<const SEGVIA*, RN_NODE_PTR> m_vias;
+
+    ///> Map that associates edges in the ratsnest model to respective tracks.
+    boost::unordered_map<const TRACK*, RN_EDGE_PTR> m_tracks;
+
+    ///> Map that associates groups of subpolygons in the ratsnest model to their respective zones.
+    boost::unordered_map<const ZONE_CONTAINER*, std::deque<RN_POLY> > m_zonePolygons;
+
+    ///> Map that associates groups of edges in the ratsnest model to their respective zones.
+    boost::unordered_map<const ZONE_CONTAINER*, std::deque<RN_EDGE_PTR> > m_zoneConnections;
+
+    ///> Visibility flag.
+    bool m_visible;
+};
+
+
+/**
+ * Class RN_DATA
+ *
+ * Stores information about unconnected items for the whole PCB.
+ */
+class RN_DATA
+{
+public:
+    /**
+     * Default constructor
+     * @param aBoard is the board to be processed in order to look for unconnected items.
+     */
+    RN_DATA( const BOARD* aBoard ) : m_board( aBoard ) {}
+
+    /**
+     * Function UpdateItem()
+     * Updates ratsnest data for an item.
+     * @param aItem is an item to be updated.
+     */
+    void Update( const BOARD_CONNECTED_ITEM* aItem );
+
+    /**
+     * Function UpdateItem()
+     * Updates ratsnest data for a module.
+     * @param aItem is a module to be updated.
+     */
+    void Update( const MODULE* aModule );
+
+    /**
+     * Function AddSimple()
+     * Sets an item to be drawn in simple mode (ie. one line per node, instead of full ratsnest).
+     * It is used for drawing temporary ratsnest, eg. while moving an item.
+     * @param aItem is an item to be drawn in simple node.
+     */
+    void AddSimple( const BOARD_CONNECTED_ITEM* aItem );
+
+    /**
+     * Function AddSimple()
+     * Sets a module to be drawn in simple mode (ie. one line per node, instead of full ratsnest).
+     * It is used for drawing temporary ratsnest, eg. while moving a module.
+     * @param aModule is a module to be drawn in simple node.
+     */
+    void AddSimple( const MODULE* aModule );
+
+    /**
+     * Function ClearSimple()
+     * Clears the list of nodes for which ratsnest is drawn in simple mode (one line per node).
+     */
+    void ClearSimple();
+
+    /**
+     * Function ProcessBoard()
+     * Prepares data for computing (computes a list of current nodes and connections). It is
+     * required to run only once after loading a board.
+     */
+    void ProcessBoard();
+
+    /**
+     * Function Recalculate()
+     * Recomputes ratsnest for selected net number or all nets that need updating.
+     * @param aNet is a net number. If it is negative, all nets that need updating are recomputed.
+     */
+    void Recalculate( int aNet = -1 );
+
+    /**
+     * Function GetNets()
+     * Returns ratsnest grouped by net numbers.
+     * @return Vector of ratsnest grouped by net numbers.
+     */
+    const std::vector<RN_NET>& GetNets() const
+    {
+        return m_nets;
+    }
+
+protected:
+    /**
+     * Function updateNet()
+     * Recomputes ratsnest for a single net.
+     * @param aNetCode is the net number to be recomputed.
+     */
+    void updateNet( int aNetCode );
+
+    ///> Board to be processed.
+    const BOARD* m_board;
+
+    ///> Stores information about ratsnest grouped by net numbers.
+    std::vector<RN_NET> m_nets;
+};
+
+#endif /* RATSNEST_DATA_H */

pcbnew/ratsnest_viewitem.cpp

+/*
+ * This program source code file is part of KICAD, a free EDA CAD application.
+ *
+ * Copyright (C) 2013 CERN
+ * @author Maciej Suminski <maciej.suminski@cern.ch>
+ *
+ * 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_viewitem.cpp
+ * @brief Class that draws missing connections on a PCB.
+ */
+
+#include <ratsnest_viewitem.h>
+#include <ratsnest_data.h>
+#include <gal/graphics_abstraction_layer.h>
+#include <layers_id_colors_and_visibility.h>
+
+#include <boost/foreach.hpp>
+
+using namespace KIGFX;
+
+RATSNEST_VIEWITEM::RATSNEST_VIEWITEM( RN_DATA* aData ) :
+        EDA_ITEM( NOT_USED ), m_data( aData )
+{
+}
+
+
+const BOX2I RATSNEST_VIEWITEM::ViewBBox() const
+{
+    // Make it always visible
+    BOX2I bbox;
+    bbox.SetMaximum();
+
+    return bbox;
+}
+
+
+void RATSNEST_VIEWITEM::ViewDraw( int aLayer, GAL* aGal ) const
+{
+    aGal->SetIsStroke( true );
+    aGal->SetIsFill( false );
+    aGal->SetLineWidth( 1.0 );
+    aGal->SetStrokeColor( COLOR4D( 1.0, 1.0, 1.0, 0.4 ) );
+
+    WITHOUT_FLAG without_flag;
+
+    // Draw the temporary ratsnest
+    BOOST_FOREACH( const RN_NET& net, m_data->GetNets() )
+    {
+        if( !net.IsVisible() )
+            continue;
+
+        // Avoid duplicate destinations for ratsnest lines by storing already used nodes
+        boost::unordered_set<RN_NODE_PTR> usedDestinations;
+
+        // Draw the "dynamic" ratsnest (ie. for objects that may be currently being moved)
+        BOOST_FOREACH( const RN_NODE_PTR& node, net.GetSimpleNodes() )
+        {
+            RN_NODE_PTR dest = net.GetClosestNode( node, without_flag );
+
+            if( dest && usedDestinations.find( dest ) == usedDestinations.end() )
+            {
+                VECTOR2D origin( node->GetX(), node->GetY() );
+                VECTOR2D end( dest->GetX(), dest->GetY() );
+
+                aGal->DrawLine( origin, end );
+                usedDestinations.insert( dest );
+            }
+        }
+
+        // Draw the "static" ratsnest
+        const std::vector<RN_EDGE_PTR>* edges = net.GetUnconnected();
+        if( edges == NULL )
+            continue;
+
+        BOOST_FOREACH( const RN_EDGE_PTR& edge, *edges )
+        {
+            const RN_NODE_PTR& sourceNode = edge->getSourceNode();
+            const RN_NODE_PTR& targetNode = edge->getTargetNode();
+            VECTOR2D source( sourceNode->GetX(), sourceNode->GetY()  );
+            VECTOR2D target( targetNode->GetX(), targetNode->GetY() );
+
+            aGal->DrawLine( source, target );
+        }
+    }
+}
+
+
+void RATSNEST_VIEWITEM::ViewGetLayers( int aLayers[], int& aCount ) const
+{
+    aCount = 1;
+    aLayers[0] = ITEM_GAL_LAYER( RATSNEST_VISIBLE );
+}

pcbnew/ratsnest_viewitem.h

+/*
+ * This program source code file is part of KICAD, a free EDA CAD application.
+ *
+ * Copyright (C) 2013 CERN
+ * @author Maciej Suminski <maciej.suminski@cern.ch>
+ *
+ * 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_viewitem.h
+ * @brief Class that draws missing connections on a PCB.
+ */
+
+#ifndef RATSNEST_VIEWITEM_H
+#define RATSNEST_VIEWITEM_H
+
+#include <base_struct.h>
+#include <math/vector2d.h>
+
+class GAL;
+class RN_DATA;
+
+namespace KIGFX
+{
+class RATSNEST_VIEWITEM : public EDA_ITEM
+{
+public:
+    RATSNEST_VIEWITEM( RN_DATA* aData );
+
+    /// @copydoc VIEW_ITEM::ViewBBox()
+    const BOX2I ViewBBox() const;
+
+    /// @copydoc VIEW_ITEM::ViewDraw()
+    void ViewDraw( int aLayer, GAL* aGal ) const;
+
+    /// @copydoc VIEW_ITEM::ViewGetLayers()
+    void ViewGetLayers( int aLayers[], int& aCount ) const;
+
+    /// @copydoc EDA_ITEM::Show()
+    void Show( int x, std::ostream& st ) const
+    {
+    }
+
+protected:
+    ///> Object containing ratsnest data.
+    RN_DATA* m_data;
+};
+
+}   // namespace KIGFX
+
+#endif /* RATSNEST_VIEWITEM_H */

pcbnew/tools/move_tool.cpp

@@ -26,8 +26,11 @@
 #include <class_module.h>
 #include <tool/tool_manager.h>
 #include <view/view_controls.h>
+#include <ratsnest_data.h>
 #include <confirm.h>
 
+#include <boost/foreach.hpp>
+
 #include "common_actions.h"
 #include "selection_tool.h"
 #include "move_tool.h"
@@ -41,11 +44,6 @@ MOVE_TOOL::MOVE_TOOL() :
 }
 
 
-MOVE_TOOL::~MOVE_TOOL()
-{
-}
-
-
 void MOVE_TOOL::Reset()
 {
     // The tool launches upon reception of action event ("pcbnew.InteractiveMove")
@@ -111,11 +109,13 @@ int MOVE_TOOL::Main( TOOL_EVENT& aEvent )
             {
                 m_state.Rotate( cursorPos, 900.0 );
                 selection.group->ViewUpdate( VIEW_ITEM::GEOMETRY );
+                updateRatsnest( true );
             }
             else if( evt->IsAction( &COMMON_ACTIONS::flip ) )        // got flip event?
             {
                 m_state.Flip( cursorPos );
                 selection.group->ViewUpdate( VIEW_ITEM::GEOMETRY );
+                updateRatsnest( true );
             }
         }
 
@@ -126,6 +126,8 @@ int MOVE_TOOL::Main( TOOL_EVENT& aEvent )
                 // Drag items to the current cursor position
                 VECTOR2D movement = ( evt->Position() - dragPosition );
                 m_state.Move( movement );
+
+                updateRatsnest( true );
             }
             else
             {
@@ -153,17 +155,50 @@ int MOVE_TOOL::Main( TOOL_EVENT& aEvent )
         // Modifications has to be rollbacked, so restore the previous state of items
         selection.group->ItemsViewUpdate( VIEW_ITEM::APPEARANCE );
         m_state.RestoreAll();
+
+        updateRatsnest( false );
     }
     else
     {
         // Changes are applied, so update the items
         selection.group->ItemsViewUpdate( m_state.GetUpdateFlag() );
         m_state.Apply();
+
     }
 
+    RN_DATA* ratsnest = static_cast<BOARD*>( m_toolMgr->GetModel() )->GetRatsnest();
+    ratsnest->Recalculate();
+
     controls->ShowCursor( false );
     controls->SetSnapping( false );
     controls->SetAutoPan( false );
 
     return 0;
 }
+
+
+void MOVE_TOOL::updateRatsnest( bool aRedraw )
+{
+    const SELECTION_TOOL::SELECTION& selection = m_selectionTool->GetSelection();
+    RN_DATA* ratsnest = static_cast<BOARD*>( m_toolMgr->GetModel() )->GetRatsnest();
+
+    ratsnest->ClearSimple();
+    BOOST_FOREACH( BOARD_ITEM* item, selection.items )
+    {
+        if( item->Type() == PCB_PAD_T || item->Type() == PCB_TRACE_T ||
+                item->Type() == PCB_VIA_T || item->Type() == PCB_ZONE_AREA_T )
+        {
+            ratsnest->Update( static_cast<BOARD_CONNECTED_ITEM*>( item ) );
+
+            if( aRedraw )
+                ratsnest->AddSimple( static_cast<BOARD_CONNECTED_ITEM*>( item ) );
+        }
+        else if( item->Type() == PCB_MODULE_T )
+        {
+            ratsnest->Update( static_cast<MODULE*>( item ) );
+
+            if( aRedraw )
+                ratsnest->AddSimple( static_cast<MODULE*>( item ) );
+        }
+    }
+}

pcbnew/tools/move_tool.h

@@ -49,7 +49,6 @@ class MOVE_TOOL : public TOOL_INTERACTIVE
 {
 public:
     MOVE_TOOL();
-    ~MOVE_TOOL();
 
     /// @copydoc TOOL_INTERACTIVE::Reset()
     void Reset();
@@ -65,6 +64,8 @@ public:
     int Main( TOOL_EVENT& aEvent );
 
 private:
+    void updateRatsnest( bool aRedraw );
+
     /// Saves the state of items and allows to restore them
     ITEM_STATE m_state;