haskal
/
kicad
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
kicad/include/ttl/ttl_constr.h

636 lines
24 KiB
C++
Raw Blame History

/*
* 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
/** \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
* <20>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
class ConstrainedTriangulation
{
public:
//------------------------------------------------------------------------------------------------
/* 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
* same_0_orbit
*/
template <class DartType>
static bool isTheConstraint(const DartType& dart, const DartType& dstart, const DartType& dend) {
DartType d0 = dart;
d0.alpha0(); // CW
if ((ttl::TriangulationHelper::same_0_orbit(dstart, dart) && ttl::TriangulationHelper::same_0_orbit(dend, d0)) ||
(ttl::TriangulationHelper::same_0_orbit(dstart, d0) && ttl::TriangulationHelper::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>
static 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
* isBoundaryNode
* positionAtNextBoundaryEdge
* TraitsType::orient2d
*/
template <class TraitsType, class DartType>
static 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::TriangulationHelper::isBoundaryNode(d_iter)) {
d_iter.alpha1(); // CW
ttl::TriangulationHelper::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>
static 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
* swappableEdge
* swapEdgeInList
* crossesConstraint
* isTheConstraint
*/
template <class TraitsType, class DartType>
void transformToConstraint(ttl::TriangulationHelper helper, DartType& dstart, DartType& dend,
std::list<DartType>& elist) const {
typename std::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::TriangulationHelper::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
helper.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 ConstrainedTriangulation class
}; // 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
* - 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 TriangulationHelper::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
std::list<DartType> elist;
DartType next_start = ttl_constr::ConstrainedTriangulation::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 (!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::ConstrainedTriangulation::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 std::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;
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_
Reference in New Issue View Git Blame Copy Permalink