KiCad PCB EDA Suite
polygon_triangulation.h
Go to the documentation of this file.
1 /*
2  * This program source code file is part of KiCad, a free EDA CAD application.
3  *
4  * Modifications Copyright (C) 2018-2021 KiCad Developers
5  *
6  * This program is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU General Public License
8  * as published by the Free Software Foundation; either version 2
9  * of the License, or (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, you may find one here:
18  * http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
19  * or you may search the http://www.gnu.org website for the version 2 license,
20  * or you may write to the Free Software Foundation, Inc.,
21  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
22  *
23  * Based on Uniform Plane Subdivision algorithm from Lamot, Marko, and Borut ┼Żalik.
24  * "A fast polygon triangulation algorithm based on uniform plane subdivision."
25  * Computers & graphics 27, no. 2 (2003): 239-253.
26  *
27  * Code derived from:
28  * K-3D which is Copyright (c) 2005-2006, Romain Behar, GPL-2, license above
29  * earcut which is Copyright (c) 2016, Mapbox, ISC
30  *
31  * ISC License:
32  * Permission to use, copy, modify, and/or distribute this software for any purpose
33  * with or without fee is hereby granted, provided that the above copyright notice
34  * and this permission notice appear in all copies.
35  *
36  * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
37  * REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
38  * FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
39  * INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
40  * OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
41  * TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
42  * THIS SOFTWARE.
43  *
44  */
45 
46 #ifndef __POLYGON_TRIANGULATION_H
47 #define __POLYGON_TRIANGULATION_H
48 
49 #include <algorithm>
50 #include <deque>
51 #include <cmath>
52 
53 #include <clipper.hpp>
56 #include <math/box2.h>
57 #include <math/vector2d.h>
58 
60 {
61 public:
63  m_result( aResult )
64  {};
65 
66  bool TesselatePolygon( const SHAPE_LINE_CHAIN& aPoly )
67  {
68  m_bbox = aPoly.BBox();
69  m_result.Clear();
70 
71  if( !m_bbox.GetWidth() || !m_bbox.GetHeight() )
72  return false;
73 
77  Vertex* firstVertex = createList( aPoly );
78  if( !firstVertex || firstVertex->prev == firstVertex->next )
79  return false;
80 
81  firstVertex->updateList();
82 
83  auto retval = earcutList( firstVertex );
84  m_vertices.clear();
85  return retval;
86  }
87 
88 private:
89  struct Vertex
90  {
91  Vertex( size_t aIndex, double aX, double aY, PolygonTriangulation* aParent ) :
92  i( aIndex ),
93  x( aX ),
94  y( aY ),
95  parent( aParent )
96  {
97  }
98 
99  Vertex& operator=( const Vertex& ) = delete;
100  Vertex& operator=( Vertex&& ) = delete;
101 
102  bool operator==( const Vertex& rhs ) const
103  {
104  return this->x == rhs.x && this->y == rhs.y;
105  }
106  bool operator!=( const Vertex& rhs ) const { return !( *this == rhs ); }
107 
108 
120  {
121  parent->m_vertices.emplace_back( i, x, y, parent );
122  Vertex* a2 = &parent->m_vertices.back();
123  parent->m_vertices.emplace_back( b->i, b->x, b->y, parent );
124  Vertex* b2 = &parent->m_vertices.back();
125  Vertex* an = next;
126  Vertex* bp = b->prev;
127 
128  next = b;
129  b->prev = this;
130 
131  a2->next = an;
132  an->prev = a2;
133 
134  b2->next = a2;
135  a2->prev = b2;
136 
137  bp->next = b2;
138  b2->prev = bp;
139 
140  return b2;
141  }
142 
146  void remove()
147  {
148  next->prev = prev;
149  prev->next = next;
150 
151  if( prevZ )
152  prevZ->nextZ = nextZ;
153 
154  if( nextZ )
155  nextZ->prevZ = prevZ;
156 
157  next = nullptr;
158  prev = nullptr;
159  nextZ = nullptr;
160  prevZ = nullptr;
161  }
162 
163  void updateOrder()
164  {
165  if( !z )
166  z = parent->zOrder( x, y );
167  }
168 
173  void updateList()
174  {
175  Vertex* p = next;
176 
177  while( p != this )
178  {
182  if( *p == *p->next )
183  {
184  p = p->prev;
185  p->next->remove();
186 
187  if( p == p->next )
188  break;
189  }
190 
191  p->updateOrder();
192  p = p->next;
193  };
194 
195  updateOrder();
196  zSort();
197  }
198 
202  void zSort()
203  {
204  std::deque<Vertex*> queue;
205 
206  queue.push_back( this );
207 
208  for( auto p = next; p && p != this; p = p->next )
209  queue.push_back( p );
210 
211  std::sort( queue.begin(), queue.end(), []( const Vertex* a, const Vertex* b )
212  {
213  return a->z < b->z;
214  } );
215 
216  Vertex* prev_elem = nullptr;
217 
218  for( auto elem : queue )
219  {
220  if( prev_elem )
221  prev_elem->nextZ = elem;
222 
223  elem->prevZ = prev_elem;
224  prev_elem = elem;
225  }
226 
227  prev_elem->nextZ = nullptr;
228  }
229 
230 
234  bool inTriangle( const Vertex& a, const Vertex& b, const Vertex& c )
235  {
236  return ( c.x - x ) * ( a.y - y ) - ( a.x - x ) * ( c.y - y ) >= 0
237  && ( a.x - x ) * ( b.y - y ) - ( b.x - x ) * ( a.y - y ) >= 0
238  && ( b.x - x ) * ( c.y - y ) - ( c.x - x ) * ( b.y - y ) >= 0;
239  }
240 
241  const size_t i;
242  const double x;
243  const double y;
245 
246  // previous and next vertices nodes in a polygon ring
247  Vertex* prev = nullptr;
248  Vertex* next = nullptr;
249 
250  // z-order curve value
251  int32_t z = 0;
252 
253  // previous and next nodes in z-order
254  Vertex* prevZ = nullptr;
255  Vertex* nextZ = nullptr;
256  };
257 
263  int32_t zOrder( const double aX, const double aY ) const
264  {
265  int32_t x = static_cast<int32_t>( 32767.0 * ( aX - m_bbox.GetX() ) / m_bbox.GetWidth() );
266  int32_t y = static_cast<int32_t>( 32767.0 * ( aY - m_bbox.GetY() ) / m_bbox.GetHeight() );
267 
268  x = ( x | ( x << 8 ) ) & 0x00FF00FF;
269  x = ( x | ( x << 4 ) ) & 0x0F0F0F0F;
270  x = ( x | ( x << 2 ) ) & 0x33333333;
271  x = ( x | ( x << 1 ) ) & 0x55555555;
272 
273  y = ( y | ( y << 8 ) ) & 0x00FF00FF;
274  y = ( y | ( y << 4 ) ) & 0x0F0F0F0F;
275  y = ( y | ( y << 2 ) ) & 0x33333333;
276  y = ( y | ( y << 1 ) ) & 0x55555555;
277 
278  return x | ( y << 1 );
279  }
280 
289  {
290  Vertex* retval = nullptr;
291  Vertex* p = aStart->next;
292 
293  while( p != aStart )
294  {
295  if( area( p->prev, p, p->next ) == 0.0 )
296  {
297  p = p->prev;
298  p->next->remove();
299  retval = aStart;
300 
301  if( p == p->next )
302  break;
303  }
304  p = p->next;
305  };
306 
307  // We needed an end point above that wouldn't be removed, so
308  // here we do the final check for this as a Steiner point
309  if( area( aStart->prev, aStart, aStart->next ) == 0.0 )
310  {
311  retval = p->next;
312  p->remove();
313  }
314 
315  return retval;
316  }
317 
321  Vertex* createList( const ClipperLib::Path& aPath )
322  {
323  Vertex* tail = nullptr;
324  double sum = 0.0;
325  auto len = aPath.size();
326 
327  // Check for winding order
328  for( size_t i = 0; i < len; i++ )
329  {
330  auto p1 = aPath.at( i );
331  auto p2 = aPath.at( ( i + 1 ) < len ? i + 1 : 0 );
332 
333  sum += ( ( p2.X - p1.X ) * ( p2.Y + p1.Y ) );
334  }
335 
336  if( sum <= 0.0 )
337  {
338  for( auto point : aPath )
339  tail = insertVertex( VECTOR2I( point.X, point.Y ), tail );
340  }
341  else
342  {
343  for( size_t i = 0; i < len; i++ )
344  {
345  auto p = aPath.at( len - i - 1 );
346  tail = insertVertex( VECTOR2I( p.X, p.Y ), tail );
347  }
348  }
349 
350  if( tail && ( *tail == *tail->next ) )
351  {
352  tail->next->remove();
353  }
354 
355  return tail;
356 
357  }
358 
363  {
364  Vertex* tail = nullptr;
365  double sum = 0.0;
366 
367  // Check for winding order
368  for( int i = 0; i < points.PointCount(); i++ )
369  {
370  VECTOR2D p1 = points.CPoint( i );
371  VECTOR2D p2 = points.CPoint( i + 1 );
372 
373  sum += ( ( p2.x - p1.x ) * ( p2.y + p1.y ) );
374  }
375 
376  if( sum > 0.0 )
377  for( int i = points.PointCount() - 1; i >= 0; i--)
378  tail = insertVertex( points.CPoint( i ), tail );
379  else
380  for( int i = 0; i < points.PointCount(); i++ )
381  tail = insertVertex( points.CPoint( i ), tail );
382 
383  if( tail && ( *tail == *tail->next ) )
384  {
385  tail->next->remove();
386  }
387 
388  return tail;
389  }
390 
404  bool earcutList( Vertex* aPoint, int pass = 0 )
405  {
406  if( !aPoint )
407  return true;
408 
409  Vertex* stop = aPoint;
410  Vertex* prev;
411  Vertex* next;
412 
413  while( aPoint->prev != aPoint->next )
414  {
415  prev = aPoint->prev;
416  next = aPoint->next;
417 
418  if( isEar( aPoint ) )
419  {
420  m_result.AddTriangle( prev->i, aPoint->i, next->i );
421  aPoint->remove();
422 
423  // Skip one vertex as the triangle will account for the prev node
424  aPoint = next->next;
425  stop = next->next;
426 
427  continue;
428  }
429 
430  Vertex* nextNext = next->next;
431 
432  if( *prev != *nextNext && intersects( prev, aPoint, next, nextNext ) &&
433  locallyInside( prev, nextNext ) &&
434  locallyInside( nextNext, prev ) )
435  {
436  m_result.AddTriangle( prev->i, aPoint->i, nextNext->i );
437 
438  // remove two nodes involved
439  next->remove();
440  aPoint->remove();
441 
442  aPoint = nextNext;
443  stop = nextNext;
444 
445  continue;
446  }
447 
448  aPoint = next;
449 
450  /*
451  * We've searched the entire polygon for available ears and there are still
452  * un-sliced nodes remaining.
453  */
454  if( aPoint == stop )
455  {
456  // First, try to remove the remaining steiner points
457  // If aPoint is a steiner, we need to re-assign both the start and stop points
458  if( auto newPoint = removeNullTriangles( aPoint ) )
459  {
460  aPoint = newPoint;
461  stop = newPoint;
462  continue;
463  }
464 
465  // If we don't have any NULL triangles left, cut the polygon in two and try again
466  splitPolygon( aPoint );
467  break;
468  }
469  }
470 
471  // Check to see if we are left with only three points in the polygon
472  if( aPoint->next && aPoint->prev == aPoint->next->next )
473  {
474  // Three concave points will never be able to be triangulated because they were
475  // created by an intersecting polygon, so just drop them.
476  if( area( aPoint->prev, aPoint, aPoint->next ) >= 0 )
477  return true;
478  }
479 
480  /*
481  * At this point, our polygon should be fully tessellated.
482  */
483  return( aPoint->prev == aPoint->next );
484  }
485 
495  bool isEar( Vertex* aEar ) const
496  {
497  const Vertex* a = aEar->prev;
498  const Vertex* b = aEar;
499  const Vertex* c = aEar->next;
500 
501  // If the area >=0, then the three points for a concave sequence
502  // with b as the reflex point
503  if( area( a, b, c ) >= 0 )
504  return false;
505 
506  // triangle bbox
507  const double minTX = std::min( a->x, std::min( b->x, c->x ) );
508  const double minTY = std::min( a->y, std::min( b->y, c->y ) );
509  const double maxTX = std::max( a->x, std::max( b->x, c->x ) );
510  const double maxTY = std::max( a->y, std::max( b->y, c->y ) );
511 
512  // z-order range for the current triangle bounding box
513  const int32_t minZ = zOrder( minTX, minTY );
514  const int32_t maxZ = zOrder( maxTX, maxTY );
515 
516  // first look for points inside the triangle in increasing z-order
517  Vertex* p = aEar->nextZ;
518 
519  while( p && p->z <= maxZ )
520  {
521  if( p != a && p != c
522  && p->inTriangle( *a, *b, *c )
523  && area( p->prev, p, p->next ) >= 0 )
524  return false;
525 
526  p = p->nextZ;
527  }
528 
529  // then look for points in decreasing z-order
530  p = aEar->prevZ;
531 
532  while( p && p->z >= minZ )
533  {
534  if( p != a && p != c
535  && p->inTriangle( *a, *b, *c )
536  && area( p->prev, p, p->next ) >= 0 )
537  return false;
538 
539  p = p->prevZ;
540  }
541 
542  return true;
543  }
544 
551  void splitPolygon( Vertex* start )
552  {
553  Vertex* origPoly = start;
554 
555  do
556  {
557  Vertex* marker = origPoly->next->next;
558 
559  while( marker != origPoly->prev )
560  {
561  // Find a diagonal line that is wholly enclosed by the polygon interior
562  if( origPoly->i != marker->i && goodSplit( origPoly, marker ) )
563  {
564  Vertex* newPoly = origPoly->split( marker );
565 
566  origPoly->updateList();
567  newPoly->updateList();
568 
569  earcutList( origPoly );
570  earcutList( newPoly );
571  return;
572  }
573 
574  marker = marker->next;
575  }
576 
577  origPoly = origPoly->next;
578  } while( origPoly != start );
579  }
580 
589  bool goodSplit( const Vertex* a, const Vertex* b ) const
590  {
591  return a->next->i != b->i &&
592  a->prev->i != b->i &&
593  !intersectsPolygon( a, b ) &&
594  locallyInside( a, b );
595  }
596 
600  double area( const Vertex* p, const Vertex* q, const Vertex* r ) const
601  {
602  return ( q->y - p->y ) * ( r->x - q->x ) - ( q->x - p->x ) * ( r->y - q->y );
603  }
604 
610  bool intersects( const Vertex* p1, const Vertex* q1, const Vertex* p2, const Vertex* q2 ) const
611  {
612  if( ( *p1 == *q1 && *p2 == *q2 ) || ( *p1 == *q2 && *p2 == *q1 ) )
613  return true;
614 
615  return ( area( p1, q1, p2 ) > 0 ) != ( area( p1, q1, q2 ) > 0 )
616  && ( area( p2, q2, p1 ) > 0 ) != ( area( p2, q2, q1 ) > 0 );
617  }
618 
625  bool intersectsPolygon( const Vertex* a, const Vertex* b ) const
626  {
627  const Vertex* p = a->next;
628 
629  do
630  {
631  if( p->i != a->i &&
632  p->next->i != a->i &&
633  p->i != b->i &&
634  p->next->i != b->i && intersects( p, p->next, a, b ) )
635  return true;
636 
637  p = p->next;
638  } while( p != a );
639 
640  return false;
641  }
642 
652  bool locallyInside( const Vertex* a, const Vertex* b ) const
653  {
654  if( area( a->prev, a, a->next ) < 0 )
655  return area( a, b, a->next ) >= 0 && area( a, a->prev, b ) >= 0;
656  else
657  return area( a, b, a->prev ) < 0 || area( a, a->next, b ) < 0;
658  }
659 
666  Vertex* insertVertex( const VECTOR2I& pt, Vertex* last )
667  {
668  m_result.AddVertex( pt );
669  m_vertices.emplace_back( m_result.GetVertexCount() - 1, pt.x, pt.y, this );
670 
671  Vertex* p = &m_vertices.back();
672 
673  if( !last )
674  {
675  p->prev = p;
676  p->next = p;
677  }
678  else
679  {
680  p->next = last->next;
681  p->prev = last;
682  last->next->prev = p;
683  last->next = p;
684  }
685  return p;
686  }
687 
688 private:
690  std::deque<Vertex> m_vertices;
692 };
693 
694 #endif //__POLYGON_TRIANGULATION_H
CITER next(CITER it)
Definition: ptree.cpp:126
bool TesselatePolygon(const SHAPE_LINE_CHAIN &aPoly)
Vertex * createList(const ClipperLib::Path &aPath)
Take a Clipper path and converts it into a circular, doubly-linked list for triangulation.
coord_type GetX() const
Definition: box2.h:173
bool operator!=(const Vertex &rhs) const
Vertex * split(Vertex *b)
Split the referenced polygon between the reference point and vertex b, assuming they are in the same ...
bool intersects(const Vertex *p1, const Vertex *q1, const Vertex *p2, const Vertex *q2) const
Check for intersection between two segments, end points included.
void remove()
Remove the node from the linked list and z-ordered linked list.
bool intersectsPolygon(const Vertex *a, const Vertex *b) const
Check whether the segment from vertex a -> vertex b crosses any of the segments of the polygon of whi...
VECTOR2< int > VECTOR2I
Definition: vector2d.h:623
int PointCount() const
Return the number of points (vertices) in this line chain.
PolygonTriangulation(SHAPE_POLY_SET::TRIANGULATED_POLYGON &aResult)
Vertex * createList(const SHAPE_LINE_CHAIN &points)
Take a SHAPE_LINE_CHAIN and links each point into a circular, doubly-linked list.
bool earcutList(Vertex *aPoint, int pass=0)
Walk through a circular linked list starting at aPoint.
const VECTOR2I & CPoint(int aIndex) const
Return a reference to a given point in the line chain.
Vertex * removeNullTriangles(Vertex *aStart)
Iterate through the list to remove NULL triangles if they exist.
int32_t zOrder(const double aX, const double aY) const
Calculate the Morton code of the Vertex http://www.graphics.stanford.edu/~seander/bithacks....
bool goodSplit(const Vertex *a, const Vertex *b) const
Check if a segment joining two vertices lies fully inside the polygon.
bool operator==(const Vertex &rhs) const
const BOX2I BBox(int aClearance=0) const override
Compute a bounding box of the shape, with a margin of aClearance a collision.
void zSort()
Sort all vertices in this vertex's list by their Morton code.
coord_type GetWidth() const
Definition: box2.h:180
SHAPE_POLY_SET::TRIANGULATED_POLYGON & m_result
E_SERIE r
Definition: eserie.cpp:41
coord_type GetY() const
Definition: box2.h:174
void updateList()
After inserting or changing nodes, this function should be called to remove duplicate vertices and en...
bool locallyInside(const Vertex *a, const Vertex *b) const
Check whether the segment from vertex a -> vertex b is inside the polygon around the immediate area o...
void AddVertex(const VECTOR2I &aP)
void AddTriangle(int a, int b, int c)
Represent a polyline (an zero-thickness chain of connected line segments).
coord_type GetHeight() const
Definition: box2.h:181
Vertex(size_t aIndex, double aX, double aY, PolygonTriangulation *aParent)
bool inTriangle(const Vertex &a, const Vertex &b, const Vertex &c)
Check to see if triangle surrounds our current vertex.
Vertex & operator=(const Vertex &)=delete
bool isEar(Vertex *aEar) const
Check whether the given vertex is in the middle of an ear.
double area(const Vertex *p, const Vertex *q, const Vertex *r) const
Return the twice the signed area of the triangle formed by vertices p, q, and r.
std::deque< Vertex > m_vertices
Vertex * insertVertex(const VECTOR2I &pt, Vertex *last)
Create an entry in the vertices lookup and optionally inserts the newly created vertex into an existi...
void splitPolygon(Vertex *start)
If we cannot find an ear to slice in the current polygon list, we use this to split the polygon into ...