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