Library of Assembled Shared Sources
planar_mesh.h
Go to the documentation of this file.
1/** @file
2 * @author Bram de Greve (bram@cocamware.com)
3 * @author Tom De Muer (tom@cocamware.com)
4 *
5 * *** BEGIN LICENSE INFORMATION ***
6 *
7 * The contents of this file are subject to the Common Public Attribution License
8 * Version 1.0 (the "License"); you may not use this file except in compliance with
9 * the License. You may obtain a copy of the License at
10 * http://lass.sourceforge.net/cpal-license. The License is based on the
11 * Mozilla Public License Version 1.1 but Sections 14 and 15 have been added to cover
12 * use of software over a computer network and provide for limited attribution for
13 * the Original Developer. In addition, Exhibit A has been modified to be consistent
14 * with Exhibit B.
15 *
16 * Software distributed under the License is distributed on an "AS IS" basis, WITHOUT
17 * WARRANTY OF ANY KIND, either express or implied. See the License for the specific
18 * language governing rights and limitations under the License.
19 *
20 * The Original Code is LASS - Library of Assembled Shared Sources.
21 *
22 * The Initial Developer of the Original Code is Bram de Greve and Tom De Muer.
23 * The Original Developer is the Initial Developer.
24 *
25 * All portions of the code written by the Initial Developer are:
26 * Copyright (C) 2004-2024 the Initial Developer.
27 * All Rights Reserved.
28 *
29 * Contributor(s):
30 *
31 * Alternatively, the contents of this file may be used under the terms of the
32 * GNU General Public License Version 2 or later (the GPL), in which case the
33 * provisions of GPL are applicable instead of those above. If you wish to allow use
34 * of your version of this file only under the terms of the GPL and not to allow
35 * others to use your version of this file under the CPAL, indicate your decision by
36 * deleting the provisions above and replace them with the notice and other
37 * provisions required by the GPL License. If you do not delete the provisions above,
38 * a recipient may use your version of this file under either the CPAL or the GPL.
39 *
40 * *** END LICENSE INFORMATION ***
41 */
42
43
44
45
46/** @class lass::spat::PlanarMesh
47 * @brief a planar mesh
48 * @author Tom De Muer [TDM]
49 * @par Credits:
50 * This code is inspired on the code of Dani Lischinski distributed under the GPL (Graphics Gems IV)
51 */
52
53#ifndef LASS_GUARDIAN_OF_INCLUSION_SPAT_PLANAR_MESH_H
54#define LASS_GUARDIAN_OF_INCLUSION_SPAT_PLANAR_MESH_H
55
56#include "spat_common.h"
57#include "quad_edge.h"
58#include "../prim/point_2d.h"
59#include "../prim/aabb_2d.h"
60#include "../prim/ray_2d.h"
61#include "../prim/line_2d.h"
62#include "../prim/line_segment_2d.h"
63#include "../prim/simple_polygon_2d.h"
64#include "../prim/triangle_2d.h"
65#include "../prim/side.h"
66#include "../util/callback_r_1.h"
67#include "../util/callback_r_2.h"
68#include "../util/callback_r_3.h"
69#include "../io/matlab_o_stream.h"
70#include "../meta/tuple.h"
71#include "../stde/extended_iterator.h"
72
73#define DEBUG_MESH 0
74
75namespace lass
76{
77namespace spat
78{
79 #define TEMPLATE_DEF template< \
80 typename T, \
81 typename PointHandle, \
82 typename EdgeHandle, \
83 typename FaceHandle >
84
85 namespace impl
86 {
87 TEMPLATE_DEF class EdgeToMatlab;
88 TEMPLATE_DEF class EdgeGatherer;
89 TEMPLATE_DEF class EdgeMarker;
90 }
91
92 namespace experimental
93 {
94
95 struct ObjectAllocator:
96 private util::AllocatorThrow<
97 util::AllocatorBinned<
98 util::AllocatorSimpleBlock<
99 8192,
100 util::AllocatorFixed<
101 util::AllocatorAlignedAlloc<LASS_SIMD_ALIGNMENT>
102 >
103 >,
104 512,
105 util::BinnerPadded<LASS_SIMD_ALIGNMENT>,
106 util::AllocatorAlignedAlloc<LASS_SIMD_ALIGNMENT>
107 >
108 >
109 {
110 template <typename T> T* make(const T& x)
111 {
112 T* const p = static_cast<T*>(allocate(sizeof(T)));
113 try
114 {
115 new (p) T(x);
116 }
117 catch (...)
118 {
119 deallocate(p, sizeof(T));
120 throw;
121 }
122 return p;
123 }
124 template <typename T> void free(T* p)
125 {
126 if (p)
127 {
128 p->~T();
129 deallocate(p, sizeof(T));
130 }
131 }
132 };
133
134 template <typename T>
135 class BitField
136 {
137 public:
138 enum { size = 8 * sizeof(T) };
139 BitField(T bits = 0x0): bits_(bits) {}
140 bool operator[](size_t i) const
141 {
142 return bits_ & mask(i) ? true : false;
143 }
144 void set(size_t i)
145 {
146 bits_ |= mask(i);
147 }
148 void set(size_t i, bool value)
149 {
150 const T m = mask(i);
151 bits_ &= ~m;
152 bits_ |= (value ? m : 0);
153 }
154 void unset(size_t i)
155 {
156 bits_ &= ~mask(i);
157 }
158 private:
159 const T mask(size_t i) const
160 {
161 LASS_ASSERT(i < size);
162 return 1 << i;
163 }
164 T bits_;
165 };
166
167 template <typename T>
168 class ResetableThreadLocalVariable
169 {
170 public:
171 ResetableThreadLocalVariable(const T& proto = T()):
172 proto_(proto)
173 {
174 // *this must be fully initizalized before constructing Impl
175 TTls* tls = reinterpret_cast<TTls*>(tls_);
176 new (tls) TTls(Impl(this));
177 }
178 ~ResetableThreadLocalVariable()
179 {
180 TTls* tls = reinterpret_cast<TTls*>(tls_);
181 tls->~TTls();
182 }
183 T& operator*()
184 {
185 return *(*reinterpret_cast<TTls*>(tls_))->it;
186 }
187 const T& operator*() const
188 {
189 return *(*reinterpret_cast<TTls*>(tls_))->it;
190 }
191 void reset(const T& proto)
192 {
193 LASS_LOCK(mutex_)
194 {
195 proto_ = proto;
196 std::fill(values_.begin(), values_.end(), proto_);
197 }
198 }
199 private:
200 struct Impl
201 {
202 Impl(ResetableThreadLocalVariable* self): self(self)
203 {
204 this->self->registerImpl(this);
205 }
206 Impl(const Impl& other): self(other.self)
207 {
208 self->registerImpl(this);
209 }
210 ~Impl()
211 {
212 self->unregisterImpl(this);
213 }
214 ResetableThreadLocalVariable* self;
215 typename std::list<T>::iterator it;
216 };
217
218 typedef util::ThreadLocalVariable<Impl> TTls;
219
220 void registerImpl(Impl* impl)
221 {
222 LASS_LOCK(mutex_)
223 {
224 values_.push_back(proto_);
225 impl->it = stde::prev(values_.end());
226 }
227 }
228 void unregisterImpl(Impl* impl)
229 {
230 LASS_LOCK(mutex_)
231 {
232 values_.erase(impl->it);
233 }
234 }
235
236 char tls_[sizeof(TTls)];
237 std::list<T> values_;
238 T proto_;
239 util::Semaphore mutex_;
240 };
241
242 }
243
244 TEMPLATE_DEF
245 class PlanarMesh: private experimental::ObjectAllocator
246 {
247 public:
248 typedef lass::prim::Point2D<T> TPoint2D;
249 typedef lass::prim::Vector2D<T> TVector2D;
250 typedef lass::prim::Line2D<T, prim::Cartesian, prim::Unnormalized > TLine2D;
251 typedef lass::prim::Ray2D<T, prim::Unnormalized, prim::Unbounded > TRay2D;
252 typedef lass::prim::LineSegment2D<T> TLineSegment2D;
253 typedef lass::prim::SimplePolygon2D<T> TSimplePolygon2D;
254 typedef lass::prim::Triangle2D<T> TTriangle2D;
255 typedef experimental::BitField<num::Tuint32> TBitField;
256
257 static const int PLANAR_MESH_STACK_DEPTH = TBitField::size - 1; /**< this determines the maximum nesting depth of forAllVertices and forAllFaces */
258 enum
259 {
260 stackDepth = TBitField::size - 1,
261 publicMarkIndex = TBitField::size - 1
262 };
263 private:
264 struct ProxyHandle: private meta::Tuple< typename meta::type_list::Make<PointHandle, EdgeHandle, FaceHandle>::Type >
265 {
266 typedef meta::Tuple< typename meta::type_list::Make<PointHandle, EdgeHandle, FaceHandle>::Type > THandles;
267
268 TPoint2D* point_;
269 TBitField internalMark_;
270
271 ProxyHandle(): point_(NULL), internalMark_(0x00) {}
272 const PointHandle& pointHandle() const { return meta::tuple::field<0>(static_cast<const THandles&>(*this)); }
273 const EdgeHandle& edgeHandle() const { return meta::tuple::field<1>(static_cast<const THandles&>(*this)); }
274 const FaceHandle& faceHandle() const { return meta::tuple::field<2>(static_cast<const THandles&>(*this)); }
275 PointHandle& pointHandle() { return meta::tuple::field<0>(static_cast<THandles&>(*this)); }
276 EdgeHandle& edgeHandle() { return meta::tuple::field<1>(static_cast<THandles&>(*this)); }
277 FaceHandle& faceHandle() { return meta::tuple::field<2>(static_cast<THandles&>(*this)); }
278 };
279 class StackIncrementer
280 {
281 size_t* stackDepth_; /**< pointer to int containing real stack depth */
282 size_t maxStackDepth_; /**< maximum allowed stack depth */
283 StackIncrementer() {}
284 public:
285 StackIncrementer( size_t* iStackVar, size_t iMaxStackDepth ) : stackDepth_(iStackVar), maxStackDepth_( iMaxStackDepth )
286 {
287 ++(*stackDepth_);
288 if (*stackDepth_>maxStackDepth_)
289 {
290 --(*stackDepth_);
291 throw std::runtime_error("PlanarMesh: stack overflow");
292 }
293 }
294 ~StackIncrementer()
295 {
296 --(*stackDepth_);
297 }
298 };
299
300 public:
301 typedef lass::spat::QuadEdge<ProxyHandle> TQuadEdge;
302 typedef typename TQuadEdge::Edge TEdge;
303 typedef std::vector< TQuadEdge* > TQuadEdgeList;
304 typedef std::vector< ProxyHandle* > TProxyHandleList;
305 typedef lass::util::CallbackR1<bool,TEdge*> TEdgeCallback;
306 typedef lass::util::CallbackR3<bool,TEdge*,const TSimplePolygon2D&, FaceHandle> TEdgePolyFaceHandleCallback;
307
308 public:
309 PlanarMesh( const TPoint2D& a, const TPoint2D& b, const TPoint2D& c);
310 PlanarMesh( const TPoint2D& a, const TPoint2D& b, const TPoint2D& c, const TPoint2D& d);
311 PlanarMesh( const prim::Aabb2D<T>& iAabb );
312 void setTolerance(const T& iTolerance) {tolerance_ = iTolerance;} //<* set the relative tolerance */
313 const T& tolerance() { return tolerance_; }
314 void setPointDistanceTolerance(const T& iPointDistanceTolerance) {pointDistanceTolerance_ = iPointDistanceTolerance;} //<* set the relative tolerance */
315 const T& pointDistanceTolerance() { return pointDistanceTolerance_; }
316 virtual ~PlanarMesh();
317
318 TEdge* startEdge() const;
319 bool isBoundingPoint( const TPoint2D& iPoint) const; /**< true for points defining the boundary */
320
321 bool forAllPrimaryEdges( const TEdgeCallback& iCallback );
322 bool forAllPrimaryUndirectedEdges( const TEdgeCallback& iCallback );
323 bool forAllPrimaryUndirectedEdgesCached( const TEdgeCallback& iCallback );
324 bool forAllDualEdges( const TEdgeCallback& iCallback );
325 bool forAllEdges( const TEdgeCallback& iCallback );
326 bool forAllVertices( const TEdgeCallback& iCallback );
327 bool forAllFaces( const TEdgeCallback& iCallback );
328 bool forAllFacesCached( const TEdgeCallback& iCallback );
329 template <typename InputPolygonIterator, typename InputFaceHandleIterator>
330 bool forAllPolygonFaces( InputPolygonIterator iFirstPolygon, InputPolygonIterator iLastPolygon, InputFaceHandleIterator iFirstFaceHandle, const TEdgePolyFaceHandleCallback& iCallback );
331
332
333 TEdge* locate( const TPoint2D& iPoint, TEdge* iHint=0 ) const; /**< locate an edge of the triangle containing iPoint */
334 TEdge* pointLocate( const TPoint2D& iPoint ) const; /**< locate an edge of which the org is the same as iPoint, useful for known degeneracy point location, possibly slower than the regular locate */
335 TEdge* shoot( const TRay2D& iRay ) const; /**< locate the edge found by shooting the ray from within the triangle containt the tail of the ray */
336 template <typename OutputIterator> OutputIterator walk( const TLineSegment2D& iSegment, OutputIterator oCrossedEdges ) const;
337 template <typename OutputIterator> OutputIterator walkIntersections( const TLineSegment2D& iSegment, OutputIterator oIntersections ) const;
338 std::pair<int, TEdge*> walkTillConstrained( const TRay2D& iRay) const;
339 TEdge* insertSite( const TPoint2D& iPoint, bool makeDelaunay = true, bool forceOnEdge = false);
340 TEdge* insertEdge( const TLineSegment2D& iSegment, EdgeHandle iLeftHandle = EdgeHandle(), EdgeHandle iRightHandle = EdgeHandle(), PointHandle iPointHandle = PointHandle(), bool forcePointHandle = false, bool makeDelaunay = true);
341 TEdge* insertPolygon( const TSimplePolygon2D& iSegment, EdgeHandle iLeftHandle = EdgeHandle(), EdgeHandle iRightHandle = EdgeHandle(), bool makeDelaunay = true);
342 void markPolygon( TEdge* iStartEdge, const TSimplePolygon2D& iPolygon, FaceHandle iFaceHandle );
343 template <typename InputPolygonIterator, typename InputFaceHandleIterator>
344 void markPolygons( InputPolygonIterator iFirstPolygon, InputPolygonIterator iLastPolygon, InputFaceHandleIterator iFirstFaceHandle);
345 void markPolygons( const std::vector<TSimplePolygon2D>& iPolygons, const std::vector<FaceHandle>& iFaceHandles);
346 bool deleteEdge( TEdge* iEdge ); /**< delete edge without using gc */
347 bool gcDeleteEdge( TEdge* iEdge ); /**< delete edge using garbage collector, useful for deletion avalanches */
348 int gc(); /**< do garbage collection after deletion avalanches, returns number of quadedge collected */
349 long edgeCount() const { return edgeCount_; }
350 void makeMaximalConvexPolygon(T iMaxSurface=T(-1));
351 void makeRectangular(T minAngle, T maxAngle);
352
353 static TTriangle2D triangle( const TEdge* iEdge);
354 static TSimplePolygon2D polygon( const TEdge* iEdge);
355 static const TPoint2D& org( const TEdge* iEdge );
356 static const TPoint2D& dest( const TEdge* iEdge );
357 static const TPoint2D& fastOrg( const TEdge* iEdge );
358 static const TPoint2D& fastDest( const TEdge* iEdge );
359 static TPoint2D along( const TEdge* iEdge, const T& iParam );
360 static TPoint2D fastAlong( const TEdge* iEdge, const T& iParam );
361 static const TVector2D direction( const TEdge* iEdge );
362 static bool rightOf( const TPoint2D& iPoint, const TEdge* iEdge );
363 static bool fastRightOf( const TPoint2D& iPoint, const TEdge* iEdge );
364 static bool leftOf( const TPoint2D& iPoint, const TEdge* iEdge );
365 static bool fastLeftOf( const TPoint2D& iPoint, const TEdge* iEdge );
366 static bool onEdge( const TPoint2D& iPoint, const TEdge* iEdge );
367 static bool hasLeftFace( const TEdge* iEdge );
368 static bool fastHasLeftFace( const TEdge* iEdge );
369 static bool hasRightFace( const TEdge* iEdge );
370 static int chainOrder( const TEdge* iEdge );
371 static int vertexOrder( const TEdge* iEdge );
372 static bool allEqualChainOrder( const TEdge* iEdge );
373 static bool inConvexCell( const TPoint2D& iPoint, const TEdge* iEdge );
374
375 // set the temp flag if you want to modify the mesh during mesh cell/face/point iteration
376 // end your algorithm with a commitTempQuadEdges
377 TEdge* makeEdge( const TPoint2D& a, const TPoint2D& b, bool makeConstrained ); /**< makes a _dangling_ edge, this is really only for ADVANCED usage, you'll have to splice it or it won't be in the mesh! */
378 TEdge* connect( TEdge* a, TEdge* b ); /**< connects the dest of a with the org of b */
379 TEdge* split( TEdge* a, T iWhere ); /**< splits an edge with iWhere in the open interval (0,1) and returns the new edge with as org the splitting point */
380 void swap( TEdge* iEdge );
381 // when set to true any allocations are done in a buffer of quadedges, when setting again to false, the buffer is added
382 void setTempQuadEdges(bool iAllocateInTemp);
383 /* move the org of iEdge to iNewLocation
384 * The new location has to be in the "neighbourhood" which is defined by any
385 * cells that have an edge incident to the org of iEdge
386 */
387 // these need to be moved in a different way into the public interface!
388 static void setOrg( const TPoint2D& iOrg, TEdge* iEdge );
389 static void setDest( const TPoint2D& iDest, TEdge* iEdge );
390 bool moveWithinNeighbourhood( TEdge* iEdge, const TPoint2D& iNewLocation );
391 /* remove a vertex from the mesh
392 * no guarantees on any mesh properties afterwards, so be careful!
393 * the most imminent property that is removed is convexness, this shall be corrected later on!
394 * [TODO] keep convexness constraint
395 */
396 bool removeVertex(TEdge* iEdge);
397
398 static bool inPrimaryMesh( const TEdge* iEdge );
399 static bool inDualMesh( const TEdge* iEdge );
400 static bool marking( const TEdge* iEdge );
401 static void setMarking( TEdge* iEdge, bool iMark );
402 static PointHandle pointHandle( const TEdge* iEdge );
403 static EdgeHandle edgeHandle( const TEdge* iEdge );
404 static FaceHandle faceHandle( const TEdge* iEdge );
405 static PointHandle& pointHandleRef( TEdge* iEdge );
406 static EdgeHandle& edgeHandleRef( TEdge* iEdge );
407 static FaceHandle& faceHandleRef( TEdge* iEdge );
408 static void setPointHandle( TEdge* iEdge, PointHandle iHandle );
409 static void setEdgeHandle( TEdge* iEdge, EdgeHandle iHandle );
410 static void setFaceHandle( TEdge* iEdge, FaceHandle iHandle );
411 static void setOrientedEdgeHandle( TEdge* iEdge, EdgeHandle iLeftHandle, EdgeHandle iRightHandle, const TVector2D& iDirection );
412 private:
413 bool allocateInTemp_;
414 TEdge* startEdge_;
415 TQuadEdgeList quadEdgeList_;
416 TQuadEdgeList tempQuadEdgeList_;
417 std::vector<TPoint2D> boundingPoints_;
418 T tolerance_;
419 T pointDistanceTolerance_;
420 long edgeCount_;
421 size_t stackDepth_;
422
423
424 mutable experimental::ResetableThreadLocalVariable<TEdge*> lastLocateEdge_;
425 mutable experimental::ResetableThreadLocalVariable<TEdge*> lastFloodEdge_;
426
427 private:
428 PlanarMesh();
429 void init4( const TPoint2D& a, const TPoint2D& b, const TPoint2D& c, const TPoint2D& d);
430 void fixEdge( TEdge* e ); /**< in quadrilateral possibly switch the diagonal for delaunay */
431 TEdge* makeEmptyEdge( bool makeConstrained );
432 void triangulate( TEdge* iEdge );
433 void splitEdge(TEdge *e, const TPoint2D& iPoint );
434 TEdge* pointShoot( const TRay2D& iRay ) const; /**< locate the edge found by shooting the ray from within the triangle with the support of the ray as a known point */
435 template <typename OutputIterator> OutputIterator pointWalk( const TLineSegment2D& iSegment, OutputIterator oCrossedEdges ) const;
436 void safeSplitEdge(TEdge *e, const TPoint2D& iPoint );
437 TPoint2D snap(const TPoint2D& a, const TPoint2D& b, const TPoint2D& c);
438
439 TEdge* bruteForceLocate( const TPoint2D& iPoint ) const;
440 TEdge* bruteForceExactLocate( const TPoint2D& iPoint ) const;
441
442 friend class impl::EdgeToMatlab<T, PointHandle, EdgeHandle, FaceHandle>;
443 friend class impl::EdgeGatherer<T, PointHandle, EdgeHandle, FaceHandle>;
444 friend class impl::EdgeMarker<T, PointHandle, EdgeHandle, FaceHandle>;
445
446 bool internalMarking( const TEdge* iEdge );
447 bool deletePoint( TEdge* e);
448 void setInternalMarking( TEdge* iEdge, bool iMark );
449 void setInternalMarkingAroundVertex( typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* iEdge, bool iMark );
450 void setInternalMarkingInFace( typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* iEdge, bool iMark );
451
452 bool floodPolygon( TEdge* iStartEdge, const TSimplePolygon2D& iPolygon, FaceHandle iFaceHandle );
453 bool floodPolygonCallback( TEdge* iStartEdge, const TSimplePolygon2D& iPolygon, FaceHandle iFaceHandle, const TEdgePolyFaceHandleCallback& iCallback);
454#ifndef NDEBUG
455 public:
456 static unsigned numSetOrientedEdgeHandleCalls;
457 static unsigned numSetOrientedEdgeHandleSwaps;
458#endif
459 };
460
461 namespace impl
462 {
463 /** fastIntersect. Returns the intersection of line iAiB with iCiD. This routine is aimed
464 at speed and less at accuracy. This routines assumes there is an intersection point and won't try to be more specific about the error
465 when there is not. It is used for calculating intersections with walking paths in meshes
466 */
467 template <typename T>
468 lass::prim::Result fastIntersect( const lass::prim::Point2D<T>& iA, const lass::prim::Point2D<T>& iB,
469 const lass::prim::Point2D<T>& iC, const lass::prim::Point2D<T>& iD, lass::prim::Point2D<T>& oP)
470 {
471 const lass::prim::Vector2D<T> dirA=iB-iA;
472 const lass::prim::Vector2D<T> difference=iC-iA;
473 const lass::prim::Vector2D<T> dirB=iD-iC;
474
475 const T denominator = lass::prim::perpDot(dirA, dirB);
476 const T denominator2= denominator*2.0; // this is (and should be) optimised by adding denominator to itself
477 oP = iA; // interwoven instruction as we assume the zero equality is seldom encountered
478 if (denominator == denominator2)
479 {
480 // we don't bother finding out which case, we only need to know that
481 // the result is not what we are looking for
482 return lass::prim::rInvalid;
483 }
484 else
485 {
486 const T oTa = perpDot(difference, dirB) / denominator;
487 oP.x += dirA.x*oTa;
488 oP.y += dirA.y*oTa;
489 return lass::prim::rOne; // intersecting
490 }
491 }
492
493 TEMPLATE_DEF
494 class EdgeMarker
495 {
496 typedef PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle> TPlanarMesh;
497 TPlanarMesh* mesh_;
498 bool marking_;
499 public:
500 EdgeMarker( TPlanarMesh* iMesh, bool iMark ) : mesh_(iMesh), marking_(iMark) {}
501 /** internalMark. This changes internal marking. This is not user code, unless
502 * your really know what you are doing, use the mark instead */
503 bool internalMark( typename TPlanarMesh::TEdge* e )
504 {
505 mesh_->setInternalMarking( e, marking_ );
506 return true;
507 }
508 bool mark( typename TPlanarMesh::TEdge* e )
509 {
510 TPlanarMesh::setMarking( e, marking_ );
511 return true;
512 }
513
514 };
515
516 TEMPLATE_DEF
517 class EdgeToMatlab
518 {
519 public:
520 typedef PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle> TPlanarMesh;
521 private:
522 TPlanarMesh* mesh_;
523 lass::io::MatlabOStream& stream_;
524 public:
525 EdgeToMatlab( TPlanarMesh* iMesh, lass::io::MatlabOStream& ioOStream ) : mesh_(iMesh), stream_( ioOStream ) {}
526 bool edgeToMatlab( typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* iEdge )
527 {
528 if (!mesh_->internalMarking( iEdge ))
529 {
530 if (iEdge->isConstrained())
531 stream_.setColor(lass::io::mcRed);
532 else
533 stream_.setColor(lass::io::mcBlack);
534 stream_ << typename TPlanarMesh::TLineSegment2D(
535 TPlanarMesh::org(iEdge), TPlanarMesh::dest(iEdge) );
536 }
537 else
538 return true;
539 mesh_->setInternalMarking( iEdge, true );
540 mesh_->setInternalMarking( iEdge->sym(), true );
541 return true;
542 }
543 bool vertexToMatlab( typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* iEdge )
544 {
545 typedef PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle> TPlanarMesh;
546 if ( !mesh_->internalMarking( iEdge ) )
547 stream_ << TPlanarMesh::org(iEdge);
548 else
549 return true;
550 mesh_->setInternalMarkingAroundVertex( iEdge, true );
551 return true;
552 }
553 };
554
555 TEMPLATE_DEF
556 class EdgeGatherer
557 {
558 typedef PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle> TPlanarMesh;
559 TPlanarMesh* mesh_;
560 T lastArea_;
561 public:
562 typedef std::list< typename TPlanarMesh::TEdge* > TEdgeList;
563 TEdgeList edgeList;
564 T angleConstraintMin; /**< min angle for which pruning is considered */
565 T angleConstraintMax; /**< max angle for which pruning is considered */
566 T maxSurface; /**< max surface of convex polygon */
567
568 EdgeGatherer( TPlanarMesh* iMesh ) : mesh_(iMesh), lastArea_(0) {}
569 virtual ~EdgeGatherer() {}
570
571 bool makeConvexPolygon( typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* e )
572 {
573 LASS_ASSERT( TPlanarMesh::inPrimaryMesh( e ) );
574 if ( mesh_->internalMarking( e ) || e->isConstrained())
575 return true;
576
577 if ( prim::weakCcw( TPlanarMesh::org(e->dNext()), TPlanarMesh::dest(e), TPlanarMesh::dest(e->lNext()))
578 && prim::weakCcw( TPlanarMesh::org(e->sym()->dNext()), TPlanarMesh::org(e), TPlanarMesh::dest(e->sym()->lNext())) )
579 {
580 T dArea = num::abs(prim::doubleTriangleArea( TPlanarMesh::org(e->dNext()), TPlanarMesh::dest(e), TPlanarMesh::dest(e->lNext())))
581 + num::abs(prim::doubleTriangleArea( TPlanarMesh::org(e->sym()->dNext()), TPlanarMesh::org(e), TPlanarMesh::dest(e->sym()->lNext())));
582 if (maxSurface>T(0))
583 {
584 // use of maximum surface criterion
585 T leftArea = TPlanarMesh::polygon(e).area();
586 T rightArea = TPlanarMesh::polygon(e->sym()).area();
587 if (leftArea+rightArea>maxSurface)
588 return true; // skip this edge
589 }
590 if (lastArea_ == T(0))
591 {
592 lastArea_ = dArea;
593 edgeList.push_back( e );
594 }
595 else
596 {
597 if (dArea<lastArea_)
598 {
599 edgeList.push_front( e );
600 lastArea_ = dArea;
601 }
602 else
603 edgeList.push_back( e );
604 }
605 return false;
606 }
607 return true;
608 }
609 bool makeRectangular( typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* e )
610 {
611 LASS_ASSERT( TPlanarMesh::inPrimaryMesh( e ) );
612 if ( mesh_->internalMarking( e ) || e->isConstrained())
613 return true;
614
615 if ( prim::ccw( TPlanarMesh::org(e->dNext()), TPlanarMesh::dest(e), TPlanarMesh::dest(e->lNext()))
616 && prim::ccw( TPlanarMesh::org(e->sym()->dNext()), TPlanarMesh::org(e), TPlanarMesh::dest(e->sym()->lNext())) )
617 {
618 // angles will always be smaller than 180deg so we can use this test
619 typename TPlanarMesh::TVector2D v1 = TPlanarMesh::fastDest(e->lNext())-TPlanarMesh::fastDest(e);
620 typename TPlanarMesh::TVector2D v2 = TPlanarMesh::fastOrg(e)-TPlanarMesh::fastDest(e->lNext());
621 typename TPlanarMesh::TVector2D v3 = TPlanarMesh::fastDest(e->sym()->lNext())-TPlanarMesh::fastOrg(e);
622 typename TPlanarMesh::TVector2D v4 = TPlanarMesh::fastDest(e)-TPlanarMesh::fastDest(e->sym()->lNext());
623 v1.normalize();
624 v2.normalize();
625 v3.normalize();
626 v4.normalize();
627 T angle1 = lass::num::acos(lass::prim::dot( v1, v2));
628 T angle2 = lass::num::acos(lass::prim::dot( v2, v3));
629 T angle3 = lass::num::acos(lass::prim::dot( v3, v4));
630 T angle4 = lass::num::acos(lass::prim::dot( v4, v1));
631 if ( lass::num::abs(angle1)>angleConstraintMin
632 && lass::num::abs(angle2)>angleConstraintMin
633 && lass::num::abs(angle3)>angleConstraintMin
634 && lass::num::abs(angle4)>angleConstraintMin
635 && lass::num::abs(angle1)<=angleConstraintMax
636 && lass::num::abs(angle2)<=angleConstraintMax
637 && lass::num::abs(angle3)<=angleConstraintMax
638 && lass::num::abs(angle4)<=angleConstraintMax )
639 {
640 edgeList.push_back( e );
641 // mark the other connecting edges as hard constraints
642 //mesh_->setInternalMarking(e->dNext(),true);
643 //mesh_->setInternalMarking(e->lNext(),true);
644 return false;
645 }
646 }
647 return true;
648 }
649 };
650
651 TEMPLATE_DEF
652 class BrutePointExactLocator
653 {
654 typedef PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle> TPlanarMesh;
655 typename TPlanarMesh::TPoint2D point_;
656 public:
657 typename TPlanarMesh::TEdge* edge;
658 BrutePointExactLocator( const typename TPlanarMesh::TPoint2D& iPoint ) : point_(iPoint), edge(NULL) {}
659 bool findPoint( typename TPlanarMesh::TEdge* e )
660 {
661 if (TPlanarMesh::org(e)==point_)
662 {
663 edge = e;
664 return false;
665 }
666 return true;
667 }
668 };
669
670 TEMPLATE_DEF
671 class BrutePointLocator
672 {
673 typedef PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle> TPlanarMesh;
674 typename TPlanarMesh::TPoint2D point_;
675 public:
676 typename TPlanarMesh::TEdge* edge;
677 BrutePointLocator( const typename TPlanarMesh::TPoint2D& iPoint ) : point_(iPoint), edge(NULL) {}
678 bool findEdge( typename TPlanarMesh::TEdge* e )
679 {
680 typename TPlanarMesh::TEdge* ce = e;
681 for (int i=0;i<3;++i)
682 {
683 if (TPlanarMesh::rightOf(point_,ce))
684 return true;
685 ce = ce->lNext();
686 }
687 // now search the closest edge
688 ce = e;
689
690 if (point_==TPlanarMesh::org(e))
691 {
692 edge = e;
693 return false;
694 }
695 if (point_==TPlanarMesh::dest(e))
696 {
697 edge = e->sym();
698 return false;
699 }
700 if (point_==TPlanarMesh::dest(e->lNext()))
701 {
702 edge = e->lNext()->sym();
703 return false;
704 }
705
706 typename TPlanarMesh::TRay2D R1(TPlanarMesh::org(e), TPlanarMesh::dest(e));
707 typename TPlanarMesh::TPoint2D p1 = R1.point( R1.t( point_ ) );
708 typename TPlanarMesh::TRay2D R2(TPlanarMesh::dest(e), TPlanarMesh::org(e->lPrev()));
709 typename TPlanarMesh::TPoint2D p2 = R2.point( R2.t( point_ ) );
710 typename TPlanarMesh::TRay2D R3(TPlanarMesh::org(e->lPrev()), TPlanarMesh::org(e));
711 typename TPlanarMesh::TPoint2D p3 = R3.point( R3.t( point_ ) );
712
713 typename TPlanarMesh::TPoint2D::TValue d1 = squaredDistance(point_,p1);
714 typename TPlanarMesh::TPoint2D::TValue d2 = squaredDistance(point_,p2);
715 typename TPlanarMesh::TPoint2D::TValue d3 = squaredDistance(point_,p3);
716
717 if ((d1<d2) && (d1<d3))
718 {
719 edge = e;
720 return false;
721 }
722 if ((d2<d3) && (d2<=d1))
723 {
724 edge = e->lNext();
725 return false;
726 }
727 edge = e->lPrev();
728 return false;
729
730 /*
731 T d1 = squaredDistance(TPlanarMesh::org(e),point_);
732 T d2 = squaredDistance(TPlanarMesh::dest(e),point_);
733 T d3 = squaredDistance(TPlanarMesh::dest(e->lNext()),point_);
734
735 if (d1<d2)
736 {
737 if (d1<d3)
738 edge=e;
739 else
740 edge=e->lNext()->sym();
741 }
742 else
743 {
744 if (d2<d3)
745 edge=e->sym();
746 else
747 edge=e->lNext()->sym();
748 }
749 return false;
750 */
751 }
752 };
753
754 TEMPLATE_DEF
755 class BrutePointLocatorVerbose
756 {
757 typedef PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle> TPlanarMesh;
758 typename TPlanarMesh::TPoint2D point_;
759 public:
760 typename TPlanarMesh::TEdge* edge;
761 std::ofstream stream;
762 BrutePointLocatorVerbose( const typename TPlanarMesh::TPoint2D& iPoint ) : point_(iPoint), edge(NULL) {}
763 bool findEdge( typename TPlanarMesh::TEdge* e )
764 {
765 typename TPlanarMesh::TEdge* ce = e;
766#if DEBUG_MESH
767 stream << "---\n";
768 stream << std::setprecision(20);
769#endif
770 for (int i=0;i<3;++i)
771 {
772#if DEBUG_MESH
773 stream << TPlanarMesh::org(ce) << "-->" << TPlanarMesh::dest(ce) << ":" << prim::doubleTriangleArea(point_,TPlanarMesh::dest(ce),TPlanarMesh::org(ce)) << "\n";
774#endif
775#pragma LASS_TODO("Get rid of the epsilon!")
776 if (prim::doubleTriangleArea(point_,TPlanarMesh::dest(ce),TPlanarMesh::org(ce))>1e-12)
777 return true;
778
779 //if (TPlanarMesh::rightOf(point_,ce))
780 // return true;
781 ce = ce->lNext();
782 }
783 // now search the closest point
784 ce = e;
785 T d1 = squaredDistance(TPlanarMesh::org(e),point_);
786 T d2 = squaredDistance(TPlanarMesh::dest(e),point_);
787 T d3 = squaredDistance(TPlanarMesh::dest(e->lNext()),point_);
788
789 if (d1<d2)
790 {
791 if (d1<d3)
792 edge=e;
793 else
794 edge=e->lNext()->sym();
795 }
796 else
797 {
798 if (d2<d3)
799 edge=e->sym();
800 else
801 edge=e->lNext()->sym();
802 }
803 return false;
804 }
805 };
806
807 }
808
809 TEMPLATE_DEF
810 lass::io::MatlabOStream& operator<<( lass::io::MatlabOStream& ioOStream, const typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge& iEdge )
811 {
812 typedef PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle> TPlanarMesh;
813 ioOStream << typename TPlanarMesh::TLineSegment2D(
814 TPlanarMesh::org(iEdge), TPlanarMesh::dest(iEdge) )
815 << std::endl;
816 return ioOStream;
817 }
818
819 TEMPLATE_DEF
820 lass::io::MatlabOStream& operator<<( lass::io::MatlabOStream& ioOStream, PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>& iMesh )
821 {
822 typedef impl::EdgeMarker<T, PointHandle, EdgeHandle, FaceHandle> TEdgeMarker;
823 typedef PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle> TPlanarMesh;
824 typedef typename TPlanarMesh::TEdgeCallback TEdgeCallback;
825 impl::EdgeToMatlab<T,PointHandle,EdgeHandle,FaceHandle> matlabWriter(&iMesh, ioOStream);
826 TEdgeMarker edgeMarker( &iMesh, false );
827 iMesh.forAllPrimaryEdges( TEdgeCallback( &edgeMarker, &TEdgeMarker::internalMark ) );
828 iMesh.forAllPrimaryEdges( TEdgeCallback( &matlabWriter,
829 &impl::EdgeToMatlab<T,PointHandle,EdgeHandle,FaceHandle>::edgeToMatlab ) );
830 iMesh.forAllPrimaryEdges( TEdgeCallback( &edgeMarker, &TEdgeMarker::internalMark ) );
831 iMesh.forAllPrimaryEdges( TEdgeCallback( &matlabWriter,
832 &impl::EdgeToMatlab<T,PointHandle,EdgeHandle,FaceHandle>::vertexToMatlab ) );
833 return ioOStream;
834 }
835
836 TEMPLATE_DEF
837 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::makeMaximalConvexPolygon(T iMaxSurface)
838 {
839 typedef impl::EdgeGatherer<T, PointHandle, EdgeHandle, FaceHandle> TEdgeGatherer;
840 typedef impl::EdgeMarker<T, PointHandle, EdgeHandle, FaceHandle> TEdgeMarker;
841
842 TEdgeMarker edgeMarker( this, false );
843 forAllPrimaryEdges( TEdgeCallback( &edgeMarker, &TEdgeMarker::internalMark ) );
844 while (true)
845 {
846 TEdgeGatherer edgeGatherer( this );
847 edgeGatherer.maxSurface = iMaxSurface;
848 edgeGatherer.edgeList.clear();
849 forAllPrimaryUndirectedEdgesCached( TEdgeCallback( &edgeGatherer, &TEdgeGatherer::makeConvexPolygon) );
850 if (edgeGatherer.edgeList.size()>0)
851 gcDeleteEdge( edgeGatherer.edgeList.front() );
852 else
853 break;
854 }
855#if DEBUG_MESH
856 std::cout << "PlanarMesh GC collected " << gc() << " edges" << std::endl;
857#endif
858 }
859
860 TEMPLATE_DEF
861 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::makeRectangular(T minAngle, T maxAngle)
862 {
863 typedef impl::EdgeGatherer<T, PointHandle, EdgeHandle, FaceHandle> TEdgeGatherer;
864 typedef impl::EdgeMarker<T, PointHandle, EdgeHandle, FaceHandle> TEdgeMarker;
865
866 TEdgeMarker edgeMarker( this, false );
867 forAllPrimaryEdges( TEdgeCallback( &edgeMarker, &TEdgeMarker::internalMark ) );
868 while (true)
869 {
870 TEdgeGatherer edgeGatherer( this );
871 edgeGatherer.edgeList.clear();
872 edgeGatherer.angleConstraintMin = minAngle;
873 edgeGatherer.angleConstraintMax = maxAngle;
874 forAllPrimaryUndirectedEdgesCached( TEdgeCallback( &edgeGatherer, &TEdgeGatherer::makeRectangular) );
875 if (edgeGatherer.edgeList.size()==0)
876 break;
877 if (edgeGatherer.edgeList.size()>0)
878 {
879 gcDeleteEdge( edgeGatherer.edgeList.back() );
880 edgeGatherer.edgeList.pop_back();
881 }
882 }
883#if DEBUG_MESH
884 std::cout << "PlanarMesh GC collected " << gc() << " edges" << std::endl;
885#endif
886 }
887
888
889 TEMPLATE_DEF
890 PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::PlanarMesh( const TPoint2D& a, const TPoint2D& b, const TPoint2D& c)
891 {
892 allocateInTemp_ = false;
893 tolerance_ = 1.0e-8;
894 pointDistanceTolerance_ = 1.0e-8;
895 TEdge* ea = makeEdge(a, b, true);
896 TEdge* eb = makeEdge(b, c, true);
897 TEdge* ec = makeEdge(c, a, true);
898
899
900 TQuadEdge::splice(ea->sym(), eb);
901 TQuadEdge::splice(eb->sym(), ec);
902 TQuadEdge::splice(ec->sym(), ea);
903
904 startEdge_ = ec;
905 lastLocateEdge_.reset(startEdge_);
906 lastFloodEdge_.reset(startEdge_);
907 edgeCount_ = 6;
908 stackDepth_ = 0;
909
910 boundingPoints_.push_back(a);
911 boundingPoints_.push_back(b);
912 boundingPoints_.push_back(c);
913 }
914
915 TEMPLATE_DEF
916 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::init4( const TPoint2D& a, const TPoint2D& b, const TPoint2D& c, const TPoint2D& d)
917 {
918 allocateInTemp_ = false;
919 tolerance_ = T(1.0e-8);
920 pointDistanceTolerance_ = T(1.0e-8);
921
922 TEdge* ea = makeEdge(a, b, true);
923 TEdge* eb = makeEdge(b, c, true);
924 TEdge* ec = makeEdge(c, d, true);
925 TEdge* ed = makeEdge(d, a, true);
926
927 TQuadEdge::splice(ea->sym(), eb);
928 TQuadEdge::splice(eb->sym(), ec);
929 TQuadEdge::splice(ec->sym(), ed);
930 TQuadEdge::splice(ed->sym(), ea);
931
932 if (prim::inCircle( c,d,a,b ))
933 connect( ec, eb );
934 else
935 connect( eb, ea );
936
937 startEdge_ = ed;
938 lastLocateEdge_.reset(startEdge_);
939 lastFloodEdge_.reset(startEdge_);
940
941 edgeCount_ = 10;
942 stackDepth_ = 0;
943
944 boundingPoints_.push_back(a);
945 boundingPoints_.push_back(b);
946 boundingPoints_.push_back(c);
947 boundingPoints_.push_back(d);
948 }
949
950 TEMPLATE_DEF
951 PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::PlanarMesh( const prim::Aabb2D<T>& iAabb )
952 {
953 allocateInTemp_ = false;
954 TPoint2D topleft(iAabb.min().x, iAabb.max().y);
955 TPoint2D topright(iAabb.max().x, iAabb.max().y);
956 TPoint2D bottomleft(iAabb.min().x, iAabb.min().y);
957 TPoint2D bottomright(iAabb.max().x, iAabb.min().y);
958 init4(topleft, bottomleft, bottomright, topright);
959 }
960
961 TEMPLATE_DEF
962 PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::PlanarMesh( const TPoint2D& a, const TPoint2D& b, const TPoint2D& c, const TPoint2D& d)
963 {
964 allocateInTemp_ = false;
965 init4(a, b, c, d);
966 }
967
968
969 TEMPLATE_DEF
970 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::deletePoint( typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* iEdge)
971 {
972 if (iEdge->handle().point_)
973 free(iEdge->handle().point_);
974 TEdge* currentEdge = iEdge;
975 do
976 {
977 currentEdge->handle().point_ = NULL;
978 currentEdge = currentEdge->oNext();
979 } while ( currentEdge != iEdge );
980 return true;
981 }
982
983 TEMPLATE_DEF
984 PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::~PlanarMesh( )
985 {
986 forAllVertices( TEdgeCallback(this, &PlanarMesh::deletePoint) );
987 typename TQuadEdgeList::iterator qIt;
988 for (qIt = quadEdgeList_.begin(); qIt != quadEdgeList_.end();++qIt)
989 {
990 (*qIt)->edgeDeconstrain();
991 (*qIt)->faceDeconstrain();
992 free(*qIt);
993 }
994 }
995
996 TEMPLATE_DEF
997 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::setTempQuadEdges( bool iAllocateInTemp )
998 {
999 allocateInTemp_ = iAllocateInTemp;
1000 if (!allocateInTemp_ )
1001 {
1002 std::copy(tempQuadEdgeList_.begin(),tempQuadEdgeList_.end(),std::back_inserter(quadEdgeList_));
1003 tempQuadEdgeList_.clear();
1004 }
1005 }
1006
1007
1008 TEMPLATE_DEF
1009 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::makeEmptyEdge( bool makeConstrained )
1010 {
1011 TQuadEdge* qe = make(TQuadEdge(makeConstrained));
1012 if (allocateInTemp_)
1013 tempQuadEdgeList_.push_back( qe );
1014 else
1015 quadEdgeList_.push_back( qe );
1016 edgeCount_ += 2;
1017
1018 TEdge* e = qe->edges();
1019 /*
1020 for (int i=0;i<4;++i)
1021 {
1022 ProxyHandle* pHandle = new ProxyHandle;
1023 e->handle() = pHandle;
1024 e = e->rot();
1025 }
1026 */
1027 return e;
1028 }
1029
1030 TEMPLATE_DEF
1031 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::connect( TEdge* a, TEdge* b)
1032 {
1033 FaceHandle fa = faceHandle(a);
1034 FaceHandle fb = faceHandle(b);
1035 if ( fa != fb )
1036 {
1037 LASS_THROW("connect of edges would violate face constraint");
1038 }
1039 PointHandle hADest = pointHandle( a->sym() );
1040 PointHandle hB = pointHandle( b );
1041 TEdge* e = makeEdge( dest(a), org(b), false);
1042 TQuadEdge::splice( e, a->lNext() );
1043 TQuadEdge::splice( e->sym(), b );
1044
1045 setFaceHandle( e, fa );
1046 setPointHandle( a->sym(), hADest );
1047 setPointHandle( b, hB );
1048 return e;
1049 }
1050
1051 TEMPLATE_DEF
1052 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::split( TEdge* iEdge, T iWhere)
1053 {
1054 LASS_ENFORCE(iWhere>T(0) && iWhere<T(1));
1055 splitEdge(iEdge,along(iEdge,iWhere));
1056 return iEdge->lNext();
1057
1058 //TPoint2D splitPoint(along(iEdge,iWhere));
1059 //TEdge* a(iEdge);
1060 //TEdge* b(iEdge->lNext());
1061 //FaceHandle fa = faceHandle(a);
1062 //FaceHandle fas = faceHandle(a->sym());
1063 //TEdge* e = makeEdge( splitPoint, dest(iEdge), iEdge->isConstrained());
1064 //*iEdge->sym()->handle().point_ = splitPoint;
1065 //TQuadEdge::splice( e, a->lNext() );
1066 //TQuadEdge::splice( e->sym(), b );
1067
1068 //setFaceHandle( e, fa );
1069 //setFaceHandle( e->sym(), fas );
1070 //return e;
1071 }
1072
1073 TEMPLATE_DEF
1074 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::removeVertex( TEdge* iEdge)
1075 {
1076 // enumerate all the incident edges
1077 int n = vertexOrder(iEdge);
1078 TEdge* currentEdge = iEdge;
1079 for (int i=0;i<n;++i)
1080 {
1081 gcDeleteEdge(currentEdge);
1082 currentEdge = currentEdge->oNext();
1083 }
1084 gc();
1085 return true;
1086 }
1087
1088 TEMPLATE_DEF
1089 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::fixEdge(TEdge *e)
1090 // A simple, but recursive way of doing things.
1091 // Flips e, if necessary, in which case calls itself
1092 // recursively on the two edges that were to the right
1093 // of e before it was flipped, to propagate the change.
1094 {
1095 if (e->isConstrained())
1096 return;
1097 TEdge *f = e->oPrev();
1098 TEdge *g = e->dNext();
1099
1100 if ( prim::inCircle( dest(e),dest(e->oNext()),org(e),dest(f)) )
1101 {
1102 swap(e);
1103 fixEdge(f);
1104 fixEdge(g);
1105 }
1106 }
1107
1108
1109 TEMPLATE_DEF
1110 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::safeSplitEdge(TEdge *e, const TPoint2D& iPoint )
1111 // checks that when the point iPoint is on the edge e that is
1112 // really between the org and the dest of e
1113 // due to numerical roundoff this could be otherwise
1114 //
1115 // only call this routine when you are splitting an edge e in two
1116 // collinear (or at least expected collinear within numerical roundoff) parts
1117 {
1118 TRay2D testRay1(org(e),dest(e));
1119 if (testRay1.t(dest(e))<=testRay1.t(iPoint))
1120 {
1121 LASS_THROW("Inconsistent splitting of constrained edge");
1122 }
1123 TRay2D testRay2(dest(e),org(e));
1124 if (testRay2.t(org(e))<=testRay2.t(iPoint))
1125 {
1126 LASS_THROW("Inconsistent splitting of constrained edge");
1127 }
1128 splitEdge(e,iPoint);
1129 }
1130
1131
1132 TEMPLATE_DEF
1133 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::splitEdge(TEdge *e, const TPoint2D& iPoint )
1134 // Shorten edge e s.t. its destination becomes x. Connect
1135 // x to the previous destination of e by a new edge. If e
1136 // is constrained, x is snapped onto the edge, and the new
1137 // edge is also marked as constrained.
1138
1139 // remember that this function is rather dumb
1140 // no snapping is performed to constrained edges
1141 // all these operations have to be performed _before_ entering this routine
1142 // the edge e is readily available so it is the responsible of the caller to
1143 // make meaningful splits
1144 {
1145 const TVector2D dir = direction(e);
1146 PointHandle iE = pointHandle(e);
1147 [[maybe_unused]] PointHandle iES = pointHandle(e->sym());
1148 EdgeHandle iLE = edgeHandle(e);
1149 EdgeHandle iRE = edgeHandle(e->sym());
1150
1151 TPoint2D thePoint = iPoint;
1152 TEdge* t = e->lNext();
1153 TQuadEdge::splice( e->sym(), t );
1154 setDest( thePoint, e );
1155 setPointHandle(e, iE);
1156 TEdge* ne = connect( e, t );
1157 if (e->isConstrained())
1158 {
1159 //setEdgeHandle(ne, iLE);
1160 //setEdgeHandle(ne->sym(), iRE);
1161 setOrientedEdgeHandle( e, iLE, iRE, dir );
1162 setOrientedEdgeHandle( ne, iLE, iRE, dir );
1163 ne->quadEdge()->edgeConstrain();
1164 }
1165
1166 // make sure the inserted edge also points to correct point information
1167 LASS_ASSERT( pointHandle(e) == iE );
1168 LASS_ASSERT( pointHandle(e->lNext()->sym()) == iES );
1169 }
1170
1171 TEMPLATE_DEF
1172 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TPoint2D PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::snap(const TPoint2D& x, const TPoint2D& a, const TPoint2D& b)
1173 {
1174//#pragma LASS_TODO("This can be more efficient")
1175 /*
1176 if (distance(x,a)<tolerance_)
1177 return a;
1178 if (distance(x,b)<tolerance_)
1179 return b;
1180 */
1181 if (x==a)
1182 return a;
1183 if (x==b)
1184 return b;
1185 T t1 = dot(x-a,b-a);
1186 T t2 = dot(x-b,a-b);
1187
1188 T t = std::max(t1,t2) / (t1 + t2);
1189
1190 if (prim::doubleTriangleArea(x,a,b)==T(0))
1191 return x;
1192
1193 if (t1>t2)
1194 {
1195 T rx1 = (T(1)-t)*a.x + t*b.x;
1196 T ry1 = (T(1)-t)*a.y + t*b.y;
1197 return TPoint2D(rx1,ry1);
1198 }
1199
1200 T rx2 = (T(1)-t)*b.x + t*a.x;
1201 T ry2 = (T(1)-t)*b.y + t*a.y;
1202
1203 return TPoint2D(rx2,ry2);
1204 }
1205
1206 TEMPLATE_DEF
1207 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::swap( TEdge* iEdge )
1208 {
1209 if ( iEdge->isConstrained() )
1210 {
1211 LASS_THROW( "cannot swap constrained edge" );
1212 }
1213 TEdge* e = iEdge;
1214 TEdge* a = e->oPrev();
1215 TEdge* b = e->sym()->oPrev();
1216
1217 /* following handles do not represent the edges, they represent the points
1218 a
1219 / \
1220 b-e-c
1221 \ /
1222 d
1223
1224 */
1225 PointHandle hA = pointHandle(e->oNext()->sym());
1226 [[maybe_unused]] PointHandle hB = pointHandle(e);
1227 [[maybe_unused]] PointHandle hC = pointHandle(e->sym());
1228 PointHandle hD = pointHandle(e->oPrev()->sym());
1229
1230 TQuadEdge::splice( e, a );
1231 TQuadEdge::splice( e->sym(), b );
1232 TQuadEdge::splice( e, a->lNext() );
1233 TQuadEdge::splice( e->sym(), b->lNext() );
1234 setOrg( dest(a), e ); // this is point a here!
1235 setDest( dest(b), e ); // this is point b here!
1236
1237 setPointHandle( e, hD );
1238 setPointHandle( e->sym(), hA );
1239
1240 LASS_ASSERT( leftOf( dest(e->lNext()), e ) );
1241 LASS_ASSERT( leftOf( dest(e->lNext()->lNext()), e->lNext() ) );
1242 LASS_ASSERT( leftOf( dest(e->lNext()->lNext()->lNext()), e->lNext()->lNext() ) );
1243 LASS_ASSERT( leftOf( dest(e->sym()->lNext()), e->sym() ) );
1244 LASS_ASSERT( leftOf( dest(e->sym()->lNext()->lNext()), e->sym()->lNext() ) );
1245 LASS_ASSERT( leftOf( dest(e->sym()->lNext()->lNext()->lNext()), e->sym()->lNext()->lNext() ) );
1246
1247 LASS_ASSERT( pointHandle( e->oNext()->sym() ) == hB );
1248 LASS_ASSERT( pointHandle( e->oPrev()->sym() ) == hC );
1249 }
1250
1251 TEMPLATE_DEF
1252 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::deleteEdge( TEdge* e )
1253 {
1254 if ( e->isConstrained() )
1255 {
1256 return false;
1257 }
1258
1259 while (startEdge_->quadEdge()==e->quadEdge())
1260 startEdge_=e->lNext();
1261
1262 // void any caching when the edge is deleted
1263 if ( *lastLocateEdge_==e
1264 || *lastLocateEdge_==e->rot()
1265 || *lastLocateEdge_==e->sym()
1266 || *lastLocateEdge_==e->sym()->rot())
1267 lastLocateEdge_.reset(startEdge_);
1268 if ( *lastFloodEdge_==e
1269 || *lastFloodEdge_==e->rot()
1270 || *lastFloodEdge_==e->sym()
1271 || *lastFloodEdge_==e->sym()->rot())
1272 lastFloodEdge_.reset(startEdge_);
1273
1274 TQuadEdge::splice( e, e->oPrev() );
1275 TQuadEdge::splice( e->sym(), e->sym()->oPrev() );
1276
1277 typename TQuadEdgeList::iterator it = std::find(quadEdgeList_.begin(),quadEdgeList_.end(),e->quadEdge());
1278 std::swap(*it, quadEdgeList_.back());
1279
1280 TQuadEdge* qe = quadEdgeList_.back();
1281 qe->edgeDeconstrain();
1282 qe->faceDeconstrain();
1283 quadEdgeList_.pop_back();
1284 free(qe);
1285
1286 edgeCount_ -= 2;
1287 return true;
1288 }
1289
1290 TEMPLATE_DEF
1291 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::gcDeleteEdge( TEdge* e )
1292 {
1293 if ( e->isConstrained() )
1294 {
1295 return false;
1296 }
1297
1298 while (startEdge_->quadEdge()==e->quadEdge())
1299 startEdge_=e->lNext();
1300 // void any caching when the edge is deleted
1301 if ( *lastLocateEdge_==e
1302 || *lastLocateEdge_==e->rot()
1303 || *lastLocateEdge_==e->sym()
1304 || *lastLocateEdge_==e->sym()->rot())
1305 lastLocateEdge_.reset(startEdge_);
1306 if ( *lastFloodEdge_==e
1307 || *lastFloodEdge_==e->rot()
1308 || *lastFloodEdge_==e->sym()
1309 || *lastFloodEdge_==e->sym()->rot())
1310 lastFloodEdge_.reset(startEdge_);
1311
1312 TQuadEdge::splice( e, e->oPrev() );
1313 TQuadEdge::splice( e->sym(), e->sym()->oPrev() );
1314
1315 TQuadEdge* qe = e->quadEdge();
1316 qe->edgeDeconstrain();
1317 qe->faceDeconstrain();
1318 qe->detach();
1319 qe->deleted = true;
1320
1321 edgeCount_ -= 2;
1322 return true;
1323 }
1324 TEMPLATE_DEF
1325 int PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::gc( )
1326 {
1327 // sort all the edges in the quadedgelist on the deleted value
1328 // because there are only two values we can do an O(n) sort, yes you read that correctly =)
1329 // we put all non-gc elements in front of the vector
1330 // this algorithm will also work on lists using iterators, we just do it on random access containers here
1331 // the iterator implementation is left as an exercise to the reader =)
1332 int lastMarker = quadEdgeList_.size()-1;
1333 int startMarker = 0;
1334 int collected = 0;
1335 int n = quadEdgeList_.size()-1;
1336
1337 while (startMarker<lastMarker)
1338 {
1339 // advance first marker till first gc element
1340 while (!quadEdgeList_[startMarker]->deleted && startMarker<n)
1341 ++startMarker;
1342 // advance last marker till first non-gc element
1343 while (quadEdgeList_[lastMarker]->deleted && lastMarker>0)
1344 --lastMarker;
1345 // double check and do termination iteration
1346 if (quadEdgeList_[startMarker]->deleted && !quadEdgeList_[lastMarker]->deleted )
1347 std::swap(quadEdgeList_[startMarker],quadEdgeList_[lastMarker]);
1348 }
1349 // lastMarker should now point at the GC starting point
1350 // we are too lazy to verify this, so we do backwards search again
1351 // the whole algorithm is O(n) anyway
1352 lastMarker = n;
1353 while (quadEdgeList_[lastMarker]->deleted && lastMarker>0)
1354 {
1355 --lastMarker;
1356 quadEdgeList_.pop_back();
1357 ++collected;
1358 }
1359 return collected;
1360 }
1361
1362
1363 TEMPLATE_DEF
1364 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::bruteForceLocate( const TPoint2D& iPoint ) const
1365 {
1366 typedef PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle> TPlanarMesh;
1367 impl::BrutePointLocator<T, PointHandle, EdgeHandle, FaceHandle> bruteLocator1( iPoint );
1368 const_cast<TPlanarMesh*>(this)->forAllFaces( TEdgeCallback( &bruteLocator1, &impl::BrutePointLocator<T, PointHandle, EdgeHandle, FaceHandle>::findEdge ) );
1369 if (bruteLocator1.edge==NULL)
1370 {
1371 impl::BrutePointLocator<T, PointHandle, EdgeHandle, FaceHandle> bruteLocator2( iPoint );
1372 const_cast<TPlanarMesh*>(this)->forAllPrimaryEdges( TEdgeCallback( &bruteLocator2, &impl::BrutePointLocator<T, PointHandle, EdgeHandle, FaceHandle>::findEdge ) );
1373#if DEBUG_MESH
1374 if (bruteLocator2.edge==NULL)
1375 {
1376 impl::BrutePointLocatorVerbose<T, PointHandle, EdgeHandle, FaceHandle> bruteLocatorVerbose( iPoint );
1377 bruteLocatorVerbose.stream.open("bruteforcelocation.txt");
1378 const_cast<TPlanarMesh*>(this)->forAllPrimaryEdges( TEdgeCallback( &bruteLocatorVerbose, &impl::BrutePointLocatorVerbose<T, PointHandle, EdgeHandle, FaceHandle>::findEdge ) );
1379 bruteLocatorVerbose.stream.close();
1380 return bruteLocatorVerbose.edge;
1381 }
1382#endif
1383 return bruteLocator2.edge;
1384 }
1385 return bruteLocator1.edge;
1386 }
1387
1388 // searches for an exact vertex
1389 TEMPLATE_DEF
1390 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::bruteForceExactLocate( const TPoint2D& iPoint ) const
1391 {
1392 typedef PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle> TPlanarMesh;
1393 impl::BrutePointExactLocator<T, PointHandle, EdgeHandle, FaceHandle> brutePointLocator( iPoint );
1394 const_cast<TPlanarMesh*>(this)->forAllVertices( TEdgeCallback( &brutePointLocator, &impl::BrutePointExactLocator<T, PointHandle, EdgeHandle, FaceHandle>::findPoint ) );
1395 return brutePointLocator.edge;
1396 }
1397
1398 TEMPLATE_DEF
1399 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::isBoundingPoint( const TPoint2D& iPoint ) const
1400 {
1401 size_t length = boundingPoints_.size();
1402 for (size_t i=0;i<length;++i)
1403 {
1404 if (boundingPoints_[i]==iPoint)
1405 return true;
1406 }
1407 return false;
1408 }
1409
1410
1411
1412 TEMPLATE_DEF
1413 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::locate( const TPoint2D& iPoint, TEdge* iHint ) const
1414 {
1415 if (iHint)
1416 *lastLocateEdge_ = iHint;
1417 TEdge*& lastLocateEdge = *lastLocateEdge_;
1418 long edgesPassed = 0;
1419
1420 TEdge* e = lastLocateEdge;
1421 while (edgesPassed<(edgeCount_+2))
1422 {
1423continueSearch:
1424 ++edgesPassed;
1425 const TPoint2D& eorg = fastOrg(e);
1426 const TPoint2D& edest = fastDest(e);
1427
1428 if ( eorg == iPoint )
1429 {
1430 lastLocateEdge = e;
1431 return e;
1432 }
1433 if ( edest == iPoint )
1434 {
1435 e = e->sym();
1436 lastLocateEdge = e;
1437 return e;
1438 }
1439 if ( fastRightOf( iPoint, e ) )
1440 e = e->sym();
1441
1442 //TEdge* elPrev= e->lPrev();
1443 bool hasALeftFace = true;
1444 if (isBoundingPoint(eorg))
1445 hasALeftFace = fastHasLeftFace(e);
1446 //bool hasALeftFace = fastLeftOf(edest,elPrev);
1447 if (hasALeftFace)
1448 {
1449 /*
1450 if ( iPoint == fastOrg(elPrev))
1451 {
1452 return e->lPrev();
1453 }
1454 */
1455 if ( fastLeftOf( iPoint, e->oNext()) )
1456 {
1457 e = e->oNext();
1458 continue;
1459 }
1460 TEdge* ce = e;
1461 // [TODO] Optimize
1462 // this for loop is introduced for point location in non-triangular, general
1463 // convex cells
1464 int loopOrder = chainOrder(e)-2;
1465 for (int i=0;i<loopOrder;++i)
1466 {
1467 if ( fastLeftOf( iPoint, ce->dPrev()) )
1468 {
1469 e = ce->dPrev();
1470 goto continueSearch;
1471 }
1472 ce = ce->lNext();
1473 }
1474
1475 TRay2D R1(fastOrg(e), fastDest(e));
1476 TPoint2D p1 = R1.point( R1.t( iPoint ) );
1477 TRay2D R2(fastDest(e), fastOrg(e->lPrev()));
1478 TPoint2D p2 = R2.point( R2.t( iPoint ) );
1479 TRay2D R3(fastOrg(e->lPrev()), fastOrg(e));
1480 TPoint2D p3 = R3.point( R3.t( iPoint ) );
1481
1482 typename TPoint2D::TValue d1 = squaredDistance(iPoint,p1);
1483 typename TPoint2D::TValue d2 = squaredDistance(iPoint,p2);
1484 typename TPoint2D::TValue d3 = squaredDistance(iPoint,p3);
1485
1486 /*if (d1*d2*d3 == T(0))
1487 {
1488 // yes, we are unlucky
1489 int a = 0;
1490 }*/
1491
1492 if ((d1<d2) && (d1<d3))
1493 {
1494 lastLocateEdge = e;
1495 return e;
1496 }
1497 if ((d2<d3) && (d2<=d1))
1498 {
1499 e = e->lNext();
1500 lastLocateEdge = e;
1501 return e;
1502 }
1503 e = e->lPrev();
1504 lastLocateEdge = e;
1505 return e;
1506 }
1507
1508 if (onEdge(iPoint,e))
1509 {
1510 e = e->sym();
1511 lastLocateEdge = e;
1512 return e;
1513 }
1514
1515 if (rightOf(iPoint, e->oPrev()))
1516 {
1517 e = e->oPrev()->sym();
1518 continue;
1519 }
1520 if (rightOf(iPoint,e->dNext()))
1521 {
1522 e = e->dNext()->sym();
1523 continue;
1524 }
1525 }
1526 //throw std::runtime_error("could not locate edge for point, probably point on edge which is not supported yet!");
1527
1528 e = bruteForceLocate(iPoint);
1529 lastLocateEdge = e;
1530 return e;
1531 }
1532
1533 TEMPLATE_DEF
1534 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::pointLocate( const TPoint2D& iPoint ) const
1535 {
1536 TEdge*& lastLocateEdge = *lastLocateEdge_;
1537 long edgesPassed = 0;
1538
1539 TEdge* e = lastLocateEdge;
1540 while (edgesPassed<(edgeCount_+2))
1541 {
1542 ++edgesPassed;
1543 if ( org( e ) == iPoint )
1544 {
1545 lastLocateEdge = e;
1546 return e;
1547 }
1548 if ( dest( e ) == iPoint )
1549 {
1550 e = e->sym();
1551 lastLocateEdge = e;
1552 return e;
1553 }
1554 if ( rightOf( iPoint, e ) )
1555 e = e->sym();
1556
1557 if (hasLeftFace(e))
1558 {
1559 if ( iPoint == org(e->lPrev()) )
1560 {
1561 e = e->lPrev();
1562 lastLocateEdge = e;
1563 return e;
1564 }
1565 if ( leftOf( iPoint, e->oNext()) )
1566 {
1567 e = e->oNext();
1568 continue;
1569 }
1570 if ( leftOf( iPoint, e->dPrev()) )
1571 {
1572 e = e->dPrev();
1573 continue;
1574 }
1575 }
1576
1577 if (rightOf(iPoint, e->oPrev()))
1578 {
1579 e = e->oPrev()->sym();
1580 continue;
1581 }
1582 if (rightOf(iPoint,e->dNext()))
1583 {
1584 e = e->dNext()->sym();
1585 continue;
1586 }
1587 }
1588 //throw std::runtime_error("could not locate edge for point, probably point on edge which is not supported yet!");
1589 e = bruteForceExactLocate(iPoint);
1590 LASS_ASSERT(e!=NULL);
1591 lastLocateEdge = e;
1592 return e;
1593 // LASS_THROW( "pointLocate: could not find point");
1594 }
1595
1596
1597 /** method is polygon-safe */
1598 TEMPLATE_DEF
1599 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::shoot( const TRay2D& iRay ) const
1600 {
1601 TEdge* locateEdge = locate(iRay.support());
1602 TEdge* startEdge = locateEdge;
1603 if (!startEdge)
1604 return NULL;
1605
1606 if (org(startEdge)==iRay.support())
1607 {
1608 TEdge* e(startEdge);
1609 //lass::prim::Side lastSide = iRay.classify(org(e));
1610 int vOrder = vertexOrder(e);
1611 for (int i=0;i<vOrder+1;++i)
1612 {
1613 /*
1614 lass::prim::Side currentSide = iRay.classify(dest(e));
1615 if (lastSide==lass::prim::sRight && lastSide!=currentSide)
1616 return e;
1617 lastSide = currentSide;
1618 */
1619 if ( num::abs(prim::doubleTriangleArea(iRay.support(),iRay.point(T(1)),dest(e)))<tolerance_
1620 && iRay.t(dest(e)) > T(0) )
1621 return e;
1622 e = e->oNext();
1623 }
1624 }
1625
1626 do
1627 {
1628 TEdge* e(startEdge);
1629 lass::prim::Side lastSide = iRay.classify(org(e));
1630 do
1631 {
1632 lass::prim::Side currentSide = iRay.classify(dest(e));
1633 if (currentSide==lass::prim::sSurface)
1634 return e->sym();
1635 if (lastSide==lass::prim::sRight && lastSide!=currentSide)
1636 return e;
1637 lastSide = currentSide;
1638 e = e->lNext();
1639 }
1640 while (e!=startEdge);
1641 startEdge = startEdge->oNext();
1642 } while (startEdge!=locateEdge);
1643
1644#if DEBUG_MESH
1645 // recreate the same situation for debugging purposes
1646 {
1647 lass::io::MatlabOStream bugMesh;
1648 bugMesh.open( "shoot_ray.m" );
1649 bugMesh << std::setprecision(15);
1650 bugMesh << *const_cast<PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>*>(this);
1651 bugMesh.close();
1652 bugMesh.open( "shoot_ray_edges.m" );
1653 bugMesh.setColor(lass::io::mcBlue);
1654 bugMesh << TLineSegment2D(iRay.support(),iRay.point(1.0));
1655 do
1656 {
1657 TEdge* e(startEdge);
1658 lass::prim::Side lastSide = iRay.classify(org(e));
1659 bugMesh.setColor(lass::io::mcRed);
1660 bugMesh << TLineSegment2D(org(startEdge),dest(startEdge));
1661 do
1662 {
1663 bugMesh.setColor(lass::io::mcGreen);
1664 bugMesh << TLineSegment2D(org(e),dest(e));
1665 lass::prim::Side currentSide = iRay.classify(dest(e));
1666 if (lastSide==lass::prim::sRight && lastSide!=currentSide)
1667 return e;
1668 lastSide = currentSide;
1669 e = e->lNext();
1670 }
1671 while (e!=startEdge);
1672 startEdge = startEdge->oNext();
1673 } while (startEdge!=locateEdge);
1674 bugMesh.close();
1675 }
1676#endif
1677
1678 return NULL; // we give up
1679 }
1680
1681 /** method is polygon-safe */
1682 TEMPLATE_DEF
1683 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::pointShoot( const TRay2D& iRay ) const
1684 {
1685 TEdge* locateEdge = pointLocate(iRay.support());
1686 LASS_ASSERT(isBoundingPoint(org(locateEdge)) || allEqualChainOrder(locateEdge));
1687 TEdge* startEdge = locateEdge;
1688 if (!startEdge)
1689 return NULL;
1690
1691 {
1692 TEdge* e(startEdge);
1693 //lass::prim::Side lastSide = iRay.classify(org(e));
1694 int vOrder = vertexOrder(e);
1695 for (int i=0;i<vOrder+1;++i)
1696 {
1697 /*
1698 lass::prim::Side currentSide = iRay.classify(dest(e));
1699 if (lastSide==lass::prim::sRight && lastSide!=currentSide)
1700 return e;
1701 lastSide = currentSide;
1702 */
1703 if ( num::abs(prim::doubleTriangleArea(iRay.support(),iRay.point(T(1)),dest(e)))<tolerance_
1704 && iRay.t(dest(e)) > T(0) )
1705 return e;
1706 e = e->oNext();
1707 }
1708 }
1709
1710 do
1711 {
1712 TEdge* e(startEdge);
1713 lass::prim::Side lastSide = iRay.classify(org(e));
1714 do
1715 {
1716 lass::prim::Side currentSide = iRay.classify(dest(e));
1717 if (lastSide==lass::prim::sRight && lastSide!=currentSide)
1718 return e;
1719 lastSide = currentSide;
1720 e = e->lNext();
1721 }
1722 while (e!=startEdge);
1723 startEdge = startEdge->oNext();
1724 } while (startEdge!=locateEdge);
1725
1726 LASS_THROW("pointShoot, could not find suitable triangle");
1727 // return NULL;
1728 }
1729
1730
1731 /* _Adds_ the intersection points of the crossed edges to the output vector, method is polygon safe */
1732 TEMPLATE_DEF
1733 template <typename OutputIterator>
1734 OutputIterator PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::walkIntersections( const TLineSegment2D& iSegment, OutputIterator oIntersections) const
1735 {
1736 std::vector<TEdge*> crossedEdges;
1737 walk(iSegment,std::back_inserter(crossedEdges));
1738 const size_t n = crossedEdges.size();
1739 for (size_t i=0;i<n;++i)
1740 {
1741 TPoint2D intersection;
1742 if (impl::fastIntersect(iSegment.tail(),iSegment.head(),fastOrg(crossedEdges[i]),fastDest(crossedEdges[i]),intersection)!=lass::prim::rOne)
1743 {
1744 // we make the bold assumption that we have parallel coinciding lines
1745 (*oIntersections++) = std::pair<TPoint2D,TEdge*>(fastOrg(crossedEdges[i]),crossedEdges[i]);
1746 if (squaredDistance(fastOrg(crossedEdges[i]),fastDest(crossedEdges[i]))<squaredDistance(iSegment.tail(),iSegment.head()))
1747 (*oIntersections++) = std::pair<TPoint2D,TEdge*>(fastDest(crossedEdges[i]),crossedEdges[i]);
1748 else
1749 (*oIntersections++) = std::pair<TPoint2D,TEdge*>(iSegment.head(),crossedEdges[i]);
1750 }
1751 else
1752 {
1753 (*oIntersections++) = std::pair<TPoint2D,TEdge*>(intersection,crossedEdges[i]);
1754 }
1755 }
1756 return oIntersections;
1757 }
1758
1759
1760 /* _Adds_ the crossed edges to the output vector, method is polygon safe */
1761 TEMPLATE_DEF
1762 template <typename OutputIterator>
1763 OutputIterator PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::walk( const TLineSegment2D& iSegment, OutputIterator crossedEdges ) const
1764 {
1765 TEdge* e=NULL;
1766 e = shoot(TRay2D(iSegment.tail(),iSegment.vector()));
1767 //LASS_ASSERT(isBoundingPoint(org(e)) || allEqualChainOrder(e));
1768 if (!e)
1769 {
1770 LASS_THROW("Could not shoot initial ray in walk");
1771 }
1772 if (inConvexCell(iSegment.head(),e))
1773 return crossedEdges;
1774 if (dest(e)==iSegment.head())
1775 {
1776 (*crossedEdges++) = e;
1777 return crossedEdges;
1778 }
1779
1780 while ( !leftOf(iSegment.head(),e)
1781 || ( (num::abs(prim::doubleTriangleArea(dest(e),iSegment.head(),iSegment.tail()))<tolerance_ )
1782 && (num::abs(prim::doubleTriangleArea(org(e),iSegment.head(),iSegment.tail()))<tolerance_ ) ) )
1783 {
1784 (*crossedEdges++) = e;
1785 /*
1786 TEdge* ne1 = e->sym()->lNext();
1787 if (cw(iSegment.tail(),iSegment.head(),dest(ne1)))
1788 e = ne1->lNext();
1789 else
1790 e = ne1;
1791 */
1792 TEdge* ce = e->sym()->lNext();
1793#pragma LASS_TODO("Optimize")
1794 // this for loop is introduced for point location in non-triangular, general
1795 // convex cells
1796 TEdge* oldE = e;
1797 for (int i=0;i<chainOrder(e);++i)
1798 {
1799 if ( weakCcw( iSegment.tail(), iSegment.head(), dest(ce) ) )
1800 {
1801 e = ce;
1802 break;
1803 }
1804 ce = ce->lNext();
1805 }
1806 if (e==oldE)
1807 {
1808 LASS_THROW("Planarmesh: stuck in walk");
1809 }
1810 }
1811 return crossedEdges;
1812 }
1813
1814 /* _Adds_ the crossed edges to the output vector, method is polygon safe */
1815 TEMPLATE_DEF
1816 std::pair<int, typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge*> PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::walkTillConstrained( const TRay2D& iRay) const
1817 {
1818 int nCrossed = 0;
1819 TEdge* e=NULL;
1820 TLineSegment2D segment(iRay.support(),iRay.support()+iRay.direction());
1821 e = shoot(iRay);
1822 //LASS_ASSERT(isBoundingPoint(org(e)) || allEqualChainOrder(e));
1823 if (!e)
1824 {
1825 LASS_THROW("Could not shoot initial ray in walk");
1826 }
1827 while ( !e->isConstrained())
1828 {
1829 /*
1830 TEdge* ne1 = e->sym()->lNext();
1831 if (cw(iSegment.tail(),iSegment.head(),dest(ne1)))
1832 e = ne1->lNext();
1833 else
1834 e = ne1;
1835 */
1836 TEdge* ce = e->sym()->lNext();
1837#pragma LASS_TODO("Optimize")
1838 // this for loop is introduced for point location in non-triangular, general
1839 // convex cells
1840 TEdge* oldE = e;
1841 int n = chainOrder(e);
1842 for (int i=0;i<n;++i)
1843 {
1844 if ( weakCcw( segment.tail(), segment.head(), dest(ce) ) )
1845 {
1846 e = ce;
1847 ++nCrossed;
1848 break;
1849 }
1850 ce = ce->lNext();
1851 }
1852 if (e==oldE)
1853 {
1854 LASS_THROW("Planarmesh: stuck in walkTillConstrained");
1855 }
1856 }
1857 return std::make_pair(nCrossed,e);
1858 }
1859
1860
1861
1862 /* _Adds_ the crossed edges to the output vector, method is polygon safe */
1863 TEMPLATE_DEF
1864 template <typename OutputIterator>
1865 OutputIterator PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::pointWalk( const TLineSegment2D& iSegment, OutputIterator crossedEdges ) const
1866 {
1867 TEdge* e=NULL;
1868 e = pointShoot(TRay2D(iSegment.tail(),iSegment.vector()));
1869 LASS_ASSERT(isBoundingPoint(org(e)) || allEqualChainOrder(e));
1870 if (!e)
1871 {
1872 LASS_THROW("Could not shoot initial ray in pointWalk");
1873 }
1874 if (dest(e)==iSegment.head())
1875 {
1876 (*crossedEdges++) = e;
1877 return crossedEdges;
1878 }
1879
1880 while ( !leftOf(iSegment.head(),e)
1881 || ( (num::abs(prim::doubleTriangleArea(dest(e),iSegment.head(),iSegment.tail()))<tolerance_ )
1882 && (num::abs(prim::doubleTriangleArea(org(e),iSegment.head(),iSegment.tail()))<tolerance_ ) ) )
1883 {
1884 (*crossedEdges++) = e;
1885 /*
1886 TEdge* ne1 = e->sym()->lNext();
1887 if (cw(iSegment.tail(),iSegment.head(),dest(ne1)))
1888 e = ne1->lNext();
1889 else
1890 e = ne1;
1891 */
1892 TEdge* ce = e->sym()->lNext();
1893#pragma LASS_TODO("Optimize")
1894 // this for loop is introduced for point location in non-triangular, general
1895 // convex cells
1896 for (int i=0;i<chainOrder(e)-1;++i)
1897 {
1898 if ( weakCcw( iSegment.tail(), iSegment.head(), dest(ce) ) )
1899 {
1900 e = ce;
1901 break;
1902 }
1903 ce = ce->lNext();
1904 }
1905 }
1906 return crossedEdges;
1907 }
1908
1909
1910 TEMPLATE_DEF
1911 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::triangulate( TEdge* iEdge )
1912 // Triangulates the left face of first, which is assumed to be closed.
1913 // It is also assumed that all the vertices of that face lie to the
1914 // left of last (the edge preceeding first).
1915 // This is NOT intended as a general simple polygon triangulation
1916 // routine. It is called by InsertEdge in order to restore a
1917 // triangulation after an edge has been inserted.
1918 // The routine consists of two loops: the outer loop continues as
1919 // long as the left face is not a triangle. The inner loop examines
1920 // successive triplets of vertices, creating a triangle whenever the
1921 // triplet is counterclockwise.
1922 {
1923 TEdge *first = iEdge;
1924 TEdge *a, *b, *e, *t1, *t2, *last = first->lPrev();
1925
1926 while (first->lNext()->lNext() != last)
1927 {
1928 e = first->lNext();
1929 t1 = first;
1930 while (e != last)
1931 {
1932 t2 = e->lNext();
1933 if (t2 == last && t1 == first)
1934 break;
1935 if (leftOf(dest(e),t1))
1936 {
1937 if (t1 == first)
1938 t1 = first = connect(e, t1)->sym();
1939 else
1940 t1 = connect(e, t1)->sym();
1941 a = t1->oPrev(), b = t1->dNext();
1942 fixEdge(a);
1943 fixEdge(b);
1944 e = t2;
1945 }
1946 else
1947 {
1948 t1 = e;
1949 e = t2;
1950 }
1951 }
1952 }
1953 a = last->lNext(), b = last->lPrev();
1954 fixEdge(a);
1955 fixEdge(b);
1956 fixEdge(last);
1957 }
1958
1959
1960 /** makeDelaunay: after insertion of the site force a _Delaunay_ triangulation
1961 * forceOnEdge : points are considered on the edge although due to numerical roundoff they
1962 * maybe not, as user you should let this default to false or you'd better
1963 * be knowning damn well what you are doing :-D
1964 */
1965 TEMPLATE_DEF
1966 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::insertSite( const TPoint2D& iPoint, bool makeDelaunay, bool forceOnEdge )
1967 {
1968 TEdge* e = locate( iPoint );
1969 if (e==NULL)
1970 {
1971 locate( iPoint );
1972 LASS_ASSERT( e!=NULL);
1973 }
1974 T sqDistOrgEP = prim::squaredDistance(org(e),iPoint);
1975 T sqDistDestEP = prim::squaredDistance(dest(e),iPoint);
1976
1977 if (sqDistOrgEP<pointDistanceTolerance_*pointDistanceTolerance_)
1978 return e;
1979 if (sqDistDestEP<pointDistanceTolerance_*pointDistanceTolerance_)
1980 return e->sym();
1981
1982 bool hasLeft = hasLeftFace(e);
1983 bool hasRight = hasRightFace(e);
1984
1985 if (!hasLeft && !hasRight)
1986 {
1987#if DEBUG_MESH
1988 lass::io::MatlabOStream bugMesh;
1989 bugMesh.open( "bugMesh.m" );
1990 bugMesh << *this;
1991 bugMesh.close();
1992#endif
1993 throw std::runtime_error("insertSite: edge does not have any faces");
1994 }
1995
1996 bool isOnEdge = onEdge(iPoint, e) || forceOnEdge;
1997 // the distance to lines is enlarged a bit for greater stability towards distances to line calculations
1998 if (!isOnEdge && TLine2D(org(e),dest(e)).squaredDistance(iPoint)<T(1.05)*pointDistanceTolerance_*pointDistanceTolerance_)
1999 {
2000 // snap the point to the line, this is a precaution to avoid numerical
2001 // instability later on
2002 isOnEdge = true;
2003 TPoint2D x = snap(iPoint, org(e), dest(e));
2004 return insertSite(x,makeDelaunay,true);
2005 }
2006
2007 bool insideLeft = hasLeft && (isOnEdge || leftOf(iPoint, e)) &&
2008 rightOf(iPoint, e->oNext()) && rightOf(iPoint, e->dPrev());
2009 bool insideRight = hasRight && (isOnEdge || rightOf(iPoint, e)) &&
2010 leftOf(iPoint, e->oPrev()) && leftOf(iPoint, e->dNext());
2011
2012 // check for numerical stability
2013 // we check if upon insertion the orientation tests still give meaningful results
2014 {
2015 TPoint2D p1 = org(e);
2016 TPoint2D p2 = dest(e);
2017 TPoint2D p3 = dest(e->lNext());
2018
2019 T s1 = prim::doubleTriangleArea(iPoint,p1,p2);
2020 T s2 = prim::doubleTriangleArea(iPoint,p2,p3);
2021 T s3 = prim::doubleTriangleArea(iPoint,p3,p1);
2022
2023 if (isOnEdge)
2024 {
2025 if ( s1*s2 < T(0) || s2*s3 < T(0) || s1*s3 < T(0) )
2026 {
2027 if (squaredDistance(org(e),iPoint)<squaredDistance(dest(e),iPoint))
2028 {
2029 if (squaredDistance(org(e),iPoint)<pointDistanceTolerance_*pointDistanceTolerance_)
2030 return insertSite(org(e),makeDelaunay,true);
2031 else
2032 {
2033 // we have come across a situation which is combinatorially not possible except for
2034 // numerical issues, we just continue as the algorithm can cope with this
2035 std::cout << "Numerical instability detected in lass::mesh\n";
2036 }
2037 }
2038 else
2039 {
2040 if (squaredDistance(dest(e),iPoint)<pointDistanceTolerance_*pointDistanceTolerance_)
2041 return insertSite(dest(e),makeDelaunay,true);
2042 else
2043 {
2044 // we have come across a situation which is combinatorially not possible except for
2045 // numerical issues, we just continue as the algorithm can cope with this
2046 std::cout << "Numerical instability detected in lass::mesh\n";
2047 }
2048 }
2049 }
2050 }
2051 else
2052 {
2053 LASS_ASSERT( s1*s2 > T(0) );
2054 LASS_ASSERT( s2*s3 > T(0) );
2055 LASS_ASSERT( s1*s3 > T(0) );
2056 }
2057 }
2058
2059
2060 if (insideLeft && iPoint == dest(e->oNext()))
2061 return e->lPrev();
2062 if (insideRight && iPoint == dest(e->oPrev()))
2063 return e->dNext();
2064
2065 if ( isOnEdge )
2066 {
2067 // snap x to e, and check for coincidence:
2068 TPoint2D x = snap(iPoint, org(e), dest(e));
2069 T sqDistOrgEX = prim::squaredDistance(org(e),x);
2070 T sqDistDestEX = prim::squaredDistance(dest(e),x);
2071
2072 if (sqDistOrgEX<pointDistanceTolerance_*pointDistanceTolerance_)
2073 return e;
2074 if (sqDistDestEX<pointDistanceTolerance_*pointDistanceTolerance_)
2075 return e->sym();
2076
2077 // bummer
2078 if (hasRight && hasLeft)
2079 {
2080 // has two faces
2081 TEdge *a, *b, *c, *d;
2082 a = e->oPrev();
2083 b = e->dNext();
2084 c = e->lNext();
2085 d = e->lPrev();
2086 safeSplitEdge(e, x);
2087 connect(e, e->lPrev());
2088 connect(e->oPrev(), e->sym());
2089 fixEdge(a);
2090 fixEdge(b);
2091 fixEdge(c);
2092 fixEdge(d);
2093 }
2094 else
2095 {
2096 // has only one face
2097 if (hasRight)
2098 e = e->sym();
2099 TEdge *c = e->lNext();
2100 TEdge *d = e->lPrev();
2101 safeSplitEdge(e, x);
2102 connect(e, e->lPrev());
2103 fixEdge(c);
2104 fixEdge(d);
2105 }
2106 LASS_ASSERT(vertexOrder(e->sym())>2);
2107 LASS_ASSERT(x==dest(e));
2108 LASS_ASSERT(allEqualChainOrder(e->sym()));
2109 return e->sym();
2110 }
2111
2112 // unused? [Bramz]
2113 // int chOrder = chainOrder(e);
2114
2115 PointHandle pH = pointHandle(e);
2116 TEdge* base = makeEdge(org(e),iPoint,false);
2117 TQuadEdge::splice(base,e);
2118 setPointHandle(e, pH);
2119 TEdge* startBase = base->sym();
2120
2121 do
2122 {
2123 base = connect( e, base->sym() );
2124 e = base->oPrev();
2125 } while (e->dPrev()!=startBase);
2126
2127 if ( makeDelaunay )
2128 {
2129 // try to build a delaunay mesh
2130 while (true)
2131 {
2132 TEdge* t=e->oPrev();
2133 if ( !e->isConstrained() && prim::inCircle(org(e),dest(t),dest(e),iPoint))
2134 {
2135 swap(e);
2136 e=e->oPrev();
2137 }
2138 else
2139 {
2140 if (e->lPrev()==startBase)
2141 return base->sym();
2142 else
2143 e = e->oNext()->lPrev();
2144 }
2145 }
2146 }
2147 LASS_ASSERT(vertexOrder(base->sym())>2);
2148 LASS_ASSERT(iPoint==dest(base));
2149 LASS_ASSERT(allEqualChainOrder(base->sym()));
2150 return base->sym();
2151 }
2152
2153
2154 /** insertEdge. Inserts a constrained edge into the planar mesh. Edges on the lefthand
2155 * side of the edge will be assigned iLeftHandle, others, iRightHandle. All inserted points
2156 * will be assigned the point handle, in case the iPointHandle is different from the NullType.
2157 */
2158 TEMPLATE_DEF
2159 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::insertEdge( const TLineSegment2D& iSegment, EdgeHandle iLeftHandle, EdgeHandle iRightHandle, PointHandle iPointHandle, bool iForcePointHandle, bool makeDelaunay )
2160 {
2161#if DEBUG_MESH
2162 std::cout << "Inserting : " << iSegment << "\n";
2163 static int passCountMesh = 0;
2164 ++passCountMesh;
2165 char buf[128];
2166 sprintf(buf,"debug_mesh_%d.m",passCountMesh);
2167 lass::io::MatlabOStream debugMesh;
2168 debugMesh.open(buf);
2169 debugMesh << *this;
2170 debugMesh.close();
2171#endif
2172 TEdge *ea, *eb;
2173 TPoint2D aa, bb, fbb, faa;
2174 if (ea = insertSite(iSegment.tail()),makeDelaunay)
2175 {
2176 aa = org(ea);
2177 if (iForcePointHandle)
2178 setPointHandle(ea, iPointHandle);
2179
2180 }
2181 if (eb = insertSite(iSegment.head()),makeDelaunay)
2182 {
2183 bb = org(eb);
2184 if (iForcePointHandle)
2185 setPointHandle(eb, iPointHandle);
2186 }
2187 fbb = bb;
2188 faa = aa;
2189
2190 if (ea == NULL || eb == NULL)
2191 {
2192 LASS_THROW("insertEdge: could not insert endpoints of edge");
2193 }
2194
2195 if (org(ea)!=aa)
2196 {
2197 ea=locate(aa);
2198 if (ea==NULL)
2199 {
2200 LASS_THROW("insertEdge: could not locate an endpoint");
2201 }
2202
2203 if (!(aa == org(ea)))
2204 ea = ea->sym();
2205 if (!(aa == org(ea)))
2206 {
2207 LASS_THROW("insertEdge: point a is not an endpoint of ea");
2208 }
2209 }
2210
2211 if (aa==bb)
2212 {
2213 return ea;
2214 }
2215
2216 // first case : edge is present, constrain it
2217 {
2218 TEdge* ce = ea;
2219 int vOrder = vertexOrder(ea);
2220 for (int i=0;i<vOrder;++i)
2221 {
2222 if (dest(ce)==fbb)
2223 {
2224 ce->quadEdge()->edgeConstrain();
2225 setOrientedEdgeHandle( ce, iLeftHandle, iRightHandle, iSegment.vector() );
2226 return ce;
2227 }
2228 ce = ce->oNext();
2229 }
2230 if (distance(aa,bb)<pointDistanceTolerance_ )
2231 {
2232 LASS_THROW("insertEdge: both ends map to same vertex within requested numerical precision");
2233 }
2234 ce = ea;
2235 for (int i=0;i<vOrder;++i)
2236 {
2237 // if is almost in line we also take it
2238 if ( prim::dot(direction(ce), iSegment.vector()) > T(0)
2239 && num::abs(lass::prim::doubleTriangleArea(dest(ce),faa,fbb))<tolerance_)
2240 {
2241 ce->quadEdge()->edgeConstrain();
2242 setOrientedEdgeHandle( ce, iLeftHandle, iRightHandle, iSegment.vector() );
2243
2244 return insertEdge( TLineSegment2D( dest(ce), fbb), iLeftHandle, iRightHandle, iPointHandle, iForcePointHandle, makeDelaunay );
2245 }
2246
2247 ce = ce->oNext();
2248 }
2249 }
2250
2251 std::vector< TEdge* > insertedEdges; // edges having as origin newly inserted points
2252 std::vector< TPoint2D > insertedPoints; // newly inserted points
2253 std::vector< TPoint2D > finalInsertedPoints; // newly inserted points
2254 std::vector< TEdge* > crossedEdges; // crossed edges by new segment
2255 std::vector< TEdge* > filteredEdges; // crossed edges by new segment
2256 //bool tookAction = false;
2257
2258 LASS_ASSERT(allEqualChainOrder(ea));
2259 LASS_ASSERT(allEqualChainOrder(eb));
2260 pointWalk(TLineSegment2D(faa,fbb),std::back_inserter(crossedEdges));
2261#if DEBUG_MESH
2262 std::cout << "Walk returned segments " << crossedEdges.size() << "\n";
2263 std::cout << std::setprecision(15);
2264#endif
2265
2266 // search all edges to split
2267 insertedPoints.push_back(faa);
2268 for (size_t i=0;i<crossedEdges.size();++i)
2269 {
2270 bool computeIntersection = (crossedEdges[i]->isConstrained());
2271 TPoint2D eorg = org(crossedEdges[i]);
2272 TPoint2D edest= dest(crossedEdges[i]);
2273 bool noCommonPoints = eorg!=faa && eorg!=fbb && edest!=faa && edest!=fbb;
2274 if ( computeIntersection
2275 && noCommonPoints)
2276 {
2277 TLine2D other(eorg,edest);
2278 TPoint2D x;
2279
2280 if (prim::intersect(TLine2D(aa,bb),other,x)==prim::rOne)
2281 {
2282 TEdge* npe = insertSite(x,true,true);
2283 TPoint2D nx = org(npe);
2284
2285 // the reverse test now for degeneration avoidance
2286 // this avoids that the crossed edges is only touching the walked path
2287 // noCommonPoints = nx != eorg && nx != edest;
2288
2289 if (nx!=insertedPoints.back())
2290 {
2291 insertedPoints.push_back(nx);
2292#if DEBUG_MESH
2293 std::cout << "-> Intersection " << x << " filtered to " << nx << "\n";
2294#endif
2295 break;
2296 }
2297 }
2298 }
2299 if (noCommonPoints)
2300 filteredEdges.push_back(crossedEdges[i]);
2301 }
2302 if (fbb!=insertedPoints.back())
2303 insertedPoints.push_back(fbb);
2304
2305 // we found some intersections, recursively insert the subparts of the line segments
2306 if (insertedPoints.size()>2)
2307 {
2308 TEdge* ce = NULL;
2309 TEdge* bce = NULL;
2310 for (size_t i=0;i<insertedPoints.size()-1;++i)
2311 {
2312 ce = insertEdge( TLineSegment2D(insertedPoints[i],insertedPoints[i+1]), iLeftHandle, iRightHandle, iPointHandle, iForcePointHandle, makeDelaunay );
2313 if (i==0)
2314 bce = ce;
2315 }
2316 if ( org(bce) != faa )
2317 {
2318 bce = pointLocate( faa );
2319 }
2320 LASS_ASSERT( org(bce)==faa);
2321
2322 return bce;
2323 }
2324
2325 std::vector< TEdge* > toSwapEdges; // edges which are not swappable yet will be hold back till the end
2326 // this will only happen once, if not we have a problem and we retry
2327 // clearing the path recursively
2328 bool madeConstraint = false;
2329 // no intersections found, we clear the path by swapping edges
2330 for (size_t i=0;i<filteredEdges.size();++i)
2331 {
2332 // is the edge swappable without creating an inside-out triangle?
2333 TPoint2D a = org(filteredEdges[i]);
2334 TPoint2D b = dest(filteredEdges[i]);
2335 TPoint2D c = dest(filteredEdges[i]->lNext());
2336 TPoint2D d = dest(filteredEdges[i]->sym()->lNext());
2337 bool swappedCcw1 = prim::ccw(d,b,c);
2338 bool swappedCcw2 = prim::ccw(a,d,c);
2339 if (!swappedCcw1 || !swappedCcw2)
2340 {
2341#if DEBUG_MESH
2342 lass::io::MatlabOStream bugMesh;
2343 bugMesh.open( "swap_constrained.m" );
2344 bugMesh << *this;
2345 for (int j=0;j<filteredEdges.size();++j)
2346 {
2347 bugMesh.setColor( lass::io::mcBlue);
2348 bugMesh << TLineSegment2D(org(filteredEdges[j]),dest(filteredEdges[j]));
2349 }
2350 bugMesh.setColor( lass::io::mcRed );
2351 bugMesh << TLineSegment2D(a,b);
2352 bugMesh.setColor( lass::io::mcGreen );
2353 bugMesh << TLineSegment2D(faa,fbb);
2354 bugMesh.close();
2355#endif
2356 toSwapEdges.push_back(filteredEdges[i]);
2357
2358 //LASS_THROW("could not swap edge without creating inside-out triangle");
2359 }
2360 else
2361 {
2362 if (filteredEdges[i]->quadEdge()->isConstrained())
2363 {
2364 LASS_THROW("found constrained edge where I should not find one!");
2365 }
2366 swap(filteredEdges[i]);
2367 TPoint2D eorg = org(filteredEdges[i]);
2368 TPoint2D edest = dest(filteredEdges[i]);
2369 T onEdgeOrg = num::abs(prim::doubleTriangleArea(eorg,edest,faa));
2370 T onEdgeDest = num::abs(prim::doubleTriangleArea(eorg,edest,fbb));
2371 //T reverseTest1 = num::abs(prim::doubleTriangleArea(eorg,faa,fbb));
2372 //T reverseTest2 = num::abs(prim::doubleTriangleArea(edest,faa,fbb));
2373 if (onEdgeOrg<tolerance_ && onEdgeDest<tolerance_)
2374 {
2375 filteredEdges[i]->quadEdge()->edgeConstrain();
2376 #if DEBUG_MESH
2377 std::cout << "@ Constraining " << org(filteredEdges[i]) << "-->" << dest(filteredEdges[i]) << "\n";
2378 #endif
2379 setOrientedEdgeHandle( filteredEdges[i], iLeftHandle, iRightHandle, iSegment.vector() );
2380 if (toSwapEdges.size()>0)
2381 {
2382 LASS_THROW("edges for delayed swapping found");
2383 }
2384 madeConstraint = true;
2385 break;
2386 }
2387 }
2388 }
2389
2390 if (toSwapEdges.size()>0)
2391 return insertEdge( iSegment, iLeftHandle, iRightHandle, iPointHandle, iForcePointHandle, makeDelaunay );
2392 if (!madeConstraint)
2393 {
2394 LASS_THROW("could not force constrained edge");
2395 }
2396 LASS_ASSERT( org(ea) == faa );
2397 return ea;
2398 }
2399
2400
2401 TEMPLATE_DEF
2402 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::insertPolygon( const TSimplePolygon2D& iPolygon, EdgeHandle iLeftHandle, EdgeHandle iRightHandle, bool /*makeDelaunay*/)
2403 {
2404 TEdge* e = NULL;
2405 for (size_t i=1;i<iPolygon.size();++i)
2406 {
2407 e=insertEdge(TLineSegment2D(iPolygon[i-1],iPolygon[i]),iLeftHandle,iRightHandle);
2408 }
2409 e=insertEdge(TLineSegment2D(iPolygon[iPolygon.size()-1],iPolygon[0]),iLeftHandle,iRightHandle);
2410 return e;
2411 }
2412
2413 /** marks all faces which have their barycentrum inside the given polygon with the provided handle */
2414 TEMPLATE_DEF
2415 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::markPolygon( TEdge* iStartEdge, const TSimplePolygon2D& iPolygon, FaceHandle iFaceHandle)
2416 {
2417 typedef impl::EdgeMarker<T, PointHandle, EdgeHandle, FaceHandle> TEdgeMarker;
2418 StackIncrementer( &stackDepth_, PLANAR_MESH_STACK_DEPTH );
2419 TEdgeMarker edgeMarker( this, false );
2420 forAllPrimaryEdges( TEdgeCallback( &edgeMarker, &TEdgeMarker::internalMark ) );
2421#pragma LASS_TODO("do a more efficient marking or allow for multiple marking")
2422
2423 int n = chainOrder(iStartEdge);
2424 for (int i=0;i<n;++i)
2425 {
2426 floodPolygon( iStartEdge->sym(), iPolygon, iFaceHandle );
2427 iStartEdge = iStartEdge->lNext();
2428 }
2429 }
2430
2431 /** marks all faces which have their barycentrum inside the given polygon with the provided handle */
2432 TEMPLATE_DEF
2433 template <typename InputPolygonIterator, typename InputFaceHandleIterator>
2434 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::markPolygons( InputPolygonIterator iFirstPolygon, InputPolygonIterator iLastPolygon, InputFaceHandleIterator iFirstFaceHandle )
2435 {
2436 typedef impl::EdgeMarker<T, PointHandle, EdgeHandle, FaceHandle> TEdgeMarker;
2437 StackIncrementer( &stackDepth_, PLANAR_MESH_STACK_DEPTH );
2438 TEdgeMarker edgeMarker( this, false );
2439 forAllPrimaryEdges( TEdgeCallback( &edgeMarker, &TEdgeMarker::internalMark ) );
2440
2441 while (iFirstPolygon != iLastPolygon)
2442 {
2443 TEdge* iStartEdge = locate((*iFirstPolygon)[0]);
2444 int n = chainOrder(iStartEdge);
2445 for (int i=0;i<n;++i)
2446 {
2447 floodPolygon( iStartEdge->sym(), *iFirstPolygon, *iFirstFaceHandle );
2448 iStartEdge = iStartEdge->lNext();
2449 }
2450 ++iFirstPolygon;
2451 ++iFirstFaceHandle;
2452 }
2453 }
2454
2455 /** marks all faces which have their barycentrum inside the given polygon with the provided handle */
2456 TEMPLATE_DEF
2457 template <typename InputPolygonIterator, typename InputFaceHandleIterator>
2458 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::forAllPolygonFaces( InputPolygonIterator iFirstPolygon, InputPolygonIterator iLastPolygon, InputFaceHandleIterator iFirstFaceHandle, const TEdgePolyFaceHandleCallback& iCallback )
2459 {
2460 typedef impl::EdgeMarker<T, PointHandle, EdgeHandle, FaceHandle> TEdgeMarker;
2461 StackIncrementer( &stackDepth_, PLANAR_MESH_STACK_DEPTH );
2462 TEdgeMarker edgeMarker( this, false );
2463 forAllPrimaryEdges( TEdgeCallback( &edgeMarker, &TEdgeMarker::internalMark ) );
2464
2465 while (iFirstPolygon != iLastPolygon)
2466 {
2467 TEdge* iStartEdge = locate((*iFirstPolygon)[0],*lastFloodEdge_);
2468 int n = chainOrder(iStartEdge);
2469 for (int i=0;i<n;++i)
2470 {
2471 if (!floodPolygonCallback( iStartEdge->sym(), *iFirstPolygon, *iFirstFaceHandle, iCallback ))
2472 {
2473 return false;
2474 }
2475 if (!floodPolygonCallback( iStartEdge, *iFirstPolygon, *iFirstFaceHandle, iCallback ))
2476 {
2477 return false;
2478 }
2479 iStartEdge = iStartEdge->lNext();
2480 }
2481 ++iFirstPolygon;
2482 ++iFirstFaceHandle;
2483 }
2484 std::cout << "Out polygon faces true" << std::endl;
2485 return true;
2486 }
2487
2488
2489 /** marks all faces which have their barycentrum inside the given polygon with the provided handle */
2490 TEMPLATE_DEF
2491 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::markPolygons( const std::vector<TSimplePolygon2D>& iPolygons, const std::vector<FaceHandle>& iFaceHandles)
2492 {
2493 if (iPolygons.size()!=iFaceHandles.size())
2494 {
2495 LASS_THROW("markPolygons: list of polygons must fit list of face handles");
2496 }
2497 markPolygons(iPolygons.begin(), iPolygons.end(), iFaceHandles.begin());
2498 }
2499
2500
2501 TEMPLATE_DEF
2502 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::floodPolygon( TEdge* iStartEdge, const TSimplePolygon2D& iPolygon, FaceHandle iFaceHandle)
2503 {
2504 TPoint2D bary = polygon(iStartEdge).surfaceCentroid().affine();
2505 if (iPolygon.contains(bary) && !internalMarking(iStartEdge))
2506 {
2507 setInternalMarking( iStartEdge, true );
2508 setFaceHandle( iStartEdge, iFaceHandle );
2509 int n = chainOrder(iStartEdge);
2510 for (int i=0;i<n;++i)
2511 {
2512 if (!floodPolygon( iStartEdge->sym(), iPolygon, iFaceHandle))
2513 return false;
2514 iStartEdge = iStartEdge->lNext();
2515 }
2516 }
2517 return true;
2518 }
2519
2520 TEMPLATE_DEF
2521 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::floodPolygonCallback( TEdge* iStartEdge, const TSimplePolygon2D& iPolygon, FaceHandle iFaceHandle, const TEdgePolyFaceHandleCallback& iCallback )
2522 {
2523 *lastFloodEdge_ = iStartEdge;
2524 TPoint2D bary = polygon(iStartEdge).surfaceCentroid().affine();
2525 if (iPolygon.contains(bary) && !internalMarking(iStartEdge))
2526 {
2527 setInternalMarking( iStartEdge, true );
2528 if (iCallback)
2529 if (!iCallback(iStartEdge, iPolygon,iFaceHandle))
2530 {
2531 return false;
2532 }
2533 int n = chainOrder(iStartEdge);
2534 for (int i=0;i<n;++i)
2535 {
2536 if (!floodPolygonCallback( iStartEdge->sym(), iPolygon, iFaceHandle, iCallback ))
2537 {
2538 return false;
2539 }
2540 iStartEdge = iStartEdge->lNext();
2541 }
2542 }
2543 return true;
2544 }
2545
2546 TEMPLATE_DEF
2547 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::makeEdge( const TPoint2D& a, const TPoint2D& b, bool makeConstrained )
2548 {
2549 TEdge* e(makeEmptyEdge( makeConstrained ));
2550
2551 TPoint2D* na = make(a);
2552 TPoint2D* nb = make(b);
2553 e->handle().point_ = na;
2554 e->sym()->handle().point_ = nb;
2555
2556 return e;
2557 }
2558
2559 TEMPLATE_DEF
2560 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TEdge* PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::startEdge() const
2561 {
2562 return startEdge_;
2563 }
2564
2565
2566 TEMPLATE_DEF
2567 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::inPrimaryMesh( const TEdge* iEdge )
2568 {
2569 return (iEdge->index() == 2) || (iEdge->index() == 0);
2570 }
2571
2572 TEMPLATE_DEF
2573 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::inDualMesh( const TEdge* iEdge )
2574 {
2575 return !inPrimaryMesh( iEdge );
2576 }
2577
2578 TEMPLATE_DEF
2579 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TTriangle2D PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::triangle( const TEdge* iEdge)
2580 {
2581 if (!inPrimaryMesh(iEdge))
2582 {
2583 LASS_THROW("PlanarMesh::triangle: edge not in primary mesh");
2584 }
2585 return TTriangle2D(org(iEdge),dest(iEdge),dest(iEdge->lNext()));
2586 }
2587
2588 TEMPLATE_DEF
2589 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TSimplePolygon2D PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::polygon( const TEdge* iEdge)
2590 {
2591 if (!inPrimaryMesh(iEdge))
2592 {
2593 LASS_THROW("PlanarMesh::polygon: edge not in primary mesh");
2594 }
2595 TSimplePolygon2D poly;
2596 const TEdge* cEdge = iEdge;
2597 do
2598 {
2599 poly.add(org(cEdge));
2600 cEdge = cEdge->lNext();
2601 } while (cEdge!=iEdge);
2602 return poly;
2603 }
2604
2605 TEMPLATE_DEF
2606 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::inConvexCell( const TPoint2D& iPoint, const TEdge* iEdge)
2607 {
2608 if (!inPrimaryMesh(iEdge))
2609 {
2610 LASS_THROW("PlanarMesh::polygon: edge not in primary mesh");
2611 }
2612 const TEdge* cEdge = iEdge;
2613 do
2614 {
2615 if (rightOf(iPoint, cEdge))
2616 return false;
2617 cEdge = cEdge->lNext();
2618 } while (cEdge!=iEdge);
2619 return true;
2620 }
2621
2622 TEMPLATE_DEF
2623 const typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TPoint2D& PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::org( const TEdge* iEdge )
2624 {
2625 if (!inPrimaryMesh( iEdge ))
2626 {
2627 LASS_THROW("PlanarMesh::org: edge not in primary mesh");
2628 }
2629 return *iEdge->handle().point_;
2630 }
2631
2632 TEMPLATE_DEF
2633 const typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TPoint2D& PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::fastOrg( const TEdge* iEdge )
2634 {
2635 return *iEdge->handle().point_;
2636 }
2637
2638 TEMPLATE_DEF
2639 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TPoint2D PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::along( const TEdge* iEdge, const T& iParam )
2640 {
2641 LASS_ASSERT(inPrimaryMesh( iEdge ));
2642 /*
2643 if (!inPrimaryMesh( iEdge ))
2644 {
2645 LASS_THROW("PlanarMesh::along: edge not in primary mesh");
2646 }
2647 */
2648 const TPoint2D& eOrg = *iEdge->handle().point_;
2649 const TPoint2D& eDest = *iEdge->sym()->handle().point_;
2650 const T oX = eDest.x*iParam + (T(1)-iParam)*eOrg.x;
2651 const T oY = eDest.y*iParam + (T(1)-iParam)*eOrg.y;
2652 return TPoint2D(oX,oY);
2653 }
2654
2655 TEMPLATE_DEF
2656 typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TPoint2D PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::fastAlong( const TEdge* iEdge, const T& iParam )
2657 {
2658 const TPoint2D& eOrg = *iEdge->handle()->point_;
2659 const TPoint2D& eDest = *iEdge->sym()->handle()->point_;
2660 const T oX = eDest.x*iParam + (T(1)-iParam)*eOrg.x;
2661 const T oY = eDest.y*iParam + (T(1)-iParam)*eOrg.y;
2662 return TPoint2D(oX,oY);
2663 }
2664
2665
2666 TEMPLATE_DEF inline
2667 const typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TPoint2D& PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::dest( const TEdge* iEdge )
2668 {
2669 return org( iEdge->sym() );
2670 }
2671
2672 TEMPLATE_DEF inline
2673 const typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TPoint2D& PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::fastDest( const TEdge* iEdge )
2674 {
2675 return fastOrg( iEdge->sym() );
2676 }
2677
2678 TEMPLATE_DEF inline
2679 const typename PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::TVector2D PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::direction( const TEdge* iEdge )
2680 {
2681 return dest( iEdge ) - org( iEdge );
2682 }
2683
2684 TEMPLATE_DEF
2685 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::setOrg( const TPoint2D& iPoint, TEdge* iEdge )
2686 {
2687 LASS_ASSERT( inPrimaryMesh( iEdge ) );
2688 *iEdge->handle().point_ = iPoint;
2689 }
2690
2691 TEMPLATE_DEF
2692 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::setDest( const TPoint2D& iPoint, TEdge* iEdge )
2693 {
2694 LASS_ASSERT( inPrimaryMesh( iEdge ) );
2695 *iEdge->sym()->handle().point_ = iPoint;
2696 }
2697
2698 TEMPLATE_DEF
2699 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::rightOf( const TPoint2D& iPoint, const TEdge* iEdge )
2700 {
2701 return lass::prim::ccw( iPoint, dest(iEdge), org(iEdge) );
2702 }
2703
2704 TEMPLATE_DEF
2705 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::fastRightOf( const TPoint2D& iPoint, const TEdge* iEdge )
2706 {
2707 return lass::prim::ccw( iPoint, fastDest(iEdge), fastOrg(iEdge) );
2708 }
2709
2710
2711 TEMPLATE_DEF
2712 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::leftOf( const TPoint2D& iPoint, const TEdge* iEdge )
2713 {
2714 return lass::prim::ccw( iPoint, org(iEdge), dest(iEdge) );
2715 }
2716
2717 TEMPLATE_DEF
2718 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::fastLeftOf( const TPoint2D& iPoint, const TEdge* iEdge )
2719 {
2720 return lass::prim::ccw( iPoint, fastOrg(iEdge), fastDest(iEdge) );
2721 }
2722
2723 TEMPLATE_DEF
2724 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::onEdge( const TPoint2D& iPoint, const TEdge* iEdge )
2725 {
2726 if (prim::doubleTriangleArea( iPoint, org(iEdge), dest(iEdge)) == T(0) )
2727 {
2728 TLineSegment2D eSeg(org(iEdge), dest(iEdge) );
2729 T rt = eSeg.t(iPoint);
2730 return (rt>=T(0)) && (rt<=T(1));
2731 }
2732 return false;
2733 }
2734 TEMPLATE_DEF
2735 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::hasLeftFace( const TEdge* e )
2736 {
2737//#pragma LASS_TODO("This only works for triangular faces... update to general convex polygonal faces")
2738// return ( org(e->lPrev()) == dest(e->lNext()) &&
2739// leftOf(org(e->lPrev()), e));
2740 return leftOf(dest(e->lNext()), e);
2741 }
2742
2743 TEMPLATE_DEF
2744 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::fastHasLeftFace( const TEdge* e )
2745 {
2746//#pragma LASS_TODO("This only works for triangular faces... update to general convex polygonal faces")
2747// return ( org(e->lPrev()) == dest(e->lNext()) &&
2748// leftOf(org(e->lPrev()), e));
2749 return fastLeftOf(fastDest(e->lNext()), e);
2750 }
2751
2752 TEMPLATE_DEF
2753 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::hasRightFace( const TEdge* iEdge )
2754 {
2755 return hasLeftFace( iEdge->sym() );
2756 }
2757
2758 /** chainOrder. returns the order of the polygonal chain starting from iEdge and
2759 * walking around the left-face of iEdge. Or in other words: the number of vertices
2760 * in the polygon on the left of iEdge.
2761 */
2762 TEMPLATE_DEF
2763 int PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::chainOrder( const TEdge* iEdge )
2764 {
2765 int order = 0;
2766 const TEdge* currentEdge = iEdge;
2767 do
2768 {
2769 ++order;
2770 currentEdge = currentEdge->lNext();
2771 } while (currentEdge!=iEdge);
2772 return order;
2773 }
2774
2775 /** allEqualchainOrder. returns true when the chainorder of all edges which share org(iEdge)
2776 * is equal, f.i. when they are all triangles/quadriliterals
2777 *
2778 */
2779 TEMPLATE_DEF
2780 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::allEqualChainOrder( const TEdge* iEdge )
2781 {
2782 int chOrd = chainOrder(iEdge);
2783 for (int i=1;i<vertexOrder(iEdge);++i)
2784 {
2785 iEdge = iEdge->oNext();
2786 if (chainOrder(iEdge)!=chOrd)
2787 return false;
2788 }
2789 return true;
2790 }
2791
2792
2793 /** vertexOrder. returns the order of the vertex defined by the origin of iEdge
2794 * This is the number of undirected edges connected to this origin.
2795 */
2796 TEMPLATE_DEF
2797 int PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::vertexOrder( const TEdge* iEdge )
2798 {
2799 int order = 0;
2800 const TEdge* currentEdge = iEdge;
2801 do
2802 {
2803 ++order;
2804 currentEdge = currentEdge->oNext();
2805 } while (currentEdge!=iEdge);
2806 return order;
2807 }
2808
2809
2810 TEMPLATE_DEF
2811 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::marking( const TEdge* iEdge )
2812 {
2813 return iEdge->handle().internalMark_[publicMarkIndex];
2814 }
2815
2816 TEMPLATE_DEF
2817 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::internalMarking( const TEdge* iEdge )
2818 {
2819 return iEdge->handle().internalMark_[stackDepth_];
2820 }
2821
2822
2823 TEMPLATE_DEF
2824 PointHandle PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::pointHandle( const TEdge* iEdge )
2825 {
2826 return inPrimaryMesh( iEdge ) ? iEdge->handle().pointHandle() : PointHandle();
2827 }
2828
2829 TEMPLATE_DEF
2830 EdgeHandle PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::edgeHandle( const TEdge* iEdge )
2831 {
2832 return inPrimaryMesh( iEdge ) ? iEdge->handle().edgeHandle() : EdgeHandle();
2833 }
2834
2835 TEMPLATE_DEF
2836 FaceHandle PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::faceHandle(const TEdge* iEdge )
2837 {
2838 if ( inPrimaryMesh( iEdge ) )
2839 {
2840 return iEdge->rot()->handle().faceHandle();
2841 }
2842 return FaceHandle();
2843 }
2844
2845
2846 TEMPLATE_DEF
2847 PointHandle& PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::pointHandleRef( TEdge* iEdge )
2848 {
2849 if (inPrimaryMesh(iEdge))
2850 return iEdge->handle().pointHandle();
2851 LASS_THROW("no point handle available");
2852 }
2853
2854 TEMPLATE_DEF
2855 EdgeHandle& PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::edgeHandleRef( TEdge* iEdge )
2856 {
2857 if (inPrimaryMesh(iEdge))
2858 return iEdge->handle().edgeHandle();
2859 LASS_THROW("no edge handle available");
2860 }
2861
2862 TEMPLATE_DEF
2863 FaceHandle& PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::faceHandleRef( TEdge* iEdge )
2864 {
2865 if ( inPrimaryMesh( iEdge ) )
2866 {
2867 return iEdge->rot()->handle().faceHandle();
2868 }
2869 LASS_THROW("no face handle available");
2870 }
2871
2872 TEMPLATE_DEF
2873 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::setMarking( TEdge* iEdge, bool iMark )
2874 {
2875 iEdge->handle().internalMark_.set(publicMarkIndex, iMark);
2876 }
2877
2878 TEMPLATE_DEF
2879 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::setInternalMarking( TEdge* iEdge, bool iMark )
2880 {
2881 iEdge->handle().internalMark_.set(stackDepth_, iMark);
2882 }
2883
2884 TEMPLATE_DEF
2885 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::setInternalMarkingAroundVertex( TEdge* iEdge, bool iMark )
2886 {
2887 TEdge* e = iEdge;
2888 do
2889 {
2890 setInternalMarking( e, iMark );
2891 e = e->oNext();
2892 } while (e!=iEdge);
2893 }
2894
2895 TEMPLATE_DEF
2896 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::setInternalMarkingInFace( TEdge* iEdge, bool iMark )
2897 {
2898 TEdge* e = iEdge;
2899 do
2900 {
2901 setInternalMarking( e, iMark );
2902 e = e->lNext();
2903 } while (e!=iEdge);
2904 }
2905
2906
2907 TEMPLATE_DEF
2908 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::setPointHandle( TEdge* iEdge, PointHandle iHandle )
2909 {
2910 if (! inPrimaryMesh( iEdge ) )
2911 throw std::runtime_error("setPointHandle : edge not in primary mesh");
2912
2913 TEdge* currentEdge = iEdge;
2914 do
2915 {
2916 currentEdge->handle().pointHandle() = iHandle;
2917 currentEdge = currentEdge->oNext();
2918 } while ( currentEdge != iEdge );
2919 }
2920
2921 TEMPLATE_DEF
2922 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::setEdgeHandle( TEdge* iEdge, EdgeHandle iHandle )
2923 {
2924 if (! inPrimaryMesh( iEdge ) )
2925 throw std::runtime_error("setEdgeHandle : edge not in primary mesh");
2926 iEdge->handle().edgeHandle() = iHandle;
2927 }
2928
2929 TEMPLATE_DEF
2930 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::setFaceHandle( TEdge* iEdge, FaceHandle iHandle )
2931 {
2932 if (! inPrimaryMesh( iEdge ) )
2933 throw std::runtime_error("setFaceHandle : edge not in dual mesh");
2934
2935 TEdge* currentEdge = iEdge;
2936 do
2937 {
2938 currentEdge->rot()->handle().faceHandle() = iHandle;
2939 if ( faceHandle( currentEdge ) != faceHandle( currentEdge->sym() ) )
2940 currentEdge->quadEdge()->faceConstrain();
2941 else
2942 currentEdge->quadEdge()->faceDeconstrain();
2943 currentEdge = currentEdge->lNext();
2944 } while ( currentEdge != iEdge );
2945 }
2946
2947 TEMPLATE_DEF
2948 void PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::setOrientedEdgeHandle(
2949 TEdge* iEdge, EdgeHandle iLeftHandle, EdgeHandle iRightHandle, const TVector2D& iOrientation )
2950 {
2951#ifndef NDEBUG
2952 ++numSetOrientedEdgeHandleCalls;
2953#endif
2954 if (prim::dot(direction(iEdge), iOrientation) < T(0))
2955 {
2956#ifndef NDEBUG
2957 ++numSetOrientedEdgeHandleSwaps;
2958#endif
2959 std::swap(iLeftHandle, iRightHandle);
2960 }
2961 setEdgeHandle(iEdge, iLeftHandle);
2962 setEdgeHandle(iEdge->sym(), iRightHandle);
2963 }
2964
2965 TEMPLATE_DEF
2966 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::forAllPrimaryEdges( const TEdgeCallback& iCallback )
2967 {
2968 typename TQuadEdgeList::iterator qIt;
2969 for (qIt = quadEdgeList_.begin(); qIt != quadEdgeList_.end();++qIt)
2970 {
2971 if (!(*qIt)->deleted)
2972 {
2973 if (inPrimaryMesh( (*qIt)->edges() ) )
2974 {
2975 if (!iCallback( (*qIt)->edges() )) return false;
2976 if (!iCallback( (*qIt)->edges()->sym() )) return false;
2977 }
2978 else
2979 {
2980 if (!iCallback( (*qIt)->edges()->rot() )) return false;
2981 if (!iCallback( (*qIt)->edges()->invRot() )) return false;
2982 }
2983 }
2984 }
2985 return true;
2986 }
2987
2988 /** forAllPrimaryUndirectedEdges. Edges sharing endpoints are only listed once */
2989 TEMPLATE_DEF
2990 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::forAllPrimaryUndirectedEdges( const TEdgeCallback& iCallback )
2991 {
2992 typename TQuadEdgeList::iterator qIt;
2993 for (qIt = quadEdgeList_.begin(); qIt != quadEdgeList_.end();++qIt)
2994 {
2995 if (!(*qIt)->deleted)
2996 {
2997 if (inPrimaryMesh( (*qIt)->edges() ) )
2998 {
2999 if (!iCallback( (*qIt)->edges() )) return false;
3000 }
3001 else
3002 {
3003 if (!iCallback( (*qIt)->edges()->rot() )) return false;
3004 }
3005 }
3006 }
3007 return true;
3008 }
3009 /** forAllPrimaryUndirectedEdges. Edges sharing endpoints are only listed once. Random starting point */
3010 TEMPLATE_DEF
3011 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::forAllPrimaryUndirectedEdgesCached( const TEdgeCallback& iCallback )
3012 {
3013 static size_t start = 0;
3014 size_t savedStart = start;
3015
3016 if (savedStart>quadEdgeList_.size()-1)
3017 {
3018 savedStart = 0;
3019 start = 0;
3020 }
3021
3022 for (size_t i = savedStart; i<quadEdgeList_.size();++i)
3023 {
3024 TQuadEdge* qe = quadEdgeList_[i];
3025 if (!qe->deleted)
3026 {
3027 if (inPrimaryMesh( (qe)->edges() ) )
3028 {
3029 if (!iCallback( (qe)->edges() )) return false;
3030 }
3031 else
3032 {
3033 if (!iCallback( (qe)->edges()->rot() )) return false;
3034 }
3035 }
3036 ++start;
3037 }
3038 start = 0;
3039 for (size_t i = 0; i<savedStart;++i)
3040 {
3041 TQuadEdge* qe = quadEdgeList_[i];
3042 if (!qe->deleted)
3043 {
3044 if (inPrimaryMesh( (qe)->edges() ) )
3045 {
3046 if (!iCallback( (qe)->edges() )) return false;
3047 }
3048 else
3049 {
3050 if (!iCallback( (qe)->edges()->rot() )) return false;
3051 }
3052 }
3053 ++start;
3054 }
3055 return true;
3056 }
3057
3058 TEMPLATE_DEF
3059 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::forAllDualEdges( const TEdgeCallback& iCallback )
3060 {
3061 typename TQuadEdgeList::iterator qIt;
3062 for (qIt = quadEdgeList_.begin(); qIt != quadEdgeList_.end();++qIt)
3063 {
3064 if (!(*qIt)->deleted)
3065 {
3066 if (inDualMesh( (*qIt)->edges() ) )
3067 {
3068 if (!iCallback( (*qIt)->edges() )) return false;
3069 if (!iCallback( (*qIt)->edges()->sym() )) return false;
3070 }
3071 else
3072 {
3073 if (!iCallback( (*qIt)->edges()->rot() )) return false;
3074 if (!iCallback( (*qIt)->edges()->invRot() )) return false;
3075 }
3076 }
3077 }
3078 return true;
3079 }
3080
3081 TEMPLATE_DEF
3082 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::forAllEdges( const TEdgeCallback& iCallback )
3083 {
3084 typename TQuadEdgeList::iterator qIt;
3085 for (qIt = quadEdgeList_.begin(); qIt != quadEdgeList_.end();++qIt)
3086 {
3087 if (!(*qIt)->deleted)
3088 {
3089 if (!iCallback( (*qIt)->edges() )) return false;
3090 if (!iCallback( (*qIt)->edges()->sym() )) return false;
3091 if (!iCallback( (*qIt)->edges()->rot() )) return false;
3092 if (!iCallback( (*qIt)->edges()->invRot() )) return false;
3093 }
3094 }
3095 return true;
3096 }
3097
3098 TEMPLATE_DEF
3099 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::forAllVertices( const TEdgeCallback& iCallback )
3100 {
3101 StackIncrementer( &stackDepth_, PLANAR_MESH_STACK_DEPTH );
3102 typedef impl::EdgeMarker<T, PointHandle, EdgeHandle, FaceHandle> TEdgeMarker;
3103 TEdgeMarker edgeMarker( this, false );
3104 forAllPrimaryEdges( TEdgeCallback( &edgeMarker, &TEdgeMarker::internalMark ) );
3105
3106 typename TQuadEdgeList::iterator qIt;
3107 for (qIt = quadEdgeList_.begin(); qIt != quadEdgeList_.end();++qIt)
3108 {
3109 if (!(*qIt)->deleted)
3110 {
3111 TEdge* baseEdge = (*qIt)->edges();
3112 if (inDualMesh( baseEdge ) )
3113 baseEdge = baseEdge->rot();
3114
3115 TEdge* e = baseEdge;
3116 for (int i=0;i<2;++i)
3117 {
3118 if ( !internalMarking( e ) )
3119 {
3120 if (!iCallback( e ))
3121 return false;
3122 setInternalMarkingAroundVertex( e, true );
3123 }
3124 e = e->sym();
3125 }
3126 }
3127 }
3128 return true;
3129 }
3130
3131 TEMPLATE_DEF
3132 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::forAllFacesCached( const TEdgeCallback& iCallback )
3133 {
3134 StackIncrementer( &stackDepth_, PLANAR_MESH_STACK_DEPTH );
3135 typedef impl::EdgeMarker<T, PointHandle, EdgeHandle, FaceHandle> TEdgeMarker;
3136 TEdgeMarker edgeMarker( this, false );
3137 forAllPrimaryEdges( TEdgeCallback( &edgeMarker, &TEdgeMarker::internalMark ) );
3138
3139 static size_t startFace = 0;
3140 size_t savedStart = startFace;
3141
3142 if (savedStart>quadEdgeList_.size()-1)
3143 {
3144 savedStart = 0;
3145 startFace = 0;
3146 }
3147
3148 for (size_t i = savedStart; i<quadEdgeList_.size();++i)
3149 {
3150 startFace = i;
3151 TQuadEdge* qe = quadEdgeList_[i];
3152 if (!qe->deleted)
3153 {
3154 TEdge* baseEdge = qe->edges();
3155 if (inDualMesh( baseEdge ) )
3156 baseEdge = baseEdge->rot();
3157
3158 TEdge* e = baseEdge;
3159 for (int i=0;i<2;++i)
3160 {
3161 if ( !internalMarking( e ) )
3162 {
3163 if (!iCallback( e ))
3164 return false;
3165 setInternalMarkingInFace( e, true );
3166 }
3167 e = e->sym();
3168 }
3169 }
3170 }
3171 for (size_t i = 0; i<savedStart;++i)
3172 {
3173 startFace = i;
3174 TQuadEdge* qe = quadEdgeList_[i];
3175 if (!qe->deleted)
3176 {
3177 TEdge* baseEdge = qe->edges();
3178 if (inDualMesh( baseEdge ) )
3179 baseEdge = baseEdge->rot();
3180
3181 TEdge* e = baseEdge;
3182 for (int i=0;i<2;++i)
3183 {
3184 if ( !internalMarking( e ) )
3185 {
3186 if (!iCallback( e ))
3187 return false;
3188 setInternalMarkingInFace( e, true );
3189 }
3190 e = e->sym();
3191 }
3192 }
3193 }
3194 return true;
3195 }
3196
3197 TEMPLATE_DEF
3198 bool PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::forAllFaces( const TEdgeCallback& iCallback )
3199 {
3200 StackIncrementer( &stackDepth_, PLANAR_MESH_STACK_DEPTH );
3201 typedef impl::EdgeMarker<T, PointHandle, EdgeHandle, FaceHandle> TEdgeMarker;
3202 TEdgeMarker edgeMarker( this, false );
3203 forAllPrimaryEdges( TEdgeCallback( &edgeMarker, &TEdgeMarker::internalMark ) );
3204
3205 typename TQuadEdgeList::iterator qIt;
3206 for (qIt = quadEdgeList_.begin(); qIt != quadEdgeList_.end();++qIt)
3207 {
3208 if (!(*qIt)->deleted)
3209 {
3210 TEdge* baseEdge = (*qIt)->edges();
3211 if (inDualMesh( baseEdge ) )
3212 baseEdge = baseEdge->rot();
3213
3214 TEdge* e = baseEdge;
3215 for (int i=0;i<2;++i)
3216 {
3217 if ( !internalMarking( e ) )
3218 {
3219 if (!iCallback( e ))
3220 return false;
3221 setInternalMarkingInFace( e, true );
3222 }
3223 e = e->sym();
3224 }
3225 }
3226 }
3227 return true;
3228 }
3229
3230
3231#ifndef NDEBUG
3232 TEMPLATE_DEF unsigned PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::numSetOrientedEdgeHandleCalls = 0;
3233 TEMPLATE_DEF unsigned PlanarMesh<T, PointHandle, EdgeHandle, FaceHandle>::numSetOrientedEdgeHandleSwaps = 0;
3234#endif
3235
3236
3237}
3238
3239}
3240
3241#endif
aabb_2d.h
callback_r_1.h
callback_r_2.h
callback_r_3.h
lass::io::MatlabOStream
Output stream for writing a selection of geometric primitives to matlab M files.
Definition matlab_o_stream.h:107
lass::prim::LineSegment2D::vector
const TVector vector() const
Return vector from tail to head.
Definition line_segment_2d.inl:144
lass::prim::Ray2D::point
const TPoint point(TParam a_t) const
Return point on ray by it's parameter.
Definition ray_2d.inl:156
lass::prim::Ray2D::t
const TValue t(const TPoint &iPoint) const
Return parameter of point on the ray that is the projection of the given point.
Definition ray_2d.inl:169
lass::prim::Ray2D::support
const TPoint & support() const
return origin of ray.
Definition ray_2d.inl:97
lass::prim::Ray2D::classify
Side classify(const TPoint &iPoint) const
Return on what side a point is located.
Definition ray_2d.inl:246
lass::prim::SimplePolygon2D
convex or concave polygon in 2D (not selfintersecting, no holes)
Definition simple_polygon_2d.h:74
lass::prim::Triangle2D
A very simple 2D polygon :)
Definition triangle_2d.h:65
lass::spat::PlanarMesh::insertEdge
TEdge * insertEdge(const TLineSegment2D &iSegment, EdgeHandle iLeftHandle=EdgeHandle(), EdgeHandle iRightHandle=EdgeHandle(), PointHandle iPointHandle=PointHandle(), bool forcePointHandle=false, bool makeDelaunay=true)
insertEdge.
Definition planar_mesh.h:2159
lass::spat::PlanarMesh::split
TEdge * split(TEdge *a, T iWhere)
splits an edge with iWhere in the open interval (0,1) and returns the new edge with as org the splitt...
Definition planar_mesh.h:1052
lass::spat::PlanarMesh::vertexOrder
static int vertexOrder(const TEdge *iEdge)
vertexOrder.
Definition planar_mesh.h:2797
lass::spat::PlanarMesh::markPolygon
void markPolygon(TEdge *iStartEdge, const TSimplePolygon2D &iPolygon, FaceHandle iFaceHandle)
marks all faces which have their barycentrum inside the given polygon with the provided handle
Definition planar_mesh.h:2415
lass::spat::PlanarMesh::deleteEdge
bool deleteEdge(TEdge *iEdge)
delete edge without using gc
Definition planar_mesh.h:1252
lass::spat::PlanarMesh::locate
TEdge * locate(const TPoint2D &iPoint, TEdge *iHint=0) const
locate an edge of the triangle containing iPoint
Definition planar_mesh.h:1413
lass::spat::PlanarMesh::chainOrder
static int chainOrder(const TEdge *iEdge)
chainOrder.
Definition planar_mesh.h:2763
lass::spat::PlanarMesh::gcDeleteEdge
bool gcDeleteEdge(TEdge *iEdge)
delete edge using garbage collector, useful for deletion avalanches
Definition planar_mesh.h:1291
lass::spat::PlanarMesh::gc
int gc()
do garbage collection after deletion avalanches, returns number of quadedge collected
Definition planar_mesh.h:1325
lass::spat::PlanarMesh::allEqualChainOrder
static bool allEqualChainOrder(const TEdge *iEdge)
allEqualchainOrder.
Definition planar_mesh.h:2780
lass::spat::PlanarMesh::insertSite
TEdge * insertSite(const TPoint2D &iPoint, bool makeDelaunay=true, bool forceOnEdge=false)
makeDelaunay: after insertion of the site force a Delaunay triangulation forceOnEdge : points are con...
Definition planar_mesh.h:1966
lass::spat::PlanarMesh::shoot
TEdge * shoot(const TRay2D &iRay) const
locate the edge found by shooting the ray from within the triangle containt the tail of the ray
Definition planar_mesh.h:1599
lass::spat::PlanarMesh::isBoundingPoint
bool isBoundingPoint(const TPoint2D &iPoint) const
true for points defining the boundary
Definition planar_mesh.h:1399
lass::spat::PlanarMesh::makeEdge
TEdge * makeEdge(const TPoint2D &a, const TPoint2D &b, bool makeConstrained)
makes a dangling edge, this is really only for ADVANCED usage, you'll have to splice it or it won't b...
Definition planar_mesh.h:2547
lass::spat::PlanarMesh::connect
TEdge * connect(TEdge *a, TEdge *b)
connects the dest of a with the org of b
Definition planar_mesh.h:1031
lass::spat::PlanarMesh::pointLocate
TEdge * pointLocate(const TPoint2D &iPoint) const
locate an edge of which the org is the same as iPoint, useful for known degeneracy point location,...
Definition planar_mesh.h:1534
lass::spat::PlanarMesh::forAllPrimaryUndirectedEdgesCached
bool forAllPrimaryUndirectedEdgesCached(const TEdgeCallback &iCallback)
forAllPrimaryUndirectedEdges.
Definition planar_mesh.h:3011
lass::spat::PlanarMesh::forAllPrimaryUndirectedEdges
bool forAllPrimaryUndirectedEdges(const TEdgeCallback &iCallback)
forAllPrimaryUndirectedEdges.
Definition planar_mesh.h:2990
lass::spat::PlanarMesh::markPolygons
void markPolygons(const std::vector< TSimplePolygon2D > &iPolygons, const std::vector< FaceHandle > &iFaceHandles)
marks all faces which have their barycentrum inside the given polygon with the provided handle
Definition planar_mesh.h:2491
extended_iterator.h
lass::num::abs
T abs(const T &x)
if x < 0 return -x, else return x.
Definition basic_ops.h:145
LASS_LOCK
#define LASS_LOCK(iLock)
Locks a iLock and starts a scope block in which it remains locked.
Definition thread.h:619
line_2d.h
line_segment_2d.h
matlab_o_stream.h
lass::prim::Side
Side
Different sides of a surface.
Definition side.h:79
lass::prim::sSurface
@ sSurface
right on the surface
Definition side.h:87
lass::prim::sRight
@ sRight
alias for sBack in 2D
Definition side.h:84
lass::prim::Result
Result
meta information on the result you have from an operation like an intersection ...
Definition result.h:74
lass::prim::rInvalid
@ rInvalid
0 is an invalid value, nothing is known.
Definition result.h:75
lass::prim::rOne
@ rOne
there's exactly one answer, 1 output argument contains the answer
Definition result.h:77
lass::spat
spatial subdivisions, quadtrees, octrees, meshes in 2D and 3D, triangulators, ...
Definition aabb8_tree.h:80
lass
Library for Assembled Shared Sources.
Definition config.h:53
lass::spat::impl::fastIntersect
lass::prim::Result fastIntersect(const lass::prim::Point2D< T > &iA, const lass::prim::Point2D< T > &iB, const lass::prim::Point2D< T > &iC, const lass::prim::Point2D< T > &iD, lass::prim::Point2D< T > &oP)
fastIntersect.
Definition planar_mesh.h:468
point_2d.h
quad_edge.h
ray_2d.h
simple_polygon_2d.h
spat_common.h
triangle_2d.h