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 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 
266  bool removeNullTriangles( Vertex* aStart )
267  {
268  bool retval = false;
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 = true;
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  p->remove();
290  retval = true;
291  }
292 
293  return retval;
294  }
295 
301  Vertex* createList( const ClipperLib::Path& aPath )
302  {
303  Vertex* tail = nullptr;
304 
305  for( auto point : aPath )
306  tail = insertVertex( VECTOR2I( point.X, point.Y ), tail );
307 
308  if( tail && ( *tail == *tail->next ) )
309  {
310  tail->next->remove();
311  }
312 
313  return tail;
314 
315  }
316 
323  {
324  Vertex* tail = nullptr;
325 
326  for( int i = 0; i < points.PointCount(); i++ )
327  tail = insertVertex( points.CPoint( i ), tail );
328 
329  if( tail && ( *tail == *tail->next ) )
330  {
331  tail->next->remove();
332  }
333 
334  return tail;
335  }
336 
348  bool earcutList( Vertex* aPoint, int pass = 0 )
349  {
350  if( !aPoint )
351  return true;
352 
353  Vertex* stop = aPoint;
354  Vertex* prev;
355  Vertex* next;
356 
357  while( aPoint->prev != aPoint->next )
358  {
359  prev = aPoint->prev;
360  next = aPoint->next;
361 
362  if( isEar( aPoint ) )
363  {
364  m_result.AddTriangle( prev->i, aPoint->i, next->i );
365  aPoint->remove();
366 
367  // Skip one vertex as the triangle will account for the prev node
368  aPoint = next->next;
369  stop = next->next;
370 
371  continue;
372  }
373 
374  Vertex* nextNext = next->next;
375 
376  if( *prev != *nextNext && intersects( prev, aPoint, next, nextNext ) &&
377  locallyInside( prev, nextNext ) &&
378  locallyInside( nextNext, prev ) )
379  {
380  m_result.AddTriangle( prev->i, aPoint->i, nextNext->i );
381 
382  // remove two nodes involved
383  next->remove();
384  aPoint->remove();
385 
386  aPoint = nextNext;
387 
388  continue;
389  }
390 
391  aPoint = next;
392 
397  if( aPoint == stop )
398  {
399  // First, try to remove the remaining steiner points
400  if( removeNullTriangles( aPoint ) )
401  continue;
402 
403  // If we don't have any NULL triangles left, cut the polygon in two and try again
404  splitPolygon( aPoint );
405  break;
406  }
407  }
408 
412  return( aPoint->prev == aPoint->next );
413  }
414 
424  bool isEar( Vertex* aEar ) const
425  {
426  const Vertex* a = aEar->prev;
427  const Vertex* b = aEar;
428  const Vertex* c = aEar->next;
429 
430  // If the area >=0, then the three points for a concave sequence
431  // with b as the reflex point
432  if( area( a, b, c ) >= 0 )
433  return false;
434 
435  // triangle bbox
436  const double minTX = std::min( a->x, std::min( b->x, c->x ) );
437  const double minTY = std::min( a->y, std::min( b->y, c->y ) );
438  const double maxTX = std::max( a->x, std::max( b->x, c->x ) );
439  const double maxTY = std::max( a->y, std::max( b->y, c->y ) );
440 
441  // z-order range for the current triangle bounding box
442  const int32_t minZ = zOrder( minTX, minTY );
443  const int32_t maxZ = zOrder( maxTX, maxTY );
444 
445  // first look for points inside the triangle in increasing z-order
446  Vertex* p = aEar->nextZ;
447 
448  while( p && p->z <= maxZ )
449  {
450  if( p != a && p != c
451  && p->inTriangle( *a, *b, *c )
452  && area( p->prev, p, p->next ) >= 0 )
453  return false;
454  p = p->nextZ;
455  }
456 
457  // then look for points in decreasing z-order
458  p = aEar->prevZ;
459 
460  while( p && p->z >= minZ )
461  {
462  if( p != a && p != c
463  && p->inTriangle( *a, *b, *c )
464  && area( p->prev, p, p->next ) >= 0 )
465  return false;
466  p = p->prevZ;
467  }
468 
469  return true;
470  }
471 
479  void splitPolygon( Vertex* start )
480  {
481  Vertex* origPoly = start;
482  do
483  {
484  Vertex* marker = origPoly->next->next;
485  while( marker != origPoly->prev )
486  {
487  // Find a diagonal line that is wholly enclosed by the polygon interior
488  if( origPoly->i != marker->i && goodSplit( origPoly, marker ) )
489  {
490  Vertex* newPoly = origPoly->split( marker );
491 
492  origPoly->updateList();
493  newPoly->updateList();
494 
495  earcutList( origPoly );
496  earcutList( newPoly );
497  return;
498  }
499  marker = marker->next;
500  }
501  origPoly = origPoly->next;
502  } while( origPoly != start );
503  }
504 
513  bool goodSplit( const Vertex* a, const Vertex* b ) const
514  {
515  return a->next->i != b->i &&
516  a->prev->i != b->i &&
517  !intersectsPolygon( a, b ) &&
518  locallyInside( a, b );
519  }
520 
526  double area( const Vertex* p, const Vertex* q, const Vertex* r ) const
527  {
528  return ( q->y - p->y ) * ( r->x - q->x ) - ( q->x - p->x ) * ( r->y - q->y );
529  }
530 
536  bool intersects( const Vertex* p1, const Vertex* q1, const Vertex* p2, const Vertex* q2 ) const
537  {
538  if( ( *p1 == *q1 && *p2 == *q2 ) || ( *p1 == *q2 && *p2 == *q1 ) )
539  return true;
540 
541  return ( area( p1, q1, p2 ) > 0 ) != ( area( p1, q1, q2 ) > 0 )
542  && ( area( p2, q2, p1 ) > 0 ) != ( area( p2, q2, q1 ) > 0 );
543  }
544 
551  bool intersectsPolygon( const Vertex* a, const Vertex* b ) const
552  {
553  const Vertex* p = a->next;
554  do
555  {
556  if( p->i != a->i &&
557  p->next->i != a->i &&
558  p->i != b->i &&
559  p->next->i != b->i && intersects( p, p->next, a, b ) )
560  return true;
561 
562  p = p->next;
563  } while( p != a );
564 
565  return false;
566  }
567 
576  bool locallyInside( const Vertex* a, const Vertex* b ) const
577  {
578  if( area( a->prev, a, a->next ) < 0 )
579  return area( a, b, a->next ) >= 0 && area( a, a->prev, b ) >= 0;
580  else
581  return area( a, b, a->prev ) < 0 || area( a, a->next, b ) < 0;
582  }
583 
590  Vertex* insertVertex( const VECTOR2I& pt, Vertex* last )
591  {
592  m_result.AddVertex( pt );
593  m_vertices.emplace_back( m_result.GetVertexCount() - 1, pt.x, pt.y, this );
594 
595  Vertex* p = &m_vertices.back();
596  if( !last )
597  {
598  p->prev = p;
599  p->next = p;
600  }
601  else
602  {
603  p->next = last->next;
604  p->prev = last;
605  last->next->prev = p;
606  last->next = p;
607  }
608  return p;
609  }
610 
611 
612 public:
613 
614  void TesselatePolygon( const SHAPE_LINE_CHAIN& aPoly )
615  {
616  ClipperLib::Clipper c;
617  m_bbox = aPoly.BBox();
618 
619  if( !m_bbox.GetWidth() || !m_bbox.GetHeight() )
620  return;
621 
622  Vertex* outerNode = createList( aPoly );
623  if( !outerNode )
624  return;
625 
626  outerNode->updateList();
627  if( !earcutList( outerNode ) )
628  {
629  m_vertices.clear();
630  m_result.Clear();
631 
632  ClipperLib::Paths simplified;
633  ClipperLib::SimplifyPolygon( aPoly.convertToClipper( true ), simplified );
634 
635  for( auto path : simplified )
636  {
637  outerNode = createList( path );
638  if( !outerNode )
639  return;
640 
641  earcutList( outerNode );
642  }
643  }
644 
645  m_vertices.clear();
646  }
647 };
648 
649 #endif //__POLYGON_TRIANGULATION_H
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.
coord_type GetY() const
Definition: box2.h:189
Vertex * createList(const ClipperLib::Path &aPath)
Function createList Takes a Clipper path and converts it into a circular, doubly-linked list for tria...
bool operator!=(const Vertex &rhs) const
ClipperLib::Path convertToClipper(bool aRequiredOrientation) const
Creates a new Clipper path from the SHAPE_LINE_CHAIN in a given orientation.
int PointCount() const
Function PointCount()
Vertex * split(Vertex *b)
Function split Splits the referenced polygon between the reference point and vertex b...
void remove()
Function remove Removes the node from the linked list and z-ordered linked list.
bool goodSplit(const Vertex *a, const Vertex *b) const
Check if a segment joining two vertices lies fully inside the polygon.
void TesselatePolygon(const SHAPE_LINE_CHAIN &aPoly)
VECTOR2< int > VECTOR2I
Definition: vector2d.h:587
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...
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, doubly-linked list.
coord_type GetWidth() const
Definition: box2.h:195
bool earcutList(Vertex *aPoint, int pass=0)
Function: earcutList Walks through a circular linked list starting at aPoint.
bool removeNullTriangles(Vertex *aStart)
Function removeNullTriangles Iterates through the list to remove NULL triangles if they exist...
const BOX2I BBox(int aClearance=0) const override
Function BBox()
void zSort()
Sort all vertices in this vertex&#39;s list by their Morton code.
SHAPE_POLY_SET::TRIANGULATED_POLYGON & m_result
coord_type GetHeight() const
Definition: box2.h:196
void updateList()
Function updateList After inserting or changing nodes, this function should be called to remove dupli...
void AddVertex(const VECTOR2I &aP)
#define max(a, b)
Definition: auxiliary.h:86
Class SHAPE_LINE_CHAIN.
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 operator==(const Vertex &rhs) const
std::deque< Vertex > m_vertices
int32_t zOrder(const double aX, const double aY) const
Calculate the Morton code of the Vertex http://www.graphics.stanford.edu/~seander/bithacks.html#InterleaveBMN.
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...
bool isEar(Vertex *aEar) const
Function isEar Checks whether the given vertex is in the middle of an ear.
coord_type GetX() const
Definition: box2.h:188
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...
const VECTOR2I & CPoint(int aIndex) const
Function CPoint()
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, we use this to split the polygon into two separate lists and slice them each independently.