polygon_triangulation.h

Go to the documentation of this file.

/*
 * This program source code file is part of KiCad, a free EDA CAD application.
 *
 * Modifications Copyright (C) 2018-2021 KiCad Developers
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, you may find one here:
 * http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
 * or you may search the http://www.gnu.org website for the version 2 license,
 * or you may write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
 *
 * Based on Uniform Plane Subdivision algorithm from Lamot, Marko, and Borut Žalik.
 * "A fast polygon triangulation algorithm based on uniform plane subdivision."
 * Computers & graphics 27, no. 2 (2003): 239-253.
 *
 * Code derived from:
 * K-3D which is Copyright (c) 2005-2006, Romain Behar, GPL-2, license above
 * earcut which is Copyright (c) 2016, Mapbox, ISC
 *
 * ISC License:
 * Permission to use, copy, modify, and/or distribute this software for any purpose
 * with or without fee is hereby granted, provided that the above copyright notice
 * and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
 * REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
 * FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
 * INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
 * OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
 * TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
 * THIS SOFTWARE.
 *
 */
 
#ifndef __POLYGON_TRIANGULATION_H
#define __POLYGON_TRIANGULATION_H
 
#include <algorithm>
#include <deque>
#include <cmath>
 
#include <clipper.hpp>
#include <geometry/shape_line_chain.h>
#include <geometry/shape_poly_set.h>
#include <math/box2.h>
#include <math/vector2d.h>
 
class PolygonTriangulation
{
public:
 PolygonTriangulation( SHAPE_POLY_SET::TRIANGULATED_POLYGON& aResult ) :
 m_result( aResult )
 {};
 
 bool TesselatePolygon( const SHAPE_LINE_CHAIN& aPoly )
 {
 m_bbox = aPoly.BBox();
 m_result.Clear();
 
 if( !m_bbox.GetWidth() || !m_bbox.GetHeight() )
 return false;
 
 Vertex* firstVertex = createList( aPoly );
 if( !firstVertex || firstVertex->prev == firstVertex->next )
 return false;
 
 firstVertex->updateList();
 
 auto retval = earcutList( firstVertex );
 m_vertices.clear();
 return retval;
 }
 
private:
 struct Vertex
 {
 Vertex( size_t aIndex, double aX, double aY, PolygonTriangulation* aParent ) :
 i( aIndex ),
 x( aX ),
 y( aY ),
 parent( aParent )
 {
 }
 
 Vertex& operator=( const Vertex& ) = delete;
 Vertex& operator=( Vertex&& ) = delete;
 
 bool operator==( const Vertex& rhs ) const
 {
 return this->x == rhs.x && this->y == rhs.y;
 }
 bool operator!=( const Vertex& rhs ) const { return !( *this == rhs ); }
 
 
 Vertex* split( Vertex* b )
 {
 parent->m_vertices.emplace_back( i, x, y, parent );
 Vertex* a2 = &parent->m_vertices.back();
 parent->m_vertices.emplace_back( b->i, b->x, b->y, parent );
 Vertex* b2 = &parent->m_vertices.back();
 Vertex* an = next;
 Vertex* bp = b->prev;
 
 next = b;
 b->prev = this;
 
 a2->next = an;
 an->prev = a2;
 
 b2->next = a2;
 a2->prev = b2;
 
 bp->next = b2;
 b2->prev = bp;
 
 return b2;
 }
 
 void remove()
 {
 next->prev = prev;
 prev->next = next;
 
 if( prevZ )
 prevZ->nextZ = nextZ;
 
 if( nextZ )
 nextZ->prevZ = prevZ;
 
 next = nullptr;
 prev = nullptr;
 nextZ = nullptr;
 prevZ = nullptr;
 }
 
 void updateOrder()
 {
 if( !z )
 z = parent->zOrder( x, y );
 }
 
 void updateList()
 {
 Vertex* p = next;
 
 while( p != this )
 {
 if( *p == *p->next )
 {
 p = p->prev;
 p->next->remove();
 
 if( p == p->next )
 break;
 }
 
 p->updateOrder();
 p = p->next;
 };
 
 updateOrder();
 zSort();
 }
 
