synfig-core/trunk/src/synfig/curve_helper.cpp

   1 /* === S Y N F I G ========================================================= */
   2 /*!     \file curve_helper.cpp
   3 **      \brief Curve Helper File
   4 **
   5 **      $Id: curve_helper.cpp,v 1.1.1.1 2005/01/04 01:23:14 darco Exp $
   6 **
   7 **      \legal
   8 **      Copyright (c) 2002 Robert B. Quattlebaum Jr.
   9 **
  10 **      This software and associated documentation
  11 **      are CONFIDENTIAL and PROPRIETARY property of
  12 **      the above-mentioned copyright holder.
  13 **
  14 **      You may not copy, print, publish, or in any
  15 **      other way distribute this software without
  16 **      a prior written agreement with
  17 **      the copyright holder.
  18 **      \endlegal
  19 */
  20 /* ========================================================================= */
  21
  22 /* === H E A D E R S ======================================================= */
  23
  24 #ifdef USING_PCH
  25 #       include "pch.h"
  26 #else
  27 #ifdef HAVE_CONFIG_H
  28 #       include <config.h>
  29 #endif
  30
  31 #include "curve_helper.h"
  32
  33 #include <algorithm>
  34 #include <vector>
  35
  36 #endif
  37
  38 /* === U S I N G =========================================================== */
  39
  40 using namespace std;
  41 using namespace etl;
  42 using namespace synfig;
  43
  44 /* === M A C R O S ========================================================= */
  45 #define ERR     1e-11
  46 const Real ERROR = 1e-11;
  47
  48 /* === G L O B A L S ======================================================= */
  49
  50 /* === P R O C E D U R E S ================================================= */
  51
  52 /* === M E T H O D S ======================================================= */
  53
  54 /* === E N T R Y P O I N T ================================================= */
  55
  56 Real synfig::find_closest(const etl::bezier<Point> &curve, const Point &point,
  57                                 float step, Real *dout, float *tout)
  58 {
  59 #if 0
  60         float time(curve.find_closest(point,4));
  61         Real dist((curve(time)-point).mag());
  62         if(dout) *dout=dist;
  63         if(tout) *tout=time;
  64         return time;
  65 #else
  66         Real d,closest = 1.0e50;
  67         float t,time,closestt = -1;
  68         Vector p0,p1,end;
  69
  70         if(dout && *dout > 0)
  71                 closest = *dout;
  72
  73         p0 = curve[0];
  74         end = curve[3];
  75
  76         for(t = step; t < 1; t+=step, p0=p1)
  77         {
  78                 p1 = curve(t);
  79                 d = line_point_distsq(p0,p1,point,time);
  80
  81                 if(d<closest)
  82                 {
  83                         closest=d;
  84                         closestt = t-step + time*step;//t+(time-1)*step; //time between [t-step,t]
  85                 }
  86         }
  87
  88         d = line_point_distsq(p0,end,point,time);
  89         if(d<closest)
  90         {
  91                 closest = d;
  92                 closestt= t-step + time*(1-t+step); //time between [t-step,1.0]
  93         }
  94
  95         //set the time value if we found a closer point
  96         if(closestt >=0)
  97         {
  98                 if(tout) *tout = closestt;
  99         }
 100
 101         return closest;
 102 #endif
 103 }
 104
 105 // Line and BezHull Definitions
 106 void BezHull::Bound(const etl::bezier<Point> &b)
 107 {
 108         #if 1
 109
 110         //with a starting vertex, find the only vertex that has all other vertices on it's right
 111         int i,j;
 112         int first,cur,last;
 113
 114         float d,ds;
 115
 116         Vector n,vi;
 117         Vector::value_type      deqn;
 118
 119         //get left most vertex
 120         d = b[0][0];
 121         first = 0;
 122         for(i = 1; i < 4; ++i)
 123         {
 124                 if(b[i][0] < d)
 125                 {
 126                         d = b[i][0];
 127                         first = i;
 128                 }
 129         }
 130         cur = last = first;
 131         size = 0;
 132
 133         //find the farthest point with all points on right
 134         ds = 0;
 135         do //should reassign cur so it won't break on first step
 136         {
 137                 for(i = 0; i < 4; ++i)
 138                 {
 139                         if(i == cur || i == last) continue;
 140
 141                         //rotate vector to right to make normal
 142                         vi = -(b[i] - b[cur]).perp();
 143                         d = vi.mag_squared();
 144
 145                         //we want only the farthest (solves the case with many points on a line)
 146                         if(d > ds)
 147                         {
 148                                 ds = d;
 149                                 deqn = n*b[cur];
 150                                 for(j = 0; j < 4; ++j)
 151                                 {
 152                                         d = n*b[i] - deqn;
 153                                         if(d < 0) break; //we're on left, nope!
 154                                 }
 155
 156                                 //everyone is on right... yay! :)
 157                                 if(d >= 0)
 158                                 {
 159                                         //advance point and add last one into hull
 160                                         p[size++] = p[last];
 161                                         last = cur;
 162                                         cur = i;
 163                                 }
 164                         }
 165                 }
 166         }while(cur != first);
 167
 168         #else
 169
 170         //will work but does not keep winding order
 171
 172         //convex hull alg.
 173         //build set of line segs which have no points on other side...
 174         //start with initial normal segments
 175
 176         //start with single triangle
 177         p[0] = b[0];
 178         p[1] = b[1];
 179         p[2] = b[2];
 180         p[3] = b[3];
 181
 182         //initial reject (if point is inside triangle don't care)
 183         {
 184                 Vector v1,v2,vp;
 185
 186                 v1 = p[1]-p[0];
 187                 v2 = p[2]-p[0];
 188
 189                 vp = p[3]-p[0];
 190
 191                 float   s = (vp*v1) / (v1*v1),
 192                                 t = (vp*v2) / (v2*v2);
 193
 194                 //if we're inside the triangle we don't this sissy point
 195                 if( s >= 0 && s <= 1 && t >= 0 && t <= 1 )
 196                 {
 197                         size = 3;
 198                         return;
 199                 }
 200         }
 201
 202         //expand triangle based on info...
 203         bool line;
 204         int index,i,j;
 205         float ds,d;
 206
 207         //distance from point to vertices
 208         line = false;
 209         index = 0;
 210         ds = (p[0]-b[3]).mag_squared();
 211         for(i = 1; i < 3; ++i)
 212         {
 213                 d = (p[3]-p[i]).mag_squared();
 214                 if(d < ds)
 215                 {
 216                         index = i;
 217                         ds = d;
 218                 }
 219         }
 220
 221         //distance to line
 222         float t;
 223         j = 2;
 224         for(i = 0; i < 3; j = i++)
 225         {
 226                 d = line_point_distsq(p[j],p[i],b[4],t);
 227                 if(d < ds)
 228                 {
 229                         index = j;
 230                         ds = d;
 231                         line = true;
 232                 }
 233         }
 234
 235         //We don't need no stinkin extra vertex, just replace
 236         if(!line)
 237         {
 238                 p[index] = p[3];
 239                 size = 3;
 240         }else
 241         {
 242                 //must expand volume to work with point...
 243                 //      after the index then
 244
 245                 /* Pattern:
 246                         0 - push 1,2 -> 2,3
 247                         1 - push 2 -> 3
 248                         2 - none
 249                 */
 250                 for(i = 3; i > index+1; --i)
 251                 {
 252                         p[i] = p[i-1];
 253                 }
 254
 255                 p[index] = b[3]; //recopy b3
 256                 size = 4;
 257         }
 258
 259         #endif
 260 }
 261
 262 //Line Intersection
 263 int
 264 synfig::intersect(const Point &p1, const Vector &v1, float &t1,
 265                                         const Point &p2, const Vector &v2, float &t2)
 266 {
 267         /* Parametric intersection:
 268                 l1 = p1 + tv1, l2 = p2 + sv2
 269
 270                 0 = p1+tv1-(p2+sv2)
 271                 group parameters: sv2 - tv1 = p1-p2
 272
 273                 ^ = transpose
 274                 invert matrix (on condition det != 0):
 275                 A[t s]^ = [p1-p2]^
 276
 277                 A = [-v1 v2]
 278
 279                 det = v1y.v2x - v1x.v2y
 280
 281                 if non 0 then A^-1 = invdet * | v2y -v2x |
 282                                                                           | v1y -v1x |
 283
 284                 [t s]^ = A^-1 [p1-p2]^
 285         */
 286
 287         Vector::value_type det = v1[1]*v2[0] - v1[0]*v2[1];
 288
 289         //is determinant valid?
 290         if(det > ERR || det < -ERR)
 291         {
 292                 Vector p_p = p1-p2;
 293
 294                 det = 1/det;
 295
 296                 t1 = det*(v2[1]*p_p[0] - v2[0]*p_p[1]);
 297                 t2 = det*(v1[1]*p_p[0] - v1[0]*p_p[1]);
 298
 299                 return 1;
 300         }
 301
 302         return 0;
 303 }
 304
 305 //Returns the true or false intersection of a rectangle and a line
 306 int intersect(const Rect &r, const Point &p, const Vector &v)
 307 {
 308         float t[4] = {0};
 309
 310         /*get horizontal intersections and then vertical intersections
 311                 and intersect them
 312
 313                 Vertical planes - n = (1,0)
 314                 Horizontal planes - n = (0,1)
 315
 316                 so if we are solving for ray with implicit line
 317         */
 318
 319         //solve horizontal
 320         if(v[0] > ERR || v[0] < -ERR)
 321         {
 322                 //solve for t0, t1
 323                 t[0] = (r.minx - p[0])/v[0];
 324                 t[1] = (r.maxx - p[0])/v[0];
 325         }else
 326         {
 327                 return (int)(p[1] >= r.miny && p[1] <= r.maxy);
 328         }
 329
 330         //solve vertical
 331         if(v[1] > ERR || v[1] < -ERR)
 332         {
 333                 //solve for t0, t1
 334                 t[2] = (r.miny - p[1])/v[1];
 335                 t[3] = (r.maxy - p[1])/v[1];
 336         }else
 337         {
 338                 return (int)(p[0] >= r.minx && p[0] <= r.maxx);
 339         }
 340
 341         return (int)(t[0] <= t[3] && t[1] >= t[2]);
 342 }
 343
 344 int synfig::intersect(const Rect &r, const Point &p)
 345 {
 346         return (p[1] < r.maxy && p[1] > r.miny) && p[0] > r.minx;
 347 }
 348
 349 //returns 0 or 1 for true or false number of intersections of a ray with a bezier convex hull
 350 int intersect(const BezHull &bh, const Point &p, const Vector &v)
 351 {
 352         float mint = 0, maxt = 1e20;
 353
 354         //polygon cliping
 355         Vector n;
 356         Vector::value_type      nv;
 357
 358         Point last = bh.p[3];
 359         for(int i = 0; i < bh.size; ++i)
 360         {
 361                 n = (bh.p[i] - last).perp(); //rotate 90 deg.
 362
 363                 /*
 364                         since rotated left
 365                         if n.v  < 0 - going in
 366                                         > 0 - going out
 367                                         = 0 - parallel
 368                 */
 369                 nv = n*v;
 370
 371                 //going OUT
 372                 if(nv > ERR)
 373                 {
 374                         maxt = min(maxt,(float)((n*(p-last))/nv));
 375                 }else
 376                 if( nv < -ERR) //going IN
 377                 {
 378                         mint = max(mint,(float)((n*(p-last))/nv));
 379                 }else
 380                 {
 381                         if( n*(p-last) > 0 ) //outside entirely
 382                         {
 383                                 return 0;
 384                         }
 385                 }
 386
 387                 last = bh.p[i];
 388         }
 389
 390         return 0;
 391 }
 392
 393 int Clip(const Rect &r, const Point &p1, const Point &p2, Point *op1, Point *op2)
 394 {
 395         float t1=0,t2=1;
 396         Vector v=p2-p1;
 397
 398         /*get horizontal intersections and then vertical intersections
 399                 and intersect them
 400
 401                 Vertical planes - n = (1,0)
 402                 Horizontal planes - n = (0,1)
 403
 404                 so if we are solving for ray with implicit line
 405         */
 406
 407         //solve horizontal
 408         if(v[0] > ERR || v[0] < -ERR)
 409         {
 410                 //solve for t0, t1
 411                 float   tt1 = (r.minx - p1[0])/v[0],
 412                                 tt2 = (r.maxx - p1[0])/v[0];
 413
 414                 //line in positive direction (normal comparisons
 415                 if(tt1 < tt2)
 416                 {
 417                         t1 = max(t1,tt1);
 418                         t2 = min(t2,tt2);
 419                 }else
 420                 {
 421                         t1 = max(t1,tt2);
 422                         t2 = min(t2,tt1);
 423                 }
 424         }else
 425         {
 426                 if(p1[1] < r.miny || p1[1] > r.maxy)
 427                         return 0;
 428         }
 429
 430         //solve vertical
 431         if(v[1] > ERR || v[1] < -ERR)
 432         {
 433                 //solve for t0, t1
 434                 float   tt1 = (r.miny - p1[1])/v[1],
 435                                 tt2 = (r.maxy - p1[1])/v[1];
 436
 437                 //line in positive direction (normal comparisons
 438                 if(tt1 < tt2)
 439                 {
 440                         t1 = max(t1,tt1);
 441                         t2 = min(t2,tt2);
 442                 }else
 443                 {
 444                         t1 = max(t1,tt2);
 445                         t2 = min(t2,tt1);
 446                 }
 447         }else
 448         {
 449                 if(p1[0] < r.minx || p1[0] > r.maxx)
 450                         return 0;
 451         }
 452
 453         if(op1) *op1 = p1 + v*t1;
 454         if(op2) *op2 = p1 + v*t2;
 455
 456         return 1;
 457 }
 458
 459 static void clean_bez(const bezier<Point> &b, bezier<Point> &out)
 460 {
 461         bezier<Point> temp;
 462
 463         temp = b;
 464         temp.set_r(0);
 465         temp.set_s(1);
 466
 467         if(b.get_r() != 0)
 468                 temp.subdivide(0,&temp,b.get_r());
 469
 470         if(b.get_s() != 1)
 471                 temp.subdivide(&temp,0,b.get_s());
 472
 473         out = temp;
 474 }
 475
 476 // CIntersect Definitions
 477
 478 CIntersect::CIntersect()
 479         : max_depth(10) //depth of 10 means timevalue parameters will have an approx. error bound of 2^-10
 480 {
 481 }
 482
 483 struct CIntersect::SCurve
 484 {
 485         bezier<Point>   b;              //the current subdivided curve
 486         float rt,st;
 487         //float                         mid,    //the midpoint time value on this section of the subdivided curve
 488         //                              scale;  //the current delta in time values this curve would be on original curve
 489
 490         float   mag;                    //approximate sum of magnitudes of each edge of control polygon
 491         Rect    aabb;                   //Axis Aligned Bounding Box for quick (albeit less accurate) collision
 492
 493         SCurve() {}
 494
 495         SCurve(const bezier<Point> &c,float rin, float sin)
 496         :b(c),rt(rin),st(sin),mag(1)
 497         {
 498                 Bound(aabb,b);
 499         }
 500
 501         void Split(SCurve &l, SCurve &r) const
 502         {
 503                 b.subdivide(&l.b,&r.b);
 504
 505                 l.rt = rt;
 506                 r.st = st;
 507                 l.st = r.rt = (rt+st)/2;
 508
 509                 Bound(l.aabb,l.b);
 510                 Bound(r.aabb,r.b);
 511         }
 512 };
 513
 514 //Curve to the left of point test
 515 static int recurse_intersect(const CIntersect::SCurve &b, const Point &p1, int depthleft = 10)
 516 {
 517         //reject when the line does not intersect the bounding box
 518         if(!intersect(b.aabb,p1)) return 0;
 519
 520         //accept curves (and perform super detailed check for intersections)
 521         //if the values are below tolerance
 522
 523         //NOTE FOR BETTERING OF ALGORITHM: SHOULD ALSO/IN-PLACE-OF CHECK MAGNITUDE OF EDGES (or approximate)
 524         if(depthleft <= 0)
 525         {
 526                 //NOTE FOR IMPROVEMENT: Polish roots based on original curve
 527                 //                                              (may be too expensive to be effective)
 528                 int turn = 0;
 529
 530                 for(int i = 0; i < 3; ++i)
 531                 {
 532                         //intersect line segmentsssss
 533
 534                         //solve for the y_value
 535                         Vector v = b.b[i+1] - b.b[i];
 536
 537                         if(v[1] > ERROR && v[1] < ERROR)
 538                         {
 539                                 Real xi = (p1[1] - b.b[i][1])/v[1];
 540
 541                                 //and add in the turn (up or down) if it's valid
 542                                 if(xi < p1[0]) turn += (v[1] > 0) ? 1 : -1;
 543                         }
 544                 }
 545
 546                 return turn;
 547         }
 548
 549         //subdivide the curve and continue
 550         CIntersect::SCurve l1,r1;
 551         b.Split(l1,r1); //subdivide left
 552
 553         //test each subdivision against the point
 554         return recurse_intersect(l1,p1) + recurse_intersect(r1,p1);
 555 }
 556
 557 int intersect(const bezier<Point> &b, const Point &p)
 558 {
 559         CIntersect::SCurve      sb;
 560         clean_bez(b,sb.b);
 561
 562         sb.rt = 0; sb.st = 1;
 563         sb.mag = 1; Bound(sb.aabb,sb.b);
 564
 565         return recurse_intersect(sb,p);
 566 }
 567
 568 //Curve curve intersection
 569 void CIntersect::recurse_intersect(const SCurve &left, const SCurve &right, int depth)
 570 {
 571         //reject curves that do not overlap with bouding boxes
 572         if(!intersect(left.aabb,right.aabb)) return;
 573
 574         //accept curves (and perform super detailed check for intersections)
 575         //if the values are below tolerance
 576
 577         //NOTE FOR BETTERING OF ALGORITHM: SHOULD ALSO/IN-PLACE-OF CHECK MAGNITUDE OF EDGES (or approximate)
 578         if(depth >= max_depth)
 579         {
 580                 //NOTE FOR IMPROVEMENT: Polish roots based on original curve with the Jacobian
 581                 //                                              (may be too expensive to be effective)
 582
 583                 //perform root approximation
 584                 //collide line segments
 585
 586                 float t,s;
 587
 588                 for(int i = 0; i < 3; ++i)
 589                 {
 590                         for(int j = 0; j < 3; ++j)
 591                         {
 592                                 //intersect line segmentsssss
 593                                 if(intersect_line_segments(left.b[i],left.b[i+1],t,right.b[j],right.b[j+1],s))
 594                                 {
 595                                         //We got one Jimmy
 596                                         times.push_back(intersect_set::value_type(t,s));
 597                                 }
 598                         }
 599                 }
 600
 601                 return;
 602         }
 603
 604         //NOTE FOR IMPROVEMENT: only subdivide one curve and choose the one that has
 605         //                                              the highest approximated length
 606         //fast approximation to curve length may be hard (accurate would
 607         // involve 3 square roots), could sum the squares which would be
 608         // quick but inaccurate
 609
 610         SCurve l1,r1,l2,r2;
 611         left.Split(l1,r1);      //subdivide left
 612         right.Split(l2,r2); //subdivide right
 613
 614         //Test each cantidate against eachother
 615         recurse_intersect(l1,l2);
 616         recurse_intersect(l1,r2);
 617         recurse_intersect(r1,l2);
 618         recurse_intersect(r1,r2);
 619 }
 620
 621
 622
 623 bool CIntersect::operator()(const bezier<Point> &c1, const bezier<Point> &c2)
 624 {
 625         times.clear();
 626
 627         //need to subdivide and check recursive bounding regions against eachother
 628         //so track a list of dirty curves and compare compare compare
 629
 630
 631         //temporary curves for subdivision
 632         CIntersect                      intersector;
 633         CIntersect::SCurve      left,right;
 634
 635         //Make sure the parameters are normalized (so we don't compare unwanted parts of the curves,
 636         //      and don't miss any for that matter)
 637
 638         //left curve
 639         //Compile information about curve
 640         clean_bez(c1,left.b);
 641         left.rt = 0; left.st = 1;
 642         Bound(left.aabb, left.b);
 643
 644         //right curve
 645         //Compile information about right curve
 646         clean_bez(c2,right.b);
 647         right.rt = 0; right.st = 1;
 648         Bound(right.aabb, right.b);
 649
 650         //Perform Curve intersection
 651         intersector.recurse_intersect(left,right);
 652
 653         //Get information about roots (yay! :P)
 654         return times.size() != 0;
 655 }
 656
 657 //point inside curve - return +/- hit up or down edge
 658 int intersect_scurve(const CIntersect::SCurve &b, const Point &p)
 659 {
 660         //initial reject/approve etc.
 661
 662         /*
 663                         *-----------*---------
 664                         |                       |
 665                         |                       |
 666                         |                       |
 667                         |         1             |    2
 668                         |                       |
 669                         |                       |
 670                         |                       |
 671                         |                       |
 672                         *-----------*--------
 673                 1,2 are only regions not rejected
 674         */
 675         if(p[0] < b.aabb.minx || p[1] < b.aabb.miny || p[1] > b.aabb.maxy)
 676                 return 0;
 677
 678         //approve only if to the right of rect around 2 end points
 679         {
 680                 Rect    r;
 681                 r.set_point(b.b[0][0],b.b[0][1]);
 682                 r.expand(b.b[3][0],b.b[3][1]);
 683
 684                 if(p[0] >= r.maxx && p[1] <= r.maxy && p[1] >= r.miny)
 685                 {
 686                         float df = b.b[3][1] - b.b[0][1];
 687
 688                         return df >= 0 ? 1 : -1;
 689                 }
 690         }
 691
 692         //subdivide and check again!
 693         CIntersect::SCurve      l,r;
 694         b.Split(l,r);
 695         return  intersect_scurve(l,p) + intersect_scurve(r,p);
 696 }
 697
 698 int synfig::intersect(const bezier<Point> &b, const Point &p)
 699 {
 700         CIntersect::SCurve      c(b,0,1);
 701
 702         return intersect_scurve(c,p);
 703 }