diff options
Diffstat (limited to 'include/ttl/ttl.h')
| -rw-r--r-- | include/ttl/ttl.h | 1901 |
1 files changed, 1901 insertions, 0 deletions
diff --git a/include/ttl/ttl.h b/include/ttl/ttl.h new file mode 100644 index 0000000..bbcf86a --- /dev/null +++ b/include/ttl/ttl.h @@ -0,0 +1,1901 @@ +/* + * 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* aMessage ) +{ + cout << "\n!!! ERROR: " << aMessage << " !!!\n" << endl; + exit(-1); +} +#endif + +// 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 DART_TYPE and \e TRAITS_TYPE 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 TRAITS_TYPE to +* support that specific function.\n +* Functionalty required in \e DART_TYPE 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 TRAITS_TYPE 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 "TRAITS_TYPE::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 get_0_orbit_interior and 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 +{ +class TRIANGULATION_HELPER +{ +#ifndef DOXYGEN_SHOULD_SKIP_THIS + +public: + TRIANGULATION_HELPER( hed::TRIANGULATION& aTriang ) : + m_triangulation( aTriang ) + { + } + + // Delaunay Triangulation + template <class TRAITS_TYPE, class DART_TYPE, class POINT_TYPE> + bool InsertNode( DART_TYPE& aDart, POINT_TYPE& aPoint ); + + template <class TRAITS_TYPE, class DART_TYPE> + void RemoveRectangularBoundary( DART_TYPE& aDart ); + + template <class TRAITS_TYPE, class DART_TYPE> + void RemoveNode( DART_TYPE& aDart ); + + template <class TRAITS_TYPE, class DART_TYPE> + void RemoveBoundaryNode( DART_TYPE& aDart ); + + template <class TRAITS_TYPE, class DART_TYPE> + void RemoveInteriorNode( DART_TYPE& aDart ); + + // Topological and Geometric Queries + // --------------------------------- + template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE> + static bool LocateFaceSimplest( const POINT_TYPE& aPoint, DART_TYPE& aDart ); + + template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE> + static bool LocateTriangle( const POINT_TYPE& aPoint, DART_TYPE& aDart ); + + template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE> + static bool InTriangle( const POINT_TYPE& aPoint, const DART_TYPE& aDart ); + + template <class DART_TYPE, class DART_LIST_TYPE> + static void GetBoundary( const DART_TYPE& aDart, DART_LIST_TYPE& aBoundary ); + + template <class DART_TYPE> + static bool IsBoundaryEdge( const DART_TYPE& aDart ); + + template <class DART_TYPE> + static bool IsBoundaryFace( const DART_TYPE& aDart ); + + template <class DART_TYPE> + static bool IsBoundaryNode( const DART_TYPE& aDart ); + + template <class DART_TYPE> + static int GetDegreeOfNode( const DART_TYPE& aDart ); + + template <class DART_TYPE, class DART_LIST_TYPE> + static void Get0OrbitInterior( const DART_TYPE& aDart, DART_LIST_TYPE& aOrbit ); + + template <class DART_TYPE, class DART_LIST_TYPE> + static void Get0OrbitBoundary( const DART_TYPE& aDart, DART_LIST_TYPE& aOrbit ); + + template <class DART_TYPE> + static bool Same0Orbit( const DART_TYPE& aD1, const DART_TYPE& aD2 ); + + template <class DART_TYPE> + static bool Same1Orbit( const DART_TYPE& aD1, const DART_TYPE& aD2 ); + + template <class DART_TYPE> + static bool Same2Orbit( const DART_TYPE& aD1, const DART_TYPE& aD2 ); + + template <class TRAITS_TYPE, class DART_TYPE> + static bool SwappableEdge( const DART_TYPE& aDart, bool aAllowDegeneracy = false ); + + template <class DART_TYPE> + static void PositionAtNextBoundaryEdge( DART_TYPE& aDart ); + + template <class TRAITS_TYPE, class DART_TYPE> + static bool ConvexBoundary( const DART_TYPE& aDart ); + + // Utilities for Delaunay Triangulation + // ------------------------------------ + template <class TRAITS_TYPE, class DART_TYPE, class DART_LIST_TYPE> + void OptimizeDelaunay( DART_LIST_TYPE& aElist ); + + template <class TRAITS_TYPE, class DART_TYPE, class DART_LIST_TYPE> + void OptimizeDelaunay( DART_LIST_TYPE& aElist, const typename DART_LIST_TYPE::iterator aEnd ); + + template <class TRAITS_TYPE, class DART_TYPE> + bool SwapTestDelaunay( const DART_TYPE& aDart, bool aCyclingCheck = false ) const; + + template <class TRAITS_TYPE, class DART_TYPE> + void RecSwapDelaunay( DART_TYPE& aDiagonal ); + + template <class TRAITS_TYPE, class DART_TYPE, class LIST_TYPE> + void SwapEdgesAwayFromInteriorNode( DART_TYPE& aDart, LIST_TYPE& aSwappedEdges ); + + template <class TRAITS_TYPE, class DART_TYPE, class LIST_TYPE> + void SwapEdgesAwayFromBoundaryNode( DART_TYPE& aDart, LIST_TYPE& aSwappedEdges ); + + template <class TRAITS_TYPE, class DART_TYPE, class DART_LIST_TYPE> + void SwapEdgeInList( const typename DART_LIST_TYPE::iterator& aIt, DART_LIST_TYPE& aElist ); + + // Constrained Triangulation + // ------------------------- + template <class TRAITS_TYPE, class DART_TYPE> + static DART_TYPE InsertConstraint( DART_TYPE& aDStart, DART_TYPE& aDEnd, bool aOptimizeDelaunay ); + +private: + hed::TRIANGULATION& m_triangulation; + + template <class TRAITS_TYPE, class FORWARD_ITERATOR, class DART_TYPE> + void insertNodes( FORWARD_ITERATOR aFirst, FORWARD_ITERATOR aLast, DART_TYPE& aDart ); + + template <class TOPOLOGY_ELEMENT_TYPE, class DART_TYPE> + static bool isMemberOfFace( const TOPOLOGY_ELEMENT_TYPE& aTopologyElement, const DART_TYPE& aDart ); + + template <class TRAITS_TYPE, class NODE_TYPE, class DART_TYPE> + static bool locateFaceWithNode( const NODE_TYPE& aNode, DART_TYPE& aDartIter ); + + template <class DART_TYPE> + static void getAdjacentTriangles( const DART_TYPE& aDart, DART_TYPE& aT1, DART_TYPE& aT2, + DART_TYPE& aT3 ); + + template <class DART_TYPE> + static void getNeighborNodes( const DART_TYPE& aDart, std::list<DART_TYPE>& aNodeList, + bool& aBoundary ); + + template <class TRAITS_TYPE, class DART_TYPE> + static bool degenerateTriangle( const DART_TYPE& aDart ); +}; + +#endif // DOXYGEN_SHOULD_SKIP_THIS + + + /** @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, + * 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 aDart + * An arbitrary CCW dart in the tringulation.\n + * Output: A CCW dart incident to the new node. + * + * \param aPoint + * 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 "TRAITS_TYPE::splitTriangle" (DART_TYPE&, const POINT_TYPE&) + * + * \using + * - locateTriangle + * - RecSwapDelaunay + * + * \note + * - For efficiency reasons \e dart should be close to the insertion \e point. + * + * \see + * removeRectangularBoundary + */ +template <class TRAITS_TYPE, class DART_TYPE, class POINT_TYPE> +bool TRIANGULATION_HELPER::InsertNode( DART_TYPE& aDart, POINT_TYPE& aPoint ) +{ + bool found = LocateTriangle<TRAITS_TYPE>( aPoint, aDart ); + + 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 + m_triangulation.splitTriangle( aDart, aPoint ); + + DART_TYPE d1 = aDart; + d1.Alpha2().Alpha1().Alpha2().Alpha0().Alpha1(); + + DART_TYPE d2 = aDart; + d2.Alpha1().Alpha0().Alpha1(); + + // Preserve a dart as output incident to the node and CCW + DART_TYPE d3 = aDart; + d3.Alpha2(); + aDart = d3; // and see below + //DART_TYPE dsav = d3; + d3.Alpha0().Alpha1(); + + //if (!TRAITS_TYPE::fixedEdge(d1) && !IsBoundaryEdge(d1)) { + if( !IsBoundaryEdge( d1 ) ) + { + d1.Alpha2(); + RecSwapDelaunay<TRAITS_TYPE>( d1 ); + } + + //if (!TRAITS_TYPE::fixedEdge(d2) && !IsBoundaryEdge(d2)) { + if( !IsBoundaryEdge( d2 ) ) + { + d2.Alpha2(); + RecSwapDelaunay<TRAITS_TYPE>( d2 ); + } + + // Preserve the incoming dart as output incident to the node and CCW + //d = dsav.Alpha2(); + aDart.Alpha2(); + //if (!TRAITS_TYPE::fixedEdge(d3) && !IsBoundaryEdge(d3)) { + if( !IsBoundaryEdge( d3 ) ) + { + d3.Alpha2(); + RecSwapDelaunay<TRAITS_TYPE>( d3 ); + } + + return true; +} + +//------------------------------------------------------------------------------------------------ +// Private/Hidden function (might change later) +template <class TRAITS_TYPE, class FORWARD_ITERATOR, class DART_TYPE> +void TRIANGULATION_HELPER::insertNodes( FORWARD_ITERATOR aFirst, FORWARD_ITERATOR aLast, + DART_TYPE& aDart ) +{ + + // Assumes that the dereferenced point objects are pointers. + // References to the point objects are then passed to TTL. + + FORWARD_ITERATOR it; + for( it = aFirst; it != aLast; ++it ) + { + InsertNode<TRAITS_TYPE>( aDart, **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 + * - RemoveBoundaryNode + * + * \note + * - This function requires that the boundary of the m_triangulation is + * a rectangle with four nodes (one in each corner). + */ +template <class TRAITS_TYPE, class DART_TYPE> +void TRIANGULATION_HELPER::RemoveRectangularBoundary( DART_TYPE& aDart ) +{ + DART_TYPE d_next = aDart; + DART_TYPE d_iter; + + for( int i = 0; i < 4; i++ ) + { + d_iter = d_next; + d_next.Alpha0(); + PositionAtNextBoundaryEdge( d_next ); + RemoveBoundaryNode<TRAITS_TYPE>( d_iter ); + } + + aDart = d_next; // Return a dart at the new boundary +} + +/** Removes the node associated with \e dart and + * updates the triangulation to be Delaunay. + * + * \using + * - RemoveBoundaryNode if \e dart represents a node at the boundary + * - 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 TRAITS_TYPE, class DART_TYPE> +void TRIANGULATION_HELPER::RemoveNode( DART_TYPE& aDart ) +{ + + if( isBoundaryNode( aDart ) ) + RemoveBoundaryNode<TRAITS_TYPE>( aDart ); + else + RemoveInteriorNode<TRAITS_TYPE>( aDart ); +} + +/** Removes the boundary node associated with \e dart and + * updates the triangulation to be Delaunay. + * + * \using + * - SwapEdgesAwayFromBoundaryNode + * - OptimizeDelaunay + * + * \require + * - \ref hed::TTLtraits::removeBoundaryTriangle "TRAITS_TYPE::removeBoundaryTriangle" (Dart&) + */ +template <class TRAITS_TYPE, class DART_TYPE> +void TRIANGULATION_HELPER::RemoveBoundaryNode( DART_TYPE& aDart ) +{ + + // ... 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 m_triangulation. + + // Position at boundary edge and CCW + if( !IsBoundaryEdge( aDart ) ) + { + aDart.Alpha1(); // ensures that next function delivers back a CCW dart (if the given dart is CCW) + PositionAtNextBoundaryEdge( aDart ); + } + + std::list<DART_TYPE> swapped_edges; + SwapEdgesAwayFromBoundaryNode<TRAITS_TYPE>( aDart, swapped_edges ); + + // Remove boundary triangles and remove the new boundary from the list + // of swapped edges, see below. + DART_TYPE d_iter = aDart; + DART_TYPE dnext = aDart; + bool bend = false; + while( bend == false ) + { + dnext.Alpha1().Alpha2(); + if( IsBoundaryEdge( dnext ) ) + bend = true; // Stop when boundary + + // Generic: Also remove the new boundary from the list of swapped edges + DART_TYPE 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 std::list<DART_TYPE>::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 + m_triangulation.removeBoundaryTriangle( d_iter ); + d_iter = dnext; + } + + // Optimize Delaunay + typedef std::list<DART_TYPE> DART_LIST_TYPE; + OptimizeDelaunay<TRAITS_TYPE, DART_TYPE, DART_LIST_TYPE>( swapped_edges ); +} + + +/** Removes the interior node associated with \e dart and + * updates the triangulation to be Delaunay. + * + * \using + * - SwapEdgesAwayFromInteriorNode + * - OptimizeDelaunay + * + * \require + * - \ref hed::TTLtraits::reverse_splitTriangle "TRAITS_TYPE::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 TRAITS_TYPE, class DART_TYPE> +void TRIANGULATION_HELPER::RemoveInteriorNode( DART_TYPE& aDart ) +{ + // ... 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 m_triangulation + // (No dart is delivered as output) + + // Assumes dart is counterclockwise + + std::list<DART_TYPE> swapped_edges; + SwapEdgesAwayFromInteriorNode<TRAITS_TYPE>( aDart, swapped_edges ); + + // The reverse operation of split triangle: + // Make one triangle of the three triangles at the node associated with dart + // TRAITS_TYPE:: + m_triangulation.reverseSplitTriangle( aDart ); + + // ???? 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 + + OptimizeDelaunay<TRAITS_TYPE, DART_TYPE>( swapped_edges ); +} + +//@} // End of Delaunay Triangulation Group + +/** @name Topological and Geometric Queries */ +//@{ +//------------------------------------------------------------------------------------------------ +// Private/Hidden function (might change later) +template <class TOPOLOGY_ELEMENT_TYPE, class DART_TYPE> +bool TRIANGULATION_HELPER::isMemberOfFace( const TOPOLOGY_ELEMENT_TYPE& aTopologyElement, + const DART_TYPE& aDart ) +{ + // Check if the given topology element (node, edge or face) is a member of the face + // Assumes: + // - DART_TYPE::isMember(TOPOLOGY_ELEMENT_TYPE) + + DART_TYPE dart_iter = aDart; + + do + { + if( dart_iter.isMember( aTopologyElement ) ) + return true; + dart_iter.Alpha0().Alpha1(); + } + while( dart_iter != aDart ); + + return false; +} + +//------------------------------------------------------------------------------------------------ +// Private/Hidden function (might change later) +template <class TRAITS_TYPE, class NODE_TYPE, class DART_TYPE> +bool TRIANGULATION_HELPER::locateFaceWithNode( const NODE_TYPE& aNode, DART_TYPE& aDartIter ) +{ + // Locate a face in the topology structure with the given node as a member + // Assumes: + // - TRAITS_TYPE::Orient2D(DART_TYPE, DART_TYPE, NODE_TYPE) + // - DART_TYPE::isMember(NODE_TYPE) + // - 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<TRAITS_TYPE>( aNode, aDartIter ); + + 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 POINT_TYPE& point, DART_TYPE& dart_iter) + + return isMemberOfFace( aNode, aDartIter ); +} + +/** 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 aPoint + * A point to be located + * + * \param aDart + * 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 "TRAITS_TYPE::Orient2D" (DART_TYPE&, DART_TYPE&, POINT_TYPE&) + * + * \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 + * LocateTriangle + */ +template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE> +bool TRIANGULATION_HELPER::LocateFaceSimplest( const POINT_TYPE& aPoint, DART_TYPE& aDart ) +{ + // 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? + // + // TRAITS_TYPE::orint2d(POINT_TYPE) 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())) + + DART_TYPE dart_start; + dart_start = aDart; + DART_TYPE dart_prev; + + DART_TYPE d0; + for( ;; ) + { + d0 = aDart; + d0.Alpha0(); + + if( TRAITS_TYPE::Orient2D( aDart, d0, aPoint ) >= 0 ) + { + aDart.Alpha0().Alpha1(); + if( aDart == dart_start ) + return true; // left to all edges in face + } + else + { + dart_prev = aDart; + aDart.Alpha2(); + + if( aDart == dart_prev ) + return false; // iteration to outside boundary + + dart_start = aDart; + dart_start.Alpha0(); + + aDart.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 m_triangulation, in which case \c false is returned, + * then the edge associated with \e dart will be at the boundary of the m_triangulation. + * + * \using + * - LocateFaceSimplest + * - InTriangle + */ +template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE> +bool TRIANGULATION_HELPER::LocateTriangle( const POINT_TYPE& aPoint, DART_TYPE& aDart ) +{ + // 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 m_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 POINT_TYPE& point, const DART_TYPE& dart). + + // Assumes: + // - CrossProduct2D, ScalarProduct2D etc., see functions called + + bool status = LocateFaceSimplest<TRAITS_TYPE>( aPoint, aDart ); + + 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<TRAITS_TYPE>( aPoint, aDart ); +} + +/** 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 aDart + * A CCW dart in the triangle + * + * \require + * - \ref hed::TTLtraits::CrossProduct2D "TRAITS_TYPE::CrossProduct2D" (DART_TYPE&, POINT_TYPE&) + * - \ref hed::TTLtraits::ScalarProduct2D "TRAITS_TYPE::ScalarProduct2D" (DART_TYPE&, POINT_TYPE&) + * + * \see + * InTriangleSimplest + */ +template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE> +bool TRIANGULATION_HELPER::InTriangle( const POINT_TYPE& aPoint, const DART_TYPE& aDart ) +{ + + // 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 aPoint coincides with all nodes; + // then false is always returned. + + typedef typename TRAITS_TYPE::REAL_TYPE REAL_TYPE; + + DART_TYPE dart_iter = aDart; + + REAL_TYPE cr1 = TRAITS_TYPE::CrossProduct2D( dart_iter, aPoint ); + if( cr1 < 0 ) + return false; + + dart_iter.Alpha0().Alpha1(); + REAL_TYPE cr2 = TRAITS_TYPE::CrossProduct2D( dart_iter, aPoint ); + + if( cr2 < 0 ) + return false; + + dart_iter.Alpha0().Alpha1(); + REAL_TYPE cr3 = TRAITS_TYPE::CrossProduct2D( dart_iter, aPoint ); + 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 aPoint is on (all) the nodes, + // then "false is returned". + + DART_TYPE dart_tmp = dart_iter; + REAL_TYPE sc1 = TRAITS_TYPE::ScalarProduct2D( dart_tmp, aPoint ); + REAL_TYPE sc2 = TRAITS_TYPE::ScalarProduct2D( dart_tmp.Alpha0(), aPoint ); + + 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 = TRAITS_TYPE::ScalarProduct2D( dart_tmp, aPoint ); + sc2 = TRAITS_TYPE::ScalarProduct2D( dart_tmp.Alpha0(), aPoint ); + + 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 = TRAITS_TYPE::ScalarProduct2D( dart_tmp, aPoint ); + sc2 = TRAITS_TYPE::ScalarProduct2D( dart_tmp.Alpha0(), aPoint ); + + 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 aPoint 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 DART_TYPE> +void TRIANGULATION_HELPER::getAdjacentTriangles( const DART_TYPE& aDart, DART_TYPE& aT1, + DART_TYPE& aT2, DART_TYPE& aT3 ) +{ + + DART_TYPE dart_iter = aDart; + + // add first + if( dart_iter.Alpha2() != aDart ) + { + aT1 = dart_iter; + dart_iter = aDart; + } + + // add second + dart_iter.Alpha0(); + dart_iter.Alpha1(); + DART_TYPE dart_prev = dart_iter; + + if( ( dart_iter.Alpha2() ) != dart_prev ) + { + aT2 = dart_iter; + dart_iter = dart_prev; + } + + // add third + dart_iter.Alpha0(); + dart_iter.Alpha1(); + dart_prev = dart_iter; + + if( ( dart_iter.Alpha2() ) != dart_prev ) + aT3 = 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 aDart + * A dart at the boundary (CCW or CW) + * + * \param aBoundary + * A sequence of darts, where the associated edges are the boundary edges + * + * \require + * - DART_LIST_TYPE::push_back (DART_TYPE&) + */ +template <class DART_TYPE, class DART_LIST_TYPE> +void TRIANGULATION_HELPER::GetBoundary( const DART_TYPE& aDart, DART_LIST_TYPE& aBoundary ) +{ + // assumes the given dart is at the boundary (by edge) + + DART_TYPE dart_iter( aDart ); + aBoundary.push_back( dart_iter ); // Given dart as first element + dart_iter.Alpha0(); + PositionAtNextBoundaryEdge( dart_iter ); + + while( dart_iter != aDart ) + { + aBoundary.push_back( dart_iter ); + dart_iter.Alpha0(); + PositionAtNextBoundaryEdge( dart_iter ); + } +} + +/** Checks if the edge associated with \e dart is at + * the boundary of the m_triangulation. + * + * \par Implements: + * \code + * DART_TYPE dart_iter = dart; + * if (dart_iter.Alpha2() == dart) + * return true; + * else + * return false; + * \endcode + */ +template <class DART_TYPE> +bool TRIANGULATION_HELPER::IsBoundaryEdge( const DART_TYPE& aDart ) +{ + DART_TYPE dart_iter = aDart; + + if( dart_iter.Alpha2() == aDart ) + return true; + else + return false; +} + +/** Checks if the face associated with \e dart is at + * the boundary of the m_triangulation. + */ +template <class DART_TYPE> +bool TRIANGULATION_HELPER::IsBoundaryFace( const DART_TYPE& aDart ) +{ + // Strategy: boundary if alpha2(d)=d + + DART_TYPE dart_iter( aDart ); + DART_TYPE 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 != aDart ); + + return false; +} + +/** Checks if the node associated with \e dart is at + * the boundary of the m_triangulation. + */ +template <class DART_TYPE> +bool TRIANGULATION_HELPER::IsBoundaryNode( const DART_TYPE& aDart ) +{ + // Strategy: boundary if alpha2(d)=d + + DART_TYPE dart_iter( aDart ); + DART_TYPE 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 != aDart ); + + 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 m_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 DART_TYPE> +int TRIANGULATION_HELPER::GetDegreeOfNode( const DART_TYPE& aDart ) +{ + DART_TYPE dart_iter( aDart ); + DART_TYPE 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 = aDart; + dart_iter.Alpha2(); + } else + return degree; + } + + } while( dart_iter != aDart ); + + 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 DART_TYPE> +void TRIANGULATION_HELPER::getNeighborNodes( const DART_TYPE& aDart, + std::list<DART_TYPE>& aNodeList, bool& aBoundary ) +{ + DART_TYPE dart_iter( aDart ); + dart_iter.Alpha0(); // position the dart at an opposite node + + DART_TYPE 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 + + DART_TYPE dart_start = dart_iter; + + do + { + aNodeList.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 + aBoundary = true; + + if( start_at_boundary == true ) + { + // add the dart which now is positioned at the opposite boundary + aNodeList.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(); + aNodeList.clear(); + getNeighborNodes( dart_iter, aNodeList, aBoundary ); + + return; // after one recursive step + } + } + } + while( dart_iter != dart_start ); + + aBoundary = false; +} + +/** Gets the 0-orbit around an interior node. + * + * \param aDart + * A dart (CCW or CW) positioned at an \e interior node. + * + * \retval aOrbit + * 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 + * - DART_LIST_TYPE::push_back (DART_TYPE&) + * + * \see + * Get0OrbitBoundary + */ +template <class DART_TYPE, class DART_LIST_TYPE> +void TRIANGULATION_HELPER::Get0OrbitInterior( const DART_TYPE& aDart, DART_LIST_TYPE& aOrbit ) +{ + DART_TYPE d_iter = aDart; + aOrbit.push_back( d_iter ); + d_iter.Alpha1().Alpha2(); + + while( d_iter != aDart ) + { + aOrbit.push_back( d_iter ); + d_iter.Alpha1().Alpha2(); + } +} + +/** Gets the 0-orbit around a node at the boundary + * + * \param aDart + * 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 + * - DART_LIST_TYPE::push_back (DART_TYPE&) + * + * \note + * - The last dart in the sequence have opposite orientation compared to the others! + * + * \see + * Get0OrbitInterior + */ +template <class DART_TYPE, class DART_LIST_TYPE> +void TRIANGULATION_HELPER::Get0OrbitBoundary( const DART_TYPE& aDart, DART_LIST_TYPE& aOrbit ) +{ + DART_TYPE dart_prev; + DART_TYPE d_iter = aDart; + + do + { + aOrbit.push_back( d_iter ); + d_iter.Alpha1(); + dart_prev = d_iter; + d_iter.Alpha2(); + } + while( d_iter != dart_prev ); + + aOrbit.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 DART_TYPE> +bool TRIANGULATION_HELPER::Same0Orbit( const DART_TYPE& aD1, const DART_TYPE& aD2 ) +{ + // Two copies of the same dart + DART_TYPE d_iter = aD2; + DART_TYPE d_end = aD2; + + if( isBoundaryNode( d_iter ) ) + { + // position at both boundary edges + PositionAtNextBoundaryEdge( d_iter ); + d_end.Alpha1(); + PositionAtNextBoundaryEdge( d_end ); + } + + for( ;; ) + { + if( d_iter == aD1 ) + return true; + + d_iter.Alpha1(); + + if( d_iter == aD1 ) + 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 DART_TYPE> +bool TRIANGULATION_HELPER::Same1Orbit( const DART_TYPE& aD1, const DART_TYPE& aD2 ) +{ + DART_TYPE d_iter = aD2; + + // (Also works at the boundary) + return ( d_iter == aD1 || d_iter.Alpha0() == aD1 || + d_iter.Alpha2() == aD1 || d_iter.Alpha0() == aD1 ); +} + +//------------------------------------------------------------------------------------------------ +/** 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 DART_TYPE> +bool TRIANGULATION_HELPER::Same2Orbit( const DART_TYPE& aD1, const DART_TYPE& aD2 ) +{ + DART_TYPE d_iter = aD2; + + return ( d_iter == aD1 || d_iter.Alpha0() == aD1 || d_iter.Alpha1() == aD1 || + d_iter.Alpha0() == aD1 || d_iter.Alpha1() == aD1 || d_iter.Alpha0() == aD1 ); +} + +// Private/Hidden function +template <class TRAITS_TYPE, class DART_TYPE> +bool TRIANGULATION_HELPER::degenerateTriangle( const DART_TYPE& aDart ) +{ + // Check if triangle is degenerate + // Assumes CCW dart + + DART_TYPE d1 = aDart; + DART_TYPE d2 = d1; + d2.Alpha1(); + + return ( TRAITS_TYPE::CrossProduct2D( d1, d2 ) == 0 ); +} + +/** 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 aAllowDegeneracy + * 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 "TRAITS_TYPE::CrossProduct2D" (Dart&, Dart&) + */ +template <class TRAITS_TYPE, class DART_TYPE> +bool TRIANGULATION_HELPER::SwappableEdge( const DART_TYPE& aDart, bool aAllowDegeneracy ) +{ + // How "safe" is it? + + if( IsBoundaryEdge( aDart ) ) + return false; + + // "angles" are at the diagonal + DART_TYPE d1 = aDart; + d1.Alpha2().Alpha1(); + DART_TYPE d2 = aDart; + d2.Alpha1(); + + if( aAllowDegeneracy ) + { + if( TRAITS_TYPE::CrossProduct2D( d1, d2 ) < 0.0 ) + return false; + } + else + { + if( TRAITS_TYPE::CrossProduct2D( d1, d2 ) <= 0.0 ) + return false; + } + + // Opposite side (still angle at the diagonal) + d1 = aDart; + d1.Alpha0(); + d2 = d1; + d1.Alpha1(); + d2.Alpha2().Alpha1(); + + if( aAllowDegeneracy ) + { + if( TRAITS_TYPE::CrossProduct2D( d1, d2 ) < 0.0 ) + return false; + } + else + { + if( TRAITS_TYPE::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 DART_TYPE> +void TRIANGULATION_HELPER::PositionAtNextBoundaryEdge( DART_TYPE& aDart ) +{ + DART_TYPE dart_prev; + + // If alpha2(d)=d, then boundary + + //old convention: dart.Alpha0(); + do + { + aDart.Alpha1(); + dart_prev = aDart; + aDart.Alpha2(); + } + while( aDart != dart_prev ); +} + +/** Checks if the boundary of a triangulation is convex. + * + * \param dart + * A CCW dart at the boundary of the m_triangulation + * + * \require + * - \ref hed::TTLtraits::CrossProduct2D "TRAITS_TYPE::CrossProduct2D" (const Dart&, const Dart&) + */ +template <class TRAITS_TYPE, class DART_TYPE> +bool TRIANGULATION_HELPER::ConvexBoundary( const DART_TYPE& aDart ) +{ + std::list<DART_TYPE> blist; + getBoundary( aDart, blist ); + + int no; + no = (int) blist.size(); + typename std::list<DART_TYPE>::const_iterator bit = blist.begin(); + DART_TYPE d1 = *bit; + ++bit; + DART_TYPE d2; + bool convex = true; + + for( ; bit != blist.end(); ++bit ) + { + d2 = *bit; + double crossProd = TRAITS_TYPE::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 = TRAITS_TYPE::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 + +/** @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 aElist + * The sequence of edges + * + * \require + * - \ref hed::TTLtraits::swapEdge "TRAITS_TYPE::swapEdge" (DART_TYPE& \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 + * - swapTestDelaunay + */ +template <class TRAITS_TYPE, class DART_TYPE, class DART_LIST_TYPE> +void TRIANGULATION_HELPER::OptimizeDelaunay( DART_LIST_TYPE& aElist ) +{ + OptimizeDelaunay<TRAITS_TYPE, DART_TYPE, DART_LIST_TYPE>( aElist, aElist.end() ); +} + +//------------------------------------------------------------------------------------------------ +template <class TRAITS_TYPE, class DART_TYPE, class DART_LIST_TYPE> +void TRIANGULATION_HELPER::OptimizeDelaunay( DART_LIST_TYPE& aElist, + const typename DART_LIST_TYPE::iterator aEnd ) +{ + // 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( aElist.empty() ) + return; + + // Avoid cycling by more extensive circumcircle test + bool cycling_check = true; + bool optimal = false; + typename DART_LIST_TYPE::iterator it; + + typename DART_LIST_TYPE::iterator end_opt = aEnd; + + // 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 DART_LIST_TYPE::iterator end_opt; + //if (*end != NULL) + // end_opt = end; + //else + // end_opt = elist.end(); + + while( !optimal ) + { + optimal = true; + for( it = aElist.begin(); it != end_opt; ++it ) + { + if( SwapTestDelaunay<TRAITS_TYPE>( *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) + + SwapEdgeInList<TRAITS_TYPE, DART_TYPE>( it, aElist ); + + 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 aCyclingCheck + * 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 "TRAITS_TYPE::ScalarProduct2D" (DART_TYPE&, DART_TYPE&) + * - \ref hed::TTLtraits::CrossProduct2D "TRAITS_TYPE::CrossProduct2D" (DART_TYPE&, DART_TYPE&) + */ +template <class TRAITS_TYPE, class DART_TYPE> +#if ((_MSC_VER > 0) && (_MSC_VER < 1300))//#ifdef _MSC_VER +bool TRIANGULATION_HELPER::SwapTestDelaunay(const DART_TYPE& aDart, bool aCyclingCheck = false) const +{ +#else +bool TRIANGULATION_HELPER::SwapTestDelaunay( const DART_TYPE& aDart, bool aCyclingCheck ) const +{ +#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 TRAITS_TYPE::REAL_TYPE REAL_TYPE; + + if( IsBoundaryEdge( aDart ) ) + return false; + + DART_TYPE v11 = aDart; + v11.Alpha1().Alpha0(); + DART_TYPE v12 = v11; + v12.Alpha1(); + + DART_TYPE v22 = aDart; + v22.Alpha2().Alpha1().Alpha0(); + DART_TYPE v21 = v22; + v21.Alpha1(); + + REAL_TYPE cos1 = TRAITS_TYPE::ScalarProduct2D( v11, v12 ); + REAL_TYPE cos2 = TRAITS_TYPE::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 = TRAITS_TYPE::CrossProduct2D( v11, v12 ); + REAL_TYPE sin2 = TRAITS_TYPE::CrossProduct2D( v21, v22 ); + REAL_TYPE sin12 = sin1 * cos2 + cos1 * sin2; + + if( sin12 >= 0 ) // equality represents a neutral case + return false; + + if( aCyclingCheck ) + { + // 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 TRAITS_TYPE::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 = TRAITS_TYPE::ScalarProduct2D( v22, v11 ); + cos2 = TRAITS_TYPE::ScalarProduct2D( v12, v21 ); + sin1 = TRAITS_TYPE::CrossProduct2D( v22, v11 ); + sin2 = TRAITS_TYPE::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 aDiagonal + * A CCW dart representing the edge where the recursion starts from. + * + * \require + * - \ref hed::TTLtraits::swapEdge "TRAITS_TYPE::swapEdge" (DART_TYPE&)\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 TRAITS_TYPE, class DART_TYPE> +void TRIANGULATION_HELPER::RecSwapDelaunay( DART_TYPE& aDiagonal ) +{ + if( !SwapTestDelaunay<TRAITS_TYPE>( aDiagonal ) ) + // ??? swapTestDelaunay also checks if boundary, so this can be optimized + return; + + // Get the other "edges" of the current triangle; see illustration above. + DART_TYPE oppEdge1 = aDiagonal; + oppEdge1.Alpha1(); + bool b1; + + if( IsBoundaryEdge( oppEdge1 ) ) + b1 = true; + else + { + b1 = false; + oppEdge1.Alpha2(); + } + + DART_TYPE oppEdge2 = aDiagonal; + oppEdge2.Alpha0().Alpha1().Alpha0(); + bool b2; + + if( IsBoundaryEdge( oppEdge2 ) ) + b2 = true; + else + { + b2 = false; + oppEdge2.Alpha2(); + } + + // Swap the given diagonal + m_triangulation.swapEdge( aDiagonal ); + + if( !b1 ) + RecSwapDelaunay<TRAITS_TYPE>( oppEdge1 ); + + if( !b2 ) + RecSwapDelaunay<TRAITS_TYPE>( 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 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 "TRAITS_TYPE::swapEdge" (DART_TYPE& \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 + * SwapEdgesAwayFromBoundaryNode + */ +template <class TRAITS_TYPE, class DART_TYPE, class LIST_TYPE> +void TRIANGULATION_HELPER::SwapEdgesAwayFromInteriorNode( DART_TYPE& aDart, + LIST_TYPE& aSwappedEdges ) +{ + + // Same iteration as in fixEdgesAtCorner, but not boundary + DART_TYPE dnext = aDart; + + // 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 = getDegreeOfNode( aDart ); + DART_TYPE d_iter; + + while( degree > 3 ) + { + d_iter = dnext; + dnext.Alpha1().Alpha2(); + + if( SwappableEdge<TRAITS_TYPE>( d_iter, allowDegeneracy ) ) + { + m_triangulation.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 + DART_TYPE swapped_edge = d_iter; // it was delivered back + swapped_edge.Alpha2().Alpha0(); // CCW (if not at boundary) + aSwappedEdges.push_back( swapped_edge ); + + degree--; + } + } + + // Output, incident to the node + aDart = 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 RemoveBoundaryNode + * + * \retval dart + * A CCW dart incident with the node + * + * \require + * - \ref hed::TTLtraits::swapEdge "TRAITS_TYPE::swapEdge" (DART_TYPE& \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 m_triangulation. + * + * \see + * SwapEdgesAwayFromInteriorNode + */ +template <class TRAITS_TYPE, class DART_TYPE, class LIST_TYPE> +void TRIANGULATION_HELPER::SwapEdgesAwayFromBoundaryNode( DART_TYPE& aDart, + LIST_TYPE& aSwappedEdges ) +{ + // 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 = getDegreeOfNode(dart); + + passes: + // Swap swappable edges that radiate from the node away + DART_TYPE d_iter = aDart; // ???? can simply use dart + d_iter.Alpha1().Alpha2(); // first not at boundary + DART_TYPE d_next = d_iter; + bool bend = false; + bool swapped_next_to_boundary = false; + bool swapped_in_pass = false; + + bool allowDegeneracy; // = true; + DART_TYPE tmp1, tmp2; + + while( !bend ) + { + d_next.Alpha1().Alpha2(); + + if( 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( IsBoundaryEdge( tmp1 ) && IsBoundaryEdge( tmp2 ) ) + allowDegeneracy = true; + else + allowDegeneracy = false; + + if( SwappableEdge<TRAITS_TYPE>( d_iter, allowDegeneracy ) ) + { + m_triangulation.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 + DART_TYPE swapped_edge = d_iter; // it was delivered back + swapped_edge.Alpha2().Alpha0(); // CCW + aSwappedEdges.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 + } + PositionAtNextBoundaryEdge( d_iter ); // CCW + + aDart = 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 TRAITS_TYPE, class DART_TYPE, class DART_LIST_TYPE> +void TRIANGULATION_HELPER::SwapEdgeInList( const typename DART_LIST_TYPE::iterator& aIt, + DART_LIST_TYPE& aElist ) +{ + + typename DART_LIST_TYPE::iterator it1, it2, it3, it4; + DART_TYPE dart( *aIt ); + + //typename TRAITS_TYPE::DART_TYPE d1 = dart; d1.Alpha2().Alpha1(); + //typename TRAITS_TYPE::DART_TYPE d2 = d1; d2.Alpha0().Alpha1(); + //typename TRAITS_TYPE::DART_TYPE d3 = dart; d3.Alpha0().Alpha1(); + //typename TRAITS_TYPE::DART_TYPE d4 = d3; d4.Alpha0().Alpha1(); + DART_TYPE d1 = dart; + d1.Alpha2().Alpha1(); + DART_TYPE d2 = d1; + d2.Alpha0().Alpha1(); + DART_TYPE d3 = dart; + d3.Alpha0().Alpha1(); + DART_TYPE 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( aElist.begin(), aElist.end(), d1 ); + it2 = find( aElist.begin(), aElist.end(), d2 ); + it3 = find( aElist.begin(), aElist.end(), d3 ); + it4 = find( aElist.begin(), aElist.end(), d4 ); + + m_triangulation.swapEdge( dart ); + // Update the current dart which may have changed + *aIt = dart; + + // Update darts that may have changed again (if they were present) + // Note that dart is delivered back after swapping + if( it1 != aElist.end() ) + { + d1 = dart; + d1.Alpha1().Alpha0(); + *it1 = d1; + } + + if( it2 != aElist.end() ) + { + d2 = dart; + d2.Alpha2().Alpha1(); + *it2 = d2; + } + + if( it3 != aElist.end() ) + { + d3 = dart; + d3.Alpha2().Alpha1().Alpha0().Alpha1(); + *it3 = d3; + } + + if( it4 != aElist.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) + +#endif // _TTL_H_ |