synfig-studio/tags/stable/src/synfigapp/blineconvert.cpp

   1 /* === S Y N F I G ========================================================= */
   2 /*!     \file blineconvert.cpp
   3 **      \brief Template File
   4 **
   5 **      $Id$
   6 **
   7 **      \legal
   8 **      Copyright (c) 2002-2005 Robert B. Quattlebaum Jr., Adrian Bentley
   9 **      Copyright (c) 2007 Chris Moore
  10 **
  11 **      This package is free software; you can redistribute it and/or
  12 **      modify it under the terms of the GNU General Public License as
  13 **      published by the Free Software Foundation; either version 2 of
  14 **      the License, or (at your option) any later version.
  15 **
  16 **      This package is distributed in the hope that it will be useful,
  17 **      but WITHOUT ANY WARRANTY; without even the implied warranty of
  18 **      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  19 **      General Public License for more details.
  20 **      \endlegal
  21 */
  22 /* ========================================================================= */
  23
  24 /* === H E A D E R S ======================================================= */
  25
  26 #ifdef USING_PCH
  27 #       include "pch.h"
  28 #else
  29 #ifdef HAVE_CONFIG_H
  30 #       include <config.h>
  31 #endif
  32
  33 #include "blineconvert.h"
  34 #include <vector>
  35 #include <ETL/gaussian>
  36 #include <ETL/hermite>
  37 #include <ETL/clock>
  38 #include <float.h>
  39 #include <algorithm>
  40 #include <synfig/general.h>
  41 #include <cassert>
  42
  43 #include "general.h"
  44
  45 #endif
  46
  47 /* === U S I N G =========================================================== */
  48
  49 using namespace std;
  50 using namespace etl;
  51 using namespace synfig;
  52
  53 /* === M A C R O S ========================================================= */
  54
  55 #define EPSILON         (1e-10)
  56
  57 /* === G L O B A L S ======================================================= */
  58
  59 /* === P R O C E D U R E S ================================================= */
  60
  61 /* === M E T H O D S ======================================================= */
  62
  63
  64 //Derivative Functions for numerical approximation
  65
  66 //bias == 0 will get F' at f3, bias < 0 will get F' at f1, and bias > 0 will get F' at f5
  67 template < class T >
  68 inline void FivePointdt(T &df, const T &f1, const T &f2, const T &f3, const T &f4, const T &f5, int bias)
  69 {
  70         if (bias == 0)                          // middle
  71                 df = (f1 - f2*8 + f4*8 - f5)*(1/12.0f);
  72         else if (bias < 0)                      // left
  73                 df = (-f1*25 + f2*48 - f3*36 + f4*16 - f5*3)*(1/12.0f);
  74         else                                            // right
  75                 df = (f1*3 - f2*16 + f3*36 - f4*48 + f5*25)*(1/12.0f);
  76 }
  77
  78 template < class T >
  79 inline void ThreePointdt(T &df, const T &f1, const T &f2, const T &f3, int bias)
  80 {
  81         if (bias == 0)                          // middle
  82                 df = (-f1 + f3)*(1/2.0f);
  83         else if (bias < 0)                      // left
  84                 df = (-f1*3 + f2*4 - f3)*(1/2.0f);
  85         else                                            // right
  86                 df = (f1 - f2*4 + f3*3)*(1/2.0f);
  87 }
  88
  89 // template < class T >
  90 // inline void ThreePointddt(T &df, const T &f1, const T &f2, const T &f3, int bias)
  91 // {
  92 //      // a 3 point approximation pretends to have constant acceleration,
  93 //      // so only one algorithm needed for left, middle, or right
  94 //      df = (f1 -f2*2 + f3)*(1/2.0f);
  95 // }
  96 //
  97 // // WARNING -- totally broken
  98 // template < class T >
  99 // inline void FivePointddt(T &df, const T &f1, const T &f2, const T &f3, int bias)
 100 // {
 101 //      if(bias == 0)
 102 //      {
 103 //              assert(0); // !?
 104 //              //middle
 105 //              //df = (- f1 + f2*16 - f3*30 +  f4*16 - f5)*(1/12.0f);
 106 //      }/*else if(bias < 0)
 107 //      {
 108 //              //left
 109 //              df = (f1*7 - f2*26*4 + f3*19*6 - f4*14*4 + f5*11)*(1/12.0f);
 110 //      }else
 111 //      {
 112 //              //right
 113 //              df = (f1*3 - f2*16 + f3*36 - f4*48 + f5*25)*(1/12.0f);
 114 //      }*/
 115 //      //side ones don't work, use 3 point
 116 // }
 117 //
 118 // //implement an arbitrary derivative
 119 // //dumb algorithm
 120 // template < class T >
 121 // void DerivativeApprox(T &df, const T f[], const Real t[], int npoints, int indexval)
 122 // {
 123 //      /*
 124 //      Lj(x) = PI_i!=j (x - xi) / PI_i!=j (xj - xi)
 125 //
 126 //      so Lj'(x) = SUM_k PI_i!=j|k (x - xi) / PI_i!=j (xj - xi)
 127 //      */
 128 //
 129 //      unsigned int i,j,k,i0,i1;
 130 //
 131 //      Real Lpj,mult,div,tj;
 132 //      Real tval = t[indexval];
 133 //
 134 //      //sum k
 135 //      for(j=0;j<npoints;++j)
 136 //      {
 137 //              Lpj = 0;
 138 //              div = 1;
 139 //              tj = t[j];
 140 //
 141 //              for(k=0;k<npoints;++k)
 142 //              {
 143 //                      if(k != j) //because there is no summand for k == j, since that term is missing from the original equation
 144 //                      {
 145 //                              //summation for k
 146 //                              for(i=0;i<npoints;++i)
 147 //                              {
 148 //                                      if(i != k)
 149 //                                      {
 150 //                                              mult *= tval - t[i];
 151 //                                      }
 152 //                              }
 153 //
 154 //                              Lpj += mult; //add into the summation
 155 //
 156 //                              //since the ks follow the exact pattern we need for the divisor (use that too)
 157 //                              div *= tj - t[k];
 158 //                      }
 159 //              }
 160 //
 161 //              //get the actual coefficient
 162 //              Lpj /= div;
 163 //
 164 //              //add it in to the equation
 165 //              df += f[j]*Lpj;
 166 //      }
 167 // }
 168
 169 //END numerical derivatives
 170
 171 // template < class T >
 172 // inline int sign(T f, T tol)
 173 // {
 174 //      if(f < -tol) return -1;
 175 //      if(f > tol) return 1;
 176 //      return 0;
 177 // }
 178
 179 void GetFirstDerivatives(const std::vector<synfig::Point> &f, unsigned int left, unsigned int right, char *out, unsigned int dfstride)
 180 {
 181         unsigned int current = left;
 182
 183         if(right - left < 2)
 184                 return;
 185         else if(right - left == 2)
 186         {
 187                 synfig::Vector v = f[left+1] - f[left];
 188
 189                 //set both to the one we want
 190                 *(synfig::Vector*)out = v;
 191                 out += dfstride;
 192                 *(synfig::Vector*)out = v;
 193                 out += dfstride;
 194         }
 195         else if(right - left < 6/*5*/) //should use 3 point
 196         {
 197                 //left then middle then right
 198                 ThreePointdt(*(synfig::Vector*)out,f[left+0], f[left+1], f[left+2], -1);
 199                 current++;
 200                 out += dfstride;
 201
 202                 for(;current < right-1; current++, out += dfstride)
 203                         ThreePointdt(*(synfig::Vector*)out,f[current-1], f[current], f[current+1], 0);
 204
 205                 ThreePointdt(*(synfig::Vector*)out,f[right-3], f[right-2], f[right-1], 1);
 206                 current++;
 207                 out += dfstride;
 208
 209         }else //can use 5 point
 210         {
 211                 //left 2 then middle bunch then right two
 212                 //may want to use 3 point for inner edge ones
 213
 214                 FivePointdt(*(synfig::Vector*)out,f[left+0], f[left+1], f[left+2], f[left+3], f[left+4], -2);
 215                 out += dfstride;
 216                 FivePointdt(*(synfig::Vector*)out,f[left+1], f[left+2], f[left+3], f[left+4], f[left+5], -1);
 217                 out += dfstride;
 218                 current += 2;
 219
 220                 for(;current < right-2; current++, out += dfstride)
 221                         FivePointdt(*(synfig::Vector*)out,f[current-2], f[current-1], f[current], f[current+1], f[current+2], 0);
 222
 223                 FivePointdt(*(synfig::Vector*)out,f[right-6], f[right-5], f[right-4], f[right-3], f[right-2], 2);
 224                 out += dfstride;
 225                 FivePointdt(*(synfig::Vector*)out,f[right-5], f[right-4], f[right-3], f[right-2], f[right-1], 1);
 226                 out += dfstride;
 227                 current += 2;
 228         }
 229 }
 230
 231 void GetSimpleDerivatives(const std::vector<synfig::Point> &f, int left, int right,
 232                                                         std::vector<synfig::Point> &df, int outleft,
 233                                                         const std::vector<synfig::Real> &/*di*/)
 234 {
 235         int i1,i2,i;
 236         int offset = 2; //df = 1/2 (f[i+o]-f[i-o])
 237
 238         assert((int)df.size() >= right-left+outleft); //must be big enough
 239
 240         for(i = left; i < right; ++i)
 241         {
 242                 //right now indices (figure out distance later)
 243                 i1 = std::max(left,i-offset);
 244                 i2 = std::max(left,i+offset);
 245
 246                 df[outleft++] = (f[i2] - f[i1])*0.5f;
 247         }
 248 }
 249
 250 //get the curve error from the double sample list of work points (hopefully that's enough)
 251 Real CurveError(const synfig::Point *pts, unsigned int n, std::vector<synfig::Point> &work, int left, int right)
 252 {
 253         if(right-left < 2) return -1;
 254
 255         int i,j;
 256
 257         //get distances to each point
 258         Real d,dtemp,dsum;
 259         //synfig::Vector v,vt;
 260         //synfig::Point p1,p2;
 261         synfig::Point pi;
 262         std::vector<synfig::Point>::const_iterator it;//,end = work.begin()+right;
 263
 264         //unsigned int size = work.size();
 265
 266         //for each line, get distance
 267         d = 0; //starts at 0
 268         for(i = 0; i < (int)n; ++i)
 269         {
 270                 pi = pts[i];
 271
 272                 dsum = FLT_MAX;
 273
 274                 it = work.begin()+left;
 275                 //p2 = *it++; //put it at left+1
 276                 for(j = left/*+1*/; j < right; ++j,++it)
 277                 {
 278                         /*p1 = p2;
 279                         p2 = *it;
 280
 281                         v = p2 - p1;
 282                         vt = pi - p1;
 283
 284                         dtemp = v.mag_squared() > 1e-12 ? (vt*v)/v.mag_squared() : 0; //get the projected time value for the current line
 285
 286                         //get distance to line segment with the time value clamped 0-1
 287                         if(dtemp >= 1)  //use p+v
 288                         {
 289                                 vt += v; //makes it pp - (p+v)
 290                         }else if(dtemp > 0)     //use vt-proj
 291                         {
 292                                 vt -= v*dtemp; // vt - proj_v(vt)       //must normalize the projection vector to work
 293                         }
 294
 295                         //else use p
 296                         dtemp = vt.mag_squared();*/
 297
 298                         dtemp = (pi - *it).mag_squared();
 299                         if(dtemp < dsum)
 300                                 dsum = dtemp;
 301                 }
 302
 303                 //accumulate the points' min distance from the curve
 304                 d += sqrt(dsum);
 305         }
 306
 307         return d;
 308 }
 309
 310 typedef synfigapp::BLineConverter::cpindex cpindex;
 311
 312 //has the index data and the tangent scale data (relevant as it may be)
 313 int tessellate_curves(const std::vector<cpindex> &inds, const std::vector<Point> &f, const std::vector<synfig::Vector> &df, std::vector<Point> &work)
 314 {
 315         if(inds.size() < 2)
 316                 return 0;
 317
 318         etl::hermite<Point>     curve;
 319         int ntess = 0;
 320
 321         std::vector<cpindex>::const_iterator j = inds.begin(),j2, end = inds.end();
 322
 323         unsigned int ibase = inds[0].curind;
 324
 325         j2 = j++;
 326         for(; j != end; j2 = j++)
 327         {
 328                 //if this curve has invalid error (in j) then retessellate its work points (requires reparametrization, etc.)
 329                 if(j->error < 0)
 330                 {
 331                         //get the stepsize etc. for the number of points in here
 332                         unsigned int n = j->curind - j2->curind + 1; //thats the number of points in the span
 333                         unsigned int k, kend, i0, i3;
 334                         //so reset the right chunk
 335
 336                         Real t, dt = 1/(Real)(n*2-2); //assuming that they own only n points
 337
 338                         //start at first intermediate
 339                         t = 0;
 340
 341                         i0 = j2->curind; i3 = j->curind;
 342                         k = (i0-ibase)*2; //start on first intermediary point (2x+1)
 343                         kend = (i3-ibase)*2; //last point to set (not intermediary)
 344
 345                         //build hermite curve, it's easier
 346                         curve.p1() = f[i0];
 347                         curve.p2() = f[i3];
 348                         curve.t1() = df[i0-ibase] * (df[i0-ibase].mag_squared() > 1e-4
 349                                                                                  ? j2->tangentscale/df[i0-ibase].mag()
 350                                                                                  : j2->tangentscale);
 351                         curve.t2() = df[i3-ibase] * (df[i3-ibase].mag_squared() > 1e-4
 352                                                                                  ? j->tangentscale/df[i3-ibase].mag()
 353                                                                                  : j->tangentscale);
 354                         curve.sync();
 355
 356                         //MUST include the end point (since we are ignoring left one)
 357                         for(; k < kend; ++k, t += dt)
 358                         {
 359                                 work[k] = curve(t);
 360                         }
 361
 362                         work[k] = curve(1); //k == kend, t == 1 -> c(t) == p2
 363                         ++ntess;
 364                 }
 365         }
 366
 367         return ntess;
 368 }
 369
 370 synfigapp::BLineConverter::BLineConverter()
 371 {
 372         pixelwidth = 1;
 373         smoothness = 0.70f;
 374         width = 0;
 375 };
 376
 377 void
 378 synfigapp::BLineConverter::clear()
 379 {
 380         point_cache.clear();
 381         width_cache.clear();
 382         ftemp.clear();
 383         deriv.clear();
 384         curvature.clear();
 385         break_tangents.clear();
 386         cum_dist.clear();
 387         this_dist.clear();
 388         work.clear();
 389         curind.clear();
 390 }
 391
 392 void
 393 synfigapp::BLineConverter::operator()(std::list<synfig::BLinePoint>  &blinepoints_out,
 394                                                                           const std::list<synfig::Point> &points_in,
 395                                                                           const std::list<synfig::Real>  &widths_in)
 396 {
 397         //Profiling information
 398         /*etl::clock::value_type initialprocess=0, curveval=0, breakeval=0, disteval=0;
 399         etl::clock::value_type preproceval=0, tesseval=0, erroreval=0, spliteval=0;
 400         unsigned int                    numpre=0, numtess=0, numerror=0, numsplit=0;
 401         etl::clock_realtime timer,total;*/
 402
 403         //total.reset();
 404         if (points_in.size() < 2)
 405                 return;
 406
 407         clear();
 408
 409         //removing digitization error harder than expected
 410
 411         //intended to fix little pen errors caused by the way people draw (tiny juts in opposite direction)
 412         //Different solutions
 413         //      Average at both end points (will probably eliminate many points at each end of the samples)
 414         //      Average after the break points are found (weird points would still affect the curve)
 415         //      Just always get rid of breaks at the beginning and end if they are a certain distance apart
 416         //              This is will be current approach so all we do now is try to remove duplicate points
 417
 418         //remove duplicate points - very bad for fitting
 419
 420         //timer.reset();
 421
 422         {
 423                 std::list<synfig::Point>::const_iterator point_iter = points_in.begin(), end = points_in.end();
 424                 std::list<synfig::Real>::const_iterator width_iter = widths_in.begin();
 425                 synfig::Point c;
 426
 427                 if (points_in.size() == widths_in.size())
 428                 {
 429                         for(bool first = true; point_iter != end; ++point_iter,++width_iter)
 430                                 if (first || *point_iter != c)          // eliminate duplicate points
 431                                 {
 432                                         first = false;
 433                                         point_cache.push_back(c = *point_iter);
 434                                         width_cache.push_back(*width_iter);
 435                                 }
 436                 }
 437                 else
 438                         for(;point_iter != end; ++point_iter)
 439                                 if(*point_iter != c)            // eliminate duplicate points
 440                                         point_cache.push_back(c = *point_iter);
 441         }
 442         //initialprocess = timer();
 443
 444         if (point_cache.size() < 7)
 445         {
 446                 info("only %d unique points - giving up", point_cache.size());
 447                 return;
 448         }
 449
 450         //get curvature information
 451         //timer.reset();
 452
 453         {
 454                 int i_this, i_prev, i_next;
 455                 synfig::Vector v_prev, v_next;
 456
 457                 curvature.resize(point_cache.size());
 458                 curvature.front() = curvature.back() = 1;
 459
 460                 for (i_this = 1; i_this < (int)point_cache.size()-1; i_this++)
 461                 {
 462                         i_prev = std::max(0, i_this-2);
 463                         i_next = std::min((int)(point_cache.size()-1), i_this+2);
 464
 465                         v_prev = point_cache[i_this] - point_cache[i_prev];
 466                         v_next = point_cache[i_next] - point_cache[i_this];
 467
 468                         curvature[i_this] = (v_prev*v_next) / (v_prev.mag()*v_next.mag());
 469                 }
 470         }
 471
 472         //curveval = timer();
 473         //synfig::info("calculated curvature");
 474
 475         //find corner points and interpolate inside those
 476         //timer.reset();
 477         {
 478                 //break at sharp derivative points
 479                 //TODO tolerance should be set based upon digitization resolution (length dependent index selection)
 480                 Real    tol = 0;                //break tolerance, for the cosine of the change in angle (really high curvature or something)
 481                 unsigned int i = 0;
 482
 483                 int             sharpest_i=-1;
 484                 int             last=0;
 485                 Real    sharpest_curvature = 1;
 486
 487                 break_tangents.push_back(0);
 488
 489                 // loop through the curvatures; in each continuous run of
 490                 // curvatures that exceed the tolerence, find the one with the
 491                 // sharpest curvature and add its index to the list of indices
 492                 // at which to split tangents
 493                 for (i = 1; i < curvature.size()-1; ++i)
 494                 {
 495                         if (curvature[i] < tol)
 496                         {
 497                                 if(curvature[i] < sharpest_curvature)
 498                                 {
 499                                         sharpest_curvature = curvature[i];
 500                                         sharpest_i = i;
 501                                 }
 502                         }
 503                         else if (sharpest_i > 0)
 504                         {
 505                                 // don't have 2 corners too close to each other
 506                                 if (sharpest_i >= last + 8) //! \todo make this configurable
 507                                 {
 508                                         //synfig::info("break: %d-%d",sharpest_i+1,curvature.size());
 509                                         break_tangents.push_back(sharpest_i);
 510                                         last = sharpest_i;
 511                                 }
 512                                 sharpest_i = -1;
 513                                 sharpest_curvature = 1;
 514                         }
 515                 }
 516
 517                 break_tangents.push_back(i);
 518
 519 // this section causes bug 1892566 if enabled
 520 #if 1
 521                 //postprocess for breaks too close to each other
 522                 Real    fixdistsq = 4*width*width; //the distance to ignore breaks at the end points (for fixing stuff)
 523                 Real d = 0;
 524                 Point p = point_cache[break_tangents.front()];
 525
 526                 //first set
 527                 for (i = 1; i < break_tangents.size()-1; ++i) //do not want to include end point...
 528                 {
 529                         d = (point_cache[break_tangents[i]] - p).mag_squared();
 530                         if(d > fixdistsq) break; //don't want to group breaks if we ever get over the dist...
 531                 }
 532                 //want to erase all points before...
 533                 if(i != 1)
 534                         break_tangents.erase(break_tangents.begin(),break_tangents.begin()+i-1);
 535
 536                 //end set
 537                 p = point_cache[break_tangents.back()];
 538                 for(i = break_tangents.size()-2; i > 0; --i) //start at one in from the end
 539                 {
 540                         d = (point_cache[break_tangents[i]] - p).mag_squared();
 541                         if(d > fixdistsq) break; //don't want to group breaks if we ever get over the dist
 542                 }
 543                 if(i != break_tangents.size()-2)
 544                         break_tangents.erase(break_tangents.begin()+i+2,break_tangents.end()); //erase all points that we found... found none if i has not advanced
 545                 //must not include the one we ended up on
 546 #endif
 547         }
 548         //breakeval = timer();
 549         //synfig::info("found break points: %d",break_tangents.size());
 550
 551         //get the distance calculation of the entire curve (for tangent scaling)
 552
 553         //timer.reset();
 554         {
 555                 synfig::Point p1,p2;
 556
 557                 p1=p2=point_cache[0];
 558
 559                 cum_dist.resize(point_cache.size()); this_dist.resize(point_cache.size());
 560                 Real d = 0;
 561                 for(unsigned int i = 0; i < point_cache.size();)
 562                 {
 563                         d += (this_dist[i] = (p2-p1).mag());
 564                         cum_dist[i] = d;
 565
 566                         p1=p2;
 567                         //! \todo is this legal?  it reads off the end of the vector
 568                         p2=point_cache[++i];
 569                 }
 570         }
 571         //disteval = timer();
 572         //synfig::info("calculated distance");
 573
 574         //now break at every point - calculate new derivatives each time
 575
 576         //TODO
 577         //must be sure that the break points are 3 or more apart
 578         //then must also store the breaks which are not smooth, etc.
 579         //and figure out tangents between there
 580
 581         //for each pair of break points (as long as they are far enough apart) recursively subdivide stuff
 582         //ignore the detected intermediate points
 583         {
 584                 unsigned int i0=0,i3=0,is=0;
 585                 int i=0,j=0;
 586
 587                 bool done = false;
 588
 589                 Real errortol = smoothness*pixelwidth; //???? what the hell should this value be
 590
 591                 BLinePoint a;
 592                 synfig::Vector v;
 593
 594                 //intemp = f; //don't want to smooth out the corners
 595
 596                 bool breaktan = false, setwidth;
 597                 a.set_split_tangent_flag(false);
 598                 //a.set_width(width);
 599                 a.set_width(1.0f);
 600
 601                 setwidth = (point_cache.size() == width_cache.size());
 602
 603                 for(j = 0; j < (int)break_tangents.size() - 1; ++j)
 604                 {
 605                         //for b[j] to b[j+1] subdivide and stuff
 606                         i0 = break_tangents[j];
 607                         i3 = break_tangents[j+1];
 608
 609                         unsigned int size = i3-i0+1; //must include the end points
 610
 611                         //new derivatives
 612                         //timer.reset();
 613                         ftemp.assign(point_cache.begin()+i0, point_cache.begin()+i3+1);
 614                         for(i=0;i<20;++i)
 615                                 gaussian_blur_3(ftemp.begin(),ftemp.end(),false);
 616
 617                         deriv.resize(size);
 618
 619                         // Wondering whether the modification of the deriv vector
 620                         // using a char* pointer and pointer arithmetric was safe,
 621                         // I looked it up...
 622                         //
 623                         // http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2369.pdf tells me:
 624                         //
 625                         //      23.2.5  Class template vector [vector]
 626                         //
 627                         //      [...] The elements of a vector are stored contiguously,
 628                         //      meaning that if v is a vector<T,Allocator> where T is
 629                         //      some type other than bool, then it obeys the identity
 630                         //      &v[n] == &v[0] + n for all 0 <= n < v.size().
 631                         //
 632                         GetFirstDerivatives(ftemp,0,size,(char*)&deriv[0],sizeof(deriv[0]));
 633
 634                         //GetSimpleDerivatives(ftemp,0,size,deriv,0,cum_dist);
 635                         //< don't have to worry about indexing stuff as it is all being taken care of right now
 636                         //preproceval += timer();
 637                         //numpre++;
 638
 639                         work.resize(size*2-1); //guarantee that all points will be tessellated correctly (one point in between every 2 adjacent points)
 640
 641                         //if size of work is size*2-1, the step size should be 1/(size*2 - 2)
 642                         //Real step = 1/(Real)(size*2 - 1);
 643
 644                         //start off with break points as indices
 645                         curind.clear();
 646                         curind.push_back(cpindex(i0,cum_dist[i3]-cum_dist[i0],0)); //0 error because no curve on the left
 647                         curind.push_back(cpindex(i3,cum_dist[i3]-cum_dist[i0],-1)); //error needs to be reevaluated
 648                         done = false; //we want to loop
 649
 650                         unsigned int dcount = 0;
 651
 652                         //while there are still enough points between us, and the error is too high subdivide (and invalidate neighbors that share tangents)
 653                         while(!done)
 654                         {
 655                                 //tessellate all curves with invalid error values
 656                                 work[0] = point_cache[i0];
 657
 658                                 //timer.reset();
 659                                 /*numtess += */tessellate_curves(curind,point_cache,deriv,work);
 660                                 //tesseval += timer();
 661
 662                                 //now get all error values
 663                                 //timer.reset();
 664                                 for(i = 1; i < (int)curind.size(); ++i)
 665                                 {
 666                                         if(curind[i].error < 0) //must have been retessellated, so now recalculate error value
 667                                         {
 668                                                 //evaluate error from points (starting at current index)
 669                                                 int size = curind[i].curind - curind[i-1].curind + 1;
 670                                                 curind[i].error = CurveError(&point_cache[curind[i-1].curind], size,
 671                                                                                                          work,(curind[i-1].curind - i0)*2,(curind[i].curind - i0)*2+1);
 672
 673                                                 /*if(curind[i].error > 1.0e5)
 674                                                 {
 675                                                         synfig::info("Holy crap %d-%d error %f",curind[i-1].curind,curind[i].curind,curind[i].error);
 676                                                         curind[i].error = -1;
 677                                                         numtess += tessellate_curves(curind,f,deriv,work);
 678                                                         curind[i].error = CurveError(&point_cache[curind[i-1].curind], size,
 679                                                                                                          work,0,work.size());//(curind[i-1].curind - i0)*2,(curind[i].curind - i0)*2+1);
 680                                                 }*/
 681                                                 //numerror++;
 682                                         }
 683                                 }
 684                                 //erroreval += timer();
 685
 686                                 //assume we're done
 687                                 done = true;
 688
 689                                 //check each error to see if it's too big, if so, then subdivide etc.
 690                                 int indsize = (int)curind.size();
 691                                 Real maxrelerror = 0;
 692                                 int maxi = -1;//, numpoints;
 693
 694                                 //timer.reset();
 695                                 //get the maximum error and split there
 696                                 for(i = 1; i < indsize; ++i)
 697                                 {
 698                                         //numpoints = curind[i].curind - curind[i-1].curind + 1;
 699
 700                                         if(curind[i].error > maxrelerror && curind[i].curind - curind[i-1].curind > 2) //only accept if it's valid
 701                                         {
 702                                                 maxrelerror = curind[i].error;
 703                                                 maxi = i;
 704                                         }
 705                                 }
 706
 707                                 //split if error is too great
 708                                 if(maxrelerror > errortol)
 709                                 {
 710                                         //add one to the left etc
 711                                         unsigned int    ibase = curind[maxi-1].curind, itop = curind[maxi].curind,
 712                                                                         ibreak = (ibase + itop)/2;
 713                                         Real scale, scale2;
 714
 715                                         assert(ibreak < point_cache.size());
 716
 717                                         //synfig::info("Split %d -%d- %d, error: %f", ibase,ibreak,itop,maxrelerror);
 718
 719                                         if(ibase != itop)
 720                                         {
 721                                                 //invalidate current error of the changed tangents and add an extra segment
 722                                                 //enforce minimum tangents property
 723                                                 curind[maxi].error = -1;
 724                                                 curind[maxi-1].error = -1;
 725                                                 if(maxi+1 < indsize) curind[maxi+1].error = -1; //if there is a curve segment beyond this it will be effected as well
 726
 727                                                 scale = cum_dist[itop] - cum_dist[ibreak];
 728                                                 scale2 = maxi+1 < indsize ? cum_dist[curind[maxi+1].curind] - cum_dist[itop] : scale; //to the right valid?
 729                                                 curind[maxi].tangentscale = std::min(scale, scale2);
 730
 731                                                 scale = cum_dist[ibreak] - cum_dist[ibase];
 732                                                 scale2 = maxi >= 2 ? cum_dist[ibase] - cum_dist[curind[maxi-2].curind] : scale; // to the left valid -2 ?
 733                                                 curind[maxi-1].tangentscale = std::min(scale, scale2);
 734
 735                                                 scale = std::min(cum_dist[ibreak] - cum_dist[ibase], cum_dist[itop] - cum_dist[ibreak]);
 736
 737                                                 curind.insert(curind.begin()+maxi,cpindex(ibreak, scale, -1));
 738                                                 //curind.push_back(cpindex(ibreak, scale, -1));
 739                                                 //std::sort(curind.begin(), curind.end());
 740
 741                                                 done = false;
 742                                                 //numsplit++;
 743                                         }
 744                                 }
 745                                 //spliteval += timer();
 746
 747                                 dcount++;
 748                         }
 749
 750                         //insert the last point too (just set tangent for now
 751                         is = curind[0].curind;
 752
 753                         //first point inherits current tangent status
 754                         v = deriv[is - i0];
 755                         if(v.mag_squared() > EPSILON)
 756                                 v *= (curind[0].tangentscale/v.mag());
 757
 758                         if(!breaktan)
 759                                 a.set_tangent(v);
 760                         else a.set_tangent2(v);
 761
 762                         a.set_vertex(point_cache[is]);
 763                         if(setwidth)a.set_width(width_cache[is]);
 764
 765                         blinepoints_out.push_back(a);
 766                         a.set_split_tangent_flag(false); //won't need to break anymore
 767                         breaktan = false;
 768
 769                         for(i = 1; i < (int)curind.size()-1; ++i)
 770                         {
 771                                 is = curind[i].curind;
 772
 773                                 //first point inherits current tangent status
 774                                 v = deriv[is-i0];
 775                                 if(v.mag_squared() > EPSILON)
 776                                         v *= (curind[i].tangentscale/v.mag());
 777
 778                                 a.set_tangent(v); // always inside, so guaranteed to be smooth
 779                                 a.set_vertex(point_cache[is]);
 780                                 if(setwidth)a.set_width(width_cache[is]);
 781
 782                                 blinepoints_out.push_back(a);
 783                         }
 784
 785                         //set the last point's data
 786                         is = curind.back().curind; //should already be this
 787
 788                         v = deriv[is-i0];
 789                         if(v.mag_squared() > EPSILON)
 790                                 v *= (curind.back().tangentscale/v.mag());
 791
 792                         a.set_tangent1(v);
 793                         a.set_split_tangent_flag(true);
 794                         breaktan = true;
 795
 796                         //will get the vertex and tangent 2 from next round
 797                 }
 798
 799                 a.set_vertex(point_cache[i3]);
 800                 a.set_split_tangent_flag(false);
 801                 if(setwidth)
 802                         a.set_width(width_cache[i3]);
 803                 blinepoints_out.push_back(a);
 804
 805                 /*etl::clock::value_type totaltime = total(),
 806                                                            misctime = totaltime - initialprocess - curveval - breakeval - disteval
 807                                                                           - preproceval - tesseval - erroreval - spliteval;
 808
 809                 synfig::info(
 810                         "Curve Convert Profile:\n"
 811                         "\tInitial Preprocess:    %f\n"
 812                         "\tCurvature Calculation: %f\n"
 813                         "\tBreak Calculation:     %f\n"
 814                         "\tDistance Calculation:  %f\n"
 815                         "  Algorithm: (numtimes,totaltime)\n"
 816                         "\tPreprocess step:      (%d,%f)\n"
 817                         "\tTessellation step:    (%d,%f)\n"
 818                         "\tError step:           (%d,%f)\n"
 819                         "\tSplit step:           (%d,%f)\n"
 820                         "  Num Input: %d, Num Output: %d\n"
 821                         "  Total time: %f, Misc time: %f\n",
 822                         initialprocess, curveval,breakeval,disteval,
 823                         numpre,preproceval,numtess,tesseval,numerror,erroreval,numsplit,spliteval,
 824                         points_in.size(),blinepoints_out.size(),
 825                         totaltime,misctime);*/
 826
 827                 return;
 828         }
 829 }
 830
 831 void synfigapp::BLineConverter::EnforceMinWidth(std::list<synfig::BLinePoint> &bline, synfig::Real min_pressure)
 832 {
 833         std::list<synfig::BLinePoint>::iterator i = bline.begin(),
 834                                                                                         end = bline.end();
 835
 836         for(i = bline.begin(); i != end; ++i)
 837                 if(i->get_width() < min_pressure)
 838                         i->set_width(min_pressure);
 839 }