1 /* -*- mode: c++; c-basic-offset: 4 -*- */
2
3 #ifndef MPL_PATH_H
4 #define MPL_PATH_H
5
6 #include <limits>
7 #include <math.h>
8 #include <vector>
9 #include <cmath>
10 #include <algorithm>
11 #include <string>
12
13 #include "agg_conv_contour.h"
14 #include "agg_conv_curve.h"
15 #include "agg_conv_stroke.h"
16 #include "agg_conv_transform.h"
17 #include "agg_path_storage.h"
18 #include "agg_trans_affine.h"
19
20 #include "path_converters.h"
21 #include "_backend_agg_basic_types.h"
22 #include "numpy_cpp.h"
23
24 /* Compatibility for PyPy3.7 before 7.3.4. */
25 #ifndef Py_DTSF_ADD_DOT_0
26 #define Py_DTSF_ADD_DOT_0 0x2
27 #endif
28
29 struct XY
30 {
31 double x;
32 double y;
33
XYXY34 XY(double x_, double y_) : x(x_), y(y_)
35 {
36 }
37
38 bool operator==(const XY& o)
39 {
40 return (x == o.x && y == o.y);
41 }
42
43 bool operator!=(const XY& o)
44 {
45 return (x != o.x || y != o.y);
46 }
47 };
48
49 typedef std::vector<XY> Polygon;
50
_finalize_polygon(std::vector<Polygon> & result,int closed_only)51 void _finalize_polygon(std::vector<Polygon> &result, int closed_only)
52 {
53 if (result.size() == 0) {
54 return;
55 }
56
57 Polygon &polygon = result.back();
58
59 /* Clean up the last polygon in the result. */
60 if (polygon.size() == 0) {
61 result.pop_back();
62 } else if (closed_only) {
63 if (polygon.size() < 3) {
64 result.pop_back();
65 } else if (polygon.front() != polygon.back()) {
66 polygon.push_back(polygon.front());
67 }
68 }
69 }
70
71 //
72 // The following function was found in the Agg 2.3 examples (interactive_polygon.cpp).
73 // It has been generalized to work on (possibly curved) polylines, rather than
74 // just polygons. The original comments have been kept intact.
75 // -- Michael Droettboom 2007-10-02
76 //
77 //======= Crossings Multiply algorithm of InsideTest ========================
78 //
79 // By Eric Haines, 3D/Eye Inc, erich@eye.com
80 //
81 // This version is usually somewhat faster than the original published in
82 // Graphics Gems IV; by turning the division for testing the X axis crossing
83 // into a tricky multiplication test this part of the test became faster,
84 // which had the additional effect of making the test for "both to left or
85 // both to right" a bit slower for triangles than simply computing the
86 // intersection each time. The main increase is in triangle testing speed,
87 // which was about 15% faster; all other polygon complexities were pretty much
88 // the same as before. On machines where division is very expensive (not the
89 // case on the HP 9000 series on which I tested) this test should be much
90 // faster overall than the old code. Your mileage may (in fact, will) vary,
91 // depending on the machine and the test data, but in general I believe this
92 // code is both shorter and faster. This test was inspired by unpublished
93 // Graphics Gems submitted by Joseph Samosky and Mark Haigh-Hutchinson.
94 // Related work by Samosky is in:
95 //
96 // Samosky, Joseph, "SectionView: A system for interactively specifying and
97 // visualizing sections through three-dimensional medical image data",
98 // M.S. Thesis, Department of Electrical Engineering and Computer Science,
99 // Massachusetts Institute of Technology, 1993.
100 //
101 // Shoot a test ray along +X axis. The strategy is to compare vertex Y values
102 // to the testing point's Y and quickly discard edges which are entirely to one
103 // side of the test ray. Note that CONVEX and WINDING code can be added as
104 // for the CrossingsTest() code; it is left out here for clarity.
105 //
106 // Input 2D polygon _pgon_ with _numverts_ number of vertices and test point
107 // _point_, returns 1 if inside, 0 if outside.
108 template <class PathIterator, class PointArray, class ResultArray>
point_in_path_impl(PointArray & points,PathIterator & path,ResultArray & inside_flag)109 void point_in_path_impl(PointArray &points, PathIterator &path, ResultArray &inside_flag)
110 {
111 uint8_t yflag1;
112 double vtx0, vty0, vtx1, vty1;
113 double tx, ty;
114 double sx, sy;
115 double x, y;
116 size_t i;
117 bool all_done;
118
119 size_t n = points.size();
120
121 std::vector<uint8_t> yflag0(n);
122 std::vector<uint8_t> subpath_flag(n);
123
124 path.rewind(0);
125
126 for (i = 0; i < n; ++i) {
127 inside_flag[i] = 0;
128 }
129
130 unsigned code = 0;
131 do {
132 if (code != agg::path_cmd_move_to) {
133 code = path.vertex(&x, &y);
134 if (code == agg::path_cmd_stop ||
135 (code & agg::path_cmd_end_poly) == agg::path_cmd_end_poly) {
136 continue;
137 }
138 }
139
140 sx = vtx0 = vtx1 = x;
141 sy = vty0 = vty1 = y;
142
143 for (i = 0; i < n; ++i) {
144 ty = points(i, 1);
145
146 if (std::isfinite(ty)) {
147 // get test bit for above/below X axis
148 yflag0[i] = (vty0 >= ty);
149
150 subpath_flag[i] = 0;
151 }
152 }
153
154 do {
155 code = path.vertex(&x, &y);
156
157 // The following cases denote the beginning on a new subpath
158 if (code == agg::path_cmd_stop ||
159 (code & agg::path_cmd_end_poly) == agg::path_cmd_end_poly) {
160 x = sx;
161 y = sy;
162 } else if (code == agg::path_cmd_move_to) {
163 break;
164 }
165
166 for (i = 0; i < n; ++i) {
167 tx = points(i, 0);
168 ty = points(i, 1);
169
170 if (!(std::isfinite(tx) && std::isfinite(ty))) {
171 continue;
172 }
173
174 yflag1 = (vty1 >= ty);
175 // Check if endpoints straddle (are on opposite sides) of
176 // X axis (i.e. the Y's differ); if so, +X ray could
177 // intersect this edge. The old test also checked whether
178 // the endpoints are both to the right or to the left of
179 // the test point. However, given the faster intersection
180 // point computation used below, this test was found to be
181 // a break-even proposition for most polygons and a loser
182 // for triangles (where 50% or more of the edges which
183 // survive this test will cross quadrants and so have to
184 // have the X intersection computed anyway). I credit
185 // Joseph Samosky with inspiring me to try dropping the
186 // "both left or both right" part of my code.
187 if (yflag0[i] != yflag1) {
188 // Check intersection of pgon segment with +X ray.
189 // Note if >= point's X; if so, the ray hits it. The
190 // division operation is avoided for the ">=" test by
191 // checking the sign of the first vertex wrto the test
192 // point; idea inspired by Joseph Samosky's and Mark
193 // Haigh-Hutchinson's different polygon inclusion
194 // tests.
195 if (((vty1 - ty) * (vtx0 - vtx1) >= (vtx1 - tx) * (vty0 - vty1)) == yflag1) {
196 subpath_flag[i] ^= 1;
197 }
198 }
199
200 // Move to the next pair of vertices, retaining info as
201 // possible.
202 yflag0[i] = yflag1;
203 }
204
205 vtx0 = vtx1;
206 vty0 = vty1;
207
208 vtx1 = x;
209 vty1 = y;
210 } while (code != agg::path_cmd_stop &&
211 (code & agg::path_cmd_end_poly) != agg::path_cmd_end_poly);
212
213 all_done = true;
214 for (i = 0; i < n; ++i) {
215 tx = points(i, 0);
216 ty = points(i, 1);
217
218 if (!(std::isfinite(tx) && std::isfinite(ty))) {
219 continue;
220 }
221
222 yflag1 = (vty1 >= ty);
223 if (yflag0[i] != yflag1) {
224 if (((vty1 - ty) * (vtx0 - vtx1) >= (vtx1 - tx) * (vty0 - vty1)) == yflag1) {
225 subpath_flag[i] = subpath_flag[i] ^ true;
226 }
227 }
228 inside_flag[i] |= subpath_flag[i];
229 if (inside_flag[i] == 0) {
230 all_done = false;
231 }
232 }
233
234 if (all_done) {
235 break;
236 }
237 } while (code != agg::path_cmd_stop);
238 }
239
240 template <class PathIterator, class PointArray, class ResultArray>
points_in_path(PointArray & points,const double r,PathIterator & path,agg::trans_affine & trans,ResultArray & result)241 inline void points_in_path(PointArray &points,
242 const double r,
243 PathIterator &path,
244 agg::trans_affine &trans,
245 ResultArray &result)
246 {
247 typedef agg::conv_transform<PathIterator> transformed_path_t;
248 typedef PathNanRemover<transformed_path_t> no_nans_t;
249 typedef agg::conv_curve<no_nans_t> curve_t;
250 typedef agg::conv_contour<curve_t> contour_t;
251
252 size_t i;
253 for (i = 0; i < points.size(); ++i) {
254 result[i] = false;
255 }
256
257 if (path.total_vertices() < 3) {
258 return;
259 }
260
261 transformed_path_t trans_path(path, trans);
262 no_nans_t no_nans_path(trans_path, true, path.has_curves());
263 curve_t curved_path(no_nans_path);
264 if (r != 0.0) {
265 contour_t contoured_path(curved_path);
266 contoured_path.width(r);
267 point_in_path_impl(points, contoured_path, result);
268 } else {
269 point_in_path_impl(points, curved_path, result);
270 }
271 }
272
273 template <class PathIterator>
point_in_path(double x,double y,const double r,PathIterator & path,agg::trans_affine & trans)274 inline bool point_in_path(
275 double x, double y, const double r, PathIterator &path, agg::trans_affine &trans)
276 {
277 npy_intp shape[] = {1, 2};
278 numpy::array_view<double, 2> points(shape);
279 points(0, 0) = x;
280 points(0, 1) = y;
281
282 int result[1];
283 result[0] = 0;
284
285 points_in_path(points, r, path, trans, result);
286
287 return result[0] != 0;
288 }
289
290 template <class PathIterator, class PointArray, class ResultArray>
points_on_path(PointArray & points,const double r,PathIterator & path,agg::trans_affine & trans,ResultArray result)291 void points_on_path(PointArray &points,
292 const double r,
293 PathIterator &path,
294 agg::trans_affine &trans,
295 ResultArray result)
296 {
297 typedef agg::conv_transform<PathIterator> transformed_path_t;
298 typedef PathNanRemover<transformed_path_t> no_nans_t;
299 typedef agg::conv_curve<no_nans_t> curve_t;
300 typedef agg::conv_stroke<curve_t> stroke_t;
301
302 size_t i;
303 for (i = 0; i < points.size(); ++i) {
304 result[i] = false;
305 }
306
307 transformed_path_t trans_path(path, trans);
308 no_nans_t nan_removed_path(trans_path, true, path.has_curves());
309 curve_t curved_path(nan_removed_path);
310 stroke_t stroked_path(curved_path);
311 stroked_path.width(r * 2.0);
312 point_in_path_impl(points, stroked_path, result);
313 }
314
315 template <class PathIterator>
point_on_path(double x,double y,const double r,PathIterator & path,agg::trans_affine & trans)316 inline bool point_on_path(
317 double x, double y, const double r, PathIterator &path, agg::trans_affine &trans)
318 {
319 npy_intp shape[] = {1, 2};
320 numpy::array_view<double, 2> points(shape);
321 points(0, 0) = x;
322 points(0, 1) = y;
323
324 int result[1];
325 result[0] = 0;
326
327 points_on_path(points, r, path, trans, result);
328
329 return result[0] != 0;
330 }
331
332 struct extent_limits
333 {
334 double x0;
335 double y0;
336 double x1;
337 double y1;
338 double xm;
339 double ym;
340 };
341
reset_limits(extent_limits & e)342 void reset_limits(extent_limits &e)
343 {
344 e.x0 = std::numeric_limits<double>::infinity();
345 e.y0 = std::numeric_limits<double>::infinity();
346 e.x1 = -std::numeric_limits<double>::infinity();
347 e.y1 = -std::numeric_limits<double>::infinity();
348 /* xm and ym are the minimum positive values in the data, used
349 by log scaling */
350 e.xm = std::numeric_limits<double>::infinity();
351 e.ym = std::numeric_limits<double>::infinity();
352 }
353
update_limits(double x,double y,extent_limits & e)354 inline void update_limits(double x, double y, extent_limits &e)
355 {
356 if (x < e.x0)
357 e.x0 = x;
358 if (y < e.y0)
359 e.y0 = y;
360 if (x > e.x1)
361 e.x1 = x;
362 if (y > e.y1)
363 e.y1 = y;
364 /* xm and ym are the minimum positive values in the data, used
365 by log scaling */
366 if (x > 0.0 && x < e.xm)
367 e.xm = x;
368 if (y > 0.0 && y < e.ym)
369 e.ym = y;
370 }
371
372 template <class PathIterator>
update_path_extents(PathIterator & path,agg::trans_affine & trans,extent_limits & extents)373 void update_path_extents(PathIterator &path, agg::trans_affine &trans, extent_limits &extents)
374 {
375 typedef agg::conv_transform<PathIterator> transformed_path_t;
376 typedef PathNanRemover<transformed_path_t> nan_removed_t;
377 double x, y;
378 unsigned code;
379
380 transformed_path_t tpath(path, trans);
381 nan_removed_t nan_removed(tpath, true, path.has_curves());
382
383 nan_removed.rewind(0);
384
385 while ((code = nan_removed.vertex(&x, &y)) != agg::path_cmd_stop) {
386 if ((code & agg::path_cmd_end_poly) == agg::path_cmd_end_poly) {
387 continue;
388 }
389 update_limits(x, y, extents);
390 }
391 }
392
393 template <class PathGenerator, class TransformArray, class OffsetArray>
get_path_collection_extents(agg::trans_affine & master_transform,PathGenerator & paths,TransformArray & transforms,OffsetArray & offsets,agg::trans_affine & offset_trans,extent_limits & extent)394 void get_path_collection_extents(agg::trans_affine &master_transform,
395 PathGenerator &paths,
396 TransformArray &transforms,
397 OffsetArray &offsets,
398 agg::trans_affine &offset_trans,
399 extent_limits &extent)
400 {
401 if (offsets.size() != 0 && offsets.dim(1) != 2) {
402 throw std::runtime_error("Offsets array must be Nx2");
403 }
404
405 size_t Npaths = paths.size();
406 size_t Noffsets = offsets.size();
407 size_t N = std::max(Npaths, Noffsets);
408 size_t Ntransforms = std::min(transforms.size(), N);
409 size_t i;
410
411 agg::trans_affine trans;
412
413 reset_limits(extent);
414
415 for (i = 0; i < N; ++i) {
416 typename PathGenerator::path_iterator path(paths(i % Npaths));
417 if (Ntransforms) {
418 size_t ti = i % Ntransforms;
419 trans = agg::trans_affine(transforms(ti, 0, 0),
420 transforms(ti, 1, 0),
421 transforms(ti, 0, 1),
422 transforms(ti, 1, 1),
423 transforms(ti, 0, 2),
424 transforms(ti, 1, 2));
425 } else {
426 trans = master_transform;
427 }
428
429 if (Noffsets) {
430 double xo = offsets(i % Noffsets, 0);
431 double yo = offsets(i % Noffsets, 1);
432 offset_trans.transform(&xo, &yo);
433 trans *= agg::trans_affine_translation(xo, yo);
434 }
435
436 update_path_extents(path, trans, extent);
437 }
438 }
439
440 template <class PathGenerator, class TransformArray, class OffsetArray>
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,e_offset_position offset_position,std::vector<int> & result)441 void point_in_path_collection(double x,
442 double y,
443 double radius,
444 agg::trans_affine &master_transform,
445 PathGenerator &paths,
446 TransformArray &transforms,
447 OffsetArray &offsets,
448 agg::trans_affine &offset_trans,
449 bool filled,
450 e_offset_position offset_position,
451 std::vector<int> &result)
452 {
453 size_t Npaths = paths.size();
454
455 if (Npaths == 0) {
456 return;
457 }
458
459 size_t Noffsets = offsets.size();
460 size_t N = std::max(Npaths, Noffsets);
461 size_t Ntransforms = std::min(transforms.size(), N);
462 size_t i;
463
464 agg::trans_affine trans;
465
466 for (i = 0; i < N; ++i) {
467 typename PathGenerator::path_iterator path = paths(i % Npaths);
468
469 if (Ntransforms) {
470 size_t ti = i % Ntransforms;
471 trans = agg::trans_affine(transforms(ti, 0, 0),
472 transforms(ti, 1, 0),
473 transforms(ti, 0, 1),
474 transforms(ti, 1, 1),
475 transforms(ti, 0, 2),
476 transforms(ti, 1, 2));
477 trans *= master_transform;
478 } else {
479 trans = master_transform;
480 }
481
482 if (Noffsets) {
483 double xo = offsets(i % Noffsets, 0);
484 double yo = offsets(i % Noffsets, 1);
485 offset_trans.transform(&xo, &yo);
486 if (offset_position == OFFSET_POSITION_DATA) {
487 trans = agg::trans_affine_translation(xo, yo) * trans;
488 } else {
489 trans *= agg::trans_affine_translation(xo, yo);
490 }
491 }
492
493 if (filled) {
494 if (point_in_path(x, y, radius, path, trans)) {
495 result.push_back(i);
496 }
497 } else {
498 if (point_on_path(x, y, radius, path, trans)) {
499 result.push_back(i);
500 }
501 }
502 }
503 }
504
505 template <class PathIterator1, class PathIterator2>
path_in_path(PathIterator1 & a,agg::trans_affine & atrans,PathIterator2 & b,agg::trans_affine & btrans)506 bool path_in_path(PathIterator1 &a,
507 agg::trans_affine &atrans,
508 PathIterator2 &b,
509 agg::trans_affine &btrans)
510 {
511 typedef agg::conv_transform<PathIterator2> transformed_path_t;
512 typedef PathNanRemover<transformed_path_t> no_nans_t;
513 typedef agg::conv_curve<no_nans_t> curve_t;
514
515 if (a.total_vertices() < 3) {
516 return false;
517 }
518
519 transformed_path_t b_path_trans(b, btrans);
520 no_nans_t b_no_nans(b_path_trans, true, b.has_curves());
521 curve_t b_curved(b_no_nans);
522
523 double x, y;
524 b_curved.rewind(0);
525 while (b_curved.vertex(&x, &y) != agg::path_cmd_stop) {
526 if (!point_in_path(x, y, 0.0, a, atrans)) {
527 return false;
528 }
529 }
530
531 return true;
532 }
533
534 /** The clip_path_to_rect code here is a clean-room implementation of
535 the Sutherland-Hodgman clipping algorithm described here:
536
537 http://en.wikipedia.org/wiki/Sutherland-Hodgman_clipping_algorithm
538 */
539
540 namespace clip_to_rect_filters
541 {
542 /* There are four different passes needed to create/remove
543 vertices (one for each side of the rectangle). The differences
544 between those passes are encapsulated in these functor classes.
545 */
546 struct bisectx
547 {
548 double m_x;
549
bisectxbisectx550 bisectx(double x) : m_x(x)
551 {
552 }
553
bisectbisectx554 inline void bisect(double sx, double sy, double px, double py, double *bx, double *by) const
555 {
556 *bx = m_x;
557 double dx = px - sx;
558 double dy = py - sy;
559 *by = sy + dy * ((m_x - sx) / dx);
560 }
561 };
562
563 struct xlt : public bisectx
564 {
xltxlt565 xlt(double x) : bisectx(x)
566 {
567 }
568
is_insidexlt569 inline bool is_inside(double x, double y) const
570 {
571 return x <= m_x;
572 }
573 };
574
575 struct xgt : public bisectx
576 {
xgtxgt577 xgt(double x) : bisectx(x)
578 {
579 }
580
is_insidexgt581 inline bool is_inside(double x, double y) const
582 {
583 return x >= m_x;
584 }
585 };
586
587 struct bisecty
588 {
589 double m_y;
590
bisectybisecty591 bisecty(double y) : m_y(y)
592 {
593 }
594
bisectbisecty595 inline void bisect(double sx, double sy, double px, double py, double *bx, double *by) const
596 {
597 *by = m_y;
598 double dx = px - sx;
599 double dy = py - sy;
600 *bx = sx + dx * ((m_y - sy) / dy);
601 }
602 };
603
604 struct ylt : public bisecty
605 {
yltylt606 ylt(double y) : bisecty(y)
607 {
608 }
609
is_insideylt610 inline bool is_inside(double x, double y) const
611 {
612 return y <= m_y;
613 }
614 };
615
616 struct ygt : public bisecty
617 {
ygtygt618 ygt(double y) : bisecty(y)
619 {
620 }
621
is_insideygt622 inline bool is_inside(double x, double y) const
623 {
624 return y >= m_y;
625 }
626 };
627 }
628
629 template <class Filter>
clip_to_rect_one_step(const Polygon & polygon,Polygon & result,const Filter & filter)630 inline void clip_to_rect_one_step(const Polygon &polygon, Polygon &result, const Filter &filter)
631 {
632 double sx, sy, px, py, bx, by;
633 bool sinside, pinside;
634 result.clear();
635
636 if (polygon.size() == 0) {
637 return;
638 }
639
640 sx = polygon.back().x;
641 sy = polygon.back().y;
642 for (Polygon::const_iterator i = polygon.begin(); i != polygon.end(); ++i) {
643 px = i->x;
644 py = i->y;
645
646 sinside = filter.is_inside(sx, sy);
647 pinside = filter.is_inside(px, py);
648
649 if (sinside ^ pinside) {
650 filter.bisect(sx, sy, px, py, &bx, &by);
651 result.push_back(XY(bx, by));
652 }
653
654 if (pinside) {
655 result.push_back(XY(px, py));
656 }
657
658 sx = px;
659 sy = py;
660 }
661 }
662
663 template <class PathIterator>
664 void
clip_path_to_rect(PathIterator & path,agg::rect_d & rect,bool inside,std::vector<Polygon> & results)665 clip_path_to_rect(PathIterator &path, agg::rect_d &rect, bool inside, std::vector<Polygon> &results)
666 {
667 double xmin, ymin, xmax, ymax;
668 if (rect.x1 < rect.x2) {
669 xmin = rect.x1;
670 xmax = rect.x2;
671 } else {
672 xmin = rect.x2;
673 xmax = rect.x1;
674 }
675
676 if (rect.y1 < rect.y2) {
677 ymin = rect.y1;
678 ymax = rect.y2;
679 } else {
680 ymin = rect.y2;
681 ymax = rect.y1;
682 }
683
684 if (!inside) {
685 std::swap(xmin, xmax);
686 std::swap(ymin, ymax);
687 }
688
689 typedef agg::conv_curve<PathIterator> curve_t;
690 curve_t curve(path);
691
692 Polygon polygon1, polygon2;
693 double x = 0, y = 0;
694 unsigned code = 0;
695 curve.rewind(0);
696
697 do {
698 // Grab the next subpath and store it in polygon1
699 polygon1.clear();
700 do {
701 if (code == agg::path_cmd_move_to) {
702 polygon1.push_back(XY(x, y));
703 }
704
705 code = curve.vertex(&x, &y);
706
707 if (code == agg::path_cmd_stop) {
708 break;
709 }
710
711 if (code != agg::path_cmd_move_to) {
712 polygon1.push_back(XY(x, y));
713 }
714 } while ((code & agg::path_cmd_end_poly) != agg::path_cmd_end_poly);
715
716 // The result of each step is fed into the next (note the
717 // swapping of polygon1 and polygon2 at each step).
718 clip_to_rect_one_step(polygon1, polygon2, clip_to_rect_filters::xlt(xmax));
719 clip_to_rect_one_step(polygon2, polygon1, clip_to_rect_filters::xgt(xmin));
720 clip_to_rect_one_step(polygon1, polygon2, clip_to_rect_filters::ylt(ymax));
721 clip_to_rect_one_step(polygon2, polygon1, clip_to_rect_filters::ygt(ymin));
722
723 // Empty polygons aren't very useful, so skip them
724 if (polygon1.size()) {
725 _finalize_polygon(results, 1);
726 results.push_back(polygon1);
727 }
728 } while (code != agg::path_cmd_stop);
729
730 _finalize_polygon(results, 1);
731 }
732
733 template <class VerticesArray, class ResultArray>
affine_transform_2d(VerticesArray & vertices,agg::trans_affine & trans,ResultArray & result)734 void affine_transform_2d(VerticesArray &vertices, agg::trans_affine &trans, ResultArray &result)
735 {
736 if (vertices.size() != 0 && vertices.dim(1) != 2) {
737 throw std::runtime_error("Invalid vertices array.");
738 }
739
740 size_t n = vertices.size();
741 double x;
742 double y;
743 double t0;
744 double t1;
745 double t;
746
747 for (size_t i = 0; i < n; ++i) {
748 x = vertices(i, 0);
749 y = vertices(i, 1);
750
751 t0 = trans.sx * x;
752 t1 = trans.shx * y;
753 t = t0 + t1 + trans.tx;
754 result(i, 0) = t;
755
756 t0 = trans.shy * x;
757 t1 = trans.sy * y;
758 t = t0 + t1 + trans.ty;
759 result(i, 1) = t;
760 }
761 }
762
763 template <class VerticesArray, class ResultArray>
affine_transform_1d(VerticesArray & vertices,agg::trans_affine & trans,ResultArray & result)764 void affine_transform_1d(VerticesArray &vertices, agg::trans_affine &trans, ResultArray &result)
765 {
766 if (vertices.dim(0) != 2) {
767 throw std::runtime_error("Invalid vertices array.");
768 }
769
770 double x;
771 double y;
772 double t0;
773 double t1;
774 double t;
775
776 x = vertices(0);
777 y = vertices(1);
778
779 t0 = trans.sx * x;
780 t1 = trans.shx * y;
781 t = t0 + t1 + trans.tx;
782 result(0) = t;
783
784 t0 = trans.shy * x;
785 t1 = trans.sy * y;
786 t = t0 + t1 + trans.ty;
787 result(1) = t;
788 }
789
790 template <class BBoxArray>
count_bboxes_overlapping_bbox(agg::rect_d & a,BBoxArray & bboxes)791 int count_bboxes_overlapping_bbox(agg::rect_d &a, BBoxArray &bboxes)
792 {
793 agg::rect_d b;
794 int count = 0;
795
796 if (a.x2 < a.x1) {
797 std::swap(a.x1, a.x2);
798 }
799 if (a.y2 < a.y1) {
800 std::swap(a.y1, a.y2);
801 }
802
803 size_t num_bboxes = bboxes.size();
804 for (size_t i = 0; i < num_bboxes; ++i) {
805 b = agg::rect_d(bboxes(i, 0, 0), bboxes(i, 0, 1), bboxes(i, 1, 0), bboxes(i, 1, 1));
806
807 if (b.x2 < b.x1) {
808 std::swap(b.x1, b.x2);
809 }
810 if (b.y2 < b.y1) {
811 std::swap(b.y1, b.y2);
812 }
813 if (!((b.x2 <= a.x1) || (b.y2 <= a.y1) || (b.x1 >= a.x2) || (b.y1 >= a.y2))) {
814 ++count;
815 }
816 }
817
818 return count;
819 }
820
821
isclose(double a,double b)822 inline bool isclose(double a, double b)
823 {
824 // relative and absolute tolerance values are chosen empirically
825 // it looks the atol value matters here because of round-off errors
826 const double rtol = 1e-10;
827 const double atol = 1e-13;
828
829 // as per python's math.isclose
830 return fabs(a-b) <= fmax(rtol * fmax(fabs(a), fabs(b)), atol);
831 }
832
833
segments_intersect(const double & x1,const double & y1,const double & x2,const double & y2,const double & x3,const double & y3,const double & x4,const double & y4)834 inline bool segments_intersect(const double &x1,
835 const double &y1,
836 const double &x2,
837 const double &y2,
838 const double &x3,
839 const double &y3,
840 const double &x4,
841 const double &y4)
842 {
843 // determinant
844 double den = ((y4 - y3) * (x2 - x1)) - ((x4 - x3) * (y2 - y1));
845
846 // If den == 0 we have two possibilities:
847 if (isclose(den, 0.0)) {
848 float t_area = (x2*y3 - x3*y2) - x1*(y3 - y2) + y1*(x3 - x2);
849 // 1 - If the area of the triangle made by the 3 first points (2 from the first segment
850 // plus one from the second) is zero, they are collinear
851 if (isclose(t_area, 0.0)) {
852 if (x1 == x2 && x2 == x3) { // segments have infinite slope (vertical lines)
853 // and lie on the same line
854 return (fmin(y1, y2) <= fmin(y3, y4) && fmin(y3, y4) <= fmax(y1, y2)) ||
855 (fmin(y3, y4) <= fmin(y1, y2) && fmin(y1, y2) <= fmax(y3, y4));
856 }
857 else {
858 return (fmin(x1, x2) <= fmin(x3, x4) && fmin(x3, x4) <= fmax(x1, x2)) ||
859 (fmin(x3, x4) <= fmin(x1, x2) && fmin(x1, x2) <= fmax(x3, x4));
860
861 }
862 }
863 // 2 - If t_area is not zero, the segments are parallel, but not collinear
864 else {
865 return false;
866 }
867 }
868
869 const double n1 = ((x4 - x3) * (y1 - y3)) - ((y4 - y3) * (x1 - x3));
870 const double n2 = ((x2 - x1) * (y1 - y3)) - ((y2 - y1) * (x1 - x3));
871
872 const double u1 = n1 / den;
873 const double u2 = n2 / den;
874
875 return ((u1 > 0.0 || isclose(u1, 0.0)) &&
876 (u1 < 1.0 || isclose(u1, 1.0)) &&
877 (u2 > 0.0 || isclose(u2, 0.0)) &&
878 (u2 < 1.0 || isclose(u2, 1.0)));
879 }
880
881 template <class PathIterator1, class PathIterator2>
path_intersects_path(PathIterator1 & p1,PathIterator2 & p2)882 bool path_intersects_path(PathIterator1 &p1, PathIterator2 &p2)
883 {
884
885 typedef PathNanRemover<py::PathIterator> no_nans_t;
886 typedef agg::conv_curve<no_nans_t> curve_t;
887
888 if (p1.total_vertices() < 2 || p2.total_vertices() < 2) {
889 return false;
890 }
891
892 no_nans_t n1(p1, true, p1.has_curves());
893 no_nans_t n2(p2, true, p2.has_curves());
894
895 curve_t c1(n1);
896 curve_t c2(n2);
897
898 double x11, y11, x12, y12;
899 double x21, y21, x22, y22;
900
901 c1.vertex(&x11, &y11);
902 while (c1.vertex(&x12, &y12) != agg::path_cmd_stop) {
903 // if the segment in path 1 is (almost) 0 length, skip to next vertex
904 if ((isclose((x11 - x12) * (x11 - x12) + (y11 - y12) * (y11 - y12), 0))){
905 continue;
906 }
907 c2.rewind(0);
908 c2.vertex(&x21, &y21);
909
910
911 while (c2.vertex(&x22, &y22) != agg::path_cmd_stop) {
912 // if the segment in path 2 is (almost) 0 length, skip to next vertex
913 if ((isclose((x21 - x22) * (x21 - x22) + (y21 - y22) * (y21 - y22), 0))){
914 continue;
915 }
916
917 if (segments_intersect(x11, y11, x12, y12, x21, y21, x22, y22)) {
918 return true;
919 }
920 x21 = x22;
921 y21 = y22;
922 }
923 x11 = x12;
924 y11 = y12;
925 }
926
927 return false;
928 }
929
930 // returns whether the segment from (x1,y1) to (x2,y2)
931 // intersects the rectangle centered at (cx,cy) with size (w,h)
932 // see doc/segment_intersects_rectangle.svg for a more detailed explanation
segment_intersects_rectangle(double x1,double y1,double x2,double y2,double cx,double cy,double w,double h)933 inline bool segment_intersects_rectangle(double x1, double y1,
934 double x2, double y2,
935 double cx, double cy,
936 double w, double h)
937 {
938 return fabs(x1 + x2 - 2.0 * cx) < fabs(x1 - x2) + w &&
939 fabs(y1 + y2 - 2.0 * cy) < fabs(y1 - y2) + h &&
940 2.0 * fabs((x1 - cx) * (y1 - y2) - (y1 - cy) * (x1 - x2)) <
941 w * fabs(y1 - y2) + h * fabs(x1 - x2);
942 }
943
944 template <class PathIterator>
path_intersects_rectangle(PathIterator & path,double rect_x1,double rect_y1,double rect_x2,double rect_y2,bool filled)945 bool path_intersects_rectangle(PathIterator &path,
946 double rect_x1, double rect_y1,
947 double rect_x2, double rect_y2,
948 bool filled)
949 {
950 typedef PathNanRemover<py::PathIterator> no_nans_t;
951 typedef agg::conv_curve<no_nans_t> curve_t;
952
953 if (path.total_vertices() == 0) {
954 return false;
955 }
956
957 no_nans_t no_nans(path, true, path.has_curves());
958 curve_t curve(no_nans);
959
960 double cx = (rect_x1 + rect_x2) * 0.5, cy = (rect_y1 + rect_y2) * 0.5;
961 double w = fabs(rect_x1 - rect_x2), h = fabs(rect_y1 - rect_y2);
962
963 double x1, y1, x2, y2;
964
965 curve.vertex(&x1, &y1);
966 if (2.0 * fabs(x1 - cx) <= w && 2.0 * fabs(y1 - cy) <= h) {
967 return true;
968 }
969
970 while (curve.vertex(&x2, &y2) != agg::path_cmd_stop) {
971 if (segment_intersects_rectangle(x1, y1, x2, y2, cx, cy, w, h)) {
972 return true;
973 }
974 x1 = x2;
975 y1 = y2;
976 }
977
978 if (filled) {
979 agg::trans_affine trans;
980 if (point_in_path(cx, cy, 0.0, path, trans)) {
981 return true;
982 }
983 }
984
985 return false;
986 }
987
988 template <class PathIterator>
convert_path_to_polygons(PathIterator & path,agg::trans_affine & trans,double width,double height,int closed_only,std::vector<Polygon> & result)989 void convert_path_to_polygons(PathIterator &path,
990 agg::trans_affine &trans,
991 double width,
992 double height,
993 int closed_only,
994 std::vector<Polygon> &result)
995 {
996 typedef agg::conv_transform<py::PathIterator> transformed_path_t;
997 typedef PathNanRemover<transformed_path_t> nan_removal_t;
998 typedef PathClipper<nan_removal_t> clipped_t;
999 typedef PathSimplifier<clipped_t> simplify_t;
1000 typedef agg::conv_curve<simplify_t> curve_t;
1001
1002 bool do_clip = width != 0.0 && height != 0.0;
1003 bool simplify = path.should_simplify();
1004
1005 transformed_path_t tpath(path, trans);
1006 nan_removal_t nan_removed(tpath, true, path.has_curves());
1007 clipped_t clipped(nan_removed, do_clip && !path.has_curves(), width, height);
1008 simplify_t simplified(clipped, simplify, path.simplify_threshold());
1009 curve_t curve(simplified);
1010
1011 result.push_back(Polygon());
1012 Polygon *polygon = &result.back();
1013 double x, y;
1014 unsigned code;
1015
1016 while ((code = curve.vertex(&x, &y)) != agg::path_cmd_stop) {
1017 if ((code & agg::path_cmd_end_poly) == agg::path_cmd_end_poly) {
1018 _finalize_polygon(result, 1);
1019 result.push_back(Polygon());
1020 polygon = &result.back();
1021 } else {
1022 if (code == agg::path_cmd_move_to) {
1023 _finalize_polygon(result, closed_only);
1024 result.push_back(Polygon());
1025 polygon = &result.back();
1026 }
1027 polygon->push_back(XY(x, y));
1028 }
1029 }
1030
1031 _finalize_polygon(result, closed_only);
1032 }
1033
1034 template <class VertexSource>
1035 void
__cleanup_path(VertexSource & source,std::vector<double> & vertices,std::vector<npy_uint8> & codes)1036 __cleanup_path(VertexSource &source, std::vector<double> &vertices, std::vector<npy_uint8> &codes)
1037 {
1038 unsigned code;
1039 double x, y;
1040 do {
1041 code = source.vertex(&x, &y);
1042 vertices.push_back(x);
1043 vertices.push_back(y);
1044 codes.push_back((npy_uint8)code);
1045 } while (code != agg::path_cmd_stop);
1046 }
1047
1048 template <class PathIterator>
cleanup_path(PathIterator & path,agg::trans_affine & trans,bool remove_nans,bool do_clip,const agg::rect_base<double> & rect,e_snap_mode snap_mode,double stroke_width,bool do_simplify,bool return_curves,SketchParams sketch_params,std::vector<double> & vertices,std::vector<unsigned char> & codes)1049 void cleanup_path(PathIterator &path,
1050 agg::trans_affine &trans,
1051 bool remove_nans,
1052 bool do_clip,
1053 const agg::rect_base<double> &rect,
1054 e_snap_mode snap_mode,
1055 double stroke_width,
1056 bool do_simplify,
1057 bool return_curves,
1058 SketchParams sketch_params,
1059 std::vector<double> &vertices,
1060 std::vector<unsigned char> &codes)
1061 {
1062 typedef agg::conv_transform<py::PathIterator> transformed_path_t;
1063 typedef PathNanRemover<transformed_path_t> nan_removal_t;
1064 typedef PathClipper<nan_removal_t> clipped_t;
1065 typedef PathSnapper<clipped_t> snapped_t;
1066 typedef PathSimplifier<snapped_t> simplify_t;
1067 typedef agg::conv_curve<simplify_t> curve_t;
1068 typedef Sketch<curve_t> sketch_t;
1069
1070 transformed_path_t tpath(path, trans);
1071 nan_removal_t nan_removed(tpath, remove_nans, path.has_curves());
1072 clipped_t clipped(nan_removed, do_clip && !path.has_curves(), rect);
1073 snapped_t snapped(clipped, snap_mode, path.total_vertices(), stroke_width);
1074 simplify_t simplified(snapped, do_simplify, path.simplify_threshold());
1075
1076 vertices.reserve(path.total_vertices() * 2);
1077 codes.reserve(path.total_vertices());
1078
1079 if (return_curves && sketch_params.scale == 0.0) {
1080 __cleanup_path(simplified, vertices, codes);
1081 } else {
1082 curve_t curve(simplified);
1083 sketch_t sketch(curve, sketch_params.scale, sketch_params.length, sketch_params.randomness);
1084 __cleanup_path(sketch, vertices, codes);
1085 }
1086 }
1087
quad2cubic(double x0,double y0,double x1,double y1,double x2,double y2,double * outx,double * outy)1088 void quad2cubic(double x0, double y0,
1089 double x1, double y1,
1090 double x2, double y2,
1091 double *outx, double *outy)
1092 {
1093
1094 outx[0] = x0 + 2./3. * (x1 - x0);
1095 outy[0] = y0 + 2./3. * (y1 - y0);
1096 outx[1] = outx[0] + 1./3. * (x2 - x0);
1097 outy[1] = outy[0] + 1./3. * (y2 - y0);
1098 outx[2] = x2;
1099 outy[2] = y2;
1100 }
1101
1102
__add_number(double val,char format_code,int precision,std::string & buffer)1103 void __add_number(double val, char format_code, int precision,
1104 std::string& buffer)
1105 {
1106 if (precision == -1) {
1107 // Special-case for compat with old ttconv code, which *truncated*
1108 // values with a cast to int instead of rounding them as printf
1109 // would do. The only point where non-integer values arise is from
1110 // quad2cubic conversion (as we already perform a first truncation
1111 // on Python's side), which can introduce additional floating point
1112 // error (by adding 2/3 delta-x and then 1/3 delta-x), so compensate by
1113 // first rounding to the closest 1/3 and then truncating.
1114 char str[255];
1115 PyOS_snprintf(str, 255, "%d", (int)(round(val * 3)) / 3);
1116 buffer += str;
1117 } else {
1118 char *str = PyOS_double_to_string(
1119 val, format_code, precision, Py_DTSF_ADD_DOT_0, NULL);
1120 // Delete trailing zeros and decimal point
1121 char *c = str + strlen(str) - 1; // Start at last character.
1122 // Rewind through all the zeros and, if present, the trailing decimal
1123 // point. Py_DTSF_ADD_DOT_0 ensures we won't go past the start of str.
1124 while (*c == '0') {
1125 --c;
1126 }
1127 if (*c == '.') {
1128 --c;
1129 }
1130 try {
1131 buffer.append(str, c + 1);
1132 } catch (std::bad_alloc& e) {
1133 PyMem_Free(str);
1134 throw e;
1135 }
1136 PyMem_Free(str);
1137 }
1138 }
1139
1140
1141 template <class PathIterator>
__convert_to_string(PathIterator & path,int precision,char ** codes,bool postfix,std::string & buffer)1142 bool __convert_to_string(PathIterator &path,
1143 int precision,
1144 char **codes,
1145 bool postfix,
1146 std::string& buffer)
1147 {
1148 const char format_code = 'f';
1149
1150 double x[3];
1151 double y[3];
1152 double last_x = 0.0;
1153 double last_y = 0.0;
1154
1155 const int sizes[] = { 1, 1, 2, 3 };
1156 int size = 0;
1157 unsigned code;
1158
1159 while ((code = path.vertex(&x[0], &y[0])) != agg::path_cmd_stop) {
1160 if (code == 0x4f) {
1161 buffer += codes[4];
1162 } else if (code < 5) {
1163 size = sizes[code - 1];
1164
1165 for (int i = 1; i < size; ++i) {
1166 unsigned subcode = path.vertex(&x[i], &y[i]);
1167 if (subcode != code) {
1168 return false;
1169 }
1170 }
1171
1172 /* For formats that don't support quad curves, convert to
1173 cubic curves */
1174 if (code == CURVE3 && codes[code - 1][0] == '\0') {
1175 quad2cubic(last_x, last_y, x[0], y[0], x[1], y[1], x, y);
1176 code++;
1177 size = 3;
1178 }
1179
1180 if (!postfix) {
1181 buffer += codes[code - 1];
1182 buffer += ' ';
1183 }
1184
1185 for (int i = 0; i < size; ++i) {
1186 __add_number(x[i], format_code, precision, buffer);
1187 buffer += ' ';
1188 __add_number(y[i], format_code, precision, buffer);
1189 buffer += ' ';
1190 }
1191
1192 if (postfix) {
1193 buffer += codes[code - 1];
1194 }
1195
1196 last_x = x[size - 1];
1197 last_y = y[size - 1];
1198 } else {
1199 // Unknown code value
1200 return false;
1201 }
1202
1203 buffer += '\n';
1204 }
1205
1206 return true;
1207 }
1208
1209 template <class PathIterator>
convert_to_string(PathIterator & path,agg::trans_affine & trans,agg::rect_d & clip_rect,bool simplify,SketchParams sketch_params,int precision,char ** codes,bool postfix,std::string & buffer)1210 bool convert_to_string(PathIterator &path,
1211 agg::trans_affine &trans,
1212 agg::rect_d &clip_rect,
1213 bool simplify,
1214 SketchParams sketch_params,
1215 int precision,
1216 char **codes,
1217 bool postfix,
1218 std::string& buffer)
1219 {
1220 size_t buffersize;
1221 typedef agg::conv_transform<py::PathIterator> transformed_path_t;
1222 typedef PathNanRemover<transformed_path_t> nan_removal_t;
1223 typedef PathClipper<nan_removal_t> clipped_t;
1224 typedef PathSimplifier<clipped_t> simplify_t;
1225 typedef agg::conv_curve<simplify_t> curve_t;
1226 typedef Sketch<curve_t> sketch_t;
1227
1228 bool do_clip = (clip_rect.x1 < clip_rect.x2 && clip_rect.y1 < clip_rect.y2);
1229
1230 transformed_path_t tpath(path, trans);
1231 nan_removal_t nan_removed(tpath, true, path.has_curves());
1232 clipped_t clipped(nan_removed, do_clip && !path.has_curves(), clip_rect);
1233 simplify_t simplified(clipped, simplify, path.simplify_threshold());
1234
1235 buffersize = path.total_vertices() * (precision + 5) * 4;
1236 if (buffersize == 0) {
1237 return true;
1238 }
1239
1240 if (sketch_params.scale != 0.0) {
1241 buffersize *= 10;
1242 }
1243
1244 buffer.reserve(buffersize);
1245
1246 if (sketch_params.scale == 0.0) {
1247 return __convert_to_string(simplified, precision, codes, postfix, buffer);
1248 } else {
1249 curve_t curve(simplified);
1250 sketch_t sketch(curve, sketch_params.scale, sketch_params.length, sketch_params.randomness);
1251 return __convert_to_string(sketch, precision, codes, postfix, buffer);
1252 }
1253
1254 }
1255
1256 template<class T>
1257 struct _is_sorted
1258 {
operator_is_sorted1259 bool operator()(PyArrayObject *array)
1260 {
1261 npy_intp size;
1262 npy_intp i;
1263 T last_value;
1264 T current_value;
1265
1266 size = PyArray_DIM(array, 0);
1267
1268 // std::isnan is only in C++11, which we don't yet require,
1269 // so we use the "self == self" trick
1270 for (i = 0; i < size; ++i) {
1271 last_value = *((T *)PyArray_GETPTR1(array, i));
1272 if (last_value == last_value) {
1273 break;
1274 }
1275 }
1276
1277 if (i == size) {
1278 // The whole array is non-finite
1279 return false;
1280 }
1281
1282 for (; i < size; ++i) {
1283 current_value = *((T *)PyArray_GETPTR1(array, i));
1284 if (current_value == current_value) {
1285 if (current_value < last_value) {
1286 return false;
1287 }
1288 last_value = current_value;
1289 }
1290 }
1291
1292 return true;
1293 }
1294 };
1295
1296
1297 template<class T>
1298 struct _is_sorted_int
1299 {
operator_is_sorted_int1300 bool operator()(PyArrayObject *array)
1301 {
1302 npy_intp size;
1303 npy_intp i;
1304 T last_value;
1305 T current_value;
1306
1307 size = PyArray_DIM(array, 0);
1308
1309 last_value = *((T *)PyArray_GETPTR1(array, 0));
1310
1311 for (i = 1; i < size; ++i) {
1312 current_value = *((T *)PyArray_GETPTR1(array, i));
1313 if (current_value < last_value) {
1314 return false;
1315 }
1316 last_value = current_value;
1317 }
1318
1319 return true;
1320 }
1321 };
1322
1323
1324 #endif
1325