klayout/src/db/db/dbPolygonTools.cc

2537 lines
70 KiB
C++

/*
KLayout Layout Viewer
Copyright (C) 2006-2018 Matthias Koefferlein
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, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
#include "dbCommon.h"
#include "dbPolygonTools.h"
#include "dbPolygonGenerators.h"
#include "tlLog.h"
#include "tlInt128Support.h"
#include <algorithm>
#include <cmath>
namespace db
{
// -------------------------------------------------------------------------
// Implementation of inside_poly_test
// internal: compare edges by higher y coordinate
template <class E>
struct inside_poly_test_edge_max_compare_f
{
bool operator() (const E &e1, const E &e2)
{
return std::max (e1.p1 ().y (), e1.p2 ().y ()) < std::max (e2.p1 ().y (), e2.p2 ().y ());
}
};
template<class P>
inside_poly_test<P>::inside_poly_test (const P &polygon)
{
m_edges.reserve (polygon.vertices ());
for (typename P::polygon_edge_iterator e = polygon.begin_edge (); ! e.at_end (); ++e) {
m_edges.push_back (*e);
}
std::sort (m_edges.begin (), m_edges.end (), inside_poly_test_edge_max_compare_f<edge_type> ());
}
template<class P>
int inside_poly_test<P>::operator() (const typename inside_poly_test<P>::point_type &pt) const
{
int wrapcount_left = 0;
typename std::vector<edge_type>::const_iterator e = std::lower_bound (m_edges.begin (), m_edges.end (), edge_type (pt, pt), inside_poly_test_edge_max_compare_f<edge_type> ());
while (e != m_edges.end () && pt.y () <= std::max (e->p1 ().y (), e->p2 ().y ())) {
if ((*e).p1 ().y () <= pt.y () && (*e).p2 ().y () > pt.y ()) {
int side = (*e).side_of (pt);
if (side < 0) {
++wrapcount_left;
} else if (side == 0) {
// "on" the line is excluded in the predicate
return 0;
}
} else if ((*e).p2 ().y () <= pt.y () && (*e).p1 ().y () > pt.y ()) {
int side = (*e).side_of (pt);
if (side > 0) {
--wrapcount_left;
} else if (side == 0) {
// "on" the line is excluded in the predicate
return 0;
}
} else if ((*e).p1 ().y () == pt.y () && (*e).p2 ().y () == pt.y () &&
(((*e).p1 ().x () <= pt.x () && (*e).p2 ().x () >= pt.x ()) ||
((*e).p2 ().x () <= pt.x () && (*e).p1 ().x () >= pt.x ()))) {
// "on" the horizontal line is excluded in the predicate
return 0;
}
++e;
}
return (wrapcount_left != 0) ? 1 : -1;
}
template class DB_PUBLIC inside_poly_test<db::SimplePolygon>;
template class DB_PUBLIC inside_poly_test<db::Polygon>;
template class DB_PUBLIC inside_poly_test<db::DSimplePolygon>;
template class DB_PUBLIC inside_poly_test<db::DPolygon>;
// -------------------------------------------------------------------------
// Implementation of cut_polygon
/**
* @brief A helper structure describing an edge cutting the cut line.
*/
template <class PointType>
struct cut_polygon_edge
{
typedef typename PointType::coord_type coord_type;
typedef typename db::edge<coord_type> edge_type;
typedef typename db::coord_traits<coord_type>::area_type projection_type;
cut_polygon_edge ()
: contour (-1), index (0), projected (0), point (), last_point ()
{ }
cut_polygon_edge (int c, unsigned int n, projection_type p, const PointType &pt, const PointType &lpt)
: contour (c), index (n), projected (p), point (pt), last_point (lpt)
{ }
int contour;
unsigned int index;
projection_type projected;
PointType point;
PointType last_point;
};
template <class CuttingEdgeType>
struct cut_polygon_segment
{
typedef typename CuttingEdgeType::edge_type edge_type;
typedef typename CuttingEdgeType::projection_type projection_type;
CuttingEdgeType leave;
CuttingEdgeType enter;
int segment;
bool hole;
cut_polygon_segment ()
: leave (), enter (), segment (-1), hole (false)
{
// .. nothing yet ..
}
bool operator< (const cut_polygon_segment &b) const
{
projection_type amin = std::min (leave.projected, enter.projected);
projection_type bmin = std::min (b.leave.projected, b.enter.projected);
projection_type amax = std::max (leave.projected, enter.projected);
projection_type bmax = std::max (b.leave.projected, b.enter.projected);
if (amin != bmin) {
return amin < bmin;
} else {
int vs = db::vprod_sign (min_edge (), b.min_edge ());
if (vs != 0) {
return vs > 0;
} else if (amax == amin) {
return bmax != bmin;
} else if (bmax == bmin) {
return false;
} else {
return amax > bmax;
}
}
}
bool after (const cut_polygon_segment &b) const
{
projection_type pb_max = std::max (b.leave.projected, b.enter.projected);
projection_type pa_min = std::min (leave.projected, enter.projected);
if (pb_max < pa_min) {
return true;
} else if (pb_max > pa_min) {
return false;
} else {
return db::vprod_sign (min_edge (), b.max_edge ()) <= 0;
}
}
edge_type max_edge () const
{
if (leave.projected > enter.projected) {
return edge_type (leave.last_point, leave.point);
} else if (leave.projected < enter.projected) {
return edge_type (enter.last_point, enter.point);
} else {
int vs = db::vprod_sign (leave.point - leave.last_point, enter.point - enter.last_point);
if (vs < 0) {
return edge_type (leave.last_point, leave.point);
} else {
return edge_type (enter.last_point, enter.point);
}
}
}
edge_type min_edge () const
{
if (leave.projected < enter.projected) {
return edge_type (leave.last_point, leave.point);
} else if (leave.projected > enter.projected) {
return edge_type (enter.last_point, enter.point);
} else {
int vs = db::vprod_sign (leave.point - leave.last_point, enter.point - enter.last_point);
if (vs > 0) {
return edge_type (leave.last_point, leave.point);
} else {
return edge_type (enter.last_point, enter.point);
}
}
}
};
template <class CuttingEdge, class CuttingEdgeSegment>
static void cut_polygon_produce_cut_points (typename std::vector<CuttingEdgeSegment>::const_iterator from, typename std::vector<CuttingEdgeSegment>::const_iterator to, std::vector<CuttingEdge> &output, bool holes)
{
typename std::vector<CuttingEdgeSegment>::const_iterator e = from;
while (e != to) {
// collect "illegal" segments (i.e. holes not being bracketed by a contour
typename std::vector<CuttingEdgeSegment>::const_iterator e1 = e;
for ( ; e != to && e->hole != holes; ++e) {
;
}
if (e == to) {
// In this case we have only wrong typed segments within the wrong bracketing segment (i.e. holes inside a hole bracket) or on the outside.
// Sort the leaving and entering edges correctly.
for (typename std::vector<CuttingEdgeSegment>::const_iterator o = e1; o != e; ++o) {
if (! holes) {
output.push_back (o->leave);
output.push_back (o->enter);
} else {
output.push_back (o->enter);
output.push_back (o->leave);
}
}
return;
}
typename std::vector<CuttingEdgeSegment>::const_iterator ee = e + 1;
for ( ; ee != to; ++ee) {
if (e->hole == holes && ee->after (*e)) {
break;
}
}
if (e->hole) {
output.push_back (e->enter);
} else {
output.push_back (e->leave);
}
// include all illegal segments into the first "good" segment, so we have a correct hierarchy at least.
if (e1 != e) {
cut_polygon_produce_cut_points<CuttingEdge, CuttingEdgeSegment> (e1, e, output, !holes);
}
if (e + 1 != ee) {
cut_polygon_produce_cut_points<CuttingEdge, CuttingEdgeSegment> (e + 1, ee, output, !holes);
}
if (e->hole) {
output.push_back (e->leave);
} else {
output.push_back (e->enter);
}
e = ee;
}
}
template <class PolygonType, class Edge>
void cut_polygon_internal (const PolygonType &input, const Edge &line, CutPolygonReceiverBase *right_of_line)
{
typedef typename PolygonType::point_type point_type;
typedef typename PolygonType::coord_type coord_type;
typedef typename PolygonType::area_type area_type;
typedef typename PolygonType::contour_type contour_type;
typedef db::edge<coord_type> edge_type;
typedef cut_polygon_edge<point_type> cut_polygon_edge_type;
typedef cut_polygon_segment<cut_polygon_edge_type> cutting_segment_type;
typedef typename db::coord_traits<coord_type>::area_type projection_type;
bool do_hole_assignment = (input.holes () > 0);
std::vector <PolygonType> hull_polygons;
std::vector <PolygonType> hole_polygons;
std::vector<cutting_segment_type> cutting_segments;
for (unsigned int nc = 0; nc < input.holes () + 1; ++nc) {
const contour_type &contour = input.contour (nc);
if (contour.size () == 0) {
continue;
}
bool any = false;
unsigned int nn = (unsigned int) contour.size () - 1;
int sc = -1;
size_t nfirst = cutting_segments.size ();
point_type last_pt = contour[nn];
for (unsigned int n = 0; n < contour.size (); ++n) {
edge_type e (last_pt, contour[n]);
last_pt = e.p2 ();
std::pair <bool, point_type> ip = line.crossed_by_point (e);
if (ip.first) {
int s1 = line.side_of (e.p1 ());
int s2 = line.side_of (e.p2 ());
projection_type p = db::sprod (ip.second - line.p1 (), line.p2 () - line.p1 ());
if (s1 < 0 && s2 >= 0) {
// right -> left or on edge
if (cutting_segments.size () == nfirst) {
cutting_segments.push_back (cutting_segment_type ());
}
cutting_segments.back ().leave = cut_polygon_edge_type (int (nc), nn, p, ip.second, e.p1 ());
any = true;
}
if (s1 >= 0 && s2 < 0) {
// left or on edge -> right
++sc;
cutting_segments.push_back (cutting_segment_type ());
cutting_segments.back ().enter = cut_polygon_edge_type (int (nc), nn, p, ip.second, e.p2 ());
cutting_segments.back ().segment = sc;
}
}
nn = n;
}
if (any) {
// tie together last and first partial segments.
if (cutting_segments[nfirst].segment < 0) {
cutting_segments[nfirst].enter = cutting_segments.back ().enter;
cutting_segments.pop_back ();
}
// assign hull/hole identities ..
for (typename std::vector<cutting_segment_type>::iterator c = cutting_segments.begin () + nfirst; c != cutting_segments.end (); ++c) {
unsigned int n = (unsigned int) contour.size ();
unsigned int n1 = c->enter.index;
unsigned int n2 = c->leave.index;
// compute area of the segment and derive the hull part/hole part nature of the segment from the sign
area_type a = 0;
point_type pl = c->enter.point;
for (unsigned int i = n1; ; ) {
a += db::vprod (pl - point_type (), contour[i] - point_type ());
pl = contour[i];
if (i == n2) {
break;
}
++i;
if (i == n) {
i = 0;
}
}
a += db::vprod (pl - point_type (), c->leave.point - point_type ());
a += db::vprod (c->leave.point - point_type (), c->enter.point - point_type ());
c->hole = (a > 0);
}
} else if (line.side_of (contour[0]) < 0) {
if (do_hole_assignment) {
if (nc == 0) {
// the hull is fully on the right side -> just output the input polygon and that's it.
right_of_line->put (&input);
return;
} else {
// remember hole contours for later assignment
hole_polygons.push_back (PolygonType ());
hole_polygons.back ().assign_hull (contour.begin (), contour.end ());
}
} else {
PolygonType poly;
poly.assign_hull (contour.begin (), contour.end ());
right_of_line->put (&poly);
}
}
}
std::sort (cutting_segments.begin (), cutting_segments.end ());
std::vector<cut_polygon_edge_type> cutting_edges;
cutting_edges.reserve (cutting_segments.size () * 2);
cut_polygon_produce_cut_points<cut_polygon_edge_type, cutting_segment_type> (typename std::vector<cutting_segment_type>::const_iterator (cutting_segments.begin ()), typename std::vector<cutting_segment_type>::const_iterator (cutting_segments.end ()), cutting_edges, false);
std::vector<point_type> contour_points;
typedef std::map <std::pair<int, int>, std::pair<typename std::vector<cut_polygon_edge_type>::iterator, typename std::vector<cut_polygon_edge_type>::iterator> > cut_point_map_type;
cut_point_map_type cut_points;
for (typename std::vector<cut_polygon_edge_type>::iterator c = cutting_edges.begin (); c != cutting_edges.end (); c += 2) {
cut_points.insert (std::make_pair (std::make_pair (c[0].contour, c[0].index), std::make_pair (c, c + 1)));
cut_points.insert (std::make_pair (std::make_pair (c[1].contour, c[1].index), std::make_pair (c + 1, c)));
}
for (typename std::vector<cut_polygon_edge_type>::iterator c = cutting_edges.begin (); c != cutting_edges.end (); c += 2) {
if (c->contour >= 0) {
typename std::vector<cut_polygon_edge_type>::iterator c1 = c;
typename std::vector<cut_polygon_edge_type>::iterator c2 = c + 1;
contour_points.clear ();
bool is_hull = false;
do {
tl_assert (c1->contour >= 0 && c2->contour >= 0);
contour_points.push_back (c1->point);
contour_points.push_back (c2->point);
int n = c2->index;
int n0 = n;
int nc = c2->contour;
const contour_type &contour = input.contour (nc);
if (nc == 0) {
is_hull = true;
}
c1->contour = -1;
c2->contour = -1;
++n;
if (n == int (contour.size ())) {
n = 0;
}
while (n != n0) {
contour_points.push_back (contour[n]);
typename cut_point_map_type::iterator cp = cut_points.find (std::make_pair (nc, n));
if (cp != cut_points.end ()) {
c1 = cp->second.first;
c2 = cp->second.second;
break;
}
++n;
if (n == int (contour.size ())) {
n = 0;
}
}
tl_assert (n != n0);
} while (c1 != c && c2 != c);
// Hint: the algorithm tends to create spikes for hole insertion edges.
// Therefore we used "remove reflected" on the assignment.
if (do_hole_assignment) {
if (is_hull) {
hull_polygons.push_back (PolygonType ());
hull_polygons.back ().assign_hull (contour_points.begin (), contour_points.end (), true, true /*remove reflected*/);
} else {
hole_polygons.push_back (PolygonType ());
hole_polygons.back ().assign_hull (contour_points.begin (), contour_points.end (), true, true /*remove reflected*/);
}
} else {
PolygonType poly;
poly.assign_hull (contour_points.begin (), contour_points.end (), true, true /*remove reflected*/);
// it might happen in some cases, that cut pieces may vanish (i.e. all points on a line). Thus we check, if that
// is the case and do not produce a polygon then.
if (poly.vertices () > 0) {
right_of_line->put (&poly);
}
}
}
}
// do hole assignment
for (typename std::vector<PolygonType>::iterator hull = hull_polygons.begin (); hull != hull_polygons.end (); ++hull) {
// it might happen in some cases, that cut pieces may vanish (i.e. all points on a line). Thus we check, if that
// is the case and do not produce a polygon then.
if (hull->vertices () > 0) {
db::inside_poly_test<PolygonType> inside_hull (*hull);
for (typename std::vector<PolygonType>::iterator hole = hole_polygons.begin (); hole != hole_polygons.end (); ++hole) {
size_t n = hole->hull ().size ();
if (n > 0) {
// look for one point "really" inside ...
int inside = 0;
for (size_t i = 0; i < n && inside == 0; ++i) {
inside = inside_hull (hole->hull ()[i]);
}
if (inside >= 0) {
hull->insert_hole (hole->hull ().begin (), hole->hull ().end ());
*hole = PolygonType ();
}
}
}
right_of_line->put (&*hull);
}
}
// use non-assigned hole (parts) as hulls
// TODO: precisely, this is possible only if the orientation is clockwise. Since we loose the orientation because
// we assign to a PolygonType, this check is not possible.
for (typename std::vector<PolygonType>::iterator hole = hole_polygons.begin (); hole != hole_polygons.end (); ++hole) {
if (hole->vertices () > 0) {
right_of_line->put (&*hole);
}
}
}
template DB_PUBLIC void cut_polygon_internal<> (const db::Polygon &polygon, const db::Polygon::edge_type &line, CutPolygonReceiverBase *right_of_line);
template DB_PUBLIC void cut_polygon_internal<> (const db::SimplePolygon &polygon, const db::SimplePolygon::edge_type &line, CutPolygonReceiverBase *right_of_line);
template DB_PUBLIC void cut_polygon_internal<> (const db::DPolygon &polygon, const db::DPolygon::edge_type &line, CutPolygonReceiverBase *right_of_line);
template DB_PUBLIC void cut_polygon_internal<> (const db::DSimplePolygon &polygon, const db::DSimplePolygon::edge_type &line, CutPolygonReceiverBase *right_of_line);
// -------------------------------------------------------------------------
// Implementation of split_polygon
template <class PolygonType>
void
split_polygon (const PolygonType &polygon, std::vector<PolygonType> &output)
{
typedef typename PolygonType::coord_type coord_type;
typedef typename PolygonType::point_type point_type;
typedef typename PolygonType::box_type box_type;
typedef db::edge<coord_type> edge_type;
box_type bbox = polygon.box ();
box_type b1, b2;
coord_type x = bbox.center ().x ();
coord_type xx = x;
bool xx_set = false;
coord_type y = bbox.center ().y ();
coord_type yy = y;
bool yy_set = false;
for (typename PolygonType::polygon_contour_iterator e = polygon.begin_hull (); e != polygon.end_hull (); ++e) {
if ((*e).x () != bbox.left () && (*e).x () != bbox.right () && (std::abs ((*e).x () - x) < std::abs (xx - x) || ! xx_set)) {
xx = (*e).x ();
xx_set = true;
}
if ((*e).y () != bbox.top () && (*e).y () != bbox.bottom () && (std::abs ((*e).y () - y) < std::abs (yy - y) || ! yy_set)) {
yy = (*e).y ();
yy_set = true;
}
}
if (! xx_set && ! yy_set) {
if (bbox.width () > bbox.height ()) {
xx_set = true;
} else {
yy_set = true;
}
} else if (xx_set && yy_set) {
// an empiric threshold for splitting polygons in one direction: don't split along the long
// axis for polygons with a aspect ratio (of the bounding box) of larger than 3
if (bbox.width () > 3 * bbox.height ()) {
yy_set = false;
} else if (bbox.height () > 3 * bbox.width ()) {
xx_set = false;
}
}
std::vector <PolygonType> xx_polygons;
size_t xx_n = std::numeric_limits<size_t>::max ();
if (xx_set) {
db::cut_polygon (polygon, edge_type (point_type (xx, 0), point_type (xx, 1)), std::back_inserter (xx_polygons));
db::cut_polygon (polygon, edge_type (point_type (xx, 1), point_type (xx, 0)), std::back_inserter (xx_polygons));
xx_n = 0;
for (typename std::vector <PolygonType>::const_iterator p = xx_polygons.begin (); p != xx_polygons.end (); ++p) {
xx_n += p->vertices ();
}
}
std::vector <PolygonType> yy_polygons;
size_t yy_n = std::numeric_limits<size_t>::max ();
if (yy_set) {
db::cut_polygon (polygon, edge_type (point_type (0, yy), point_type (1, yy)), std::back_inserter (yy_polygons));
db::cut_polygon (polygon, edge_type (point_type (1, yy), point_type (0, yy)), std::back_inserter (yy_polygons));
yy_n = 0;
for (typename std::vector <PolygonType>::const_iterator p = yy_polygons.begin (); p != yy_polygons.end (); ++p) {
yy_n += p->vertices ();
}
}
if (xx_n < yy_n) {
output.swap (xx_polygons);
} else {
output.swap (yy_polygons);
}
}
template DB_PUBLIC void split_polygon<> (const db::Polygon &polygon, std::vector<db::Polygon> &output);
template DB_PUBLIC void split_polygon<> (const db::SimplePolygon &polygon, std::vector<db::SimplePolygon> &output);
template DB_PUBLIC void split_polygon<> (const db::DPolygon &polygon, std::vector<db::DPolygon> &output);
template DB_PUBLIC void split_polygon<> (const db::DSimplePolygon &polygon, std::vector<db::DSimplePolygon> &output);
// -------------------------------------------------------------------------
// Smoothing tools
void
smooth_contour (db::Polygon::polygon_contour_iterator from, db::Polygon::polygon_contour_iterator to, std::vector <db::Point> &points, db::Coord d)
{
points.clear ();
points.reserve (std::distance (from, to));
std::vector<db::Point> org_points;
std::vector<size_t> point_indexes;
point_indexes.reserve (std::distance (from, to));
// collect the points into a vector
size_t pi = 0;
for (db::Polygon::polygon_contour_iterator p = from; p != to; ++p, ++pi) {
points.push_back (*p);
point_indexes.push_back (pi);
}
org_points = points;
// proceed until there is nothing to do
bool even = false;
int cont = 2;
while (points.size () >= 3 && cont > 0) {
std::vector<db::Point> new_points;
new_points.reserve (points.size ());
std::vector<size_t> new_point_indexes;
new_point_indexes.reserve (points.size ());
bool any = false;
size_t i;
bool first_point_deleted = false;
for (i = (even ? 0 : 1); i < points.size (); i += 2) {
if (i == points.size () - 1 && first_point_deleted) {
break;
}
db::Point pm1 = points [(i + points.size () - 2) % points.size ()];
db::Point p0 = points [(i + points.size () - 1) % points.size ()];
db::Point p1 = points [i];
db::Point p2 = points [(i + 1) % points.size ()];
size_t pi0 = point_indexes [(i + points.size () - 1) % points.size ()];
size_t pi1 = point_indexes [i];
size_t pi2 = point_indexes [(i + 1) % points.size ()];
if (i > 0) {
new_points.push_back (p0);
new_point_indexes.push_back (pi0);
}
bool can_drop = false;
if (db::Coord (p1.distance(p0)) <= d && db::sprod_sign (p2 - p1, p0 - pm1) > 0 && std::abs (db::vprod (p2 - p1, p0 - pm1)) < 0.8 * p2.distance (p1) * p0.distance (pm1)) {
// jog configurations with small edges are candidates
can_drop = true;
} else if (db::vprod_sign (p2 - p1, p1 - p0) < 0) {
// concave corners are always candidates
can_drop = true;
} else {
// convex corners enclosing a little more than 45 degree are candidates too
can_drop = (db::sprod_sign (p2 - p1, p1 - p0) > 0 && std::abs (db::vprod (p2 - p1, p1 - p0)) < 0.8 * p2.distance (p1) * p1.distance (p0));
}
for (size_t j = pi0; can_drop; ) {
if (std::abs (db::Edge (p0, p2).distance (org_points [j])) > d) {
can_drop = false;
}
if (j == pi2) {
break;
}
++j;
if (j == org_points.size ()) {
j = 0;
}
}
if (can_drop) {
// drop this point
any = true;
if (i == 0) {
first_point_deleted = true;
}
} else {
new_points.push_back (p1);
new_point_indexes.push_back (pi1);
}
}
if (any) {
cont = 2;
} else {
cont -= 1;
}
while (i <= points.size ()) {
new_points.push_back (points [i - 1]);
new_point_indexes.push_back (point_indexes [i - 1]);
++i;
}
points.swap (new_points);
point_indexes.swap (new_point_indexes);
even = !even;
}
}
db::Polygon
smooth (const db::Polygon &polygon, db::Coord d)
{
db::Polygon new_poly;
std::vector <db::Point> new_pts;
smooth_contour (polygon.begin_hull (), polygon.end_hull (), new_pts, d);
if (new_pts.size () >= 3) {
new_poly.assign_hull (new_pts.begin (), new_pts.end (), false /*don't compress*/);
for (unsigned int h = 0; h < polygon.holes (); ++h) {
new_pts.clear ();
smooth_contour (polygon.begin_hole (h), polygon.end_hole (h), new_pts, d);
if (new_pts.size () >= 3) {
new_poly.insert_hole (new_pts.begin (), new_pts.end (), false /*don't compress*/);
}
}
}
return new_poly;
}
// -------------------------------------------------------------------------
// Rounding tools
template <class C>
static bool
do_extract_rad_from_contour (typename db::polygon<C>::polygon_contour_iterator from, typename db::polygon<C>::polygon_contour_iterator to, double &rinner, double &router, unsigned int &n, std::vector <db::point<C> > *new_pts, bool fallback)
{
if (from == to) {
return false;
}
typename db::polygon<C>::polygon_contour_iterator p0 = from;
typename db::polygon<C>::polygon_contour_iterator p1 = p0;
++p1;
if (p1 == to) {
p1 = from;
}
typename db::polygon<C>::polygon_contour_iterator p2 = p1;
const double cos_thr = 0.8;
const double circle_segment_thr = 2.5;
// search for the first circle segment (where cos(a) > cos_thr)
double ls_inner = 0.0, ls_outer = 0.0;
unsigned long n_ls_inner = 0, n_ls_outer = 0;
if (! fallback) {
do {
++p2;
if (p2 == to) {
p2 = from;
}
db::edge<C> ep (*p0, *p1);
db::edge<C> e (*p1, *p2);
bool inner = (db::vprod_sign (ep, e) > 0);
double &ls = inner ? ls_inner : ls_outer;
unsigned long &n_ls = inner ? n_ls_inner : n_ls_outer;
if (db::sprod (ep, e) > cos_thr * e.double_length () * ep.double_length ()) {
ls += std::min (e.double_length (), ep.double_length ());
++n_ls;
}
p0 = p1;
p1 = p2;
} while (p0 != from);
if (n_ls_inner > 0) {
ls_inner /= n_ls_inner;
}
if (n_ls_outer > 0) {
ls_outer /= n_ls_outer;
}
}
bool found = false;
// search for the first circle segment (where cos(a) > cos_thr)
// or a long segment is followed by a short one or the curvature changes.
typename db::polygon<C>::polygon_contour_iterator pm1 = from;
p0 = pm1;
++p0;
if (p0 == to) {
p0 = from;
}
p1 = p0;
++p1;
if (p1 == to) {
p1 = from;
}
p2 = p1;
++p2;
if (p2 == to) {
p2 = from;
}
typename db::polygon<C>::polygon_contour_iterator p3 = p2;
do {
++p3;
if (p3 == to) {
p3 = from;
}
db::edge<C> em (*pm1, *p0);
db::edge<C> ep (*p0, *p1);
db::edge<C> e (*p1, *p2);
db::edge<C> en (*p2, *p3);
bool first_or_last = fallback || (e.double_length () > circle_segment_thr * ep.double_length () || ep.double_length () > circle_segment_thr * e.double_length ())
|| (db::vprod_sign (em, ep) * db::vprod_sign (ep, e) < 0 || db::vprod_sign (ep, e) * db::vprod_sign (e, en) < 0);
if (first_or_last && db::sprod (ep, e) > cos_thr * e.double_length () * ep.double_length ()) {
double ls = (db::vprod_sign (ep, e) > 0 ? ls_inner : ls_outer);
if (! fallback && ((e.double_length () < circle_segment_thr * ls && ep.double_length () > circle_segment_thr * ls)
|| db::vprod_sign (em, ep) * db::vprod_sign (ep, e) < 0)) {
found = true;
break;
} else if (fallback && (ep.dx () == 0 || ep.dy () == 0)) {
found = true;
break;
}
}
pm1 = p0;
p0 = p1;
p1 = p2;
p2 = p3;
} while (pm1 != from);
if (! found) {
return false;
}
// create a list of new points without the rounded corners and compute rounding radii
if (new_pts) {
new_pts->clear ();
}
typename db::polygon<C>::polygon_contour_iterator pfirst = p0;
bool in_corner = false;
double ls_corner = 0.0;
db::edge<C> elast;
typename db::polygon<C>::polygon_contour_iterator plast;
double asum = 0.0;
unsigned int nseg = 0;
double rxi_sum = 0.0;
double rxo_sum = 0.0;
double da_sum = 0.0;
int n_corners = 0;
int ni_corners = 0;
int no_corners = 0;
p3 = p2;
do {
++p3;
if (p3 == to) {
p3 = from;
}
db::edge<C> em (*pm1, *p0);
db::edge<C> ep (*p0, *p1);
db::edge<C> e (*p1, *p2);
db::edge<C> en (*p2, *p3);
// Heuristic detection of a new circle segment:
// In fallback mode vertical or horizontal edges separate circle segments.
// In non-fallback mode either a long edge followed by a short one indicates the beginning of a new
// circle segment or a circle segment is detected when the curvature changes.
// The latter case detects situations where two circle segments directly attach to each other
// with different bending direction.
bool first_or_last = fallback || (e.double_length () > circle_segment_thr * ep.double_length () || ep.double_length () > circle_segment_thr * e.double_length ())
|| (db::vprod_sign (em, ep) * db::vprod_sign (ep, e) < 0 || db::vprod_sign (ep, e) * db::vprod_sign (e, en) < 0);
if (db::sprod (ep, e) > cos_thr * e.double_length () * ep.double_length ()) {
double ls = db::vprod_sign (ep, e) > 0 ? ls_inner : ls_outer;
if ((! fallback && first_or_last && ((e.double_length () < circle_segment_thr * ls && ep.double_length () > circle_segment_thr * ls) || db::vprod_sign (em, ep) * db::vprod_sign (ep, e) < 0)) ||
(fallback && (ep.dx () == 0 || ep.dy () == 0))) {
if (! in_corner) {
elast = ep;
plast = p1;
asum = db::vprod (*p1 - db::point<C> (), *p2 - db::point<C> ());
nseg = 1;
ls_corner = ls;
}
in_corner = true;
} else if ((!fallback && first_or_last && ((e.double_length () > circle_segment_thr * ls_corner && ep.double_length () < circle_segment_thr * ls_corner) || db::vprod_sign (ep, e) * db::vprod_sign (e, en) < 0)) ||
(fallback && (e.dx () == 0 || e.dy () == 0))) {
if (in_corner) {
std::pair<bool, db::point<C> > cp = elast.cut_point (e);
if (! cp.first) {
// We have a full 180 degree bend without a stop (actually two corners).
// Use the segment in between that is perpendicular to the start and end segment as stop edge.
typename db::polygon<C>::polygon_contour_iterator pp1 = plast;
typename db::polygon<C>::polygon_contour_iterator pp2 = pp1;
double asum_part = 0.0;
unsigned int nseg_part = 0;
while (pp1 != p1) {
++pp2;
if (pp2 == to) {
pp2 = from;
}
e = db::edge<C> (*pp1, *pp2);
if (db::sprod_sign (elast, e) == 0) {
break;
}
asum_part += db::vprod (*pp1 - db::point<C> (), *pp2 - db::point<C> ());
++nseg_part;
pp1 = pp2;
}
cp = elast.cut_point (e);
if (! cp.first) {
return false;
}
if (new_pts) {
new_pts->push_back (cp.second);
}
++nseg_part;
asum -= asum_part;
asum -= db::vprod (e.p1 () - db::point<C> (), e.p2 () - db::point<C> ());
nseg -= nseg_part;
asum_part += db::vprod (cp.second - db::point<C> (), elast.p2 () - db::point<C> ());
asum_part += db::vprod (*pp1 - db::point<C> (), cp.second - db::point<C> ());
double sin_atot = db::vprod (elast, e);
double cos_atot = db::sprod (elast, e);
double atot = fabs (atan2 (sin_atot, cos_atot));
double rx = sqrt (fabs (asum_part) * 0.5 / (tan (atot * 0.5) - tan (atot * 0.5 / nseg_part) * nseg_part));
double da = atot / nseg_part;
if (sin_atot > 0.0) {
rxi_sum += rx;
++ni_corners;
} else {
rxo_sum += rx;
++no_corners;
}
da_sum += da;
++n_corners;
elast = e;
e = db::edge<C> (*p1, *p2);
cp = elast.cut_point (e);
if (! cp.first) {
return false;
}
}
if (new_pts) {
new_pts->push_back (cp.second);
}
asum += db::vprod (cp.second - db::point<C> (), elast.p2 () - db::point<C> ());
asum += db::vprod (*p1 - db::point<C> (), cp.second - db::point<C> ());
++nseg;
double sin_atot = db::vprod (elast, e);
double cos_atot = db::sprod (elast, e);
double atot = fabs (atan2 (sin_atot, cos_atot));
double rx = sqrt (fabs (asum) * 0.5 / (tan (atot * 0.5) - tan (atot * 0.5 / nseg) * nseg));
double da = atot / nseg;
if (sin_atot > 0.0) {
rxi_sum += rx;
++ni_corners;
} else {
rxo_sum += rx;
++no_corners;
}
da_sum += da;
++n_corners;
}
in_corner = false;
} else if (in_corner) {
asum += db::vprod (*p1 - db::point<C> (), *p2 - db::point<C> ());
++nseg;
} else {
if (new_pts) {
new_pts->push_back (*p1);
}
}
} else {
if (new_pts) {
new_pts->push_back (*p1);
}
}
pm1 = p0;
p0 = p1;
p1 = p2;
p2 = p3;
} while (p0 != pfirst);
if (n_corners < 2) {
return false;
} else {
n = (unsigned int) floor (2.0 * M_PI / (da_sum / n_corners) + 0.5);
if (ni_corners > 0) {
rinner = floor ((rxi_sum / ni_corners * 0.5) + 0.5) * 2;
}
if (no_corners > 0) {
router = floor ((rxo_sum / no_corners * 0.5) + 0.5) * 2;
}
return true;
}
}
bool
extract_rad_from_contour (db::Polygon::polygon_contour_iterator from, db::Polygon::polygon_contour_iterator to, double &rinner, double &router, unsigned int &n, std::vector <db::Point> *new_pts, bool fallback)
{
return do_extract_rad_from_contour (from, to, rinner, router, n, new_pts, fallback);
}
bool
extract_rad_from_contour (db::DPolygon::polygon_contour_iterator from, db::DPolygon::polygon_contour_iterator to, double &rinner, double &router, unsigned int &n, std::vector <db::DPoint> *new_pts, bool fallback)
{
return do_extract_rad_from_contour (from, to, rinner, router, n, new_pts, fallback);
}
template <class C>
static bool
do_extract_rad (const db::polygon<C> &polygon, double &rinner, double &router, unsigned int &n, db::polygon<C> *new_polygon)
{
if (new_polygon) {
std::vector <db::point<C> > new_pts;
if (! do_extract_rad_from_contour (polygon.begin_hull (), polygon.end_hull (), rinner, router, n, &new_pts, false) &&
! do_extract_rad_from_contour (polygon.begin_hull (), polygon.end_hull (), rinner, router, n, &new_pts, true)) {
// no radius found
return false;
} else {
new_polygon->assign_hull (new_pts.begin (), new_pts.end ());
}
for (unsigned int h = 0; h < polygon.holes (); ++h) {
new_pts.clear ();
if (! do_extract_rad_from_contour (polygon.begin_hole (h), polygon.end_hole (h), rinner, router, n, &new_pts, false) &&
! do_extract_rad_from_contour (polygon.begin_hole (h), polygon.end_hole (h), rinner, router, n, &new_pts, true)) {
// no radius found
return false;
} else {
new_polygon->insert_hole (new_pts.begin (), new_pts.end ());
}
}
} else {
if (! do_extract_rad_from_contour (polygon.begin_hull (), polygon.end_hull (), rinner, router, n, (std::vector<db::point<C> > *) 0, false)) {
if (! do_extract_rad_from_contour (polygon.begin_hull (), polygon.end_hull (), rinner, router, n, (std::vector<db::point<C> > *) 0, true)) {
return false;
}
}
for (unsigned int h = 0; h < polygon.holes (); ++h) {
if (! do_extract_rad_from_contour (polygon.begin_hole (h), polygon.end_hole (h), rinner, router, n, (std::vector<db::point<C> > *) 0, false)) {
if (! do_extract_rad_from_contour (polygon.begin_hole (h), polygon.end_hole (h), rinner, router, n, (std::vector<db::point<C> > *) 0, true)) {
return false;
}
}
}
}
return true;
}
bool
extract_rad (const db::Polygon &polygon, double &rinner, double &router, unsigned int &n, db::Polygon *new_polygon)
{
return do_extract_rad (polygon, rinner, router, n, new_polygon);
}
bool
extract_rad (const db::DPolygon &polygon, double &rinner, double &router, unsigned int &n, db::DPolygon *new_polygon)
{
return do_extract_rad (polygon, rinner, router, n, new_polygon);
}
template <class C>
static void
do_compute_rounded_contour (typename db::polygon<C>::polygon_contour_iterator from, typename db::polygon<C>::polygon_contour_iterator to, std::vector <db::point<C> > &new_pts, double rinner, double router, unsigned int n)
{
std::vector<db::point<C> > points;
// collect the points into a vector
if (from != to) {
typename db::polygon<C>::polygon_contour_iterator p0 = from;
typename db::polygon<C>::polygon_contour_iterator p1 = p0;
++p1;
if (p1 == to) {
p1 = from;
}
typename db::polygon<C>::polygon_contour_iterator p2 = p1;
do {
++p2;
if (p2 == to) {
p2 = from;
}
if (! db::edge<C> (*p0, *p1).parallel (db::edge<C> (*p1, *p2))) {
points.push_back (*p1);
}
p0 = p1;
p1 = p2;
} while (p0 != from);
}
// compute the radii and segment length
std::vector<double> rad (points.size (), 0.0);
std::vector<double> seg (points.size (), 0.0);
for (size_t i = 0; i < points.size (); ++i) {
db::point<C> p0 = points [(i + points.size () - 1) % points.size ()];
db::point<C> p1 = points [i];
db::point<C> p2 = points [(i + points.size () + 1) % points.size ()];
db::DVector e1 = (db::DPoint (p1) - db::DPoint (p0)) * (1.0 / p0.double_distance (p1));
db::DVector e2 = (db::DPoint (p2) - db::DPoint (p1)) * (1.0 / p1.double_distance (p2));
double sin_a = db::vprod (e1, e2);
double cos_a = db::sprod (e1, e2);
double a = fabs (atan2 (sin_a, cos_a));
double r = sin_a > 0.0 ? rinner : router;
double s = r * fabs (sin (a * 0.5) / cos (a * 0.5));
rad [i] = r;
seg [i] = s;
}
// compute the rounded points
for (size_t i = 0; i < points.size (); ++i) {
db::point<C> p0 = points [(i + points.size () - 1) % points.size ()];
db::point<C> p1 = points [i];
db::point<C> p2 = points [(i + points.size () + 1) % points.size ()];
db::DVector e1 = (db::DPoint (p1) - db::DPoint (p0)) * (1.0 / p0.double_distance (p1));
db::DVector e2 = (db::DPoint (p2) - db::DPoint (p1)) * (1.0 / p1.double_distance (p2));
double sin_a = db::vprod (e1, e2);
double cos_a = db::sprod (e1, e2);
double a = fabs (atan2 (sin_a, cos_a));
double s0 = seg [(i + points.size () - 1) % points.size ()];
double s1 = seg [i];
double s2 = seg [(i + points.size () + 1) % points.size ()];
double f0 = std::min (1.0, p0.double_distance (p1) / (s0 + s1));
double f1 = std::min (1.0, p1.double_distance (p2) / (s1 + s2));
double r = std::min (f0, f1) * rad [i];
if (r > 0.0) {
db::DPoint q0 = db::DPoint (p1) - e1 * (tan (a * 0.5) * r);
db::DVector n1;
if (sin_a > 0) {
n1 = db::DVector (e1.y (), -e1.x ());
} else {
n1 = db::DVector (-e1.y (), e1.x ());
}
db::DPoint pr = q0 - n1 * r;
double ares = (2.0 * M_PI) / double (n);
unsigned int nseg = (unsigned int) floor (a / ares + 0.5);
if (nseg == 0) {
new_pts.push_back (p1);
} else {
double da = a / floor (a / ares + 0.5);
for (double aa = 0.0; aa < a - 1e-6; aa += da) {
db::DPoint q1 = pr + n1 * (r * cos (aa + da)) + e1 * (r * sin (aa + da));
// do a interpolation by computing the crossing point of the tangents of the
// circle at a and a+da. This scheme guarantees a low distortion of the original
// polygon and enables reverting back to the original polygon to some degree.
db::DPoint qm = q0 + (q1 - q0) * 0.5;
db::DPoint q = qm + (qm - pr) * (q0.sq_distance (qm) / pr.sq_distance (qm));
new_pts.push_back (db::point<C> (q));
q0 = q1;
}
}
} else {
new_pts.push_back (p1);
}
}
}
void
compute_rounded_contour (db::Polygon::polygon_contour_iterator from, db::Polygon::polygon_contour_iterator to, std::vector <db::Point> &new_pts, double rinner, double router, unsigned int n)
{
do_compute_rounded_contour (from, to, new_pts, rinner, router, n);
}
void
compute_rounded_contour (db::DPolygon::polygon_contour_iterator from, db::DPolygon::polygon_contour_iterator to, std::vector <db::DPoint> &new_pts, double rinner, double router, unsigned int n)
{
do_compute_rounded_contour (from, to, new_pts, rinner, router, n);
}
template <class C>
static db::polygon<C>
do_compute_rounded (const db::polygon<C> &polygon, double rinner, double router, unsigned int n)
{
db::polygon<C> new_poly;
std::vector <db::point<C> > new_pts;
compute_rounded_contour (polygon.begin_hull (), polygon.end_hull (), new_pts, rinner, router, n);
new_poly.assign_hull (new_pts.begin (), new_pts.end (), false /*don't compress*/);
for (unsigned int h = 0; h < polygon.holes (); ++h) {
new_pts.clear ();
compute_rounded_contour (polygon.begin_hole (h), polygon.end_hole (h), new_pts, rinner, router, n);
new_poly.insert_hole (new_pts.begin (), new_pts.end (), false /*don't compress*/);
}
return new_poly;
}
db::Polygon
compute_rounded (const db::Polygon &polygon, double rinner, double router, unsigned int n)
{
return do_compute_rounded (polygon, rinner, router, n);
}
db::DPolygon
compute_rounded (const db::DPolygon &polygon, double rinner, double router, unsigned int n)
{
return do_compute_rounded (polygon, rinner, router, n);
}
// -------------------------------------------------------------------------
// Implementation of AreaMap
AreaMap::AreaMap ()
: m_nx (0), m_ny (0)
{
mp_av = 0;
}
AreaMap::AreaMap (const db::Point &p0, const db::Vector &d, size_t nx, size_t ny)
: m_p0 (p0), m_d (d), m_nx (nx), m_ny (ny)
{
mp_av = new area_type [nx * ny];
clear ();
}
AreaMap::~AreaMap ()
{
if (mp_av) {
delete[] mp_av;
}
mp_av = 0;
}
void
AreaMap::reinitialize (const db::Point &p0, const db::Vector &d, size_t nx, size_t ny)
{
m_p0 = p0;
m_d = d;
m_nx = nx;
m_ny = ny;
if (mp_av) {
delete mp_av;
}
mp_av = new area_type [nx * ny];
clear ();
}
void
AreaMap::clear ()
{
if (mp_av) {
area_type *a = mp_av;
for (size_t n = m_nx * m_ny; n > 0; --n) {
*a++ = 0;
}
}
}
void
AreaMap::swap (AreaMap &other)
{
std::swap (m_p0, other.m_p0);
std::swap (m_d, other.m_d);
std::swap (m_nx, other.m_nx);
std::swap (m_ny, other.m_ny);
std::swap (mp_av, other.mp_av);
}
AreaMap::area_type
AreaMap::total_area () const
{
area_type asum = 0;
if (mp_av) {
const area_type *a = mp_av;
for (size_t n = m_nx * m_ny; n > 0; --n) {
asum += *a++;
}
}
return asum;
}
// -------------------------------------------------------------------------
// Implementation of rasterize
void
rasterize (const db::Polygon &polygon, db::AreaMap &am)
{
typedef db::AreaMap::area_type area_type;
db::Box box = am.bbox ();
db::Box pbox = polygon.box ();
// check if the polygon overlaps the rasterization area. Otherwise, we simply do nothing.
if (! pbox.overlaps (box)) {
return;
}
db::Coord ymin = box.bottom (), ymax = box.top ();
db::Coord dy = am.d ().y (), dx = am.d ().x ();
db::Coord y0 = am.p0 ().y (), x0 = am.p0 ().x ();
size_t ny = am.ny (), nx = am.nx ();
size_t iy0 = std::min (ny, size_t (std::max (db::Coord (0), (pbox.bottom () - am.p0 ().y ()) / am.d ().y ())));
size_t iy1 = std::min (ny, size_t (std::max (db::Coord (0), (pbox.top () - am.p0 ().y () + am.d ().y () - 1) / am.d ().y ())));
size_t ix0 = std::min (nx, size_t (std::max (db::Coord (0), (pbox.left () - am.p0 ().x ()) / am.d ().x ())));
size_t ix1 = std::min (nx, size_t (std::max (db::Coord (0), (pbox.right () - am.p0 ().x () + am.d ().x () - 1) / am.d ().x ())));
// no scanning required (i.e. degenerated polygon) -> do nothing
if (iy0 == iy1 || ix0 == ix1) {
return;
}
// collect edges
size_t n = 0;
for (db::Polygon::polygon_edge_iterator e = polygon.begin_edge (); ! e.at_end (); ++e) {
if ((*e).dy () != 0 && db::edge_ymax (*e) > ymin && db::edge_ymin (*e) < ymax) {
++n;
}
}
std::vector <db::Edge> edges;
edges.reserve (n);
for (db::Polygon::polygon_edge_iterator e = polygon.begin_edge (); ! e.at_end (); ++e) {
if ((*e).dy () != 0 && db::edge_ymax (*e) > ymin && db::edge_ymin (*e) < ymax) {
edges.push_back (*e);
}
}
// sort edges
std::sort (edges.begin (), edges.end (), db::edge_ymin_compare<db::Coord> ());
std::vector <db::Edge>::iterator c = edges.begin ();
db::Coord y = y0 + dy * db::Coord (iy0);
while (c != edges.end () && db::edge_ymax (*c) <= y) {
++c;
}
std::vector <db::Edge>::iterator f = c;
for (size_t iy = iy0; iy < iy1; ++iy) {
db::Coord yy = y + dy;
while (f != edges.end () && db::edge_ymin (*f) < yy) {
++f;
}
std::sort (c, f, db::edge_xmin_compare <db::Coord> ());
db::Coord x = x0 + dx * db::Coord (ix0);
db::Coord xl = pbox.left ();
area_type a = 0;
std::vector <db::Edge>::iterator cc = c;
while (cc != edges.end () && db::edge_xmax (*cc) <= x) {
db::Coord y1 = std::max (y, std::min (yy, cc->p1 ().y ()));
db::Coord y2 = std::max (y, std::min (yy, cc->p2 ().y ()));
a += area_type (dx) * area_type (y2 - y1);
++cc;
}
std::vector <db::Edge>::iterator ff = cc;
for (size_t ix = ix0; ix < ix1; ++ix) {
db::Coord xx = x + dx;
// TODO: edge_xmin_at_interval(y, yy) and edge_xmax.. would be more efficient in the
// all-angle case. However, it is crucial that the edge clipping produces
// connected edge segments and it is questionable whether the at_interval
// functions produce a sorting/filter criterion compatible with the clip.
while (ff != f && db::edge_xmin (*ff) < xx) {
++ff;
}
if (xl < x) {
// consider all edges or parts of those left of the first cell
db::Box left (xl, y, x, yy);
for (std::vector <db::Edge>::iterator e = cc; e != ff; ++e) {
std::pair<bool, db::Edge> ec = e->clipped (left);
if (ec.first && db::edge_xmin (ec.second) < x) {
a += area_type (ec.second.dy ()) * area_type (dx);
}
}
}
db::Box cell (x, y, xx, yy);
area_type aa = a;
for (std::vector <db::Edge>::iterator e = cc; e != ff; ++e) {
std::pair<bool, db::Edge> ec = e->clipped (cell);
if (ec.first && db::edge_xmin (ec.second) < xx) {
aa += area_type (ec.second.dy ()) * area_type (2 * xx - (ec.second.p2 ().x () + ec.second.p1 ().x ())) / 2;
a += area_type (ec.second.dy ()) * area_type (dx);
}
}
am.get (ix, iy) += aa;
x = xx;
xl = xx;
for (std::vector <db::Edge>::iterator ccx = cc; ccx != ff; ++ccx) {
if (db::edge_xmax (*ccx) <= x) {
std::swap (*ccx, *cc);
++cc;
}
}
}
y = yy;
for (std::vector <db::Edge>::iterator cx = c; cx != f; ++cx) {
if (db::edge_ymax (*cx) <= y) {
std::swap (*cx, *c);
++c;
}
}
}
}
// -------------------------------------------------------------------------
// Implementation of minkowsky_sum
/**
* @brief A helper class that produces edges into an EdgeProcessor from a sequence of points
*/
class EdgeInputIterator
{
public:
EdgeInputIterator (db::EdgeProcessor &ep, bool inverse = false)
: m_last_set (false), m_last (), mp_ep (&ep), m_inverse (inverse)
{ }
~EdgeInputIterator ()
{
// close the polygon
if (m_last_set && m_last != m_first) {
if (!m_inverse) {
mp_ep->insert (db::Edge (m_last, m_first));
} else {
mp_ep->insert (db::Edge (m_first, m_last));
}
m_last_set = false;
}
mp_ep = 0;
}
void operator+= (const db::Point &p)
{
if (m_last_set) {
if (!m_inverse) {
mp_ep->insert (db::Edge (m_last, p));
} else {
mp_ep->insert (db::Edge (p, m_last));
}
} else {
m_first = p;
}
m_last = p;
m_last_set = true;
}
private:
bool m_last_set;
db::Point m_last, m_first;
db::EdgeProcessor *mp_ep;
bool m_inverse;
};
/**
* @brief Produce edges for the partial Minkowsky sum of an edge with an input polygon
*/
static void
ms_production (const db::Polygon &a, const db::Point &p1, const db::Point &p2, db::EdgeProcessor &ep)
{
double d12 = p2.double_distance (p1);
db::DPoint d (-double (p2.y () - p1.y ()) / d12, double (p2.x () - p1.x ()) / d12);
db::EdgeInputIterator e (ep);
db::Polygon::polygon_contour_iterator ci = a.begin_hull ();
db::Polygon::polygon_contour_iterator cf = a.end_hull ();
db::Polygon::polygon_contour_iterator cmin = cf, cmax = cf;
double pmin = 0.0, pmax = 0.0;
// Look for the points in the contour bounding the partial sum perpendicular to the edge
for (db::Polygon::polygon_contour_iterator c = ci; c != cf; ++c) {
double p = (*c).x() * d.x () + (*c).y () * d.y ();
if (cmin == cf || pmin > p) {
pmin = p;
cmin = c;
}
if (cmax == cf || pmax < p) {
pmax = p;
cmax = c;
}
}
tl_assert (cmin != cf && cmax != cf);
db::Polygon::polygon_contour_iterator c = cmin;
db::Polygon::polygon_contour_iterator cc = cf;
db::Polygon::polygon_contour_iterator cl = cf;
bool pcc_set = false;
double pcc = 0.0;
while (true) {
double pc = (*c).x() * d.x () + (*c).y () * d.y ();
// detect inversion due to a convex pattern and create a cover polygon for that case
if (pcc_set) {
cc = c;
if (cc == ci) {
cc = cf;
}
--cc;
if (pcc > pc + 1e-6) {
if (cl == cf) {
cl = cc;
}
} else if (cl != cf) {
EdgeInputIterator ee (ep, true);
// create the cover polygon
db::Polygon::polygon_contour_iterator k = cl;
while (k != cc) {
ee += p1 + (*k - db::Point ());
if (++k == cf) {
k = ci;
}
}
ee += p1 + (*k - db::Point ());
while (k != cl) {
ee += p2 + (*k - db::Point ());
if (k == ci) {
k = cf;
}
--k;
}
ee += p2 + (*cl - db::Point ());
cl = cf;
}
}
// produce a new edge
e += p1 + (*c - db::Point ());
if (c == cmax) {
break;
}
if (++c == cf) {
c = ci;
}
pcc = pc;
pcc_set = true;
}
cl = cf;
pcc_set = false;
while (true) {
double pc = (*c).x() * d.x () + (*c).y () * d.y ();
// detect inversion due to a convex pattern and create a cover polygon for that case
if (pcc_set) {
cc = c;
if (cc == ci) {
cc = cf;
}
--cc;
pcc = (*cc).x() * d.x () + (*cc).y () * d.y ();
if (pcc < pc - 1e-6) {
if (cl == cf) {
cl = cc;
}
} else if (cl != cf) {
EdgeInputIterator ee (ep, true);
// create the cover polygon
db::Polygon::polygon_contour_iterator k = cl;
// create the cover polygon
while (k != cc) {
ee += p2 + (*k - db::Point ());
if (++k == cf) {
k = ci;
}
}
ee += p2 + (*cc - db::Point ());
while (k != cl) {
ee += p1 + (*k - db::Point ());
if (k == ci) {
k = cf;
}
--k;
}
ee += p1 + (*cl - db::Point ());
cl = cf;
}
}
e += p2 + (*c - db::Point ());
if (c == cmin) {
break;
}
if (++c == cf) {
c = ci;
}
pcc = pc;
pcc_set = true;
}
}
static db::Polygon
ms_extraction (db::EdgeProcessor &ep, bool resolve_holes)
{
db::SimpleMerge op (-1);
std::vector <db::Polygon> polygons;
db::PolygonContainer pc (polygons);
db::PolygonGenerator out (pc, resolve_holes, false);
ep.process (out, op);
if (polygons.empty ()) {
return db::Polygon ();
} else {
tl_assert (polygons.size () == 1);
return polygons [0];
}
}
static db::Polygon
do_minkowsky_sum (const db::Polygon &a, const db::Edge &b, bool resolve_holes)
{
db::EdgeProcessor ep;
db::ms_production (a, b.p1 (), b.p2 (), ep);
return db::ms_extraction (ep, resolve_holes);
}
db::Polygon
minkowsky_sum (const db::Polygon &a, const db::Edge &b, bool rh)
{
if (a.holes () > 0) {
return do_minkowsky_sum (db::resolve_holes (a), b, rh);
} else {
return do_minkowsky_sum (a, b, rh);
}
}
static db::Polygon
do_minkowsky_sum (const db::Polygon &a, const db::Polygon &b, bool resolve_holes)
{
tl_assert (a.begin_hull () != a.end_hull ());
db::Vector p0 = *a.begin_hull () - db::Point ();
db::EdgeProcessor ep;
for (db::Polygon::polygon_edge_iterator e = b.begin_edge (); ! e.at_end (); ++e) {
ep.insert (db::Edge ((*e).p1 () + p0, (*e).p2 () + p0));
db::ms_production (a, (*e).p1 (), (*e).p2 (), ep);
}
return db::ms_extraction (ep, resolve_holes);
}
db::Polygon
minkowsky_sum (const db::Polygon &a, const db::Polygon &b, bool rh)
{
if (a.holes () > 0) {
return do_minkowsky_sum (db::resolve_holes (a), b, rh);
} else {
return do_minkowsky_sum (a, b, rh);
}
}
static db::Polygon
do_minkowsky_sum (const db::Polygon &a, const db::Box &b, bool resolve_holes)
{
return minkowsky_sum (a, db::Polygon (b), resolve_holes);
}
db::Polygon
minkowsky_sum (const db::Polygon &a, const db::Box &b, bool rh)
{
if (a.holes () > 0) {
return do_minkowsky_sum (db::resolve_holes (a), b, rh);
} else {
return do_minkowsky_sum (a, b, rh);
}
}
static db::Polygon
do_minkowsky_sum (const db::Polygon &a, const std::vector<db::Point> &c, bool resolve_holes)
{
db::EdgeProcessor ep;
for (size_t i = 1; i < c.size (); ++i) {
db::ms_production (a, c[i - 1], c[i], ep);
}
return db::ms_extraction (ep, resolve_holes);
}
db::Polygon
minkowsky_sum (const db::Polygon &a, const std::vector<db::Point> &c, bool rh)
{
if (a.holes () > 0) {
return do_minkowsky_sum (db::resolve_holes (a), c, rh);
} else {
return do_minkowsky_sum (a, c, rh);
}
}
// -------------------------------------------------------------------------
// Implementation of hole resolution and polygon to simple polygon conversion
db::Polygon
resolve_holes (const db::Polygon &p)
{
db::EdgeProcessor ep;
ep.insert_sequence (p.begin_edge ());
std::vector<db::Polygon> polygons;
db::PolygonContainer pc (polygons);
db::PolygonGenerator out (pc, true /*resolve holes*/, false /*max coherence to get one polygon*/);
db::SimpleMerge op;
ep.process (out, op);
if (polygons.empty ()) {
return db::Polygon ();
} else {
tl_assert (polygons.size () == 1);
return polygons [0];
}
}
db::Polygon
simple_polygon_to_polygon (const db::SimplePolygon &sp)
{
db::Polygon p;
p.assign_hull (sp.begin_hull (), sp.end_hull ());
return p;
}
db::SimplePolygon
polygon_to_simple_polygon (const db::Polygon &p)
{
if (p.holes () > 0) {
db::Polygon pp = resolve_holes (p);
db::SimplePolygon sp;
sp.assign_hull (pp.begin_hull (), pp.end_hull ());
return sp;
} else {
db::SimplePolygon sp;
sp.assign_hull (p.begin_hull (), p.end_hull ());
return sp;
}
}
static void decompose_convex_helper (int depth, PreferredOrientation po, const db::SimplePolygon &sp, SimplePolygonSink &sink)
{
size_t n = sp.hull ().size ();
if (n < 4 || depth <= 0) {
if (n > 2) {
sink.put (sp);
}
return;
}
db::Box bbox = sp.box ();
db::coord_traits<db::Coord>::area_type atot = 0;
for (size_t i = 0; i < n; ++i) {
db::Edge ep (sp.hull ()[(i + n - 1) % n], sp.hull ()[i]);
atot += db::vprod (ep.p2 () - db::Point (), ep.p1 () - db::Point ());
}
std::set<db::Point> skipped;
while (true) {
// Look for the convex corner closed to the median
size_t imed = std::numeric_limits<size_t>::max ();
db::Coord dmin = 0;
for (size_t i = 0; i < n; ++i) {
db::Edge ep (sp.hull ()[(i + n - 1) % n], sp.hull ()[i]);
db::Edge ec (sp.hull ()[i], sp.hull ()[(i + 1) % n]);
if (db::vprod_sign (ep, ec) > 0 && skipped.find (ep.p2 ()) == skipped.end ()) {
db::Vector v = sp.hull ()[i] - bbox.center ();
db::Coord d = std::min (std::abs (v.x ()), std::abs (v.y ()));
if (imed == std::numeric_limits<size_t>::max () || d < dmin) {
imed = i;
dmin = d;
}
}
}
// is convex already
if (imed == std::numeric_limits<size_t>::max ()) {
if (! skipped.empty ()) {
tl::warn << "sp=" << sp << tl::endl << "po=" << int(po);
tl_assert (false);
}
sink.put (sp);
return;
}
db::Point p (sp.hull ()[imed]);
db::Edge ep (sp.hull ()[(imed + n - 1) % n], p);
db::Edge ec (p, sp.hull ()[(imed + 1) % n]);
// convex corner
std::set<db::Vector> cuts;
db::Vector pv = db::Vector (ep.dy (), -ep.dx ()) + db::Vector (ec.dy (), -ec.dx ());
bool ortho = ((ep.dx () == 0 || ep.dy () == 0) && (ec.dx () == 0 || ec.dy () == 0));
if (po == PO_any || po == PO_horizontal || po == PO_htrapezoids || (po == PO_vtrapezoids && ortho)) {
if (pv.x () >= 0) {
cuts.insert (db::Vector (1, 0));
}
if (pv.x () <= 0) {
cuts.insert (db::Vector (-1, 0));
}
}
if (po == PO_any || po == PO_vertical || po == PO_vtrapezoids || (po == PO_htrapezoids && ortho)) {
if (pv.y () >= 0) {
cuts.insert (db::Vector (0, 1));
}
if (pv.y () <= 0) {
cuts.insert (db::Vector (0, -1));
}
}
int cut_rating = 0;
size_t jmin = 0;
db::Point xmin;
db::coord_traits<db::Coord>::area_type acutoff = 0;
for (std::set<db::Vector>::const_iterator c = cuts.begin (); c != cuts.end (); ++c) {
db::Vector nv = *c;
int cut_rating_inner = 0;
size_t jmin_inner = 0;
db::Point xmin_inner;
db::coord_traits<db::Coord>::area_type acutoff_inner = 0;
db::coord_traits<db::Coord>::area_type asum = 0;
db::coord_traits<db::Coord>::area_type min_dist = std::numeric_limits<db::coord_traits<db::Coord>::area_type>::max ();
for (size_t j = 1; j != n - 1; ++j) {
db::Edge efc (sp.hull ()[(imed + j) % n], sp.hull ()[(imed + j + 1) % n]);
db::Edge efp (sp.hull ()[(imed + j + n - 1) % n], sp.hull ()[(imed + j) % n]);
asum += db::vprod (efp.p2 () - db::Point (), efp.p1 () - db::Point ());
std::pair <bool, db::Point> x;
if (db::vprod_sign (nv, efc.d ()) == 0) {
if (efc.side_of (p) == 0) {
db::coord_traits<db::Coord>::area_type d1 = db::sprod (efc.p1 () - p, nv);
db::coord_traits<db::Coord>::area_type d2 = db::sprod (efc.p2 () - p, nv);
if (d1 <= 0 && d2 >= 0) {
x = std::make_pair (true, p);
} else if (d2 <= 0 && d1 >= 0) {
x = std::make_pair (true, p);
} else if (d1 >= 0 && d2 >= 0) {
if (d1 < d2) {
x = std::make_pair (true, efc.p1 ());
} else {
x = std::make_pair (true, efc.p2 ());
}
}
}
} else {
x = db::Edge (p, nv).crossed_by_point (efc);
}
if (x.first && x.second != efc.p2 ()) {
db::coord_traits<db::Coord>::area_type dist = db::sprod (x.second - p, nv);
if (dist >= 0 && dist < min_dist) {
// a will be the area of the half we cut off
db::coord_traits<db::Coord>::area_type a = asum +
db::vprod (x.second - db::Point (), efp.p2 () - db::Point ()) +
db::vprod (p - db::Point (), x.second - db::Point ());
// due to rounding, the new cut point will modify the total
// area. We compute the new total area now:
db::coord_traits<db::Coord>::area_type atot_eff = atot +
db::vprod (efc.p2 () - db::Point (), x.second - db::Point ()) +
db::vprod (x.second - db::Point (), efc.p1 () - db::Point ()) +
db::vprod (efc.p1 () - db::Point (), efc.p2 () - db::Point ());
db::coord_traits<db::Coord>::area_type ac;
if (a > atot_eff / 2) {
ac = a - atot_eff / 2;
} else {
ac = atot_eff / 2 - a;
}
if (db::vprod_sign (nv, efc.d ()) <= 0 && a >= 0 && a <= atot_eff) {
// compute rating
int cr = 0;
if (x.second == efc.p1 ()) {
if (db::vprod (efc, efp) < 0) {
cr = 2; // cut terminates at another concave corner
} else {
cr = 1; // cut terminates at a convex corner
}
}
jmin_inner = j;
cut_rating_inner = cr;
xmin_inner = x.second;
acutoff_inner = ac;
min_dist = dist;
} else if (db::vprod_sign (nv, efc.d ()) < 0) {
min_dist = dist;
jmin_inner = 0;
}
}
}
}
if (jmin_inner > 0 && (!jmin || cut_rating_inner > cut_rating || (cut_rating_inner == cut_rating && acutoff_inner < acutoff))) {
jmin = jmin_inner;
cut_rating = cut_rating_inner;
xmin = xmin_inner;
acutoff = acutoff_inner;
}
}
if (jmin > 0) {
std::vector<db::Point> pts;
pts.reserve (n);
db::SimplePolygon sp_out;
for (size_t i = imed; i <= imed + jmin; ++i) {
pts.push_back (sp.hull ()[i % n]);
}
if (pts.back () != xmin) {
pts.push_back (xmin);
}
sp_out.assign_hull (pts.begin (), pts.end (), true, true);
decompose_convex_helper (depth - 1, po, sp_out, sink);
pts.clear ();
for (size_t i = imed + jmin + 1; i <= imed + n; ++i) {
pts.push_back (sp.hull ()[i % n]);
}
if (pts.front () != xmin) {
pts.push_back (xmin);
}
sp_out.assign_hull (pts.begin (), pts.end (), true, true);
decompose_convex_helper (depth - 1, po, sp_out, sink);
break;
} else {
// no decomposition found -> next try
skipped.insert (p);
}
}
}
class ConvexDecompositionFilter
: public db::SimplePolygonSink
{
public:
ConvexDecompositionFilter (PreferredOrientation po, bool swap_xy, db::SimplePolygonSink *out)
: mp_out (out), m_po (po), m_swap_xy (swap_xy)
{
// .. nothing yet ..
}
virtual void put (const db::SimplePolygon &polygon)
{
if (m_swap_xy) {
db::SimplePolygon p (polygon);
p.transform (db::FTrans (db::FTrans::m45));
decompose_convex_helper (std::numeric_limits<int>::max (), m_po, p, *mp_out);
} else {
decompose_convex_helper (std::numeric_limits<int>::max (), m_po, polygon, *mp_out);
}
}
private:
db::SimplePolygonSink *mp_out;
PreferredOrientation m_po;
bool m_swap_xy;
};
void
decompose_convex (const db::Polygon &p, PreferredOrientation po, db::SimplePolygonSink &sink)
{
if (p.is_box ()) {
sink.put (db::SimplePolygon (p.box ()));
} else {
// Because the hole resolution strategy favours horizontal lines we need to swap x and y
// for the PO_vertical and PO_vtrapezoids case
bool swap_xy = (po == PO_vertical || po == PO_vtrapezoids);
ConvexDecompositionFilter cd (po, swap_xy, &sink);
db::PolygonGenerator pg (cd, true /*min coherence*/);
// Does some pre-decomposition and avoids self-interacting polygons:
pg.open_contours (true);
db::EdgeProcessor ep;
if (swap_xy) {
for (db::Polygon::polygon_edge_iterator e = p.begin_edge (); ! e.at_end (); ++e) {
ep.insert ((*e).transformed (db::FTrans (db::FTrans::m45)));
}
} else {
ep.insert_sequence (p.begin_edge ());
}
db::SimpleMerge op;
ep.process (pg, op);
}
}
void
decompose_convex (const db::SimplePolygon &sp, PreferredOrientation po, db::SimplePolygonSink &sink)
{
if (sp.is_box ()) {
sink.put (sp);
} else {
decompose_convex_helper (std::numeric_limits<int>::max (), po, sp, sink);
}
}
template <class P>
bool
is_convex_helper (const P &p)
{
size_t n = p.hull ().size ();
if (n < 4) {
return true;
}
for (size_t i = 0; i < n; ++i) {
db::Edge ep (p.hull ()[(i + n - 1) % n], p.hull ()[i]);
db::Edge ec (p.hull ()[i], p.hull ()[(i + 1) % n]);
if (db::vprod_sign (ep, ec) > 0) {
return false;
}
}
return true;
}
bool
is_convex (const db::SimplePolygon &p)
{
return is_convex_helper (p);
}
bool
is_convex (const db::Polygon &p)
{
if (p.holes () > 0) {
return false;
} else {
return is_convex_helper (p);
}
}
static void
decompose_convex_to_trapezoids (const db::SimplePolygon &sp, bool horizontal, db::SimplePolygonSink &sink)
{
if (sp.hull ().size () < 3) {
return;
}
std::vector<db::Edge> edges;
edges.reserve (sp.hull ().size ());
for (db::SimplePolygon::polygon_edge_iterator e = sp.begin_edge (); ! e.at_end (); ++e) {
db::Edge ee = *e;
if (! horizontal) {
ee.transform (db::FTrans (db::FTrans::m45));
}
if (ee.dy () != 0) {
edges.push_back (ee);
}
}
std::sort (edges.begin (), edges.end (), db::edge_ymin_compare<db::Coord> ());
db::Coord y = db::edge_ymin (edges.front ());
std::vector<db::Edge>::iterator c = edges.begin ();
while (c != edges.end ()) {
std::vector<db::Edge>::iterator cc = c;
while (cc != edges.end () && db::edge_ymin (*cc) <= y) {
++cc;
}
// this condition will be fulfilled always if the input polygon is convex
tl_assert (cc - c == 2);
db::Coord x1 = db::edge_xaty (c[0], y);
db::Coord x2 = db::edge_xaty (c[1], y);
if (x1 > x2) {
std::swap (x1, x2);
}
db::Coord yy;
if (cc == edges.end ()) {
yy = db::edge_ymax (*c);
tl_assert (db::edge_ymax (c[1]) == db::edge_ymax (*c));
} else {
yy = db::edge_ymin (*cc);
}
db::Coord xx1 = db::edge_xaty (c[0], yy);
db::Coord xx2 = db::edge_xaty (c[1], yy);
if (xx1 > xx2) {
std::swap (xx1, xx2);
}
db::SimplePolygon sp;
if (x1 == x2) {
db::Point pts[] = {
db::Point (x1, y),
db::Point (xx1, yy),
db::Point (xx2, yy)
};
sp.assign_hull (&pts[0], &pts[sizeof (pts) / sizeof (pts[0])]);
} else if (xx1 == xx2) {
db::Point pts[] = {
db::Point (x1, y),
db::Point (xx1, yy),
db::Point (x2, y)
};
sp.assign_hull (&pts[0], &pts[sizeof (pts) / sizeof (pts[0])]);
} else {
db::Point pts[] = {
db::Point (x1, y),
db::Point (xx1, yy),
db::Point (xx2, yy),
db::Point (x2, y)
};
sp.assign_hull (&pts[0], &pts[sizeof (pts) / sizeof (pts[0])]);
}
if (! horizontal) {
sp.transform (db::FTrans (db::FTrans::m45));
}
sink.put (sp);
std::vector<db::Edge>::iterator c0 = c;
for (std::vector<db::Edge>::iterator i = c; i != cc; ++i) {
if (db::edge_ymax (*i) <= yy) {
if (c != i) {
std::swap (*c, *i);
}
++c;
}
}
y = yy;
tl_assert (c0 != c);
}
}
namespace {
struct TrapezoidConverter
: public db::SimplePolygonSink
{
TrapezoidConverter (bool horizontal, db::SimplePolygonSink *target)
: m_horizontal (horizontal), mp_target (target)
{
// .. nothing yet ..
}
virtual void put (const db::SimplePolygon &polygon)
{
decompose_convex_to_trapezoids (polygon, m_horizontal, *mp_target);
}
public:
bool m_horizontal;
db::SimplePolygonSink *mp_target;
};
}
void
decompose_trapezoids (const db::Polygon &p, TrapezoidDecompositionMode mode, db::SimplePolygonSink &sink)
{
if (mode == TD_htrapezoids || mode == TD_vtrapezoids) {
// Implementation uses convex decomposition and trapezoid decomposition
if (p.is_box ()) {
sink.put (db::SimplePolygon (p.box ()));
} else {
bool swap_xy = (mode == TD_vtrapezoids);
TrapezoidConverter trap_maker (mode == TD_htrapezoids, &sink);
ConvexDecompositionFilter cd (mode == TD_htrapezoids ? PO_htrapezoids : PO_vtrapezoids, swap_xy, &trap_maker);
db::PolygonGenerator pg (cd, true /*min coherence*/);
// Does some pre-decomposition and avoids self-interacting polygons:
pg.open_contours (true);
db::EdgeProcessor ep;
if (swap_xy) {
for (db::Polygon::polygon_edge_iterator e = p.begin_edge (); ! e.at_end (); ++e) {
ep.insert ((*e).transformed (db::FTrans (db::FTrans::m45)));
}
} else {
ep.insert_sequence (p.begin_edge ());
}
db::SimpleMerge op;
ep.process (pg, op);
}
} else {
// Implementation uses trapezoid generator
if (p.is_box ()) {
sink.put (db::SimplePolygon (p.box ()));
} else {
db::TrapezoidGenerator pg (sink);
db::EdgeProcessor ep;
db::SimpleMerge op;
ep.insert_sequence (p.begin_edge ());
ep.process (pg, op);
}
}
}
void
decompose_trapezoids (const db::SimplePolygon &sp, TrapezoidDecompositionMode mode, db::SimplePolygonSink &sink)
{
if (mode == TD_htrapezoids || mode == TD_vtrapezoids) {
if (sp.is_box ()) {
sink.put (sp);
} else {
TrapezoidConverter trap_maker (mode == TD_htrapezoids, &sink);
decompose_convex_helper (std::numeric_limits<int>::max (), mode == TD_htrapezoids ? PO_htrapezoids : PO_vtrapezoids, sp, trap_maker);
}
} else {
// This implementation uses trapezoid generator
if (sp.is_box ()) {
sink.put (db::SimplePolygon (sp.box ()));
} else {
db::TrapezoidGenerator pg (sink);
db::EdgeProcessor ep;
db::SimpleMerge op;
ep.insert_sequence (sp.begin_edge ());
ep.process (pg, op);
}
}
}
}