ttl_constr.h 23.3 KB
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 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 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 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 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 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 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 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632
/*
 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT,
 * Applied Mathematics, Norway.
 *
 * Contact information: E-mail: tor.dokken@sintef.no                      
 * SINTEF ICT, Department of Applied Mathematics,                         
 * P.O. Box 124 Blindern,                                                 
 * 0314 Oslo, Norway.                                                     
 *
 * This file is part of TTL.
 *
 * TTL is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version. 
 *
 * TTL is distributed in the hope that it will be useful,        
 * but WITHOUT ANY WARRANTY; without even the implied warranty of         
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          
 * GNU Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public
 * License along with TTL. If not, see
 * <http://www.gnu.org/licenses/>.
 *
 * In accordance with Section 7(b) of the GNU Affero General Public
 * License, a covered work must retain the producer line in every data
 * file that is created or manipulated using TTL.
 *
 * Other Usage
 * You can be released from the requirements of the license by purchasing
 * a commercial license. Buying such a license is mandatory as soon as you
 * develop commercial activities involving the TTL library without
 * disclosing the source code of your own applications.
 *
 * This file may be used in accordance with the terms contained in a
 * written agreement between you and SINTEF ICT. 
 */

#ifndef _TTL_CONSTR_H_
#define _TTL_CONSTR_H_


#include <list>
#include <cmath>


// Debugging
#ifdef DEBUG_TTL_CONSTR_PLOT
  #include <fstream>
  static ofstream ofile_constr("qweCons.dat");
#endif


//using namespace std;

/** \brief Constrained Delaunay triangulation
*
*   Basic generic algorithms in TTL for inserting a constrained edge between two existing nodes.\n
*
*   See documentation for the namespace ttl for general requirements and assumptions.
*
*   \author 
*   Øyvind Hjelle, oyvindhj@ifi.uio.no
*/

namespace ttl_constr {

  //  ??? A constant used to evluate a numerical expression against a user spesified
  //  roundoff-zero number
#ifdef DEBUG_TTL_CONSTR
  static const double ROUNDOFFZERO = 0.0; // 0.1e-15;
#endif


  //------------------------------------------------------------------------------------------------
  /* Checks if \e dart has start and end points in \e dstart and \e dend.
  *
  *   \param dart 
  *   The dart that should be controlled to see if it's the constraint
  *
  *   \param dstart 
  *   A CCW dart with the startnode of the constraint as the startnode
  *
  *   \param dend 
  *   A CCW dart with the endnode of the constraint as the startnode
  *
  *   \retval bool 
  *   A bool confirming that it's the constraint or not 
  *
  *   \using
  *   ttl::same_0_orbit
  */
  template <class DartType>
    bool isTheConstraint(const DartType& dart, const DartType& dstart, const DartType& dend) {
    DartType d0 = dart;
    d0.alpha0(); // CW
    if ((ttl::same_0_orbit(dstart, dart) && ttl::same_0_orbit(dend, d0)) ||
      (ttl::same_0_orbit(dstart, d0)  && ttl::same_0_orbit(dend, dart))) {
      return true;
    }
    return false;
  }


  //------------------------------------------------------------------------------------------------
  /* Checks if \e d1 and \e d2 are on the same side of the line between \e dstart and \e dend.
  *   (The start nodes of \e d1 and \e d2 represent an edge).
  *
  *   \param dstart
  *   A CCW dart with the start node of the constraint as the source node of the dart.
  *
  *   \param dend
  *    A CCW dart with the end node of the constraint as the source node of the dart.
  * 
  *   \param d1
  *    A CCW dart with the first node as the start node of the dart.
  * 
  *   \param d2
  *   A CCW dart with the other node as the start node of the dart.
  *
  *   \using
  *   TraitsType::orient2d
  */
  template <class TraitsType, class DartType>
    bool crossesConstraint(DartType& dstart, DartType& dend, DartType& d1, DartType& d2) {
    
    typename TraitsType::real_type orient_1 = TraitsType::orient2d(dstart,d1,dend);
    typename TraitsType::real_type orient_2 = TraitsType::orient2d(dstart,d2,dend);
    // ??? Should we refine this? e.g. find if (dstart,dend) (d1,d2) represent the same edge
    if ((orient_1 <= 0 && orient_2 <= 0) || (orient_1 >= 0 && orient_2 >= 0))
      return false;
    
    return true;
  }


