1 #define INKSCAPE_HELPER_GEOM_CPP
3 /**
4 * Specific geometry functions for Inkscape, not provided my lib2geom.
5 *
6 * Author:
7 * Johan Engelen <goejendaagh@zonnet.nl>
8 *
9 * Copyright (C) 2008 Johan Engelen
10 *
11 * Released under GNU GPL
12 */
14 #include "helper/geom.h"
15 #include <typeinfo>
16 #include <2geom/pathvector.h>
17 #include <2geom/path.h>
18 #include <2geom/bezier-curve.h>
19 #include <2geom/hvlinesegment.h>
20 #include <2geom/transforms.h>
21 #include <2geom/rect.h>
22 #include <2geom/coord.h>
23 #include <2geom/sbasis-to-bezier.h>
24 #include <libnr/nr-convert2geom.h>
25 #include <glibmm.h>
27 using Geom::X;
28 using Geom::Y;
30 #define NR_HUGE 1e18
32 //#################################################################################
33 // BOUNDING BOX CALCULATIONS
35 /* Fast bbox calculation */
36 /* Thanks to Nathan Hurst for suggesting it */
37 static void
38 cubic_bbox (Geom::Coord x000, Geom::Coord y000, Geom::Coord x001, Geom::Coord y001, Geom::Coord x011, Geom::Coord y011, Geom::Coord x111, Geom::Coord y111, Geom::Rect &bbox)
39 {
40 Geom::Coord a, b, c, D;
42 bbox[0].extendTo(x111);
43 bbox[1].extendTo(y111);
45 // It already contains (x000,y000) and (x111,y111)
46 // All points of the Bezier lie in the convex hull of (x000,y000), (x001,y001), (x011,y011) and (x111,y111)
47 // So, if it also contains (x001,y001) and (x011,y011) we don't have to compute anything else!
48 // Note that we compute it for the X and Y range separately to make it easier to use them below
49 bool containsXrange = bbox[0].contains(x001) && bbox[0].contains(x011);
50 bool containsYrange = bbox[1].contains(y001) && bbox[1].contains(y011);
52 /*
53 * xttt = s * (s * (s * x000 + t * x001) + t * (s * x001 + t * x011)) + t * (s * (s * x001 + t * x011) + t * (s * x011 + t * x111))
54 * xttt = s * (s2 * x000 + s * t * x001 + t * s * x001 + t2 * x011) + t * (s2 * x001 + s * t * x011 + t * s * x011 + t2 * x111)
55 * xttt = s * (s2 * x000 + 2 * st * x001 + t2 * x011) + t * (s2 * x001 + 2 * st * x011 + t2 * x111)
56 * xttt = s3 * x000 + 2 * s2t * x001 + st2 * x011 + s2t * x001 + 2st2 * x011 + t3 * x111
57 * xttt = s3 * x000 + 3s2t * x001 + 3st2 * x011 + t3 * x111
58 * xttt = s3 * x000 + (1 - s) 3s2 * x001 + (1 - s) * (1 - s) * 3s * x011 + (1 - s) * (1 - s) * (1 - s) * x111
59 * xttt = s3 * x000 + (3s2 - 3s3) * x001 + (3s - 6s2 + 3s3) * x011 + (1 - 2s + s2 - s + 2s2 - s3) * x111
60 * xttt = (x000 - 3 * x001 + 3 * x011 - x111) * s3 +
61 * ( 3 * x001 - 6 * x011 + 3 * x111) * s2 +
62 * ( 3 * x011 - 3 * x111) * s +
63 * ( x111)
64 * xttt' = (3 * x000 - 9 * x001 + 9 * x011 - 3 * x111) * s2 +
65 * ( 6 * x001 - 12 * x011 + 6 * x111) * s +
66 * ( 3 * x011 - 3 * x111)
67 */
69 if (!containsXrange) {
70 a = 3 * x000 - 9 * x001 + 9 * x011 - 3 * x111;
71 b = 6 * x001 - 12 * x011 + 6 * x111;
72 c = 3 * x011 - 3 * x111;
74 /*
75 * s = (-b +/- sqrt (b * b - 4 * a * c)) / 2 * a;
76 */
77 if (fabs (a) < Geom::EPSILON) {
78 /* s = -c / b */
79 if (fabs (b) > Geom::EPSILON) {
80 double s, t, xttt;
81 s = -c / b;
82 if ((s > 0.0) && (s < 1.0)) {
83 t = 1.0 - s;
84 xttt = s * s * s * x000 + 3 * s * s * t * x001 + 3 * s * t * t * x011 + t * t * t * x111;
85 bbox[0].extendTo(xttt);
86 }
87 }
88 } else {
89 /* s = (-b +/- sqrt (b * b - 4 * a * c)) / 2 * a; */
90 D = b * b - 4 * a * c;
91 if (D >= 0.0) {
92 Geom::Coord d, s, t, xttt;
93 /* Have solution */
94 d = sqrt (D);
95 s = (-b + d) / (2 * a);
96 if ((s > 0.0) && (s < 1.0)) {
97 t = 1.0 - s;
98 xttt = s * s * s * x000 + 3 * s * s * t * x001 + 3 * s * t * t * x011 + t * t * t * x111;
99 bbox[0].extendTo(xttt);
100 }
101 s = (-b - d) / (2 * a);
102 if ((s > 0.0) && (s < 1.0)) {
103 t = 1.0 - s;
104 xttt = s * s * s * x000 + 3 * s * s * t * x001 + 3 * s * t * t * x011 + t * t * t * x111;
105 bbox[0].extendTo(xttt);
106 }
107 }
108 }
109 }
111 if (!containsYrange) {
112 a = 3 * y000 - 9 * y001 + 9 * y011 - 3 * y111;
113 b = 6 * y001 - 12 * y011 + 6 * y111;
114 c = 3 * y011 - 3 * y111;
116 if (fabs (a) < Geom::EPSILON) {
117 /* s = -c / b */
118 if (fabs (b) > Geom::EPSILON) {
119 double s, t, yttt;
120 s = -c / b;
121 if ((s > 0.0) && (s < 1.0)) {
122 t = 1.0 - s;
123 yttt = s * s * s * y000 + 3 * s * s * t * y001 + 3 * s * t * t * y011 + t * t * t * y111;
124 bbox[1].extendTo(yttt);
125 }
126 }
127 } else {
128 /* s = (-b +/- sqrt (b * b - 4 * a * c)) / 2 * a; */
129 D = b * b - 4 * a * c;
130 if (D >= 0.0) {
131 Geom::Coord d, s, t, yttt;
132 /* Have solution */
133 d = sqrt (D);
134 s = (-b + d) / (2 * a);
135 if ((s > 0.0) && (s < 1.0)) {
136 t = 1.0 - s;
137 yttt = s * s * s * y000 + 3 * s * s * t * y001 + 3 * s * t * t * y011 + t * t * t * y111;
138 bbox[1].extendTo(yttt);
139 }
140 s = (-b - d) / (2 * a);
141 if ((s > 0.0) && (s < 1.0)) {
142 t = 1.0 - s;
143 yttt = s * s * s * y000 + 3 * s * s * t * y001 + 3 * s * t * t * y011 + t * t * t * y111;
144 bbox[1].extendTo(yttt);
145 }
146 }
147 }
148 }
149 }
151 Geom::Rect
152 bounds_fast_transformed(Geom::PathVector const & pv, Geom::Matrix const & t)
153 {
154 return bounds_exact_transformed(pv, t); //use this as it is faster for now! :)
155 // return Geom::bounds_fast(pv * t);
156 }
158 Geom::Rect
159 bounds_exact_transformed(Geom::PathVector const & pv, Geom::Matrix const & t)
160 {
161 Geom::Rect bbox;
163 if (pv.empty())
164 return bbox;
166 Geom::Point initial = pv.front().initialPoint() * t;
167 bbox = Geom::Rect(initial, initial); // obtain well defined bbox as starting point to unionWith
169 for (Geom::PathVector::const_iterator it = pv.begin(); it != pv.end(); ++it) {
170 bbox.expandTo(it->initialPoint() * t);
172 // don't loop including closing segment, since that segment can never increase the bbox
173 for (Geom::Path::const_iterator cit = it->begin(); cit != it->end_open(); ++cit) {
174 Geom::Curve const &c = *cit;
176 if( dynamic_cast<Geom::LineSegment const*>(&c) ||
177 dynamic_cast<Geom::HLineSegment const*>(&c) ||
178 dynamic_cast<Geom::VLineSegment const*>(&c) )
179 {
180 bbox.expandTo( c.finalPoint() * t );
181 }
182 else if(Geom::CubicBezier const *cubic_bezier = dynamic_cast<Geom::CubicBezier const *>(&c))
183 {
184 Geom::Point c0 = (*cubic_bezier)[0] * t;
185 Geom::Point c1 = (*cubic_bezier)[1] * t;
186 Geom::Point c2 = (*cubic_bezier)[2] * t;
187 Geom::Point c3 = (*cubic_bezier)[3] * t;
188 cubic_bbox( c0[0], c0[1],
189 c1[0], c1[1],
190 c2[0], c2[1],
191 c3[0], c3[1],
192 bbox );
193 }
194 else
195 {
196 // should handle all not-so-easy curves:
197 Geom::Curve *ctemp = cit->transformed(t);
198 bbox.unionWith( ctemp->boundsExact());
199 delete ctemp;
200 }
201 }
202 }
203 //return Geom::bounds_exact(pv * t);
204 return bbox;
205 }
209 static void
210 geom_line_wind_distance (Geom::Coord x0, Geom::Coord y0, Geom::Coord x1, Geom::Coord y1, Geom::Point const &pt, int *wind, Geom::Coord *best)
211 {
212 Geom::Coord Ax, Ay, Bx, By, Dx, Dy, s;
213 Geom::Coord dist2;
215 /* Find distance */
216 Ax = x0;
217 Ay = y0;
218 Bx = x1;
219 By = y1;
220 Dx = x1 - x0;
221 Dy = y1 - y0;
222 const Geom::Coord Px = pt[X];
223 const Geom::Coord Py = pt[Y];
225 if (best) {
226 s = ((Px - Ax) * Dx + (Py - Ay) * Dy) / (Dx * Dx + Dy * Dy);
227 if (s <= 0.0) {
228 dist2 = (Px - Ax) * (Px - Ax) + (Py - Ay) * (Py - Ay);
229 } else if (s >= 1.0) {
230 dist2 = (Px - Bx) * (Px - Bx) + (Py - By) * (Py - By);
231 } else {
232 Geom::Coord Qx, Qy;
233 Qx = Ax + s * Dx;
234 Qy = Ay + s * Dy;
235 dist2 = (Px - Qx) * (Px - Qx) + (Py - Qy) * (Py - Qy);
236 }
238 if (dist2 < (*best * *best)) *best = sqrt (dist2);
239 }
241 if (wind) {
242 /* Find wind */
243 if ((Ax >= Px) && (Bx >= Px)) return;
244 if ((Ay >= Py) && (By >= Py)) return;
245 if ((Ay < Py) && (By < Py)) return;
246 if (Ay == By) return;
247 /* Ctach upper y bound */
248 if (Ay == Py) {
249 if (Ax < Px) *wind -= 1;
250 return;
251 } else if (By == Py) {
252 if (Bx < Px) *wind += 1;
253 return;
254 } else {
255 Geom::Coord Qx;
256 /* Have to calculate intersection */
257 Qx = Ax + Dx * (Py - Ay) / Dy;
258 if (Qx < Px) {
259 *wind += (Dy > 0.0) ? 1 : -1;
260 }
261 }
262 }
263 }
265 static void
266 geom_cubic_bbox_wind_distance (Geom::Coord x000, Geom::Coord y000,
267 Geom::Coord x001, Geom::Coord y001,
268 Geom::Coord x011, Geom::Coord y011,
269 Geom::Coord x111, Geom::Coord y111,
270 Geom::Point const &pt,
271 Geom::Rect *bbox, int *wind, Geom::Coord *best,
272 Geom::Coord tolerance)
273 {
274 Geom::Coord x0, y0, x1, y1, len2;
275 int needdist, needwind, needline;
277 const Geom::Coord Px = pt[X];
278 const Geom::Coord Py = pt[Y];
280 needdist = 0;
281 needwind = 0;
282 needline = 0;
284 if (bbox) cubic_bbox (x000, y000, x001, y001, x011, y011, x111, y111, *bbox);
286 x0 = MIN (x000, x001);
287 x0 = MIN (x0, x011);
288 x0 = MIN (x0, x111);
289 y0 = MIN (y000, y001);
290 y0 = MIN (y0, y011);
291 y0 = MIN (y0, y111);
292 x1 = MAX (x000, x001);
293 x1 = MAX (x1, x011);
294 x1 = MAX (x1, x111);
295 y1 = MAX (y000, y001);
296 y1 = MAX (y1, y011);
297 y1 = MAX (y1, y111);
299 if (best) {
300 /* Quickly adjust to endpoints */
301 len2 = (x000 - Px) * (x000 - Px) + (y000 - Py) * (y000 - Py);
302 if (len2 < (*best * *best)) *best = (Geom::Coord) sqrt (len2);
303 len2 = (x111 - Px) * (x111 - Px) + (y111 - Py) * (y111 - Py);
304 if (len2 < (*best * *best)) *best = (Geom::Coord) sqrt (len2);
306 if (((x0 - Px) < *best) && ((y0 - Py) < *best) && ((Px - x1) < *best) && ((Py - y1) < *best)) {
307 /* Point is inside sloppy bbox */
308 /* Now we have to decide, whether subdivide */
309 /* fixme: (Lauris) */
310 if (((y1 - y0) > 5.0) || ((x1 - x0) > 5.0)) {
311 needdist = 1;
312 } else {
313 needline = 1;
314 }
315 }
316 }
317 if (!needdist && wind) {
318 if ((y1 >= Py) && (y0 < Py) && (x0 < Px)) {
319 /* Possible intersection at the left */
320 /* Now we have to decide, whether subdivide */
321 /* fixme: (Lauris) */
322 if (((y1 - y0) > 5.0) || ((x1 - x0) > 5.0)) {
323 needwind = 1;
324 } else {
325 needline = 1;
326 }
327 }
328 }
330 if (needdist || needwind) {
331 Geom::Coord x00t, x0tt, xttt, x1tt, x11t, x01t;
332 Geom::Coord y00t, y0tt, yttt, y1tt, y11t, y01t;
333 Geom::Coord s, t;
335 t = 0.5;
336 s = 1 - t;
338 x00t = s * x000 + t * x001;
339 x01t = s * x001 + t * x011;
340 x11t = s * x011 + t * x111;
341 x0tt = s * x00t + t * x01t;
342 x1tt = s * x01t + t * x11t;
343 xttt = s * x0tt + t * x1tt;
345 y00t = s * y000 + t * y001;
346 y01t = s * y001 + t * y011;
347 y11t = s * y011 + t * y111;
348 y0tt = s * y00t + t * y01t;
349 y1tt = s * y01t + t * y11t;
350 yttt = s * y0tt + t * y1tt;
352 geom_cubic_bbox_wind_distance (x000, y000, x00t, y00t, x0tt, y0tt, xttt, yttt, pt, NULL, wind, best, tolerance);
353 geom_cubic_bbox_wind_distance (xttt, yttt, x1tt, y1tt, x11t, y11t, x111, y111, pt, NULL, wind, best, tolerance);
354 } else if (1 || needline) {
355 geom_line_wind_distance (x000, y000, x111, y111, pt, wind, best);
356 }
357 }
359 static void
360 geom_curve_bbox_wind_distance(Geom::Curve const & c, Geom::Matrix const &m,
361 Geom::Point const &pt,
362 Geom::Rect *bbox, int *wind, Geom::Coord *dist,
363 Geom::Coord tolerance, Geom::Rect const *viewbox,
364 Geom::Point &p0) // pass p0 through as it represents the last endpoint added (the finalPoint of last curve)
365 {
366 if( dynamic_cast<Geom::LineSegment const*>(&c) ||
367 dynamic_cast<Geom::HLineSegment const*>(&c) ||
368 dynamic_cast<Geom::VLineSegment const*>(&c) )
369 {
370 Geom::Point pe = c.finalPoint() * m;
371 if (bbox) {
372 bbox->expandTo(pe);
373 }
374 if (dist || wind) {
375 if (wind) { // we need to pick fill, so do what we're told
376 geom_line_wind_distance (p0[X], p0[Y], pe[X], pe[Y], pt, wind, dist);
377 } else { // only stroke is being picked; skip this segment if it's totally outside the viewbox
378 Geom::Rect swept(p0, pe);
379 if (!viewbox || swept.intersects(*viewbox))
380 geom_line_wind_distance (p0[X], p0[Y], pe[X], pe[Y], pt, wind, dist);
381 }
382 }
383 p0 = pe;
384 }
385 else if(Geom::CubicBezier const *cubic_bezier = dynamic_cast<Geom::CubicBezier const *>(&c)) {
386 Geom::Point p1 = (*cubic_bezier)[1] * m;
387 Geom::Point p2 = (*cubic_bezier)[2] * m;
388 Geom::Point p3 = (*cubic_bezier)[3] * m;
390 // get approximate bbox from handles (convex hull property of beziers):
391 Geom::Rect swept(p0, p3);
392 swept.expandTo(p1);
393 swept.expandTo(p2);
395 if (!viewbox || swept.intersects(*viewbox)) { // we see this segment, so do full processing
396 geom_cubic_bbox_wind_distance ( p0[X], p0[Y],
397 p1[X], p1[Y],
398 p2[X], p2[Y],
399 p3[X], p3[Y],
400 pt,
401 bbox, wind, dist, tolerance);
402 } else {
403 if (wind) { // if we need fill, we can just pretend it's a straight line
404 geom_line_wind_distance (p0[X], p0[Y], p3[X], p3[Y], pt, wind, dist);
405 } else { // otherwise, skip it completely
406 }
407 }
408 p0 = p3;
409 } else {
410 //this case handles sbasis as well as all other curve types
411 Geom::Path sbasis_path = Geom::cubicbezierpath_from_sbasis(c.toSBasis(), 0.1);
413 //recurse to convert the new path resulting from the sbasis to svgd
414 for(Geom::Path::iterator iter = sbasis_path.begin(); iter != sbasis_path.end(); ++iter) {
415 geom_curve_bbox_wind_distance(*iter, m, pt, bbox, wind, dist, tolerance, viewbox, p0);
416 }
417 }
418 }
420 /* Calculates...
421 and returns ... in *wind and the distance to ... in *dist.
422 Returns bounding box in *bbox if bbox!=NULL.
423 */
424 void
425 pathv_matrix_point_bbox_wind_distance (Geom::PathVector const & pathv, Geom::Matrix const &m, Geom::Point const &pt,
426 Geom::Rect *bbox, int *wind, Geom::Coord *dist,
427 Geom::Coord tolerance, Geom::Rect const *viewbox)
428 {
429 if (pathv.empty()) {
430 if (wind) *wind = 0;
431 if (dist) *dist = NR_HUGE;
432 return;
433 }
435 // remember last point of last curve
436 Geom::Point p0(0,0);
438 // remembering the start of subpath
439 Geom::Point p_start(0,0);
440 bool start_set = false;
442 for (Geom::PathVector::const_iterator it = pathv.begin(); it != pathv.end(); ++it) {
444 if (start_set) { // this is a new subpath
445 if (wind && (p0 != p_start)) // for correct fill picking, each subpath must be closed
446 geom_line_wind_distance (p0[X], p0[Y], p_start[X], p_start[Y], pt, wind, dist);
447 }
448 p0 = it->initialPoint() * m;
449 p_start = p0;
450 start_set = true;
451 if (bbox) {
452 bbox->expandTo(p0);
453 }
455 // loop including closing segment if path is closed
456 for (Geom::Path::const_iterator cit = it->begin(); cit != it->end_default(); ++cit) {
457 geom_curve_bbox_wind_distance(*cit, m, pt, bbox, wind, dist, tolerance, viewbox, p0);
458 }
459 }
461 if (start_set) {
462 if (wind && (p0 != p_start)) // for correct picking, each subpath must be closed
463 geom_line_wind_distance (p0[X], p0[Y], p_start[X], p_start[Y], pt, wind, dist);
464 }
465 }
467 // temporary wrapper
468 void
469 pathv_matrix_point_bbox_wind_distance (Geom::PathVector const & pathv, NR::Matrix const &m, NR::Point const &pt,
470 NR::Rect *bbox, int *wind, NR::Coord *dist,
471 NR::Coord tolerance, NR::Rect const *viewbox)
472 {
473 Geom::Rect _bbox;
474 if (bbox)
475 _bbox = to_2geom(*bbox);
476 Geom::Coord _dist;
477 if (dist)
478 _dist = *dist;
479 Geom::Rect _viewbox;
480 if (viewbox)
481 _viewbox = to_2geom(*viewbox);
483 pathv_matrix_point_bbox_wind_distance( pathv, to_2geom(m), to_2geom(pt),
484 bbox ? &_bbox : NULL,
485 wind,
486 dist ? &_dist : NULL,
487 tolerance,
488 viewbox ? &_viewbox : NULL );
490 if (bbox)
491 *bbox = from_2geom(_bbox);
492 if (dist)
493 *dist = _dist;
494 }
495 //#################################################################################
497 /*
498 * Converts all segments in all paths to Geom::LineSegment or Geom::HLineSegment or
499 * Geom::VLineSegment or Geom::CubicBezier.
500 */
501 Geom::PathVector
502 pathv_to_linear_and_cubic_beziers( Geom::PathVector const &pathv )
503 {
504 Geom::PathVector output;
506 for (Geom::PathVector::const_iterator pit = pathv.begin(); pit != pathv.end(); ++pit) {
507 output.push_back( Geom::Path() );
508 output.back().start( pit->initialPoint() );
509 output.back().close( pit->closed() );
511 for (Geom::Path::const_iterator cit = pit->begin(); cit != pit->end_open(); ++cit) {
512 if( dynamic_cast<Geom::LineSegment const*>(&*cit) ||
513 dynamic_cast<Geom::HLineSegment const*>(&*cit) ||
514 dynamic_cast<Geom::VLineSegment const*>(&*cit) )
515 {
516 output.back().append(*cit);
517 }
518 else if(Geom::CubicBezier const *cubic_bezier = dynamic_cast<Geom::CubicBezier const*>(&*cit)) {
519 (void)cubic_bezier;
520 output.back().append(*cit);
521 }
522 else {
523 // convert all other curve types to cubicbeziers
524 Geom::Path cubicbezier_path = Geom::cubicbezierpath_from_sbasis(cit->toSBasis(), 0.1);
526 for(Geom::Path::iterator iter = cubicbezier_path.begin(); iter != cubicbezier_path.end(); ++iter) {
527 output.back().append(*iter);
528 }
529 }
530 }
531 }
533 return output;
534 }
539 /*
540 Local Variables:
541 mode:c++
542 c-file-style:"stroustrup"
543 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
544 indent-tabs-mode:nil
545 fill-column:99
546 End:
547 */
548 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :