ttl.h 63 KB
Newer Older
Maciej Suminski's avatar
Maciej Suminski committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
/*
 * 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
48 49 50
static void errorAndExit( char* aMessage )
{
    cout << "\n!!! ERROR: " << aMessage << " !!!\n" << endl;
Maciej Suminski's avatar
Maciej Suminski committed
51
    exit(-1);
52
}
Maciej Suminski's avatar
Maciej Suminski committed
53 54 55 56 57 58 59 60 61 62 63 64
#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.

65 66
/**
 * \brief Main interface to TTL
Maciej Suminski's avatar
Maciej Suminski committed
67
*
68 69
* This namespace contains the basic generic algorithms for the TTL,
* the Triangulation Template Library.\n
Maciej Suminski's avatar
Maciej Suminski committed
70
*
71 72 73 74 75 76
* 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.
Maciej Suminski's avatar
Maciej Suminski committed
77
*
78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
* \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.
Maciej Suminski's avatar
Maciej Suminski committed
95
*
96 97 98 99 100 101 102 103
* \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
Maciej Suminski's avatar
Maciej Suminski committed
104
*
105 106
* \see
* \ref ttl_util and \ref api
Maciej Suminski's avatar
Maciej Suminski committed
107
*
108 109
* \author
* �yvind Hjelle, oyvindhj@ifi.uio.no
Maciej Suminski's avatar
Maciej Suminski committed
110 111 112
*/