  //------------------------------------------------------------------------------------------------
  /* Return the dart \e d making the smallest non-negative angle,
  *   as calculated with: orient2d(dstart, d.alpha0(), dend),
  *   at the 0-orbit of dstart.
  *   If (dstart,dend) is a CCW boundary edge \e d will be CW, otherwise CCW (since CCW in)
  *   at the 0-orbit of dstart.
  * 
  *   \par Assumes:
  *   - CCW dstart and dend, but returned dart can be CW at the boundary.
  *   - Boundary is convex?
  *
  *   \param dstart
  *   A CCW dart dstart
  *
  *   \param dend 
  *   A CCW dart dend
  *
  *   \retval DartType
  *   The dart \e d making the smallest positive (or == 0) angle
  *
  *   \using
  *   ttl::isBoundaryNode
  *   ttl::positionAtNextBoundaryEdge
  *   TraitsType::orient2d
  */
  template <class TraitsType, class DartType>
    DartType getAtSmallestAngle(const DartType& dstart, const DartType& dend) {
    
    // - Must boundary be convex???
    // - Handle the case where the constraint is already present???
    // - Handle dstart and/or dend at the boundary
    //   (dstart and dend may define a boundary edge)
    
    DartType d_iter = dstart;
    if (ttl::isBoundaryNode(d_iter)) {
      d_iter.alpha1(); // CW
      ttl::positionAtNextBoundaryEdge(d_iter); // CCW (was rotated CW to the boundary)
    }
    
    // assume convex boundary; see comments
    
    DartType d0 = d_iter;
    d0.alpha0();
    bool ccw = true; // the rotation later
    typename TraitsType::real_type o_iter = TraitsType::orient2d(d_iter, d0, dend);
    if (o_iter == 0) { // collinear BUT can be on "back side"
      d0.alpha1().alpha0(); // CW
      if (TraitsType::orient2d(dstart, dend, d0) > 0)
        return d_iter; //(=dstart) collinear
      else {
        // collinear on "back side"
        d_iter.alpha1().alpha2(); // assume convex boundary
        ccw = true;
      }
    }
    else if (o_iter < 0) {
      // Prepare for rotating CW and with d_iter CW
      d_iter.alpha1();
      ccw = false;
    }
    
    // Set first angle
    d0 = d_iter; d0.alpha0();
    o_iter = TraitsType::orient2d(dstart, d0, dend);
    
    typename TraitsType::real_type o_next;
    
    // Rotate towards the constraint CCW or CW.
    // Here we assume that the boundary is convex.
    DartType d_next = d_iter;
    for (;;) {
      d_next.alpha1(); // CW !!! (if ccw == true)
      d0 = d_next; d0.alpha0();
      o_next = TraitsType::orient2d(dstart, d0, dend);
      
      if (ccw && o_next < 0) // and o_iter > 0
        return d_iter;
      else if (!ccw && o_next > 0)
        return d_next; // CCW
      else if (o_next == 0) {
        if (ccw)
          return d_next.alpha2(); // also ok if boundary
        else
          return d_next;
      }
      
      // prepare next
      d_next.alpha2(); // CCW if ccw
      d_iter = d_next; // also ok if boundary CCW if ccw == true
    }
  }


