src/display/bezier-utils-test.cpp

   1 #include "../utest/utest.h"
   2 #include <glib.h>
   3 #include <libnr/nr-macros.h> /* NR_DF_TEST_CLOSE */
   4
   5 /* mental disclaims all responsibility for this evil idea for testing
   6    static functions.  The main disadvantages are that we retain the
   7    #define's and `using' directives of the included file. */
   8 #include "bezier-utils.cpp"
   9
  10 using NR::Point;
  11
  12 static bool range_approx_equal(double const a[], double const b[], unsigned len);
  13
  14 /* (Returns false if NaN encountered.) */
  15 template<class T>
  16 static bool range_equal(T const a[], T const b[], unsigned len) {
  17     for (unsigned i = 0; i < len; ++i) {
  18         if ( a[i] != b[i] ) {
  19             return false;
  20         }
  21     }
  22     return true;
  23 }
  24
  25 inline bool point_approx_equal(NR::Point const &a, NR::Point const &b, double const eps)
  26 {
  27     using NR::X; using NR::Y;
  28     return ( NR_DF_TEST_CLOSE(a[X], b[X], eps) &&
  29              NR_DF_TEST_CLOSE(a[Y], b[Y], eps) );
  30 }
  31
  32 static inline double square(double const x) {
  33     return x * x;
  34 }
  35
  36 /** Determine whether the found control points are the same as previously found on some developer's
  37     machine.  Doesn't call utest__fail, just writes a message to stdout for diagnostic purposes:
  38     the most important test is that the root-mean-square of errors in the estimation are low rather
  39     than that the control points found are the same.
  40 **/
  41 static void compare_ctlpts(Point const est_b[], Point const exp_est_b[])
  42 {
  43     unsigned diff_mask = 0;
  44     for (unsigned i = 0; i < 4; ++i) {
  45         for (unsigned d = 0; d < 2; ++d) {
  46             if ( fabs( est_b[i][d] - exp_est_b[i][d] ) > 1.1e-5 ) {
  47                 diff_mask |= 1 << ( i * 2 + d );
  48             }
  49         }
  50     }
  51     if ( diff_mask != 0 ) {
  52         printf("Warning: got different control points from previously-coded (diffs=0x%x).\n",
  53                diff_mask);
  54         printf(" Previous:");
  55         for (unsigned i = 0; i < 4; ++i) {
  56             printf(" (%g, %g)", exp_est_b[i][0], exp_est_b[i][1]); // localizing ok
  57         }
  58         putchar('\n');
  59         printf(" Found:   ");
  60         for (unsigned i = 0; i < 4; ++i) {
  61             printf(" (%g, %g)", est_b[i][0], est_b[i][1]); // localizing ok
  62         }
  63         putchar('\n');
  64     }
  65 }
  66
  67 static void compare_rms(Point const est_b[], double const t[], Point const d[], unsigned const n,
  68                         double const exp_rms_error)
  69 {
  70     double sum_errsq = 0.0;
  71     for (unsigned i = 0; i < n; ++i) {
  72         Point const fit_pt = bezier_pt(3, est_b, t[i]);
  73         Point const diff = fit_pt - d[i];
  74         sum_errsq += dot(diff, diff);
  75     }
  76     double const rms_error = sqrt( sum_errsq / n );
  77     UTEST_ASSERT( rms_error <= exp_rms_error + 1.1e-6 );
  78     if ( rms_error < exp_rms_error - 1.1e-6 ) {
  79         /* The fitter code appears to have improved [or the floating point calculations differ
  80            on this machine from the machine where exp_rms_error was calculated]. */
  81         printf("N.B. rms_error regression requirement can be decreased: have rms_error=%g.\n", rms_error); // localizing ok
  82     }
  83 }
  84
  85 int main(int argc, char *argv[]) {
  86     utest_start("bezier-utils.cpp");
  87
  88     UTEST_TEST("copy_without_nans_or_adjacent_duplicates") {
  89         NR::Point const src[] = {
  90             Point(2., 3.),
  91             Point(2., 3.),
  92             Point(0., 0.),
  93             Point(2., 3.),
  94             Point(2., 3.),
  95             Point(1., 9.),
  96             Point(1., 9.)
  97         };
  98         Point const exp_dest[] = {
  99             Point(2., 3.),
 100             Point(0., 0.),
 101             Point(2., 3.),
 102             Point(1., 9.)
 103         };
 104         g_assert( G_N_ELEMENTS(src) == 7 );
 105         Point dest[7];
 106         struct tst {
 107             unsigned src_ix0;
 108             unsigned src_len;
 109             unsigned exp_dest_ix0;
 110             unsigned exp_dest_len;
 111         } const test_data[] = {
 112             /* src start ix, src len, exp_dest start ix, exp dest len */
 113             {0, 0, 0, 0},
 114             {2, 1, 1, 1},
 115             {0, 1, 0, 1},
 116             {0, 2, 0, 1},
 117             {0, 3, 0, 2},
 118             {1, 3, 0, 3},
 119             {0, 5, 0, 3},
 120             {0, 6, 0, 4},
 121             {0, 7, 0, 4}
 122         };
 123         for (unsigned i = 0 ; i < G_N_ELEMENTS(test_data) ; ++i) {
 124             tst const &t = test_data[i];
 125             UTEST_ASSERT( t.exp_dest_len
 126                           == copy_without_nans_or_adjacent_duplicates(src + t.src_ix0,
 127                                                                       t.src_len,
 128                                                                       dest) );
 129             UTEST_ASSERT(range_equal(dest,
 130                                      exp_dest + t.exp_dest_ix0,
 131                                      t.exp_dest_len));
 132         }
 133     }
 134
 135     UTEST_TEST("bezier_pt(1)") {
 136         Point const a[] = {Point(2.0, 4.0),
 137                            Point(1.0, 8.0)};
 138         UTEST_ASSERT( bezier_pt(1, a, 0.0) == a[0] );
 139         UTEST_ASSERT( bezier_pt(1, a, 1.0) == a[1] );
 140         UTEST_ASSERT( bezier_pt(1, a, 0.5) == Point(1.5, 6.0) );
 141         double const t[] = {0.5, 0.25, 0.3, 0.6};
 142         for (unsigned i = 0; i < G_N_ELEMENTS(t); ++i) {
 143             double const ti = t[i], si = 1.0 - ti;
 144             UTEST_ASSERT( bezier_pt(1, a, ti) == si * a[0] + ti * a[1] );
 145         }
 146     }
 147
 148     UTEST_TEST("bezier_pt(2)") {
 149         Point const b[] = {Point(1.0, 2.0),
 150                            Point(8.0, 4.0),
 151                            Point(3.0, 1.0)};
 152         UTEST_ASSERT( bezier_pt(2, b, 0.0) == b[0] );
 153         UTEST_ASSERT( bezier_pt(2, b, 1.0) == b[2] );
 154         UTEST_ASSERT( bezier_pt(2, b, 0.5) == Point(5.0, 2.75) );
 155         double const t[] = {0.5, 0.25, 0.3, 0.6};
 156         for (unsigned i = 0; i < G_N_ELEMENTS(t); ++i) {
 157             double const ti = t[i], si = 1.0 - ti;
 158             Point const exp_pt( si*si * b[0] + 2*si*ti * b[1] + ti*ti * b[2] );
 159             Point const pt(bezier_pt(2, b, ti));
 160             UTEST_ASSERT(point_approx_equal(pt, exp_pt, 1e-11));
 161         }
 162     }
 163
 164     Point const c[] = {Point(1.0, 2.0),
 165                        Point(8.0, 4.0),
 166                        Point(3.0, 1.0),
 167                        Point(-2.0, -4.0)};
 168     UTEST_TEST("bezier_pt(3)") {
 169         UTEST_ASSERT( bezier_pt(3, c, 0.0) == c[0] );
 170         UTEST_ASSERT( bezier_pt(3, c, 1.0) == c[3] );
 171         UTEST_ASSERT( bezier_pt(3, c, 0.5) == Point(4.0, 13.0/8.0) );
 172         double const t[] = {0.5, 0.25, 0.3, 0.6};
 173         for (unsigned i = 0; i < G_N_ELEMENTS(t); ++i) {
 174             double const ti = t[i], si = 1.0 - ti;
 175             UTEST_ASSERT( LInfty( bezier_pt(3, c, ti)
 176                                   - ( si*si*si * c[0] +
 177                                       3*si*si*ti * c[1] +
 178                                       3*si*ti*ti * c[2] +
 179                                       ti*ti*ti * c[3] ) )
 180                           < 1e-4 );
 181         }
 182     }
 183
 184     struct Err_tst {
 185         Point pt;
 186         double u;
 187         double err;
 188     } const err_tst[] = {
 189         {c[0], 0.0, 0.0},
 190         {Point(4.0, 13.0/8.0), 0.5, 0.0},
 191         {Point(4.0, 2.0), 0.5, 9.0/64.0},
 192         {Point(3.0, 2.0), 0.5, 1.0 + 9.0/64.0},
 193         {Point(6.0, 2.0), 0.5, 4.0 + 9.0/64.0},
 194         {c[3], 1.0, 0.0},
 195     };
 196
 197     UTEST_TEST("compute_max_error_ratio") {
 198         Point d[G_N_ELEMENTS(err_tst)];
 199         double u[G_N_ELEMENTS(err_tst)];
 200         for (unsigned i = 0; i < G_N_ELEMENTS(err_tst); ++i) {
 201             Err_tst const &t = err_tst[i];
 202             d[i] = t.pt;
 203             u[i] = t.u;
 204         }
 205         g_assert( G_N_ELEMENTS(u) == G_N_ELEMENTS(d) );
 206         unsigned max_ix = ~0u;
 207         double const err_ratio = compute_max_error_ratio(d, u, G_N_ELEMENTS(d), c, 1.0, &max_ix);
 208         UTEST_ASSERT( fabs( sqrt(err_tst[4].err) - err_ratio ) < 1e-12 );
 209         UTEST_ASSERT( max_ix == 4 );
 210     }
 211
 212     UTEST_TEST("chord_length_parameterize") {
 213         /* n == 2 */
 214         {
 215             Point const d[] = {Point(2.9415, -5.8149),
 216                                Point(23.021, 4.9814)};
 217             double u[G_N_ELEMENTS(d)];
 218             double const exp_u[] = {0.0, 1.0};
 219             g_assert( G_N_ELEMENTS(u) == G_N_ELEMENTS(exp_u) );
 220             chord_length_parameterize(d, u, G_N_ELEMENTS(d));
 221             UTEST_ASSERT(range_equal(u, exp_u, G_N_ELEMENTS(exp_u)));
 222         }
 223
 224         /* Straight line. */
 225         {
 226             double const exp_u[] = {0.0, 0.1829, 0.2105, 0.2105, 0.619, 0.815, 0.999, 1.0};
 227             unsigned const n = G_N_ELEMENTS(exp_u);
 228             Point d[n];
 229             double u[n];
 230             Point const a(-23.985, 4.915), b(4.9127, 5.203);
 231             for (unsigned i = 0; i < n; ++i) {
 232                 double bi = exp_u[i], ai = 1.0 - bi;
 233                 d[i] = ai * a  +  bi * b;
 234             }
 235             chord_length_parameterize(d, u, n);
 236             UTEST_ASSERT(range_approx_equal(u, exp_u, n));
 237         }
 238     }
 239
 240     /* Feed it some points that can be fit exactly with a single bezier segment, and see how
 241        well it manages. */
 242     Point const src_b[4] = {Point(5., -3.),
 243                             Point(8., 0.),
 244                             Point(4., 2.),
 245                             Point(3., 3.)};
 246     double const t[] = {0.0, .001, .03, .05, .09, .13, .18, .25, .29, .33, .39, .44,
 247                         .51, .57, .62, .69, .75, .81, .91, .93, .97, .98, .999, 1.0};
 248     unsigned const n = G_N_ELEMENTS(t);
 249     Point d[n];
 250     for (unsigned i = 0; i < n; ++i) {
 251         d[i] = bezier_pt(3, src_b, t[i]);
 252     }
 253     Point const tHat1(unit_vector( src_b[1] - src_b[0] ));
 254     Point const tHat2(unit_vector( src_b[2] - src_b[3] ));
 255
 256     UTEST_TEST("generate_bezier") {
 257         Point est_b[4];
 258         generate_bezier(est_b, d, t, n, tHat1, tHat2, 1.0);
 259
 260         compare_ctlpts(est_b, src_b);
 261
 262         /* We're being unfair here in using our t[] rather than best t[] for est_b: we
 263            may over-estimate RMS of errors. */
 264         compare_rms(est_b, t, d, n, 1e-8);
 265     }
 266
 267     UTEST_TEST("sp_bezier_fit_cubic_full") {
 268         Point est_b[4];
 269         int splitpoints[2];
 270         gint const succ = sp_bezier_fit_cubic_full(est_b, splitpoints, d, n, tHat1, tHat2, square(1.2), 1);
 271         UTEST_ASSERT( succ == 1 );
 272
 273         Point const exp_est_b[4] = {
 274             Point(5.000000, -3.000000),
 275             Point(7.5753, -0.4247),
 276             Point(4.77533, 1.22467),
 277             Point(3, 3)
 278         };
 279         compare_ctlpts(est_b, exp_est_b);
 280
 281         /* We're being unfair here in using our t[] rather than best t[] for est_b: we
 282            may over-estimate RMS of errors. */
 283         compare_rms(est_b, t, d, n, .307911);
 284     }
 285
 286     UTEST_TEST("sp_bezier_fit_cubic") {
 287         Point est_b[4];
 288         gint const succ = sp_bezier_fit_cubic(est_b, d, n, square(1.2));
 289         UTEST_ASSERT( succ == 1 );
 290
 291         Point const exp_est_b[4] = {
 292             Point(5.000000, -3.000000),
 293             Point(7.57134, -0.423509),
 294             Point(4.77929, 1.22426),
 295             Point(3, 3)
 296         };
 297         compare_ctlpts(est_b, exp_est_b);
 298
 299 #if 1 /* A change has been made to right_tangent.  I believe that usually this change
 300          will result in better fitting, but it won't do as well for this example where
 301          we happen to be feeding a t=0.999 point to the fitter. */
 302         printf("TODO: Update this test case for revised right_tangent implementation.\n");
 303         /* In particular, have a test case to show whether the new implementation
 304            really is likely to be better on average. */
 305 #else
 306         /* We're being unfair here in using our t[] rather than best t[] for est_b: we
 307            may over-estimate RMS of errors. */
 308         compare_rms(est_b, t, d, n, .307983);
 309 #endif
 310     }
 311
 312     return !utest_end();
 313 }
 314
 315 /* (Returns false if NaN encountered.) */
 316 static bool range_approx_equal(double const a[], double const b[], unsigned const len) {
 317     for (unsigned i = 0; i < len; ++i) {
 318         if (!( fabs( a[i] - b[i] ) < 1e-4 )) {
 319             return false;
 320         }
 321     }
 322     return true;
 323 }
 324
 325 /*
 326   Local Variables:
 327   mode:c++
 328   c-file-style:"stroustrup"
 329   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 330   indent-tabs-mode:nil
 331   fill-column:99
 332   End:
 333 */
 334 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :