src/display/bezier-utils-test.h

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