113 114 115
namespace ttl
{
class TRIANGULATION_HELPER
116
{
Maciej Suminski's avatar
Maciej Suminski committed
117 118
#ifndef DOXYGEN_SHOULD_SKIP_THIS

119
public:
120 121 122 123
    TRIANGULATION_HELPER( hed::TRIANGULATION& aTriang ) :
        m_triangulation( aTriang )
    {
    }
Maciej Suminski's avatar
Maciej Suminski committed
124

125 126 127
    // Delaunay Triangulation
    template <class TRAITS_TYPE, class DART_TYPE, class POINT_TYPE>
    bool InsertNode( DART_TYPE& aDart, POINT_TYPE& aPoint );
Maciej Suminski's avatar
Maciej Suminski committed
128

129 130
    template <class TRAITS_TYPE, class DART_TYPE>
    void RemoveRectangularBoundary( DART_TYPE& aDart );
Maciej Suminski's avatar
Maciej Suminski committed
131

132 133
    template <class TRAITS_TYPE, class DART_TYPE>
    void RemoveNode( DART_TYPE& aDart );
Maciej Suminski's avatar
Maciej Suminski committed
134

135 136
    template <class TRAITS_TYPE, class DART_TYPE>
    void RemoveBoundaryNode( DART_TYPE& aDart );
Maciej Suminski's avatar
Maciej Suminski committed
137

138 139
    template <class TRAITS_TYPE, class DART_TYPE>
    void RemoveInteriorNode( DART_TYPE& aDart );
Maciej Suminski's avatar
Maciej Suminski committed
140

141 142 143 144
    // Topological and Geometric Queries
    // ---------------------------------
    template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE>
    static bool LocateFaceSimplest( const POINT_TYPE& aPoint, DART_TYPE& aDart );
Maciej Suminski's avatar
Maciej Suminski committed
145

146 147
    template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE>
    static bool LocateTriangle( const POINT_TYPE& aPoint, DART_TYPE& aDart );
Maciej Suminski's avatar
Maciej Suminski committed
148

149 150
    template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE>
    static bool InTriangleSimplest( const POINT_TYPE& aPoint, const DART_TYPE& aDart );
Maciej Suminski's avatar
Maciej Suminski committed
151

152 153
    template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE>
    static bool InTriangle( const POINT_TYPE& aPoint, const DART_TYPE& aDart );
Maciej Suminski's avatar
Maciej Suminski committed
154

155 156
    template <class DART_TYPE, class DART_LIST_TYPE>
    static void GetBoundary( const DART_TYPE& aDart, DART_LIST_TYPE& aBoundary );
Maciej Suminski's avatar
Maciej Suminski committed
157

158 159
    template <class DART_TYPE>
    static bool IsBoundaryEdge( const DART_TYPE& aDart );
Maciej Suminski's avatar
Maciej Suminski committed
160

161 162
    template <class DART_TYPE>
    static bool IsBoundaryFace( const DART_TYPE& aDart );
Maciej Suminski's avatar
Maciej Suminski committed
163

164 165
    template <class DART_TYPE>
    static bool IsBoundaryNode( const DART_TYPE& aDart );
Maciej Suminski's avatar
Maciej Suminski committed
166

167 168
    template <class DART_TYPE>
    static int GetDegreeOfNode( const DART_TYPE& aDart );
Maciej Suminski's avatar
Maciej Suminski committed
169

170 171
    template <class DART_TYPE, class DART_LIST_TYPE>
    static void Get0OrbitInterior( const DART_TYPE& aDart, DART_LIST_TYPE& aOrbit );
Maciej Suminski's avatar
Maciej Suminski committed
172

173 174
    template <class DART_TYPE, class DART_LIST_TYPE>
    static void Get0OrbitBoundary( const DART_TYPE& aDart, DART_LIST_TYPE& aOrbit );
Maciej Suminski's avatar
Maciej Suminski committed
175

176 177
    template <class DART_TYPE>
    static bool Same0Orbit( const DART_TYPE& aD1, const DART_TYPE& aD2 );
Maciej Suminski's avatar
Maciej Suminski committed
178

179 180
    template <class DART_TYPE>
    static bool Same1Orbit( const DART_TYPE& aD1, const DART_TYPE& aD2 );
Maciej Suminski's avatar
Maciej Suminski committed
181

182 183
    template <class DART_TYPE>
    static bool Same2Orbit( const DART_TYPE& aD1, const DART_TYPE& aD2 );
Maciej Suminski's avatar
Maciej Suminski committed
184

185 186
    template <class TRAITS_TYPE, class DART_TYPE>
    static bool SwappableEdge( const DART_TYPE& aDart, bool aAllowDegeneracy = false );
Maciej Suminski's avatar
Maciej Suminski committed
187

188 189
    template <class DART_TYPE>
    static void PositionAtNextBoundaryEdge( DART_TYPE& aDart );
Maciej Suminski's avatar
Maciej Suminski committed
190

191 192
    template <class TRAITS_TYPE, class DART_TYPE>
    static bool ConvexBoundary( const DART_TYPE& aDart );
Maciej Suminski's avatar
Maciej Suminski committed
193

194 195 196 197
    // Utilities for Delaunay Triangulation
    // ------------------------------------
    template <class TRAITS_TYPE, class DART_TYPE, class DART_LIST_TYPE>
    void OptimizeDelaunay( DART_LIST_TYPE& aElist );
Maciej Suminski's avatar
Maciej Suminski committed
198

199 200
    template <class TRAITS_TYPE, class DART_TYPE, class DART_LIST_TYPE>
    void OptimizeDelaunay( DART_LIST_TYPE& aElist, const typename DART_LIST_TYPE::iterator aEnd );
Maciej Suminski's avatar
Maciej Suminski committed
201

202 203
    template <class TRAITS_TYPE, class DART_TYPE>
    bool SwapTestDelaunay( const DART_TYPE& aDart, bool aCyclingCheck = false ) const;
Maciej Suminski's avatar
Maciej Suminski committed
204

205 206
    template <class TRAITS_TYPE, class DART_TYPE>
    void RecSwapDelaunay( DART_TYPE& aDiagonal );
Maciej Suminski's avatar
Maciej Suminski committed
207

208 209
    template <class TRAITS_TYPE, class DART_TYPE, class LIST_TYPE>
    void SwapEdgesAwayFromInteriorNode( DART_TYPE& aDart, LIST_TYPE& aSwappedEdges );
Maciej Suminski's avatar
Maciej Suminski committed
210

211 212
    template <class TRAITS_TYPE, class DART_TYPE, class LIST_TYPE>
    void SwapEdgesAwayFromBoundaryNode( DART_TYPE& aDart, LIST_TYPE& aSwappedEdges );
Maciej Suminski's avatar
Maciej Suminski committed
213

214 215
    template <class TRAITS_TYPE, class DART_TYPE, class DART_LIST_TYPE>
    void SwapEdgeInList( const typename DART_LIST_TYPE::iterator& aIt, DART_LIST_TYPE& aElist );
Maciej Suminski's avatar
Maciej Suminski committed
216

217 218 219 220
    // Constrained Triangulation
    // -------------------------
    template <class TRAITS_TYPE, class DART_TYPE>
    static DART_TYPE InsertConstraint( DART_TYPE& aDStart, DART_TYPE& aDEnd, bool aOptimizeDelaunay );
221 222

private:
223
    hed::TRIANGULATION& m_triangulation;
224

225 226
    template <class TRAITS_TYPE, class FORWARD_ITERATOR, class DART_TYPE>
    void insertNodes( FORWARD_ITERATOR aFirst, FORWARD_ITERATOR aLast, DART_TYPE& aDart );
227

228 229
    template <class TOPOLOGY_ELEMENT_TYPE, class DART_TYPE>
    static bool isMemberOfFace( const TOPOLOGY_ELEMENT_TYPE& aTopologyElement, const DART_TYPE& aDart );
230

231 232
    template <class TRAITS_TYPE, class NODE_TYPE, class DART_TYPE>
    static bool locateFaceWithNode( const NODE_TYPE& aNode, DART_TYPE& aDartIter );
233

234 235 236
    template <class DART_TYPE>
    static void getAdjacentTriangles( const DART_TYPE& aDart, DART_TYPE& aT1, DART_TYPE& aT2,
                                      DART_TYPE& aT3 );
237

238 239 240
    template <class DART_TYPE>
    static void getNeighborNodes( const DART_TYPE& aDart, std::list<DART_TYPE>& aNodeList,
                                  bool& aBoundary );
241

242 243
    template <class TRAITS_TYPE, class DART_TYPE>
    static bool degenerateTriangle( const DART_TYPE& aDart );
244
};
Maciej Suminski's avatar
Maciej Suminski committed
245 246 247 248 249

#endif // DOXYGEN_SHOULD_SKIP_THIS


