Code

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