 void zSort()
 {
 std::deque<Vertex*> queue;
 
 queue.push_back( this );
 
 for( auto p = next; p && p != this; p = p->next )
 queue.push_back( p );
 
 std::sort( queue.begin(), queue.end(), []( const Vertex* a, const Vertex* b )
 {
 if( a->z != b->z )
 return a->z < b->z;
 
 if( a->x != b->x )
 return a->x < b->x;
 
 if( a->y != b->y )
 return a->y < b->y;
 
 return a->i < b->i;
 } );
 
 Vertex* prev_elem = nullptr;
 
 for( auto elem : queue )
 {
 if( prev_elem )
 prev_elem->nextZ = elem;
 
 elem->prevZ = prev_elem;
 prev_elem = elem;
 }
 
 prev_elem->nextZ = nullptr;
 }
 
 
 bool inTriangle( const Vertex& a, const Vertex& b, const Vertex& c )
 {
 return ( c.x - x ) * ( a.y - y ) - ( a.x - x ) * ( c.y - y ) >= 0
 && ( a.x - x ) * ( b.y - y ) - ( b.x - x ) * ( a.y - y ) >= 0
 && ( b.x - x ) * ( c.y - y ) - ( c.x - x ) * ( b.y - y ) >= 0;
 }
 
 const size_t i;
 const double x;
 const double y;
 PolygonTriangulation* parent;
 
 // previous and next vertices nodes in a polygon ring
 Vertex* prev = nullptr;
 Vertex* next = nullptr;
 
 // z-order curve value
 int32_t z = 0;
 
 // previous and next nodes in z-order
 Vertex* prevZ = nullptr;
 Vertex* nextZ = nullptr;
 };
 
 int32_t zOrder( const double aX, const double aY ) const
 {
 int32_t x = static_cast<int32_t>( 32767.0 * ( aX - m_bbox.GetX() ) / m_bbox.GetWidth() );
 int32_t y = static_cast<int32_t>( 32767.0 * ( aY - m_bbox.GetY() ) / m_bbox.GetHeight() );
 
 x = ( x | ( x << 8 ) ) & 0x00FF00FF;
 x = ( x | ( x << 4 ) ) & 0x0F0F0F0F;
 x = ( x | ( x << 2 ) ) & 0x33333333;
 x = ( x | ( x << 1 ) ) & 0x55555555;
 
 y = ( y | ( y << 8 ) ) & 0x00FF00FF;
 y = ( y | ( y << 4 ) ) & 0x0F0F0F0F;
 y = ( y | ( y << 2 ) ) & 0x33333333;
 y = ( y | ( y << 1 ) ) & 0x55555555;
 
 return x | ( y << 1 );
 }
 
 Vertex* removeNullTriangles( Vertex* aStart )
 {
 Vertex* retval = nullptr;
 Vertex* p = aStart->next;
 
 while( p != aStart )
 {
 if( area( p->prev, p, p->next ) == 0.0 )
 {
 p = p->prev;
 p->next->remove();
 retval = aStart;
 
 if( p == p->next )
 break;
 }
 p = p->next;
 };
 
 // We needed an end point above that wouldn't be removed, so
 // here we do the final check for this as a Steiner point
 if( area( aStart->prev, aStart, aStart->next ) == 0.0 )
 {
 retval = p->next;
 p->remove();
 }
 
 return retval;
 }
 
 Vertex* createList( const ClipperLib::Path& aPath )
 {
 Vertex* tail = nullptr;
 double sum = 0.0;
 auto len = aPath.size();
 
 // Check for winding order
 for( size_t i = 0; i < len; i++ )
 {
 auto p1 = aPath.at( i );
 auto p2 = aPath.at( ( i + 1 ) < len ? i + 1 : 0 );
 
 sum += ( ( p2.X - p1.X ) * ( p2.Y + p1.Y ) );
 }
 
 if( sum <= 0.0 )
 {
 for( auto point : aPath )
 tail = insertVertex( VECTOR2I( point.X, point.Y ), tail );
 }
 else
 {
 for( size_t i = 0; i < len; i++ )
 {
 auto p = aPath.at( len - i - 1 );
 tail = insertVertex( VECTOR2I( p.X, p.Y ), tail );
 }
 }
 
 if( tail && ( *tail == *tail->next ) )
 {
 tail->next->remove();
 }
 
 return tail;
 
 }
 
