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