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 <deque>
51 #include <cmath>
52 
53 #include <clipper.hpp>
56 #include <math/box2.h>
57 #include <math/vector2d.h>
58 
60 {
61 
62 public:
63 
65  m_result( aResult )
66  {};
67 
68 private:
69  struct Vertex
70  {
71  Vertex( size_t aIndex, double aX, double aY, PolygonTriangulation* aParent ) :
72  i( aIndex ), x( aX ), y( aY ), parent( aParent )
73  {
74  }
75  Vertex& operator=( const Vertex& ) = delete;
76  Vertex& operator=( Vertex&& ) = delete;
77 
78  bool operator==( const Vertex& rhs ) const
79  {
80  return this->x == rhs.x && this->y == rhs.y;
81  }
82  bool operator!=( const Vertex& rhs ) const { return !( *this == rhs ); }
83 
84 
97  {
98  parent->m_vertices.emplace_back( i, x, y, parent );
99  Vertex* a2 = &parent->m_vertices.back();
100  parent->m_vertices.emplace_back( b->i, b->x, b->y, parent );
101  Vertex* b2 = &parent->m_vertices.back();
102  Vertex* an = next;
103  Vertex* bp = b->prev;
104 
105  next = b;
106  b->prev = this;
107 
108  a2->next = an;
109  an->prev = a2;
110 
111  b2->next = a2;
112  a2->prev = b2;
113 
114  bp->next = b2;
115  b2->prev = bp;
116 
117  return b2;
118  }
119 
124  void remove()
125  {
126  next->prev = prev;
127  prev->next = next;
128 
129  if( prevZ )
130  prevZ->nextZ = nextZ;
131  if( nextZ )
132  nextZ->prevZ = prevZ;
133  next = NULL;
134  prev = NULL;
135  nextZ = NULL;
136  prevZ = NULL;
137  }
138 
139 
140  void updateOrder()
141  {
142  if( !z )
143  z = parent->zOrder( x, y );
144  }
145 
151  void updateList()
152  {
153  Vertex* p = next;
154 
155  while( p != this )
156  {
160  if( *p == *p->next )
161  {
162  p = p->prev;
163  p->next->remove();
164 
165  if( p == p->next )
166  break;
167  }
168 
169  p->updateOrder();
170  p = p->next;
171  };
172 
173  updateOrder();
174  zSort();
175  }
176 
180  void zSort()
181  {
182  std::deque<Vertex*> queue;
183 
184  queue.push_back( this );
185 
186  for( auto p = next; p && p != this; p = p->next )
187  queue.push_back( p );
188 
189  std::sort( queue.begin(), queue.end(), []( const Vertex* a, const Vertex* b)
190  {
191  return a->z < b->z;
192  } );
193 
194  Vertex* prev_elem = nullptr;
195  for( auto elem : queue )
196  {
197  if( prev_elem )
198  prev_elem->nextZ = elem;
199 
200  elem->prevZ = prev_elem;
201  prev_elem = elem;
202  }
203 
204  prev_elem->nextZ = nullptr;
205  }
206 
207 
211  bool inTriangle( const Vertex& a, const Vertex& b, const Vertex& c )
212  {
213  return ( c.x - x ) * ( a.y - y ) - ( a.x - x ) * ( c.y - y ) >= 0
214  && ( a.x - x ) * ( b.y - y ) - ( b.x - x ) * ( a.y - y ) >= 0
215  && ( b.x - x ) * ( c.y - y ) - ( c.x - x ) * ( b.y - y ) >= 0;
216  }
217 
218  const size_t i;
219  const double x;
220  const double y;
222 
223  // previous and next vertices nodes in a polygon ring
224  Vertex* prev = nullptr;
225  Vertex* next = nullptr;
226 
227  // z-order curve value
228  int32_t z = 0;
229 
230  // previous and next nodes in z-order
231  Vertex* prevZ = nullptr;
232  Vertex* nextZ = nullptr;
233  };
234 
236  std::deque<Vertex> m_vertices;
238 
244  int32_t zOrder( const double aX, const double aY ) const
245  {
246  int32_t x = static_cast<int32_t>( 32767.0 * ( aX - m_bbox.GetX() ) / m_bbox.GetWidth() );
247  int32_t y = static_cast<int32_t>( 32767.0 * ( aY - m_bbox.GetY() ) / m_bbox.GetHeight() );
248 
249  x = ( x | ( x << 8 ) ) & 0x00FF00FF;
250  x = ( x | ( x << 4 ) ) & 0x0F0F0F0F;
251  x = ( x | ( x << 2 ) ) & 0x33333333;
252  x = ( x | ( x << 1 ) ) & 0x55555555;
253 
254  y = ( y | ( y << 8 ) ) & 0x00FF00FF;
255  y = ( y | ( y << 4 ) ) & 0x0F0F0F0F;
256  y = ( y | ( y << 2 ) ) & 0x33333333;
257  y = ( y | ( y << 1 ) ) & 0x55555555;
258 
259  return x | ( y << 1 );
260  }
261 
270  {
271  Vertex* retval = nullptr;
272  Vertex* p = aStart->next;
273 
274  while( p != aStart )
275  {
276  if( area( p->prev, p, p->next ) == 0.0 )
277  {
278  p = p->prev;
279  p->next->remove();
280  retval = aStart;
281 
282  if( p == p->next )
283  break;
284  }
285  p = p->next;
286  };
287 
288  // We needed an end point above that wouldn't be removed, so
289  // here we do the final check for this as a Steiner point
290  if( area( aStart->prev, aStart, aStart->next ) == 0.0 )
291  {
292  retval = p->next;
293  p->remove();
294  }
295 
296  return retval;
297  }
298 
304  Vertex* createList( const ClipperLib::Path& aPath )
305  {
306  Vertex* tail = nullptr;
307  double sum = 0.0;
308  auto len = aPath.size();
309 
310  // Check for winding order
311  for( size_t i = 0; i < len; i++ )
312  {
313  auto p1 = aPath.at( i );
314  auto p2 = aPath.at( ( i + 1 ) < len ? i + 1 : 0 );
315 
316  sum += ( ( p2.X - p1.X ) * ( p2.Y + p1.Y ) );
317  }
318 
319  if( sum <= 0.0 )
320  {
321  for( auto point : aPath )
322  tail = insertVertex( VECTOR2I( point.X, point.Y ), tail );
323  }
324  else
325  {
326  for( size_t i = 0; i < len; i++ )
327  {
328  auto p = aPath.at( len - i - 1 );
329  tail = insertVertex( VECTOR2I( p.X, p.Y ), tail );
330  }
331  }
332 
333  if( tail && ( *tail == *tail->next ) )
334  {
335  tail->next->remove();
336  }
337 
338  return tail;
339 
340  }
341 
348  {
349  Vertex* tail = nullptr;
350  double sum = 0.0;
351 
352  // Check for winding order
353  for( int i = 0; i < points.PointCount(); i++ )
354  {
355  VECTOR2D p1 = points.CPoint( i );
356  VECTOR2D p2 = points.CPoint( i + 1 );
357 
358  sum += ( ( p2.x - p1.x ) * ( p2.y + p1.y ) );
359  }
360 
361  if( sum > 0.0 )
362  for( int i = points.PointCount() - 1; i >= 0; i--)
363  tail = insertVertex( points.CPoint( i ), tail );
364  else
365  for( int i = 0; i < points.PointCount(); i++ )
366  tail = insertVertex( points.CPoint( i ), tail );
367 
368  if( tail && ( *tail == *tail->next ) )
369  {
370  tail->next->remove();
371  }
372 
373  return tail;
374  }
375 
387  bool earcutList( Vertex* aPoint, int pass = 0 )
388  {
389  if( !aPoint )
390  return true;
391 
392  Vertex* stop = aPoint;
393  Vertex* prev;
394  Vertex* next;
395 
396  while( aPoint->prev != aPoint->next )
397  {
398  prev = aPoint->prev;
399  next = aPoint->next;
400 
401  if( isEar( aPoint ) )
402  {
403  m_result.AddTriangle( prev->i, aPoint->i, next->i );
404  aPoint->remove();
405 
406  // Skip one vertex as the triangle will account for the prev node
407  aPoint = next->next;
408  stop = next->next;
409 
410  continue;
411  }
412 
413  Vertex* nextNext = next->next;
414 
415  if( *prev != *nextNext && intersects( prev, aPoint, next, nextNext ) &&
416  locallyInside( prev, nextNext ) &&
417  locallyInside( nextNext, prev ) )
418  {
419  m_result.AddTriangle( prev->i, aPoint->i, nextNext->i );
420 
421  // remove two nodes involved
422  next->remove();
423  aPoint->remove();
424 
425  aPoint = nextNext;
426  stop = nextNext;
427 
428  continue;
429  }
430 
431  aPoint = next;
432 
437  if( aPoint == stop )
438  {
439  // First, try to remove the remaining steiner points
440  // If aPoint is a steiner, we need to re-assign both the start and stop points
441  if( auto newPoint = removeNullTriangles( aPoint ) )
442  {
443  aPoint = newPoint;
444  stop = newPoint;
445  continue;
446  }
447 
448  // If we don't have any NULL triangles left, cut the polygon in two and try again
449  splitPolygon( aPoint );
450  break;
451  }
452  }
453 
457  return( aPoint->prev == aPoint->next );
458  }
459 
469  bool isEar( Vertex* aEar ) const
470  {
471  const Vertex* a = aEar->prev;
472  const Vertex* b = aEar;
473  const Vertex* c = aEar->next;
474 
475  // If the area >=0, then the three points for a concave sequence
476  // with b as the reflex point
477  if( area( a, b, c ) >= 0 )
478  return false;
479 
480  // triangle bbox
481  const double minTX = std::min( a->x, std::min( b->x, c->x ) );
482  const double minTY = std::min( a->y, std::min( b->y, c->y ) );
483  const double maxTX = std::max( a->x, std::max( b->x, c->x ) );
484  const double maxTY = std::max( a->y, std::max( b->y, c->y ) );
485 
486  // z-order range for the current triangle bounding box
487  const int32_t minZ = zOrder( minTX, minTY );
488  const int32_t maxZ = zOrder( maxTX, maxTY );
489 
490  // first look for points inside the triangle in increasing z-order
491  Vertex* p = aEar->nextZ;
492 
493  while( p && p->z <= maxZ )
494  {
495  if( p != a && p != c
496  && p->inTriangle( *a, *b, *c )
497  && area( p->prev, p, p->next ) >= 0 )
498  return false;
499  p = p->nextZ;
500  }
501 
502  // then look for points in decreasing z-order
503  p = aEar->prevZ;
504 
505  while( p && p->z >= minZ )
506  {
507  if( p != a && p != c
508  && p->inTriangle( *a, *b, *c )
509  && area( p->prev, p, p->next ) >= 0 )
510  return false;
511  p = p->prevZ;
512  }
513 
514  return true;
515  }
516 
524  void splitPolygon( Vertex* start )
525  {
526  Vertex* origPoly = start;
527  do
528  {
529  Vertex* marker = origPoly->next->next;
530  while( marker != origPoly->prev )
531  {
532  // Find a diagonal line that is wholly enclosed by the polygon interior
533  if( origPoly->i != marker->i && goodSplit( origPoly, marker ) )
534  {
535  Vertex* newPoly = origPoly->split( marker );
536 
537  origPoly->updateList();
538  newPoly->updateList();
539 
540  earcutList( origPoly );
541  earcutList( newPoly );
542  return;
543  }
544  marker = marker->next;
545  }
546  origPoly = origPoly->next;
547  } while( origPoly != start );
548  }
549 
558  bool goodSplit( const Vertex* a, const Vertex* b ) const
559  {
560  return a->next->i != b->i &&
561  a->prev->i != b->i &&
562  !intersectsPolygon( a, b ) &&
563  locallyInside( a, b );
564  }
565 
571  double area( const Vertex* p, const Vertex* q, const Vertex* r ) const
572  {
573  return ( q->y - p->y ) * ( r->x - q->x ) - ( q->x - p->x ) * ( r->y - q->y );
574  }
575 
581  bool intersects( const Vertex* p1, const Vertex* q1, const Vertex* p2, const Vertex* q2 ) const
582  {
583  if( ( *p1 == *q1 && *p2 == *q2 ) || ( *p1 == *q2 && *p2 == *q1 ) )
584  return true;
585 
586  return ( area( p1, q1, p2 ) > 0 ) != ( area( p1, q1, q2 ) > 0 )
587  && ( area( p2, q2, p1 ) > 0 ) != ( area( p2, q2, q1 ) > 0 );
588  }
589 
596  bool intersectsPolygon( const Vertex* a, const Vertex* b ) const
597  {
598  const Vertex* p = a->next;
599  do
600  {
601  if( p->i != a->i &&
602  p->next->i != a->i &&
603  p->i != b->i &&
604  p->next->i != b->i && intersects( p, p->next, a, b ) )
605  return true;
606 
607  p = p->next;
608  } while( p != a );
609 
610  return false;
611  }
612 
621  bool locallyInside( const Vertex* a, const Vertex* b ) const
622  {
623  if( area( a->prev, a, a->next ) < 0 )
624  return area( a, b, a->next ) >= 0 && area( a, a->prev, b ) >= 0;
625  else
626  return area( a, b, a->prev ) < 0 || area( a, a->next, b ) < 0;
627  }
628 
635  Vertex* insertVertex( const VECTOR2I& pt, Vertex* last )
636  {
637  m_result.AddVertex( pt );
638  m_vertices.emplace_back( m_result.GetVertexCount() - 1, pt.x, pt.y, this );
639 
640  Vertex* p = &m_vertices.back();
641  if( !last )
642  {
643  p->prev = p;
644  p->next = p;
645  }
646  else
647  {
648  p->next = last->next;
649  p->prev = last;
650  last->next->prev = p;
651  last->next = p;
652  }
653  return p;
654  }
655 
656 
657 public:
658 
659  bool TesselatePolygon( const SHAPE_LINE_CHAIN& aPoly )
660  {
661  m_bbox = aPoly.BBox();
662  m_result.Clear();
663 
664  if( !m_bbox.GetWidth() || !m_bbox.GetHeight() )
665  return false;
666 
670  Vertex* firstVertex = createList( aPoly );
671  if( !firstVertex || firstVertex->prev == firstVertex->next )
672  return false;
673 
674  firstVertex->updateList();
675 
676  auto retval = earcutList( firstVertex );
677  m_vertices.clear();
678  return retval;
679  }
680 };
681 
682 #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:189
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:594
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 Point()
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()
#define NULL
void zSort()
Sort all vertices in this vertex's list by their Morton code.
coord_type GetWidth() const
Definition: box2.h:196
SHAPE_POLY_SET::TRIANGULATED_POLYGON & m_result
coord_type GetY() const
Definition: box2.h:190
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)
SHAPE_LINE_CHAIN.
coord_type GetHeight() const
Definition: box2.h:197
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...
void splitPolygon(Vertex *start)
Function splitPolygon If we cannot find an ear to slice in the current polygon list,...