 Vertex* createList( const SHAPE_LINE_CHAIN& points )
 {
 Vertex* tail = nullptr;
 double sum = 0.0;
 
 // Check for winding order
 for( int i = 0; i < points.PointCount(); i++ )
 {
 VECTOR2D p1 = points.CPoint( i );
 VECTOR2D p2 = points.CPoint( i + 1 );
 
 sum += ( ( p2.x - p1.x ) * ( p2.y + p1.y ) );
 }
 
 if( sum > 0.0 )
 for( int i = points.PointCount() - 1; i >= 0; i--)
 tail = insertVertex( points.CPoint( i ), tail );
 else
 for( int i = 0; i < points.PointCount(); i++ )
 tail = insertVertex( points.CPoint( i ), tail );
 
 if( tail && ( *tail == *tail->next ) )
 {
 tail->next->remove();
 }
 
 return tail;
 }
 
 bool earcutList( Vertex* aPoint, int pass = 0 )
 {
 if( !aPoint )
 return true;
 
 Vertex* stop = aPoint;
 Vertex* prev;
 Vertex* next;
 
 while( aPoint->prev != aPoint->next )
 {
 prev = aPoint->prev;
 next = aPoint->next;
 
 if( isEar( aPoint ) )
  {
 m_result.AddTriangle( prev->i, aPoint->i, next->i );
 aPoint->remove();
 
 // Skip one vertex as the triangle will account for the prev node
 aPoint = next->next;
 stop = next->next;
 
 continue;
 }
 
 Vertex* nextNext = next->next;
 
 if( *prev != *nextNext && intersects( prev, aPoint, next, nextNext ) &&
 locallyInside( prev, nextNext ) &&
 locallyInside( nextNext, prev ) )
 {
 m_result.AddTriangle( prev->i, aPoint->i, nextNext->i );
 
 // remove two nodes involved
 next->remove();
 aPoint->remove();
 
 aPoint = nextNext;
 stop = nextNext;
 
 continue;
 }
 
 aPoint = next;
 
 /*
 * We've searched the entire polygon for available ears and there are still
 * un-sliced nodes remaining.
 */
 if( aPoint == stop )
 {
 // First, try to remove the remaining steiner points
 // If aPoint is a steiner, we need to re-assign both the start and stop points
 if( auto newPoint = removeNullTriangles( aPoint ) )
 {
 aPoint = newPoint;
 stop = newPoint;
 continue;
 }
 
 // If we don't have any NULL triangles left, cut the polygon in two and try again
 splitPolygon( aPoint );
 break;
 }
 }
 
 // Check to see if we are left with only three points in the polygon
 if( aPoint->next && aPoint->prev == aPoint->next->next )
 {
 // Three concave points will never be able to be triangulated because they were
 // created by an intersecting polygon, so just drop them.
 if( area( aPoint->prev, aPoint, aPoint->next ) >= 0 )
 return true;
 }
 
 /*
 * At this point, our polygon should be fully tessellated.
 */
 return( aPoint->prev == aPoint->next );
 }
 