  //------------------------------------------------------------------------------------------------
  /* This function finds all the edges in the triangulation crossing
  *   the spesified constraint and puts them in a list.
  *   In the case of collinearity, an attempt is made to detect this.
  *   The first collinear node between dstart and dend is then returned.
  *
  *   Strategy:
  *   - Iterate such that \e d_iter is always strictly "below" the constraint
  *     as seen with \e dstart to the left and \e dend to the right.
  *   - Add CCW darts, whose edges intersect the constrait, to a list.
  *     These edges are found by the orient2d predicate:
  *     If two nodes of an edge are on opposite sides of the constraint,
  *     the edge between them intersect.
  *   - Must handle collinnear cases, i.e., if a node falls on the constraint,
  *     and possibly restarting collection of edges. Detecting collinearity
  *     heavily relies on the orient2d predicate which is provided by the
  *     traits class.
  *
  *   Action:
  *   1) Find cone/opening angle containing \e dstart and \e dend
  *   2) Find first edge from the first 0-orbit that intersects
  *   3) Check which of the two opposite that intersects
  *
  *   1)  
  *   Rotate CCW and find the (only) case where \e d_iter and \e d_next satisfy:
  *   - orient2d(d_iter, d_iter.alpha0(), dend) > 0
  *   - orient2d(d_next, d_next.alpha0(), dend) < 0
  *      
  *   - check if we are done, i.e., if (d_next.alpha0() == my_dend)
  *   - Note also the situation if, e.g., the constraint is a boundary edge in which case
  *     \e my_dend wil be CW
  *
  *   \param dstart 
  *   A CCW dart with the startnode of the constraint as the startnode
  *
  *   \param dend
  *   A CCW dart with the endnode of the constraint as the startnode
  *
  *   \param elist
  *   A list where all the edges crossing the spesified constraint will be put
  * 
  *   \retval dartType 
  *   Returns the next "collinear" starting node such that dend is returned when done.
  */
  template <class TraitsType, class DartType, class ListType>
    DartType findCrossingEdges(const DartType& dstart, const DartType& dend, ListType& elist) {
    
    const DartType my_start = getAtSmallestAngle<TraitsType>(dstart, dend);
    DartType my_end   = getAtSmallestAngle<TraitsType>(dend, dstart);
    
    DartType d_iter = my_start;
    if (d_iter.alpha0().alpha2() == my_end)
      return d_iter; // The constraint is an existing edge and we are done
    
    // Facts/status so far:
    // - my_start and my_end are now both CCW and the constraint is not a boundary edge.
    // - Further, the constraint is not one single existing edge, but it might be a collection
    //   of collinear edges in which case we return the current collinear edge
    //   and calling this function until all are collected.
    
    my_end.alpha1(); // CW! // ??? this is probably ok for testing now?
    
    d_iter = my_start;
    d_iter.alpha0().alpha1(); // alpha0 is downwards or along the constraint
    
    // Facts:
    // - d_iter is guaranteed to intersect, but can be in start point.
    // - d_iter.alpha0() is not at dend yet
    typename TraitsType::real_type orient = TraitsType::orient2d(dstart, d_iter, dend);

    // Use round-off error/tolerance or rely on the orient2d predicate ???
    // Make a warning message if orient != exact 0
    if (orient == 0)
      return d_iter;

#ifdef DEBUG_TTL_CONSTR
    else if (fabs(orient) <= ROUNDOFFZERO) {
      cout << "The darts are not exactly colinear, but |d1 x d2| <= " << ROUNDOFFZERO << endl;
      return d_iter; // collinear, not done (and not collect in the list)
    }
#endif

    // Collect intersecting edges
    // --------------------------
    elist.push_back(d_iter); // The first with interior intersection point
    
    // Facts, status so far:
    // - The first intersecting edge is now collected
    // (- d_iter.alpha0() is still not at dend)
    
    // d_iter should always be the edge that intersects and be below or on the constraint
    // One of the two edges opposite to d_iter must intersect, or we have collinearity
    
    // Note: Almost collinear cases can be handled on the
    // application side with orient2d. We should probably
    // return an int and the application will set it to zero
    for(;;) {
      // assume orient have been calc. and collinearity has been tested,
      // above the first time and below later
      d_iter.alpha2().alpha1(); // 2a same node
      
      DartType d0 = d_iter;
      d0.alpha0(); // CW
      if (d0 == my_end)
        return dend; // WE ARE DONE (but can we enter an endless loop???)
      
      // d_iter or d_iter.alpha0().alpha1() must intersect
      orient = TraitsType::orient2d(dstart, d0, dend);
      
      if (orient == 0) 
        return d0.alpha1();

#ifdef DEBUG_TTL_CONSTR
      else if (fabs(orient) <= ROUNDOFFZERO) {
        return d0.alpha1(); // CCW, collinear
      }
#endif

      else if (orient > 0) { // orient > 0 and still below
        // This one must intersect!
        d_iter = d0.alpha1();
      }
      elist.push_back(d_iter);
    }
  }


