1 /* === S Y N F I G ========================================================= */
2 /*! \file curve_helper.cpp
3 ** \brief Curve Helper File
5 ** $Id: curve_helper.cpp,v 1.1.1.1 2005/01/04 01:23:14 darco Exp $
8 ** Copyright (c) 2002-2005 Robert B. Quattlebaum Jr., Adrian Bentley
10 ** This package is free software; you can redistribute it and/or
11 ** modify it under the terms of the GNU General Public License as
12 ** published by the Free Software Foundation; either version 2 of
13 ** the License, or (at your option) any later version.
15 ** This package is distributed in the hope that it will be useful,
16 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
17 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 ** General Public License for more details.
21 /* ========================================================================= */
23 /* === H E A D E R S ======================================================= */
32 #include "curve_helper.h"
39 /* === U S I N G =========================================================== */
43 using namespace synfig;
45 /* === M A C R O S ========================================================= */
47 const Real ERROR = 1e-11;
49 /* === G L O B A L S ======================================================= */
51 /* === P R O C E D U R E S ================================================= */
53 /* === M E T H O D S ======================================================= */
55 /* === E N T R Y P O I N T ================================================= */
57 Real synfig::find_closest(const etl::bezier<Point> &curve, const Point &point,
58 float step, Real *dout, float *tout)
61 float time(curve.find_closest(point,4));
62 Real dist((curve(time)-point).mag());
67 Real d,closest = 1.0e50;
68 float t,time,closestt = -1;
77 for(t = step; t < 1; t+=step, p0=p1)
80 d = line_point_distsq(p0,p1,point,time);
85 closestt = t-step + time*step;//t+(time-1)*step; //time between [t-step,t]
89 d = line_point_distsq(p0,end,point,time);
93 closestt= t-step + time*(1-t+step); //time between [t-step,1.0]
96 //set the time value if we found a closer point
99 if(tout) *tout = closestt;
106 // Line and BezHull Definitions
107 void BezHull::Bound(const etl::bezier<Point> &b)
111 //with a starting vertex, find the only vertex that has all other vertices on it's right
118 Vector::value_type deqn;
120 //get left most vertex
123 for(i = 1; i < 4; ++i)
134 //find the farthest point with all points on right
136 do //should reassign cur so it won't break on first step
138 for(i = 0; i < 4; ++i)
140 if(i == cur || i == last) continue;
142 //rotate vector to right to make normal
143 vi = -(b[i] - b[cur]).perp();
144 d = vi.mag_squared();
146 //we want only the farthest (solves the case with many points on a line)
151 for(j = 0; j < 4; ++j)
154 if(d < 0) break; //we're on left, nope!
157 //everyone is on right... yay! :)
160 //advance point and add last one into hull
167 }while(cur != first);
171 //will work but does not keep winding order
174 //build set of line segs which have no points on other side...
175 //start with initial normal segments
177 //start with single triangle
183 //initial reject (if point is inside triangle don't care)
192 float s = (vp*v1) / (v1*v1),
193 t = (vp*v2) / (v2*v2);
195 //if we're inside the triangle we don't this sissy point
196 if( s >= 0 && s <= 1 && t >= 0 && t <= 1 )
203 //expand triangle based on info...
208 //distance from point to vertices
211 ds = (p[0]-b[3]).mag_squared();
212 for(i = 1; i < 3; ++i)
214 d = (p[3]-p[i]).mag_squared();
225 for(i = 0; i < 3; j = i++)
227 d = line_point_distsq(p[j],p[i],b[4],t);
236 //We don't need no stinkin extra vertex, just replace
243 //must expand volume to work with point...
244 // after the index then
251 for(i = 3; i > index+1; --i)
256 p[index] = b[3]; //recopy b3
265 synfig::intersect(const Point &p1, const Vector &v1, float &t1,
266 const Point &p2, const Vector &v2, float &t2)
268 /* Parametric intersection:
269 l1 = p1 + tv1, l2 = p2 + sv2
272 group parameters: sv2 - tv1 = p1-p2
275 invert matrix (on condition det != 0):
280 det = v1y.v2x - v1x.v2y
282 if non 0 then A^-1 = invdet * | v2y -v2x |
285 [t s]^ = A^-1 [p1-p2]^
288 Vector::value_type det = v1[1]*v2[0] - v1[0]*v2[1];
290 //is determinant valid?
291 if(det > ERR || det < -ERR)
297 t1 = det*(v2[1]*p_p[0] - v2[0]*p_p[1]);
298 t2 = det*(v1[1]*p_p[0] - v1[0]*p_p[1]);
306 //Returns the true or false intersection of a rectangle and a line
307 int intersect(const Rect &r, const Point &p, const Vector &v)
311 /*get horizontal intersections and then vertical intersections
314 Vertical planes - n = (1,0)
315 Horizontal planes - n = (0,1)
317 so if we are solving for ray with implicit line
321 if(v[0] > ERR || v[0] < -ERR)
324 t[0] = (r.minx - p[0])/v[0];
325 t[1] = (r.maxx - p[0])/v[0];
328 return (int)(p[1] >= r.miny && p[1] <= r.maxy);
332 if(v[1] > ERR || v[1] < -ERR)
335 t[2] = (r.miny - p[1])/v[1];
336 t[3] = (r.maxy - p[1])/v[1];
339 return (int)(p[0] >= r.minx && p[0] <= r.maxx);
342 return (int)(t[0] <= t[3] && t[1] >= t[2]);
345 int synfig::intersect(const Rect &r, const Point &p)
347 return (p[1] < r.maxy && p[1] > r.miny) && p[0] > r.minx;
350 //returns 0 or 1 for true or false number of intersections of a ray with a bezier convex hull
351 int intersect(const BezHull &bh, const Point &p, const Vector &v)
353 float mint = 0, maxt = 1e20;
357 Vector::value_type nv;
359 Point last = bh.p[3];
360 for(int i = 0; i < bh.size; ++i)
362 n = (bh.p[i] - last).perp(); //rotate 90 deg.
366 if n.v < 0 - going in
375 maxt = min(maxt,(float)((n*(p-last))/nv));
377 if( nv < -ERR) //going IN
379 mint = max(mint,(float)((n*(p-last))/nv));
382 if( n*(p-last) > 0 ) //outside entirely
394 int Clip(const Rect &r, const Point &p1, const Point &p2, Point *op1, Point *op2)
399 /*get horizontal intersections and then vertical intersections
402 Vertical planes - n = (1,0)
403 Horizontal planes - n = (0,1)
405 so if we are solving for ray with implicit line
409 if(v[0] > ERR || v[0] < -ERR)
412 float tt1 = (r.minx - p1[0])/v[0],
413 tt2 = (r.maxx - p1[0])/v[0];
415 //line in positive direction (normal comparisons
427 if(p1[1] < r.miny || p1[1] > r.maxy)
432 if(v[1] > ERR || v[1] < -ERR)
435 float tt1 = (r.miny - p1[1])/v[1],
436 tt2 = (r.maxy - p1[1])/v[1];
438 //line in positive direction (normal comparisons
450 if(p1[0] < r.minx || p1[0] > r.maxx)
454 if(op1) *op1 = p1 + v*t1;
455 if(op2) *op2 = p1 + v*t2;
460 static void clean_bez(const bezier<Point> &b, bezier<Point> &out)
469 temp.subdivide(0,&temp,b.get_r());
472 temp.subdivide(&temp,0,b.get_s());
477 // CIntersect Definitions
479 CIntersect::CIntersect()
480 : max_depth(10) //depth of 10 means timevalue parameters will have an approx. error bound of 2^-10
484 struct CIntersect::SCurve
486 bezier<Point> b; //the current subdivided curve
488 //float mid, //the midpoint time value on this section of the subdivided curve
489 // scale; //the current delta in time values this curve would be on original curve
491 float mag; //approximate sum of magnitudes of each edge of control polygon
492 Rect aabb; //Axis Aligned Bounding Box for quick (albeit less accurate) collision
496 SCurve(const bezier<Point> &c,float rin, float sin)
497 :b(c),rt(rin),st(sin),mag(1)
502 void Split(SCurve &l, SCurve &r) const
504 b.subdivide(&l.b,&r.b);
508 l.st = r.rt = (rt+st)/2;
515 //Curve to the left of point test
516 static int recurse_intersect(const CIntersect::SCurve &b, const Point &p1, int depthleft = 10)
518 //reject when the line does not intersect the bounding box
519 if(!intersect(b.aabb,p1)) return 0;
521 //accept curves (and perform super detailed check for intersections)
522 //if the values are below tolerance
524 //NOTE FOR BETTERING OF ALGORITHM: SHOULD ALSO/IN-PLACE-OF CHECK MAGNITUDE OF EDGES (or approximate)
527 //NOTE FOR IMPROVEMENT: Polish roots based on original curve
528 // (may be too expensive to be effective)
531 for(int i = 0; i < 3; ++i)
533 //intersect line segmentsssss
535 //solve for the y_value
536 Vector v = b.b[i+1] - b.b[i];
538 if(v[1] > ERROR && v[1] < ERROR)
540 Real xi = (p1[1] - b.b[i][1])/v[1];
542 //and add in the turn (up or down) if it's valid
543 if(xi < p1[0]) turn += (v[1] > 0) ? 1 : -1;
550 //subdivide the curve and continue
551 CIntersect::SCurve l1,r1;
552 b.Split(l1,r1); //subdivide left
554 //test each subdivision against the point
555 return recurse_intersect(l1,p1) + recurse_intersect(r1,p1);
558 int intersect(const bezier<Point> &b, const Point &p)
560 CIntersect::SCurve sb;
563 sb.rt = 0; sb.st = 1;
564 sb.mag = 1; Bound(sb.aabb,sb.b);
566 return recurse_intersect(sb,p);
569 //Curve curve intersection
570 void CIntersect::recurse_intersect(const SCurve &left, const SCurve &right, int depth)
572 //reject curves that do not overlap with bouding boxes
573 if(!intersect(left.aabb,right.aabb)) return;
575 //accept curves (and perform super detailed check for intersections)
576 //if the values are below tolerance
578 //NOTE FOR BETTERING OF ALGORITHM: SHOULD ALSO/IN-PLACE-OF CHECK MAGNITUDE OF EDGES (or approximate)
579 if(depth >= max_depth)
581 //NOTE FOR IMPROVEMENT: Polish roots based on original curve with the Jacobian
582 // (may be too expensive to be effective)
584 //perform root approximation
585 //collide line segments
589 for(int i = 0; i < 3; ++i)
591 for(int j = 0; j < 3; ++j)
593 //intersect line segmentsssss
594 if(intersect_line_segments(left.b[i],left.b[i+1],t,right.b[j],right.b[j+1],s))
597 times.push_back(intersect_set::value_type(t,s));
605 //NOTE FOR IMPROVEMENT: only subdivide one curve and choose the one that has
606 // the highest approximated length
607 //fast approximation to curve length may be hard (accurate would
608 // involve 3 square roots), could sum the squares which would be
609 // quick but inaccurate
612 left.Split(l1,r1); //subdivide left
613 right.Split(l2,r2); //subdivide right
615 //Test each cantidate against eachother
616 recurse_intersect(l1,l2);
617 recurse_intersect(l1,r2);
618 recurse_intersect(r1,l2);
619 recurse_intersect(r1,r2);
624 bool CIntersect::operator()(const bezier<Point> &c1, const bezier<Point> &c2)
628 //need to subdivide and check recursive bounding regions against eachother
629 //so track a list of dirty curves and compare compare compare
632 //temporary curves for subdivision
633 CIntersect intersector;
634 CIntersect::SCurve left,right;
636 //Make sure the parameters are normalized (so we don't compare unwanted parts of the curves,
637 // and don't miss any for that matter)
640 //Compile information about curve
641 clean_bez(c1,left.b);
642 left.rt = 0; left.st = 1;
643 Bound(left.aabb, left.b);
646 //Compile information about right curve
647 clean_bez(c2,right.b);
648 right.rt = 0; right.st = 1;
649 Bound(right.aabb, right.b);
651 //Perform Curve intersection
652 intersector.recurse_intersect(left,right);
654 //Get information about roots (yay! :P)
655 return times.size() != 0;
658 //point inside curve - return +/- hit up or down edge
659 int intersect_scurve(const CIntersect::SCurve &b, const Point &p)
661 //initial reject/approve etc.
664 *-----------*---------
673 *-----------*--------
674 1,2 are only regions not rejected
676 if(p[0] < b.aabb.minx || p[1] < b.aabb.miny || p[1] > b.aabb.maxy)
679 //approve only if to the right of rect around 2 end points
682 r.set_point(b.b[0][0],b.b[0][1]);
683 r.expand(b.b[3][0],b.b[3][1]);
685 if(p[0] >= r.maxx && p[1] <= r.maxy && p[1] >= r.miny)
687 float df = b.b[3][1] - b.b[0][1];
689 return df >= 0 ? 1 : -1;
693 //subdivide and check again!
694 CIntersect::SCurve l,r;
696 return intersect_scurve(l,p) + intersect_scurve(r,p);
699 int synfig::intersect(const bezier<Point> &b, const Point &p)
701 CIntersect::SCurve c(b,0,1);
703 return intersect_scurve(c,p);