remove bezier-utils. now we use 2geom's one.

remove bezier-utils-test as well.

diff --git a/src/display/CMakeLists.txt b/src/display/CMakeLists.txt
index f8dc72afe6d0b6e613c2e514bbd500ea97e9835d..a1d61de9479664c976fc6ad5503b64307a5d881a 100644 (file)
--- a/src/display/CMakeLists.txt
+++ b/src/display/CMakeLists.txt
@@ -1,6 +1,4 @@
 SET(display_SRC
-bezier-utils.cpp
-#bezier-utils-test.cpp
 canvas-arena.cpp
 canvas-axonomgrid.cpp
 canvas-bpath.cpp

diff --git a/src/display/Makefile_insert b/src/display/Makefile_insert
index e6877ccfb61dc20fa9b2de946ce7fc20da45b1bb..ea1994040bdd17f5ad1623ac84b16369fe316479 100644 (file)
--- a/src/display/Makefile_insert
+++ b/src/display/Makefile_insert
@@ -32,8 +32,6 @@ display_libspdisplay_a_SOURCES = \
        display/nr-arena-glyphs.h       \
        display/canvas-arena.cpp        \
        display/canvas-arena.h  \
-       display/bezier-utils.cpp        \
-       display/bezier-utils.h  \
        display/canvas-bpath.cpp        \
        display/canvas-bpath.h  \
        display/canvas-grid.cpp \
@@ -132,7 +130,6 @@ display_libspdisplay_a_SOURCES = \
 # ### CxxTest stuff ####
 # ######################
 ddislay_testsuites = \
-       $(srcdir)/display/bezier-utils-test.h \
        $(srcdir)/display/curve-test.h
 
 display_test_display_SOURCES = \

diff --git a/src/display/bezier-utils-test.cpp b/src/display/bezier-utils-test.cpp
deleted file mode 100644 (file)
index 4d63138..0000000
--- a/src/display/bezier-utils-test.cpp
+++ /dev/null
@@ -1,334 +0,0 @@
-#include "../utest/utest.h"
-#include <glib.h>
-#include <libnr/nr-macros.h> /* NR_DF_TEST_CLOSE */
-
-/* mental disclaims all responsibility for this evil idea for testing
-   static functions.  The main disadvantages are that we retain the
-   #define's and `using' directives of the included file. */
-#include "bezier-utils.cpp"
-
-using Geom::Point;
-
-static bool range_approx_equal(double const a[], double const b[], unsigned len);
-
-/* (Returns false if NaN encountered.) */
-template<class T>
-static bool range_equal(T const a[], T const b[], unsigned len) {
-    for (unsigned i = 0; i < len; ++i) {
-        if ( a[i] != b[i] ) {
-            return false;
-        }
-    }
-    return true;
-}
-
-inline bool point_approx_equal(Geom::Point const &a, Geom::Point const &b, double const eps)
-{
-    using Geom::X; using Geom::Y;
-    return ( NR_DF_TEST_CLOSE(a[X], b[X], eps) &&
-             NR_DF_TEST_CLOSE(a[Y], b[Y], eps) );
-}
-
-static inline double square(double const x) {
-    return x * x;
-}
-
-/** Determine whether the found control points are the same as previously found on some developer's
-    machine.  Doesn't call utest__fail, just writes a message to stdout for diagnostic purposes:
-    the most important test is that the root-mean-square of errors in the estimation are low rather
-    than that the control points found are the same.
-**/
-static void compare_ctlpts(Point const est_b[], Point const exp_est_b[])
-{
-    unsigned diff_mask = 0;
-    for (unsigned i = 0; i < 4; ++i) {
-        for (unsigned d = 0; d < 2; ++d) {
-            if ( fabs( est_b[i][d] - exp_est_b[i][d] ) > 1.1e-5 ) {
-                diff_mask |= 1 << ( i * 2 + d );
-            }
-        }
-    }
-    if ( diff_mask != 0 ) {
-        printf("Warning: got different control points from previously-coded (diffs=0x%x).\n",
-               diff_mask);
-        printf(" Previous:");
-        for (unsigned i = 0; i < 4; ++i) {
-            printf(" (%g, %g)", exp_est_b[i][0], exp_est_b[i][1]); // localizing ok
-        }
-        putchar('\n');
-        printf(" Found:   ");
-        for (unsigned i = 0; i < 4; ++i) {
-            printf(" (%g, %g)", est_b[i][0], est_b[i][1]); // localizing ok
-        }
-        putchar('\n');
-    }
-}
-
-static void compare_rms(Point const est_b[], double const t[], Point const d[], unsigned const n,
-                        double const exp_rms_error)
-{
-    double sum_errsq = 0.0;
-    for (unsigned i = 0; i < n; ++i) {
-        Point const fit_pt = bezier_pt(3, est_b, t[i]);
-        Point const diff = fit_pt - d[i];
-        sum_errsq += dot(diff, diff);
-    }
-    double const rms_error = sqrt( sum_errsq / n );
-    UTEST_ASSERT( rms_error <= exp_rms_error + 1.1e-6 );
-    if ( rms_error < exp_rms_error - 1.1e-6 ) {
-        /* The fitter code appears to have improved [or the floating point calculations differ
-           on this machine from the machine where exp_rms_error was calculated]. */
-        printf("N.B. rms_error regression requirement can be decreased: have rms_error=%g.\n", rms_error); // localizing ok
-    }
-}
-
-int main(int /*argc*/, char */*argv*/[]) {
-    utest_start("bezier-utils.cpp");
-
-    UTEST_TEST("copy_without_nans_or_adjacent_duplicates") {
-        Geom::Point const src[] = {
-            Point(2., 3.),
-            Point(2., 3.),
-            Point(0., 0.),
-            Point(2., 3.),
-            Point(2., 3.),
-            Point(1., 9.),
-            Point(1., 9.)
-        };
-        Point const exp_dest[] = {
-            Point(2., 3.),
-            Point(0., 0.),
-            Point(2., 3.),
-            Point(1., 9.)
-        };
-        g_assert( G_N_ELEMENTS(src) == 7 );
-        Point dest[7];
-        struct tst {
-            unsigned src_ix0;
-            unsigned src_len;
-            unsigned exp_dest_ix0;
-            unsigned exp_dest_len;
-        } const test_data[] = {
-            /* src start ix, src len, exp_dest start ix, exp dest len */
-            {0, 0, 0, 0},
-            {2, 1, 1, 1},
-            {0, 1, 0, 1},
-            {0, 2, 0, 1},
-            {0, 3, 0, 2},
-            {1, 3, 0, 3},
-            {0, 5, 0, 3},
-            {0, 6, 0, 4},
-            {0, 7, 0, 4}
-        };
-        for (unsigned i = 0 ; i < G_N_ELEMENTS(test_data) ; ++i) {
-            tst const &t = test_data[i];
-            UTEST_ASSERT( t.exp_dest_len
-                          == copy_without_nans_or_adjacent_duplicates(src + t.src_ix0,
-                                                                      t.src_len,
-                                                                      dest) );
-            UTEST_ASSERT(range_equal(dest,
-                                     exp_dest + t.exp_dest_ix0,
-                                     t.exp_dest_len));
-        }
-    }
-
-    UTEST_TEST("bezier_pt(1)") {
-        Point const a[] = {Point(2.0, 4.0),
-                           Point(1.0, 8.0)};
-        UTEST_ASSERT( bezier_pt(1, a, 0.0) == a[0] );
-        UTEST_ASSERT( bezier_pt(1, a, 1.0) == a[1] );
-        UTEST_ASSERT( bezier_pt(1, a, 0.5) == Point(1.5, 6.0) );
-        double const t[] = {0.5, 0.25, 0.3, 0.6};
-        for (unsigned i = 0; i < G_N_ELEMENTS(t); ++i) {
-            double const ti = t[i], si = 1.0 - ti;
-            UTEST_ASSERT( bezier_pt(1, a, ti) == si * a[0] + ti * a[1] );
-        }
-    }
-
-    UTEST_TEST("bezier_pt(2)") {
-        Point const b[] = {Point(1.0, 2.0),
-                           Point(8.0, 4.0),
-                           Point(3.0, 1.0)};
-        UTEST_ASSERT( bezier_pt(2, b, 0.0) == b[0] );
-        UTEST_ASSERT( bezier_pt(2, b, 1.0) == b[2] );
-        UTEST_ASSERT( bezier_pt(2, b, 0.5) == Point(5.0, 2.75) );
-        double const t[] = {0.5, 0.25, 0.3, 0.6};
-        for (unsigned i = 0; i < G_N_ELEMENTS(t); ++i) {
-            double const ti = t[i], si = 1.0 - ti;
-            Point const exp_pt( si*si * b[0] + 2*si*ti * b[1] + ti*ti * b[2] );
-            Point const pt(bezier_pt(2, b, ti));
-            UTEST_ASSERT(point_approx_equal(pt, exp_pt, 1e-11));
-        }
-    }
-
-    Point const c[] = {Point(1.0, 2.0),
-                       Point(8.0, 4.0),
-                       Point(3.0, 1.0),
-                       Point(-2.0, -4.0)};
-    UTEST_TEST("bezier_pt(3)") {
-        UTEST_ASSERT( bezier_pt(3, c, 0.0) == c[0] );
-        UTEST_ASSERT( bezier_pt(3, c, 1.0) == c[3] );
-        UTEST_ASSERT( bezier_pt(3, c, 0.5) == Point(4.0, 13.0/8.0) );
-        double const t[] = {0.5, 0.25, 0.3, 0.6};
-        for (unsigned i = 0; i < G_N_ELEMENTS(t); ++i) {
-            double const ti = t[i], si = 1.0 - ti;
-            UTEST_ASSERT( LInfty( bezier_pt(3, c, ti)
-                                  - ( si*si*si * c[0] +
-                                      3*si*si*ti * c[1] +
-                                      3*si*ti*ti * c[2] +
-                                      ti*ti*ti * c[3] ) )
-                          < 1e-4 );
-        }
-    }
-
-    struct Err_tst {
-        Point pt;
-        double u;
-        double err;
-    } const err_tst[] = {
-        {c[0], 0.0, 0.0},
-        {Point(4.0, 13.0/8.0), 0.5, 0.0},
-        {Point(4.0, 2.0), 0.5, 9.0/64.0},
-        {Point(3.0, 2.0), 0.5, 1.0 + 9.0/64.0},
-        {Point(6.0, 2.0), 0.5, 4.0 + 9.0/64.0},
-        {c[3], 1.0, 0.0},
-    };
-
-    UTEST_TEST("compute_max_error_ratio") {
-        Point d[G_N_ELEMENTS(err_tst)];
-        double u[G_N_ELEMENTS(err_tst)];
-        for (unsigned i = 0; i < G_N_ELEMENTS(err_tst); ++i) {
-            Err_tst const &t = err_tst[i];
-            d[i] = t.pt;
-            u[i] = t.u;
-        }
-        g_assert( G_N_ELEMENTS(u) == G_N_ELEMENTS(d) );
-        unsigned max_ix = ~0u;
-        double const err_ratio = compute_max_error_ratio(d, u, G_N_ELEMENTS(d), c, 1.0, &max_ix);
-        UTEST_ASSERT( fabs( sqrt(err_tst[4].err) - err_ratio ) < 1e-12 );
-        UTEST_ASSERT( max_ix == 4 );
-    }
-
-    UTEST_TEST("chord_length_parameterize") {
-        /* n == 2 */
-        {
-            Point const d[] = {Point(2.9415, -5.8149),
-                               Point(23.021, 4.9814)};
-            double u[G_N_ELEMENTS(d)];
-            double const exp_u[] = {0.0, 1.0};
-            g_assert( G_N_ELEMENTS(u) == G_N_ELEMENTS(exp_u) );
-            chord_length_parameterize(d, u, G_N_ELEMENTS(d));
-            UTEST_ASSERT(range_equal(u, exp_u, G_N_ELEMENTS(exp_u)));
-        }
-
-        /* Straight line. */
-        {
-            double const exp_u[] = {0.0, 0.1829, 0.2105, 0.2105, 0.619, 0.815, 0.999, 1.0};
-            unsigned const n = G_N_ELEMENTS(exp_u);
-            Point d[n];
-            double u[n];
-            Point const a(-23.985, 4.915), b(4.9127, 5.203);
-            for (unsigned i = 0; i < n; ++i) {
-                double bi = exp_u[i], ai = 1.0 - bi;
-                d[i] = ai * a  +  bi * b;
-            }
-            chord_length_parameterize(d, u, n);
-            UTEST_ASSERT(range_approx_equal(u, exp_u, n));
-        }
-    }
-
-    /* Feed it some points that can be fit exactly with a single bezier segment, and see how
-       well it manages. */
-    Point const src_b[4] = {Point(5., -3.),
-                            Point(8., 0.),
-                            Point(4., 2.),
-                            Point(3., 3.)};
-    double const t[] = {0.0, .001, .03, .05, .09, .13, .18, .25, .29, .33, .39, .44,
-                        .51, .57, .62, .69, .75, .81, .91, .93, .97, .98, .999, 1.0};
-    unsigned const n = G_N_ELEMENTS(t);
-    Point d[n];
-    for (unsigned i = 0; i < n; ++i) {
-        d[i] = bezier_pt(3, src_b, t[i]);
-    }
-    Point const tHat1(unit_vector( src_b[1] - src_b[0] ));
-    Point const tHat2(unit_vector( src_b[2] - src_b[3] ));
-
-    UTEST_TEST("generate_bezier") {
-        Point est_b[4];
-        generate_bezier(est_b, d, t, n, tHat1, tHat2, 1.0);
-
-        compare_ctlpts(est_b, src_b);
-
-        /* We're being unfair here in using our t[] rather than best t[] for est_b: we
-           may over-estimate RMS of errors. */
-        compare_rms(est_b, t, d, n, 1e-8);
-    }
-
-    UTEST_TEST("sp_bezier_fit_cubic_full") {
-        Point est_b[4];
-        int splitpoints[2];
-        gint const succ = sp_bezier_fit_cubic_full(est_b, splitpoints, d, n, tHat1, tHat2, square(1.2), 1);
-        UTEST_ASSERT( succ == 1 );
-
-        Point const exp_est_b[4] = {
-            Point(5.000000, -3.000000),
-            Point(7.5753, -0.4247),
-            Point(4.77533, 1.22467),
-            Point(3, 3)
-        };
-        compare_ctlpts(est_b, exp_est_b);
-
-        /* We're being unfair here in using our t[] rather than best t[] for est_b: we
-           may over-estimate RMS of errors. */
-        compare_rms(est_b, t, d, n, .307911);
-    }
-
-    UTEST_TEST("sp_bezier_fit_cubic") {
-        Point est_b[4];
-        gint const succ = sp_bezier_fit_cubic(est_b, d, n, square(1.2));
-        UTEST_ASSERT( succ == 1 );
-
-        Point const exp_est_b[4] = {
-            Point(5.000000, -3.000000),
-            Point(7.57134, -0.423509),
-            Point(4.77929, 1.22426),
-            Point(3, 3)
-        };
-        compare_ctlpts(est_b, exp_est_b);
-
-#if 1 /* A change has been made to right_tangent.  I believe that usually this change
-         will result in better fitting, but it won't do as well for this example where
-         we happen to be feeding a t=0.999 point to the fitter. */
-        printf("TODO: Update this test case for revised right_tangent implementation.\n");
-        /* In particular, have a test case to show whether the new implementation
-           really is likely to be better on average. */
-#else
-        /* We're being unfair here in using our t[] rather than best t[] for est_b: we
-           may over-estimate RMS of errors. */
-        compare_rms(est_b, t, d, n, .307983);
-#endif
-    }
-
-    return !utest_end();
-}
-
-/* (Returns false if NaN encountered.) */
-static bool range_approx_equal(double const a[], double const b[], unsigned const len) {
-    for (unsigned i = 0; i < len; ++i) {
-        if (!( fabs( a[i] - b[i] ) < 1e-4 )) {
-            return false;
-        }
-    }
-    return true;
-}
-
-/*
-  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 :

diff --git a/src/display/bezier-utils-test.h b/src/display/bezier-utils-test.h
deleted file mode 100644 (file)
index 3be6435..0000000
--- a/src/display/bezier-utils-test.h
+++ /dev/null
@@ -1,349 +0,0 @@
-#include <cxxtest/TestSuite.h>
-
-#include <glib.h>
-#include <libnr/nr-macros.h> /* NR_DF_TEST_CLOSE */
-#include <sstream>
-
-/* mental disclaims all responsibility for this evil idea for testing
-   static functions.  The main disadvantages are that we retain the
-   #define's and `using' directives of the included file. */
-#include "bezier-utils.cpp"
-
-/* (Returns false if NaN encountered.) */
-static bool range_approx_equal(double const a[], double const b[], unsigned const len) {
-    for (unsigned i = 0; i < len; ++i) {
-        if (!( fabs( a[i] - b[i] ) < 1e-4 )) {
-            return false;
-        }
-    }
-    return true;
-}
-
-static inline bool point_approx_equal(Geom::Point const &a, Geom::Point const &b, double const eps)
-{
-    return ( NR_DF_TEST_CLOSE(a[Geom::X], b[Geom::X], eps) &&
-             NR_DF_TEST_CLOSE(a[Geom::Y], b[Geom::Y], eps) );
-}
-
-static inline double square(double const x) {
-    return x * x;
-}
-
-/** Determine whether the found control points are the same as previously found on some developer's
-    machine.  Doesn't call utest__fail, just writes a message to stdout for diagnostic purposes:
-    the most important test is that the root-mean-square of errors in the estimation are low rather
-    than that the control points found are the same.
-**/
-static void compare_ctlpts(Geom::Point const est_b[], Geom::Point const exp_est_b[])
-{
-    unsigned diff_mask = 0;
-    for (unsigned i = 0; i < 4; ++i) {
-        for (unsigned d = 0; d < 2; ++d) {
-            if ( fabs( est_b[i][d] - exp_est_b[i][d] ) > 1.1e-5 ) {
-                diff_mask |= 1 << ( i * 2 + d );
-            }
-        }
-    }
-    if ( diff_mask != 0 ) {
-        std::stringstream msg;
-        msg << "Got different control points from previously-coded (diffs=0x" << std::hex << diff_mask << "\n";
-        msg << " Previous:";
-        for (unsigned i = 0; i < 4; ++i) {
-            msg << " (" << exp_est_b[i][0] << ", " << exp_est_b[i][1] << ")"; // localizing ok
-        }
-        msg << "\n";
-        msg << " Found:   ";
-        for (unsigned i = 0; i < 4; ++i) {
-            msg << " (" << est_b[i][0] << ", " << est_b[i][1] << ")"; // localizing ok
-        }
-        msg << "\n";
-        TS_WARN(msg.str().c_str());
-    }
-}
-
-static void compare_rms(Geom::Point const est_b[], double const t[], Geom::Point const d[], unsigned const n,
-                        double const exp_rms_error)
-{
-    double sum_errsq = 0.0;
-    for (unsigned i = 0; i < n; ++i) {
-        Geom::Point const fit_pt = bezier_pt(3, est_b, t[i]);
-        Geom::Point const diff = fit_pt - d[i];
-        sum_errsq += dot(diff, diff);
-    }
-    double const rms_error = sqrt( sum_errsq / n );
-    TS_ASSERT_LESS_THAN_EQUALS( rms_error , exp_rms_error + 1.1e-6 );
-    if ( rms_error < exp_rms_error - 1.1e-6 ) {
-        /* The fitter code appears to have improved [or the floating point calculations differ
-           on this machine from the machine where exp_rms_error was calculated]. */
-        char msg[200];
-        sprintf(msg, "N.B. rms_error regression requirement can be decreased: have rms_error=%g.", rms_error); // localizing ok
-        TS_TRACE(msg);
-    }
-}
-
-class BezierUtilsTest : public CxxTest::TestSuite {
-public:
-    static Geom::Point const c[4];
-    static double const t[24];
-    static unsigned const n;
-    Geom::Point d[24];
-    static Geom::Point const src_b[4];
-    static Geom::Point const tHat1;
-    static Geom::Point const tHat2;
-
-    BezierUtilsTest()
-    {
-        /* Feed it some points that can be fit exactly with a single bezier segment, and see how
-        well it manages. */
-        for (unsigned i = 0; i < n; ++i) {
-            d[i] = bezier_pt(3, src_b, t[i]);
-        }
-    }
-    virtual ~BezierUtilsTest() {}
-
-// createSuite and destroySuite get us per-suite setup and teardown
-// without us having to worry about static initialization order, etc.
-    static BezierUtilsTest *createSuite() { return new BezierUtilsTest(); }
-    static void destroySuite( BezierUtilsTest *suite ) { delete suite; }
-
-    void testCopyWithoutNansOrAdjacentDuplicates()
-    {
-        Geom::Point const src[] = {
-            Geom::Point(2., 3.),
-            Geom::Point(2., 3.),
-            Geom::Point(0., 0.),
-            Geom::Point(2., 3.),
-            Geom::Point(2., 3.),
-            Geom::Point(1., 9.),
-            Geom::Point(1., 9.)
-        };
-        Geom::Point const exp_dest[] = {
-            Geom::Point(2., 3.),
-            Geom::Point(0., 0.),
-            Geom::Point(2., 3.),
-            Geom::Point(1., 9.)
-        };
-        g_assert( G_N_ELEMENTS(src) == 7 );
-        Geom::Point dest[7];
-        struct tst {
-            unsigned src_ix0;
-            unsigned src_len;
-            unsigned exp_dest_ix0;
-            unsigned exp_dest_len;
-        } const test_data[] = {
-            /* src start ix, src len, exp_dest start ix, exp dest len */
-            {0, 0, 0, 0},
-            {2, 1, 1, 1},
-            {0, 1, 0, 1},
-            {0, 2, 0, 1},
-            {0, 3, 0, 2},
-            {1, 3, 0, 3},
-            {0, 5, 0, 3},
-            {0, 6, 0, 4},
-            {0, 7, 0, 4}
-        };
-        for (unsigned i = 0 ; i < G_N_ELEMENTS(test_data) ; ++i) {
-            tst const &t = test_data[i];
-            TS_ASSERT_EQUALS( t.exp_dest_len,
-                              copy_without_nans_or_adjacent_duplicates(src + t.src_ix0,
-                                                                       t.src_len,
-                                                                       dest) );
-            TS_ASSERT_SAME_DATA(dest,
-                                exp_dest + t.exp_dest_ix0,
-                                t.exp_dest_len);
-        }
-    }
-
-    void testBezierPt1()
-    {
-        Geom::Point const a[] = {Geom::Point(2.0, 4.0),
-                               Geom::Point(1.0, 8.0)};
-        TS_ASSERT_EQUALS( bezier_pt(1, a, 0.0) , a[0] );
-        TS_ASSERT_EQUALS( bezier_pt(1, a, 1.0) , a[1] );
-        TS_ASSERT_EQUALS( bezier_pt(1, a, 0.5) , Geom::Point(1.5, 6.0) );
-        double const t[] = {0.5, 0.25, 0.3, 0.6};
-        for (unsigned i = 0; i < G_N_ELEMENTS(t); ++i) {
-            double const ti = t[i], si = 1.0 - ti;
-            TS_ASSERT_EQUALS( bezier_pt(1, a, ti) , si * a[0] + ti * a[1] );
-        }
-    }
-
-    void testBezierPt2()
-    {
-        Geom::Point const b[] = {Geom::Point(1.0, 2.0),
-                               Geom::Point(8.0, 4.0),
-                               Geom::Point(3.0, 1.0)};
-        TS_ASSERT_EQUALS( bezier_pt(2, b, 0.0) , b[0] );
-        TS_ASSERT_EQUALS( bezier_pt(2, b, 1.0) , b[2] );
-        TS_ASSERT_EQUALS( bezier_pt(2, b, 0.5) , Geom::Point(5.0, 2.75) );
-        double const t[] = {0.5, 0.25, 0.3, 0.6};
-        for (unsigned i = 0; i < G_N_ELEMENTS(t); ++i) {
-            double const ti = t[i], si = 1.0 - ti;
-            Geom::Point const exp_pt( si*si * b[0] + 2*si*ti * b[1] + ti*ti * b[2] );
-            Geom::Point const pt(bezier_pt(2, b, ti));
-            TS_ASSERT(point_approx_equal(pt, exp_pt, 1e-11));
-        }
-    }
-
-    void testBezierPt3()
-    {
-        TS_ASSERT_EQUALS( bezier_pt(3, c, 0.0) , c[0] );
-        TS_ASSERT_EQUALS( bezier_pt(3, c, 1.0) , c[3] );
-        TS_ASSERT_EQUALS( bezier_pt(3, c, 0.5) , Geom::Point(4.0, 13.0/8.0) );
-        double const t[] = {0.5, 0.25, 0.3, 0.6};
-        for (unsigned i = 0; i < G_N_ELEMENTS(t); ++i) {
-            double const ti = t[i], si = 1.0 - ti;
-            TS_ASSERT( LInfty( bezier_pt(3, c, ti)
-                                  - ( si*si*si * c[0] +
-                                      3*si*si*ti * c[1] +
-                                      3*si*ti*ti * c[2] +
-                                      ti*ti*ti * c[3] ) )
-                          < 1e-4 );
-        }
-    }
-
-    void testComputeMaxErrorRatio()
-    {
-        struct Err_tst {
-            Geom::Point pt;
-            double u;
-            double err;
-        } const err_tst[] = {
-            {c[0], 0.0, 0.0},
-            {Geom::Point(4.0, 13.0/8.0), 0.5, 0.0},
-            {Geom::Point(4.0, 2.0), 0.5, 9.0/64.0},
-            {Geom::Point(3.0, 2.0), 0.5, 1.0 + 9.0/64.0},
-            {Geom::Point(6.0, 2.0), 0.5, 4.0 + 9.0/64.0},
-            {c[3], 1.0, 0.0},
-        };
-        Geom::Point d[G_N_ELEMENTS(err_tst)];
-        double u[G_N_ELEMENTS(err_tst)];
-        for (unsigned i = 0; i < G_N_ELEMENTS(err_tst); ++i) {
-            Err_tst const &t = err_tst[i];
-            d[i] = t.pt;
-            u[i] = t.u;
-        }
-        g_assert( G_N_ELEMENTS(u) == G_N_ELEMENTS(d) );
-        unsigned max_ix = ~0u;
-        double const err_ratio = compute_max_error_ratio(d, u, G_N_ELEMENTS(d), c, 1.0, &max_ix);
-        TS_ASSERT_LESS_THAN( fabs( sqrt(err_tst[4].err) - err_ratio ) , 1e-12 );
-        TS_ASSERT_EQUALS( max_ix , 4u );
-    }
-
-    void testChordLengthParameterize()
-    {
-        /* n == 2 */
-        {
-            Geom::Point const d[] = {Geom::Point(2.9415, -5.8149),
-                                   Geom::Point(23.021, 4.9814)};
-            double u[G_N_ELEMENTS(d)];
-            double const exp_u[] = {0.0, 1.0};
-            g_assert( G_N_ELEMENTS(u) == G_N_ELEMENTS(exp_u) );
-            chord_length_parameterize(d, u, G_N_ELEMENTS(d));
-            TS_ASSERT_SAME_DATA(u, exp_u, G_N_ELEMENTS(exp_u));
-        }
-
-        /* Straight line. */
-        {
-            double const exp_u[] = {0.0, 0.1829, 0.2105, 0.2105, 0.619, 0.815, 0.999, 1.0};
-            unsigned const n = G_N_ELEMENTS(exp_u);
-            Geom::Point d[n];
-            double u[n];
-            Geom::Point const a(-23.985, 4.915), b(4.9127, 5.203);
-            for (unsigned i = 0; i < n; ++i) {
-                double bi = exp_u[i], ai = 1.0 - bi;
-                d[i] = ai * a  +  bi * b;
-            }
-            chord_length_parameterize(d, u, n);
-            TS_ASSERT(range_approx_equal(u, exp_u, n));
-        }
-    }
-
-    void testGenerateBezier()
-    {
-        Geom::Point est_b[4];
-        generate_bezier(est_b, d, t, n, tHat1, tHat2, 1.0);
-
-        compare_ctlpts(est_b, src_b);
-
-        /* We're being unfair here in using our t[] rather than best t[] for est_b: we
-           may over-estimate RMS of errors. */
-        compare_rms(est_b, t, d, n, 1e-8);
-    }
-
-    void testSpBezierFitCubicFull()
-    {
-        Geom::Point est_b[4];
-        int splitpoints[2];
-        gint const succ = sp_bezier_fit_cubic_full(est_b, splitpoints, d, n, tHat1, tHat2, square(1.2), 1);
-        TS_ASSERT_EQUALS( succ , 1 );
-
-        Geom::Point const exp_est_b[4] = {
-            Geom::Point(5.000000, -3.000000),
-            Geom::Point(7.5753, -0.4247),
-            Geom::Point(4.77533, 1.22467),
-            Geom::Point(3, 3)
-        };
-        compare_ctlpts(est_b, exp_est_b);
-
-        /* We're being unfair here in using our t[] rather than best t[] for est_b: we
-           may over-estimate RMS of errors. */
-        compare_rms(est_b, t, d, n, .307911);
-    }
-
-    void testSpBezierFitCubic()
-    {
-        Geom::Point est_b[4];
-        gint const succ = sp_bezier_fit_cubic(est_b, d, n, square(1.2));
-        TS_ASSERT_EQUALS( succ , 1 );
-
-        Geom::Point const exp_est_b[4] = {
-            Geom::Point(5.000000, -3.000000),
-            Geom::Point(7.57134, -0.423509),
-            Geom::Point(4.77929, 1.22426),
-            Geom::Point(3, 3)
-        };
-        compare_ctlpts(est_b, exp_est_b);
-
-#if 1 /* A change has been made to right_tangent.  I believe that usually this change
-         will result in better fitting, but it won't do as well for this example where
-         we happen to be feeding a t=0.999 point to the fitter. */
-        TS_WARN("TODO: Update this test case for revised right_tangent implementation.");
-        /* In particular, have a test case to show whether the new implementation
-           really is likely to be better on average. */
-#else
-        /* We're being unfair here in using our t[] rather than best t[] for est_b: we
-           may over-estimate RMS of errors. */
-        compare_rms(est_b, t, d, n, .307983);
-#endif
-    }
-};
-
-// This is not very neat, but since we know this header is only included by the generated CxxTest file it shouldn't give any problems
-Geom::Point const BezierUtilsTest::c[4] = {
-    Geom::Point(1.0, 2.0),
-    Geom::Point(8.0, 4.0),
-    Geom::Point(3.0, 1.0),
-    Geom::Point(-2.0, -4.0)};
-double const BezierUtilsTest::t[24] = {
-    0.0, .001, .03, .05, .09, .13, .18, .25, .29, .33,  .39, .44,
-    .51,  .57, .62, .69, .75, .81, .91, .93, .97, .98, .999, 1.0};
-unsigned const BezierUtilsTest::n = G_N_ELEMENTS(BezierUtilsTest::t);
-Geom::Point const BezierUtilsTest::src_b[4] = {
-    Geom::Point(5., -3.),
-    Geom::Point(8., 0.),
-    Geom::Point(4., 2.),
-    Geom::Point(3., 3.)};
-Geom::Point const BezierUtilsTest::tHat1(unit_vector( BezierUtilsTest::src_b[1] - BezierUtilsTest::src_b[0] ));
-Geom::Point const BezierUtilsTest::tHat2(unit_vector( BezierUtilsTest::src_b[2] - BezierUtilsTest::src_b[3] ));
-
-/*
-  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 :

diff --git a/src/display/bezier-utils-work.txt b/src/display/bezier-utils-work.txt
deleted file mode 100644 (file)
index 75ec5be..0000000
--- a/src/display/bezier-utils-work.txt
+++ /dev/null
@@ -1,34 +0,0 @@
-min .5 * sum_i lensq(bez_pt(b, u[i]) - d[i])
-
-lensq(d)=dot(d, d) = d.x * d.x + d.y * d.y
-
-sum_i (f(i) + g(i)) = sum_i f(i) + sum_i g(i), so
-we can separate into x,y parts.  Since they are the same, we write `z' in the below
-to mean either x or y.
-
-.5 * sum_i (bez_pt(b, u[i]) - d[i]).z^2
-
-= .5 * sum_i (B0(u[i]) * b[0] +
-              B1(u[i]) * b[1] +
-              B2(u[i]) * b[2] +
-              B3(u[i]) * b[3]
-              - d[i]           ).z^2
-
-= H.
-
-Suppose that b[0,1,3] are fixed (with b[1] perhaps being calculated
-from a prior call to existing generate_bezier).
-
-d H / d b[2].z = sum_i B2(u[i]) * (bez_pt(b, u[i]) - d[i]).z
-
-Solve for dH/db[2].z==0:
-
--sum_i B2(u[i]) B2(u[i]) * b[2].z = sum_i B2(u[i]) * (B0(u[i]) * b[0] +
-                                                      B1(u[i]) * b[1] +
-                                                      B3(u[i]) * b[3]
-                                                      - d[i]           ).z
-b[2].z = ((sum_i B2(u[i]) * (B0(u[i]) * b[0] +
-                            B1(u[i]) * b[1] +
-                            B3(u[i]) * b[3]
-                            - d[i]           ).z)
-          / -sum_i (B2(u[i]))^2)

diff --git a/src/display/bezier-utils.cpp b/src/display/bezier-utils.cpp
deleted file mode 100644 (file)
index b538f7f..0000000
--- a/src/display/bezier-utils.cpp
+++ /dev/null
@@ -1,984 +0,0 @@
-#define __SP_BEZIER_UTILS_C__
-
-/** \file
- * Bezier interpolation for inkscape drawing code.
- */
-/*
- * Original code published in:
- *   An Algorithm for Automatically Fitting Digitized Curves
- *   by Philip J. Schneider
- *  "Graphics Gems", Academic Press, 1990
- *
- * Authors:
- *   Philip J. Schneider
- *   Lauris Kaplinski <lauris@kaplinski.com>
- *   Peter Moulder <pmoulder@mail.csse.monash.edu.au>
- *
- * Copyright (C) 1990 Philip J. Schneider
- * Copyright (C) 2001 Lauris Kaplinski
- * Copyright (C) 2001 Ximian, Inc.
- * Copyright (C) 2003,2004 Monash University
- *
- * Released under GNU GPL, read the file 'COPYING' for more information
- */
-
-#define SP_HUGE 1e5
-#define noBEZIER_DEBUG
-
-#ifdef HAVE_CONFIG_H
-# include <config.h>
-#endif
-
-#ifdef HAVE_IEEEFP_H
-# include <ieeefp.h>
-#endif
-
-#include <glib.h> // g_assert()
-#include <glib/gmessages.h>
-#include <glib/gmem.h>
-#include "bezier-utils.h"
-#include <libnr/nr-point-fns.h>
-
-#include "2geom/isnan.h"
-
-
-typedef Geom::Point BezierCurve[];
-
-/* Forward declarations */
-static void generate_bezier(Geom::Point b[], Geom::Point const d[], gdouble const u[], unsigned len,
-                            Geom::Point const &tHat1, Geom::Point const &tHat2, double tolerance_sq);
-static void estimate_lengths(Geom::Point bezier[],
-                             Geom::Point const data[], gdouble const u[], unsigned len,
-                             Geom::Point const &tHat1, Geom::Point const &tHat2);
-static void estimate_bi(Geom::Point b[4], unsigned ei,
-                        Geom::Point const data[], double const u[], unsigned len);
-static void reparameterize(Geom::Point const d[], unsigned len, double u[], BezierCurve const bezCurve);
-static gdouble NewtonRaphsonRootFind(BezierCurve const Q, Geom::Point const &P, gdouble u);
-static Geom::Point sp_darray_center_tangent(Geom::Point const d[], unsigned center, unsigned length);
-static Geom::Point sp_darray_right_tangent(Geom::Point const d[], unsigned const len);
-static unsigned copy_without_nans_or_adjacent_duplicates(Geom::Point const src[], unsigned src_len, Geom::Point dest[]);
-static void chord_length_parameterize(Geom::Point const d[], gdouble u[], unsigned len);
-static double compute_max_error_ratio(Geom::Point const d[], double const u[], unsigned len,
-                                      BezierCurve const bezCurve, double tolerance,
-                                      unsigned *splitPoint);
-static double compute_hook(Geom::Point const &a, Geom::Point const &b, double const u, BezierCurve const bezCurve,
-                           double const tolerance);
-
-
-static Geom::Point const unconstrained_tangent(0, 0);
-
-
-/*
- *  B0, B1, B2, B3 : Bezier multipliers
- */
-
-#define B0(u) ( ( 1.0 - u )  *  ( 1.0 - u )  *  ( 1.0 - u ) )
-#define B1(u) ( 3 * u  *  ( 1.0 - u )  *  ( 1.0 - u ) )
-#define B2(u) ( 3 * u * u  *  ( 1.0 - u ) )
-#define B3(u) ( u * u * u )
-
-#ifdef BEZIER_DEBUG
-# define DOUBLE_ASSERT(x) g_assert( ( (x) > -SP_HUGE ) && ( (x) < SP_HUGE ) )
-# define BEZIER_ASSERT(b) do { \
-           DOUBLE_ASSERT((b)[0][Geom::X]); DOUBLE_ASSERT((b)[0][Geom::Y]);  \
-           DOUBLE_ASSERT((b)[1][Geom::X]); DOUBLE_ASSERT((b)[1][Geom::Y]);  \
-           DOUBLE_ASSERT((b)[2][Geom::X]); DOUBLE_ASSERT((b)[2][Geom::Y]);  \
-           DOUBLE_ASSERT((b)[3][Geom::X]); DOUBLE_ASSERT((b)[3][Geom::Y]);  \
-         } while(0)
-#else
-# define DOUBLE_ASSERT(x) do { } while(0)
-# define BEZIER_ASSERT(b) do { } while(0)
-#endif
-
-
-/**
- * Fit a single-segment Bezier curve to a set of digitized points.
- *
- * \return Number of segments generated, or -1 on error.
- */
-gint
-sp_bezier_fit_cubic(Geom::Point *bezier, Geom::Point const *data, gint len, gdouble error)
-{
-    return sp_bezier_fit_cubic_r(bezier, data, len, error, 1);
-}
-
-/**
- * Fit a multi-segment Bezier curve to a set of digitized points, with
- * possible weedout of identical points and NaNs.
- *
- * \param max_beziers Maximum number of generated segments
- * \param Result array, must be large enough for n. segments * 4 elements.
- *
- * \return Number of segments generated, or -1 on error.
- */
-gint
-sp_bezier_fit_cubic_r(Geom::Point bezier[], Geom::Point const data[], gint const len, gdouble const error, unsigned const max_beziers)
-{
-    g_return_val_if_fail(bezier != NULL, -1);
-    g_return_val_if_fail(data != NULL, -1);
-    g_return_val_if_fail(len > 0, -1);
-    g_return_val_if_fail(max_beziers < (1ul << (31 - 2 - 1 - 3)), -1);
-
-    Geom::Point *uniqued_data = g_new(Geom::Point, len);
-    unsigned uniqued_len = copy_without_nans_or_adjacent_duplicates(data, len, uniqued_data);
-
-    if ( uniqued_len < 2 ) {
-        g_free(uniqued_data);
-        return 0;
-    }
-
-    /* Call fit-cubic function with recursion. */
-    gint const ret = sp_bezier_fit_cubic_full(bezier, NULL, uniqued_data, uniqued_len,
-                                              unconstrained_tangent, unconstrained_tangent,
-                                              error, max_beziers);
-    g_free(uniqued_data);
-    return ret;
-}
-
-/** 
- * Copy points from src to dest, filter out points containing NaN and
- * adjacent points with equal x and y.
- * \return length of dest
- */
-static unsigned
-copy_without_nans_or_adjacent_duplicates(Geom::Point const src[], unsigned src_len, Geom::Point dest[])
-{
-    unsigned si = 0;
-    for (;;) {
-        if ( si == src_len ) {
-            return 0;
-        }
-        if (!IS_NAN(src[si][Geom::X]) &&
-            !IS_NAN(src[si][Geom::Y])) {
-            dest[0] = Geom::Point(src[si]);
-            ++si;
-            break;
-        }
-        si ++;
-    }
-    unsigned di = 0;
-    for (; si < src_len; ++si) {
-        Geom::Point const src_pt = Geom::Point(src[si]);
-        if ( src_pt != dest[di]
-             && !IS_NAN(src_pt[Geom::X])
-             && !IS_NAN(src_pt[Geom::Y])) {
-            dest[++di] = src_pt;
-        }
-    }
-    unsigned dest_len = di + 1;
-    g_assert( dest_len <= src_len );
-    return dest_len;
-}
-
-/**
- * Fit a multi-segment Bezier curve to a set of digitized points, without
- * possible weedout of identical points and NaNs.
- * 
- * \pre data is uniqued, i.e. not exist i: data[i] == data[i + 1].
- * \param max_beziers Maximum number of generated segments
- * \param Result array, must be large enough for n. segments * 4 elements.
- */
-gint
-sp_bezier_fit_cubic_full(Geom::Point bezier[], int split_points[],
-                         Geom::Point const data[], gint const len,
-                         Geom::Point const &tHat1, Geom::Point const &tHat2,
-                         double const error, unsigned const max_beziers)
-{
-    int const maxIterations = 4;   /* Max times to try iterating */
-
-    g_return_val_if_fail(bezier != NULL, -1);
-    g_return_val_if_fail(data != NULL, -1);
-    g_return_val_if_fail(len > 0, -1);
-    g_return_val_if_fail(max_beziers >= 1, -1);
-    g_return_val_if_fail(error >= 0.0, -1);
-
-    if ( len < 2 ) return 0;
-
-    if ( len == 2 ) {
-        /* We have 2 points, which can be fitted trivially. */
-        bezier[0] = data[0];
-        bezier[3] = data[len - 1];
-        double const dist = ( L2( data[len - 1]
-                                  - data[0] )
-                              / 3.0 );
-        if (IS_NAN(dist)) {
-            /* Numerical problem, fall back to straight line segment. */
-            bezier[1] = bezier[0];
-            bezier[2] = bezier[3];
-        } else {
-            bezier[1] = ( is_zero(tHat1)
-                          ? ( 2 * bezier[0] + bezier[3] ) / 3.
-                          : bezier[0] + dist * tHat1 );
-            bezier[2] = ( is_zero(tHat2)
-                          ? ( bezier[0] + 2 * bezier[3] ) / 3.
-                          : bezier[3] + dist * tHat2 );
-        }
-        BEZIER_ASSERT(bezier);
-        return 1;
-    }
-
-    /*  Parameterize points, and attempt to fit curve */
-    unsigned splitPoint;   /* Point to split point set at. */
-    bool is_corner;
-    {
-        double *u = g_new(double, len);
-        chord_length_parameterize(data, u, len);
-        if ( u[len - 1] == 0.0 ) {
-            /* Zero-length path: every point in data[] is the same.
-             *
-             * (Clients aren't allowed to pass such data; handling the case is defensive
-             * programming.)
-             */
-            g_free(u);
-            return 0;
-        }
-
-        generate_bezier(bezier, data, u, len, tHat1, tHat2, error);
-        reparameterize(data, len, u, bezier);
-
-        /* Find max deviation of points to fitted curve. */
-        double const tolerance = sqrt(error + 1e-9);
-        double maxErrorRatio = compute_max_error_ratio(data, u, len, bezier, tolerance, &splitPoint);
-
-        if ( fabs(maxErrorRatio) <= 1.0 ) {
-            BEZIER_ASSERT(bezier);
-            g_free(u);
-            return 1;
-        }
-
-        /* If error not too large, then try some reparameterization and iteration. */
-        if ( 0.0 <= maxErrorRatio && maxErrorRatio <= 3.0 ) {
-            for (int i = 0; i < maxIterations; i++) {
-                generate_bezier(bezier, data, u, len, tHat1, tHat2, error);
-                reparameterize(data, len, u, bezier);
-                maxErrorRatio = compute_max_error_ratio(data, u, len, bezier, tolerance, &splitPoint);
-                if ( fabs(maxErrorRatio) <= 1.0 ) {
-                    BEZIER_ASSERT(bezier);
-                    g_free(u);
-                    return 1;
-                }
-            }
-        }
-        g_free(u);
-        is_corner = (maxErrorRatio < 0);
-    }
-
-    if (is_corner) {
-        g_assert(splitPoint < unsigned(len));
-        if (splitPoint == 0) {
-            if (is_zero(tHat1)) {
-                /* Got spike even with unconstrained initial tangent. */
-                ++splitPoint;
-            } else {
-                return sp_bezier_fit_cubic_full(bezier, split_points, data, len, unconstrained_tangent, tHat2,
-                                                error, max_beziers);
-            }
-        } else if (splitPoint == unsigned(len - 1)) {
-            if (is_zero(tHat2)) {
-                /* Got spike even with unconstrained final tangent. */
-                --splitPoint;
-            } else {
-                return sp_bezier_fit_cubic_full(bezier, split_points, data, len, tHat1, unconstrained_tangent,
-                                                error, max_beziers);
-            }
-        }
-    }
-
-    if ( 1 < max_beziers ) {
-        /*
-         *  Fitting failed -- split at max error point and fit recursively
-         */
-        unsigned const rec_max_beziers1 = max_beziers - 1;
-
-        Geom::Point recTHat2, recTHat1;
-        if (is_corner) {
-            g_return_val_if_fail(0 < splitPoint && splitPoint < unsigned(len - 1), -1);
-            recTHat1 = recTHat2 = unconstrained_tangent;
-        } else {
-            /* Unit tangent vector at splitPoint. */
-            recTHat2 = sp_darray_center_tangent(data, splitPoint, len);
-            recTHat1 = -recTHat2;
-        }
-        gint const nsegs1 = sp_bezier_fit_cubic_full(bezier, split_points, data, splitPoint + 1,
-                                                     tHat1, recTHat2, error, rec_max_beziers1);
-        if ( nsegs1 < 0 ) {
-#ifdef BEZIER_DEBUG
-            g_print("fit_cubic[1]: recursive call failed\n");
-#endif
-            return -1;
-        }
-        g_assert( nsegs1 != 0 );
-        if (split_points != NULL) {
-            split_points[nsegs1 - 1] = splitPoint;
-        }
-        unsigned const rec_max_beziers2 = max_beziers - nsegs1;
-        gint const nsegs2 = sp_bezier_fit_cubic_full(bezier + nsegs1*4,
-                                                     ( split_points == NULL
-                                                       ? NULL
-                                                       : split_points + nsegs1 ),
-                                                     data + splitPoint, len - splitPoint,
-                                                     recTHat1, tHat2, error, rec_max_beziers2);
-        if ( nsegs2 < 0 ) {
-#ifdef BEZIER_DEBUG
-            g_print("fit_cubic[2]: recursive call failed\n");
-#endif
-            return -1;
-        }
-
-#ifdef BEZIER_DEBUG
-        g_print("fit_cubic: success[nsegs: %d+%d=%d] on max_beziers:%u\n",
-                nsegs1, nsegs2, nsegs1 + nsegs2, max_beziers);
-#endif
-        return nsegs1 + nsegs2;
-    } else {
-        return -1;
-    }
-}
-
-
-/**
- * Fill in \a bezier[] based on the given data and tangent requirements, using
- * a least-squares fit.
- *
- * Each of tHat1 and tHat2 should be either a zero vector or a unit vector.
- * If it is zero, then bezier[1 or 2] is estimated without constraint; otherwise,
- * it bezier[1 or 2] is placed in the specified direction from bezier[0 or 3].
- *
- * \param tolerance_sq Used only for an initial guess as to tangent directions
- *   when \a tHat1 or \a tHat2 is zero.
- */
-static void
-generate_bezier(Geom::Point bezier[],
-                Geom::Point const data[], gdouble const u[], unsigned const len,
-                Geom::Point const &tHat1, Geom::Point const &tHat2,
-                double const tolerance_sq)
-{
-    bool const est1 = is_zero(tHat1);
-    bool const est2 = is_zero(tHat2);
-    Geom::Point est_tHat1( est1
-                         ? sp_darray_left_tangent(data, len, tolerance_sq)
-                         : tHat1 );
-    Geom::Point est_tHat2( est2
-                         ? sp_darray_right_tangent(data, len, tolerance_sq)
-                         : tHat2 );
-    estimate_lengths(bezier, data, u, len, est_tHat1, est_tHat2);
-    /* We find that sp_darray_right_tangent tends to produce better results
-       for our current freehand tool than full estimation. */
-    if (est1) {
-        estimate_bi(bezier, 1, data, u, len);
-        if (bezier[1] != bezier[0]) {
-            est_tHat1 = unit_vector(bezier[1] - bezier[0]);
-        }
-        estimate_lengths(bezier, data, u, len, est_tHat1, est_tHat2);
-    }
-}
-
-
-static void
-estimate_lengths(Geom::Point bezier[],
-                 Geom::Point const data[], gdouble const uPrime[], unsigned const len,
-                 Geom::Point const &tHat1, Geom::Point const &tHat2)
-{
-    double C[2][2];   /* Matrix C. */
-    double X[2];      /* Matrix X. */
-
-    /* Create the C and X matrices. */
-    C[0][0] = 0.0;
-    C[0][1] = 0.0;
-    C[1][0] = 0.0;
-    C[1][1] = 0.0;
-    X[0]    = 0.0;
-    X[1]    = 0.0;
-
-    /* First and last control points of the Bezier curve are positioned exactly at the first and
-       last data points. */
-    bezier[0] = data[0];
-    bezier[3] = data[len - 1];
-
-    for (unsigned i = 0; i < len; i++) {
-        /* Bezier control point coefficients. */
-        double const b0 = B0(uPrime[i]);
-        double const b1 = B1(uPrime[i]);
-        double const b2 = B2(uPrime[i]);
-        double const b3 = B3(uPrime[i]);
-
-        /* rhs for eqn */
-        Geom::Point const a1 = b1 * tHat1;
-        Geom::Point const a2 = b2 * tHat2;
-
-        C[0][0] += dot(a1, a1);
-        C[0][1] += dot(a1, a2);
-        C[1][0] = C[0][1];
-        C[1][1] += dot(a2, a2);
-
-        /* Additional offset to the data point from the predicted point if we were to set bezier[1]
-           to bezier[0] and bezier[2] to bezier[3]. */
-        Geom::Point const shortfall
-            = ( data[i]
-                - ( ( b0 + b1 ) * bezier[0] )
-                - ( ( b2 + b3 ) * bezier[3] ) );
-        X[0] += dot(a1, shortfall);
-        X[1] += dot(a2, shortfall);
-    }
-
-    /* We've constructed a pair of equations in the form of a matrix product C * alpha = X.
-       Now solve for alpha. */
-    double alpha_l, alpha_r;
-
-    /* Compute the determinants of C and X. */
-    double const det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
-    if ( det_C0_C1 != 0 ) {
-        /* Apparently Kramer's rule. */
-        double const det_C0_X  = C[0][0] * X[1]    - C[0][1] * X[0];
-        double const det_X_C1  = X[0]    * C[1][1] - X[1]    * C[0][1];
-        alpha_l = det_X_C1 / det_C0_C1;
-        alpha_r = det_C0_X / det_C0_C1;
-    } else {
-        /* The matrix is under-determined.  Try requiring alpha_l == alpha_r.
-         *
-         * One way of implementing the constraint alpha_l == alpha_r is to treat them as the same
-         * variable in the equations.  We can do this by adding the columns of C to form a single
-         * column, to be multiplied by alpha to give the column vector X.
-         *
-         * We try each row in turn.
-         */
-        double const c0 = C[0][0] + C[0][1];
-        if (c0 != 0) {
-            alpha_l = alpha_r = X[0] / c0;
-        } else {
-            double const c1 = C[1][0] + C[1][1];
-            if (c1 != 0) {
-                alpha_l = alpha_r = X[1] / c1;
-            } else {
-                /* Let the below code handle this. */
-                alpha_l = alpha_r = 0.;
-            }
-        }
-    }
-
-    /* If alpha negative, use the Wu/Barsky heuristic (see text).  (If alpha is 0, you get
-       coincident control points that lead to divide by zero in any subsequent
-       NewtonRaphsonRootFind() call.) */
-    /// \todo Check whether this special-casing is necessary now that 
-    /// NewtonRaphsonRootFind handles non-positive denominator.
-    if ( alpha_l < 1.0e-6 ||
-         alpha_r < 1.0e-6   )
-    {
-        alpha_l = alpha_r = ( L2( data[len - 1]
-                                  - data[0] )
-                              / 3.0 );
-    }
-
-    /* Control points 1 and 2 are positioned an alpha distance out on the tangent vectors, left and
-       right, respectively. */
-    bezier[1] = alpha_l * tHat1 + bezier[0];
-    bezier[2] = alpha_r * tHat2 + bezier[3];
-
-    return;
-}
-
-static double lensq(Geom::Point const p) {
-    return dot(p, p);
-}
-
-static void
-estimate_bi(Geom::Point bezier[4], unsigned const ei,
-            Geom::Point const data[], double const u[], unsigned const len)
-{
-    g_return_if_fail(1 <= ei && ei <= 2);
-    unsigned const oi = 3 - ei;
-    double num[2] = {0., 0.};
-    double den = 0.;
-    for (unsigned i = 0; i < len; ++i) {
-        double const ui = u[i];
-        double const b[4] = {
-            B0(ui),
-            B1(ui),
-            B2(ui),
-            B3(ui)
-        };
-
-        for (unsigned d = 0; d < 2; ++d) {
-            num[d] += b[ei] * (b[0]  * bezier[0][d] +
-                               b[oi] * bezier[oi][d] +
-                               b[3]  * bezier[3][d] +
-                               - data[i][d]);
-        }
-        den -= b[ei] * b[ei];
-    }
-
-    if (den != 0.) {
-        for (unsigned d = 0; d < 2; ++d) {
-            bezier[ei][d] = num[d] / den;
-        }
-    } else {
-        bezier[ei] = ( oi * bezier[0] + ei * bezier[3] ) / 3.;
-    }
-}
-
-/**
- * Given set of points and their parameterization, try to find a better assignment of parameter
- * values for the points.
- *
- *  \param d  Array of digitized points.
- *  \param u  Current parameter values.
- *  \param bezCurve  Current fitted curve.
- *  \param len  Number of values in both d and u arrays.
- *              Also the size of the array that is allocated for return.
- */
-static void
-reparameterize(Geom::Point const d[],
-               unsigned const len,
-               double u[],
-               BezierCurve const bezCurve)
-{
-    g_assert( 2 <= len );
-
-    unsigned const last = len - 1;
-    g_assert( bezCurve[0] == d[0] );
-    g_assert( bezCurve[3] == d[last] );
-    g_assert( u[0] == 0.0 );
-    g_assert( u[last] == 1.0 );
-    /* Otherwise, consider including 0 and last in the below loop. */
-
-    for (unsigned i = 1; i < last; i++) {
-        u[i] = NewtonRaphsonRootFind(bezCurve, d[i], u[i]);
-    }
-}
-
-/**
- *  Use Newton-Raphson iteration to find better root.
- *  
- *  \param Q  Current fitted curve
- *  \param P  Digitized point
- *  \param u  Parameter value for "P"
- *  
- *  \return Improved u
- */
-static gdouble
-NewtonRaphsonRootFind(BezierCurve const Q, Geom::Point const &P, gdouble const u)
-{
-    g_assert( 0.0 <= u );
-    g_assert( u <= 1.0 );
-
-    /* Generate control vertices for Q'. */
-    Geom::Point Q1[3];
-    for (unsigned i = 0; i < 3; i++) {
-        Q1[i] = 3.0 * ( Q[i+1] - Q[i] );
-    }
-
-    /* Generate control vertices for Q''. */
-    Geom::Point Q2[2];
-    for (unsigned i = 0; i < 2; i++) {
-        Q2[i] = 2.0 * ( Q1[i+1] - Q1[i] );
-    }
-
-    /* Compute Q(u), Q'(u) and Q''(u). */
-    Geom::Point const Q_u  = bezier_pt(3, Q, u);
-    Geom::Point const Q1_u = bezier_pt(2, Q1, u);
-    Geom::Point const Q2_u = bezier_pt(1, Q2, u);
-
-    /* Compute f(u)/f'(u), where f is the derivative wrt u of distsq(u) = 0.5 * the square of the
-       distance from P to Q(u).  Here we're using Newton-Raphson to find a stationary point in the
-       distsq(u), hopefully corresponding to a local minimum in distsq (and hence a local minimum
-       distance from P to Q(u)). */
-    Geom::Point const diff = Q_u - P;
-    double numerator = dot(diff, Q1_u);
-    double denominator = dot(Q1_u, Q1_u) + dot(diff, Q2_u);
-
-    double improved_u;
-    if ( denominator > 0. ) {
-        /* One iteration of Newton-Raphson:
-           improved_u = u - f(u)/f'(u) */
-        improved_u = u - ( numerator / denominator );
-    } else {
-        /* Using Newton-Raphson would move in the wrong direction (towards a local maximum rather
-           than local minimum), so we move an arbitrary amount in the right direction. */
-        if ( numerator > 0. ) {
-            improved_u = u * .98 - .01;
-        } else if ( numerator < 0. ) {
-            /* Deliberately asymmetrical, to reduce the chance of cycling. */
-            improved_u = .031 + u * .98;
-        } else {
-            improved_u = u;
-        }
-    }
-
-    if (!IS_FINITE(improved_u)) {
-        improved_u = u;
-    } else if ( improved_u < 0.0 ) {
-        improved_u = 0.0;
-    } else if ( improved_u > 1.0 ) {
-        improved_u = 1.0;
-    }
-
-    /* Ensure that improved_u isn't actually worse. */
-    {
-        double const diff_lensq = lensq(diff);
-        for (double proportion = .125; ; proportion += .125) {
-            if ( lensq( bezier_pt(3, Q, improved_u) - P ) > diff_lensq ) {
-                if ( proportion > 1.0 ) {
-                    //g_warning("found proportion %g", proportion);
-                    improved_u = u;
-                    break;
-                }
-                improved_u = ( ( 1 - proportion ) * improved_u  +
-                               proportion         * u            );
-            } else {
-                break;
-            }
-        }
-    }
-
-    DOUBLE_ASSERT(improved_u);
-    return improved_u;
-}
-
-/** 
- * Evaluate a Bezier curve at parameter value \a t.
- * 
- * \param degree The degree of the Bezier curve: 3 for cubic, 2 for quadratic etc.
- * \param V The control points for the Bezier curve.  Must have (\a degree+1)
- *    elements.
- * \param t The "parameter" value, specifying whereabouts along the curve to
- *    evaluate.  Typically in the range [0.0, 1.0].
- *
- * Let s = 1 - t.
- * BezierII(1, V) gives (s, t) * V, i.e. t of the way
- * from V[0] to V[1].
- * BezierII(2, V) gives (s**2, 2*s*t, t**2) * V.
- * BezierII(3, V) gives (s**3, 3 s**2 t, 3s t**2, t**3) * V.
- *
- * The derivative of BezierII(i, V) with respect to t
- * is i * BezierII(i-1, V'), where for all j, V'[j] =
- * V[j + 1] - V[j].
- */
-Geom::Point
-bezier_pt(unsigned const degree, Geom::Point const V[], gdouble const t)
-{
-    /** Pascal's triangle. */
-    static int const pascal[4][4] = {{1},
-                                     {1, 1},
-                                     {1, 2, 1},
-                                     {1, 3, 3, 1}};
-    g_assert( degree < G_N_ELEMENTS(pascal) );
-    gdouble const s = 1.0 - t;
-
-    /* Calculate powers of t and s. */
-    double spow[4];
-    double tpow[4];
-    spow[0] = 1.0; spow[1] = s;
-    tpow[0] = 1.0; tpow[1] = t;
-    for (unsigned i = 1; i < degree; ++i) {
-        spow[i + 1] = spow[i] * s;
-        tpow[i + 1] = tpow[i] * t;
-    }
-
-    Geom::Point ret = spow[degree] * V[0];
-    for (unsigned i = 1; i <= degree; ++i) {
-        ret += pascal[degree][i] * spow[degree - i] * tpow[i] * V[i];
-    }
-    return ret;
-}
-
-/*
- * ComputeLeftTangent, ComputeRightTangent, ComputeCenterTangent :
- * Approximate unit tangents at endpoints and "center" of digitized curve
- */
-
-/** 
- * Estimate the (forward) tangent at point d[first + 0.5].
- *
- * Unlike the center and right versions, this calculates the tangent in 
- * the way one might expect, i.e., wrt increasing index into d.
- * \pre (2 \<= len) and (d[0] != d[1]).
- **/
-Geom::Point
-sp_darray_left_tangent(Geom::Point const d[], unsigned const len)
-{
-    g_assert( len >= 2 );
-    g_assert( d[0] != d[1] );
-    return unit_vector( d[1] - d[0] );
-}
-
-/** 
- * Estimates the (backward) tangent at d[last - 0.5].
- *
- * \note The tangent is "backwards", i.e. it is with respect to 
- * decreasing index rather than increasing index.
- *
- * \pre 2 \<= len.
- * \pre d[len - 1] != d[len - 2].
- * \pre all[p in d] in_svg_plane(p).
- */
-static Geom::Point
-sp_darray_right_tangent(Geom::Point const d[], unsigned const len)
-{
-    g_assert( 2 <= len );
-    unsigned const last = len - 1;
-    unsigned const prev = last - 1;
-    g_assert( d[last] != d[prev] );
-    return unit_vector( d[prev] - d[last] );
-}
-
-/** 
- * Estimate the (forward) tangent at point d[0].
- *
- * Unlike the center and right versions, this calculates the tangent in 
- * the way one might expect, i.e., wrt increasing index into d.
- *
- * \pre 2 \<= len.
- * \pre d[0] != d[1].
- * \pre all[p in d] in_svg_plane(p).
- * \post is_unit_vector(ret).
- **/
-Geom::Point
-sp_darray_left_tangent(Geom::Point const d[], unsigned const len, double const tolerance_sq)
-{
-    g_assert( 2 <= len );
-    g_assert( 0 <= tolerance_sq );
-    for (unsigned i = 1;;) {
-        Geom::Point const pi(d[i]);
-        Geom::Point const t(pi - d[0]);
-        double const distsq = dot(t, t);
-        if ( tolerance_sq < distsq ) {
-            return unit_vector(t);
-        }
-        ++i;
-        if (i == len) {
-            return ( distsq == 0
-                     ? sp_darray_left_tangent(d, len)
-                     : unit_vector(t) );
-        }
-    }
-}
-
-/** 
- * Estimates the (backward) tangent at d[last].
- *
- * \note The tangent is "backwards", i.e. it is with respect to 
- * decreasing index rather than increasing index.
- *
- * \pre 2 \<= len.
- * \pre d[len - 1] != d[len - 2].
- * \pre all[p in d] in_svg_plane(p).
- */
-Geom::Point
-sp_darray_right_tangent(Geom::Point const d[], unsigned const len, double const tolerance_sq)
-{
-    g_assert( 2 <= len );
-    g_assert( 0 <= tolerance_sq );
-    unsigned const last = len - 1;
-    for (unsigned i = last - 1;; i--) {
-        Geom::Point const pi(d[i]);
-        Geom::Point const t(pi - d[last]);
-        double const distsq = dot(t, t);
-        if ( tolerance_sq < distsq ) {
-            return unit_vector(t);
-        }
-        if (i == 0) {
-            return ( distsq == 0
-                     ? sp_darray_right_tangent(d, len)
-                     : unit_vector(t) );
-        }
-    }
-}
-
-/** 
- * Estimates the (backward) tangent at d[center], by averaging the two 
- * segments connected to d[center] (and then normalizing the result).
- *
- * \note The tangent is "backwards", i.e. it is with respect to 
- * decreasing index rather than increasing index.
- *
- * \pre (0 \< center \< len - 1) and d is uniqued (at least in 
- * the immediate vicinity of \a center).
- */
-static Geom::Point
-sp_darray_center_tangent(Geom::Point const d[],
-                         unsigned const center,
-                         unsigned const len)
-{
-    g_assert( center != 0 );
-    g_assert( center < len - 1 );
-
-    Geom::Point ret;
-    if ( d[center + 1] == d[center - 1] ) {
-        /* Rotate 90 degrees in an arbitrary direction. */
-        Geom::Point const diff = d[center] - d[center - 1];
-        ret = Geom::rot90(diff);
-    } else {
-        ret = d[center - 1] - d[center + 1];
-    }
-    ret.normalize();
-    return ret;
-}
-
-
-/**
- *  Assign parameter values to digitized points using relative distances between points.
- *
- *  \pre Parameter array u must have space for \a len items.
- */
-static void
-chord_length_parameterize(Geom::Point const d[], gdouble u[], unsigned const len)
-{
-    g_return_if_fail( 2 <= len );
-
-    /* First let u[i] equal the distance travelled along the path from d[0] to d[i]. */
-    u[0] = 0.0;
-    for (unsigned i = 1; i < len; i++) {
-        double const dist = L2( d[i] - d[i-1] );
-        u[i] = u[i-1] + dist;
-    }
-
-    /* Then scale to [0.0 .. 1.0]. */
-    gdouble tot_len = u[len - 1];
-    g_return_if_fail( tot_len != 0 );
-    if (IS_FINITE(tot_len)) {
-        for (unsigned i = 1; i < len; ++i) {
-            u[i] /= tot_len;
-        }
-    } else {
-        /* We could do better, but this probably never happens anyway. */
-        for (unsigned i = 1; i < len; ++i) {
-            u[i] = i / (gdouble) ( len - 1 );
-        }
-    }
-
-    /** \todo
-     * It's been reported that u[len - 1] can differ from 1.0 on some 
-     * systems (amd64), despite it having been calculated as x / x where x 
-     * is IS_FINITE and non-zero.
-     */
-    if (u[len - 1] != 1) {
-        double const diff = u[len - 1] - 1;
-        if (fabs(diff) > 1e-13) {
-            g_warning("u[len - 1] = %19g (= 1 + %19g), expecting exactly 1",
-                      u[len - 1], diff);
-        }
-        u[len - 1] = 1;
-    }
-
-#ifdef BEZIER_DEBUG
-    g_assert( u[0] == 0.0 && u[len - 1] == 1.0 );
-    for (unsigned i = 1; i < len; i++) {
-        g_assert( u[i] >= u[i-1] );
-    }
-#endif
-}
-
-
-
-
-/**
- * Find the maximum squared distance of digitized points to fitted curve, and (if this maximum
- * error is non-zero) set \a *splitPoint to the corresponding index.
- *
- * \pre 2 \<= len.
- * \pre u[0] == 0.
- * \pre u[len - 1] == 1.0.
- * \post ((ret == 0.0)
- *        || ((*splitPoint \< len - 1)
- *            \&\& (*splitPoint != 0 || ret \< 0.0))).
- */
-static gdouble
-compute_max_error_ratio(Geom::Point const d[], double const u[], unsigned const len,
-                        BezierCurve const bezCurve, double const tolerance,
-                        unsigned *const splitPoint)
-{
-    g_assert( 2 <= len );
-    unsigned const last = len - 1;
-    g_assert( bezCurve[0] == d[0] );
-    g_assert( bezCurve[3] == d[last] );
-    g_assert( u[0] == 0.0 );
-    g_assert( u[last] == 1.0 );
-    /* I.e. assert that the error for the first & last points is zero.
-     * Otherwise we should include those points in the below loop.
-     * The assertion is also necessary to ensure 0 < splitPoint < last.
-     */
-
-    double maxDistsq = 0.0; /* Maximum error */
-    double max_hook_ratio = 0.0;
-    unsigned snap_end = 0;
-    Geom::Point prev = bezCurve[0];
-    for (unsigned i = 1; i <= last; i++) {
-        Geom::Point const curr = bezier_pt(3, bezCurve, u[i]);
-        double const distsq = lensq( curr - d[i] );
-        if ( distsq > maxDistsq ) {
-            maxDistsq = distsq;
-            *splitPoint = i;
-        }
-        double const hook_ratio = compute_hook(prev, curr, .5 * (u[i - 1] + u[i]), bezCurve, tolerance);
-        if (max_hook_ratio < hook_ratio) {
-            max_hook_ratio = hook_ratio;
-            snap_end = i;
-        }
-        prev = curr;
-    }
-
-    double const dist_ratio = sqrt(maxDistsq) / tolerance;
-    double ret;
-    if (max_hook_ratio <= dist_ratio) {
-        ret = dist_ratio;
-    } else {
-        g_assert(0 < snap_end);
-        ret = -max_hook_ratio;
-        *splitPoint = snap_end - 1;
-    }
-    g_assert( ret == 0.0
-              || ( ( *splitPoint < last )
-                   && ( *splitPoint != 0 || ret < 0. ) ) );
-    return ret;
-}
-
-/** 
- * Whereas compute_max_error_ratio() checks for itself that each data point 
- * is near some point on the curve, this function checks that each point on 
- * the curve is near some data point (or near some point on the polyline 
- * defined by the data points, or something like that: we allow for a
- * "reasonable curviness" from such a polyline).  "Reasonable curviness" 
- * means we draw a circle centred at the midpoint of a..b, of radius 
- * proportional to the length |a - b|, and require that each point on the 
- * segment of bezCurve between the parameters of a and b be within that circle.
- * If any point P on the bezCurve segment is outside of that allowable 
- * region (circle), then we return some metric that increases with the 
- * distance from P to the circle.
- *
- *  Given that this is a fairly arbitrary criterion for finding appropriate 
- *  places for sharp corners, we test only one point on bezCurve, namely 
- *  the point on bezCurve with parameter halfway between our estimated 
- *  parameters for a and b.  (Alternatives are taking the farthest of a
- *  few parameters between those of a and b, or even using a variant of 
- *  NewtonRaphsonFindRoot() for finding the maximum rather than minimum 
- *  distance.)
- */
-static double
-compute_hook(Geom::Point const &a, Geom::Point const &b, double const u, BezierCurve const bezCurve,
-             double const tolerance)
-{
-    Geom::Point const P = bezier_pt(3, bezCurve, u);
-    Geom::Point const diff = .5 * (a + b) - P;
-    double const dist = Geom::L2(diff);
-    if (dist < tolerance) {
-        return 0;
-    }
-    double const allowed = Geom::L2(b - a) + tolerance;
-    return dist / allowed;
-    /** \todo 
-     * effic: Hooks are very rare.  We could start by comparing 
-     * distsq, only resorting to the more expensive L2 in cases of 
-     * uncertainty.
-     */
-}
-
-/*
-  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 :

diff --git a/src/display/bezier-utils.h b/src/display/bezier-utils.h
deleted file mode 100644 (file)
index fc5fa61..0000000
--- a/src/display/bezier-utils.h
+++ /dev/null
@@ -1,50 +0,0 @@
-#ifndef __SP_BEZIER_UTILS_H__
-#define __SP_BEZIER_UTILS_H__
-
-/*
- * An Algorithm for Automatically Fitting Digitized Curves
- * by Philip J. Schneider
- * from "Graphics Gems", Academic Press, 1990
- *
- * Authors:
- *   Philip J. Schneider
- *   Lauris Kaplinski <lauris@ximian.com>
- *
- * Copyright (C) 1990 Philip J. Schneider
- * Copyright (C) 2001 Lauris Kaplinski and Ximian, Inc.
- *
- * Released under GNU GPL
- */
-
-#include <2geom/forward.h>
-#include <glib/gtypes.h>
-
-/* Bezier approximation utils */
-Geom::Point bezier_pt(unsigned degree, Geom::Point const V[], gdouble t);
-
-gint sp_bezier_fit_cubic(Geom::Point bezier[], Geom::Point const data[], gint len, gdouble error);
-
-gint sp_bezier_fit_cubic_r(Geom::Point bezier[], Geom::Point const data[], gint len, gdouble error,
-                           unsigned max_beziers);
-
-gint sp_bezier_fit_cubic_full(Geom::Point bezier[], int split_points[], Geom::Point const data[], gint len,
-                              Geom::Point const &tHat1, Geom::Point const &tHat2,
-                              gdouble error, unsigned max_beziers);
-
-Geom::Point sp_darray_left_tangent(Geom::Point const d[], unsigned const len);
-Geom::Point sp_darray_left_tangent(Geom::Point const d[], unsigned const len, double const tolerance_sq);
-Geom::Point sp_darray_right_tangent(Geom::Point const d[], unsigned const length, double const tolerance_sq);
-
-
-#endif /* __SP_BEZIER_UTILS_H__ */
-
-/*
-  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 :