/* -*- mode: c++; c-basic-offset: 4 -*- */
#ifndef MPL_PATH_H
#define MPL_PATH_H
#include
#include
#include
#include
#include
#include
#include "agg_conv_contour.h"
#include "agg_conv_curve.h"
#include "agg_conv_stroke.h"
#include "agg_conv_transform.h"
#include "agg_trans_affine.h"
#include "path_converters.h"
#include "_backend_agg_basic_types.h"
const size_t NUM_VERTICES[] = { 1, 1, 1, 2, 3 };
struct XY
{
double x;
double y;
XY() : x(0), y(0) {}
XY(double x_, double y_) : x(x_), y(y_)
{
}
bool operator==(const XY& o)
{
return (x == o.x && y == o.y);
}
bool operator!=(const XY& o)
{
return (x != o.x || y != o.y);
}
};
typedef std::vector Polygon;
inline void
_finalize_polygon(std::vector &result, bool closed_only)
{
if (result.size() == 0) {
return;
}
Polygon &polygon = result.back();
/* Clean up the last polygon in the result. */
if (polygon.size() == 0) {
result.pop_back();
} else if (closed_only) {
if (polygon.size() < 3) {
result.pop_back();
} else if (polygon.front() != polygon.back()) {
polygon.push_back(polygon.front());
}
}
}
//
// The following function was found in the Agg 2.3 examples (interactive_polygon.cpp).
// It has been generalized to work on (possibly curved) polylines, rather than
// just polygons. The original comments have been kept intact.
// -- Michael Droettboom 2007-10-02
//
//======= Crossings Multiply algorithm of InsideTest ========================
//
// By Eric Haines, 3D/Eye Inc, erich@eye.com
//
// This version is usually somewhat faster than the original published in
// Graphics Gems IV; by turning the division for testing the X axis crossing
// into a tricky multiplication test this part of the test became faster,
// which had the additional effect of making the test for "both to left or
// both to right" a bit slower for triangles than simply computing the
// intersection each time. The main increase is in triangle testing speed,
// which was about 15% faster; all other polygon complexities were pretty much
// the same as before. On machines where division is very expensive (not the
// case on the HP 9000 series on which I tested) this test should be much
// faster overall than the old code. Your mileage may (in fact, will) vary,
// depending on the machine and the test data, but in general I believe this
// code is both shorter and faster. This test was inspired by unpublished
// Graphics Gems submitted by Joseph Samosky and Mark Haigh-Hutchinson.
// Related work by Samosky is in:
//
// Samosky, Joseph, "SectionView: A system for interactively specifying and
// visualizing sections through three-dimensional medical image data",
// M.S. Thesis, Department of Electrical Engineering and Computer Science,
// Massachusetts Institute of Technology, 1993.
//
// Shoot a test ray along +X axis. The strategy is to compare vertex Y values
// to the testing point's Y and quickly discard edges which are entirely to one
// side of the test ray. Note that CONVEX and WINDING code can be added as
// for the CrossingsTest() code; it is left out here for clarity.
//
// Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
// _point_, returns 1 if inside, 0 if outside.
template
void point_in_path_impl(PointArray &points, PathIterator &path, ResultArray &inside_flag)
{
uint8_t yflag1;
double vtx0, vty0, vtx1, vty1;
double tx, ty;
double sx, sy;
double x, y;
size_t i;
bool all_done;
size_t n = safe_first_shape(points);
std::vector yflag0(n);
std::vector subpath_flag(n);
path.rewind(0);
for (i = 0; i < n; ++i) {
inside_flag[i] = 0;
}
unsigned code = 0;
do {
if (code != agg::path_cmd_move_to) {
code = path.vertex(&x, &y);
if (code == agg::path_cmd_stop ||
(code & agg::path_cmd_end_poly) == agg::path_cmd_end_poly) {
continue;
}
}
sx = vtx0 = vtx1 = x;
sy = vty0 = vty1 = y;
for (i = 0; i < n; ++i) {
ty = points(i, 1);
if (std::isfinite(ty)) {
// get test bit for above/below X axis
yflag0[i] = (vty0 >= ty);
subpath_flag[i] = 0;
}
}
do {
code = path.vertex(&x, &y);
// The following cases denote the beginning on a new subpath
if (code == agg::path_cmd_stop ||
(code & agg::path_cmd_end_poly) == agg::path_cmd_end_poly) {
x = sx;
y = sy;
} else if (code == agg::path_cmd_move_to) {
break;
}
for (i = 0; i < n; ++i) {
tx = points(i, 0);
ty = points(i, 1);
if (!(std::isfinite(tx) && std::isfinite(ty))) {
continue;
}
yflag1 = (vty1 >= ty);
// Check if endpoints straddle (are on opposite sides) of
// X axis (i.e. the Y's differ); if so, +X ray could
// intersect this edge. The old test also checked whether
// the endpoints are both to the right or to the left of
// the test point. However, given the faster intersection
// point computation used below, this test was found to be
// a break-even proposition for most polygons and a loser
// for triangles (where 50% or more of the edges which
// survive this test will cross quadrants and so have to
// have the X intersection computed anyway). I credit
// Joseph Samosky with inspiring me to try dropping the
// "both left or both right" part of my code.
if (yflag0[i] != yflag1) {
// Check intersection of pgon segment with +X ray.
// Note if >= point's X; if so, the ray hits it. The
// division operation is avoided for the ">=" test by
// checking the sign of the first vertex wrto the test
// point; idea inspired by Joseph Samosky's and Mark
// Haigh-Hutchinson's different polygon inclusion
// tests.
if (((vty1 - ty) * (vtx0 - vtx1) >= (vtx1 - tx) * (vty0 - vty1)) == yflag1) {
subpath_flag[i] ^= 1;
}
}
// Move to the next pair of vertices, retaining info as
// possible.
yflag0[i] = yflag1;
}
vtx0 = vtx1;
vty0 = vty1;
vtx1 = x;
vty1 = y;
} while (code != agg::path_cmd_stop &&
(code & agg::path_cmd_end_poly) != agg::path_cmd_end_poly);
all_done = true;
for (i = 0; i < n; ++i) {
tx = points(i, 0);
ty = points(i, 1);
if (!(std::isfinite(tx) && std::isfinite(ty))) {
continue;
}
yflag1 = (vty1 >= ty);
if (yflag0[i] != yflag1) {
if (((vty1 - ty) * (vtx0 - vtx1) >= (vtx1 - tx) * (vty0 - vty1)) == yflag1) {
subpath_flag[i] = subpath_flag[i] ^ true;
}
}
inside_flag[i] |= subpath_flag[i];
if (inside_flag[i] == 0) {
all_done = false;
}
}
if (all_done) {
break;
}
} while (code != agg::path_cmd_stop);
}
template
inline void points_in_path(PointArray &points,
const double r,
PathIterator &path,
agg::trans_affine &trans,
ResultArray &result)
{
for (auto i = 0; i < safe_first_shape(points); ++i) {
result[i] = false;
}
if (path.total_vertices() < 3) {
return;
}
auto trans_path = agg::conv_transform{path, trans};
auto no_nans_path = PathNanRemover{trans_path, true, path.has_codes()};
auto curved_path = agg::conv_curve{no_nans_path};
if (r != 0.0) {
auto contoured_path = agg::conv_contour{curved_path};
contoured_path.width(r);
point_in_path_impl(points, contoured_path, result);
} else {
point_in_path_impl(points, curved_path, result);
}
}
template
inline bool point_in_path(
double x, double y, const double r, PathIterator &path, agg::trans_affine &trans)
{
py::ssize_t shape[] = {1, 2};
py::array_t points_arr(shape);
*points_arr.mutable_data(0, 0) = x;
*points_arr.mutable_data(0, 1) = y;
auto points = points_arr.mutable_unchecked();
int result[1];
result[0] = 0;
points_in_path(points, r, path, trans, result);
return result[0] != 0;
}
template
inline bool point_on_path(
double x, double y, const double r, PathIterator &path, agg::trans_affine &trans)
{
py::ssize_t shape[] = {1, 2};
py::array_t points_arr(shape);
*points_arr.mutable_data(0, 0) = x;
*points_arr.mutable_data(0, 1) = y;
auto points = points_arr.mutable_unchecked();
int result[1];
result[0] = 0;
auto trans_path = agg::conv_transform{path, trans};
auto nan_removed_path = PathNanRemover{trans_path, true, path.has_codes()};
auto curved_path = agg::conv_curve{nan_removed_path};
auto stroked_path = agg::conv_stroke{curved_path};
stroked_path.width(r * 2.0);
point_in_path_impl(points, stroked_path, result);
return result[0] != 0;
}
struct extent_limits
{
XY start;
XY end;
/* minpos is the minimum positive values in the data; used by log scaling. */
XY minpos;
extent_limits() : start{0,0}, end{0,0}, minpos{0,0} {
reset();
}
void reset()
{
start.x = std::numeric_limits::infinity();
start.y = std::numeric_limits::infinity();
end.x = -std::numeric_limits::infinity();
end.y = -std::numeric_limits::infinity();
minpos.x = std::numeric_limits::infinity();
minpos.y = std::numeric_limits::infinity();
}
void update(double x, double y)
{
start.x = std::min(start.x, x);
start.y = std::min(start.y, y);
end.x = std::max(end.x, x);
end.y = std::max(end.y, y);
if (x > 0.0) {
minpos.x = std::min(minpos.x, x);
}
if (y > 0.0) {
minpos.y = std::min(minpos.y, y);
}
}
};
template
void update_path_extents(PathIterator &path, agg::trans_affine &trans, extent_limits &extents)
{
double x, y;
unsigned code;
auto tpath = agg::conv_transform{path, trans};
auto nan_removed = PathNanRemover{tpath, true, path.has_codes()};
nan_removed.rewind(0);
while ((code = nan_removed.vertex(&x, &y)) != agg::path_cmd_stop) {
if ((code & agg::path_cmd_end_poly) == agg::path_cmd_end_poly) {
continue;
}
extents.update(x, y);
}
}
template
void get_path_collection_extents(agg::trans_affine &master_transform,
PathGenerator &paths,
TransformArray &transforms,
OffsetArray &offsets,
agg::trans_affine &offset_trans,
extent_limits &extent)
{
if (offsets.size() != 0 && offsets.shape(1) != 2) {
throw std::runtime_error("Offsets array must have shape (N, 2)");
}
auto Npaths = paths.size();
auto Noffsets = safe_first_shape(offsets);
auto N = std::max(Npaths, Noffsets);
auto Ntransforms = std::min(safe_first_shape(transforms), N);
agg::trans_affine trans;
extent.reset();
for (auto i = 0; i < N; ++i) {
typename PathGenerator::path_iterator path(paths(i % Npaths));
if (Ntransforms) {
py::ssize_t ti = i % Ntransforms;
trans = agg::trans_affine(transforms(ti, 0, 0),
transforms(ti, 1, 0),
transforms(ti, 0, 1),
transforms(ti, 1, 1),
transforms(ti, 0, 2),
transforms(ti, 1, 2));
} else {
trans = master_transform;
}
if (Noffsets) {
double xo = offsets(i % Noffsets, 0);
double yo = offsets(i % Noffsets, 1);
offset_trans.transform(&xo, &yo);
trans *= agg::trans_affine_translation(xo, yo);
}
update_path_extents(path, trans, extent);
}
}
template
void point_in_path_collection(double x,
double y,
double radius,
agg::trans_affine &master_transform,
PathGenerator &paths,
TransformArray &transforms,
OffsetArray &offsets,
agg::trans_affine &offset_trans,
bool filled,
std::vector &result)
{
auto Npaths = paths.size();
if (Npaths == 0) {
return;
}
auto Noffsets = safe_first_shape(offsets);
auto N = std::max(Npaths, Noffsets);
auto Ntransforms = std::min(safe_first_shape(transforms), N);
agg::trans_affine trans;
for (auto i = 0; i < N; ++i) {
typename PathGenerator::path_iterator path = paths(i % Npaths);
if (Ntransforms) {
auto ti = i % Ntransforms;
trans = agg::trans_affine(transforms(ti, 0, 0),
transforms(ti, 1, 0),
transforms(ti, 0, 1),
transforms(ti, 1, 1),
transforms(ti, 0, 2),
transforms(ti, 1, 2));
trans *= master_transform;
} else {
trans = master_transform;
}
if (Noffsets) {
double xo = offsets(i % Noffsets, 0);
double yo = offsets(i % Noffsets, 1);
offset_trans.transform(&xo, &yo);
trans *= agg::trans_affine_translation(xo, yo);
}
if (filled) {
if (point_in_path(x, y, radius, path, trans)) {
result.push_back(i);
}
} else {
if (point_on_path(x, y, radius, path, trans)) {
result.push_back(i);
}
}
}
}
template
bool path_in_path(PathIterator1 &a,
agg::trans_affine &atrans,
PathIterator2 &b,
agg::trans_affine &btrans)
{
if (a.total_vertices() < 3) {
return false;
}
auto b_path_trans = agg::conv_transform{b, btrans};
auto b_no_nans = PathNanRemover{b_path_trans, true, b.has_codes()};
auto b_curved = agg::conv_curve{b_no_nans};
double x, y;
b_curved.rewind(0);
while (b_curved.vertex(&x, &y) != agg::path_cmd_stop) {
if (!point_in_path(x, y, 0.0, a, atrans)) {
return false;
}
}
return true;
}
/** The clip_path_to_rect code here is a clean-room implementation of
the Sutherland-Hodgman clipping algorithm described here:
https://en.wikipedia.org/wiki/Sutherland-Hodgman_clipping_algorithm
*/
namespace clip_to_rect_filters
{
/* There are four different passes needed to create/remove
vertices (one for each side of the rectangle). The differences
between those passes are encapsulated in these functor classes.
*/
struct bisectx
{
double m_x;
bisectx(double x) : m_x(x)
{
}
inline XY bisect(const XY s, const XY p) const
{
double dx = p.x - s.x;
double dy = p.y - s.y;
return {
m_x,
s.y + dy * ((m_x - s.x) / dx),
};
}
};
struct xlt : public bisectx
{
xlt(double x) : bisectx(x)
{
}
inline bool is_inside(const XY point) const
{
return point.x = m_x;
}
};
struct bisecty
{
double m_y;
bisecty(double y) : m_y(y)
{
}
inline XY bisect(const XY s, const XY p) const
{
double dx = p.x - s.x;
double dy = p.y - s.y;
return {
s.x + dx * ((m_y - s.y) / dy),
m_y,
};
}
};
struct ylt : public bisecty
{
ylt(double y) : bisecty(y)
{
}
inline bool is_inside(const XY point) const
{
return point.y = m_y;
}
};
}
template
inline void clip_to_rect_one_step(const Polygon &polygon, Polygon &result, const Filter &filter)
{
bool sinside, pinside;
result.clear();
if (polygon.size() == 0) {
return;
}
auto s = polygon.back();
for (auto p : polygon) {
sinside = filter.is_inside(s);
pinside = filter.is_inside(p);
if (sinside ^ pinside) {
result.emplace_back(filter.bisect(s, p));
}
if (pinside) {
result.emplace_back(p);
}
s = p;
}
}
template
auto
clip_path_to_rect(PathIterator &path, agg::rect_d &rect, bool inside)
{
rect.normalize();
auto xmin = rect.x1, xmax = rect.x2;
auto ymin = rect.y1, ymax = rect.y2;
if (!inside) {
std::swap(xmin, xmax);
std::swap(ymin, ymax);
}
auto curve = agg::conv_curve{path};
Polygon polygon1, polygon2;
XY point;
unsigned code = 0;
curve.rewind(0);
std::vector results;
do {
// Grab the next subpath and store it in polygon1
polygon1.clear();
do {
if (code == agg::path_cmd_move_to) {
polygon1.emplace_back(point);
}
code = curve.vertex(&point.x, &point.y);
if (code == agg::path_cmd_stop) {
break;
}
if (code != agg::path_cmd_move_to) {
polygon1.emplace_back(point);
}
} while ((code & agg::path_cmd_end_poly) != agg::path_cmd_end_poly);
// The result of each step is fed into the next (note the
// swapping of polygon1 and polygon2 at each step).
clip_to_rect_one_step(polygon1, polygon2, clip_to_rect_filters::xlt(xmax));
clip_to_rect_one_step(polygon2, polygon1, clip_to_rect_filters::xgt(xmin));
clip_to_rect_one_step(polygon1, polygon2, clip_to_rect_filters::ylt(ymax));
clip_to_rect_one_step(polygon2, polygon1, clip_to_rect_filters::ygt(ymin));
// Empty polygons aren't very useful, so skip them
if (polygon1.size()) {
_finalize_polygon(results, true);
results.push_back(polygon1);
}
} while (code != agg::path_cmd_stop);
_finalize_polygon(results, true);
return results;
}
template
void affine_transform_2d(VerticesArray &vertices, agg::trans_affine &trans, ResultArray &result)
{
if (vertices.size() != 0 && vertices.shape(1) != 2) {
throw std::runtime_error("Invalid vertices array.");
}
size_t n = vertices.shape(0);
double x;
double y;
double t0;
double t1;
double t;
for (size_t i = 0; i < n; ++i) {
x = vertices(i, 0);
y = vertices(i, 1);
t0 = trans.sx * x;
t1 = trans.shx * y;
t = t0 + t1 + trans.tx;
result(i, 0) = t;
t0 = trans.shy * x;
t1 = trans.sy * y;
t = t0 + t1 + trans.ty;
result(i, 1) = t;
}
}
template
void affine_transform_1d(VerticesArray &vertices, agg::trans_affine &trans, ResultArray &result)
{
if (vertices.shape(0) != 2) {
throw std::runtime_error("Invalid vertices array.");
}
double x;
double y;
double t0;
double t1;
double t;
x = vertices(0);
y = vertices(1);
t0 = trans.sx * x;
t1 = trans.shx * y;
t = t0 + t1 + trans.tx;
result(0) = t;
t0 = trans.shy * x;
t1 = trans.sy * y;
t = t0 + t1 + trans.ty;
result(1) = t;
}
template
int count_bboxes_overlapping_bbox(agg::rect_d &a, BBoxArray &bboxes)
{
agg::rect_d b;
int count = 0;
if (a.x2 < a.x1) {
std::swap(a.x1, a.x2);
}
if (a.y2 < a.y1) {
std::swap(a.y1, a.y2);
}
size_t num_bboxes = safe_first_shape(bboxes);
for (size_t i = 0; i < num_bboxes; ++i) {
b = agg::rect_d(bboxes(i, 0, 0), bboxes(i, 0, 1), bboxes(i, 1, 0), bboxes(i, 1, 1));
if (b.x2 < b.x1) {
std::swap(b.x1, b.x2);
}
if (b.y2 < b.y1) {
std::swap(b.y1, b.y2);
}
if (!((b.x2 = a.y2))) {
++count;
}
}
return count;
}
inline bool isclose(double a, double b)
{
// relative and absolute tolerance values are chosen empirically
// it looks the atol value matters here because of round-off errors
const double rtol = 1e-10;
const double atol = 1e-13;
// as per python's math.isclose
return fabs(a-b)