  //------------------------------------------------------------------------------------------------
  /* This function recives a constrained edge and a list of all the edges crossing a constraint.
  *   It then swaps the crossing edges away from the constraint. This is done according to a
  *   scheme suggested by Dyn, Goren & Rippa (slightly modified).
  *   The resulting triangulation is a constrained one, but not necessarily constrained Delaunay.
  *   In other to run optimization later to obtain a constrained Delaunay triangulation,
  *   the swapped edges are maintained in a list.
  *
  *   Strategy :
  *   - Situation A: Run through the list and swap crossing edges away from the constraint.
  *     All the swapped edges are moved to the end of the list, and are "invisible" to this procedure.
  *   - Situation B: We may come in a situation where none of the crossing edges can be swapped away
  *     from the constraint.
  *     Then we follow the strategy of Dyn, Goren & Rippa and allow edges to be swapped,
  *     even if they are not swapped away from the constraint.
  *     These edges are NOT moved to the end of the list. They are later swapped to none-crossing
  *     edges when the locked situation is solved.
  *   - We keep on swapping edges in Situation B until we have iterated trough the list.
  *     We then resume Situation A.
  *   - This is done until the list is virtualy empty. The resulting \c elist has the constraint
  *     as the last element.
  *
  *   \param dstart 
  *   A CCW dart dstart
  *
  *   \param dend
  *   A CCW dart dend
  *
  *   \param elist
  *   A list containing all the edges crossing the spesified constraint
  *
  *   \using
  *   ttl::swappableEdge
  *   ttl::swapEdgeInList
  *   ttl::crossesConstraint
  *   ttl::isTheConstraint
  */
  template <class TraitsType, class DartType>
    void transformToConstraint(DartType& dstart, DartType& dend, std::list<DartType>& elist) {
    
    typename list<DartType>::iterator it, used;
    
    // We may enter in a situation where dstart and dend are altered because of a swap.
    // (The general rule is that darts inside the actual quadrilateral can be changed,
    // but not those outside.)
    // So, we need some look-ahead strategies for dstart and dend and change these
    // after a swap if necessary.
    
    int dartsInList = (int)elist.size();
    if (dartsInList == 0)
      return;

    bool erase; // indicates if an edge is swapped away from the constraint such that it can be
                // moved to the back of the list
    bool locked = false;
    do {
      int noswap = 0;
      it = elist.begin();

      // counts how many edges that have been swapped per list-cycle
      int counter = 1;
      while(it != elist.end()) { // ??? change this test with counter > dartsInList
        erase = false;
        // Check if our virtual end of the list has been crossed. It breaks the
        // while and starts all over again in the do-while loop
        if (counter > dartsInList)
          break;
        
        if (ttl::swappableEdge<TraitsType, DartType>(*it, true)) {
          // Dyn & Goren & Rippa 's notation:
          // The node assosiated with dart *it is denoted u_m. u_m has edges crossing the constraint
          // named w_1, ... , w_r . The other node to the edge assosiated with dart *it is w_s.
          // We want to swap from edge u_m<->w_s to edge w_{s-1}<->w_{s+1}.
          DartType op1 = *it;
          DartType op2 = op1;
          op1.alpha1().alpha0(); //finds dart with node w_{s-1}
          op2.alpha2().alpha1().alpha0(); // (CW) finds dart with node w_{s+1}
          DartType tmp = *it; tmp.alpha0(); // Dart with assosiated node opposite to node of *it allong edge
          // If there is a locked situation we swap, even if the result is crossing the constraint
          // If there is a looked situation, but we do an ordinary swap, it should be treated as
          // if we were not in a locked situation!!

          // The flag swap_away indicates if the edge is swapped away from the constraint such that
          // it does not cross the constraint.
          bool swap_away = (crossesConstraint<TraitsType>(dstart, dend, *it, tmp) &&
                            !crossesConstraint<TraitsType>(dstart, dend, op1, op2));
          if (swap_away || locked) {
            // Do a look-ahead to see if dstart and/or dend are in the quadrilateral
            // If so, we mark it with a flag to make sure we update them after the swap
            // (they may have been changed during the swap according to the general rule!)
            bool start = false;
            bool end = false;
            
            DartType d = *it;
            if (d.alpha1().alpha0() == dstart)
              start = true;
            d = *it;
            if (d.alpha2().alpha1().alpha0().alpha1() == dend)
              end = true;
            
            // This is the only place swapping is called when inserting a constraint
            ttl::swapEdgeInList<TraitsType, DartType>(it,elist);
            
            // If we, during look-ahead, found that dstart and/or dend were in the quadrilateral,
            // we update them.
            if (end)
              dend = *it;
            if (start) {
              dstart = *it;
              dstart.alpha0().alpha2();
            }
            
            if (swap_away) { // !locked || //it should be sufficient with swap_away ???
              noswap++;
              erase = true;
            }
            
            if (isTheConstraint(*it, dstart, dend)) {
              // Move the constraint to the end of the list
              DartType the_constraint = *it;
              elist.erase(it);
              elist.push_back(the_constraint);
              return;
            } //endif
          } //endif
        } //endif "swappable edge"
        
                
        // Move the edge to the end of the list if it was swapped away from the constraint
        if (erase) {
          used = it;
          elist.push_back(*it);
          ++it;
          elist.erase(used);
          --dartsInList;
        } 
        else {
          ++it;
          ++counter;
        }
        
      } //end while
      
      if (noswap == 0)
        locked = true;
      
    } while (dartsInList != 0);
    
    
#ifdef DEBUG_TTL_CONSTR
    // We will never enter here. (If elist is empty, we return above).
    cout << "??????? ERROR 2, should never enter here ????????????????????????? SKIP ???? " << endl;
    exit(-1);
#endif

  }

}; // End of ttl_constr namespace scope


