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-2019 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 <cmath>
51 #include <vector>
52 #include <math/box2.h>
53 
54 #include "clipper.hpp"
55 
57 {
58 
59 public:
60 
62  m_result( aResult )
63  {};
64 
65 private:
66  struct Vertex
67  {
68  Vertex( size_t aIndex, double aX, double aY, PolygonTriangulation* aParent ) :
69  i( aIndex ), x( aX ), y( aY ), parent( aParent )
70  {
71  }
72  Vertex& operator=( const Vertex& ) = delete;
73  Vertex& operator=( Vertex&& ) = delete;
74 
75  bool operator==( const Vertex& rhs ) const
76  {
77  return this->x == rhs.x && this->y == rhs.y;
78  }
79  bool operator!=( const Vertex& rhs ) const { return !( *this == rhs ); }
80 
81 
94  {
95  parent->m_vertices.emplace_back( i, x, y, parent );
96  Vertex* a2 = &parent->m_vertices.back();
97  parent->m_vertices.emplace_back( b->i, b->x, b->y, parent );
98  Vertex* b2 = &parent->m_vertices.back();
99  Vertex* an = next;
100  Vertex* bp = b->prev;
101 
102  next = b;
103  b->prev = this;
104 
105  a2->next = an;
106  an->prev = a2;
107 
108  b2->next = a2;
109  a2->prev = b2;
110 
111  bp->next = b2;
112  b2->prev = bp;
113 
114  return b2;
115  }
116 
121  void remove()
122  {
123  next->prev = prev;
124  prev->next = next;
125 
126  if( prevZ )
127  prevZ->nextZ = nextZ;
128  if( nextZ )
129  nextZ->prevZ = prevZ;
130  next = NULL;
131  prev = NULL;
132  nextZ = NULL;
133  prevZ = NULL;
134  }
135 
136 
137  void updateOrder()
138  {
139  if( !z )
140  z = parent->zOrder( x, y );
141  }
142 
148  void updateList()
149  {
150  Vertex* p = next;
151 
152  while( p != this )
153  {
157  if( *p == *p->next )
158  {
159  p = p->prev;
160  p->next->remove();
161 
162  if( p == p->next )
163  break;
164  }
165 
166  p->updateOrder();
167  p = p->next;
168  };
169 
170  updateOrder();
171  zSort();
172  }
173 
177  void zSort()
178  {
179  std::deque<Vertex*> queue;
180 
181  queue.push_back( this );
182 
183  for( auto p = next; p && p != this; p = p->next )
184  queue.push_back( p );
185 
186  std::sort( queue.begin(), queue.end(), []( const Vertex* a, const Vertex* b)
187  {
188  return a->z < b->z;
189  } );
190 
191  Vertex* prev_elem = nullptr;
192  for( auto elem : queue )
193  {
194  if( prev_elem )
195  prev_elem->nextZ = elem;
196 
197  elem->prevZ = prev_elem;
198  prev_elem = elem;
199  }
200 
201  prev_elem->nextZ = nullptr;
202  }
203 
204 
208  bool inTriangle( const Vertex& a, const Vertex& b, const Vertex& c )
209  {
210  return ( c.x - x ) * ( a.y - y ) - ( a.x - x ) * ( c.y - y ) >= 0
211  && ( a.x - x ) * ( b.y - y ) - ( b.x - x ) * ( a.y - y ) >= 0
212  && ( b.x - x ) * ( c.y - y ) - ( c.x - x ) * ( b.y - y ) >= 0;
213  }
214 
215  const size_t i;
216  const double x;
217  const double y;
219 
220  // previous and next vertices nodes in a polygon ring
221  Vertex* prev = nullptr;
222  Vertex* next = nullptr;
223 
224  // z-order curve value
225  int32_t z = 0;
226 
227  // previous and next nodes in z-order
228  Vertex* prevZ = nullptr;
229  Vertex* nextZ = nullptr;
230  };
231 
233  std::deque<Vertex> m_vertices;
235 
241  int32_t zOrder( const double aX, const double aY ) const
242  {
243  int32_t x = static_cast<int32_t>( 32767.0 * ( aX - m_bbox.GetX() ) / m_bbox.GetWidth() );
244  int32_t y = static_cast<int32_t>( 32767.0 * ( aY - m_bbox.GetY() ) / m_bbox.GetHeight() );
245 
246  x = ( x | ( x << 8 ) ) & 0x00FF00FF;
247  x = ( x | ( x << 4 ) ) & 0x0F0F0F0F;
248  x = ( x | ( x << 2 ) ) & 0x33333333;
249  x = ( x | ( x << 1 ) ) & 0x55555555;
250 
251  y = ( y | ( y << 8 ) ) & 0x00FF00FF;
252  y = ( y | ( y << 4 ) ) & 0x0F0F0F0F;
253  y = ( y | ( y << 2 ) ) & 0x33333333;
254  y = ( y | ( y << 1 ) ) & 0x55555555;
255 
256  return x | ( y << 1 );
257  }
258 
267  {
268  Vertex* retval = nullptr;
269  Vertex* p = aStart->next;
270 
271  while( p != aStart )
272  {
273  if( area( p->prev, p, p->next ) == 0.0 )
274  {
275  p = p->prev;
276  p->next->remove();
277  retval = aStart;
278 
279  if( p == p->next )
280  break;
281  }
282  p = p->next;
283  };
284 
285  // We needed an end point above that wouldn't be removed, so
286  // here we do the final check for this as a Steiner point
287  if( area( aStart->prev, aStart, aStart->next ) == 0.0 )
288  {
289  retval = p->next;
290  p->remove();
291  }
292 
293  return retval;
294  }
295 
301  Vertex* createList( const ClipperLib::Path& aPath )
302  {
303  Vertex* tail = nullptr;
304  double sum = 0.0;
305  auto len = aPath.size();
306 
307  // Check for winding order
308  for( size_t i = 0; i < len; i++ )
309  {
310  auto p1 = aPath.at( i );
311  auto p2 = aPath.at( ( i + 1 ) < len ? i + 1 : 0 );
312 
313  sum += ( ( p2.X - p1.X ) * ( p2.Y + p1.Y ) );
314  }
315 
316  if( sum <= 0.0 )
317  {
318  for( auto point : aPath )
319  tail = insertVertex( VECTOR2I( point.X, point.Y ), tail );
320  }
321  else
322  {
323  for( size_t i = 0; i < len; i++ )
324  {
325  auto p = aPath.at( len - i - 1 );
326  tail = insertVertex( VECTOR2I( p.X, p.Y ), tail );
327  }
328  }
329 
330  if( tail && ( *tail == *tail->next ) )
331  {
332  tail->next->remove();
333  }
334 
335  return tail;
336 
337  }
338 
345  {
346  Vertex* tail = nullptr;
347  double sum = 0.0;
348 
349  // Check for winding order
350  for( int i = 0; i < points.PointCount(); i++ )
351  {
352  VECTOR2D p1 = points.CPoint( i );
353  VECTOR2D p2 = points.CPoint( i + 1 );
354 
355  sum += ( ( p2.x - p1.x ) * ( p2.y + p1.y ) );
356  }
357 
358  if( sum > 0.0 )
359  for( int i = points.PointCount() - 1; i >= 0; i--)
360  tail = insertVertex( points.CPoint( i ), tail );
361  else
362  for( int i = 0; i < points.PointCount(); i++ )
363  tail = insertVertex( points.CPoint( i ), tail );
364 
365  if( tail && ( *tail == *tail->next ) )
366  {
367  tail->next->remove();
368  }
369 
370  return tail;
371  }
372 
384  bool earcutList( Vertex* aPoint, int pass = 0 )
385  {
386  if( !aPoint )
387  return true;
388 
389  Vertex* stop = aPoint;
390  Vertex* prev;
391  Vertex* next;
392 
393  while( aPoint->prev != aPoint->next )
394  {
395  prev = aPoint->prev;
396  next = aPoint->next;
397 
398  if( isEar( aPoint ) )
399  {
400  m_result.AddTriangle( prev->i, aPoint->i, next->i );
401  aPoint->remove();
402 
403  // Skip one vertex as the triangle will account for the prev node
404  aPoint = next->next;
405  stop = next->next;
406 
407  continue;
408  }
409 
410  Vertex* nextNext = next->next;
411 
412  if( *prev != *nextNext && intersects( prev, aPoint, next, nextNext ) &&
413  locallyInside( prev, nextNext ) &&
414  locallyInside( nextNext, prev ) )
415  {
416  m_result.AddTriangle( prev->i, aPoint->i, nextNext->i );
417 
418  // remove two nodes involved
419  next->remove();
420  aPoint->remove();
421 
422  aPoint = nextNext;
423  stop = nextNext;
424 
425  continue;
426  }
427 
428  aPoint = next;
429 
434  if( aPoint == stop )
435  {
436  // First, try to remove the remaining steiner points
437  // If aPoint is a steiner, we need to re-assign both the start and stop points
438  if( auto newPoint = removeNullTriangles( aPoint ) )
439  {
440  aPoint = newPoint;
441  stop = newPoint;
442  continue;
443  }
444 
445  // If we don't have any NULL triangles left, cut the polygon in two and try again
446  splitPolygon( aPoint );
447  break;
448  }
449  }
450 
454  return( aPoint->prev == aPoint->next );
455  }
456 
466  bool isEar( Vertex* aEar ) const
467  {
468  const Vertex* a = aEar->prev;
469  const Vertex* b = aEar;
470  const Vertex* c = aEar->next;
471 
472  // If the area >=0, then the three points for a concave sequence
473  // with b as the reflex point
474  if( area( a, b, c ) >= 0 )
475  return false;
476 
477  // triangle bbox
478  const double minTX = std::min( a->x, std::min( b->x, c->x ) );
479  const double minTY = std::min( a->y, std::min( b->y, c->y ) );
480  const double maxTX = std::max( a->x, std::max( b->x, c->x ) );
481  const double maxTY = std::max( a->y, std::max( b->y, c->y ) );
482 
483  // z-order range for the current triangle bounding box
484  const int32_t minZ = zOrder( minTX, minTY );
485  const int32_t maxZ = zOrder( maxTX, maxTY );
486 
487  // first look for points inside the triangle in increasing z-order
488  Vertex* p = aEar->nextZ;
489 
490  while( p && p->z <= maxZ )
491  {
492  if( p != a && p != c
493  && p->inTriangle( *a, *b, *c )
494  && area( p->prev, p, p->next ) >= 0 )
495  return false;
496  p = p->nextZ;
497  }
498 
499  // then look for points in decreasing z-order
500  p = aEar->prevZ;
501 
502  while( p && p->z >= minZ )
503  {
504  if( p != a && p != c
505  && p->inTriangle( *a, *b, *c )
506  && area( p->prev, p, p->next ) >= 0 )
507  return false;
508  p = p->prevZ;
509  }
510 
511  return true;
512  }
513 
521  void splitPolygon( Vertex* start )
522  {
523  Vertex* origPoly = start;
524  do
525  {
526  Vertex* marker = origPoly->next->next;
527  while( marker != origPoly->prev )
528  {
529  // Find a diagonal line that is wholly enclosed by the polygon interior
530  if( origPoly->i != marker->i && goodSplit( origPoly, marker ) )
531  {
532  Vertex* newPoly = origPoly->split( marker );
533 
534  origPoly->updateList();
535  newPoly->updateList();
536 
537  earcutList( origPoly );
538  earcutList( newPoly );
539  return;
540  }
541  marker = marker->next;
542  }
543  origPoly = origPoly->next;
544  } while( origPoly != start );
545  }
546 
555  bool goodSplit( const Vertex* a, const Vertex* b ) const
556  {
557  return a->next->i != b->i &&
558  a->prev->i != b->i &&
559  !intersectsPolygon( a, b ) &&
560  locallyInside( a, b );
561  }
562 
568  double area( const Vertex* p, const Vertex* q, const Vertex* r ) const
569  {
570  return ( q->y - p->y ) * ( r->x - q->x ) - ( q->x - p->x ) * ( r->y - q->y );
571  }
572 
578  bool intersects( const Vertex* p1, const Vertex* q1, const Vertex* p2, const Vertex* q2 ) const
579  {
580  if( ( *p1 == *q1 && *p2 == *q2 ) || ( *p1 == *q2 && *p2 == *q1 ) )
581  return true;
582 
583  return ( area( p1, q1, p2 ) > 0 ) != ( area( p1, q1, q2 ) > 0 )
584  && ( area( p2, q2, p1 ) > 0 ) != ( area( p2, q2, q1 ) > 0 );
585  }
586 
593  bool intersectsPolygon( const Vertex* a, const Vertex* b ) const
594  {
595  const Vertex* p = a->next;
596  do
597  {
598  if( p->i != a->i &&
599  p->next->i != a->i &&
600  p->i != b->i &&
601  p->next->i != b->i && intersects( p, p->next, a, b ) )
602  return true;
603 
604  p = p->next;
605  } while( p != a );
606 
607  return false;
608  }
609 
618  bool locallyInside( const Vertex* a, const Vertex* b ) const
619  {
620  if( area( a->prev, a, a->next ) < 0 )
621  return area( a, b, a->next ) >= 0 && area( a, a->prev, b ) >= 0;
622  else
623  return area( a, b, a->prev ) < 0 || area( a, a->next, b ) < 0;
624  }
625 
632  Vertex* insertVertex( const VECTOR2I& pt, Vertex* last )
633  {
634  m_result.AddVertex( pt );
635  m_vertices.emplace_back( m_result.GetVertexCount() - 1, pt.x, pt.y, this );
636 
637  Vertex* p = &m_vertices.back();
638  if( !last )
639  {
640  p->prev = p;
641  p->next = p;
642  }
643  else
644  {
645  p->next = last->next;
646  p->prev = last;
647  last->next->prev = p;
648  last->next = p;
649  }
650  return p;
651  }
652 
653 
654 public:
655 
656  bool TesselatePolygon( const SHAPE_LINE_CHAIN& aPoly )
657  {
658  m_bbox = aPoly.BBox();
659  m_result.Clear();
660 
661  if( !m_bbox.GetWidth() || !m_bbox.GetHeight() )
662  return false;
663 
667  Vertex* firstVertex = createList( aPoly );
668  if( !firstVertex || firstVertex->prev == firstVertex->next )
669  return false;
670 
671  firstVertex->updateList();
672 
673  auto retval = earcutList( firstVertex );
674  m_vertices.clear();
675  return retval;
676  }
677 };
678 
679 #endif //__POLYGON_TRIANGULATION_H
CITER next(CITER it)
Definition: ptree.cpp:130
bool TesselatePolygon(const SHAPE_LINE_CHAIN &aPoly)
Vertex * createList(const ClipperLib::Path &aPath)
Function createList Takes a Clipper path and converts it into a circular, doubly-linked list for tria...
coord_type GetX() const
Definition: box2.h:188
bool operator!=(const Vertex &rhs) const
Vertex * split(Vertex *b)
Function split Splits the referenced polygon between the reference point and vertex b,...
bool intersects(const Vertex *p1, const Vertex *q1, const Vertex *p2, const Vertex *q2) const
Function intersects Checks for intersection between two segments, end points included.
void remove()
Function remove Removes the node from the linked list and z-ordered linked list.
bool intersectsPolygon(const Vertex *a, const Vertex *b) const
Function intersectsPolygon Checks whether the segment from vertex a -> vertex b crosses any of the se...
VECTOR2< int > VECTOR2I
Definition: vector2d.h:587
int PointCount() const
Function PointCount()
PolygonTriangulation(SHAPE_POLY_SET::TRIANGULATED_POLYGON &aResult)
Vertex * createList(const SHAPE_LINE_CHAIN &points)
Function createList Takes the SHAPE_LINE_CHAIN and links each point into a circular,...
bool earcutList(Vertex *aPoint, int pass=0)
Function: earcutList Walks through a circular linked list starting at aPoint.
const VECTOR2I & CPoint(int aIndex) const
Function CPoint()
Vertex * removeNullTriangles(Vertex *aStart)
Function removeNullTriangles Iterates 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
Function BBox()
void zSort()
Sort all vertices in this vertex's list by their Morton code.
coord_type GetWidth() const
Definition: box2.h:195
SHAPE_POLY_SET::TRIANGULATED_POLYGON & m_result
coord_type GetY() const
Definition: box2.h:189
void updateList()
Function updateList After inserting or changing nodes, this function should be called to remove dupli...
bool locallyInside(const Vertex *a, const Vertex *b) const
Function locallyInside Checks whether the segment from vertex a -> vertex b is inside the polygon aro...
void AddVertex(const VECTOR2I &aP)
#define max(a, b)
Definition: auxiliary.h:86
Class SHAPE_LINE_CHAIN.
size_t i
Definition: json11.cpp:597
coord_type GetHeight() const
Definition: box2.h:196
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
Function isEar Checks whether the given vertex is in the middle of an ear.
double area(const Vertex *p, const Vertex *q, const Vertex *r) const
Function area Returns the twice the signed area of the triangle formed by vertices p,...
std::deque< Vertex > m_vertices
Vertex * insertVertex(const VECTOR2I &pt, Vertex *last)
Function insertVertex Creates an entry in the vertices lookup and optionally inserts the newly create...
#define min(a, b)
Definition: auxiliary.h:85
void splitPolygon(Vertex *start)
Function splitPolygon If we cannot find an ear to slice in the current polygon list,...