From 38b5dc2ffd01219d3c9241b404ac1ae4f9ec2f33 Mon Sep 17 00:00:00 2001 From: johanengelen Date: Sat, 13 Dec 2008 21:21:10 +0000 Subject: [PATCH] remove bezier-utils. now we use 2geom's one. remove bezier-utils-test as well. --- src/display/CMakeLists.txt | 2 - src/display/Makefile_insert | 3 - src/display/bezier-utils-test.cpp | 334 ---------- src/display/bezier-utils-test.h | 349 ----------- src/display/bezier-utils-work.txt | 34 -- src/display/bezier-utils.cpp | 984 ------------------------------ src/display/bezier-utils.h | 50 -- 7 files changed, 1756 deletions(-) delete mode 100644 src/display/bezier-utils-test.cpp delete mode 100644 src/display/bezier-utils-test.h delete mode 100644 src/display/bezier-utils-work.txt delete mode 100644 src/display/bezier-utils.cpp delete mode 100644 src/display/bezier-utils.h diff --git a/src/display/CMakeLists.txt b/src/display/CMakeLists.txt index f8dc72afe..a1d61de94 100644 --- 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 e6877ccfb..ea1994040 100644 --- 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 index 4d6313886..000000000 --- a/src/display/bezier-utils-test.cpp +++ /dev/null @@ -1,334 +0,0 @@ -#include "../utest/utest.h" -#include -#include /* 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 -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 index 3be643568..000000000 --- a/src/display/bezier-utils-test.h +++ /dev/null @@ -1,349 +0,0 @@ -#include - -#include -#include /* NR_DF_TEST_CLOSE */ -#include - -/* 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 index 75ec5be86..000000000 --- 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 index b538f7fba..000000000 --- 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 - * Peter Moulder - * - * 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 -#endif - -#ifdef HAVE_IEEEFP_H -# include -#endif - -#include // g_assert() -#include -#include -#include "bezier-utils.h" -#include - -#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 index fc5fa61a4..000000000 --- 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 - * - * Copyright (C) 1990 Philip J. Schneider - * Copyright (C) 2001 Lauris Kaplinski and Ximian, Inc. - * - * Released under GNU GPL - */ - -#include <2geom/forward.h> -#include - -/* 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 : -- 2.30.2