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>
44 /* === U S I N G =========================================================== */
48 using namespace synfig;
50 /* === M A C R O S ========================================================= */
52 #define EPSILON (1e-10)
54 /* === G L O B A L S ======================================================= */
56 /* === P R O C E D U R E S ================================================= */
58 /* === M E T H O D S ======================================================= */
61 //Derivative Functions for numerical approximation
63 //bias == 0 will get F' at f3, bias < 0 will get F' at f1, and bias > 0 will get F' at f5
65 inline void FivePointdt(T &df, const T &f1, const T &f2, const T &f3, const T &f4, const T &f5, int bias)
67 if (bias == 0) // middle
68 df = (f1 - f2*8 + f4*8 - f5)*(1/12.0f);
69 else if (bias < 0) // left
70 df = (-f1*25 + f2*48 - f3*36 + f4*16 - f5*3)*(1/12.0f);
72 df = (f1*3 - f2*16 + f3*36 - f4*48 + f5*25)*(1/12.0f);
76 inline void ThreePointdt(T &df, const T &f1, const T &f2, const T &f3, int bias)
78 if (bias == 0) // middle
79 df = (-f1 + f3)*(1/2.0f);
80 else if (bias < 0) // left
81 df = (-f1*3 + f2*4 - f3)*(1/2.0f);
83 df = (f1 - f2*4 + f3*3)*(1/2.0f);
86 // template < class T >
87 // inline void ThreePointddt(T &df, const T &f1, const T &f2, const T &f3, int bias)
89 // // a 3 point approximation pretends to have constant acceleration,
90 // // so only one algorithm needed for left, middle, or right
91 // df = (f1 -f2*2 + f3)*(1/2.0f);
94 // // WARNING -- totally broken
95 // template < class T >
96 // inline void FivePointddt(T &df, const T &f1, const T &f2, const T &f3, int bias)
102 // //df = (- f1 + f2*16 - f3*30 + f4*16 - f5)*(1/12.0f);
103 // }/*else if(bias < 0)
106 // df = (f1*7 - f2*26*4 + f3*19*6 - f4*14*4 + f5*11)*(1/12.0f);
110 // df = (f1*3 - f2*16 + f3*36 - f4*48 + f5*25)*(1/12.0f);
112 // //side ones don't work, use 3 point
115 // //implement an arbitrary derivative
117 // template < class T >
118 // void DerivativeApprox(T &df, const T f[], const Real t[], int npoints, int indexval)
121 // Lj(x) = PI_i!=j (x - xi) / PI_i!=j (xj - xi)
123 // so Lj'(x) = SUM_k PI_i!=j|k (x - xi) / PI_i!=j (xj - xi)
126 // unsigned int i,j,k,i0,i1;
128 // Real Lpj,mult,div,tj;
129 // Real tval = t[indexval];
132 // for(j=0;j<npoints;++j)
138 // for(k=0;k<npoints;++k)
140 // if(k != j) //because there is no summand for k == j, since that term is missing from the original equation
143 // for(i=0;i<npoints;++i)
147 // mult *= tval - t[i];
151 // Lpj += mult; //add into the summation
153 // //since the ks follow the exact pattern we need for the divisor (use that too)
158 // //get the actual coefficient
161 // //add it in to the equation
166 //END numerical derivatives
168 // template < class T >
169 // inline int sign(T f, T tol)
171 // if(f < -tol) return -1;
172 // if(f > tol) return 1;
176 void GetFirstDerivatives(const std::vector<synfig::Point> &f, unsigned int left, unsigned int right, char *out, unsigned int dfstride)
178 unsigned int current = left;
182 else if(right - left == 2)
184 synfig::Vector v = f[left+1] - f[left];
186 //set both to the one we want
187 *(synfig::Vector*)out = v;
189 *(synfig::Vector*)out = v;
192 else if(right - left < 6/*5*/) //should use 3 point
194 //left then middle then right
195 ThreePointdt(*(synfig::Vector*)out,f[left+0], f[left+1], f[left+2], -1);
199 for(;current < right-1; current++, out += dfstride)
200 ThreePointdt(*(synfig::Vector*)out,f[current-1], f[current], f[current+1], 0);
202 ThreePointdt(*(synfig::Vector*)out,f[right-3], f[right-2], f[right-1], 1);
206 }else //can use 5 point
208 //left 2 then middle bunch then right two
209 //may want to use 3 point for inner edge ones
211 FivePointdt(*(synfig::Vector*)out,f[left+0], f[left+1], f[left+2], f[left+3], f[left+4], -2);
213 FivePointdt(*(synfig::Vector*)out,f[left+1], f[left+2], f[left+3], f[left+4], f[left+5], -1);
217 for(;current < right-2; current++, out += dfstride)
218 FivePointdt(*(synfig::Vector*)out,f[current-2], f[current-1], f[current], f[current+1], f[current+2], 0);
220 FivePointdt(*(synfig::Vector*)out,f[right-6], f[right-5], f[right-4], f[right-3], f[right-2], 2);
222 FivePointdt(*(synfig::Vector*)out,f[right-5], f[right-4], f[right-3], f[right-2], f[right-1], 1);
228 void GetSimpleDerivatives(const std::vector<synfig::Point> &f, int left, int right,
229 std::vector<synfig::Point> &df, int outleft,
230 const std::vector<synfig::Real> &/*di*/)
233 int offset = 2; //df = 1/2 (f[i+o]-f[i-o])
235 assert((int)df.size() >= right-left+outleft); //must be big enough
237 for(i = left; i < right; ++i)
239 //right now indices (figure out distance later)
240 i1 = std::max(left,i-offset);
241 i2 = std::max(left,i+offset);
243 df[outleft++] = (f[i2] - f[i1])*0.5f;
247 //get the curve error from the double sample list of work points (hopefully that's enough)
248 Real CurveError(const synfig::Point *pts, unsigned int n, std::vector<synfig::Point> &work, int left, int right)
250 if(right-left < 2) return -1;
254 //get distances to each point
256 //synfig::Vector v,vt;
257 //synfig::Point p1,p2;
259 std::vector<synfig::Point>::const_iterator it;//,end = work.begin()+right;
261 //unsigned int size = work.size();
263 //for each line, get distance
265 for(i = 0; i < (int)n; ++i)
271 it = work.begin()+left;
272 //p2 = *it++; //put it at left+1
273 for(j = left/*+1*/; j < right; ++j,++it)
281 dtemp = v.mag_squared() > 1e-12 ? (vt*v)/v.mag_squared() : 0; //get the projected time value for the current line
283 //get distance to line segment with the time value clamped 0-1
284 if(dtemp >= 1) //use p+v
286 vt += v; //makes it pp - (p+v)
287 }else if(dtemp > 0) //use vt-proj
289 vt -= v*dtemp; // vt - proj_v(vt) //must normalize the projection vector to work
293 dtemp = vt.mag_squared();*/
295 dtemp = (pi - *it).mag_squared();
300 //accumulate the points' min distance from the curve
307 typedef synfigapp::BLineConverter::cpindex cpindex;
309 //has the index data and the tangent scale data (relevant as it may be)
310 int tessellate_curves(const std::vector<cpindex> &inds, const std::vector<Point> &f, const std::vector<synfig::Vector> &df, std::vector<Point> &work)
315 etl::hermite<Point> curve;
318 std::vector<cpindex>::const_iterator j = inds.begin(),j2, end = inds.end();
320 unsigned int ibase = inds[0].curind;
323 for(; j != end; j2 = j++)
325 //if this curve has invalid error (in j) then retessellate its work points (requires reparametrization, etc.)
328 //get the stepsize etc. for the number of points in here
329 unsigned int n = j->curind - j2->curind + 1; //thats the number of points in the span
330 unsigned int k, kend, i0, i3;
331 //so reset the right chunk
333 Real t, dt = 1/(Real)(n*2-2); //assuming that they own only n points
335 //start at first intermediate
338 i0 = j2->curind; i3 = j->curind;
339 k = (i0-ibase)*2; //start on first intermediary point (2x+1)
340 kend = (i3-ibase)*2; //last point to set (not intermediary)
342 //build hermite curve, it's easier
345 curve.t1() = df[i0-ibase] * (df[i0-ibase].mag_squared() > 1e-4
346 ? j2->tangentscale/df[i0-ibase].mag()
348 curve.t2() = df[i3-ibase] * (df[i3-ibase].mag_squared() > 1e-4
349 ? j->tangentscale/df[i3-ibase].mag()
353 //MUST include the end point (since we are ignoring left one)
354 for(; k < kend; ++k, t += dt)
359 work[k] = curve(1); //k == kend, t == 1 -> c(t) == p2
367 synfigapp::BLineConverter::BLineConverter()
375 synfigapp::BLineConverter::clear()
390 synfigapp::BLineConverter::operator () (std::list<synfig::BLinePoint> &out, const std::list<synfig::Point> &in,const std::list<synfig::Real> &in_w)
392 //Profiling information
393 /*etl::clock::value_type initialprocess=0, curveval=0, breakeval=0, disteval=0;
394 etl::clock::value_type preproceval=0, tesseval=0, erroreval=0, spliteval=0;
395 unsigned int numpre=0, numtess=0, numerror=0, numsplit=0;
396 etl::clock_realtime timer,total;*/
404 //removing digitization error harder than expected
406 //intended to fix little pen errors caused by the way people draw (tiny juts in opposite direction)
407 //Different solutions
408 // Average at both end points (will probably eliminate many points at each end of the samples)
409 // Average after the break points are found (weird points would still affect the curve)
410 // Just always get rid of breaks at the beginning and end if they are a certain distance apart
411 // This is will be current approach so all we do now is try to remove duplicate points
413 //remove duplicate points - very bad for fitting
418 std::list<synfig::Point>::const_iterator i = in.begin(), end = in.end();
419 std::list<synfig::Real>::const_iterator iw = in_w.begin();
422 if(in.size() == in_w.size())
424 for(;i != end; ++i,++iw)
425 if(*i != c) // eliminate duplicate points
434 if(*i != c) // eliminate duplicate points
438 //initialprocess = timer();
443 //get curvature information
448 synfig::Vector v1,v2;
450 cvt.resize(f.size());
455 for(i = 1; i < (int)f.size()-1; ++i)
457 i0 = std::max(0,i - 2);
458 i1 = std::min((int)(f.size()-1),i + 2);
463 cvt[i] = (v1*v2)/(v1.mag()*v2.mag());
467 //curveval = timer();
468 //synfig::info("calculated curvature");
470 //find corner points and interpolate inside those
473 //break at sharp derivative points
474 //TODO tolerance should be set based upon digitization resolution (length dependent index selection)
475 Real tol = 0; //break tolerance, for the cosine of the change in angle (really high curvature or something)
476 Real fixdistsq = 4*width*width; //the distance to ignore breaks at the end points (for fixing stuff)
479 int maxi = -1, last=0;
484 for(i = 1; i < cvt.size()-1; ++i)
486 //insert if too sharp (we need to break the tangents to insert onto the break list)
500 //synfig::info("break: %d-%d",maxi+1,cvt.size());
511 //postprocess for breaks too close to each other
513 Point p = f[brk.front()];
516 for(i = 1; i < brk.size()-1; ++i) //do not want to include end point...
518 d = (f[brk[i]] - p).mag_squared();
519 if(d > fixdistsq) break; //don't want to group breaks if we ever get over the dist...
521 //want to erase all points before...
523 brk.erase(brk.begin(),brk.begin()+i-1);
527 for(i = brk.size()-2; i > 0; --i) //start at one in from the end
529 d = (f[brk[i]] - p).mag_squared();
530 if(d > fixdistsq) break; //don't want to group breaks if we ever get over the dist
532 if(i != brk.size()-2)
533 brk.erase(brk.begin()+i+2,brk.end()); //erase all points that we found... found none if i has not advanced
534 //must not include the one we ended up on
536 //breakeval = timer();
537 //synfig::info("found break points: %d",brk.size());
539 //get the distance calculation of the entire curve (for tangent scaling)
547 di.resize(f.size()); d_i.resize(f.size());
549 for(unsigned int i = 0; i < f.size();)
551 d += (d_i[i] = (p2-p1).mag());
558 //disteval = timer();
559 //synfig::info("calculated distance");
561 //now break at every point - calculate new derivatives each time
564 //must be sure that the break points are 3 or more apart
565 //then must also store the breaks which are not smooth, etc.
566 //and figure out tangents between there
568 //for each pair of break points (as long as they are far enough apart) recursively subdivide stuff
569 //ignore the detected intermediate points
571 unsigned int i0=0,i3=0,is=0;
576 Real errortol = smoothness*pixelwidth; //???? what the hell should this value be
581 //intemp = f; //don't want to smooth out the corners
583 bool breaktan = false, setwidth;
584 a.set_split_tangent_flag(false);
585 //a.set_width(width);
588 setwidth = (f.size() == f_w.size());
590 for(j = 0; j < (int)brk.size() - 1; ++j)
592 //for b[j] to b[j+1] subdivide and stuff
596 unsigned int size = i3-i0+1; //must include the end points
600 ftemp.assign(f.begin()+i0, f.begin()+i3+1);
602 gaussian_blur_3(ftemp.begin(),ftemp.end(),false);
606 // Wondering whether the modification of the df vector
607 // using a char* pointer and pointer arithmetric was safe,
610 // http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2369.pdf tells me:
612 // 23.2.5 Class template vector [vector]
614 // [...] The elements of a vector are stored contiguously,
615 // meaning that if v is a vector<T,Allocator> where T is
616 // some type other than bool, then it obeys the identity
617 // &v[n] == &v[0] + n for all 0 <= n < v.size().
619 GetFirstDerivatives(ftemp,0,size,(char*)&df[0],sizeof(df[0]));
621 //GetSimpleDerivatives(ftemp,0,size,df,0,di);
622 //< don't have to worry about indexing stuff as it is all being taken care of right now
623 //preproceval += timer();
626 work.resize(size*2-1); //guarantee that all points will be tessellated correctly (one point in between every 2 adjacent points)
628 //if size of work is size*2-1, the step size should be 1/(size*2 - 2)
629 //Real step = 1/(Real)(size*2 - 1);
631 //start off with break points as indices
633 curind.push_back(cpindex(i0,di[i3]-di[i0],0)); //0 error because no curve on the left
634 curind.push_back(cpindex(i3,di[i3]-di[i0],-1)); //error needs to be reevaluated
635 done = false; //we want to loop
637 unsigned int dcount = 0;
639 //while there are still enough points between us, and the error is too high subdivide (and invalidate neighbors that share tangents)
642 //tessellate all curves with invalid error values
646 /*numtess += */tessellate_curves(curind,f,df,work);
647 //tesseval += timer();
649 //now get all error values
651 for(i = 1; i < (int)curind.size(); ++i)
653 if(curind[i].error < 0) //must have been retessellated, so now recalculate error value
655 //evaluate error from points (starting at current index)
656 int size = curind[i].curind - curind[i-1].curind + 1;
657 curind[i].error = CurveError(&f[curind[i-1].curind], size,
658 work,(curind[i-1].curind - i0)*2,(curind[i].curind - i0)*2+1);
660 /*if(curind[i].error > 1.0e5)
662 synfig::info("Holy crap %d-%d error %f",curind[i-1].curind,curind[i].curind,curind[i].error);
663 curind[i].error = -1;
664 numtess += tessellate_curves(curind,f,df,work);
665 curind[i].error = CurveError(&f[curind[i-1].curind], size,
666 work,0,work.size());//(curind[i-1].curind - i0)*2,(curind[i].curind - i0)*2+1);
671 //erroreval += timer();
676 //check each error to see if it's too big, if so, then subdivide etc.
677 int indsize = (int)curind.size();
678 Real maxrelerror = 0;
679 int maxi = -1;//, numpoints;
682 //get the maximum error and split there
683 for(i = 1; i < indsize; ++i)
685 //numpoints = curind[i].curind - curind[i-1].curind + 1;
687 if(curind[i].error > maxrelerror && curind[i].curind - curind[i-1].curind > 2) //only accept if it's valid
689 maxrelerror = curind[i].error;
694 //split if error is too great
695 if(maxrelerror > errortol)
697 //add one to the left etc
698 unsigned int ibase = curind[maxi-1].curind, itop = curind[maxi].curind,
699 ibreak = (ibase + itop)/2;
702 assert(ibreak < f.size());
704 //synfig::info("Split %d -%d- %d, error: %f", ibase,ibreak,itop,maxrelerror);
708 //invalidate current error of the changed tangents and add an extra segment
709 //enforce minimum tangents property
710 curind[maxi].error = -1;
711 curind[maxi-1].error = -1;
712 if(maxi+1 < indsize) curind[maxi+1].error = -1; //if there is a curve segment beyond this it will be effected as well
714 scale = di[itop] - di[ibreak];
715 scale2 = maxi+1 < indsize ? di[curind[maxi+1].curind] - di[itop] : scale; //to the right valid?
716 curind[maxi].tangentscale = std::min(scale, scale2);
718 scale = di[ibreak] - di[ibase];
719 scale2 = maxi >= 2 ? di[ibase] - di[curind[maxi-2].curind] : scale; // to the left valid -2 ?
720 curind[maxi-1].tangentscale = std::min(scale, scale2);
722 scale = std::min(di[ibreak] - di[ibase], di[itop] - di[ibreak]);
724 curind.insert(curind.begin()+maxi,cpindex(ibreak, scale, -1));
725 //curind.push_back(cpindex(ibreak, scale, -1));
726 //std::sort(curind.begin(), curind.end());
732 //spliteval += timer();
737 //insert the last point too (just set tangent for now
738 is = curind[0].curind;
740 //first point inherits current tangent status
742 if(v.mag_squared() > EPSILON)
743 v *= (curind[0].tangentscale/v.mag());
747 else a.set_tangent2(v);
750 if(setwidth)a.set_width(f_w[is]);
753 a.set_split_tangent_flag(false); //won't need to break anymore
756 for(i = 1; i < (int)curind.size()-1; ++i)
758 is = curind[i].curind;
760 //first point inherits current tangent status
762 if(v.mag_squared() > EPSILON)
763 v *= (curind[i].tangentscale/v.mag());
765 a.set_tangent(v); // always inside, so guaranteed to be smooth
767 if(setwidth)a.set_width(f_w[is]);
772 //set the last point's data
773 is = curind.back().curind; //should already be this
776 if(v.mag_squared() > EPSILON)
777 v *= (curind.back().tangentscale/v.mag());
780 a.set_split_tangent_flag(true);
783 //will get the vertex and tangent 2 from next round
787 a.set_split_tangent_flag(false);
789 a.set_width(f_w[i3]);
792 /*etl::clock::value_type totaltime = total(),
793 misctime = totaltime - initialprocess - curveval - breakeval - disteval
794 - preproceval - tesseval - erroreval - spliteval;
797 "Curve Convert Profile:\n"
798 "\tInitial Preprocess: %f\n"
799 "\tCurvature Calculation: %f\n"
800 "\tBreak Calculation: %f\n"
801 "\tDistance Calculation: %f\n"
802 " Algorithm: (numtimes,totaltime)\n"
803 "\tPreprocess step: (%d,%f)\n"
804 "\tTessellation step: (%d,%f)\n"
805 "\tError step: (%d,%f)\n"
806 "\tSplit step: (%d,%f)\n"
807 " Num Input: %d, Num Output: %d\n"
808 " Total time: %f, Misc time: %f\n",
809 initialprocess, curveval,breakeval,disteval,
810 numpre,preproceval,numtess,tesseval,numerror,erroreval,numsplit,spliteval,
811 in.size(),out.size(),
812 totaltime,misctime);*/
818 void synfigapp::BLineConverter::EnforceMinWidth(std::list<synfig::BLinePoint> &bline, synfig::Real min_pressure)
820 std::list<synfig::BLinePoint>::iterator i = bline.begin(),
823 for(i = bline.begin(); i != end; ++i)
824 if(i->get_width() < min_pressure)
825 i->set_width(min_pressure);