Adding upstream version 1.3.upstream/1.3upstream

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

1 files changed, 510 insertions, 0 deletions

diff --git a/src/toys/implicit-toy.cpp b/src/toys/implicit-toy.cpp
new file mode 100644
index 0000000..c90c082
--- /dev/null
+++ b/src/toys/implicit-toy.cpp
@@ -0,0 +1,510 @@
+#include <toys/path-cairo.h>
+#include <toys/toy-framework-2.h>
+
+#include <2geom/geom.h>
+#include <2geom/d2.h>
+#include <2geom/polynomial.h>
+#include <2geom/sbasis-poly.h>
+#include <2geom/transforms.h>
+
+#include <2geom/symbolic/implicit.h>
+
+#include <aa.h>
+
+#include <algorithm>
+#include <ctime>
+#include <functional>
+
+
+using namespace Geom;
+
+
+
+struct PtLexCmp{
+    bool operator()(const Point &a, const Point &b) {
+        return (a[0] < b[0]) || ((a[0] == b[0]) and (a[1] < b[1]));
+    }
+};
+
+//typedef AAF (*implicit_curve_t)(AAF, AAF);
+typedef std::function<AAF (AAF const&, AAF const&)> implicit_curve_t;
+
+// draw ax + by + c = 0
+void draw_line_in_rect(cairo_t*cr, Rect &r, Point n, double c)
+{
+    vector<Geom::Point> result;
+    Point resultp;
+    if(intersects == line_intersection(Point(1, 0), r.left(),
+                       n, c,
+                       resultp) && r[1].contains(resultp[1]))
+    result.push_back(resultp);
+    if(intersects == line_intersection(Point(1, 0), r.right(),
+                       n, c,
+                       resultp) && r[1].contains(resultp[1]))
+    result.push_back(resultp);
+    if(intersects == line_intersection(Point(0, 1), r.top(),
+                       n, c,
+                       resultp) && r[0].contains(resultp[0]))
+    result.push_back(resultp);
+    if(intersects == line_intersection(Point(0, 1), r.bottom(),
+                       n, c,
+                       resultp) && r[0].contains(resultp[0]))
+    result.push_back(resultp);
+    if(result.size() > 2) {
+        std::sort(result.begin(), result.end(), PtLexCmp());
+        vector<Geom::Point>::iterator new_end = std::unique(result.begin(), result.end());
+        result.resize(new_end-result.begin());
+    }
+    if(result.size() == 2)
+    {
+        cairo_move_to(cr, result[0]);
+        cairo_line_to(cr, result[1]);
+        cairo_stroke(cr);
+    }
+}
+
+OptRect tighten(Rect const&r, Point n, Interval lu)
+{
+    vector<Geom::Point> result;
+    Point resultp;
+    for(int i = 0; i < 4; i++)
+    {
+        Point cnr = r.corner(i);
+        double z = dot(cnr, n);
+        if ((z > lu[0]) && (z < lu[1]))
+            result.push_back(cnr);
+    }
+    for(int i = 0; i < 2; i++)
+    {
+        double c = lu[i];
+        if(intersects == line_intersection(Point(1, 0), r.left(),
+                                           n, c,
+                                           resultp) && r[1].contains(resultp[1]))
+            result.push_back(resultp);
+        if(intersects == line_intersection(Point(1, 0), r.right(),
+                                           n, c,
+                                           resultp) && r[1].contains(resultp[1]))
+            result.push_back(resultp);
+        if(intersects == line_intersection(Point(0, 1), r.top(),
+                                           n, c,
+                                           resultp) && r[0].contains(resultp[0]))
+            result.push_back(resultp);
+        if(intersects == line_intersection(Point(0, 1), r.bottom(),
+                                           n, c,
+                                           resultp) && r[0].contains(resultp[0]))
+            result.push_back(resultp);
+    }
+    if(result.size() < 2)
+        return OptRect;
+    Rect nr(result[0], result[1]);
+    for(size_t i = 2; i < result.size(); i++)
+    {
+        nr.expandTo(result[i]);
+    }
+    return intersect(nr, r);
+}
+
+static const unsigned int DEG = 5;
+double bvp[DEG+1][DEG+1]
+    = {{-1, 0.00945115, -4.11799e-05, 1.01365e-07, -1.35037e-10, 7.7868e-14},
+       {0.00837569, -6.24676e-05, 1.96093e-07, -3.09683e-10, 1.95681e-13, 0},
+       {-2.39448e-05, 1.3331e-07, -2.65787e-10, 1.96698e-13, 0, 0},
+       {2.76173e-08, -1.01069e-10, 9.88596e-14, 0, 0, 0},
+       {-1.43584e-11, 2.48433e-14, 0, 0, 0, 0}, {2.49723e-15, 0, 0, 0, 0, 0}};
+
+
+AAF trial_eval(AAF const& x, AAF const& y) {
+//    AAF x = _x/100;
+//    AAF y = _y/100;
+    //return x*x - 1;
+    //return y - pow(x,3);
+    //return y - pow_sample_based(x,2.5);
+    //return y - log_sample_based(x);
+    //return y - log(x);
+    //return y - exp_sample_based(x*log(x));
+    //return y - sqrt(sin(x));
+    //return sqrt(y)*x - sqrt(x) - y - 1;
+    //return y-1/x;
+    //return exp(x)-y;
+    //return sin(x)-y;
+    //return exp_sample_based(x)-y;
+    //return atan(x)-y;
+    //return atan_sample_based(x)-y;
+    //return atanh(x)-y;
+    //return x*y;
+    //return 4*x+3*y-1;
+    //return x*x + y*y - 1;
+    //return sin(x*y) + cos(pow(x, 3)) - atan(x);
+    //return pow((x*x + y*y), 2) - (x*x-y*y);
+    //return x*x-y;
+    //return (x*x*x-y*x)*sin(x) + (x-y*y)*cos(y)-0.5;
+//    return -120.75 +(-64.4688 +(-16.6875 +(0.53125 -0.00390625*y)*y)*y)*y
+//            + (-15.9375 + ( 1.5 +( 4.375 -0.0625*y)*y)*y
+//            + (17 +( 9.5 -0.375*y)*y + (2 + -1*y -1*x)*x)*x)*x;
+
+//    AAF v(0);
+//    for (size_t i = DEG; i > 0; --i)
+//    {
+//        AAF vy(0);
+//        for (size_t j = DEG - i; j > 0; --j)
+//        {
+//            vy += bvp[i][j];
+//            vy *= y;
+//        }
+//        vy += bvp[i][0];
+//        v += vy;
+//        v *= x;
+//    }
+//    AAF vy(0);
+//    for (size_t j = DEG; j > 0; --j)
+//    {
+//        vy += bvp[0][j];
+//        vy *= y;
+//    }
+//    vy += bvp[0][0];
+//    v += vy;
+//    return v;
+
+    int i = DEG;
+    int j = DEG - i;
+    AAF vy(bvp[i][j]);
+    --j;
+    for (; j >= 0; --j)
+    {
+        vy *= y;
+        vy += bvp[DEG][j];
+    }
+    AAF v(vy);
+    --i;
+    for (; i >= 0; --i)
+    {
+        int j = DEG - i;
+        AAF vy(bvp[i][j]);
+        --j;
+        for (; j >= 0; --j)
+        {
+            vy *= y;
+            vy += bvp[i][j];
+        }
+        v *= x;
+        v += vy;
+    }
+    return v;
+
+//    return
+//      -1 +( 0.00945115 +( -4.11799e-05 +( 1.01365e-07 +( -1.35037e-10 + 7.7868e-14*y)*y)*y)*y)*y
+//        + (0.00837569 +( -6.24676e-05 +( 1.96093e-07 +( -3.09683e-10 + 1.95681e-13*y)*y)*y)*y
+//        + (-2.39448e-05 +( 1.3331e-07 +( -2.65787e-10 + 1.96698e-13*y)*y)*y
+//        + (2.76173e-08 +( -1.01069e-10 + 9.88596e-14*y)*y
+//        + (-1.43584e-11 + 2.48433e-14*y  + 2.49723e-15*x)*x)*x)*x)*x;
+}
+
+
+
+double max_modulus (SL::MVPoly2 const& p)
+{
+    double a, m = 1;
+
+    for (size_t i = 0; i < p.get_poly().size(); ++i)
+        for (double j : p)
+        {
+            a = std::abs(j);
+            if (m < a) m = a;
+        }
+    return m;
+}
+
+void poly_to_mvpoly1(SL::MVPoly1& p, Geom::Poly const& q)
+{
+    for (size_t i = 0; i < q.size(); ++i)
+    {
+        p.coefficient(i, q[i]);
+    }
+    p.normalize();
+}
+
+void make_implicit_curve (SL::MVPoly2& ic, D2<SBasis> const& pc)
+{
+    Geom::Poly pc0 = sbasis_to_poly(pc[0]);
+    Geom::Poly pc1 = sbasis_to_poly(pc[1]);
+
+//    std::cerr << "parametrization: \n";
+//    std::cerr << "pc0 = " << pc0 << std::endl;
+//    std::cerr << "pc1 = " << pc1 << "\n\n";
+
+    SL::MVPoly1 f, g;
+    poly_to_mvpoly1(f, pc0);
+    poly_to_mvpoly1(g, pc1);
+
+//    std::cerr << "parametrization: \n";
+//    std::cerr << "f = " << f << std::endl;
+//    std::cerr << "g = " << g << "\n\n";
+
+    Geom::SL::basis_type b;
+    microbasis(b, f, g);
+
+    Geom::SL::MVPoly3 p, q;
+    basis_to_poly(p, b[0]);
+    basis_to_poly(q, b[1]);
+
+//    std::cerr << "generators as polynomial in R[t,x,y] : \n";
+//    std::cerr << "p = " << p << std::endl;
+//    std::cerr << "q = " << q << "\n\n";
+
+
+    Geom::SL::Matrix<Geom::SL::MVPoly2> B = make_bezout_matrix(p, q);
+    ic = determinant_minor(B);
+    ic.normalize();
+    double m = max_modulus(ic);
+    ic /= m;
+
+//    std::cerr << "Bezout matrix: (entries are bivariate polynomials) \n";
+//    std::cerr << "B = " << B << "\n\n";
+//    std::cerr << "determinant: \n";
+//    std::cerr << "r(x, y) = " << ic << "\n\n";
+
+}
+
+//namespace Geom{ namespace SL{
+//
+//template<>
+//struct zero<AAF, false>
+//{
+//    AAF operator() () const
+//    {
+//        return AAF(0);
+//    }
+//};
+//
+//} }
+
+class ImplicitToy : public Toy
+{
+    bool contains_zero (implicit_curve_t const& eval,
+                        Rect r, double w=1e-5)
+    {
+        ++iters;
+        AAF x(interval(r.left(), r.right()));
+        AAF y(interval(r.top(), r.bottom()));
+        AAF f = eval(x, y);
+        double a = f.index_coeff(x.get_index(0)) / x.index_coeff(x.get_index(0));
+        double b = f.index_coeff(y.get_index(0)) / y.index_coeff(y.get_index(0));
+        AAF d = a*x + b*y - f;
+        interval ivl(d);
+        Point n(a,b);
+        OptRect out = tighten(r, n, Interval(ivl.min(), ivl.max()));
+        if (f.straddles_zero())
+        {
+            if ((r.width() > w) || (r.height() > w))
+            {
+                Point c = r.midpoint();
+                Rect oldr = r;
+                if (out)  r = *out;
+                // Three possibilities:
+                // 1) the trim operation buys us enough that we should just iterate
+                if (1 && (r.area() < oldr.area()*0.25))
+                {
+                    return contains_zero(eval, r,  w);
+                }
+                // 2) one dimension is significantly smaller
+                else if (1 && (r[1].extent() < oldr[1].extent()*0.5))
+                {
+                    return contains_zero (eval,
+                                          Rect(Interval(r.left(), r.right()),
+                                               Interval(r.top(), c[1])), w)
+                        || contains_zero (eval,
+                                          Rect(Interval(r.left(), r.right()),
+                                               Interval(c[1], r.bottom())), w);
+                }
+                else if (1 && (r[0].extent() < oldr[0].extent()*0.5))
+                {
+                    return contains_zero (eval,
+                                          Rect(Interval(r.left(), c[0]),
+                                               Interval(r.top(), r.bottom())), w)
+                        || contains_zero (eval,
+                                          Rect(Interval(c[0], r.right()),
+                                               Interval(r.top(), r.bottom())), w);
+                }
+                // 3) to ensure progress we must do a four way split
+                else
+                {
+                    return contains_zero (eval,
+                                          Rect(Interval(r.left(), c[0]),
+                                               Interval(r.top(), c[1])), w)
+                        || contains_zero (eval, 
+                                          Rect(Interval(c[0], r.right()),
+                                               Interval(r.top(), c[1])), w)
+                        || contains_zero (eval, 
+                                          Rect(Interval(r.left(), c[0]),
+                                               Interval(c[1], r.bottom())), w)
+                        || contains_zero (eval, 
+                                          Rect(Interval(c[0], r.right()),
+                                               Interval(c[1], r.bottom())), w);
+                }
+            }
+            //std::cout << w << " < " << r.width() << " , " << r.height() << std::endl;
+            //std::cout << r.min() << " - " << r.max() << std::endl;
+            return true;
+        }
+        return false;
+    }  // end recursive_implicit
+
+    
+    void draw_implicit_curve (cairo_t*cr, implicit_curve_t const& eval,
+                              Point const& origin, Rect r, double w)
+    {
+        ++iters;
+        AAF x(interval(r.left(), r.right()));
+        AAF y(interval(r.top(), r.bottom()));
+        //assert(x.rad() > 0);
+        //assert(y.rad() > 0);
+//        time(&t0);
+        AAF f = eval(x-origin[X], y-origin[Y]);
+//        time(&t1);
+//        d1 += std::difftime(t1, t0);
+        // pivot
+//        time(&t2);
+        double a = f.index_coeff(x.get_index(0)) / x.index_coeff(x.get_index(0));
+        double b = f.index_coeff(y.get_index(0)) / y.index_coeff(y.get_index(0));
+        AAF d = a*x + b*y - f;
+        interval ivl(d);
+        Point n(a,b);
+        OptRect out = tighten(r, n, Interval(ivl.min(), ivl.max()));
+        if (ivl.extent() < 0.5*L2(n))
+        {
+            draw_line_in_rect(cr, r, n, ivl.middle());
+            return;
+        }
+//        time(&t3);
+//        d2 += std::difftime(t3, t2);
+        if ((r.width() > w) || (r.height() > w))
+        {
+            if (f.straddles_zero())
+            {
+                Point c = r.midpoint();
+                Rect oldr = r;
+                if (out)  r = *out;
+                // Three possibilities:
+                // 1) the trim operation buys us enough that we should just iterate
+                if (1 && (r.area() < oldr.area()*0.25))
+                {
+                    draw_implicit_curve(cr, eval, origin, r,  w);
+                }
+                // 2) one dimension is significantly smaller
+                else if (1 && (r[1].extent() < oldr[1].extent()*0.5))
+                {
+                    draw_implicit_curve (cr, eval, origin,
+                                         Rect(Interval(r.left(), r.right()),
+                                              Interval(r.top(), c[1])), w);
+                    draw_implicit_curve (cr, eval, origin,
+                                         Rect(Interval(r.left(), r.right()),
+                                              Interval(c[1], r.bottom())), w);
+                }
+                else if (1 && (r[0].extent() < oldr[0].extent()*0.5))
+                {
+                    draw_implicit_curve (cr, eval, origin,
+                                         Rect(Interval(r.left(), c[0]),
+                                              Interval(r.top(), r.bottom())), w);
+                    draw_implicit_curve (cr, eval, origin,
+                                         Rect(Interval(c[0], r.right()),
+                                              Interval(r.top(), r.bottom())), w);
+                }
+                // 3) to ensure progress we must do a four way split
+                else
+                {
+                    draw_implicit_curve (cr, eval, origin,
+                                         Rect(Interval(r.left(), c[0]),
+                                              Interval(r.top(), c[1])), w);
+                    draw_implicit_curve (cr, eval, origin,
+                                         Rect(Interval(c[0], r.right()),
+                                              Interval(r.top(), c[1])), w);
+                    draw_implicit_curve (cr, eval, origin,
+                                         Rect(Interval(r.left(), c[0]),
+                                              Interval(c[1], r.bottom())), w);
+                    draw_implicit_curve (cr, eval, origin,
+                                         Rect(Interval(c[0], r.right()),
+                                              Interval(c[1], r.bottom())), w);
+                }
+            }
+        } else {
+            if(contains_zero(eval, r*Geom::Translate(-origin))) {
+                cairo_save(cr);
+                cairo_set_source_rgb(cr, 0,0.5,0);
+                cairo_rectangle(cr, r);
+                cairo_fill(cr);
+                cairo_restore(cr);
+            }
+        }
+    }  // end recursive_implicit
+
+    void draw( cairo_t *cr, std::ostringstream *notify,
+               int width, int height, bool save, std::ostringstream *timer_stream) override
+    {
+        iters = 0;
+        d1 = d2 = 0;
+        cairo_set_line_width (cr, 0.3);
+        D2<SBasis> A = pshA.asBezier();
+        cairo_d2_sb(cr, A);
+        cairo_stroke(cr);
+
+        SL::MVPoly2 ic;
+        make_implicit_curve(ic, A);
+
+        cairo_set_source_rgba (cr, 0., 0., 0, 1);
+        cairo_set_line_width (cr, 0.8);
+        draw_implicit_curve (cr, ic, orig_handle.pos,
+                             Rect(Interval(0,width), Interval(0, height)), 1);
+        cairo_stroke(cr);
+
+//        std::cerr << "D1 = " << d1 << std::endl;
+//        std::cerr << "D2 = " << d2 << std::endl;
+
+        *notify << "iter: " << iters;
+        Toy::draw(cr, notify, width, height, save,timer_stream);
+    }
+
+
+public:
+    ImplicitToy(unsigned int _A_bez_ord)
+        : A_bez_ord(_A_bez_ord)
+    {
+        handles.push_back(&orig_handle);
+        orig_handle.pos = Point(0,0); //Point(300,300);
+
+        handles.push_back(&pshA);
+        for (unsigned int i = 0; i < A_bez_ord; ++i)
+            pshA.push_back(Geom::Point(uniform()*400, uniform()*400));
+    }
+
+private:
+    unsigned int A_bez_ord;
+    PointHandle orig_handle;
+    PointSetHandle pshA;
+    time_t t0, t1, t2, t3;
+    double d1, d2;
+    unsigned int iters;
+};
+
+
+int main(int argc, char **argv)
+{
+    unsigned int A_bez_ord=5;
+    if(argc > 1)
+        sscanf(argv[1], "%d", &A_bez_ord);
+
+
+    init( argc, argv, new ImplicitToy(A_bez_ord));
+    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:fileencoding=utf-8:textwidth=99 :