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  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 
424  continue;
425  }
426 
427  aPoint = next;
428 
433  if( aPoint == stop )
434  {
435  // First, try to remove the remaining steiner points
436  if( removeNullTriangles( aPoint ) )
437  continue;
438 
439  // If we don't have any NULL triangles left, cut the polygon in two and try again
440  splitPolygon( aPoint );
441  break;
442  }
443  }
444 
448  return( aPoint->prev == aPoint->next );
449  }
450 
460  bool isEar( Vertex* aEar ) const
461  {
462  const Vertex* a = aEar->prev;
463  const Vertex* b = aEar;
464  const Vertex* c = aEar->next;
465 
466  // If the area >=0, then the three points for a concave sequence
467  // with b as the reflex point
468  if( area( a, b, c ) >= 0 )
469  return false;
470 
471  // triangle bbox
472  const double minTX = std::min( a->x, std::min( b->x, c->x ) );
473  const double minTY = std::min( a->y, std::min( b->y, c->y ) );
474  const double maxTX = std::max( a->x, std::max( b->x, c->x ) );
475  const double maxTY = std::max( a->y, std::max( b->y, c->y ) );
476 
477  // z-order range for the current triangle bounding box
478  const int32_t minZ = zOrder( minTX, minTY );
479  const int32_t maxZ = zOrder( maxTX, maxTY );
480 
481  // first look for points inside the triangle in increasing z-order
482  Vertex* p = aEar->nextZ;
483 
484  while( p && p->z <= maxZ )
485  {
486  if( p != a && p != c
487  && p->inTriangle( *a, *b, *c )
488  && area( p->prev, p, p->next ) >= 0 )
489  return false;
490  p = p->nextZ;
491  }
492 
493  // then look for points in decreasing z-order
494  p = aEar->prevZ;
495 
496  while( p && p->z >= minZ )
497  {
498  if( p != a && p != c
499  && p->inTriangle( *a, *b, *c )
500  && area( p->prev, p, p->next ) >= 0 )
501  return false;
502  p = p->prevZ;
503  }
504 
505  return true;
506  }
507 
515  void splitPolygon( Vertex* start )
516  {
517  Vertex* origPoly = start;
518  do
519  {
520  Vertex* marker = origPoly->next->next;
521  while( marker != origPoly->prev )
522  {
523  // Find a diagonal line that is wholly enclosed by the polygon interior
524  if( origPoly->i != marker->i && goodSplit( origPoly, marker ) )
525  {
526  Vertex* newPoly = origPoly->split( marker );
527 
528  origPoly->updateList();
529  newPoly->updateList();
530 
531  earcutList( origPoly );
532  earcutList( newPoly );
533  return;
534  }
535  marker = marker->next;
536  }
537  origPoly = origPoly->next;
538  } while( origPoly != start );
539  }
540 
549  bool goodSplit( const Vertex* a, const Vertex* b ) const
550  {
551  return a->next->i != b->i &&
552  a->prev->i != b->i &&
553  !intersectsPolygon( a, b ) &&
554  locallyInside( a, b );
555  }
556 
562  double area( const Vertex* p, const Vertex* q, const Vertex* r ) const
563  {
564  return ( q->y - p->y ) * ( r->x - q->x ) - ( q->x - p->x ) * ( r->y - q->y );
565  }
566 
572  bool intersects( const Vertex* p1, const Vertex* q1, const Vertex* p2, const Vertex* q2 ) const
573  {
574  if( ( *p1 == *q1 && *p2 == *q2 ) || ( *p1 == *q2 && *p2 == *q1 ) )
575  return true;
576 
577  return ( area( p1, q1, p2 ) > 0 ) != ( area( p1, q1, q2 ) > 0 )
578  && ( area( p2, q2, p1 ) > 0 ) != ( area( p2, q2, q1 ) > 0 );
579  }
580 
587  bool intersectsPolygon( const Vertex* a, const Vertex* b ) const
588  {
589  const Vertex* p = a->next;
590  do
591  {
592  if( p->i != a->i &&
593  p->next->i != a->i &&
594  p->i != b->i &&
595  p->next->i != b->i && intersects( p, p->next, a, b ) )
596  return true;
597 
598  p = p->next;
599  } while( p != a );
600 
601  return false;
602  }
603 
612  bool locallyInside( const Vertex* a, const Vertex* b ) const
613  {
614  if( area( a->prev, a, a->next ) < 0 )
615  return area( a, b, a->next ) >= 0 && area( a, a->prev, b ) >= 0;
616  else
617  return area( a, b, a->prev ) < 0 || area( a, a->next, b ) < 0;
618  }
619 
626  Vertex* insertVertex( const VECTOR2I& pt, Vertex* last )
627  {
628  m_result.AddVertex( pt );
629  m_vertices.emplace_back( m_result.GetVertexCount() - 1, pt.x, pt.y, this );
630 
631  Vertex* p = &m_vertices.back();
632  if( !last )
633  {
634  p->prev = p;
635  p->next = p;
636  }
637  else
638  {
639  p->next = last->next;
640  p->prev = last;
641  last->next->prev = p;
642  last->next = p;
643  }
644  return p;
645  }
646 
647 
648 public:
649 
650  bool TesselatePolygon( const SHAPE_LINE_CHAIN& aPoly )
651  {
652  m_bbox = aPoly.BBox();
653  m_result.Clear();
654 
655  if( !m_bbox.GetWidth() || !m_bbox.GetHeight() )
656  return false;
657 
661  Vertex* firstVertex = createList( aPoly );
662  if( !firstVertex || firstVertex->prev == firstVertex->next )
663  return false;
664 
665  firstVertex->updateList();
666 
667  auto retval = earcutList( firstVertex );
668  m_vertices.clear();
669  return retval;
670  }
671 };
672 
673 #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
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...
bool operator!=(const Vertex &rhs) const
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.
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.