Code-Eli  0.3.6
bezier_curve_test_suite.hpp
Go to the documentation of this file.
1 /*********************************************************************************
2 * Copyright (c) 2013 David D. Marshall <ddmarsha@calpoly.edu>
3 *
4 * All rights reserved. This program and the accompanying materials
5 * are made available under the terms of the Eclipse Public License v1.0
6 * which accompanies this distribution, and is available at
7 * http://www.eclipse.org/legal/epl-v10.html
8 *
9 * Contributors:
10 * David D. Marshall - initial code and implementation
11 ********************************************************************************/
12 
13 #ifndef bezier_curve_test_suite_hpp
14 #define bezier_curve_test_suite_hpp
15 
16 #include <cmath> // std::pow, std::exp
17 
18 #include <typeinfo> // typeid
19 #include <string> // std::string
20 #include <sstream> // std::stringstream
21 #include <iomanip> // std::setw
22 #include <limits> // std::numeric_limits
23 
24 #include "eli/constants/math.hpp"
29 
30 template<typename data__>
31 class bezier_curve_test_suite : public Test::Suite
32 {
33  private:
40 
41  protected:
42  void AddTests(const float &)
43  {
44  // add the tests
62  }
63  void AddTests(const double &)
64  {
65  // add the tests
83  }
84  void AddTests(const long double &)
85  {
86  // add the tests
104  }
105 
106  public:
108  {
109  AddTests(data__());
110  }
112  {
113  }
114 
115  private:
116  void octave_print(int figno, const std::vector<point_type> &pts, const bezier_type &bez) const
117  {
118  size_t i;
119 
120  std::cout << "figure(" << figno << ");" << std::endl;
121  std::cout << "xpts=[" << pts[0].x();
122  for (i=1; i<pts.size(); ++i)
123  std::cout << ", " << pts[i].x();
124  std::cout << "];" << std::endl;
125  std::cout << "ypts=[" << pts[0].y();
126  for (i=1; i<pts.size(); ++i)
127  std::cout << ", " << pts[i].y();
128  std::cout << "];" << std::endl;
129 
130  std::vector<data__> t(101);
131  for (i=0; i<t.size(); ++i)
132  t[i]=static_cast<data__>(i)/(t.size()-1);
133 
134  std::cout << "xint=[" << bez.f(t[0])(0);
135  for (i=1; i<t.size(); ++i)
136  std::cout << ", " << bez.f(t[i])(0);
137  std::cout << "];" << std::endl;
138  std::cout << "yint=[" << bez.f(t[0])(1);
139  for (i=1; i<t.size(); ++i)
140  std::cout << ", " << bez.f(t[i])(1);
141  std::cout << "];" << std::endl;
142 
143  std::cout << "plot(xpts, ypts, 'bo', xint, yint, 'k-');" << std::endl;
144  }
145 
146  void create_circle(std::vector<point_type> &pts)
147  {
148  // NOTE: This will not create a closed circle, the last point will be
149  // the point just prior to the 2*pi point
150  size_t n=pts.size();
151  for (size_t i=0; i<n; ++i)
152  {
153  data__ theta(eli::constants::math<data__>::two_pi()*static_cast<data__>(i)/n);
154  pts[i](0)=std::cos(theta);
155  pts[i](1)=std::sin(theta);
156  pts[i](2)=0;
157  }
158  }
159 
161  {
162  bezier_type bc1, bc2;
163  point_type cntrl_in[4];
164  index_type i;
165 
166  // test default constructor then set control points
167  cntrl_in[0] << 2, 2, 0;
168  cntrl_in[1] << 1, 1, 0;
169  cntrl_in[2] << 3, 0, 0;
170  cntrl_in[3] << 4, 1, 0;
171  bc1.resize(3);
172  for (i=0; i<4; ++i)
173  {
174  bc1.set_control_point(cntrl_in[i], i);
175  }
176  for (i=0; i<4; ++i)
177  {
178  TEST_ASSERT(cntrl_in[i]==bc1.get_control_point(i));
179  }
180 
181  // test copy ctr
182  bezier_type bcc1(bc1);
183  for (i=0; i<4; ++i)
184  {
185  TEST_ASSERT(bc1.get_control_point(i)==bcc1.get_control_point(i));
186  }
187 
188  // test assignment operator
189  bc2=bc1;
190  for (i=0; i<4; ++i)
191  {
192  TEST_ASSERT(bc1.get_control_point(i)==bc2.get_control_point(i));
193  }
194 
195  // test order
196  TEST_ASSERT(bc2.degree()==4-1);
197  }
198 
200  {
201  bezier_type bc;
202  point_type cntrl_in[4];
203  index_type i;
204 
205  // test default constructor then set control points
206  cntrl_in[0] << 2, 2, 0;
207  cntrl_in[1] << 1, 1, 0;
208  cntrl_in[2] << 3, 0, 0;
209  cntrl_in[3] << 4, 1, 0;
210  bc.resize(3);
211  for (i=0; i<4; ++i)
212  {
213  bc.set_control_point(cntrl_in[i], i);
214  }
215 
216  // test bounding box
217  typename bezier_type::bounding_box_type bb;
218  point_type pmin_ref, pmax_ref;
219 
220  bc.get_bounding_box(bb);
221  pmin_ref << 1, 0, 0;
222  pmax_ref << 4, 2, 0;
223  TEST_ASSERT(bb.get_min()==pmin_ref);
224  TEST_ASSERT(bb.get_max()==pmax_ref);
225  }
226 
228  {
229  point_type cntrl_in[4];
230  data_type eps(std::numeric_limits<data__>::epsilon());
231 
232  // set control points
233  cntrl_in[0] << 2.0, 2.0, 0.0;
234  cntrl_in[1] << 1.0, 1.5, 0.0;
235  cntrl_in[2] << 3.5, 0.0, 0.0;
236  cntrl_in[3] << 4.0, 1.0, 0.0;
237 
238  bezier_type bc1(3), bc2;
239  point_type eval_out, eval_ref;
240  data_type t;
241 
242  // set control points
243  for (index_type i=0; i<4; ++i)
244  {
245  bc1.set_control_point(cntrl_in[i], i);
246  }
247 
248  // test translation
249  {
250  point_type trans;
251 
252  // set up translation vector and apply
253  bc2=bc1;
254  trans << 2, 1, 3;
255  bc2.translate(trans);
256 
257  // test evaluations
258  t=0;
259  eval_out=bc2.f(t);
260  eval_ref=trans+cntrl_in[0];
261  TEST_ASSERT(eval_out==eval_ref);
262  t=1;
263  eval_out=bc2.f(t);
264  eval_ref=trans+cntrl_in[3];
265  TEST_ASSERT(eval_out==eval_ref);
266 
267  // test evaluation at interior point
268  t=static_cast<data__>(0.45);
269  eval_out=bc2.f(t);
270  eval_ref << static_cast<data__>(2.2750625), static_cast<data__>(1.0364375), static_cast<data__>(0);
271  eval_ref+=trans;
272  TEST_ASSERT((eval_out-eval_ref).norm()<5e3*eps);
273  }
274 
275  // test rotation about origin
276  {
277  typename bezier_type::rotation_matrix_type rmat;
278 
279  // set up rotation and apply
280  data_type one(1);
281  bc2=bc1;
282  rmat << std::cos(one), 0, -std::sin(one),
283  0, one, 0,
284  std::sin(one), 0, std::cos(one);
285  bc2.rotate(rmat);
286 
287  // test evaluations
288  t=0;
289  eval_out=bc2.f(t);
290  eval_ref=(cntrl_in[0]*rmat.transpose());
291  TEST_ASSERT(eval_out==eval_ref);
292  t=1;
293  eval_out=bc2.f(t);
294  eval_ref=(cntrl_in[3]*rmat.transpose());
295  TEST_ASSERT(eval_out==eval_ref);
296 
297  // test evaluation at interior point
298  t=static_cast<data__>(0.45);
299  eval_out=bc2.f(t);
300  eval_ref << static_cast<data__>(2.2750625), static_cast<data__>(1.0364375), static_cast<data__>(0);
301  eval_ref*=rmat.transpose();
302  TEST_ASSERT((eval_out-eval_ref).norm()<5e3*eps);
303  }
304 
305  // test rotation about point
306  {
307  point_type rorig;
308  typename bezier_type::rotation_matrix_type rmat;
309 
310  // set up rotation and apply
311  data_type one(1);
312  bc2=bc1;
313  rorig << 2, 1, 3;
314  rmat << std::cos(one), 0, -std::sin(one),
315  0, one, 0,
316  std::sin(one), 0, std::cos(one);
317  bc2.rotate(rmat, rorig);
318 
319  // test evaluations
320  t=0;
321  eval_out=bc2.f(t);
322  eval_ref=rorig+(cntrl_in[0]-rorig)*rmat.transpose();
323  TEST_ASSERT(eval_out==eval_ref);
324  t=1;
325  eval_out=bc2.f(t);
326  eval_ref=rorig+(cntrl_in[3]-rorig)*rmat.transpose();
327  TEST_ASSERT(eval_out==eval_ref);
328 
329  // test evaluation at interior point
330  t=static_cast<data__>(0.45);
331  eval_out=bc2.f(t);
332  eval_ref << static_cast<data__>(2.2750625), static_cast<data__>(1.0364375), static_cast<data__>(0);
333  eval_ref=rorig+(eval_ref-rorig)*rmat.transpose();
334  TEST_ASSERT((eval_out-eval_ref).norm()<5e3*eps);
335  }
336  }
337 
339  {
340  point_type cntrl_in[4];
341  data_type eps(std::numeric_limits<data__>::epsilon());
342 
343  // set control points
344  cntrl_in[0] << 2.0, 2.0, 0.0;
345  cntrl_in[1] << 1.0, 1.5, 0.0;
346  cntrl_in[2] << 3.5, 0.0, 0.0;
347  cntrl_in[3] << 4.0, 1.0, 0.0;
348 
349  bezier_type bc1(3);
350  point_type eval_out, eval_ref;
351  data_type t;
352 
353  // set control points
354  for (index_type i=0; i<4; ++i)
355  {
356  bc1.set_control_point(cntrl_in[i], i);
357  }
358 
359  // test evaluation at end points
360  t=0;
361  eval_out=bc1.f(t);
362  TEST_ASSERT(eval_out==cntrl_in[0]);
363  t=1;
364  eval_out=bc1.f(t);
365  TEST_ASSERT(eval_out==cntrl_in[3]);
366 
367  // test evaluation at interior point
368  t=static_cast<data__>(0.45);
369  eval_out=bc1.f(t);
370  eval_ref << static_cast<data__>(2.2750625), static_cast<data__>(1.0364375), static_cast<data__>(0);
371  TEST_ASSERT((eval_out-eval_ref).norm()<5e3*eps);
372  }
373 
375  {
376  point_type cntrl_in[4];
377  data_type eps(std::numeric_limits<data__>::epsilon());
378 
379  // set control points
380  cntrl_in[0] << 2.0, 2.0, 0.0;
381  cntrl_in[1] << 1.0, 1.5, 0.0;
382  cntrl_in[2] << 3.5, 0.0, 0.0;
383  cntrl_in[3] << 4.0, 1.0, 0.0;
384 
385  bezier_type bc1(3);
386  point_type eval_out, eval_ref;
387  data_type t;
388 
389  // set control points
390  for (index_type i=0; i<4; ++i)
391  {
392  bc1.set_control_point(cntrl_in[i], i);
393  }
394 
395  // test 1st derivative at end points
396  t=0;
397  eval_out=bc1.fp(t);
398  eval_ref << -3, -1.5, 0;
399  TEST_ASSERT(eval_out==eval_ref);
400  t=1;
401  eval_out=bc1.fp(t);
402  eval_ref << 1.5, 3, 0;
403  TEST_ASSERT(eval_out==eval_ref);
404 
405  // test 1st derivative at interior point
406  t=static_cast<data__>(0.45);
407  eval_out=bc1.fp(t);
408  eval_ref << static_cast<data__>(3.10875), static_cast<data__>(-2.07375), static_cast<data__>(0);
409  TEST_ASSERT((eval_out-eval_ref).norm() < 5e3*eps);
410  }
411 
413  {
414  point_type cntrl_in[4];
415  data_type eps(std::numeric_limits<data__>::epsilon());
416 
417  // set control points
418  cntrl_in[0] << 2.0, 2.0, 0.0;
419  cntrl_in[1] << 1.0, 1.5, 0.0;
420  cntrl_in[2] << 3.5, 0.0, 0.0;
421  cntrl_in[3] << 4.0, 1.0, 0.0;
422 
423  bezier_type bc1(3);
424  point_type eval_out, eval_ref;
425  data_type t;
426 
427  // set control points
428  for (index_type i=0; i<4; ++i)
429  {
430  bc1.set_control_point(cntrl_in[i], i);
431  }
432 
433  // test 2nd derivative at end points
434  t=0;
435  eval_out=bc1.fpp(t);
436  eval_ref << 21, -6, 0;
437  TEST_ASSERT(eval_out==eval_ref);
438  t=1;
439  eval_out=bc1.fpp(t);
440  eval_ref << -12, 15, 0;
441  TEST_ASSERT(eval_out==eval_ref);
442 
443  // test 2nd derivative at interior point
444  t=static_cast<data__>(0.45);
445  eval_out=bc1.fpp(t);
446  eval_ref << static_cast<data__>(6.15), static_cast<data__>(3.45), static_cast<data__>(0);
447  TEST_ASSERT_DELTA((eval_out-eval_ref).norm(), 0, 1e4*eps);
448  }
449 
451  {
452  point_type cntrl_in[4];
453  data_type eps(std::numeric_limits<data__>::epsilon());
454 
455  // set control points
456  cntrl_in[0] << 2.0, 2.0, 0.0;
457  cntrl_in[1] << 1.0, 1.5, 0.0;
458  cntrl_in[2] << 3.5, 0.0, 0.0;
459  cntrl_in[3] << 4.0, 1.0, 0.0;
460 
461  bezier_type bc1(3);
462  point_type eval_out, eval_ref;
463  data_type t;
464 
465  // set control points
466  for (index_type i=0; i<4; ++i)
467  {
468  bc1.set_control_point(cntrl_in[i], i);
469  }
470 
471  // test 3rd derivative at end points
472  t=0;
473  eval_out=bc1.fppp(t);
474  eval_ref << -33, 21, 0;
475  TEST_ASSERT(eval_out==eval_ref);
476  t=1;
477  eval_out=bc1.fppp(t);
478  eval_ref << -33, 21, 0;
479  TEST_ASSERT(eval_out==eval_ref);
480 
481  // test 3rd derivative at interior point
482  t=static_cast<data__>(0.45);
483  eval_out=bc1.fppp(t);
484  eval_ref << -33, 21, 0;
485  TEST_ASSERT(eval_out==eval_ref);
486 
487  // test curvature at end points
488  data_type curv_out, curv_ref;
489  t=0;
490  eli::geom::curve::curvature(curv_out, bc1, t);
491  curv_ref=static_cast<data__>(1.31182654679988);
492  TEST_ASSERT_DELTA(curv_out, curv_ref, std::max(static_cast<data__>(1e-13), static_cast<data__>(1e2)*eps));
493  t=1;
494  eli::geom::curve::curvature(curv_out, bc1, t);
495  curv_ref=static_cast<data__>(1.55034046439985);
496  TEST_ASSERT_DELTA(curv_out, curv_ref, std::max(static_cast<data__>(1e-13), static_cast<data__>(1e2)*eps));
497 
498  // test curvature at interior point
499  t=static_cast<data__>(0.45);
500  eli::geom::curve::curvature(curv_out, bc1, t);
501  curv_ref=static_cast<data__>(0.449908807121445);
502  TEST_ASSERT_DELTA(curv_out, curv_ref, std::max(static_cast<data__>(1e-13), static_cast<data__>(1e2)*eps));
503  }
504 
506  {
507  point_type cntrl_in[4];
508  data_type eps(std::numeric_limits<data__>::epsilon());
509 
510  // set control points
511  cntrl_in[0] << 2.0, 2.0, 0.0;
512  cntrl_in[1] << 1.0, 1.5, 0.0;
513  cntrl_in[2] << 3.5, 0.0, 0.0;
514  cntrl_in[3] << 4.0, 1.0, 0.0;
515 
516  bezier_type bc(3);
517  point_type t_dir, n_dir, b_dir, t_dir_ref;
518  data_type t;
519 
520  // set control points
521  for (index_type i=0; i<4; ++i)
522  {
523  bc.set_control_point(cntrl_in[i], i);
524  }
525 
526  t=0.5;
527  t_dir=bc.tangent(t);
528  t_dir_ref << static_cast<data_type>(0.874157), static_cast<data_type>(-0.485643), 0;
529  TEST_ASSERT((t_dir-t_dir_ref).norm()<1e-5);
530 
531  bc.frenet_serret_frame(t_dir, n_dir, b_dir, t);
532  TEST_ASSERT((t_dir.cross(n_dir)-b_dir).norm()<2*eps);
533  }
534 
536  {
537  point_type cntrl_in[4];
538 
539  // set control points
540  cntrl_in[0] << 2.0, 2.0, 0.0;
541  cntrl_in[1] << 1.0, 1.5, 0.0;
542  cntrl_in[2] << 3.5, 0.0, 0.0;
543  cntrl_in[3] << 4.0, 1.0, 0.0;
544 
545  bezier_type bc1(3);
546 
547  // set control points
548  for (index_type i=0; i<4; ++i)
549  {
550  bc1.set_control_point(cntrl_in[i], i);
551  }
552 
553  // reverse curve parameterization
554  bc1.reverse();
555 
556  // get control points and check them
557  for (index_type i=0; i<4; ++i)
558  {
559  TEST_ASSERT(bc1.get_control_point(i)==cntrl_in[3-i]);
560  }
561  }
562 
564  {
565  point_type cntrl_in[4], cntrl_out[4];
566  index_type i;
567 
568  // set control points
569  cntrl_in[0] << 2.0, 2.0, 1.0;
570  cntrl_in[1] << 1.0, 1.5, 2.0;
571  cntrl_in[2] << 3.5, 0.0, 3.0;
572  cntrl_in[3] << 4.0, 1.0, 4.0;
573 
574  bezier_type bc(3), bcc;
575 
576  // set control points
577  for (index_type i=0; i<4; ++i)
578  {
579  bc.set_control_point(cntrl_in[i], i);
580  }
581  bcc=bc;
582 
583  // reflect about xy-plane
584  {
585  bc.reflect_xy();
586  for (i=0; i<4; ++i)
587  {
588  cntrl_out[i].x()= cntrl_in[i].x();
589  cntrl_out[i].y()= cntrl_in[i].y();
590  cntrl_out[i].z()=-cntrl_in[i].z();
591  }
592 
593  // check control points
594  for (i=0; i<4; ++i)
595  {
596  TEST_ASSERT(bc.get_control_point(i)==cntrl_out[i]);
597  }
598  }
599 
600  // reflect about xz-plane
601  bc=bcc;
602  {
603  bc.reflect_xz();
604  for (i=0; i<4; ++i)
605  {
606  cntrl_out[i].x()= cntrl_in[i].x();
607  cntrl_out[i].y()=-cntrl_in[i].y();
608  cntrl_out[i].z()= cntrl_in[i].z();
609  }
610 
611  // check control points
612  for (i=0; i<4; ++i)
613  {
614  TEST_ASSERT(bc.get_control_point(i)==cntrl_out[i]);
615  }
616  }
617 
618  // reflect about yz-plane
619  bc=bcc;
620  {
621  bc.reflect_yz();
622  for (i=0; i<4; ++i)
623  {
624  cntrl_out[i].x()=-cntrl_in[i].x();
625  cntrl_out[i].y()= cntrl_in[i].y();
626  cntrl_out[i].z()= cntrl_in[i].z();
627  }
628 
629  // check control points
630  for (i=0; i<4; ++i)
631  {
632  TEST_ASSERT(bc.get_control_point(i)==cntrl_out[i]);
633  }
634  }
635 
636  // reflect about arbitrary plane
637  bc=bcc;
638  {
639  data_type eps(std::numeric_limits<data__>::epsilon());
640  point_type normal, n;
641 
642  normal << 1, 2, 3;
643  n=normal;
644  n.normalize();
645 
646  bc.reflect(normal);
647  for (i=0; i<4; ++i)
648  {
649  cntrl_out[i]=cntrl_in[i]-2*cntrl_in[i].dot(n)*n;
650  }
651 
652  // check control points
653  for (i=0; i<4; ++i)
654  {
655  TEST_ASSERT((bc.get_control_point(i)-cntrl_out[i]).norm()<10*eps);
656  }
657  }
658 
659  // reflect about arbitrary plane not through origin
660  bc=bcc;
661  {
662  data_type eps(std::numeric_limits<data__>::epsilon());
663  point_type normal, n;
664  data_type d(4);
665 
666  normal << 1, 2, 3;
667  n=normal;
668  n.normalize();
669 
670  bc.reflect(normal, d);
671  for (i=0; i<4; ++i)
672  {
673  cntrl_out[i]=cntrl_in[i]-2*(cntrl_in[i].dot(n)-d)*n;
674  }
675 
676  // check control points
677  for (i=0; i<4; ++i)
678  {
679  TEST_ASSERT((bc.get_control_point(i)-cntrl_out[i]).norm()<10*eps);
680  }
681  }
682  }
683 
685  {
686  point_type cntrl_in[5];
687  data_type eps(std::numeric_limits<data__>::epsilon());
688 
689  // set control points
690  cntrl_in[0] << 0, 0, 0;
691  cntrl_in[1] << 0, 4, 0;
692  cntrl_in[2] << 2, 4, 0;
693  cntrl_in[3] << 2, 3, 0;
694  cntrl_in[4] << 1.5, 3, 0;
695 
696  bezier_type bc1(4), bc2;
697  point_type eval_out, eval_ref;
698  data_type t, curv_out, curv_ref;
699 
700  // set control points
701  for (index_type i=0; i<5; ++i)
702  {
703  bc1.set_control_point(cntrl_in[i], i);
704  }
705  bc2=bc1;
706 
707  // promote curve
708  bc2.degree_promote();
709 
710  // test to see if degree has increased
711  TEST_ASSERT(bc1.degree()+1==bc2.degree());
712 
713  // test to see if get expected polygon
714  if (bc1.degree()+1==bc2.degree())
715  {
716  point_type cntrl_ref[6];
717 
718  cntrl_ref[0] << 0.0, 0.0, 0.0;
719  cntrl_ref[1] << 0.0, static_cast<data__>(3.2), 0.0;
720  cntrl_ref[2] << static_cast<data__>(1.2), 4.0, 0.0;
721  cntrl_ref[3] << 2.0, static_cast<data__>(3.6), 0.0;
722  cntrl_ref[4] << static_cast<data__>(1.9), 3.0, 0.0;
723  cntrl_ref[5] << static_cast<data__>(1.5), 3.0, 0.0;
724  for (index_type i=0; i<6; ++i)
725  {
726  if (typeid(data__)==typeid(long double))
727  {
728  TEST_ASSERT((bc2.get_control_point(i)-cntrl_ref[i]).norm()<std::sqrt(eps));
729  }
730  else
731  {
732  TEST_ASSERT(bc2.get_control_point(i)==cntrl_ref[i]);
733  }
734  }
735  }
736 
737  // test evaluation at end points
738  t=0;
739  eval_out=bc2.f(t);
740  eval_ref=bc1.f(t);
741  TEST_ASSERT(eval_out==eval_ref);
742  t=1;
743  eval_out=bc2.f(t);
744  eval_ref=bc1.f(t);
745  TEST_ASSERT(eval_out==eval_ref);
746 
747  // test evaluation at interior point
748  t=static_cast<data__>(0.45);
749  eval_out=bc2.f(t);
750  eval_ref=bc1.f(t);
751  TEST_ASSERT((eval_out-eval_ref).norm()<5*eps);
752 
753  // test 1st derivative at end points
754  t=0;
755  eval_out=bc2.fp(t);
756  eval_ref=bc1.fp(t);
757  TEST_ASSERT(eval_out==eval_ref);
758  t=1;
759  eval_out=bc2.fp(t);
760  eval_ref=bc1.fp(t);
761  TEST_ASSERT((eval_out-eval_ref).norm()<5*eps);
762 
763  // test 1st derivative at interior point
764  t=static_cast<data__>(0.45);
765  eval_out=bc2.fp(t);
766  eval_ref=bc1.fp(t);
767  TEST_ASSERT((eval_out-eval_ref).norm()<5*eps);
768 
769  // test 2nd derivative at end points
770  t=0;
771  eval_out=bc2.fpp(t);
772  eval_ref=bc1.fpp(t);
773  TEST_ASSERT((eval_out-eval_ref).norm()<33*eps);
774  t=1;
775  eval_out=bc2.fpp(t);
776  eval_ref=bc1.fpp(t);
777  TEST_ASSERT((eval_out-eval_ref).norm()<18*eps);
778 
779  // test 2nd derivative at interior point
780  t=static_cast<data__>(0.45);
781  eval_out=bc2.fpp(t);
782  eval_ref=bc1.fpp(t);
783  TEST_ASSERT((eval_out-eval_ref).norm()<33*eps);
784 
785  // test 3rd derivative at end points
786  t=0;
787  eval_out=bc2.fppp(t);
788  eval_ref=bc1.fppp(t);
789  TEST_ASSERT((eval_out-eval_ref).norm()<390*eps);
790  t=1;
791  eval_out=bc2.fppp(t);
792  eval_ref=bc1.fppp(t);
793  TEST_ASSERT((eval_out-eval_ref).norm()<390*eps);
794 
795  // test 3rd derivative at interior point
796  t=static_cast<data__>(0.45);
797  eval_out=bc2.fppp(t);
798  eval_ref=bc1.fppp(t);
799  TEST_ASSERT((eval_out-eval_ref).norm()<203*eps);
800 
801  // test curvature at end points
802  t=0;
803  eli::geom::curve::curvature(curv_out, bc2, t);
804  eli::geom::curve::curvature(curv_ref, bc1, t);
805  TEST_ASSERT(curv_out==curv_ref);
806  t=1;
807  eli::geom::curve::curvature(curv_out, bc2, t);
808  eli::geom::curve::curvature(curv_ref, bc1, t);
809  TEST_ASSERT(std::abs(curv_out-curv_ref)<203*eps);
810 
811  // test curvature at interior point
812  t=static_cast<data__>(0.45);
813  eli::geom::curve::curvature(curv_out, bc2, t);
814  eli::geom::curve::curvature(curv_ref, bc1, t);
815  TEST_ASSERT(std::abs(curv_out-curv_ref)<203*eps);
816  }
817 
819  {
820  point_type cntrl_in[5];
821  data_type eps(std::numeric_limits<data__>::epsilon());
822 
823  // set control points
824  cntrl_in[0] << 0, 0, 0;
825  cntrl_in[1] << 0, 4, 0;
826  cntrl_in[2] << 2, 4, 0;
827  cntrl_in[3] << 2, 3, 0;
828  cntrl_in[4] << 1.5, 3, 0;
829 
830  bezier_type bc1(4), bc2, bc3, bc4;
831  point_type eval_out, eval_ref;
832  data_type t, curv_out, curv_ref;
833 
834  // set control points
835  for (index_type i=0; i<5; ++i)
836  {
837  bc1.set_control_point(cntrl_in[i], i);
838  }
839  bc2=bc1;
840  bc3=bc1;
841  bc4=bc1;
842 
843  // promote curve single degree
844  bc2.degree_promote();
845  // promote_to curve
846  bc3.degree_promote_to(bc1.degree()+1);
847 
848  // test to see if degree has increased
849  TEST_ASSERT(bc3.degree()==bc2.degree());
850 
851  // test to see if get exact behavior of single-degree promotion
852  if (bc3.degree()==bc2.degree())
853  {
854  for (index_type i=0; i<6; ++i)
855  {
856  TEST_ASSERT(bc3.get_control_point(i)==bc2.get_control_point(i));
857  }
858  }
859 
860  // promote_to curve high order
861  bc4.degree_promote_to(bc1.degree()+5);
862 
863  // test to see if degree has increased
864  TEST_ASSERT(bc4.degree()==bc1.degree()+5);
865 
866  // test evaluation at end points
867  t=0;
868  eval_out=bc4.f(t);
869  eval_ref=bc1.f(t);
870  TEST_ASSERT(eval_out==eval_ref);
871  t=1;
872  eval_out=bc4.f(t);
873  eval_ref=bc1.f(t);
874  TEST_ASSERT(eval_out==eval_ref);
875 
876  // test evaluation at interior point
877  t=static_cast<data__>(0.45);
878  eval_out=bc4.f(t);
879  eval_ref=bc1.f(t);
880  TEST_ASSERT((eval_out-eval_ref).norm()<5*eps);
881 
882  // test 1st derivative at end points
883  t=0;
884  eval_out=bc4.fp(t);
885  eval_ref=bc1.fp(t);
886  TEST_ASSERT((eval_out-eval_ref).norm()<20*eps);
887  t=1;
888  eval_out=bc4.fp(t);
889  eval_ref=bc1.fp(t);
890  TEST_ASSERT((eval_out-eval_ref).norm()<5*4*eps);
891 
892  // test 1st derivative at interior point
893  t=static_cast<data__>(0.45);
894  eval_out=bc4.fp(t);
895  eval_ref=bc1.fp(t);
896  TEST_ASSERT((eval_out-eval_ref).norm()<5*eps);
897 
898  // test 2nd derivative at end points
899  t=0;
900  eval_out=bc4.fpp(t);
901  eval_ref=bc1.fpp(t);
902  TEST_ASSERT((eval_out-eval_ref).norm()<200*eps);
903  t=1;
904  eval_out=bc4.fpp(t);
905  eval_ref=bc1.fpp(t);
906  TEST_ASSERT((eval_out-eval_ref).norm()<18*16*eps);
907 
908  // test 2nd derivative at interior point
909  t=static_cast<data__>(0.45);
910  eval_out=bc4.fpp(t);
911  eval_ref=bc1.fpp(t);
912  TEST_ASSERT((eval_out-eval_ref).norm()<35*eps);
913 
914  // test 3rd derivative at end points
915  t=0;
916  eval_out=bc4.fppp(t);
917  eval_ref=bc1.fppp(t);
918  TEST_ASSERT((eval_out-eval_ref).norm()<390*16*eps);
919  t=1;
920  eval_out=bc4.fppp(t);
921  eval_ref=bc1.fppp(t);
922  TEST_ASSERT((eval_out-eval_ref).norm()<390*4*eps);
923 
924  // test 3rd derivative at interior point
925  t=static_cast<data__>(0.45);
926  eval_out=bc4.fppp(t);
927  eval_ref=bc1.fppp(t);
928  TEST_ASSERT((eval_out-eval_ref).norm()<203*4*eps);
929 
930  // test curvature at end points
931  t=0;
932  eli::geom::curve::curvature(curv_out, bc4, t);
933  eli::geom::curve::curvature(curv_ref, bc1, t);
934  TEST_ASSERT(std::abs(curv_out-curv_ref)<5*eps);
935  t=1;
936  eli::geom::curve::curvature(curv_out, bc4, t);
937  eli::geom::curve::curvature(curv_ref, bc1, t);
938  TEST_ASSERT(std::abs(curv_out-curv_ref)<203*eps);
939 
940  // test curvature at interior point
941  t=static_cast<data__>(0.45);
942  eli::geom::curve::curvature(curv_out, bc4, t);
943  eli::geom::curve::curvature(curv_ref, bc1, t);
944  TEST_ASSERT(std::abs(curv_out-curv_ref)<203*eps);
945  }
946 
948  {
949  data_type eps(std::numeric_limits<data__>::epsilon());
950 
951  {
952  // hand checked with "Degree Reduction of Bezier Curves" by Dave Morgan
953  point_type cntrl_in[7], cntrl_ref[6];
954 
955  cntrl_in[0] << 0, 0, 0;
956  cntrl_in[1] << 2, 6, 0;
957  cntrl_in[2] << 3, 0, 0;
958  cntrl_in[3] << 5, 4, 0;
959  cntrl_in[4] << 7, 1, 0;
960  cntrl_in[5] << 5, 5, 0;
961  cntrl_in[6] << 10, 6, 0;
962 
963  bezier_type bc1(6), bc2;
964 
965  // set control points
966  for (index_type i=0; i<7; ++i)
967  {
968  bc1.set_control_point(cntrl_in[i], i);
969  }
970 
971  bc2=bc1;
973 
974  // test if the degree is correct
975  TEST_ASSERT(bc1.degree()+1==bc2.degree());
976 
977  if (bc1.degree()+1==bc2.degree())
978  {
979  // test of get the correct control points
980  cntrl_ref[0] << static_cast<data__>(-0.00878906), static_cast<data__>(0.0610352), 0;
981  cntrl_ref[1] << static_cast<data__>( 2.5177734), static_cast<data__>(6.3821289), 0;
982  cntrl_ref[2] << static_cast<data__>( 2.8060547), static_cast<data__>(-0.16982422), 0;
983  cntrl_ref[3] << static_cast<data__>( 8.0060547), static_cast<data__>( 2.5301758), 0;
984  cntrl_ref[4] << static_cast<data__>( 4.1177734), static_cast<data__>( 3.9821289), 0;
985  cntrl_ref[5] << static_cast<data__>( 9.9912109), static_cast<data__>( 6.0610352), 0;
986  for (index_type i=0; i<6; ++i)
987  {
988  TEST_ASSERT_DELTA((bc1.get_control_point(i)-cntrl_ref[i]).norm(), 0, 1e-6);
989  }
990  }
991  }
992 
993  // test whether can promote and demote and get the same control points
994  {
995  point_type cntrl_in[5];
996 
997  // set control points
998  cntrl_in[0] << 0, 0, 0;
999  cntrl_in[1] << 0, 4, 0;
1000  cntrl_in[2] << 2, 4, 0;
1001  cntrl_in[3] << 2, 3, 0;
1002  cntrl_in[4] << 1.5, 3, 0;
1003 
1004  bezier_type bc1(4), bc2;
1005 
1006  // set control points
1007  for (index_type i=0; i<5; ++i)
1008  {
1009  bc1.set_control_point(cntrl_in[i], i);
1010  }
1011  bc2=bc1;
1012  bc2.degree_promote();
1014 
1015  // test to see if degree has changed
1016  TEST_ASSERT(bc1.degree()==bc2.degree());
1017 
1018  // test to see if get expected polygon
1019  for (index_type i=0; i<5; ++i)
1020  {
1021  if (typeid(data__)==typeid(long double))
1022  {
1023  TEST_ASSERT((bc2.get_control_point(i)-cntrl_in[i]).norm()<std::sqrt(eps));
1024  }
1025  else
1026  {
1027  TEST_ASSERT((bc2.get_control_point(i)-cntrl_in[i]).norm()<4*eps);
1028  }
1029  }
1030  }
1031 
1032  // test whether can demote and maintain continuity at end points
1033  {
1034  point_type cntrl_in[10];
1035 
1036  cntrl_in[0] << 0, 0, 0,
1037  cntrl_in[1] << 2, 6, 0,
1038  cntrl_in[2] << 3, 0, 0,
1039  cntrl_in[3] << 5, 4, 0,
1040  cntrl_in[4] << 7, 1, 0,
1041  cntrl_in[5] << 5, 5, 0,
1042  cntrl_in[6] << 6, 6, 0,
1043  cntrl_in[7] << 7, 5, 0,
1044  cntrl_in[8] << 5, 5, 0,
1045  cntrl_in[9] << 10, 6, 0;
1046 
1047  bezier_type bc1(9);
1048 
1049  // set control points
1050  for (index_type i=0; i<10; ++i)
1051  {
1052  bc1.set_control_point(cntrl_in[i], i);
1053  }
1054 
1055  // No End Constraint
1056  {
1057  bezier_type bc2(bc1);
1058  point_type eval_out[3], eval_ref[3];
1059  data_type d[3], tv[3]={0, 1, static_cast<data__>(0.45)};
1060  bool success;
1061 
1063  TEST_ASSERT(success);
1064 
1065  for (size_t k=0; k<3; ++k)
1066  {
1067  data_type t(tv[k]);
1068  eval_out[0]=bc2.f(t);
1069  eval_out[1]=bc2.fp(t);
1070  eval_out[2]=bc2.fpp(t);
1071  eval_ref[0]=bc1.f(t);
1072  eval_ref[1]=bc1.fp(t);
1073  eval_ref[2]=bc1.fpp(t);
1074  d[0]=eli::geom::point::distance(eval_out[0], eval_ref[0]);
1075  d[1]=eli::geom::point::distance(eval_out[1], eval_ref[1]);
1076  d[2]=eli::geom::point::distance(eval_out[2], eval_ref[2]);
1077 
1078  if (k<2)
1079  {
1080  TEST_ASSERT_DELTA(d[0], 0.00435372, 5e-4);
1081  TEST_ASSERT_DELTA(d[1], 0.70530, 5e-3);
1082  TEST_ASSERT_DELTA(d[2], 37.6161, 5e-2);
1083  }
1084  else
1085  {
1086  TEST_ASSERT_DELTA(d[0], 0.003414, 1e-5);
1087  TEST_ASSERT_DELTA(d[1], 0.04887, 5e-4);
1088  TEST_ASSERT_DELTA(d[2], 1.10759, 5e-4);
1089  }
1090  }
1091  }
1092 
1093  // Point End Constraint
1094  {
1095  bezier_type bc2(bc1);
1096  point_type eval_out[3], eval_ref[3];
1097  data_type d[3], tv[3]={0, 1, static_cast<data__>(0.45)};
1098  bool success;
1099 
1101  TEST_ASSERT(success);
1102 
1103  for (size_t k=0; k<3; ++k)
1104  {
1105  data_type t(tv[k]);
1106  eval_out[0]=bc2.f(t);
1107  eval_out[1]=bc2.fp(t);
1108  eval_out[2]=bc2.fpp(t);
1109  eval_ref[0]=bc1.f(t);
1110  eval_ref[1]=bc1.fp(t);
1111  eval_ref[2]=bc1.fpp(t);
1112  d[0]=eli::geom::point::distance(eval_out[0], eval_ref[0]);
1113  d[1]=eli::geom::point::distance(eval_out[1], eval_ref[1]);
1114  d[2]=eli::geom::point::distance(eval_out[2], eval_ref[2]);
1115 
1116  if (k<2)
1117  {
1118  TEST_ASSERT_DELTA(d[0], 0.0, 5e-6);
1119  TEST_ASSERT_DELTA(d[1], 0.64553, 5e-5);
1120  TEST_ASSERT_DELTA(d[2], 37.4409, 5e-4);
1121  }
1122  else
1123  {
1124  TEST_ASSERT_DELTA(d[0], 0.00305362, 5e-6);
1125  TEST_ASSERT_DELTA(d[1], 0.042752, 5e-5);
1126  TEST_ASSERT_DELTA(d[2], 1.043406, 1e-4);
1127  }
1128  }
1129  }
1130 
1131  // Slope End Constraint
1132  {
1133  bezier_type bc2(bc1);
1134  point_type eval_out[3], eval_ref[3];
1135  data_type d[3], tv[3]={0, 1, static_cast<data__>(0.45)};
1136  bool success;
1137 
1139  TEST_ASSERT(success);
1140 
1141  for (size_t k=0; k<3; ++k)
1142  {
1143  data_type t(tv[k]);
1144  eval_out[0]=bc2.f(t);
1145  eval_out[1]=bc2.fp(t);
1146  eval_out[2]=bc2.fpp(t);
1147  eval_ref[0]=bc1.f(t);
1148  eval_ref[1]=bc1.fp(t);
1149  eval_ref[2]=bc1.fpp(t);
1150  d[0]=eli::geom::point::distance(eval_out[0], eval_ref[0]);
1151  d[1]=eli::geom::point::distance(eval_out[1], eval_ref[1]);
1152  d[2]=eli::geom::point::distance(eval_out[2], eval_ref[2]);
1153 
1154  if (k<2)
1155  {
1156  TEST_ASSERT(d[0]==0);
1157  TEST_ASSERT(d[1]==0);
1158  TEST_ASSERT_DELTA(d[2], 16.7838, 5e-4);
1159  }
1160  else
1161  {
1162  TEST_ASSERT_DELTA(d[0], 0.0057941, 5e-6);
1163  TEST_ASSERT_DELTA(d[1], 0.086370, 5e-5);
1164  TEST_ASSERT_DELTA(d[2], 1.69369, 5e-4);
1165  }
1166  }
1167  }
1168 
1169  // Second Derivative Constraint
1170  {
1171  bezier_type bc2(bc1);
1172  point_type eval_out[3], eval_ref[3];
1173  data_type d[3], tv[3]={0, 1, static_cast<data__>(0.45)};
1174  bool success;
1175 
1177  TEST_ASSERT(success);
1178 
1179  for (size_t k=0; k<3; ++k)
1180  {
1181  data_type t(tv[k]);
1182  eval_out[0]=bc2.f(t);
1183  eval_out[1]=bc2.fp(t);
1184  eval_out[2]=bc2.fpp(t);
1185  eval_ref[0]=bc1.f(t);
1186  eval_ref[1]=bc1.fp(t);
1187  eval_ref[2]=bc1.fpp(t);
1188  d[0]=eli::geom::point::distance(eval_out[0], eval_ref[0]);
1189  d[1]=eli::geom::point::distance(eval_out[1], eval_ref[1]);
1190  d[2]=eli::geom::point::distance(eval_out[2], eval_ref[2]);
1191 
1192  if (k<2)
1193  {
1194  TEST_ASSERT(d[0]==0);
1195  TEST_ASSERT(d[1]==0);
1196  TEST_ASSERT(d[2] < 129*eps);
1197  }
1198  else
1199  {
1200  TEST_ASSERT_DELTA(d[0], 0.0180029, 5e-5);
1201  TEST_ASSERT_DELTA(d[1], 0.294979, 5e-5);
1202  TEST_ASSERT_DELTA(d[2], 3.78229, 5e-4);
1203  }
1204  }
1205  }
1206  }
1207  }
1208 
1209  void degree_to_cubic_test() // degree_to_cubic();
1210  {
1211  data_type eps(std::numeric_limits<data__>::epsilon());
1212 
1213  // Test linear promotion case.
1214  {
1215  point_type cntrl_in[2];
1216 
1217  cntrl_in[0] << 0, 0, 0;
1218  cntrl_in[1] << 2, 6, 0;
1219 
1220  bezier_type bc1(1), bc2;
1221  point_type eval_out, eval_ref;
1222  data_type t, curv_out, curv_ref;
1223 
1224  // set control points
1225  for (index_type i=0; i<2; ++i)
1226  {
1227  bc1.set_control_point(cntrl_in[i], i);
1228  }
1229 
1230  bc2=bc1;
1231  bc1.degree_to_cubic();
1232 
1233  // test if the degree is correct
1234  TEST_ASSERT(bc1.degree()==3);
1235 
1236  // test evaluation at end points
1237  t=0;
1238  eval_out=bc2.f(t);
1239  eval_ref=bc1.f(t);
1240  TEST_ASSERT(eval_out==eval_ref);
1241  t=1;
1242  eval_out=bc2.f(t);
1243  eval_ref=bc1.f(t);
1244  TEST_ASSERT(eval_out==eval_ref);
1245 
1246  // test evaluation at interior point
1247  t=static_cast<data__>(0.45);
1248  eval_out=bc2.f(t);
1249  eval_ref=bc1.f(t);
1250  TEST_ASSERT((eval_out-eval_ref).norm()<5*eps);
1251 
1252  // test 1st derivative at end points
1253  t=0;
1254  eval_out=bc2.fp(t);
1255  eval_ref=bc1.fp(t);
1256  TEST_ASSERT((eval_out-eval_ref).norm()<5*eps);
1257  t=1;
1258  eval_out=bc2.fp(t);
1259  eval_ref=bc1.fp(t);
1260  TEST_ASSERT((eval_out-eval_ref).norm()<5*eps);
1261 
1262  // test 1st derivative at interior point
1263  t=static_cast<data__>(0.45);
1264  eval_out=bc2.fp(t);
1265  eval_ref=bc1.fp(t);
1266  TEST_ASSERT((eval_out-eval_ref).norm()<5*eps);
1267 
1268  // test 2nd derivative at end points
1269  t=0;
1270  eval_out=bc2.fpp(t);
1271  eval_ref=bc1.fpp(t);
1272  TEST_ASSERT((eval_out-eval_ref).norm()<33*eps);
1273  t=1;
1274  eval_out=bc2.fpp(t);
1275  eval_ref=bc1.fpp(t);
1276  TEST_ASSERT((eval_out-eval_ref).norm()<18*eps);
1277 
1278  // test 2nd derivative at interior point
1279  t=static_cast<data__>(0.45);
1280  eval_out=bc2.fpp(t);
1281  eval_ref=bc1.fpp(t);
1282  TEST_ASSERT((eval_out-eval_ref).norm()<33*eps);
1283 
1284  // test 3rd derivative at end points
1285  t=0;
1286  eval_out=bc2.fppp(t);
1287  eval_ref=bc1.fppp(t);
1288  TEST_ASSERT((eval_out-eval_ref).norm()<390*eps);
1289  t=1;
1290  eval_out=bc2.fppp(t);
1291  eval_ref=bc1.fppp(t);
1292  TEST_ASSERT((eval_out-eval_ref).norm()<390*eps);
1293 
1294  // test 3rd derivative at interior point
1295  t=static_cast<data__>(0.45);
1296  eval_out=bc2.fppp(t);
1297  eval_ref=bc1.fppp(t);
1298  TEST_ASSERT((eval_out-eval_ref).norm()<203*eps);
1299 
1300  // test curvature at end points
1301  t=0;
1302  eli::geom::curve::curvature(curv_out, bc2, t);
1303  eli::geom::curve::curvature(curv_ref, bc1, t);
1304  TEST_ASSERT(curv_out==curv_ref);
1305  t=1;
1306  eli::geom::curve::curvature(curv_out, bc2, t);
1307  eli::geom::curve::curvature(curv_ref, bc1, t);
1308  TEST_ASSERT(std::abs(curv_out-curv_ref)<203*eps);
1309 
1310  // test curvature at interior point
1311  t=static_cast<data__>(0.45);
1312  eli::geom::curve::curvature(curv_out, bc2, t);
1313  eli::geom::curve::curvature(curv_ref, bc1, t);
1314  TEST_ASSERT(std::abs(curv_out-curv_ref)<203*eps);
1315  }
1316 
1317  // Test quadratic promotion case.
1318  {
1319  point_type cntrl_in[3];
1320 
1321  cntrl_in[0] << 0, 0, 0;
1322  cntrl_in[1] << 2, 6, 0;
1323  cntrl_in[2] << 3, 0, 0;
1324 
1325  bezier_type bc1(2), bc2;
1326  point_type eval_out, eval_ref;
1327  data_type t, curv_out, curv_ref;
1328 
1329  // set control points
1330  for (index_type i=0; i<3; ++i)
1331  {
1332  bc1.set_control_point(cntrl_in[i], i);
1333  }
1334 
1335  bc2=bc1;
1336  bc1.degree_to_cubic();
1337 
1338  // test if the degree is correct
1339  TEST_ASSERT(bc1.degree()==3);
1340 
1341  // test evaluation at end points
1342  t=0;
1343  eval_out=bc2.f(t);
1344  eval_ref=bc1.f(t);
1345  TEST_ASSERT(eval_out==eval_ref);
1346  t=1;
1347  eval_out=bc2.f(t);
1348  eval_ref=bc1.f(t);
1349  TEST_ASSERT(eval_out==eval_ref);
1350 
1351  // test evaluation at interior point
1352  t=static_cast<data__>(0.45);
1353  eval_out=bc2.f(t);
1354  eval_ref=bc1.f(t);
1355  TEST_ASSERT((eval_out-eval_ref).norm()<5*eps);
1356 
1357  // test 1st derivative at end points
1358  t=0;
1359  eval_out=bc2.fp(t);
1360  eval_ref=bc1.fp(t);
1361  TEST_ASSERT((eval_out-eval_ref).norm()<5*eps);
1362  t=1;
1363  eval_out=bc2.fp(t);
1364  eval_ref=bc1.fp(t);
1365  TEST_ASSERT((eval_out-eval_ref).norm()<5*eps);
1366 
1367  // test 1st derivative at interior point
1368  t=static_cast<data__>(0.45);
1369  eval_out=bc2.fp(t);
1370  eval_ref=bc1.fp(t);
1371  TEST_ASSERT((eval_out-eval_ref).norm()<5*eps);
1372 
1373  // test 2nd derivative at end points
1374  t=0;
1375  eval_out=bc2.fpp(t);
1376  eval_ref=bc1.fpp(t);
1377  TEST_ASSERT((eval_out-eval_ref).norm()<33*eps);
1378  t=1;
1379  eval_out=bc2.fpp(t);
1380  eval_ref=bc1.fpp(t);
1381  TEST_ASSERT((eval_out-eval_ref).norm()<18*eps);
1382 
1383  // test 2nd derivative at interior point
1384  t=static_cast<data__>(0.45);
1385  eval_out=bc2.fpp(t);
1386  eval_ref=bc1.fpp(t);
1387  TEST_ASSERT((eval_out-eval_ref).norm()<33*eps);
1388 
1389  // test 3rd derivative at end points
1390  t=0;
1391  eval_out=bc2.fppp(t);
1392  eval_ref=bc1.fppp(t);
1393  TEST_ASSERT((eval_out-eval_ref).norm()<390*eps);
1394  t=1;
1395  eval_out=bc2.fppp(t);
1396  eval_ref=bc1.fppp(t);
1397  TEST_ASSERT((eval_out-eval_ref).norm()<390*eps);
1398 
1399  // test 3rd derivative at interior point
1400  t=static_cast<data__>(0.45);
1401  eval_out=bc2.fppp(t);
1402  eval_ref=bc1.fppp(t);
1403  TEST_ASSERT((eval_out-eval_ref).norm()<203*eps);
1404 
1405  // test curvature at end points
1406  t=0;
1407  eli::geom::curve::curvature(curv_out, bc2, t);
1408  eli::geom::curve::curvature(curv_ref, bc1, t);
1409  TEST_ASSERT((curv_out-curv_ref)<10*eps);
1410  t=1;
1411  eli::geom::curve::curvature(curv_out, bc2, t);
1412  eli::geom::curve::curvature(curv_ref, bc1, t);
1413  TEST_ASSERT(std::abs(curv_out-curv_ref)<203*eps);
1414 
1415  // test curvature at interior point
1416  t=static_cast<data__>(0.45);
1417  eli::geom::curve::curvature(curv_out, bc2, t);
1418  eli::geom::curve::curvature(curv_ref, bc1, t);
1419  TEST_ASSERT(std::abs(curv_out-curv_ref)<203*eps);
1420  }
1421 
1422  // Test do-nothing cubic case.
1423  {
1424  point_type cntrl_in[4];
1425 
1426  cntrl_in[0] << 0, 0, 0;
1427  cntrl_in[1] << 2, 6, 0;
1428  cntrl_in[2] << 3, 0, 0;
1429  cntrl_in[3] << 5, 4, 0;
1430 
1431  bezier_type bc1(3), bc2;
1432 
1433  // set control points
1434  for (index_type i=0; i<4; ++i)
1435  {
1436  bc1.set_control_point(cntrl_in[i], i);
1437  }
1438 
1439  bc2=bc1;
1440  bc1.degree_to_cubic();
1441 
1442  // test if the degree is correct
1443  TEST_ASSERT(bc1.degree()==3);
1444 
1445  // test to see that nothing happened
1446  if (bc1.degree()==3)
1447  {
1448  for (index_type i=0; i<4; ++i)
1449  {
1450  TEST_ASSERT(bc1.get_control_point(i)==bc2.get_control_point(i));
1451  }
1452  }
1453  }
1454 
1455  // Test demotion case.
1456  {
1457  point_type cntrl_in[7];
1458 
1459  cntrl_in[0] << 0, 0, 0;
1460  cntrl_in[1] << 2, 6, 0;
1461  cntrl_in[2] << 3, 0, 0;
1462  cntrl_in[3] << 5, 4, 0;
1463  cntrl_in[4] << 7, 1, 0;
1464  cntrl_in[5] << 5, 5, 0;
1465  cntrl_in[6] << 10, 6, 0;
1466 
1467  bezier_type bc1(6), bc2;
1468 
1469  // set control points
1470  for (index_type i=0; i<7; ++i)
1471  {
1472  bc1.set_control_point(cntrl_in[i], i);
1473  }
1474 
1475  bc2=bc1;
1476  bc1.degree_to_cubic();
1477 
1478  // test if the degree is correct
1479  TEST_ASSERT(bc1.degree()==3);
1480 
1481  // test that start/end point/derivatives are equal.
1482  data_type t;
1483  t = 0;
1484  TEST_ASSERT(eli::geom::point::distance(bc2.f(t), bc1.f(t))==0);
1485  TEST_ASSERT(eli::geom::point::distance(bc2.fp(t), bc1.fp(t))==0);
1486 
1487  t = 1;
1488  TEST_ASSERT(eli::geom::point::distance(bc2.f(t), bc1.f(t))==0);
1489  TEST_ASSERT(eli::geom::point::distance(bc2.fp(t), bc1.fp(t))==0);
1490  }
1491  }
1492 
1494  {
1495  // Test that curve has zero distance itself
1496  {
1497  point_type cntrl_in[5];
1498 
1499  // set control points
1500  cntrl_in[0] << 0, 0, 0;
1501  cntrl_in[1] << 0, 4, 0;
1502  cntrl_in[2] << 2, 4, 0;
1503  cntrl_in[3] << 2, 3, 0;
1504  cntrl_in[4] << 1.5, 3, 0;
1505 
1506  bezier_type bc1(4);
1507 
1508  // set control points
1509  for (index_type i=0; i<5; ++i)
1510  {
1511  bc1.set_control_point(cntrl_in[i], i);
1512  }
1513 
1514  data_type d1 = bc1.eqp_distance_bound(bc1);
1515  TEST_ASSERT(d1==0);
1516  }
1517 
1518  // Test that curve has zero distance from promoted-self (and vis-versa)
1519  {
1520  point_type cntrl_in[5];
1521 
1522  // set control points
1523  cntrl_in[0] << 0, 0, 0;
1524  cntrl_in[1] << 0, 4, 0;
1525  cntrl_in[2] << 2, 4, 0;
1526  cntrl_in[3] << 2, 3, 0;
1527  cntrl_in[4] << 1.5, 3, 0;
1528 
1529  bezier_type bc1(4), bc2;
1530 
1531  // set control points
1532  for (index_type i=0; i<5; ++i)
1533  {
1534  bc1.set_control_point(cntrl_in[i], i);
1535  }
1536  bc2=bc1;
1537 
1538  // promote_to curve high order
1539  bc2.degree_promote_to(bc1.degree()+5);
1540 
1541  data_type d1 = bc2.eqp_distance_bound(bc1);
1542  TEST_ASSERT(d1==0);
1543 
1544  data_type d2 = bc1.eqp_distance_bound(bc2);
1545  TEST_ASSERT(d2==0);
1546  }
1547 
1548  // Test that curve has known distance from offset & promoted-self (and vis-versa)
1549  {
1550  point_type cntrl_in[5];
1551 
1552  // set control points
1553  cntrl_in[0] << 0, 0, 0;
1554  cntrl_in[1] << 0, 4, 0;
1555  cntrl_in[2] << 2, 4, 0;
1556  cntrl_in[3] << 2, 3, 0;
1557  cntrl_in[4] << 1.5, 3, 0;
1558 
1559  data_type dz(3);
1560  point_type offset;
1561  offset << 0, 0, dz;
1562 
1563  bezier_type bc1(4), bc2(4);
1564 
1565  // set control points
1566  for (index_type i=0; i<5; ++i)
1567  {
1568  bc1.set_control_point(cntrl_in[i], i);
1569  bc2.set_control_point(cntrl_in[i]+offset, i);
1570  }
1571 
1572  // promote_to curve high order
1573  bc2.degree_promote_to(bc1.degree()+5);
1574 
1575  data_type d1 = bc2.eqp_distance_bound(bc1);
1576  TEST_ASSERT(d1==dz);
1577 
1578  data_type d2 = bc1.eqp_distance_bound(bc2);
1579  TEST_ASSERT(d2==dz);
1580  }
1581 
1582  // Test that curve has known distance from reversed & promoted-self (and vis-versa)
1583  {
1584  point_type cntrl_in[5];
1585 
1586  // set control points
1587  // curve constructed such that greatest distance from reverse will be the endpoints.
1588  cntrl_in[0] << 0, 0, 0;
1589  cntrl_in[1] << 1, 1, 0;
1590  cntrl_in[2] << 2, 0.5, 0;
1591  cntrl_in[3] << 3, 2.5, 0;
1592  cntrl_in[4] << 4, 6.0, 0;
1593 
1594  bezier_type bc1(4), bc2(4);
1595 
1596  // set control points
1597  for (index_type i=0; i<5; ++i)
1598  {
1599  bc1.set_control_point(cntrl_in[i], i);
1600  }
1601  bc2=bc1;
1602 
1603  bc2.reverse();
1604 
1605  data_type dendpts = (bc1.get_control_point(4)-bc1.get_control_point(0) ).norm();
1606 
1607  // promote_to curve high order
1608  bc2.degree_promote_to(bc1.degree()+5);
1609 
1610  data_type d1 = bc2.eqp_distance_bound(bc1);
1611  TEST_ASSERT(d1==dendpts);
1612 
1613  data_type d2 = bc1.eqp_distance_bound(bc2);
1614  TEST_ASSERT(d2==dendpts);
1615  }
1616  }
1617 
1618  void split_test()
1619  {
1620  point_type cntrl_in[4], cntrl_ref[4];
1621  data_type eps(std::numeric_limits<data__>::epsilon());
1622 
1623  // set control points
1624  cntrl_in[0] << 0, 0, 0;
1625  cntrl_in[1] << 0, 2, 0;
1626  cntrl_in[2] << 8, 2, 0;
1627  cntrl_in[3] << 4, 0, 0;
1628 
1629  bezier_type bc1(3), bc1l, bc1r;
1630  point_type eval_out, eval_ref;
1631  data_type t;
1632 
1633  // set control points
1634  for (index_type i=0; i<4; ++i)
1635  {
1636  bc1.set_control_point(cntrl_in[i], i);
1637  }
1638 
1639  // test split with known control points
1640  t=0.5;
1641  bc1.split(bc1l, bc1r, t);
1642  cntrl_ref[0] << 0.0, 0.0, 0.0;
1643  cntrl_ref[1] << 0.0, 1.0, 0.0;
1644  cntrl_ref[2] << 2.0, 1.5, 0.0;
1645  cntrl_ref[3] << 3.5, 1.5, 0.0;
1646  for (index_type i=0; i<4; ++i)
1647  {
1648  TEST_ASSERT(bc1l.get_control_point(i)==cntrl_ref[i]);
1649  }
1650  cntrl_ref[0] << 3.5, 1.5, 0.0;
1651  cntrl_ref[1] << 5.0, 1.5, 0.0;
1652  cntrl_ref[2] << 6.0, 1.0, 0.0;
1653  cntrl_ref[3] << 4.0, 0.0, 0.0;
1654  for (index_type i=0; i<4; ++i)
1655  {
1656  TEST_ASSERT(bc1r.get_control_point(i)==cntrl_ref[i]);
1657  }
1658 
1659  // split the curve and check the evaluations
1660  data_type tl, tr, ts;
1661  tl=static_cast<data__>(0.3);
1662  tr=static_cast<data__>(0.87);
1663  ts=static_cast<data__>(0.586);
1664 
1665  bc1.split(bc1l, bc1r, ts);
1666 
1667  // check the left curve
1668  t=tl*ts;
1669  eval_out=bc1l.f(tl);
1670  eval_ref=bc1.f(t);
1671  TEST_ASSERT((eval_out-eval_ref).norm()<1e2*eps);
1672  eval_out=bc1l.fp(tl);
1673  eval_ref=bc1.fp(t)*ts;
1674  TEST_ASSERT((eval_out-eval_ref).norm()<1e2*eps);
1675  eval_out=bc1l.fpp(tl);
1676  eval_ref=bc1.fpp(t)*ts*ts;
1677  TEST_ASSERT((eval_out-eval_ref).norm()<1e2*eps);
1678  eval_out=bc1l.fppp(tl);
1679  eval_ref=bc1.fppp(t)*ts*ts*ts;
1680  TEST_ASSERT((eval_out-eval_ref).norm()<1e2*eps);
1681 
1682  // check the right curve
1683  t=ts+tr*(1-ts);
1684  eval_out=bc1r.f(tr);
1685  eval_ref=bc1.f(t);
1686  TEST_ASSERT((eval_out-eval_ref).norm()<1e2*eps);
1687  eval_out=bc1r.fp(tr);
1688  eval_ref=bc1.fp(t)*(1-ts);
1689  TEST_ASSERT((eval_out-eval_ref).norm()<1e2*eps);
1690  eval_out=bc1r.fpp(tr);
1691  eval_ref=bc1.fpp(t)*(1-ts)*(1-ts);
1692  TEST_ASSERT((eval_out-eval_ref).norm()<1e2*eps);
1693  eval_out=bc1r.fppp(tr);
1694  eval_ref=bc1.fppp(t)*(1-ts)*(1-ts)*(1-ts);
1695  TEST_ASSERT((eval_out-eval_ref).norm()<1.21e2*eps);
1696  }
1697 
1699  {
1700  point_type cntrl_in[4];
1701  data_type eps(std::numeric_limits<data__>::epsilon());
1702 
1703  // set control points
1704  cntrl_in[0] << 0, 0, 0;
1705  cntrl_in[1] << 0, 2, 0;
1706  cntrl_in[2] << 8, 2, 0;
1707  cntrl_in[3] << 4, 0, 0;
1708 
1709  bezier_type bc1(3);
1710  data_type length_cal, length_ref;
1711 
1712  // set control points
1713  for (index_type i=0; i<4; ++i)
1714  {
1715  bc1.set_control_point(cntrl_in[i], i);
1716  }
1717 
1718  // calculate the reference length using simpson's rule
1719  size_t i, npts(1001), n(npts-1);
1720  std::vector<data__> speed(npts);
1721  for (i=0; i<npts; ++i)
1722  {
1723  data_type t(i/static_cast<data__>(n));
1724  speed[i]=bc1.fp(t).norm();
1725  }
1726  length_ref=speed[0];
1727  for (i=1; i<=n-1; i+=2)
1728  {
1729  length_ref+=4*speed[i];
1730  length_ref+=2*speed[i+1];
1731  }
1732  length_ref-=speed[n];
1733  length_ref*=static_cast<data__>(1-0)/n/3;
1734 
1735  // compute length and compare
1736  data_type tol(std::sqrt(eps));
1737  length(length_cal, bc1, tol);
1738  TEST_ASSERT_DELTA(1, length_cal/length_ref, tol);
1739 
1740  // test computing some segment length
1741  data_type length01_cal, length12_cal, t0, t1, t2;
1742  t0=0;
1743  t1=static_cast<data__>(0.3);
1744  t2=1;
1745 
1746  length(length01_cal, bc1, t0, t1, tol);
1747  length(length12_cal, bc1, t1, t2, tol);
1748  TEST_ASSERT_DELTA(1, (length01_cal+length12_cal)/length_cal, tol);
1749  }
1750 };
1751 
1752 #endif
1753 
Definition: continuity.hpp:26
void evaluation_test()
Definition: bezier_curve_test_suite.hpp:338
void get_bounding_box(bounding_box_type &bb) const
Definition: bezier.hpp:270
void reflect_xy()
Definition: bezier.hpp:225
Eigen::Matrix< data_type, dim__, dim__ > rotation_matrix_type
Definition: bezier.hpp:119
Definition: continuity.hpp:28
bezier_type::point_type point_type
Definition: bezier_curve_test_suite.hpp:38
void promotion_to_test()
Definition: bezier_curve_test_suite.hpp:818
Derived1__::Scalar distance(const Eigen::MatrixBase< Derived1__ > &p1, const Eigen::MatrixBase< Derived2__ > &p2)
Definition: distance.hpp:33
void reflect_xz()
Definition: bezier.hpp:230
Definition: bounding_box.hpp:27
void derivative_1_test()
Definition: bezier_curve_test_suite.hpp:374
void rotate(const rotation_matrix_type &rmat)
Definition: bezier.hpp:281
constraint_info::point_type constraint_point_type
Definition: fit_container.hpp:117
Eigen::Matrix< data_type, 1, dim__ > point_type
Definition: bezier.hpp:114
void create_circle(std::vector< point_type > &pts)
Definition: bezier_curve_test_suite.hpp:146
void bounding_box_test()
Definition: bezier_curve_test_suite.hpp:199
void degree_to_cubic_test()
Definition: bezier_curve_test_suite.hpp:1209
void AddTests(const double &)
Definition: bezier_curve_test_suite.hpp:63
void split_test()
Definition: bezier_curve_test_suite.hpp:1618
void demotion_test()
Definition: bezier_curve_test_suite.hpp:947
point_type get_max() const
Definition: bounding_box.hpp:104
void degree_promote()
Definition: bezier.hpp:436
void derivative_3_test()
Definition: bezier_curve_test_suite.hpp:450
point_type f(const data_type &t) const
Definition: bezier.hpp:324
void derivative_2_test()
Definition: bezier_curve_test_suite.hpp:412
Definition: fit_container.hpp:34
void degree_promote_to(const index_type target_degree)
Definition: bezier.hpp:448
void AddTests(const long double &)
Definition: bezier_curve_test_suite.hpp:84
point_type tangent(const data_type &t) const
Definition: bezier.hpp:419
void translate(const point_type &trans)
Definition: bezier.hpp:293
void frenet_serret_frame(point_type &t, point_type &n, point_type &b, const data_type &t0)
Definition: bezier.hpp:427
Definition: math.hpp:25
bezier_type::fit_container_type fit_container_type
Definition: bezier_curve_test_suite.hpp:35
~bezier_curve_test_suite()
Definition: bezier_curve_test_suite.hpp:111
index_type degree() const
Definition: bezier.hpp:191
point_type fpp(const data_type &t) const
Definition: bezier.hpp:369
control_point_type get_control_point(const index_type &i) const
Definition: bezier.hpp:213
Definition: continuity.hpp:27
fit_container_type::constraint_point_type constraint_point_type
Definition: bezier_curve_test_suite.hpp:36
void promotion_test()
Definition: bezier_curve_test_suite.hpp:684
void degree_to_cubic()
Definition: bezier.hpp:494
void curvature(typename curve__::data_type &rho, const curve__ &c, const typename curve__::data_type &t)
Definition: curvature.hpp:25
void length(typename piecewise< curve__, data__, dim__, tol__ >::data_type &len, const piecewise< curve__, data__, dim__, tol__ > &pc, const typename piecewise< curve__, data__, dim__, tol__ >::data_type &tol)
Definition: length.hpp:43
void frenet_serret_test()
Definition: bezier_curve_test_suite.hpp:505
point_type get_min() const
Definition: bounding_box.hpp:93
bezier_type::data_type data_type
Definition: bezier_curve_test_suite.hpp:39
void assignment_test()
Definition: bezier_curve_test_suite.hpp:160
bezier_curve_test_suite()
Definition: bezier_curve_test_suite.hpp:107
void resize(const index_type &t_dim)
Definition: bezier.hpp:186
void AddTests(const float &)
Definition: bezier_curve_test_suite.hpp:42
void set_control_point(const control_point_type &cp, const index_type &i)
Definition: bezier.hpp:201
void octave_print(int figno, const std::vector< point_type > &pts, const bezier_type &bez) const
Definition: bezier_curve_test_suite.hpp:116
eli::geom::curve::bezier< data__, 3 > bezier_type
Definition: bezier_curve_test_suite.hpp:34
void length_test()
Definition: bezier_curve_test_suite.hpp:1698
void transformation_test()
Definition: bezier_curve_test_suite.hpp:227
point_type fppp(const data_type &t) const
Definition: bezier.hpp:394
point_type fp(const data_type &t) const
Definition: bezier.hpp:344
data__ data_type
Definition: bezier.hpp:113
void reverse_test()
Definition: bezier_curve_test_suite.hpp:535
void distance_bound_test()
Definition: bezier_curve_test_suite.hpp:1493
void reflection_test()
Definition: bezier_curve_test_suite.hpp:563
Definition: bezier.hpp:109
void reverse()
Definition: bezier.hpp:256
data_type eqp_distance_bound(const bezier< data_type, dim__, tolerance_type > &bc) const
Definition: bezier.hpp:174
Definition: continuity.hpp:29
Definition: bezier_curve_test_suite.hpp:31
void reflect_yz()
Definition: bezier.hpp:235
bool degree_demote(const geom::general::continuity &continuity_degree=geom::general::C0)
Definition: bezier.hpp:460
bezier_type::index_type index_type
Definition: bezier_curve_test_suite.hpp:37
point_type::Index index_type
Definition: bezier.hpp:116