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
9 ** Copyright (c) 2007 Chris Moore
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.
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.
22 /* ========================================================================= */
24 /* === H E A D E R S ======================================================= */
33 #include "blineconvert.h"
35 #include <ETL/gaussian>
36 #include <ETL/hermite>
40 #include <synfig/general.h>
47 /* === U S I N G =========================================================== */
51 using namespace synfig;
53 /* === M A C R O S ========================================================= */
55 #define EPSILON (1e-10)
57 /* === G L O B A L S ======================================================= */
59 /* === P R O C E D U R E S ================================================= */
61 /* === M E T H O D S ======================================================= */
64 //Derivative Functions for numerical approximation
66 //bias == 0 will get F' at f3, bias < 0 will get F' at f1, and bias > 0 will get F' at f5
68 inline void FivePointdt(T &df, const T &f1, const T &f2, const T &f3, const T &f4, const T &f5, int bias)
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);
75 df = (f1*3 - f2*16 + f3*36 - f4*48 + f5*25)*(1/12.0f);
79 inline void ThreePointdt(T &df, const T &f1, const T &f2, const T &f3, int bias)
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);
86 df = (f1 - f2*4 + f3*3)*(1/2.0f);
89 // template < class T >
90 // inline void ThreePointddt(T &df, const T &f1, const T &f2, const T &f3, int bias)
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);
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)
105 // //df = (- f1 + f2*16 - f3*30 + f4*16 - f5)*(1/12.0f);
106 // }/*else if(bias < 0)
109 // df = (f1*7 - f2*26*4 + f3*19*6 - f4*14*4 + f5*11)*(1/12.0f);
113 // df = (f1*3 - f2*16 + f3*36 - f4*48 + f5*25)*(1/12.0f);
115 // //side ones don't work, use 3 point
118 // //implement an arbitrary derivative
120 // template < class T >
121 // void DerivativeApprox(T &df, const T f[], const Real t[], int npoints, int indexval)
124 // Lj(x) = PI_i!=j (x - xi) / PI_i!=j (xj - xi)
126 // so Lj'(x) = SUM_k PI_i!=j|k (x - xi) / PI_i!=j (xj - xi)
129 // unsigned int i,j,k,i0,i1;
131 // Real Lpj,mult,div,tj;
132 // Real tval = t[indexval];
135 // for(j=0;j<npoints;++j)
141 // for(k=0;k<npoints;++k)
143 // if(k != j) //because there is no summand for k == j, since that term is missing from the original equation
146 // for(i=0;i<npoints;++i)
150 // mult *= tval - t[i];
154 // Lpj += mult; //add into the summation
156 // //since the ks follow the exact pattern we need for the divisor (use that too)
161 // //get the actual coefficient
164 // //add it in to the equation
169 //END numerical derivatives
171 // template < class T >
172 // inline int sign(T f, T tol)
174 // if(f < -tol) return -1;
175 // if(f > tol) return 1;
179 void GetFirstDerivatives(const std::vector<synfig::Point> &f, unsigned int left, unsigned int right, char *out, unsigned int dfstride)
181 unsigned int current = left;
185 else if(right - left == 2)
187 synfig::Vector v = f[left+1] - f[left];
189 //set both to the one we want
190 *(synfig::Vector*)out = v;
192 *(synfig::Vector*)out = v;
195 else if(right - left < 6/*5*/) //should use 3 point
197 //left then middle then right
198 ThreePointdt(*(synfig::Vector*)out,f[left+0], f[left+1], f[left+2], -1);
202 for(;current < right-1; current++, out += dfstride)
203 ThreePointdt(*(synfig::Vector*)out,f[current-1], f[current], f[current+1], 0);
205 ThreePointdt(*(synfig::Vector*)out,f[right-3], f[right-2], f[right-1], 1);
209 }else //can use 5 point
211 //left 2 then middle bunch then right two
212 //may want to use 3 point for inner edge ones
214 FivePointdt(*(synfig::Vector*)out,f[left+0], f[left+1], f[left+2], f[left+3], f[left+4], -2);
216 FivePointdt(*(synfig::Vector*)out,f[left+1], f[left+2], f[left+3], f[left+4], f[left+5], -1);
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);
223 FivePointdt(*(synfig::Vector*)out,f[right-6], f[right-5], f[right-4], f[right-3], f[right-2], 2);
225 FivePointdt(*(synfig::Vector*)out,f[right-5], f[right-4], f[right-3], f[right-2], f[right-1], 1);
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*/)
236 int offset = 2; //df = 1/2 (f[i+o]-f[i-o])
238 assert((int)df.size() >= right-left+outleft); //must be big enough
240 for(i = left; i < right; ++i)
242 //right now indices (figure out distance later)
243 i1 = std::max(left,i-offset);
244 i2 = std::max(left,i+offset);
246 df[outleft++] = (f[i2] - f[i1])*0.5f;
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)
253 if(right-left < 2) return -1;
257 //get distances to each point
259 //synfig::Vector v,vt;
260 //synfig::Point p1,p2;
262 std::vector<synfig::Point>::const_iterator it;//,end = work.begin()+right;
264 //unsigned int size = work.size();
266 //for each line, get distance
268 for(i = 0; i < (int)n; ++i)
274 it = work.begin()+left;
275 //p2 = *it++; //put it at left+1
276 for(j = left/*+1*/; j < right; ++j,++it)
284 dtemp = v.mag_squared() > 1e-12 ? (vt*v)/v.mag_squared() : 0; //get the projected time value for the current line
286 //get distance to line segment with the time value clamped 0-1
287 if(dtemp >= 1) //use p+v
289 vt += v; //makes it pp - (p+v)
290 }else if(dtemp > 0) //use vt-proj
292 vt -= v*dtemp; // vt - proj_v(vt) //must normalize the projection vector to work
296 dtemp = vt.mag_squared();*/
298 dtemp = (pi - *it).mag_squared();
303 //accumulate the points' min distance from the curve
310 typedef synfigapp::BLineConverter::cpindex cpindex;
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)
318 etl::hermite<Point> curve;
321 std::vector<cpindex>::const_iterator j = inds.begin(),j2, end = inds.end();
323 unsigned int ibase = inds[0].curind;
326 for(; j != end; j2 = j++)
328 //if this curve has invalid error (in j) then retessellate its work points (requires reparametrization, etc.)
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
336 Real t, dt = 1/(Real)(n*2-2); //assuming that they own only n points
338 //start at first intermediate
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)
345 //build hermite curve, it's easier
348 curve.t1() = df[i0-ibase] * (df[i0-ibase].mag_squared() > 1e-4
349 ? j2->tangentscale/df[i0-ibase].mag()
351 curve.t2() = df[i3-ibase] * (df[i3-ibase].mag_squared() > 1e-4
352 ? j->tangentscale/df[i3-ibase].mag()
356 //MUST include the end point (since we are ignoring left one)
357 for(; k < kend; ++k, t += dt)
362 work[k] = curve(1); //k == kend, t == 1 -> c(t) == p2
370 synfigapp::BLineConverter::BLineConverter()
378 synfigapp::BLineConverter::clear()
385 break_tangents.clear();
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)
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;*/
404 if (points_in.size() < 2)
409 //removing digitization error harder than expected
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
418 //remove duplicate points - very bad for fitting
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();
427 if (points_in.size() == widths_in.size())
429 for(bool first = true; point_iter != end; ++point_iter,++width_iter)
430 if (first || *point_iter != c) // eliminate duplicate points
433 point_cache.push_back(c = *point_iter);
434 width_cache.push_back(*width_iter);
438 for(;point_iter != end; ++point_iter)
439 if(*point_iter != c) // eliminate duplicate points
440 point_cache.push_back(c = *point_iter);
442 //initialprocess = timer();
444 if (point_cache.size() < 7)
446 info("only %d unique points - giving up", point_cache.size());
450 //get curvature information
454 int i_this, i_prev, i_next;
455 synfig::Vector v_prev, v_next;
457 curvature.resize(point_cache.size());
458 curvature.front() = curvature.back() = 1;
460 for (i_this = 1; i_this < (int)point_cache.size()-1; i_this++)
462 i_prev = std::max(0, i_this-2);
463 i_next = std::min((int)(point_cache.size()-1), i_this+2);
465 v_prev = point_cache[i_this] - point_cache[i_prev];
466 v_next = point_cache[i_next] - point_cache[i_this];
468 curvature[i_this] = (v_prev*v_next) / (v_prev.mag()*v_next.mag());
472 //curveval = timer();
473 //synfig::info("calculated curvature");
475 //find corner points and interpolate inside those
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)
485 Real sharpest_curvature = 1;
487 break_tangents.push_back(0);
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)
495 if (curvature[i] < tol)
497 if(curvature[i] < sharpest_curvature)
499 sharpest_curvature = curvature[i];
503 else if (sharpest_i > 0)
505 // don't have 2 corners too close to each other
506 if (sharpest_i >= last + 8) //! \todo make this configurable
508 //synfig::info("break: %d-%d",sharpest_i+1,curvature.size());
509 break_tangents.push_back(sharpest_i);
513 sharpest_curvature = 1;
517 break_tangents.push_back(i);
519 // this section causes bug 1892566 if enabled
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)
524 Point p = point_cache[break_tangents.front()];
527 for (i = 1; i < break_tangents.size()-1; ++i) //do not want to include end point...
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...
532 //want to erase all points before...
534 break_tangents.erase(break_tangents.begin(),break_tangents.begin()+i-1);
537 p = point_cache[break_tangents.back()];
538 for(i = break_tangents.size()-2; i > 0; --i) //start at one in from the end
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
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
548 //breakeval = timer();
549 //synfig::info("found break points: %d",break_tangents.size());
551 //get the distance calculation of the entire curve (for tangent scaling)
557 p1=p2=point_cache[0];
559 cum_dist.resize(point_cache.size()); this_dist.resize(point_cache.size());
561 for(unsigned int i = 0; i < point_cache.size();)
563 d += (this_dist[i] = (p2-p1).mag());
567 //! \todo is this legal? it reads off the end of the vector
571 //disteval = timer();
572 //synfig::info("calculated distance");
574 //now break at every point - calculate new derivatives each time
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
581 //for each pair of break points (as long as they are far enough apart) recursively subdivide stuff
582 //ignore the detected intermediate points
584 unsigned int i0=0,i3=0,is=0;
589 Real errortol = smoothness*pixelwidth; //???? what should this value be
594 //intemp = f; //don't want to smooth out the corners
596 bool breaktan = false, setwidth;
597 a.set_split_tangent_flag(false);
598 //a.set_width(width);
601 setwidth = (point_cache.size() == width_cache.size());
603 for(j = 0; j < (int)break_tangents.size() - 1; ++j)
605 //for b[j] to b[j+1] subdivide and stuff
606 i0 = break_tangents[j];
607 i3 = break_tangents[j+1];
609 unsigned int size = i3-i0+1; //must include the end points
613 ftemp.assign(point_cache.begin()+i0, point_cache.begin()+i3+1);
615 gaussian_blur_3(ftemp.begin(),ftemp.end(),false);
619 // Wondering whether the modification of the deriv vector
620 // using a char* pointer and pointer arithmetric was safe,
623 // http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2369.pdf tells me:
625 // 23.2.5 Class template vector [vector]
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().
632 GetFirstDerivatives(ftemp,0,size,(char*)&deriv[0],sizeof(deriv[0]));
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();
639 work.resize(size*2-1); //guarantee that all points will be tessellated correctly (one point in between every 2 adjacent points)
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);
644 //start off with break points as indices
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
650 unsigned int dcount = 0;
652 //while there are still enough points between us, and the error is too high subdivide (and invalidate neighbors that share tangents)
655 //tessellate all curves with invalid error values
656 work[0] = point_cache[i0];
659 /*numtess += */tessellate_curves(curind,point_cache,deriv,work);
660 //tesseval += timer();
662 //now get all error values
664 for(i = 1; i < (int)curind.size(); ++i)
666 if(curind[i].error < 0) //must have been retessellated, so now recalculate error value
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);
673 /*if(curind[i].error > 1.0e5)
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);
684 //erroreval += timer();
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;
695 //get the maximum error and split there
696 for(i = 1; i < indsize; ++i)
698 //numpoints = curind[i].curind - curind[i-1].curind + 1;
700 if(curind[i].error > maxrelerror && curind[i].curind - curind[i-1].curind > 2) //only accept if it's valid
702 maxrelerror = curind[i].error;
707 //split if error is too great
708 if(maxrelerror > errortol)
710 //add one to the left etc
711 unsigned int ibase = curind[maxi-1].curind, itop = curind[maxi].curind,
712 ibreak = (ibase + itop)/2;
715 assert(ibreak < point_cache.size());
717 //synfig::info("Split %d -%d- %d, error: %f", ibase,ibreak,itop,maxrelerror);
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
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);
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);
735 scale = std::min(cum_dist[ibreak] - cum_dist[ibase], cum_dist[itop] - cum_dist[ibreak]);
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());
745 //spliteval += timer();
750 //insert the last point too (just set tangent for now
751 is = curind[0].curind;
753 //first point inherits current tangent status
755 if(v.mag_squared() > EPSILON)
756 v *= (curind[0].tangentscale/v.mag());
760 else a.set_tangent2(v);
762 a.set_vertex(point_cache[is]);
763 if(setwidth)a.set_width(width_cache[is]);
765 blinepoints_out.push_back(a);
766 a.set_split_tangent_flag(false); //won't need to break anymore
769 for(i = 1; i < (int)curind.size()-1; ++i)
771 is = curind[i].curind;
773 //first point inherits current tangent status
775 if(v.mag_squared() > EPSILON)
776 v *= (curind[i].tangentscale/v.mag());
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]);
782 blinepoints_out.push_back(a);
785 //set the last point's data
786 is = curind.back().curind; //should already be this
789 if(v.mag_squared() > EPSILON)
790 v *= (curind.back().tangentscale/v.mag());
793 a.set_split_tangent_flag(true);
796 //will get the vertex and tangent 2 from next round
799 a.set_vertex(point_cache[i3]);
800 a.set_split_tangent_flag(false);
802 a.set_width(width_cache[i3]);
803 blinepoints_out.push_back(a);
805 /*etl::clock::value_type totaltime = total(),
806 misctime = totaltime - initialprocess - curveval - breakeval - disteval
807 - preproceval - tesseval - erroreval - spliteval;
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);*/
831 void synfigapp::BLineConverter::EnforceMinWidth(std::list<synfig::BLinePoint> &bline, synfig::Real min_pressure)
833 std::list<synfig::BLinePoint>::iterator i = bline.begin(),
836 for(i = bline.begin(); i != end; ++i)
837 if(i->get_width() < min_pressure)
838 i->set_width(min_pressure);