namespace ttl { // (extension)
  
  /** @name Constrained (Delaunay) Triangulation */
  //@{

  //------------------------------------------------------------------------------------------------
  /** Inserts a constrained edge between two existing nodes in a triangulation.
  *   If the constraint falls on one or more existing nodes and this is detected by the
  *   predicate \c TraitsType::orient2d, which should return zero in this case, the
  *   constraint is split. Otherwise a degenerate triangle will be made along
  *   the constraint. 
  *
  *   \param dstart 
  *   A CCW dart with the start node of the constraint as the source node
  *
  *   \param dend 
  *   A CCW dart with the end node of the constraint as the source node
  *
  *   \param optimize_delaunay
  *   If set to \c true, the resulting triangulation will be
  *   a \e constrained \e Delaunay \e triangulation. If set to \c false, the resulting
  *   triangulation will not necessarily be of constrained Delaunay type.
  *
  *   \retval DartType 
  *   A dart representing the constrained edge.
  *
  *   \require
  *   - \ref hed::TTLtraits::orient2d "TraitsType::orient2d" (DartType&, DartType&, PointType&)
  *   - \ref hed::TTLtraits::swapEdge "TraitsType::swapEdge" (DartType&)
  *
  *   \using
  *   - ttl::optimizeDelaunay if \e optimize_delaunay is set to \c true
  *
  *   \par Assumes:
  *   - The constrained edge must be inside the existing triangulation (and it cannot
  *     cross the boundary of the triangulation).
  */
  template <class TraitsType, class DartType>
    DartType insertConstraint(DartType& dstart, DartType& dend, bool optimize_delaunay) {
    
    // Assumes:
    // - It is the users responsibility to avoid crossing constraints
    // - The constraint cannot cross the boundary, i.e., the boundary must be
    //   convex in the area of crossing edges.
    // - dtart and dend are preserved (same node associated.)
    

    // Find edges crossing the constraint and put them in elist.
    // If findCrossingEdges reaches a Node lying on the constraint, this function 
    // calls itself recursively. 

    // RECURSION
    list<DartType> elist;
    DartType next_start = ttl_constr::findCrossingEdges<TraitsType>(dstart, dend, elist);

    // If there are no crossing edges (elist is empty), we assume that the constraint
    // is an existing edge.
    // In this case, findCrossingEdges returns the constraint.
    // Put the constraint in the list to fit with the procedures below
    // (elist can also be empty in the case of invalid input data (the constraint is in
    // a non-convex area) but this is the users responsibility.)

    //by Thomas Sevaldrud   if (elist.size() == 0)
    //by Thomas Sevaldrud   elist.push_back(next_start);

    // findCrossingEdges stops if it finds a node lying on the constraint.
    // A dart with this node as start node is returned
    // We call insertConstraint recursivly until the received dart is dend
    if (!ttl::same_0_orbit(next_start, dend)) {

#ifdef DEBUG_TTL_CONSTR_PLOT
      cout << "RECURSION due to collinearity along constraint" << endl;
#endif

      insertConstraint<TraitsType,DartType>(next_start, dend, optimize_delaunay);
    }
    
    // Swap edges such that the constraint edge is present in the transformed triangulation.
    if (elist.size() > 0) // by Thomas Sevaldrud
       ttl_constr::transformToConstraint<TraitsType>(dstart, next_start, elist);
    
#ifdef DEBUG_TTL_CONSTR_PLOT
    cout << "size of elist = " << elist.size() << endl;
    if (elist.size() > 0) {
      DartType the_constraint = elist.back();
      ofile_constr << the_constraint.x() << " " << the_constraint.y() << " " << 0 << endl;
      the_constraint.alpha0();
      ofile_constr << the_constraint.x() << " " << the_constraint.y() << " " << 0 << endl << endl;
    }
#endif

    // Optimize to constrained Delaunay triangulation if required.
    typename list<DartType>::iterator end_opt = elist.end();
    if (optimize_delaunay) {

      // Indicate that the constrained edge, which is the last element in the list,
      // should not be swapped
      --end_opt;
      ttl::optimizeDelaunay<TraitsType, DartType>(elist, end_opt);
    }

    if(elist.size() == 0) // by Thomas Sevaldrud
       return next_start; // by Thomas Sevaldrud
     
    // Return the constraint, which is still the last element in the list
    end_opt = elist.end();
    --end_opt;
    return *end_opt;
  }

  //@} // End of Constrained Triangulation Group

}; // End of ttl namespace scope (extension)

#endif // _TTL_CONSTR_H_