 bool isEar( Vertex* aEar ) const
 {
 const Vertex* a = aEar->prev;
 const Vertex* b = aEar;
 const Vertex* c = aEar->next;
 
 // If the area >=0, then the three points for a concave sequence
 // with b as the reflex point
 if( area( a, b, c ) >= 0 )
 return false;
 
 // triangle bbox
 const double minTX = std::min( a->x, std::min( b->x, c->x ) );
 const double minTY = std::min( a->y, std::min( b->y, c->y ) );
 const double maxTX = std::max( a->x, std::max( b->x, c->x ) );
 const double maxTY = std::max( a->y, std::max( b->y, c->y ) );
 
 // z-order range for the current triangle bounding box
 const int32_t minZ = zOrder( minTX, minTY );
 const int32_t maxZ = zOrder( maxTX, maxTY );
 
 // first look for points inside the triangle in increasing z-order
 Vertex* p = aEar->nextZ;
 
 while( p && p->z <= maxZ )
 {
 if( p != a && p != c
 && p->inTriangle( *a, *b, *c )
 && area( p->prev, p, p->next ) >= 0 )
 return false;
 
 p = p->nextZ;
 }
 
 // then look for points in decreasing z-order
 p = aEar->prevZ;
 
 while( p && p->z >= minZ )
 {
 if( p != a && p != c
 && p->inTriangle( *a, *b, *c )
 && area( p->prev, p, p->next ) >= 0 )
 return false;
 
 p = p->prevZ;
 }
 
 return true;
 }
 
 void splitPolygon( Vertex* start )
 {
 Vertex* origPoly = start;
 
 do
 {
 Vertex* marker = origPoly->next->next;
 
 while( marker != origPoly->prev )
 {
 // Find a diagonal line that is wholly enclosed by the polygon interior
 if( origPoly->i != marker->i && goodSplit( origPoly, marker ) )
 {
 Vertex* newPoly = origPoly->split( marker );
 
 origPoly->updateList();
 newPoly->updateList();
 
 earcutList( origPoly );
 earcutList( newPoly );
 return;
 }
 
 marker = marker->next;
 }
 
 origPoly = origPoly->next;
 } while( origPoly != start );
 }
 
 bool goodSplit( const Vertex* a, const Vertex* b ) const
 {
 return a->next->i != b->i &&
 a->prev->i != b->i &&
 !intersectsPolygon( a, b ) &&
 locallyInside( a, b );
 }
 
 double area( const Vertex* p, const Vertex* q, const Vertex* r ) const
 {
 return ( q->y - p->y ) * ( r->x - q->x ) - ( q->x - p->x ) * ( r->y - q->y );
 }
 
 bool intersects( const Vertex* p1, const Vertex* q1, const Vertex* p2, const Vertex* q2 ) const
 {
 if( ( *p1 == *q1 && *p2 == *q2 ) || ( *p1 == *q2 && *p2 == *q1 ) )
 return true;
 
 return ( area( p1, q1, p2 ) > 0 ) != ( area( p1, q1, q2 ) > 0 )
 && ( area( p2, q2, p1 ) > 0 ) != ( area( p2, q2, q1 ) > 0 );
 }
 
 bool intersectsPolygon( const Vertex* a, const Vertex* b ) const
 {
 const Vertex* p = a->next;
 
 do
 {
 if( p->i != a->i &&
 p->next->i != a->i &&
 p->i != b->i &&
 p->next->i != b->i && intersects( p, p->next, a, b ) )
 return true;
 
 p = p->next;
 } while( p != a );
 
 return false;
 }
 
 bool locallyInside( const Vertex* a, const Vertex* b ) const
 {
 if( area( a->prev, a, a->next ) < 0 )
 return area( a, b, a->next ) >= 0 && area( a, a->prev, b ) >= 0;
 else
 return area( a, b, a->prev ) < 0 || area( a, a->next, b ) < 0;
 }
 
 Vertex* insertVertex( const VECTOR2I& pt, Vertex* last )
 {
 m_result.AddVertex( pt );
 m_vertices.emplace_back( m_result.GetVertexCount() - 1, pt.x, pt.y, this );
 
 Vertex* p = &m_vertices.back();
 
 if( !last )
 {
 p->prev = p;
 p->next = p;
 }
 else
 {
 p->next = last->next;
 p->prev = last;
 last->next->prev = p;
 last->next = p;
 }
 return p;
 }
 
private:
 BOX2I m_bbox;
 std::deque<Vertex> m_vertices;
 SHAPE_POLY_SET::TRIANGULATED_POLYGON& m_result;
};
 
#endif //__POLYGON_TRIANGULATION_H