Adding upstream version 1.3.upstream/1.3upstream

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

1 files changed, 404 insertions, 0 deletions

diff --git a/src/toys/sb-to-bez.cpp b/src/toys/sb-to-bez.cpp
new file mode 100644
index 0000000..337180f
--- /dev/null
+++ b/src/toys/sb-to-bez.cpp
@@ -0,0 +1,404 @@
+/* 
+ * sb-to-bez Toy - Tests conversions from sbasis to cubic bezier.
+ *
+ * Copyright 2007 jf barraud.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ *
+ */
+
+// mainly experimental atm...
+// do not expect to find anything understandable here atm. 
+
+#include <2geom/d2.h>
+#include <2geom/sbasis.h>
+#include <2geom/sbasis-geometric.h>
+#include <2geom/basic-intersection.h>
+
+#include <toys/path-cairo.h>
+#include <toys/toy-framework-2.h>
+
+#define ZERO 1e-7
+
+using std::vector;
+using namespace Geom;
+using namespace std;
+
+#include <stdio.h>
+#include <gsl/gsl_poly.h>
+
+void cairo_pw(cairo_t *cr, Piecewise<SBasis> p, double hscale=1., double vscale=1.) {
+    for(unsigned i = 0; i < p.size(); i++) {
+        D2<SBasis> B;
+        B[0] = Linear(150+p.cuts[i]*hscale, 150+p.cuts[i+1]*hscale);
+        B[1] = Linear(450) - p[i]*vscale;
+        cairo_d2_sb(cr, B);
+    }
+}
+
+//===================================================================================
+
+D2<SBasis>
+naive_sb_seg_to_bez(Piecewise<D2<SBasis> > const &M,double t0,double t1){
+
+    Piecewise<D2<SBasis> > dM = derivative(M); 
+    Point M0  = M(t0);
+    Point dM0 = dM(t0)*(t1-t0);
+    Point M1  = M(t1);
+    Point dM1 = dM(t1)*(t1-t0);
+    D2<SBasis> result;
+    for (unsigned dim=0; dim<2; dim++){
+        SBasis r(2, Linear());
+        r[0] = Linear(M0[dim],M1[dim]);
+        r[1] = Linear(M0[dim]-M1[dim]+dM0[dim],-(M0[dim]-M1[dim]+dM1[dim]));
+        result[dim] = r;
+    }
+    return result;
+}
+
+D2<SBasis>
+sb_seg_to_bez(Piecewise<D2<SBasis> > const &M,double t0,double t1){
+    Point M0,dM0,d2M0,M1,dM1,d2M1,A0,V0,A1,V1;
+    Piecewise<D2<SBasis> > dM,d2M;
+    dM=derivative(M);
+    d2M=derivative(dM);
+    M0  =M(t0);
+    M1  =M(t1);
+    dM0 =dM(t0);
+    dM1 =dM(t1);
+    d2M0=d2M(t0);
+    d2M1=d2M(t1);
+    A0=M(t0);
+    A1=M(t1);
+    
+    std::vector<D2<SBasis> > candidates = cubics_fitting_curvature(M0,M1,dM0,dM1,d2M0,d2M1);
+    if (candidates.empty()){
+        return D2<SBasis>(SBasis(M0[X],M1[X]),SBasis(M0[Y],M1[Y])) ;
+    }
+    double maxlength = -1;
+    unsigned best = 0;
+    for (unsigned i=0; i<candidates.size(); i++){
+        double l = length(candidates[i]);
+        if ( l < maxlength || maxlength < 0 ){
+            maxlength = l;
+            best = i;
+        }
+    }
+    return candidates[best];
+}
+#include <2geom/sbasis-to-bezier.h>
+
+int goal_function_type = 0;
+
+double goal_function(Piecewise<D2<SBasis> >const &A, 
+                     Piecewise<D2<SBasis> >const&B) {
+    if(goal_function_type) {
+        OptInterval bnds = bounds_fast(dot(derivative(A), rot90(derivative(B))));
+        //double h_dist = bnds.dimensions().length();
+//0 is in the rect!, TODO:gain factor ~2 for free.
+// njh: not really, the benefit is actually rather small.
+        double h_dist = 0;
+        if(bnds)
+            h_dist = bnds->extent();
+        return h_dist ;
+    } else {
+        Rect bnds = *bounds_fast(A - B);
+        return max(bnds.min().length(), bnds.max().length());
+    }  
+}
+
+int recursive_curvature_fitter(cairo_t* cr, Piecewise<D2<SBasis> > const &f, double t0, double t1, double precision) {
+      if (t0>=t1) return 0;//TODO: fix me...
+      if (t0+0.001>=t1) return 0;//TODO: fix me...
+      
+      //TODO: don't re-compute derivative(f) at each try!!
+      D2<SBasis> k_bez = sb_seg_to_bez(f,t0,t1);      
+
+      if(k_bez[0].size() > 1 and k_bez[1].size() > 1) {
+          Piecewise<SBasis> s = arcLengthSb(k_bez);
+          s *= (t1-t0)/arcLengthSb(k_bez).segs.back().at1();
+          s += t0;
+          double h_dist = goal_function(compose(f,s), Piecewise<D2<SBasis> >(k_bez));
+          if(h_dist < precision) {
+              cairo_save(cr);
+              cairo_set_line_width (cr, 0.93);
+              cairo_set_source_rgba (cr, 0.7, 0.0, 0.0, 1);
+              draw_handle(cr, k_bez.at0());
+              cairo_d2_sb(cr, k_bez);
+              cairo_stroke(cr);
+              cairo_restore(cr);
+              return 1;
+          }
+      }
+      //TODO: find a better place where to cut (at the worst fit?).
+      return recursive_curvature_fitter(cr, f, t0, (t0+t1)/2, precision) +
+          recursive_curvature_fitter(cr, f, (t0+t1)/2, t1, precision);
+}
+
+double single_curvature_fitter(Piecewise<D2<SBasis> > const &f, double t0, double t1) {
+      if (t0>=t1) return 0;//TODO: fix me...
+      if (t0+0.001>=t1) return 0;//TODO: fix me...
+      
+      D2<SBasis> k_bez = sb_seg_to_bez(f,t0,t1);      
+
+      if(k_bez[0].size() > 1 and k_bez[1].size() > 1) {
+          Piecewise<SBasis> s = arcLengthSb(k_bez);
+          s *= (t1-t0)/arcLengthSb(k_bez).segs.back().at1();
+          s += t0;
+          return goal_function(compose(f,s), Piecewise<D2<SBasis> >(k_bez));
+      }
+      return 1e100;
+}
+
+struct quadratic_params
+{
+    Piecewise<D2<SBasis> > const *f;
+    double t0, precision;
+};
+     
+
+double quadratic (double x, void *params) {
+    struct quadratic_params *p 
+        = (struct quadratic_params *) params;
+     
+    return single_curvature_fitter(*p->f, p->t0, x) - p->precision;
+}
+     
+#include <stdio.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_roots.h>
+     
+     
+int sequential_curvature_fitter(cairo_t* cr, Piecewise<D2<SBasis> > const &f, double t0, double t1, double precision) {
+    if(t0 >= t1) return 0;
+    
+    double r = t1;
+    if(single_curvature_fitter(f, t0, t1) > precision) {
+        int status;
+        int iter = 0, max_iter = 100;
+        const gsl_root_fsolver_type *T;
+        gsl_root_fsolver *s;
+        gsl_function F;
+        struct quadratic_params params = {&f, t0, precision};
+     
+        F.function = &quadratic;
+        F.params = &params;
+     
+        T = gsl_root_fsolver_brent;
+        s = gsl_root_fsolver_alloc (T);
+        gsl_root_fsolver_set (s, &F, t0, t1);
+     
+        do
+        {
+            iter++;
+            status = gsl_root_fsolver_iterate (s);
+            r = gsl_root_fsolver_root (s);
+            double x_lo = gsl_root_fsolver_x_lower (s);
+            double x_hi = gsl_root_fsolver_x_upper (s);
+            status = gsl_root_test_interval (x_lo, x_hi,
+                                             0, 0.001);
+     
+     
+        }
+        while (status == GSL_CONTINUE && iter < max_iter);
+    
+        double x_lo = gsl_root_fsolver_x_lower (s);
+        double x_hi = gsl_root_fsolver_x_upper (s);
+        printf ("%5d [%.7f, %.7f] %.7f %.7f\n",
+                iter, x_lo, x_hi,
+                r, 
+                x_hi - x_lo);
+        gsl_root_fsolver_free (s);
+    }    
+    D2<SBasis> k_bez = sb_seg_to_bez(f,t0,r);      
+    
+    cairo_save(cr);
+    cairo_set_line_width (cr, 0.93);
+    cairo_set_source_rgba (cr, 0.7, 0.0, 0.0, 1);
+    draw_handle(cr, k_bez.at0());
+    cairo_d2_sb(cr, k_bez);
+    cairo_stroke(cr);
+    cairo_restore(cr);
+    
+    if(r < t1)
+        return sequential_curvature_fitter(cr, f, r, t1, precision) + 1;
+    return 1;
+}
+
+
+class SbToBezierTester: public Toy {
+    //std::vector<Slider> sliders;
+    std::vector<PointSetHandle*> path_psh;
+    PointHandle adjuster, adjuster2;
+    std::vector<Toggle> toggles;
+
+  void draw(cairo_t *cr, std::ostringstream *notify, int width, int height, bool save, std::ostringstream *timer_stream) override {
+      cairo_save(cr);
+      for(unsigned i = 1; i < path_psh.size(); i++)
+          path_psh[i-1]->pts.back() = path_psh[i]->pts[0];
+      Piecewise<D2<SBasis> > f_as_pw(path_psh[0]->asBezier());
+    for(unsigned i = 1; i < path_psh.size(); i++) {
+        f_as_pw.push_seg(path_psh[i]->asBezier());
+    }
+      //f=handles_to_sbasis(handles.begin(), SIZE-1);
+      adjuster.pos[1]=450;
+      adjuster.pos[0]=std::max(adjuster.pos[0],150.);
+      adjuster.pos[0]=std::min(adjuster.pos[0],450.);
+      double t0=0;//(adjuster.pos[0]-150)/300;
+      double t1=(adjuster.pos[0]-150)/300;
+      //if (t0>t1) {double temp=t0;t0=t1;t1=temp;}
+
+      cairo_set_source_rgba (cr, 0., 0., 0., 1);
+      cairo_set_line_width (cr, 0.5);
+      cairo_pw_d2_sb(cr, f_as_pw);
+      cairo_stroke(cr);
+      if (t0==t1) return;//TODO: fix me...
+#if 0
+      if(0) {
+      Piecewise<D2<SBasis> > g = f_as_pw;
+      cairo_set_line_width (cr, 1);
+      cairo_set_source_rgba (cr, 0., 0., 0.9, .7);
+      double error=0;
+
+      cairo_set_line_width (cr, 1);
+      cairo_set_source_rgba (cr, 0.9, 0., 0., .7);
+      D2<SBasis> naive_bez = naive_sb_seg_to_bez(g,0,t1);
+      cairo_d2_sb(cr, naive_bez);
+      cairo_stroke(cr);
+
+      adjuster2.pos[0]=150;
+      adjuster2.pos[1]=std::min(std::max(adjuster2.pos[1],150.),450.);
+
+      double scale0=(450-adjuster2.pos[1])/150;
+
+      cairo_set_line_width (cr, 1);
+      cairo_set_source_rgba (cr, 0.7, 0., 0.7, .7);
+      D2<SBasis> k_bez = sb_seg_to_bez(g,t0,t1);
+      cairo_d2_sb(cr, k_bez);
+      cairo_stroke(cr);
+      double h_a_t = 0, h_b_t = 0;
+      
+      double h_dist = hausdorfl( k_bez, f, 1e-6, &h_a_t, &h_b_t);
+      {
+          Point At = k_bez(h_a_t);
+          Point Bu = f(h_b_t);
+          cairo_move_to(cr, At);
+          cairo_line_to(cr, Bu);
+          draw_handle(cr, At);
+          draw_handle(cr, Bu);
+          cairo_save(cr);
+          cairo_set_line_width (cr, 0.3);
+          cairo_set_source_rgba (cr, 0.7, 0.0, 0.0, 1);
+          cairo_stroke(cr);
+          cairo_restore(cr);
+      }
+      *notify << "Move handle 6 to set the segment to be approximated by cubic bezier.\n";
+      *notify << " -red:  bezier approx derived from parametrization.\n";
+      *notify << " -blue: bezier approx derived from curvature.\n";
+      *notify << "      max distance (to original): "<<h_dist<<"\n";
+      }
+#endif
+
+
+      f_as_pw = arc_length_parametrization(f_as_pw);
+      adjuster2.pos[0]=150;
+      adjuster2.pos[1]=std::min(std::max(adjuster2.pos[1],150.),450.);
+      cairo_move_to(cr, 150, 150);
+      cairo_line_to(cr, 150, 450);
+      cairo_stroke(cr);
+      ostringstream val_s;
+      double scale0=(450-adjuster2.pos[1])/300;
+      double curve_precision = pow(10, scale0*5-2);
+      val_s << curve_precision;
+      draw_text(cr, adjuster2.pos, val_s.str().c_str());
+      
+    int segs = 0;
+    goal_function_type = toggles[1].on;
+      if(toggles[0].on)
+          segs = sequential_curvature_fitter(cr, f_as_pw, 0, f_as_pw.cuts.back(), curve_precision);
+      else {
+          segs = recursive_curvature_fitter(cr, f_as_pw, 0, f_as_pw.cuts.back(),curve_precision);
+      }
+      Geom::PathVector vpt = path_from_piecewise(f_as_pw, curve_precision, true);
+      unsigned default_number_curves = 0;
+      for(const auto & i : vpt) {
+          default_number_curves += i.size();
+      }
+      
+      *notify << "      segments from default algorithm: "<< default_number_curves <<"\n";
+      *notify << "      total segments: "<< segs <<"\n";
+      cairo_restore(cr);
+      Point p(25, height - 100), d(50,25);
+      toggles[0].bounds = Rect(p,     p + d);
+      p+= Point(75, 0);
+      toggles[1].bounds = Rect(p,     p + d);
+      draw_toggles(cr, toggles);
+      Toy::draw(cr, notify, width, height, save,timer_stream);
+  }
+  
+public:
+    void key_hit(GdkEventKey *e) override {
+        if(e->keyval == 's') toggles[0].toggle();
+        redraw();
+    }
+    void mouse_pressed(GdkEventButton* e) override {
+        toggle_events(toggles, e);
+        Toy::mouse_pressed(e);
+    }
+    SbToBezierTester() {
+      //if(handles.empty()) {
+        for(int j = 0; j < 3; j++) {
+            path_psh.push_back(new PointSetHandle());
+            for(unsigned i = 0; i < 6; i++)
+                path_psh.back()->push_back(150+300*uniform(),150+300*uniform());
+            handles.push_back(path_psh.back());
+        }
+      adjuster.pos = Geom::Point(150+300*uniform(),150+300*uniform());
+      handles.push_back(&adjuster);
+      adjuster2.pos = Geom::Point(150,300);
+      handles.push_back(&adjuster2);
+      toggles.emplace_back("Seq", true);
+      toggles.emplace_back("Linfty", true);
+      //}
+    //sliders.push_back(Slider(0.0, 1.0, 0.0, 0.0, "t"));
+    //handles.push_back(&(sliders[0]));
+  }
+};
+
+int main(int argc, char **argv) {
+    init(argc, argv, new SbToBezierTester);
+    return 0;
+}
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype = cpp:expandtab:shiftwidth = 4:tabstop = 8:softtabstop = 4:encoding = utf-8:textwidth = 99 :