Updated towards 0.61.09 release
[synfig.git] / synfig-studio / trunk / 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 }