  /** @name Delaunay Triangulation */
250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294
//@{
/**
 * 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 );
Maciej Suminski's avatar
Maciej Suminski committed
295
    
296 297
    if( !found )
    {
Maciej Suminski's avatar
Maciej Suminski committed
298
#ifdef DEBUG_TTL
299
        cout << "ERROR: Triangulation::insertNode: NO triangle found. /n";
Maciej Suminski's avatar
Maciej Suminski committed
300
#endif
301
        return false;
Maciej Suminski's avatar
Maciej Suminski committed
302
    }
303

Maciej Suminski's avatar
Maciej Suminski committed
304
    // ??? can we hide the dart? this is not possible if one triangle only
305 306 307 308 309 310 311 312
    m_triangulation.splitTriangle( aDart, aPoint );

    DART_TYPE d1 = aDart;
    d1.Alpha2().Alpha1().Alpha2().Alpha0().Alpha1();

    DART_TYPE d2 = aDart;
    d2.Alpha1().Alpha0().Alpha1();

Maciej Suminski's avatar
Maciej Suminski committed
313
    // Preserve a dart as output incident to the node and CCW
314 315 316 317 318 319 320 321 322 323 324
    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 );
Maciej Suminski's avatar
Maciej Suminski committed
325
    }
326 327 328 329 330 331

    //if (!TRAITS_TYPE::fixedEdge(d2) && !IsBoundaryEdge(d2)) {
    if( !IsBoundaryEdge( d2 ) )
    {
        d2.Alpha2();
        RecSwapDelaunay<TRAITS_TYPE>( d2 );
Maciej Suminski's avatar
Maciej Suminski committed
332
    }
333

Maciej Suminski's avatar
Maciej Suminski committed
334
    // Preserve the incoming dart as output incident to the node and CCW
335 336 337 338 339 340 341
    //d = dsav.Alpha2();
    aDart.Alpha2();
    //if (!TRAITS_TYPE::fixedEdge(d3) && !IsBoundaryEdge(d3)) {
    if( !IsBoundaryEdge( d3 ) )
    {
        d3.Alpha2();
        RecSwapDelaunay<TRAITS_TYPE>( d3 );
Maciej Suminski's avatar
Maciej Suminski committed
342
    }
343

Maciej Suminski's avatar
Maciej Suminski committed
344
    return true;
345
}
Maciej Suminski's avatar
Maciej Suminski committed
346

347 348 349 350 351 352
//------------------------------------------------------------------------------------------------
// 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 )
{
Maciej Suminski's avatar
Maciej Suminski committed
353 354 355

    // Assumes that the dereferenced point objects are pointers.
    // References to the point objects are then passed to TTL.
356 357 358 359 360

    FORWARD_ITERATOR it;
    for( it = aFirst; it != aLast; ++it )
    {
        InsertNode<TRAITS_TYPE>( aDart, **it );
Maciej Suminski's avatar
Maciej Suminski committed
361
    }
362
}
Maciej Suminski's avatar
Maciej Suminski committed
363 364


365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392
/** 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 );
Maciej Suminski's avatar
Maciej Suminski committed
393 394
    }

395 396
    aDart = d_next; // Return a dart at the new boundary
}
Maciej Suminski's avatar
Maciej Suminski committed
397

398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414
/** 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 );
Maciej Suminski's avatar
Maciej Suminski committed
415
    else
416 417
        RemoveInteriorNode<TRAITS_TYPE>( aDart );
}
Maciej Suminski's avatar
Maciej Suminski committed
418

419 420 421 422 423 424 425 426 427 428 429 430 431
/** 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 )
{
Maciej Suminski's avatar
Maciej Suminski committed
432 433 434 435

    // ... and update Delaunay
    // - CCW dart must be given (for remove)
    // - No dart is delivered back now (but this is possible if
436 437
    //   we assume that there is not only one triangle left in the m_triangulation.

Maciej Suminski's avatar
Maciej Suminski committed
438
    // Position at boundary edge and CCW
439 440 441 442
    if( !IsBoundaryEdge( aDart ) )
    {
        aDart.Alpha1(); // ensures that next function delivers back a CCW dart (if the given dart is CCW)
        PositionAtNextBoundaryEdge( aDart );
Maciej Suminski's avatar
Maciej Suminski committed
443 444
    }

445 446 447
    std::list<DART_TYPE> swapped_edges;
    SwapEdgesAwayFromBoundaryNode<TRAITS_TYPE>( aDart, swapped_edges );

Maciej Suminski's avatar
Maciej Suminski committed
448 449
    // Remove boundary triangles and remove the new boundary from the list
    // of swapped edges, see below.
450 451
    DART_TYPE d_iter = aDart;
    DART_TYPE dnext = aDart;
Maciej Suminski's avatar
Maciej Suminski committed
452
    bool bend = false;
453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472
    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;
Maciej Suminski's avatar
Maciej Suminski committed
473
    }
474

Maciej Suminski's avatar
Maciej Suminski committed
475
    // Optimize Delaunay
476 477 478
    typedef std::list<DART_TYPE> DART_LIST_TYPE;
    OptimizeDelaunay<TRAITS_TYPE, DART_TYPE, DART_LIST_TYPE>( swapped_edges );
}
Maciej Suminski's avatar
Maciej Suminski committed
479 480


481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497
/** 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 )
{
Maciej Suminski's avatar
Maciej Suminski committed
498 499 500 501 502
    // ... 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.
503

Maciej Suminski's avatar
Maciej Suminski committed
504 505 506
    // 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
507
    // 3) Runs LOP on the platelet to obtain a Delaunay m_triangulation
Maciej Suminski's avatar
Maciej Suminski committed
508
    // (No dart is delivered as output)
509

Maciej Suminski's avatar
Maciej Suminski committed
510
    // Assumes dart is counterclockwise
511 512 513 514

    std::list<DART_TYPE> swapped_edges;
    SwapEdgesAwayFromInteriorNode<TRAITS_TYPE>( aDart, swapped_edges );

Maciej Suminski's avatar
Maciej Suminski committed
515 516
    // The reverse operation of split triangle:
    // Make one triangle of the three triangles at the node associated with dart
517 518 519
    // TRAITS_TYPE::
    m_triangulation.reverseSplitTriangle( aDart );

Maciej Suminski's avatar
Maciej Suminski committed
520 521 522 523 524 525
    // ???? 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???
526

Maciej Suminski's avatar
Maciej Suminski committed
527 528 529 530
    // Make the swapped-away edges Delaunay.
    // Note the theoretical result: if there are no edges in the list,
    // the triangulation is Delaunay already

531 532
    OptimizeDelaunay<TRAITS_TYPE, DART_TYPE>( swapped_edges );
}
Maciej Suminski's avatar
Maciej Suminski committed
533

534
//@} // End of Delaunay Triangulation Group
Maciej Suminski's avatar
Maciej Suminski committed
535

536 537 538 539 540 541 542 543
/** @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 )
{
Maciej Suminski's avatar
Maciej Suminski committed
544 545
    // Check if the given topology element (node, edge or face) is a member of the face
    // Assumes:
546 547 548
    // - DART_TYPE::isMember(TOPOLOGY_ELEMENT_TYPE)

    DART_TYPE dart_iter = aDart;
Maciej Suminski's avatar
Maciej Suminski committed
549
    
550 551 552 553 554 555 556
    do
    {
        if( dart_iter.isMember( aTopologyElement ) )
            return true;
        dart_iter.Alpha0().Alpha1();
    }
    while( dart_iter != aDart );
Maciej Suminski's avatar
Maciej Suminski committed
557

558 559
    return false;
}
Maciej Suminski's avatar
Maciej Suminski committed
560

561 562 563 564 565
//------------------------------------------------------------------------------------------------
// 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 )
{
Maciej Suminski's avatar
Maciej Suminski committed
566 567
    // Locate a face in the topology structure with the given node as a member
    // Assumes:
568 569
    // - TRAITS_TYPE::Orient2D(DART_TYPE, DART_TYPE, NODE_TYPE)
    // - DART_TYPE::isMember(NODE_TYPE)
Maciej Suminski's avatar
Maciej Suminski committed
570 571 572
    // - 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;
573 574 575 576 577 578 579 580
    //   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
Maciej Suminski's avatar
Maciej Suminski committed
581 582 583
    // 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
584
    // and in locateTriangle(const POINT_TYPE& point, DART_TYPE& dart_iter)
Maciej Suminski's avatar
Maciej Suminski committed
585

586 587
    return isMemberOfFace( aNode, aDartIter );
}
Maciej Suminski's avatar
Maciej Suminski committed
588

589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615
/** 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 )
{
Maciej Suminski's avatar
Maciej Suminski committed
616 617 618 619 620 621 622
    // 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?
    //
623
    // TRAITS_TYPE::orint2d(POINT_TYPE) is the half open half-plane defined
Maciej Suminski's avatar
Maciej Suminski committed
624 625
    // by the dart:
    // n1 = dart.node()
626
    // n2 = dart.Alpha0().node
Maciej Suminski's avatar
Maciej Suminski committed
627 628
    // Only the following gives true:
    // ((n2->x()-n1->x())*(point.y()-n1->y()) >= (point.x()-n1->x())*(n2->y()-n1->y()))
629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658

    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
        }
Maciej Suminski's avatar
Maciej Suminski committed
659
    }
660
}
Maciej Suminski's avatar
Maciej Suminski committed
661 662


663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687
/** 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 )
{
Maciej Suminski's avatar
Maciej Suminski committed
688 689 690
    // 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
691

Maciej Suminski's avatar
Maciej Suminski committed
692 693 694
    // 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.
695
    // ii) The boundary of the m_triangulation is not convex.
Maciej Suminski's avatar
Maciej Suminski committed
696 697 698
    // 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.
699

Maciej Suminski's avatar
Maciej Suminski committed
700 701 702 703 704
    // 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:
705 706
    // bool inTriangle(const POINT_TYPE& point, const DART_TYPE& dart).

Maciej Suminski's avatar
Maciej Suminski committed
707
    // Assumes:
708 709 710
    // - CrossProduct2D, ScalarProduct2D etc., see functions called

    bool status = LocateFaceSimplest<TRAITS_TYPE>( aPoint, aDart );
Maciej Suminski's avatar
Maciej Suminski committed
711
    
712 713 714
    if( status == false )
        return status;

Maciej Suminski's avatar
Maciej Suminski committed
715 716
    // There may be degeneracy, i.e., the point might be outside the triangle
    // on the extension of the edges of a degenerate triangle.
717

Maciej Suminski's avatar
Maciej Suminski committed
718 719 720
    // 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.
721 722
    return InTriangle<TRAITS_TYPE>( aPoint, aDart );
}
Maciej Suminski's avatar
Maciej Suminski committed
723

724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739
//------------------------------------------------------------------------------------------------
/** Checks if \e point is inside the triangle associated with \e dart.
 *   A fast and simple function that does not deal with degeneracy.
 *
 *   \param aDart
 *   A CCW dart in the triangle
 *
 *   \require
 *   - \ref hed::TTLtraits::Orient2D "TRAITS_TYPE::Orient2D" (DART_TYPE&, DART_TYPE&, POINT_TYPE&)
 *
 *   \see
 *   InTriangle for a more robust function
 */
template <class TRAITS_TYPE, class POINT_TYPE, class DART_TYPE>
bool TRIANGULATION_HELPER::InTriangleSimplest( const POINT_TYPE& aPoint, const DART_TYPE& aDart )
{
Maciej Suminski's avatar
Maciej Suminski committed
740 741 742
    // 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
743 744 745 746

    DART_TYPE d_iter = aDart;
    DART_TYPE d0 = d_iter;
    d0.Alpha0();
Maciej Suminski's avatar
Maciej Suminski committed
747
    
748 749 750 751
    if( !TRAITS_TYPE::Orient2D( d_iter, d0, aPoint ) >= 0 )
        return false;

    d_iter.Alpha0().Alpha1();
Maciej Suminski's avatar
Maciej Suminski committed
752
    d0 = d_iter;
753
    d0.Alpha0();
Maciej Suminski's avatar
Maciej Suminski committed
754
    
755 756 757 758
    if( !TRAITS_TYPE::Orient2D( d_iter, d0, aPoint ) >= 0 )
        return false;

    d_iter.Alpha0().Alpha1();
Maciej Suminski's avatar
Maciej Suminski committed
759
    d0 = d_iter;
760
    d0.Alpha0();
Maciej Suminski's avatar
Maciej Suminski committed
761
    
762 763 764
    if( !TRAITS_TYPE::Orient2D( d_iter, d0, aPoint ) >= 0 )
        return false;

Maciej Suminski's avatar
Maciej Suminski committed
765
    return true;
766
}
Maciej Suminski's avatar
Maciej Suminski committed
767

768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784
/** 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 )
{
Maciej Suminski's avatar
Maciej Suminski committed
785 786 787 788

    // 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 ?
789

Maciej Suminski's avatar
Maciej Suminski committed
790 791
    // Returns true if inside (but not necessarily strictly inside)
    // Works for degenerate triangles, but not when all edges are degenerate,
792
    // and the aPoint coincides with all nodes;
Maciej Suminski's avatar
Maciej Suminski committed
793
    // then false is always returned.
794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813

    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;

Maciej Suminski's avatar
Maciej Suminski committed
814 815
    // All cross products are >= 0
    // Check for degeneracy
816 817 818
    if( cr1 != 0 || cr2 != 0 || cr3 != 0 )
        return true; // inside non-degenerate face

Maciej Suminski's avatar
Maciej Suminski committed
819
    // All cross-products are zero, i.e. degenerate triangle, check if inside
820 821
    // 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,
Maciej Suminski's avatar
Maciej Suminski committed
822
    // then "false is returned".
823 824 825 826 827 828 829 830 831 832

    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)
Maciej Suminski's avatar
Maciej Suminski committed
833
    }
834 835 836 837

    dart_tmp = dart_iter.Alpha0().Alpha1();
    sc1 = TRAITS_TYPE::ScalarProduct2D( dart_tmp, aPoint );
    sc2 = TRAITS_TYPE::ScalarProduct2D( dart_tmp.Alpha0(), aPoint );
Maciej Suminski's avatar
Maciej Suminski committed
838
    
839 840 841 842 843
    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)
Maciej Suminski's avatar
Maciej Suminski committed
844
    }
845 846 847 848

    dart_tmp = dart_iter.Alpha1();
    sc1 = TRAITS_TYPE::ScalarProduct2D( dart_tmp, aPoint );
    sc2 = TRAITS_TYPE::ScalarProduct2D( dart_tmp.Alpha0(), aPoint );
Maciej Suminski's avatar
Maciej Suminski committed
849
    
850 851 852 853 854
    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)
Maciej Suminski's avatar
Maciej Suminski committed
855
    }
856

Maciej Suminski's avatar
Maciej Suminski committed
857
    // Not on any of the edges of the degenerate triangle.
858
    // The only possibility for the aPoint to be "inside" is that all edges are degenerate
Maciej Suminski's avatar
Maciej Suminski committed
859
    // and the point coincide with all nodes. So false is returned in this case.
860

Maciej Suminski's avatar
Maciej Suminski committed
861
    return false;
862
}
Maciej Suminski's avatar
Maciej Suminski committed
863 864 865


  //------------------------------------------------------------------------------------------------
866 867 868 869 870 871 872 873
// 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;

Maciej Suminski's avatar
Maciej Suminski committed
874
    // add first
875 876 877 878
    if( dart_iter.Alpha2() != aDart )
    {
        aT1 = dart_iter;
        dart_iter = aDart;
Maciej Suminski's avatar
Maciej Suminski committed
879
    }
880

Maciej Suminski's avatar
Maciej Suminski committed
881
    // add second
882 883 884 885 886 887 888 889
    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;
Maciej Suminski's avatar
Maciej Suminski committed
890
    }
891

Maciej Suminski's avatar
Maciej Suminski committed
892
    // add third
893 894
    dart_iter.Alpha0();
    dart_iter.Alpha1();
Maciej Suminski's avatar
Maciej Suminski committed
895 896
    dart_prev = dart_iter;

897 898 899
    if( ( dart_iter.Alpha2() ) != dart_prev )
        aT3 = dart_iter;
}
Maciej Suminski's avatar
Maciej Suminski committed
900

901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919
//------------------------------------------------------------------------------------------------
/** 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 )
{
Maciej Suminski's avatar
Maciej Suminski committed
920 921
    // assumes the given dart is at the boundary (by edge)

922 923 924 925
    DART_TYPE dart_iter( aDart );
    aBoundary.push_back( dart_iter ); // Given dart as first element
    dart_iter.Alpha0();
    PositionAtNextBoundaryEdge( dart_iter );
Maciej Suminski's avatar
Maciej Suminski committed
926

927 928 929 930 931 932 933
    while( dart_iter != aDart )
    {
        aBoundary.push_back( dart_iter );
        dart_iter.Alpha0();
        PositionAtNextBoundaryEdge( dart_iter );
    }
}
Maciej Suminski's avatar
Maciej Suminski committed
934

935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950
/** 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;
Maciej Suminski's avatar
Maciej Suminski committed
951
    
952 953
    if( dart_iter.Alpha2() == aDart )
        return true;
Maciej Suminski's avatar
Maciej Suminski committed
954
    else
955 956
        return false;
}
Maciej Suminski's avatar
Maciej Suminski committed
957

958 959 960 961 962 963
/** 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 )
{
Maciej Suminski's avatar
Maciej Suminski committed
964 965
    // Strategy: boundary if alpha2(d)=d

966 967
    DART_TYPE dart_iter( aDart );
    DART_TYPE dart_prev;
Maciej Suminski's avatar
Maciej Suminski committed
968

969 970 971
    do
    {
        dart_prev = dart_iter;
Maciej Suminski's avatar
Maciej Suminski committed
972

973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025
        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;
Maciej Suminski's avatar
Maciej Suminski committed
1026 1027 1028

    // If input dart is reached again, then interior node
    // If alpha2(d)=d, then boundary
1029

Maciej Suminski's avatar
Maciej Suminski committed
1030 1031
    int degree = 0;
    bool boundaryVisited = false;
1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050
    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;
Maciej Suminski's avatar
Maciej Suminski committed
1051
        }
1052 1053 1054

    } while( dart_iter != aDart );

Maciej Suminski's avatar
Maciej Suminski committed
1055
    return degree;
1056
}
Maciej Suminski's avatar
Maciej Suminski committed
1057

1058 1059 1060 1061
// 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.
Maciej Suminski's avatar
Maciej Suminski committed
1062

1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087
// 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;
Maciej Suminski's avatar
Maciej Suminski committed
1088
    bool start_at_boundary = false;
1089 1090 1091 1092
    dart_iter.Alpha2();

    if( dart_iter == dart_prev )
        start_at_boundary = true;
Maciej Suminski's avatar
Maciej Suminski committed
1093
    else
1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127
        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
            }
Maciej Suminski's avatar
Maciej Suminski committed
1128
        }
1129 1130
    }
    while( dart_iter != dart_start );
Maciej Suminski's avatar
Maciej Suminski committed
1131

1132 1133
    aBoundary = false;
}
Maciej Suminski's avatar
Maciej Suminski committed
1134

1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161
/** 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();
Maciej Suminski's avatar
Maciej Suminski committed
1162
    }
1163
}
Maciej Suminski's avatar
Maciej Suminski committed
1164

1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197
/** 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 );
Maciej Suminski's avatar
Maciej Suminski committed
1198

1199 1200
    aOrbit.push_back( d_iter ); // the last one with opposite orientation
}
Maciej Suminski's avatar
Maciej Suminski committed
1201

1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214
/** 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 )
{
Maciej Suminski's avatar
Maciej Suminski committed
1215
    // Two copies of the same dart
1216 1217 1218 1219 1220 1221 1222 1223 1224
    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 );
Maciej Suminski's avatar
Maciej Suminski committed
1225 1226
    }

1227 1228 1229 1230
    for( ;; )
    {
        if( d_iter == aD1 )
            return true;
Maciej Suminski's avatar
Maciej Suminski committed
1231

1232
        d_iter.Alpha1();
Maciej Suminski's avatar
Maciej Suminski committed
1233

1234 1235 1236 1237 1238 1239 1240 1241
        if( d_iter == aD1 )
            return true;

        d_iter.Alpha2();

        if( d_iter == d_end )
            break;
    }
Maciej Suminski's avatar
Maciej Suminski committed
1242 1243

    return false;
1244
}
Maciej Suminski's avatar
Maciej Suminski committed
1245

1246 1247 1248 1249 1250 1251 1252 1253
/** 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;
Maciej Suminski's avatar
Maciej Suminski committed
1254
    
1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268
    // (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;
Maciej Suminski's avatar
Maciej Suminski committed
1269
    
1270 1271 1272
    return ( d_iter == aD1 || d_iter.Alpha0() == aD1 || d_iter.Alpha1() == aD1 ||
            d_iter.Alpha0() == aD1 || d_iter.Alpha1() == aD1 || d_iter.Alpha0() == aD1 );
}
Maciej Suminski's avatar
Maciej Suminski committed
1273

1274 1275 1276 1277 1278 1279
// 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
Maciej Suminski's avatar
Maciej Suminski committed
1280

1281 1282 1283
    DART_TYPE d1 = aDart;
    DART_TYPE d2 = d1;
    d2.Alpha1();
Maciej Suminski's avatar
Maciej Suminski committed
1284
    
1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301
    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 )
{
Maciej Suminski's avatar
Maciej Suminski committed
1302
    // How "safe" is it?
1303 1304

    if( IsBoundaryEdge( aDart ) )
Maciej Suminski's avatar
Maciej Suminski committed
1305
        return false;
1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316

    // "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;
Maciej Suminski's avatar
Maciej Suminski committed
1317
    }
1318 1319 1320 1321
    else
    {
        if( TRAITS_TYPE::CrossProduct2D( d1, d2 ) <= 0.0 )
            return false;
Maciej Suminski's avatar
Maciej Suminski committed
1322
    }
1323

Maciej Suminski's avatar
Maciej Suminski committed
1324
    // Opposite side (still angle at the diagonal)
1325 1326
    d1 = aDart;
    d1.Alpha0();
Maciej Suminski's avatar
Maciej Suminski committed
1327
    d2 = d1;
1328 1329 1330 1331 1332 1333 1334
    d1.Alpha1();
    d2.Alpha2().Alpha1();

    if( aAllowDegeneracy )
    {
        if( TRAITS_TYPE::CrossProduct2D( d1, d2 ) < 0.0 )
            return false;
Maciej Suminski's avatar
Maciej Suminski committed
1335
    }
1336 1337 1338 1339
    else
    {
        if( TRAITS_TYPE::CrossProduct2D( d1, d2 ) <= 0.0 )
            return false;
Maciej Suminski's avatar
Maciej Suminski committed
1340 1341
    }

1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359
    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;
Maciej Suminski's avatar
Maciej Suminski committed
1360 1361 1362

    // If alpha2(d)=d, then boundary

1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385
    //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 );
Maciej Suminski's avatar
Maciej Suminski committed
1386 1387

    int no;
1388 1389 1390
    no = (int) blist.size();
    typename std::list<DART_TYPE>::const_iterator bit = blist.begin();
    DART_TYPE d1 = *bit;
Maciej Suminski's avatar
Maciej Suminski committed
1391
    ++bit;
1392
    DART_TYPE d2;
Maciej Suminski's avatar
Maciej Suminski committed
1393
    bool convex = true;
1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407

    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;
Maciej Suminski's avatar
Maciej Suminski committed
1408
    }
1409

Maciej Suminski's avatar
Maciej Suminski committed
1410 1411
    // Check the last angle
    d2 = *blist.begin();
1412 1413 1414 1415 1416 1417
    double crossProd = TRAITS_TYPE::CrossProduct2D( d1, d2 );

    if( crossProd < 0.0 )
    {
        //cout << "!!! Boundary is NOT convex: crossProd = " << crossProd << endl;
        convex = false;
Maciej Suminski's avatar
Maciej Suminski committed
1418
    }
1419

Maciej Suminski's avatar
Maciej Suminski committed
1420 1421 1422 1423
    //if (convex)
    //  cout << "\n---> Boundary is convex\n" << endl;
    //cout << endl;
    return convex;
1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452
}

//@} // 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() );
}
Maciej Suminski's avatar
Maciej Suminski committed
1453

1454 1455 1456 1457 1458
//------------------------------------------------------------------------------------------------
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 )
{
Maciej Suminski's avatar
Maciej Suminski committed
1459 1460 1461 1462 1463 1464
    // 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
1465

Maciej Suminski's avatar
Maciej Suminski committed
1466 1467 1468 1469
    // 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.
1470

Maciej Suminski's avatar
Maciej Suminski committed
1471 1472
    // Requires that swap retuns a dart in the "same position when rotated CCW"
    // (A vector instead of a list may be better.)
1473

Maciej Suminski's avatar
Maciej Suminski committed
1474
    // First check that elist is not empty
1475
    if( aElist.empty() )
Maciej Suminski's avatar
Maciej Suminski committed
1476 1477 1478 1479 1480
        return;

    // Avoid cycling by more extensive circumcircle test
    bool cycling_check = true;
    bool optimal = false;
1481
    typename DART_LIST_TYPE::iterator it;
Maciej Suminski's avatar
Maciej Suminski committed
1482

1483
    typename DART_LIST_TYPE::iterator end_opt = aEnd;
Maciej Suminski's avatar
Maciej Suminski committed
1484 1485 1486 1487 1488 1489 1490

    // 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"
1491
    //typename DART_LIST_TYPE::iterator end_opt;
Maciej Suminski's avatar
Maciej Suminski committed
1492 1493 1494 1495 1496
    //if (*end != NULL)
    //  end_opt = end;
    //else
    //  end_opt = elist.end();

1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513
    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
Maciej Suminski's avatar
Maciej Suminski committed
1514
    } // end pass
1515
}
Maciej Suminski's avatar
Maciej Suminski committed
1516

1517 1518 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530
/** 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>
Maciej Suminski's avatar
Maciej Suminski committed
1531
#if ((_MSC_VER > 0) && (_MSC_VER < 1300))//#ifdef _MSC_VER
1532 1533
bool TRIANGULATION_HELPER::SwapTestDelaunay(const DART_TYPE& aDart, bool aCyclingCheck = false) const
{
Maciej Suminski's avatar
Maciej Suminski committed
1534
#else
1535 1536
bool TRIANGULATION_HELPER::SwapTestDelaunay( const DART_TYPE& aDart, bool aCyclingCheck ) const
{
Maciej Suminski's avatar
Maciej Suminski committed
1537 1538 1539 1540 1541 1542 1543
#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.
1544 1545 1546 1547 1548 1549 1550 1551 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562

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

Maciej Suminski's avatar
Maciej Suminski committed
1563 1564 1565 1566 1567
    // "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.
1568
    if( cos1 >= 0 && cos2 >= 0 ) // both angles are grater or equual 90
Maciej Suminski's avatar
Maciej Suminski committed
1569
        return false;
1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581 1582 1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594 1595 1596 1597 1598 1599 1600 1601 1602 1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 1613 1614 1615 1616 1617 1618

    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;
        }
Maciej Suminski's avatar
Maciej Suminski committed
1619
    }
1620

Maciej Suminski's avatar
Maciej Suminski committed
1621
    return true;
1622 1623 1624 1625 1626 1627 1628 1629 1630 1631 1632 1633 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656
}

//-----------------------------------------------------------------------
//
//        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;

Maciej Suminski's avatar
Maciej Suminski committed
1657
    // Get the other "edges" of the current triangle; see illustration above.
1658 1659
    DART_TYPE oppEdge1 = aDiagonal;
    oppEdge1.Alpha1();
Maciej Suminski's avatar
Maciej Suminski committed
1660
    bool b1;
1661 1662 1663 1664 1665 1666 1667

    if( IsBoundaryEdge( oppEdge1 ) )
        b1 = true;
    else
    {
        b1 = false;
        oppEdge1.Alpha2();
Maciej Suminski's avatar
Maciej Suminski committed
1668
    }
1669 1670 1671

    DART_TYPE oppEdge2 = aDiagonal;
    oppEdge2.Alpha0().Alpha1().Alpha0();
Maciej Suminski's avatar
Maciej Suminski committed
1672
    bool b2;
1673 1674 1675 1676 1677 1678 1679

    if( IsBoundaryEdge( oppEdge2 ) )
        b2 = true;
    else
    {
        b2 = false;
        oppEdge2.Alpha2();
Maciej Suminski's avatar
Maciej Suminski committed
1680
    }
1681

Maciej Suminski's avatar
Maciej Suminski committed
1682
    // Swap the given diagonal
1683
    m_triangulation.swapEdge( aDiagonal );
Maciej Suminski's avatar
Maciej Suminski committed
1684

1685 1686 1687 1688 1689 1690 1691 1692 1693 1694 1695 1696 1697 1698 1699 1700 1701 1702 1703 1704 1705 1706 1707 1708 1709 1710 1711 1712 1713 1714 1715 1716 1717 1718 1719 1720
    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 )
{
Maciej Suminski's avatar
Maciej Suminski committed
1721 1722

    // Same iteration as in fixEdgesAtCorner, but not boundary    
1723 1724
    DART_TYPE dnext = aDart;

Maciej Suminski's avatar
Maciej Suminski committed
1725 1726 1727 1728 1729 1730 1731 1732
    // 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;
1733 1734 1735 1736 1737 1738 1739 1740 1741 1742 1743 1744 1745 1746 1747 1748 1749 1750 1751 1752 1753

    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--;
        }
Maciej Suminski's avatar
Maciej Suminski committed
1754 1755
    }

1756 1757 1758
    // Output, incident to the node
    aDart = dnext;
}
Maciej Suminski's avatar
Maciej Suminski committed
1759

1760 1761 1762 1763 1764 1765 1766 1767 1768 1769 1770 1771 1772 1773 1774 1775 1776 1777 1778 1779 1780 1781 1782
/** 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 )
{
Maciej Suminski's avatar
Maciej Suminski committed
1783 1784 1785 1786 1787 1788 1789 1790 1791
    // 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
1792

1793
    //int degree = getDegreeOfNode(dart);
1794 1795 1796 1797 1798 1799 1800 1801 1802 1803 1804 1805 1806 1807 1808 1809 1810 1811 1812 1813 1814 1815 1816 1817 1818 1819 1820 1821 1822 1823 1824 1825 1826 1827 1828 1829 1830 1831 1832 1833 1834 1835 1836 1837 1838 1839 1840 1841 1842 1843 1844 1845 1846 1847 1848 1849 1850 1851 1852 1853

    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
    }
Maciej Suminski's avatar
Maciej Suminski committed
1854
    else
1855 1856
    {
        d_iter.Alpha1(); // CW and see below
Maciej Suminski's avatar
Maciej Suminski committed
1857
    }
1858
    PositionAtNextBoundaryEdge( d_iter ); // CCW
Maciej Suminski's avatar
Maciej Suminski committed
1859

1860 1861 1862 1863 1864 1865 1866 1867 1868 1869 1870 1871 1872 1873 1874 1875 1876 1877 1878 1879 1880 1881 1882 1883 1884 1885 1886 1887 1888 1889 1890 1891 1892
    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();
Maciej Suminski's avatar
Maciej Suminski committed
1893 1894 1895 1896 1897 1898 1899 1900 1901 1902

    // 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?
1903 1904 1905 1906 1907 1908
    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 );
Maciej Suminski's avatar
Maciej Suminski committed
1909
    // Update the current dart which may have changed
1910 1911
    *aIt = dart;

Maciej Suminski's avatar
Maciej Suminski committed
1912 1913
    // Update darts that may have changed again (if they were present)
    // Note that dart is delivered back after swapping
1914 1915 1916 1917 1918
    if( it1 != aElist.end() )
    {
        d1 = dart;
        d1.Alpha1().Alpha0();
        *it1 = d1;
Maciej Suminski's avatar
Maciej Suminski committed
1919
    }
1920 1921 1922 1923 1924 1925

    if( it2 != aElist.end() )
    {
        d2 = dart;
        d2.Alpha2().Alpha1();
        *it2 = d2;
Maciej Suminski's avatar
Maciej Suminski committed
1926
    }
1927 1928 1929 1930 1931 1932

    if( it3 != aElist.end() )
    {
        d3 = dart;
        d3.Alpha2().Alpha1().Alpha0().Alpha1();
        *it3 = d3;
Maciej Suminski's avatar
Maciej Suminski committed
1933
    }
1934 1935 1936 1937 1938 1939

    if( it4 != aElist.end() )
    {
        d4 = dart;
        d4.Alpha0().Alpha1();
        *it4 = d4;
Maciej Suminski's avatar
Maciej Suminski committed
1940
    }
1941 1942 1943 1944 1945 1946
}

//@} // End of Utilities for Delaunay Triangulation Group

}
// End of ttl namespace scope (but other files may also contain functions for ttl)
Maciej Suminski's avatar
Maciej Suminski committed
1947 1948

#endif // _TTL_H_