1 /* === S Y N F I G ========================================================= */
2 /*! \file blineconvert.cpp
3 ** \brief Template File
8 ** Copyright (c) 2002-2005 Robert B. Quattlebaum Jr., Adrian Bentley
10 ** This package is free software; you can redistribute it and/or
11 ** modify it under the terms of the GNU General Public License as
12 ** published by the Free Software Foundation; either version 2 of
13 ** the License, or (at your option) any later version.
15 ** This package is distributed in the hope that it will be useful,
16 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
17 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 ** General Public License for more details.
21 /* ========================================================================= */
23 /* === H E A D E R S ======================================================= */
32 #include "blineconvert.h"
34 #include <ETL/gaussian>
35 #include <ETL/hermite>
39 #include <synfig/general.h>
46 /* === U S I N G =========================================================== */
50 using namespace synfig;
52 /* === M A C R O S ========================================================= */
54 #define EPSILON (1e-10)
56 /* === G L O B A L S ======================================================= */
58 /* === P R O C E D U R E S ================================================= */
60 /* === M E T H O D S ======================================================= */
63 //Derivative Functions for numerical approximation
65 //bias == 0 will get F' at f3, bias < 0 will get F' at f1, and bias > 0 will get F' at f5
67 inline void FivePointdt(T &df, const T &f1, const T &f2, const T &f3, const T &f4, const T &f5, int bias)
69 if (bias == 0) // middle
70 df = (f1 - f2*8 + f4*8 - f5)*(1/12.0f);
71 else if (bias < 0) // left
72 df = (-f1*25 + f2*48 - f3*36 + f4*16 - f5*3)*(1/12.0f);
74 df = (f1*3 - f2*16 + f3*36 - f4*48 + f5*25)*(1/12.0f);
78 inline void ThreePointdt(T &df, const T &f1, const T &f2, const T &f3, int bias)
80 if (bias == 0) // middle
81 df = (-f1 + f3)*(1/2.0f);
82 else if (bias < 0) // left
83 df = (-f1*3 + f2*4 - f3)*(1/2.0f);
85 df = (f1 - f2*4 + f3*3)*(1/2.0f);
88 // template < class T >
89 // inline void ThreePointddt(T &df, const T &f1, const T &f2, const T &f3, int bias)
91 // // a 3 point approximation pretends to have constant acceleration,
92 // // so only one algorithm needed for left, middle, or right
93 // df = (f1 -f2*2 + f3)*(1/2.0f);
96 // // WARNING -- totally broken
97 // template < class T >
98 // inline void FivePointddt(T &df, const T &f1, const T &f2, const T &f3, int bias)
104 // //df = (- f1 + f2*16 - f3*30 + f4*16 - f5)*(1/12.0f);
105 // }/*else if(bias < 0)
108 // df = (f1*7 - f2*26*4 + f3*19*6 - f4*14*4 + f5*11)*(1/12.0f);
112 // df = (f1*3 - f2*16 + f3*36 - f4*48 + f5*25)*(1/12.0f);
114 // //side ones don't work, use 3 point
117 // //implement an arbitrary derivative
119 // template < class T >
120 // void DerivativeApprox(T &df, const T f[], const Real t[], int npoints, int indexval)
123 // Lj(x) = PI_i!=j (x - xi) / PI_i!=j (xj - xi)
125 // so Lj'(x) = SUM_k PI_i!=j|k (x - xi) / PI_i!=j (xj - xi)
128 // unsigned int i,j,k,i0,i1;
130 // Real Lpj,mult,div,tj;
131 // Real tval = t[indexval];
134 // for(j=0;j<npoints;++j)
140 // for(k=0;k<npoints;++k)
142 // if(k != j) //because there is no summand for k == j, since that term is missing from the original equation
145 // for(i=0;i<npoints;++i)
149 // mult *= tval - t[i];
153 // Lpj += mult; //add into the summation
155 // //since the ks follow the exact pattern we need for the divisor (use that too)
160 // //get the actual coefficient
163 // //add it in to the equation
168 //END numerical derivatives
170 // template < class T >
171 // inline int sign(T f, T tol)
173 // if(f < -tol) return -1;
174 // if(f > tol) return 1;
178 void GetFirstDerivatives(const std::vector<synfig::Point> &f, unsigned int left, unsigned int right, char *out, unsigned int dfstride)
180 unsigned int current = left;
184 else if(right - left == 2)
186 synfig::Vector v = f[left+1] - f[left];
188 //set both to the one we want
189 *(synfig::Vector*)out = v;
191 *(synfig::Vector*)out = v;
194 else if(right - left < 6/*5*/) //should use 3 point
196 //left then middle then right
197 ThreePointdt(*(synfig::Vector*)out,f[left+0], f[left+1], f[left+2], -1);
201 for(;current < right-1; current++, out += dfstride)
202 ThreePointdt(*(synfig::Vector*)out,f[current-1], f[current], f[current+1], 0);
204 ThreePointdt(*(synfig::Vector*)out,f[right-3], f[right-2], f[right-1], 1);
208 }else //can use 5 point
210 //left 2 then middle bunch then right two
211 //may want to use 3 point for inner edge ones
213 FivePointdt(*(synfig::Vector*)out,f[left+0], f[left+1], f[left+2], f[left+3], f[left+4], -2);
215 FivePointdt(*(synfig::Vector*)out,f[left+1], f[left+2], f[left+3], f[left+4], f[left+5], -1);
219 for(;current < right-2; current++, out += dfstride)
220 FivePointdt(*(synfig::Vector*)out,f[current-2], f[current-1], f[current], f[current+1], f[current+2], 0);
222 FivePointdt(*(synfig::Vector*)out,f[right-6], f[right-5], f[right-4], f[right-3], f[right-2], 2);
224 FivePointdt(*(synfig::Vector*)out,f[right-5], f[right-4], f[right-3], f[right-2], f[right-1], 1);
230 void GetSimpleDerivatives(const std::vector<synfig::Point> &f, int left, int right,
231 std::vector<synfig::Point> &df, int outleft,
232 const std::vector<synfig::Real> &/*di*/)
235 int offset = 2; //df = 1/2 (f[i+o]-f[i-o])
237 assert((int)df.size() >= right-left+outleft); //must be big enough
239 for(i = left; i < right; ++i)
241 //right now indices (figure out distance later)
242 i1 = std::max(left,i-offset);
243 i2 = std::max(left,i+offset);
245 df[outleft++] = (f[i2] - f[i1])*0.5f;
249 //get the curve error from the double sample list of work points (hopefully that's enough)
250 Real CurveError(const synfig::Point *pts, unsigned int n, std::vector<synfig::Point> &work, int left, int right)
252 if(right-left < 2) return -1;
256 //get distances to each point
258 //synfig::Vector v,vt;
259 //synfig::Point p1,p2;
261 std::vector<synfig::Point>::const_iterator it;//,end = work.begin()+right;
263 //unsigned int size = work.size();
265 //for each line, get distance
267 for(i = 0; i < (int)n; ++i)
273 it = work.begin()+left;
274 //p2 = *it++; //put it at left+1
275 for(j = left/*+1*/; j < right; ++j,++it)
283 dtemp = v.mag_squared() > 1e-12 ? (vt*v)/v.mag_squared() : 0; //get the projected time value for the current line
285 //get distance to line segment with the time value clamped 0-1
286 if(dtemp >= 1) //use p+v
288 vt += v; //makes it pp - (p+v)
289 }else if(dtemp > 0) //use vt-proj
291 vt -= v*dtemp; // vt - proj_v(vt) //must normalize the projection vector to work
295 dtemp = vt.mag_squared();*/
297 dtemp = (pi - *it).mag_squared();
302 //accumulate the points' min distance from the curve
309 typedef synfigapp::BLineConverter::cpindex cpindex;
311 //has the index data and the tangent scale data (relevant as it may be)
312 int tessellate_curves(const std::vector<cpindex> &inds, const std::vector<Point> &f, const std::vector<synfig::Vector> &df, std::vector<Point> &work)
317 etl::hermite<Point> curve;
320 std::vector<cpindex>::const_iterator j = inds.begin(),j2, end = inds.end();
322 unsigned int ibase = inds[0].curind;
325 for(; j != end; j2 = j++)
327 //if this curve has invalid error (in j) then retessellate its work points (requires reparametrization, etc.)
330 //get the stepsize etc. for the number of points in here
331 unsigned int n = j->curind - j2->curind + 1; //thats the number of points in the span
332 unsigned int k, kend, i0, i3;
333 //so reset the right chunk
335 Real t, dt = 1/(Real)(n*2-2); //assuming that they own only n points
337 //start at first intermediate
340 i0 = j2->curind; i3 = j->curind;
341 k = (i0-ibase)*2; //start on first intermediary point (2x+1)
342 kend = (i3-ibase)*2; //last point to set (not intermediary)
344 //build hermite curve, it's easier
347 curve.t1() = df[i0-ibase] * (df[i0-ibase].mag_squared() > 1e-4
348 ? j2->tangentscale/df[i0-ibase].mag()
350 curve.t2() = df[i3-ibase] * (df[i3-ibase].mag_squared() > 1e-4
351 ? j->tangentscale/df[i3-ibase].mag()
355 //MUST include the end point (since we are ignoring left one)
356 for(; k < kend; ++k, t += dt)
361 work[k] = curve(1); //k == kend, t == 1 -> c(t) == p2
369 synfigapp::BLineConverter::BLineConverter()
377 synfigapp::BLineConverter::clear()
392 synfigapp::BLineConverter::operator () (std::list<synfig::BLinePoint> &out, const std::list<synfig::Point> &in,const std::list<synfig::Real> &in_w)
394 //Profiling information
395 /*etl::clock::value_type initialprocess=0, curveval=0, breakeval=0, disteval=0;
396 etl::clock::value_type preproceval=0, tesseval=0, erroreval=0, spliteval=0;
397 unsigned int numpre=0, numtess=0, numerror=0, numsplit=0;
398 etl::clock_realtime timer,total;*/
406 //removing digitization error harder than expected
408 //intended to fix little pen errors caused by the way people draw (tiny juts in opposite direction)
409 //Different solutions
410 // Average at both end points (will probably eliminate many points at each end of the samples)
411 // Average after the break points are found (weird points would still affect the curve)
412 // Just always get rid of breaks at the beginning and end if they are a certain distance apart
413 // This is will be current approach so all we do now is try to remove duplicate points
415 //remove duplicate points - very bad for fitting
420 std::list<synfig::Point>::const_iterator i = in.begin(), end = in.end();
421 std::list<synfig::Real>::const_iterator iw = in_w.begin();
424 if(in.size() == in_w.size())
426 for(;i != end; ++i,++iw)
427 if(*i != c) // eliminate duplicate points
436 if(*i != c) // eliminate duplicate points
440 //initialprocess = timer();
445 //get curvature information
450 synfig::Vector v1,v2;
452 cvt.resize(f.size());
457 for(i = 1; i < (int)f.size()-1; ++i)
459 i0 = std::max(0,i - 2);
460 i1 = std::min((int)(f.size()-1),i + 2);
465 cvt[i] = (v1*v2)/(v1.mag()*v2.mag());
469 //curveval = timer();
470 //synfig::info("calculated curvature");
472 //find corner points and interpolate inside those
475 //break at sharp derivative points
476 //TODO tolerance should be set based upon digitization resolution (length dependent index selection)
477 Real tol = 0; //break tolerance, for the cosine of the change in angle (really high curvature or something)
478 Real fixdistsq = 4*width*width; //the distance to ignore breaks at the end points (for fixing stuff)
481 int maxi = -1, last=0;
486 for(i = 1; i < cvt.size()-1; ++i)
488 //insert if too sharp (we need to break the tangents to insert onto the break list)
502 //synfig::info("break: %d-%d",maxi+1,cvt.size());
513 //postprocess for breaks too close to each other
515 Point p = f[brk.front()];
518 for(i = 1; i < brk.size()-1; ++i) //do not want to include end point...
520 d = (f[brk[i]] - p).mag_squared();
521 if(d > fixdistsq) break; //don't want to group breaks if we ever get over the dist...
523 //want to erase all points before...
525 brk.erase(brk.begin(),brk.begin()+i-1);
529 for(i = brk.size()-2; i > 0; --i) //start at one in from the end
531 d = (f[brk[i]] - p).mag_squared();
532 if(d > fixdistsq) break; //don't want to group breaks if we ever get over the dist
534 if(i != brk.size()-2)
535 brk.erase(brk.begin()+i+2,brk.end()); //erase all points that we found... found none if i has not advanced
536 //must not include the one we ended up on
538 //breakeval = timer();
539 //synfig::info("found break points: %d",brk.size());
541 //get the distance calculation of the entire curve (for tangent scaling)
549 di.resize(f.size()); d_i.resize(f.size());
551 for(unsigned int i = 0; i < f.size();)
553 d += (d_i[i] = (p2-p1).mag());
560 //disteval = timer();
561 //synfig::info("calculated distance");
563 //now break at every point - calculate new derivatives each time
566 //must be sure that the break points are 3 or more apart
567 //then must also store the breaks which are not smooth, etc.
568 //and figure out tangents between there
570 //for each pair of break points (as long as they are far enough apart) recursively subdivide stuff
571 //ignore the detected intermediate points
573 unsigned int i0=0,i3=0,is=0;
578 Real errortol = smoothness*pixelwidth; //???? what the hell should this value be
583 //intemp = f; //don't want to smooth out the corners
585 bool breaktan = false, setwidth;
586 a.set_split_tangent_flag(false);
587 //a.set_width(width);
590 setwidth = (f.size() == f_w.size());
592 for(j = 0; j < (int)brk.size() - 1; ++j)
594 //for b[j] to b[j+1] subdivide and stuff
598 unsigned int size = i3-i0+1; //must include the end points
602 ftemp.assign(f.begin()+i0, f.begin()+i3+1);
604 gaussian_blur_3(ftemp.begin(),ftemp.end(),false);
608 // Wondering whether the modification of the df vector
609 // using a char* pointer and pointer arithmetric was safe,
612 // http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2369.pdf tells me:
614 // 23.2.5 Class template vector [vector]
616 // [...] The elements of a vector are stored contiguously,
617 // meaning that if v is a vector<T,Allocator> where T is
618 // some type other than bool, then it obeys the identity
619 // &v[n] == &v[0] + n for all 0 <= n < v.size().
621 GetFirstDerivatives(ftemp,0,size,(char*)&df[0],sizeof(df[0]));
623 //GetSimpleDerivatives(ftemp,0,size,df,0,di);
624 //< don't have to worry about indexing stuff as it is all being taken care of right now
625 //preproceval += timer();
628 work.resize(size*2-1); //guarantee that all points will be tessellated correctly (one point in between every 2 adjacent points)
630 //if size of work is size*2-1, the step size should be 1/(size*2 - 2)
631 //Real step = 1/(Real)(size*2 - 1);
633 //start off with break points as indices
635 curind.push_back(cpindex(i0,di[i3]-di[i0],0)); //0 error because no curve on the left
636 curind.push_back(cpindex(i3,di[i3]-di[i0],-1)); //error needs to be reevaluated
637 done = false; //we want to loop
639 unsigned int dcount = 0;
641 //while there are still enough points between us, and the error is too high subdivide (and invalidate neighbors that share tangents)
644 //tessellate all curves with invalid error values
648 /*numtess += */tessellate_curves(curind,f,df,work);
649 //tesseval += timer();
651 //now get all error values
653 for(i = 1; i < (int)curind.size(); ++i)
655 if(curind[i].error < 0) //must have been retessellated, so now recalculate error value
657 //evaluate error from points (starting at current index)
658 int size = curind[i].curind - curind[i-1].curind + 1;
659 curind[i].error = CurveError(&f[curind[i-1].curind], size,
660 work,(curind[i-1].curind - i0)*2,(curind[i].curind - i0)*2+1);
662 /*if(curind[i].error > 1.0e5)
664 synfig::info("Holy crap %d-%d error %f",curind[i-1].curind,curind[i].curind,curind[i].error);
665 curind[i].error = -1;
666 numtess += tessellate_curves(curind,f,df,work);
667 curind[i].error = CurveError(&f[curind[i-1].curind], size,
668 work,0,work.size());//(curind[i-1].curind - i0)*2,(curind[i].curind - i0)*2+1);
673 //erroreval += timer();
678 //check each error to see if it's too big, if so, then subdivide etc.
679 int indsize = (int)curind.size();
680 Real maxrelerror = 0;
681 int maxi = -1;//, numpoints;
684 //get the maximum error and split there
685 for(i = 1; i < indsize; ++i)
687 //numpoints = curind[i].curind - curind[i-1].curind + 1;
689 if(curind[i].error > maxrelerror && curind[i].curind - curind[i-1].curind > 2) //only accept if it's valid
691 maxrelerror = curind[i].error;
696 //split if error is too great
697 if(maxrelerror > errortol)
699 //add one to the left etc
700 unsigned int ibase = curind[maxi-1].curind, itop = curind[maxi].curind,
701 ibreak = (ibase + itop)/2;
704 assert(ibreak < f.size());
706 //synfig::info("Split %d -%d- %d, error: %f", ibase,ibreak,itop,maxrelerror);
710 //invalidate current error of the changed tangents and add an extra segment
711 //enforce minimum tangents property
712 curind[maxi].error = -1;
713 curind[maxi-1].error = -1;
714 if(maxi+1 < indsize) curind[maxi+1].error = -1; //if there is a curve segment beyond this it will be effected as well
716 scale = di[itop] - di[ibreak];
717 scale2 = maxi+1 < indsize ? di[curind[maxi+1].curind] - di[itop] : scale; //to the right valid?
718 curind[maxi].tangentscale = std::min(scale, scale2);
720 scale = di[ibreak] - di[ibase];
721 scale2 = maxi >= 2 ? di[ibase] - di[curind[maxi-2].curind] : scale; // to the left valid -2 ?
722 curind[maxi-1].tangentscale = std::min(scale, scale2);
724 scale = std::min(di[ibreak] - di[ibase], di[itop] - di[ibreak]);
726 curind.insert(curind.begin()+maxi,cpindex(ibreak, scale, -1));
727 //curind.push_back(cpindex(ibreak, scale, -1));
728 //std::sort(curind.begin(), curind.end());
734 //spliteval += timer();
739 //insert the last point too (just set tangent for now
740 is = curind[0].curind;
742 //first point inherits current tangent status
744 if(v.mag_squared() > EPSILON)
745 v *= (curind[0].tangentscale/v.mag());
749 else a.set_tangent2(v);
752 if(setwidth)a.set_width(f_w[is]);
755 a.set_split_tangent_flag(false); //won't need to break anymore
758 for(i = 1; i < (int)curind.size()-1; ++i)
760 is = curind[i].curind;
762 //first point inherits current tangent status
764 if(v.mag_squared() > EPSILON)
765 v *= (curind[i].tangentscale/v.mag());
767 a.set_tangent(v); // always inside, so guaranteed to be smooth
769 if(setwidth)a.set_width(f_w[is]);
774 //set the last point's data
775 is = curind.back().curind; //should already be this
778 if(v.mag_squared() > EPSILON)
779 v *= (curind.back().tangentscale/v.mag());
782 a.set_split_tangent_flag(true);
785 //will get the vertex and tangent 2 from next round
789 a.set_split_tangent_flag(false);
791 a.set_width(f_w[i3]);
794 /*etl::clock::value_type totaltime = total(),
795 misctime = totaltime - initialprocess - curveval - breakeval - disteval
796 - preproceval - tesseval - erroreval - spliteval;
799 "Curve Convert Profile:\n"
800 "\tInitial Preprocess: %f\n"
801 "\tCurvature Calculation: %f\n"
802 "\tBreak Calculation: %f\n"
803 "\tDistance Calculation: %f\n"
804 " Algorithm: (numtimes,totaltime)\n"
805 "\tPreprocess step: (%d,%f)\n"
806 "\tTessellation step: (%d,%f)\n"
807 "\tError step: (%d,%f)\n"
808 "\tSplit step: (%d,%f)\n"
809 " Num Input: %d, Num Output: %d\n"
810 " Total time: %f, Misc time: %f\n",
811 initialprocess, curveval,breakeval,disteval,
812 numpre,preproceval,numtess,tesseval,numerror,erroreval,numsplit,spliteval,
813 in.size(),out.size(),
814 totaltime,misctime);*/
820 void synfigapp::BLineConverter::EnforceMinWidth(std::list<synfig::BLinePoint> &bline, synfig::Real min_pressure)
822 std::list<synfig::BLinePoint>::iterator i = bline.begin(),
825 for(i = bline.begin(); i != end; ++i)
826 if(i->get_width() < min_pressure)
827 i->set_width(min_pressure);