Code-Eli  0.3.6
bezier_surface_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_surface_test_suite_hpp
14 #define bezier_surface_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"
28 
29 template<typename data__>
30 class bezier_surface_test_suite : public Test::Suite
31 {
32  private:
40 
41  tolerance_type tol;
42 
43  protected:
44  void AddTests(const float &)
45  {
46  // add the tests
65  }
66  void AddTests(const double &)
67  {
68  // add the tests
87  }
88  void AddTests(const long double &)
89  {
90  // add the tests
109  }
110 
111  public:
113  {
114  AddTests(data__());
115  }
117  {
118  }
119 
120  private:
121  void octave_print(int figno, /*const point_type pts[][], */const bezier_type &bez) const
122  {
123  index_type i, j;
124 
125  std::cout << "figure(" << figno << ");" << std::endl;
126  std::cout << "cp_x=[";
127  for (i=0; i<=bez.degree_u(); ++i)
128  {
129  std::cout << bez.get_control_point(i, 0).x();
130  for (j=1; j<bez.degree_v(); ++j)
131  {
132  std::cout << ", " << bez.get_control_point(i, j).x();
133  }
134  j=bez.degree_v();
135  std::cout << ", " << bez.get_control_point(i, j).x();
136  if (i<bez.degree_u())
137  std::cout << "; ";
138  }
139  std::cout << "];" << std::endl;
140 
141  std::cout << "cp_y=[";
142  for (i=0; i<=bez.degree_u(); ++i)
143  {
144  std::cout << bez.get_control_point(i, 0).y();
145  for (j=1; j<bez.degree_v(); ++j)
146  {
147  std::cout << ", " << bez.get_control_point(i, j).y();
148  }
149  j=bez.degree_v();
150  std::cout << ", " << bez.get_control_point(i, j).y();
151  if (i<bez.degree_u())
152  std::cout << "; ";
153  }
154  std::cout << "];" << std::endl;
155 
156  std::cout << "cp_z=[";
157  for (i=0; i<=bez.degree_u(); ++i)
158  {
159  std::cout << bez.get_control_point(i, 0).z();
160  for (j=1; j<bez.degree_v(); ++j)
161  {
162  std::cout << ", " << bez.get_control_point(i, j).z();
163  }
164  j=bez.degree_v();
165  std::cout << ", " << bez.get_control_point(i, j).z();
166  if (i<bez.degree_u())
167  std::cout << "; ";
168  }
169  std::cout << "];" << std::endl;
170 
171  // initialize the u & v parameters
172  std::vector<data__> u(51), v(51);
173  for (i=0; i<static_cast<index_type>(u.size()); ++i)
174  {
175  u[i]=static_cast<data__>(i)/(u.size()-1);
176  }
177  for (j=0; j<static_cast<index_type>(v.size()); ++j)
178  {
179  v[j]=static_cast<data__>(j)/(v.size()-1);
180  }
181 
182  // set the surface points
183  std::cout << "surf_x=[";
184  for (i=0; i<static_cast<index_type>(u.size()); ++i)
185  {
186  std::cout << bez.f(u[i], v[0]).x();
187  for (j=1; j<static_cast<index_type>(v.size()-1); ++j)
188  {
189  std::cout << ", " << bez.f(u[i], v[j]).x();
190  }
191  j=static_cast<index_type>(v.size()-1);
192  std::cout << ", " << bez.f(u[i], v[j]).x();
193  if (i<static_cast<index_type>(u.size()-1))
194  std::cout << "; " << std::endl;
195  }
196  std::cout << "];" << std::endl;
197 
198  std::cout << "surf_y=[";
199  for (i=0; i<static_cast<index_type>(u.size()); ++i)
200  {
201  std::cout << bez.f(u[i], v[0]).y();
202  for (j=1; j<static_cast<index_type>(v.size()-1); ++j)
203  {
204  std::cout << ", " << bez.f(u[i], v[j]).y();
205  }
206  j=static_cast<index_type>(v.size()-1);
207  std::cout << ", " << bez.f(u[i], v[j]).y();
208  if (i<static_cast<index_type>(u.size()-1))
209  std::cout << "; " << std::endl;
210  }
211  std::cout << "];" << std::endl;
212 
213  std::cout << "surf_z=[";
214  for (i=0; i<static_cast<index_type>(u.size()); ++i)
215  {
216  std::cout << bez.f(u[i], v[0]).z();
217  for (j=1; j<static_cast<index_type>(v.size()-1); ++j)
218  {
219  std::cout << ", " << bez.f(u[i], v[j]).z();
220  }
221  j=static_cast<index_type>(v.size()-1);
222  std::cout << ", " << bez.f(u[i], v[j]).z();
223  if (i<static_cast<index_type>(u.size()-1))
224  std::cout << "; " << std::endl;
225  }
226  std::cout << "];" << std::endl;
227 
228  std::cout << "setenv('GNUTERM', 'x11');" << std::endl;
229  std::cout << "mesh(surf_x, surf_y, surf_z, zeros(size(surf_z)), 'EdgeColor', [0 0 0]);" << std::endl;
230  std::cout << "hold on;" << std::endl;
231  std::cout << "plot3(cp_x, cp_y, cp_z, '-ok', 'MarkerFaceColor', [0 0 0]);" << std::endl;
232  std::cout << "plot3(cp_x', cp_y', cp_z', '-k');" << std::endl;
233  std::cout << "hold off;" << std::endl;
234  }
235 
236  void create_circle(std::vector<point_type> &/*pts*/)
237  {
238 #if 0
239  // NOTE: This will not create a closed circle, the last point will be
240  // the point just prior to the 2*pi point
241  size_t n=pts.size();
242  for (size_t i=0; i<n; ++i)
243  {
244  data__ theta(eli::constants::math<data__>::two_pi()*static_cast<data__>(i)/n);
245  pts[i](0)=std::cos(theta);
246  pts[i](1)=std::sin(theta);
247  pts[i](2)=0;
248  }
249 #endif
250  }
251 
253  {
254  index_type i, j, n(3), m(3);
255  point_type pt[3+1][3+1], pt_out;
256 
257  // create surface with specified control points
258  pt[0][0] << -15, 0, 15;
259  pt[1][0] << -5, 5, 15;
260  pt[2][0] << 5, 5, 15;
261  pt[3][0] << 15, 0, 15;
262  pt[0][1] << -15, 5, 5;
263  pt[1][1] << -5, 5, 5;
264  pt[2][1] << 5, 5, 5;
265  pt[3][1] << 15, 5, 5;
266  pt[0][2] << -15, 5, -5;
267  pt[1][2] << -5, 5, -5;
268  pt[2][2] << 5, 5, -5;
269  pt[3][2] << 15, 5, -5;
270  pt[0][3] << -15, 0, -15;
271  pt[1][3] << -5, 5, -15;
272  pt[2][3] << 5, 5, -15;
273  pt[3][3] << 15, 0, -15;
274 
275  // create surface with specified dimensions and set control points
276  bezier_type bez1(n, m);
277 
278  for (i=0; i<=n; ++i)
279  {
280  for (j=0; j<=m; ++j)
281  {
282  bez1.set_control_point(pt[i][j], i, j);
283  }
284  }
285  for (i=0; i<=n; ++i)
286  {
287  for (j=0; j<=m; ++j)
288  {
289  TEST_ASSERT(bez1.get_control_point(i, j) == pt[i][j]);
290  }
291  }
292 
293  // test order
294  TEST_ASSERT(bez1.degree_u()==n);
295  TEST_ASSERT(bez1.degree_v()==m);
296 
297  // test dimension
298  TEST_ASSERT(bez1.dimension()==3);
299 
300  // test copy ctr
301  bezier_type bez2(bez1);
302 
303  for (i=0; i<=n; ++i)
304  {
305  for (j=0; j<=m; ++j)
306  {
307  TEST_ASSERT(bez2.get_control_point(i, j) == bez1.get_control_point(i, j));
308  }
309  }
310 
311  // test equivalence operator
312  TEST_ASSERT(bez2==bez1);
313 
314  // test about equivalent
315  TEST_ASSERT(bez2.abouteq(bez1, 1e-6));
316 
317  // test assignment operator
318  bezier_type bez3;
319 
320  bez3=bez1;
321  TEST_ASSERT(bez3==bez1);
322 
323  // create surface then resize to needed dimensions
324  bezier_type bez4;
325 
326  bez4.resize(n, m);
327  for (i=0; i<=n; ++i)
328  {
329  for (j=0; j<=m; ++j)
330  {
331  bez4.set_control_point(pt[i][j], i, j);
332  }
333  }
334  for (i=0; i<=n; ++i)
335  {
336  for (j=0; j<=m; ++j)
337  {
338  TEST_ASSERT(bez4.get_control_point(i, j) == pt[i][j]);
339  }
340  }
341  }
342 
344  {
345  index_type i, j, n(3), m(3);
346  point_type pt[3+1][3+1], pt_out;
347 
348  // create surface with specified control points
349  pt[0][0] << -15, 0, 15;
350  pt[1][0] << -5, 5, 15;
351  pt[2][0] << 5, 5, 15;
352  pt[3][0] << 15, 0, 15;
353  pt[0][1] << -15, 5, 5;
354  pt[1][1] << -5, 5, 5;
355  pt[2][1] << 5, 5, 5;
356  pt[3][1] << 15, 5, 5;
357  pt[0][2] << -15, 5, -5;
358  pt[1][2] << -5, 5, -5;
359  pt[2][2] << 5, 5, -5;
360  pt[3][2] << 15, 5, -5;
361  pt[0][3] << -15, 0, -15;
362  pt[1][3] << -5, 5, -15;
363  pt[2][3] << 5, 5, -15;
364  pt[3][3] << 15, 0, -15;
365 
366  // create surface with specified dimensions and set control points
367  bezier_type bez1(n, m);
368 
369  for (i=0; i<=n; ++i)
370  {
371  for (j=0; j<=m; ++j)
372  {
373  bez1.set_control_point(pt[i][j], i, j);
374  }
375  }
376  for (i=0; i<=n; ++i)
377  {
378  for (j=0; j<=m; ++j)
379  {
380  TEST_ASSERT(bez1.get_control_point(i, j) == pt[i][j]);
381  }
382  }
383 
384  // test bounding box
385  typename bezier_type::bounding_box_type bb;
386  point_type pmin_ref, pmax_ref;
387 
388  bez1.get_bounding_box(bb);
389  pmin_ref << -15, 0, -15;
390  pmax_ref << 15, 5, 15;
391  TEST_ASSERT(bb.get_min()==pmin_ref);
392  TEST_ASSERT(bb.get_max()==pmax_ref);
393  }
394 
396  {
397  index_type i, j, n(3), m(3);
398  point_type pt[3+1][3+1], pt_out;
399 
400  // create surface with specified control points
401  pt[0][0] << -15, 0, 15;
402  pt[1][0] << -5, 5, 15;
403  pt[2][0] << 5, 5, 15;
404  pt[3][0] << 15, 0, 15;
405  pt[0][1] << -15, 5, 5;
406  pt[1][1] << -5, 5, 5;
407  pt[2][1] << 5, 5, 5;
408  pt[3][1] << 15, 5, 5;
409  pt[0][2] << -15, 5, -5;
410  pt[1][2] << -5, 5, -5;
411  pt[2][2] << 5, 5, -5;
412  pt[3][2] << 15, 5, -5;
413  pt[0][3] << -15, 0, -15;
414  pt[1][3] << -5, 5, -15;
415  pt[2][3] << 5, 5, -15;
416  pt[3][3] << 15, 0, -15;
417 
418  // create surface with specified dimensions and set control points
419  bezier_type bez1(n, m), bez2;
420 
421  for (i=0; i<=n; ++i)
422  {
423  for (j=0; j<=m; ++j)
424  {
425  bez1.set_control_point(pt[i][j], i, j);
426  }
427  }
428 
429  // reverse u-direction
430  bez2=bez1;
431  bez2.reverse_u();
432  for (i=0; i<=n; ++i)
433  {
434  for (j=0; j<=m; ++j)
435  {
436  TEST_ASSERT(bez2.get_control_point(i, j) == bez1.get_control_point(n-i, j));
437  }
438  }
439 
440  // reverse v-direction
441  bez2=bez1;
442  bez2.reverse_v();
443  for (i=0; i<=n; ++i)
444  {
445  for (j=0; j<=m; ++j)
446  {
447  TEST_ASSERT(bez2.get_control_point(i, j) == bez1.get_control_point(i, m-j));
448  }
449  }
450  }
451 
452  void swap_test()
453  {
454  index_type i, j, n(3), m(3);
455  point_type pt[3+1][3+1], pt_out;
456 
457  // create surface with specified control points
458  pt[0][0] << -15, 0, 15;
459  pt[1][0] << -5, 5, 15;
460  pt[2][0] << 5, 5, 15;
461  pt[3][0] << 15, 0, 15;
462  pt[0][1] << -15, 5, 5;
463  pt[1][1] << -5, 5, 5;
464  pt[2][1] << 5, 5, 5;
465  pt[3][1] << 15, 5, 5;
466  pt[0][2] << -15, 5, -5;
467  pt[1][2] << -5, 5, -5;
468  pt[2][2] << 5, 5, -5;
469  pt[3][2] << 15, 5, -5;
470  pt[0][3] << -15, 0, -15;
471  pt[1][3] << -5, 5, -15;
472  pt[2][3] << 5, 5, -15;
473  pt[3][3] << 15, 0, -15;
474 
475  // create surface with specified dimensions and set control points
476  bezier_type bez1(n, m), bez2;
477 
478  for (i=0; i<=n; ++i)
479  {
480  for (j=0; j<=m; ++j)
481  {
482  bez1.set_control_point(pt[i][j], i, j);
483  }
484  }
485 
486  // reverse swap u- & v-directions
487  bez2=bez1;
488  bez2.swap_uv();
489  TEST_ASSERT(bez1.degree_u()==bez2.degree_v());
490  TEST_ASSERT(bez1.degree_v()==bez2.degree_u());
491  for (i=0; i<=n; ++i)
492  {
493  for (j=0; j<=m; ++j)
494  {
495  TEST_ASSERT(bez1.get_control_point(i, j) == bez2.get_control_point(j, i));
496  }
497  }
498  }
499 
501  {
502  index_type n(3), m(3);
503  point_type pt[3+1][3+1], pt_out, pt_ref;
504  data_type u, v;
505 
506  // create surface with specified control points
507  pt[0][0] << -15, 0, 15;
508  pt[1][0] << -5, 5, 15;
509  pt[2][0] << 5, 5, 15;
510  pt[3][0] << 15, 0, 15;
511  pt[0][1] << -15, 5, 5;
512  pt[1][1] << -5, 5, 5;
513  pt[2][1] << 5, 5, 5;
514  pt[3][1] << 15, 5, 5;
515  pt[0][2] << -15, 5, -5;
516  pt[1][2] << -5, 5, -5;
517  pt[2][2] << 5, 5, -5;
518  pt[3][2] << 15, 5, -5;
519  pt[0][3] << -15, 0, -15;
520  pt[1][3] << -5, 5, -15;
521  pt[2][3] << 5, 5, -15;
522  pt[3][3] << 15, 0, -15;
523 
524  // create surface with specified dimensions and set control points
525  bezier_type bez(n, m), bez2;
526 
527  for (index_type i=0; i<=n; ++i)
528  {
529  for (index_type j=0; j<=m; ++j)
530  {
531  bez.set_control_point(pt[i][j], i, j);
532  }
533  }
534 
535  // test translation
536  {
537  point_type trans;
538 
539  // set up translation vector and apply
540  bez2=bez;
541  trans << 2, 1, 3;
542  bez2.translate(trans);
543 
544  // test evaluation at corners
545  u=0; v=0;
546  pt_out=bez2.f(u, v);
547  pt_ref=trans+pt[0][0];
548  TEST_ASSERT(pt_out==pt_ref);
549  u=1; v=0;
550  pt_out=bez2.f(u, v);
551  pt_ref=trans+pt[n][0];
552  TEST_ASSERT(pt_out==pt_ref);
553  u=0; v=1;
554  pt_out=bez2.f(u, v);
555  pt_ref=trans+pt[0][m];
556  TEST_ASSERT(pt_out==pt_ref);
557  u=1; v=1;
558  pt_out=bez2.f(u, v);
559  pt_ref=trans+pt[n][m];
560  TEST_ASSERT(pt_out==pt_ref);
561 
562  // test evaluation at interior point u=v=1/2
563  u=0.5; v=0.5;
564  pt_ref << 0, 4.6875 , 0;
565  pt_out=bez2.f(u, v);
566  pt_ref=trans+pt_ref;
567  TEST_ASSERT(pt_out==pt_ref);
568 
569  // test evaluation at interior point u=1/4, v=3/4
570  u=0.25; v=0.75;
571  pt_ref << -7.5, 4.04296875, -7.5;
572  pt_out=bez2.f(u, v);
573  pt_ref=trans+pt_ref;
574  TEST_ASSERT(pt_out==pt_ref);
575  }
576 
577  // test rotation about origin
578  {
579  typename bezier_type::rotation_matrix_type rmat;
580 
581  // set up rotation and apply
582  data_type one(1);
583  bez2=bez;
584  rmat << std::cos(one), 0, -std::sin(one),
585  0, one, 0,
586  std::sin(one), 0, std::cos(one);
587  bez2.rotate(rmat);
588 
589  // test evaluation at corners
590  u=0; v=0;
591  pt_out=bez2.f(u, v);
592  pt_ref=pt[0][0]*rmat.transpose();
593  TEST_ASSERT(pt_out==pt_ref);
594  u=1; v=0;
595  pt_out=bez2.f(u, v);
596  pt_ref=pt[n][0]*rmat.transpose();
597  TEST_ASSERT(pt_out==pt_ref);
598  u=0; v=1;
599  pt_out=bez2.f(u, v);
600  pt_ref=pt[0][m]*rmat.transpose();
601  TEST_ASSERT(pt_out==pt_ref);
602  u=1; v=1;
603  pt_out=bez2.f(u, v);
604  pt_ref=pt[n][m]*rmat.transpose();
605  TEST_ASSERT(pt_out==pt_ref);
606 
607  // test evaluation at interior point u=v=1/2
608  u=0.5; v=0.5;
609  pt_ref << 0, 4.6875 , 0;
610  pt_out=bez2.f(u, v);
611  pt_ref=pt_ref*rmat.transpose();
612  TEST_ASSERT(pt_out==pt_ref);
613 
614  // test evaluation at interior point u=1/4, v=3/4
615  u=0.25; v=0.75;
616  pt_ref << -7.5, 4.04296875, -7.5;
617  pt_out=bez2.f(u, v);
618  pt_ref=pt_ref*rmat.transpose();
619  TEST_ASSERT(tol.approximately_equal(pt_out,pt_ref));
620  }
621 
622  // test rotation about point
623  {
624  point_type rorig;
625  typename bezier_type::rotation_matrix_type rmat;
626 
627  // set up rotation and apply
628  data_type one(1);
629  bez2=bez;
630  rorig << 2, 1, 3;
631  rmat << std::cos(one), 0, -std::sin(one),
632  0, one, 0,
633  std::sin(one), 0, std::cos(one);
634  bez2.rotate(rmat, rorig);
635 
636  // test evaluation at corners
637  u=0; v=0;
638  pt_out=bez2.f(u, v);
639  pt_ref=rorig+(pt[0][0]-rorig)*rmat.transpose();
640  TEST_ASSERT(pt_out==pt_ref);
641  u=1; v=0;
642  pt_out=bez2.f(u, v);
643  pt_ref=rorig+(pt[n][0]-rorig)*rmat.transpose();
644  TEST_ASSERT(pt_out==pt_ref);
645  u=0; v=1;
646  pt_out=bez2.f(u, v);
647  pt_ref=rorig+(pt[0][m]-rorig)*rmat.transpose();
648  if (typeid(data_type)==typeid(float))
649  {
650  TEST_ASSERT(tol.approximately_equal(pt_out, pt_ref));
651  }
652  else
653  {
654  TEST_ASSERT(pt_out==pt_ref);
655  }
656  u=1; v=1;
657  pt_out=bez2.f(u, v);
658  pt_ref=rorig+(pt[n][m]-rorig)*rmat.transpose();
659  TEST_ASSERT(pt_out==pt_ref);
660 
661  // test evaluation at interior point u=v=1/2
662  u=0.5; v=0.5;
663  pt_ref << 0, 4.6875 , 0;
664  pt_out=bez2.f(u, v);
665  pt_ref=rorig+(pt_ref-rorig)*rmat.transpose();
666  TEST_ASSERT(tol.approximately_equal(pt_out,pt_ref));
667 
668  // test evaluation at interior point u=1/4, v=3/4
669  u=0.25; v=0.75;
670  pt_ref << -7.5, 4.04296875, -7.5;
671  pt_out=bez2.f(u, v);
672  pt_ref=rorig+(pt_ref-rorig)*rmat.transpose();
673  TEST_ASSERT(tol.approximately_equal(pt_out,pt_ref));
674  }
675  }
676 
678  {
679  index_type n(3), m(3);
680  point_type pt[3+1][3+1], pt_out, pt_ref;
681  data_type u, v;
682 
683  // create surface with specified control points
684  pt[0][0] << -15, 0, 15;
685  pt[1][0] << -5, 5, 15;
686  pt[2][0] << 5, 5, 15;
687  pt[3][0] << 15, 0, 15;
688  pt[0][1] << -15, 5, 5;
689  pt[1][1] << -5, 5, 5;
690  pt[2][1] << 5, 5, 5;
691  pt[3][1] << 15, 5, 5;
692  pt[0][2] << -15, 5, -5;
693  pt[1][2] << -5, 5, -5;
694  pt[2][2] << 5, 5, -5;
695  pt[3][2] << 15, 5, -5;
696  pt[0][3] << -15, 0, -15;
697  pt[1][3] << -5, 5, -15;
698  pt[2][3] << 5, 5, -15;
699  pt[3][3] << 15, 0, -15;
700 
701  // create surface with specified dimensions and set control points
702  bezier_type bez(n, m);
703 
704  for (index_type i=0; i<=n; ++i)
705  {
706  for (index_type j=0; j<=m; ++j)
707  {
708  bez.set_control_point(pt[i][j], i, j);
709  }
710  }
711 
712  // test evaluation at corners
713  u=0; v=0;
714  pt_out=bez.f(u, v);
715  TEST_ASSERT(pt_out==pt[0][0]);
716  u=1; v=0;
717  pt_out=bez.f(u, v);
718  TEST_ASSERT(pt_out==pt[n][0]);
719  u=0; v=1;
720  pt_out=bez.f(u, v);
721  TEST_ASSERT(pt_out==pt[0][m]);
722  u=1; v=1;
723  pt_out=bez.f(u, v);
724  TEST_ASSERT(pt_out==pt[n][m]);
725 
726  // test evaluation at interior point u=v=1/2
727  u=0.5; v=0.5;
728  pt_ref << 0, 4.6875 , 0;
729  pt_out=bez.f(u, v);
730  TEST_ASSERT(pt_out==pt_ref);
731 
732  // test evaluation at interior point u=1/4, v=3/4
733  u=0.25; v=0.75;
734  pt_ref << -7.5, 4.04296875, -7.5;
735  pt_out=bez.f(u, v);
736  TEST_ASSERT(pt_out==pt_ref);
737  }
738 
740  {
741  index_type n(3), m(3);
742  point_type pt[3+1][3+1], pt_out, pt_ref;
743  data_type u, v;
744 
745  // create surface with specified control points
746  pt[0][0] << -15, 0, 15;
747  pt[1][0] << -5, 5, 15;
748  pt[2][0] << 5, 5, 15;
749  pt[3][0] << 15, 0, 15;
750  pt[0][1] << -15, 5, 5;
751  pt[1][1] << -5, 5, 5;
752  pt[2][1] << 5, 5, 5;
753  pt[3][1] << 15, 5, 5;
754  pt[0][2] << -15, 5, -5;
755  pt[1][2] << -5, 5, -5;
756  pt[2][2] << 5, 5, -5;
757  pt[3][2] << 15, 5, -5;
758  pt[0][3] << -15, 0, -15;
759  pt[1][3] << -5, 5, -15;
760  pt[2][3] << 5, 5, -15;
761  pt[3][3] << 15, 0, -15;
762 
763  // create surface with specified dimensions and set control points
764  bezier_type bez(n, m);
765 
766  for (index_type i=0; i<=n; ++i)
767  {
768  for (index_type j=0; j<=m; ++j)
769  {
770  bez.set_control_point(pt[i][j], i, j);
771  }
772  }
773 
774  // test evaluation at corners
775  u=0; v=0;
776  pt_out=bez.f_u(u, v);
777  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(pt[1][0]-pt[0][0]));
778  pt_out=bez.f_v(u, v);
779  TEST_ASSERT(pt_out==static_cast<data_type>(m)*(pt[0][1]-pt[0][0]));
780  u=1; v=0;
781  pt_out=bez.f_u(u, v);
782  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(pt[3][0]-pt[2][0]));
783  pt_out=bez.f_v(u, v);
784  TEST_ASSERT(pt_out==static_cast<data_type>(m)*(pt[3][1]-pt[3][0]));
785  u=0; v=1;
786  pt_out=bez.f_u(u, v);
787  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(pt[1][3]-pt[0][3]));
788  pt_out=bez.f_v(u, v);
789  TEST_ASSERT(pt_out==static_cast<data_type>(m)*(pt[0][3]-pt[0][2]));
790  u=1; v=1;
791  pt_out=bez.f_u(u, v);
792  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(pt[3][3]-pt[2][3]));
793  pt_out=bez.f_v(u, v);
794  TEST_ASSERT(pt_out==static_cast<data_type>(m)*(pt[3][3]-pt[3][2]));
795 
796  // test evaluation at interior point u=v=1/2
797  u=0.5; v=0.5;
798  pt_ref << 30, 0, 0;
799  pt_out=bez.f_u(u, v);
800  TEST_ASSERT(pt_out==pt_ref);
801  pt_ref << 0, 0, -30;
802  pt_out=bez.f_v(u, v);
803  TEST_ASSERT(pt_out==pt_ref);
804 
805  // test evaluation at interior point u=1/4, v=3/4
806  u=0.25; v=0.75;
807  pt_ref << 30, 3.28125, 0;
808  pt_out=bez.f_u(u, v);
809  TEST_ASSERT(pt_out==pt_ref);
810  pt_ref << 0, -3.28125, -30;
811  pt_out=bez.f_v(u, v);
812  TEST_ASSERT(pt_out==pt_ref);
813  }
814 
816  {
817  index_type n(3), m(3);
818  point_type pt[3+1][3+1], pt_out, pt_ref;
819  data_type u, v;
820 
821  // create surface with specified control points
822  pt[0][0] << -15, 0, 15;
823  pt[1][0] << -5, 5, 15;
824  pt[2][0] << 5, 5, 15;
825  pt[3][0] << 15, 0, 15;
826  pt[0][1] << -15, 5, 5;
827  pt[1][1] << -5, 5, 5;
828  pt[2][1] << 5, 5, 5;
829  pt[3][1] << 15, 5, 5;
830  pt[0][2] << -15, 5, -5;
831  pt[1][2] << -5, 5, -5;
832  pt[2][2] << 5, 5, -5;
833  pt[3][2] << 15, 5, -5;
834  pt[0][3] << -15, 0, -15;
835  pt[1][3] << -5, 5, -15;
836  pt[2][3] << 5, 5, -15;
837  pt[3][3] << 15, 0, -15;
838 
839  // create surface with specified dimensions and set control points
840  bezier_type bez(n, m);
841 
842  for (index_type i=0; i<=n; ++i)
843  {
844  for (index_type j=0; j<=m; ++j)
845  {
846  bez.set_control_point(pt[i][j], i, j);
847  }
848  }
849 
850  // test evaluation at corners
851  u=0; v=0;
852  pt_out=bez.f_uu(u, v);
853  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(n-1)*(pt[2][0]-2*pt[1][0]+pt[0][0]));
854  pt_out=bez.f_uv(u, v);
855  TEST_ASSERT(pt_out==static_cast<data_type>(n)*m*(pt[1][1]-pt[1][0]-pt[0][1]+pt[0][0]));
856  pt_out=bez.f_vv(u, v);
857  TEST_ASSERT(pt_out==static_cast<data_type>(m)*(m-1)*(pt[0][2]-2*pt[0][1]+pt[0][0]));
858  u=1; v=0;
859  pt_out=bez.f_uu(u, v);
860  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(n-1)*(pt[3][0]-2*pt[2][0]+pt[1][0]));
861  pt_out=bez.f_uv(u, v);
862  TEST_ASSERT(pt_out==static_cast<data_type>(n)*m*(pt[3][1]-pt[3][0]-pt[2][1]+pt[2][0]));
863  pt_out=bez.f_vv(u, v);
864  TEST_ASSERT(pt_out==static_cast<data_type>(m)*(m-1)*(pt[3][2]-2*pt[3][1]+pt[3][0]));
865  u=0; v=1;
866  pt_out=bez.f_uu(u, v);
867  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(n-1)*(pt[2][3]-2*pt[1][3]+pt[0][3]));
868  pt_out=bez.f_uv(u, v);
869  TEST_ASSERT(pt_out==static_cast<data_type>(n)*m*(pt[1][3]-pt[1][2]-pt[0][3]+pt[0][2]));
870  pt_out=bez.f_vv(u, v);
871  TEST_ASSERT(pt_out==static_cast<data_type>(m)*(m-1)*(pt[0][3]-2*pt[0][2]+pt[0][1]));
872  u=1; v=1;
873  pt_out=bez.f_uu(u, v);
874  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(n-1)*(pt[3][3]-2*pt[2][3]+pt[1][3]));
875  pt_out=bez.f_uv(u, v);
876  TEST_ASSERT(pt_out==static_cast<data_type>(n)*m*(pt[3][3]-pt[3][2]-pt[2][3]+pt[2][2]));
877  pt_out=bez.f_vv(u, v);
878  TEST_ASSERT(pt_out==static_cast<data_type>(m)*(m-1)*(pt[3][3]-2*pt[3][2]+pt[3][1]));
879 
880  // test evaluation at interior point u=v=1/2
881  u=0.5; v=0.5;
882  pt_ref << 0, -7.5, 0;
883  pt_out=bez.f_uu(u, v);
884  TEST_ASSERT(pt_out==pt_ref);
885  pt_ref << 0, 0, 0;
886  pt_out=bez.f_uv(u, v);
887  TEST_ASSERT(pt_out==pt_ref);
888  pt_ref << 0, -7.5, 0;
889  pt_out=bez.f_vv(u, v);
890  TEST_ASSERT(pt_out==pt_ref);
891 
892  // test evaluation at interior point u=1/4, v=3/4
893  u=0.25; v=0.75;
894  pt_ref << 0, -13.125, 0;
895  pt_out=bez.f_uu(u, v);
896  TEST_ASSERT(pt_out==pt_ref);
897  pt_ref << 0, 11.25, 0;
898  pt_out=bez.f_uv(u, v);
899  TEST_ASSERT(pt_out==pt_ref);
900  pt_ref << 0, -13.125, 0;
901  pt_out=bez.f_vv(u, v);
902  TEST_ASSERT(pt_out==pt_ref);
903  }
904 
906  {
907  index_type n(3), m(3);
908  point_type pt[3+1][3+1], pt_out, pt_ref;
909  data_type u, v;
910 
911  // create surface with specified control points
912  pt[0][0] << -15, 0, 15;
913  pt[1][0] << -5, 5, 15;
914  pt[2][0] << 5, 5, 15;
915  pt[3][0] << 15, 0, 15;
916  pt[0][1] << -15, 5, 5;
917  pt[1][1] << -5, 5, 5;
918  pt[2][1] << 5, 5, 5;
919  pt[3][1] << 15, 5, 5;
920  pt[0][2] << -15, 5, -5;
921  pt[1][2] << -5, 5, -5;
922  pt[2][2] << 5, 5, -5;
923  pt[3][2] << 15, 5, -5;
924  pt[0][3] << -15, 0, -15;
925  pt[1][3] << -5, 5, -15;
926  pt[2][3] << 5, 5, -15;
927  pt[3][3] << 15, 0, -15;
928 
929  // create surface with specified dimensions and set control points
930  bezier_type bez(n, m);
931 
932  for (index_type i=0; i<=n; ++i)
933  {
934  for (index_type j=0; j<=m; ++j)
935  {
936  bez.set_control_point(pt[i][j], i, j);
937  }
938  }
939 
940  // test evaluation at corners
941  u=0; v=0;
942  pt_out=bez.f_uuu(u, v);
943  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(n-1)*(n-2)*(pt[3][0]-3*pt[2][0]+3*pt[1][0]-pt[0][0]));
944  pt_out=bez.f_uuv(u, v);
945  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(n-1)*m*((pt[2][1]-pt[2][0])-2*(pt[1][1]-pt[1][0])+(pt[0][1]-pt[0][0])));
946  pt_out=bez.f_uvv(u, v);
947  TEST_ASSERT(pt_out==static_cast<data_type>(n)*m*(m-1)*((pt[1][2]-pt[0][2])-2*(pt[1][1]-pt[0][1])+(pt[1][0]-pt[0][0])));
948  pt_out=bez.f_vvv(u, v);
949  TEST_ASSERT(pt_out==static_cast<data_type>(m)*(m-1)*(m-2)*(pt[0][3]-3*pt[0][2]+3*pt[0][1]-pt[0][0]));
950  u=1; v=0;
951  pt_out=bez.f_uuu(u, v);
952  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(n-1)*(n-2)*(pt[3][0]-3*pt[2][0]+3*pt[1][0]-pt[0][0]));
953  pt_out=bez.f_uuv(u, v);
954  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(n-1)*m*((pt[3][1]-pt[3][0])-2*(pt[2][1]-pt[2][0])+(pt[1][1]-pt[1][0])));
955  pt_out=bez.f_uvv(u, v);
956  TEST_ASSERT(pt_out==static_cast<data_type>(n)*m*(m-1)*((pt[3][2]-pt[2][2])-2*(pt[3][1]-pt[2][1])+(pt[3][0]-pt[2][0])));
957  pt_out=bez.f_vvv(u, v);
958  TEST_ASSERT(pt_out==static_cast<data_type>(m)*(m-1)*(m-2)*(pt[3][3]-3*pt[3][2]+3*pt[3][1]-pt[3][0]));
959  u=0; v=1;
960  pt_out=bez.f_uuu(u, v);
961  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(n-1)*(n-2)*(pt[3][3]-3*pt[2][3]+3*pt[1][3]-pt[0][3]));
962  pt_out=bez.f_uuv(u, v);
963  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(n-1)*m*((pt[2][3]-pt[2][2])-2*(pt[1][3]-pt[1][2])+(pt[0][3]-pt[0][2])));
964  pt_out=bez.f_uvv(u, v);
965  TEST_ASSERT(pt_out==static_cast<data_type>(n)*m*(m-1)*((pt[1][3]-pt[0][3])-2*(pt[1][2]-pt[0][2])+(pt[1][1]-pt[0][1])));
966  pt_out=bez.f_vvv(u, v);
967  TEST_ASSERT(pt_out==static_cast<data_type>(m)*(m-1)*(m-2)*(pt[0][3]-3*pt[0][2]+3*pt[0][1]-pt[0][0]));
968  u=1; v=1;
969  pt_out=bez.f_uuu(u, v);
970  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(n-1)*(n-2)*(pt[3][3]-3*pt[2][3]+3*pt[1][3]-pt[0][3]));
971  pt_out=bez.f_uuv(u, v);
972  TEST_ASSERT(pt_out==static_cast<data_type>(n)*(n-1)*m*((pt[3][3]-pt[3][2])-2*(pt[2][3]-pt[2][2])+(pt[1][3]-pt[1][2])));
973  pt_out=bez.f_uvv(u, v);
974  TEST_ASSERT(pt_out==static_cast<data_type>(n)*m*(m-1)*((pt[3][3]-pt[2][3])-2*(pt[3][2]-pt[2][2])+(pt[3][1]-pt[2][1])));
975  pt_out=bez.f_vvv(u, v);
976  TEST_ASSERT(pt_out==static_cast<data_type>(m)*(m-1)*(m-2)*(pt[3][3]-3*pt[3][2]+3*pt[3][1]-pt[3][0]));
977 
978  // test evaluation at interior point u=v=1/2
979  u=0.5; v=0.5;
980  pt_ref << 0, 0, 0;
981  pt_out=bez.f_uuu(u, v);
982  TEST_ASSERT(pt_out==pt_ref);
983  pt_ref << 0, 0, 0;
984  pt_out=bez.f_uuv(u, v);
985  TEST_ASSERT(pt_out==pt_ref);
986  pt_ref << 0, 0, 0;
987  pt_out=bez.f_uvv(u, v);
988  TEST_ASSERT(pt_out==pt_ref);
989  pt_ref << 0, 0, 0;
990  pt_out=bez.f_vvv(u, v);
991  TEST_ASSERT(pt_out==pt_ref);
992 
993  // test evaluation at interior point u=1/4, v=3/4
994  u=0.25; v=0.75;
995  pt_ref << 0, 0, 0;
996  pt_out=bez.f_uuu(u, v);
997  TEST_ASSERT(pt_out==pt_ref);
998  pt_ref << 0, -45, 0;
999  pt_out=bez.f_uuv(u, v);
1000  TEST_ASSERT(pt_out==pt_ref);
1001  pt_ref << 0, 45, 0;
1002  pt_out=bez.f_uvv(u, v);
1003  TEST_ASSERT(pt_out==pt_ref);
1004  pt_ref << 0, 0, 0;
1005  pt_out=bez.f_vvv(u, v);
1006  TEST_ASSERT(pt_out==pt_ref);
1007  }
1008 
1010  {
1011  index_type n(3), m(3);
1012  point_type pt[3+1][3+1], pt_out, pt_ref;
1013  data_type u, v;
1014 
1015  // create surface with specified control points
1016  pt[0][0] << -15, 0, 15;
1017  pt[1][0] << -5, 5, 15;
1018  pt[2][0] << 5, 5, 15;
1019  pt[3][0] << 15, 0, 15;
1020  pt[0][1] << -15, 5, 5;
1021  pt[1][1] << -5, 5, 5;
1022  pt[2][1] << 5, 5, 5;
1023  pt[3][1] << 15, 5, 5;
1024  pt[0][2] << -15, 5, -5;
1025  pt[1][2] << -5, 5, -5;
1026  pt[2][2] << 5, 5, -5;
1027  pt[3][2] << 15, 5, -5;
1028  pt[0][3] << -15, 0, -15;
1029  pt[1][3] << -5, 5, -15;
1030  pt[2][3] << 5, 5, -15;
1031  pt[3][3] << 15, 0, -15;
1032 
1033  // create surface with specified dimensions and set control points
1034  bezier_type bez(n, m);
1035 
1036  for (index_type i=0; i<=n; ++i)
1037  {
1038  for (index_type j=0; j<=m; ++j)
1039  {
1040  bez.set_control_point(pt[i][j], i, j);
1041  }
1042  }
1043 
1044  // test mean curvature
1045  {
1046  data_type curv_out, curv_ref;
1047 
1048  // test evaluation at corners
1049  u=0; v=0;
1050  curv_ref=static_cast<data_type>(-0.015876322406928);
1051  eli::geom::surface::mean_curvature(curv_out, bez, u, v);
1052  TEST_ASSERT(std::abs(curv_ref-curv_out)<1e-6);
1053  u=1; v=0;
1054  curv_ref=static_cast<data_type>(-0.015876322406928);
1055  eli::geom::surface::mean_curvature(curv_out, bez, u, v);
1056  TEST_ASSERT(std::abs(curv_ref-curv_out)<1e-6);
1057  u=0; v=1;
1058  curv_ref=static_cast<data_type>(-0.015876322406928);
1059  eli::geom::surface::mean_curvature(curv_out, bez, u, v);
1060  TEST_ASSERT(std::abs(curv_ref-curv_out)<1e-6);
1061  u=1; v=1;
1062  curv_ref=static_cast<data_type>(-0.015876322406928);
1063  eli::geom::surface::mean_curvature(curv_out, bez, u, v);
1064  TEST_ASSERT(std::abs(curv_ref-curv_out)<1e-6);
1065 
1066  // test at interior point u=v=1/2
1067  u=0.5; v=0.5;
1068  curv_ref=static_cast<data_type>(-0.008333333333333);
1069  eli::geom::surface::mean_curvature(curv_out, bez, u, v);
1070  TEST_ASSERT(std::abs(curv_ref-curv_out)<1e-6);
1071 
1072  // test at interior point u=1/4 & v=3/4
1073  u=0.25; v=0.75;
1074  curv_ref=static_cast<data_type>(-0.014099238414151);
1075  eli::geom::surface::mean_curvature(curv_out, bez, u, v);
1076  TEST_ASSERT(std::abs(curv_ref-curv_out)<1e-6);
1077  }
1078 
1079  // test Gaussian curvature
1080  {
1081  data_type curv_out, curv_ref;
1082 
1083  // test evaluation at corners
1084  u=0; v=0;
1085  curv_ref=static_cast<data_type>(-0.00061728395);
1086  eli::geom::surface::gaussian_curvature(curv_out, bez, u, v);
1087  TEST_ASSERT(std::abs(curv_ref-curv_out)<1e-10);
1088  u=1; v=0;
1089  curv_ref=static_cast<data_type>(-0.00061728395);
1090  eli::geom::surface::gaussian_curvature(curv_out, bez, u, v);
1091  TEST_ASSERT(std::abs(curv_ref-curv_out)<1e-10);
1092  u=0; v=1;
1093  curv_ref=static_cast<data_type>(-0.00061728395);
1094  eli::geom::surface::gaussian_curvature(curv_out, bez, u, v);
1095  TEST_ASSERT(std::abs(curv_ref-curv_out)<1e-10);
1096  u=1; v=1;
1097  curv_ref=static_cast<data_type>(-0.00061728395);
1098  eli::geom::surface::gaussian_curvature(curv_out, bez, u, v);
1099  TEST_ASSERT(std::abs(curv_ref-curv_out)<1e-10);
1100 
1101  // test at interior point u=v=1/2
1102  u=0.5; v=0.5;
1103  curv_ref=static_cast<data_type>(6.944444444444444e-5);
1104  eli::geom::surface::gaussian_curvature(curv_out, bez, u, v);
1105  TEST_ASSERT(std::abs(curv_ref-curv_out)<1e-10);
1106 
1107  // test at interior point u=1/4 & v=3/4
1108  u=0.25; v=0.75;
1109  curv_ref=static_cast<data_type>(5.381754978392456e-05);
1110  eli::geom::surface::gaussian_curvature(curv_out, bez, u, v);
1111  TEST_ASSERT(std::abs(curv_ref-curv_out)<1e-10);
1112  }
1113 
1114  // test principal curvatures without directions
1115  {
1116  data_type kmax_out, kmin_out, H_out, K_out;
1117 
1118  // test evaluation at corners
1119  u=0; v=0;
1120  eli::geom::surface::principal_curvature(kmax_out, kmin_out, bez, u, v);
1121  eli::geom::surface::mean_curvature(H_out, bez, u, v);
1122  eli::geom::surface::gaussian_curvature(K_out, bez, u, v);
1123  TEST_ASSERT(std::abs((kmax_out+kmin_out)/2-H_out)<std::numeric_limits<data_type>::epsilon());
1124  TEST_ASSERT(std::abs((kmax_out*kmin_out)-K_out)<std::numeric_limits<data_type>::epsilon());
1125  u=1; v=0;
1126  eli::geom::surface::principal_curvature(kmax_out, kmin_out, bez, u, v);
1127  eli::geom::surface::mean_curvature(H_out, bez, u, v);
1128  eli::geom::surface::gaussian_curvature(K_out, bez, u, v);
1129  TEST_ASSERT(std::abs((kmax_out+kmin_out)/2-H_out)<std::numeric_limits<data_type>::epsilon());
1130  TEST_ASSERT(std::abs((kmax_out*kmin_out)-K_out)<std::numeric_limits<data_type>::epsilon());
1131  u=0; v=1;
1132  eli::geom::surface::principal_curvature(kmax_out, kmin_out, bez, u, v);
1133  eli::geom::surface::mean_curvature(H_out, bez, u, v);
1134  eli::geom::surface::gaussian_curvature(K_out, bez, u, v);
1135  TEST_ASSERT(std::abs((kmax_out+kmin_out)/2-H_out)<std::numeric_limits<data_type>::epsilon());
1136  TEST_ASSERT(std::abs((kmax_out*kmin_out)-K_out)<std::numeric_limits<data_type>::epsilon());
1137  u=1; v=1;
1138  eli::geom::surface::principal_curvature(kmax_out, kmin_out, bez, u, v);
1139  eli::geom::surface::mean_curvature(H_out, bez, u, v);
1140  eli::geom::surface::gaussian_curvature(K_out, bez, u, v);
1141  TEST_ASSERT(std::abs((kmax_out+kmin_out)/2-H_out)<std::numeric_limits<data_type>::epsilon());
1142  TEST_ASSERT(std::abs((kmax_out*kmin_out)-K_out)<std::numeric_limits<data_type>::epsilon());
1143 
1144  // test at interior point u=v=1/2
1145  u=0.5; v=0.5;
1146  eli::geom::surface::principal_curvature(kmax_out, kmin_out, bez, u, v);
1147  eli::geom::surface::mean_curvature(H_out, bez, u, v);
1148  eli::geom::surface::gaussian_curvature(K_out, bez, u, v);
1149  TEST_ASSERT(std::abs((kmax_out+kmin_out)/2-H_out)<std::numeric_limits<data_type>::epsilon());
1150  TEST_ASSERT(std::abs((kmax_out*kmin_out)-K_out)<std::numeric_limits<data_type>::epsilon());
1151 
1152  // test at interior point u=1/4 & v=3/4
1153  u=0.25; v=0.75;
1154  eli::geom::surface::principal_curvature(kmax_out, kmin_out, bez, u, v);
1155  eli::geom::surface::mean_curvature(H_out, bez, u, v);
1156  eli::geom::surface::gaussian_curvature(K_out, bez, u, v);
1157  TEST_ASSERT(std::abs((kmax_out+kmin_out)/2-H_out)<std::numeric_limits<data_type>::epsilon());
1158  TEST_ASSERT(std::abs((kmax_out*kmin_out)-K_out)<std::numeric_limits<data_type>::epsilon());
1159  }
1160 
1161  // test principal curvatures with directions
1162  {
1163  data_type kmax_out, kmin_out, H_out, K_out;
1164  point_type kmax_dir_out, kmin_dir_out, n_out, kmax_dir_ref, kmin_dir_ref, n_ref;
1165 
1166  // test evaluation at corners
1167  u=0; v=0;
1168  eli::geom::surface::principal_curvature(kmax_out, kmin_out, kmax_dir_out, kmin_dir_out, n_out, bez, u, v);
1169  eli::geom::surface::mean_curvature(H_out, bez, u, v);
1170  eli::geom::surface::gaussian_curvature(K_out, bez, u, v);
1171  kmax_dir_ref << static_cast<data_type>( 0.70710678), static_cast<data_type>(0), static_cast<data_type>( 0.70710678);
1172  kmin_dir_ref << static_cast<data_type>( 0.57735027), static_cast<data_type>(0.57735027), static_cast<data_type>(-0.57735027);
1173  n_ref << static_cast<data_type>(-0.40824829), static_cast<data_type>(0.81649658), static_cast<data_type>( 0.40824829);
1174  TEST_ASSERT(std::abs((kmax_out+kmin_out)/2-H_out)<std::numeric_limits<data_type>::epsilon());
1175  TEST_ASSERT(std::abs((kmax_out*kmin_out)-K_out)<std::numeric_limits<data_type>::epsilon());
1176  TEST_ASSERT((kmax_dir_out - kmax_dir_ref).norm()<1e-7);
1177  TEST_ASSERT((kmin_dir_out - kmin_dir_ref).norm()<1e-7);
1178  TEST_ASSERT((n_out - n_ref).norm()<1e-7);
1179  TEST_ASSERT(kmax_dir_out.dot(kmin_dir_out)<std::numeric_limits<data_type>::epsilon());
1180  TEST_ASSERT(kmax_dir_out.dot(n_out)<std::numeric_limits<data_type>::epsilon());
1181  TEST_ASSERT(kmin_dir_out.dot(n_out)<std::numeric_limits<data_type>::epsilon());
1182  u=1; v=0;
1183  eli::geom::surface::principal_curvature(kmax_out, kmin_out, kmax_dir_out, kmin_dir_out, n_out, bez, u, v);
1184  eli::geom::surface::mean_curvature(H_out, bez, u, v);
1185  eli::geom::surface::gaussian_curvature(K_out, bez, u, v);
1186  kmax_dir_ref << static_cast<data_type>( 0.70710678), static_cast<data_type>(0), static_cast<data_type>(-0.70710678);
1187  kmin_dir_ref << static_cast<data_type>(-0.57735027), static_cast<data_type>(0.57735027), static_cast<data_type>(-0.57735027);
1188  n_ref << static_cast<data_type>( 0.40824829), static_cast<data_type>(0.81649658), static_cast<data_type>( 0.40824829);
1189  TEST_ASSERT(std::abs((kmax_out+kmin_out)/2-H_out)<std::numeric_limits<data_type>::epsilon());
1190  TEST_ASSERT(std::abs((kmax_out*kmin_out)-K_out)<std::numeric_limits<data_type>::epsilon());
1191  TEST_ASSERT((kmax_dir_out - kmax_dir_ref).norm()<1e-7);
1192  TEST_ASSERT((kmin_dir_out - kmin_dir_ref).norm()<1e-7);
1193  TEST_ASSERT((n_out - n_ref).norm()<1e-7);
1194  TEST_ASSERT(kmax_dir_out.dot(kmin_dir_out)<std::numeric_limits<data_type>::epsilon());
1195  TEST_ASSERT(kmax_dir_out.dot(n_out)<std::numeric_limits<data_type>::epsilon());
1196  TEST_ASSERT(kmin_dir_out.dot(n_out)<std::numeric_limits<data_type>::epsilon());
1197  u=0; v=1;
1198  eli::geom::surface::principal_curvature(kmax_out, kmin_out, kmax_dir_out, kmin_dir_out, n_out, bez, u, v);
1199  eli::geom::surface::mean_curvature(H_out, bez, u, v);
1200  eli::geom::surface::gaussian_curvature(K_out, bez, u, v);
1201  kmax_dir_ref << static_cast<data_type>( 0.70710678), static_cast<data_type>( 0), static_cast<data_type>(-0.70710678);
1202  kmin_dir_ref << static_cast<data_type>(-0.57735027), static_cast<data_type>(-0.57735027), static_cast<data_type>(-0.57735027);
1203  n_ref << static_cast<data_type>(-0.40824829), static_cast<data_type>( 0.81649658), static_cast<data_type>(-0.40824829);
1204  TEST_ASSERT(std::abs((kmax_out+kmin_out)/2-H_out)<std::numeric_limits<data_type>::epsilon());
1205  TEST_ASSERT(std::abs((kmax_out*kmin_out)-K_out)<std::numeric_limits<data_type>::epsilon());
1206  TEST_ASSERT((kmax_dir_out - kmax_dir_ref).norm()<1e-7);
1207  TEST_ASSERT((kmin_dir_out - kmin_dir_ref).norm()<1e-7);
1208  TEST_ASSERT((n_out - n_ref).norm()<1e-7);
1209  TEST_ASSERT(kmax_dir_out.dot(kmin_dir_out)<std::numeric_limits<data_type>::epsilon());
1210  TEST_ASSERT(kmax_dir_out.dot(n_out)<std::numeric_limits<data_type>::epsilon());
1211  TEST_ASSERT(kmin_dir_out.dot(n_out)<std::numeric_limits<data_type>::epsilon());
1212  u=1; v=1;
1213  eli::geom::surface::principal_curvature(kmax_out, kmin_out, kmax_dir_out, kmin_dir_out, n_out, bez, u, v);
1214  eli::geom::surface::mean_curvature(H_out, bez, u, v);
1215  eli::geom::surface::gaussian_curvature(K_out, bez, u, v);
1216  kmax_dir_ref << static_cast<data_type>( 0.70710678), static_cast<data_type>( 0), static_cast<data_type>( 0.70710678);
1217  kmin_dir_ref << static_cast<data_type>( 0.57735027), static_cast<data_type>(-0.57735027), static_cast<data_type>(-0.57735027);
1218  n_ref << static_cast<data_type>( 0.40824829), static_cast<data_type>( 0.81649658), static_cast<data_type>(-0.40824829);
1219  TEST_ASSERT(std::abs((kmax_out+kmin_out)/2-H_out)<std::numeric_limits<data_type>::epsilon());
1220  TEST_ASSERT(std::abs((kmax_out*kmin_out)-K_out)<std::numeric_limits<data_type>::epsilon());
1221  TEST_ASSERT((kmax_dir_out - kmax_dir_ref).norm()<1e-7);
1222  TEST_ASSERT((kmin_dir_out - kmin_dir_ref).norm()<1e-7);
1223  TEST_ASSERT((n_out - n_ref).norm()<1e-7);
1224  TEST_ASSERT(kmax_dir_out.dot(kmin_dir_out)<std::numeric_limits<data_type>::epsilon());
1225  TEST_ASSERT(kmax_dir_out.dot(n_out)<std::numeric_limits<data_type>::epsilon());
1226  TEST_ASSERT(kmin_dir_out.dot(n_out)<std::numeric_limits<data_type>::epsilon());
1227 
1228  // test at interior point u=v=1/2
1229  u=0.5; v=0.5;
1230  eli::geom::surface::principal_curvature(kmax_out, kmin_out, kmax_dir_out, kmin_dir_out, n_out, bez, u, v);
1231  eli::geom::surface::mean_curvature(H_out, bez, u, v);
1232  eli::geom::surface::gaussian_curvature(K_out, bez, u, v);
1233  kmax_dir_ref << 1, 0, 0;
1234  kmin_dir_ref << 0, 0, -1;
1235  n_ref << 0, 1, 0;
1236  TEST_ASSERT(std::abs((kmax_out+kmin_out)/2-H_out)<std::numeric_limits<data_type>::epsilon());
1237  TEST_ASSERT(std::abs((kmax_out*kmin_out)-K_out)<std::numeric_limits<data_type>::epsilon());
1238  TEST_ASSERT((kmax_dir_out - kmax_dir_ref).norm()<1e-7);
1239  TEST_ASSERT((kmin_dir_out - kmin_dir_ref).norm()<1e-7);
1240  TEST_ASSERT((n_out - n_ref).norm()<1e-7);
1241  TEST_ASSERT(kmax_dir_out.dot(kmin_dir_out)<std::numeric_limits<data_type>::epsilon());
1242  TEST_ASSERT(kmax_dir_out.dot(n_out)<std::numeric_limits<data_type>::epsilon());
1243  TEST_ASSERT(kmin_dir_out.dot(n_out)<std::numeric_limits<data_type>::epsilon());
1244 
1245  // test at interior point u=1/4 & v=3/4
1246  u=0.25; v=0.75;
1247  eli::geom::surface::principal_curvature(kmax_out, kmin_out, kmax_dir_out, kmin_dir_out, n_out, bez, u, v);
1248  eli::geom::surface::mean_curvature(H_out, bez, u, v);
1249  eli::geom::surface::gaussian_curvature(K_out, bez, u, v);
1250  kmax_dir_ref << static_cast<data_type>( 0.70710678), static_cast<data_type>(0), static_cast<data_type>(-0.70710678);
1251  kmin_dir_ref << static_cast<data_type>(-0.69879657), static_cast<data_type>(-0.15286175), static_cast<data_type>(-0.69879657);
1252  n_ref << static_cast<data_type>(-0.10808958), static_cast<data_type>( 0.98824758), static_cast<data_type>(-0.10808958);
1253  TEST_ASSERT(std::abs((kmax_out+kmin_out)/2-H_out)<std::numeric_limits<data_type>::epsilon());
1254  TEST_ASSERT(std::abs((kmax_out*kmin_out)-K_out)<std::numeric_limits<data_type>::epsilon());
1255  TEST_ASSERT((kmax_dir_out - kmax_dir_ref).norm()<1e-7);
1256  TEST_ASSERT((kmin_dir_out - kmin_dir_ref).norm()<1e-7);
1257  TEST_ASSERT((n_out - n_ref).norm()<1e-7);
1258  TEST_ASSERT(kmax_dir_out.dot(kmin_dir_out)<std::numeric_limits<data_type>::epsilon());
1259  TEST_ASSERT(kmax_dir_out.dot(n_out)<std::numeric_limits<data_type>::epsilon());
1260  TEST_ASSERT(kmin_dir_out.dot(n_out)<std::numeric_limits<data_type>::epsilon());
1261  }
1262  }
1263 
1265  {
1266  index_type n(3), m(3);
1267  point_type pt[3+1][3+1];
1268  data_type u, v;
1269 
1270  // create surface with specified control points
1271  pt[0][0] << -15, 0, 15;
1272  pt[1][0] << -5, 5, 15;
1273  pt[2][0] << 5, 5, 15;
1274  pt[3][0] << 15, 0, 15;
1275  pt[0][1] << -15, 5, 5;
1276  pt[1][1] << -5, 5, 5;
1277  pt[2][1] << 5, 5, 5;
1278  pt[3][1] << 15, 5, 5;
1279  pt[0][2] << -15, 5, -5;
1280  pt[1][2] << -5, 5, -5;
1281  pt[2][2] << 5, 5, -5;
1282  pt[3][2] << 15, 5, -5;
1283  pt[0][3] << -15, 0, -15;
1284  pt[1][3] << -5, 5, -15;
1285  pt[2][3] << 5, 5, -15;
1286  pt[3][3] << 15, 0, -15;
1287 
1288  // create surface with specified dimensions and set control points
1289  bezier_type bez(n, m);
1290 
1291  for (index_type i=0; i<=n; ++i)
1292  {
1293  for (index_type j=0; j<=m; ++j)
1294  {
1295  bez.set_control_point(pt[i][j], i, j);
1296  }
1297  }
1298 
1299  // test promote in u-direction
1300  {
1301  bezier_type bez2(bez);
1302 
1303  // promote
1304  bez2.promote_u();
1305 
1306  // test if degree increased
1307  TEST_ASSERT(bez2.degree_u()==bez.degree_u()+1);
1308 
1309  // test evaluation at u=v=1/2
1310  u=0.5; v=0.5;
1311  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
1312  TEST_ASSERT(bez.f_u(u, v)==bez2.f_u(u, v));
1313  TEST_ASSERT(bez.f_v(u, v)==bez2.f_v(u, v));
1314  TEST_ASSERT(bez.f_uu(u, v)==bez2.f_uu(u, v));
1315  TEST_ASSERT(bez.f_uv(u, v)==bez2.f_uv(u, v));
1316  TEST_ASSERT(bez.f_vv(u, v)==bez2.f_vv(u, v));
1317  TEST_ASSERT(bez.f_uuu(u, v)==bez2.f_uuu(u, v));
1318  TEST_ASSERT(bez.f_uuv(u, v)==bez2.f_uuv(u, v));
1319  TEST_ASSERT(bez.f_uvv(u, v)==bez2.f_uvv(u, v));
1320  TEST_ASSERT(bez.f_vvv(u, v)==bez2.f_vvv(u, v));
1321 
1322  // test evaluation at interior point u=1/4, v=3/4
1323  u=0.25; v=0.75;
1324  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
1325  TEST_ASSERT(bez.f_u(u, v)==bez2.f_u(u, v));
1326  TEST_ASSERT(bez.f_v(u, v)==bez2.f_v(u, v));
1327  TEST_ASSERT(bez.f_uu(u, v)==bez2.f_uu(u, v));
1328  TEST_ASSERT(bez.f_uv(u, v)==bez2.f_uv(u, v));
1329  TEST_ASSERT(bez.f_vv(u, v)==bez2.f_vv(u, v));
1330  TEST_ASSERT(bez.f_uuu(u, v)==bez2.f_uuu(u, v));
1331  TEST_ASSERT(bez.f_uuv(u, v)==bez2.f_uuv(u, v));
1332  TEST_ASSERT(bez.f_uvv(u, v)==bez2.f_uvv(u, v));
1333  TEST_ASSERT(bez.f_vvv(u, v)==bez2.f_vvv(u, v));
1334  }
1335 
1336  // test promote in v-direction
1337  {
1338  bezier_type bez2(bez);
1339 
1340  // promote
1341  bez2.promote_v();
1342 
1343  // test if degree increased
1344  TEST_ASSERT(bez2.degree_v()==bez.degree_v()+1);
1345 
1346  // test evaluation at u=v=1/2
1347  u=0.5; v=0.5;
1348  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
1349  TEST_ASSERT(bez.f_u(u, v)==bez2.f_u(u, v));
1350  TEST_ASSERT(bez.f_v(u, v)==bez2.f_v(u, v));
1351  TEST_ASSERT(bez.f_uu(u, v)==bez2.f_uu(u, v));
1352  TEST_ASSERT(bez.f_uv(u, v)==bez2.f_uv(u, v));
1353  TEST_ASSERT(bez.f_vv(u, v)==bez2.f_vv(u, v));
1354  TEST_ASSERT(bez.f_uuu(u, v)==bez2.f_uuu(u, v));
1355  TEST_ASSERT(bez.f_uuv(u, v)==bez2.f_uuv(u, v));
1356  TEST_ASSERT(bez.f_uvv(u, v)==bez2.f_uvv(u, v));
1357  TEST_ASSERT(bez.f_vvv(u, v)==bez2.f_vvv(u, v));
1358 
1359  // test evaluation at interior point u=1/4, v=3/4
1360  u=0.25; v=0.75;
1361  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
1362  TEST_ASSERT(bez.f_u(u, v)==bez2.f_u(u, v));
1363  TEST_ASSERT(bez.f_v(u, v)==bez2.f_v(u, v));
1364  TEST_ASSERT(bez.f_uu(u, v)==bez2.f_uu(u, v));
1365  TEST_ASSERT(bez.f_uv(u, v)==bez2.f_uv(u, v));
1366  TEST_ASSERT(bez.f_vv(u, v)==bez2.f_vv(u, v));
1367  TEST_ASSERT(bez.f_uuu(u, v)==bez2.f_uuu(u, v));
1368  TEST_ASSERT(bez.f_uuv(u, v)==bez2.f_uuv(u, v));
1369  TEST_ASSERT(bez.f_uvv(u, v)==bez2.f_uvv(u, v));
1370  TEST_ASSERT(bez.f_vvv(u, v)==bez2.f_vvv(u, v));
1371  }
1372  }
1373 
1375  {
1376  index_type n(3), m(3);
1377  point_type pt[3+1][3+1];
1378  data_type u, v;
1379 
1380  // create surface with specified control points
1381  pt[0][0] << -15, 0, 15;
1382  pt[1][0] << -5, 5, 15;
1383  pt[2][0] << 5, 5, 15;
1384  pt[3][0] << 15, 0, 15;
1385  pt[0][1] << -15, 5, 5;
1386  pt[1][1] << -5, 5, 5;
1387  pt[2][1] << 5, 5, 5;
1388  pt[3][1] << 15, 5, 5;
1389  pt[0][2] << -15, 5, -5;
1390  pt[1][2] << -5, 5, -5;
1391  pt[2][2] << 5, 5, -5;
1392  pt[3][2] << 15, 5, -5;
1393  pt[0][3] << -15, 0, -15;
1394  pt[1][3] << -5, 5, -15;
1395  pt[2][3] << 5, 5, -15;
1396  pt[3][3] << 15, 0, -15;
1397 
1398  // create surface with specified dimensions and set control points
1399  bezier_type bez(n, m);
1400 
1401  for (index_type i=0; i<=n; ++i)
1402  {
1403  for (index_type j=0; j<=m; ++j)
1404  {
1405  bez.set_control_point(pt[i][j], i, j);
1406  }
1407  }
1408 
1409  // test promote in u-direction
1410  {
1411  bezier_type bez2(bez);
1412  bezier_type bez3(bez);
1413 
1414  // promote
1415  bez2.promote_u();
1416 
1417  // promote to +1
1418  bez3.promote_u_to(bez.degree_u()+1);
1419 
1420  // test if degree increased
1421  TEST_ASSERT(bez3.degree_u()==bez.degree_u()+1);
1422 
1423  for (index_type i=0; i<=n+1; ++i)
1424  {
1425  for (index_type j=0; j<=m; ++j)
1426  {
1427  TEST_ASSERT(bez2.get_control_point(i, j)==bez3.get_control_point(i, j));
1428  }
1429  }
1430 
1431  }
1432 
1433  // test promote in v-direction
1434  {
1435  bezier_type bez2(bez);
1436  bezier_type bez3(bez);
1437 
1438  // promote
1439  bez2.promote_v();
1440 
1441  bez3.promote_v_to(bez.degree_v()+1);
1442 
1443  // test if degree increased
1444  TEST_ASSERT(bez3.degree_v()==bez.degree_v()+1);
1445 
1446  for (index_type i=0; i<=n; ++i)
1447  {
1448  for (index_type j=0; j<=m+1; ++j)
1449  {
1450  TEST_ASSERT(bez2.get_control_point(i, j)==bez3.get_control_point(i, j));
1451  }
1452  }
1453  }
1454 
1455  // test promote in u-direction
1456  {
1457  bezier_type bez2(bez);
1458 
1459  index_type inc = 3;
1460 
1461  // promote
1462  bez2.promote_u_to(bez.degree_u()+inc);
1463 
1464  // test if degree increased
1465  TEST_ASSERT(bez2.degree_u()==bez.degree_u()+inc);
1466 
1467  // test evaluation at u=v=1/2
1468  u=0.5; v=0.5;
1469  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
1470  TEST_ASSERT(bez.f_u(u, v)==bez2.f_u(u, v));
1471  TEST_ASSERT(bez.f_v(u, v)==bez2.f_v(u, v));
1472  TEST_ASSERT(bez.f_uu(u, v)==bez2.f_uu(u, v));
1473  TEST_ASSERT(bez.f_uv(u, v)==bez2.f_uv(u, v));
1474  TEST_ASSERT(bez.f_vv(u, v)==bez2.f_vv(u, v));
1475  TEST_ASSERT(bez.f_uuu(u, v)==bez2.f_uuu(u, v));
1476  TEST_ASSERT(bez.f_uuv(u, v)==bez2.f_uuv(u, v));
1477  TEST_ASSERT(bez.f_uvv(u, v)==bez2.f_uvv(u, v));
1478  TEST_ASSERT(bez.f_vvv(u, v)==bez2.f_vvv(u, v));
1479 
1480  // test evaluation at interior point u=1/4, v=3/4
1481  u=0.25; v=0.75;
1482  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
1483  TEST_ASSERT(bez.f_u(u, v)==bez2.f_u(u, v));
1484  TEST_ASSERT(bez.f_v(u, v)==bez2.f_v(u, v));
1485  TEST_ASSERT(bez.f_uu(u, v)==bez2.f_uu(u, v));
1486  TEST_ASSERT(bez.f_uv(u, v)==bez2.f_uv(u, v));
1487  TEST_ASSERT(bez.f_vv(u, v)==bez2.f_vv(u, v));
1488  TEST_ASSERT(bez.f_uuu(u, v)==bez2.f_uuu(u, v));
1489  TEST_ASSERT(bez.f_uuv(u, v)==bez2.f_uuv(u, v));
1490  TEST_ASSERT(bez.f_uvv(u, v)==bez2.f_uvv(u, v));
1491  TEST_ASSERT(bez.f_vvv(u, v)==bez2.f_vvv(u, v));
1492  }
1493 
1494  // test promote in v-direction
1495  {
1496  bezier_type bez2(bez);
1497 
1498  index_type inc = 3;
1499 
1500  // promote
1501  bez2.promote_v_to(bez.degree_v()+inc);
1502 
1503  // test if degree increased
1504  TEST_ASSERT(bez2.degree_v()==bez.degree_v()+inc);
1505 
1506  // test evaluation at u=v=1/2
1507  u=0.5; v=0.5;
1508  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
1509  TEST_ASSERT(bez.f_u(u, v)==bez2.f_u(u, v));
1510  TEST_ASSERT(bez.f_v(u, v)==bez2.f_v(u, v));
1511  TEST_ASSERT(bez.f_uu(u, v)==bez2.f_uu(u, v));
1512  TEST_ASSERT(bez.f_uv(u, v)==bez2.f_uv(u, v));
1513  TEST_ASSERT(bez.f_vv(u, v)==bez2.f_vv(u, v));
1514  TEST_ASSERT(bez.f_uuu(u, v)==bez2.f_uuu(u, v));
1515  TEST_ASSERT(bez.f_uuv(u, v)==bez2.f_uuv(u, v));
1516  TEST_ASSERT(bez.f_uvv(u, v)==bez2.f_uvv(u, v));
1517  TEST_ASSERT(bez.f_vvv(u, v)==bez2.f_vvv(u, v));
1518 
1519  // test evaluation at interior point u=1/4, v=3/4
1520  u=0.25; v=0.75;
1521  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
1522  TEST_ASSERT(bez.f_u(u, v)==bez2.f_u(u, v));
1523  TEST_ASSERT(bez.f_v(u, v)==bez2.f_v(u, v));
1524  TEST_ASSERT(bez.f_uu(u, v)==bez2.f_uu(u, v));
1525  TEST_ASSERT(bez.f_uv(u, v)==bez2.f_uv(u, v));
1526  TEST_ASSERT(bez.f_vv(u, v)==bez2.f_vv(u, v));
1527  TEST_ASSERT(bez.f_uuu(u, v)==bez2.f_uuu(u, v));
1528  TEST_ASSERT(bez.f_uuv(u, v)==bez2.f_uuv(u, v));
1529  TEST_ASSERT(bez.f_uvv(u, v)==bez2.f_uvv(u, v));
1530  TEST_ASSERT(bez.f_vvv(u, v)==bez2.f_vvv(u, v));
1531  }
1532  }
1533 
1535  {
1536  // hand checked with "Degree Reduction of Bezier Curves" by Dave Morgan
1537  {
1538  index_type n(6), m(4);
1539  point_type pt[6+1][4+1], pt_ref[5+1][4+1];
1540  bool demoted;
1541 
1542  // create surface with specified control points
1543  pt[0][0] << 0, 0, 15;
1544  pt[1][0] << 2, 6, 15;
1545  pt[2][0] << 3, 0, 15;
1546  pt[3][0] << 5, 4, 15;
1547  pt[4][0] << 7, 1, 15;
1548  pt[5][0] << 5, 5, 15;
1549  pt[6][0] << 10, 6, 15;
1550  pt[0][1] << 0, 0, 11;
1551  pt[1][1] << 2, 6, 11;
1552  pt[2][1] << 3, 0, 11;
1553  pt[3][1] << 5, 4, 11;
1554  pt[4][1] << 7, 1, 11;
1555  pt[5][1] << 5, 5, 11;
1556  pt[6][1] << 10, 6, 11;
1557  pt[0][2] << 0, 0, 3;
1558  pt[1][2] << 2, 6, 3;
1559  pt[2][2] << 3, 0, 3;
1560  pt[3][2] << 5, 4, 3;
1561  pt[4][2] << 7, 1, 3;
1562  pt[5][2] << 5, 5, 3;
1563  pt[6][2] << 10, 6, 3;
1564  pt[0][3] << 0, 0, 0;
1565  pt[1][3] << 2, 6, 0;
1566  pt[2][3] << 3, 0, 0;
1567  pt[3][3] << 5, 4, 0;
1568  pt[4][3] << 7, 1, 0;
1569  pt[5][3] << 5, 5, 0;
1570  pt[6][3] << 10, 6, 0;
1571  pt[0][4] << 0, 0, -5;
1572  pt[1][4] << 2, 6, -5;
1573  pt[2][4] << 3, 0, -5;
1574  pt[3][4] << 5, 4, -5;
1575  pt[4][4] << 7, 1, -5;
1576  pt[5][4] << 5, 5, -5;
1577  pt[6][4] << 10, 6, -5;
1578 
1579  // create surface with specified dimensions and set control points
1580  bezier_type bez(n, m), bez2;
1581 
1582  for (index_type i=0; i<=n; ++i)
1583  {
1584  for (index_type j=0; j<=m; ++j)
1585  {
1586  bez.set_control_point(pt[i][j], i, j);
1587  }
1588  }
1589  bez2=bez;
1590 
1592 
1593  TEST_ASSERT(demoted);
1594 
1595  // test if the degree is correct
1596  TEST_ASSERT(bez.degree_u()==bez2.degree_u()+1);
1597  TEST_ASSERT(bez.degree_v()==bez2.degree_v());
1598 
1599  // set the reference control points
1600  pt_ref[0][0] << static_cast<data_type>(-0.00878906), static_cast<data_type>(0.0610352), 15;
1601  pt_ref[1][0] << static_cast<data_type>( 2.5177734), static_cast<data_type>( 6.3821289), 15;
1602  pt_ref[2][0] << static_cast<data_type>( 2.8060547), static_cast<data_type>(-0.16982422), 15;
1603  pt_ref[3][0] << static_cast<data_type>( 8.0060547), static_cast<data_type>( 2.5301758), 15;
1604  pt_ref[4][0] << static_cast<data_type>( 4.1177734), static_cast<data_type>( 3.9821289), 15;
1605  pt_ref[5][0] << static_cast<data_type>( 9.9912109), static_cast<data_type>( 6.0610352), 15;
1606  pt_ref[0][1] << static_cast<data_type>(-0.00878906), static_cast<data_type>( 0.0610352), 11;
1607  pt_ref[1][1] << static_cast<data_type>( 2.5177734), static_cast<data_type>( 6.3821289), 11;
1608  pt_ref[2][1] << static_cast<data_type>( 2.8060547), static_cast<data_type>(-0.16982422), 11;
1609  pt_ref[3][1] << static_cast<data_type>( 8.0060547), static_cast<data_type>( 2.5301758), 11;
1610  pt_ref[4][1] << static_cast<data_type>( 4.1177734), static_cast<data_type>( 3.9821289), 11;
1611  pt_ref[5][1] << static_cast<data_type>( 9.9912109), static_cast<data_type>( 6.0610352), 11;
1612  pt_ref[0][2] << static_cast<data_type>(-0.00878906), static_cast<data_type>( 0.0610352), 3;
1613  pt_ref[1][2] << static_cast<data_type>( 2.5177734), static_cast<data_type>( 6.3821289), 3;
1614  pt_ref[2][2] << static_cast<data_type>( 2.8060547), static_cast<data_type>(-0.16982422), 3;
1615  pt_ref[3][2] << static_cast<data_type>( 8.0060547), static_cast<data_type>( 2.5301758), 3;
1616  pt_ref[4][2] << static_cast<data_type>( 4.1177734), static_cast<data_type>( 3.9821289), 3;
1617  pt_ref[5][2] << static_cast<data_type>( 9.9912109), static_cast<data_type>( 6.0610352), 3;
1618  pt_ref[0][3] << static_cast<data_type>(-0.00878906), static_cast<data_type>( 0.0610352), 0;
1619  pt_ref[1][3] << static_cast<data_type>( 2.5177734), static_cast<data_type>( 6.3821289), 0;
1620  pt_ref[2][3] << static_cast<data_type>( 2.8060547), static_cast<data_type>(-0.16982422), 0;
1621  pt_ref[3][3] << static_cast<data_type>( 8.0060547), static_cast<data_type>( 2.5301758), 0;
1622  pt_ref[4][3] << static_cast<data_type>( 4.1177734), static_cast<data_type>( 3.9821289), 0;
1623  pt_ref[5][3] << static_cast<data_type>( 9.9912109), static_cast<data_type>( 6.0610352), 0;
1624  pt_ref[0][4] << static_cast<data_type>(-0.00878906), static_cast<data_type>( 0.0610352), -5;
1625  pt_ref[1][4] << static_cast<data_type>( 2.5177734), static_cast<data_type>( 6.3821289), -5;
1626  pt_ref[2][4] << static_cast<data_type>( 2.8060547), static_cast<data_type>(-0.16982422), -5;
1627  pt_ref[3][4] << static_cast<data_type>( 8.0060547), static_cast<data_type>( 2.5301758), -5;
1628  pt_ref[4][4] << static_cast<data_type>( 4.1177734), static_cast<data_type>( 3.9821289), -5;
1629  pt_ref[5][4] << static_cast<data_type>( 9.9912109), static_cast<data_type>( 6.0610352), -5;
1630 
1631  // test if get the correct control points
1632  for (index_type j=0; j<=m; ++j)
1633  {
1634  for (index_type i=0; i<=n-1; ++i)
1635  {
1636  TEST_ASSERT((bez2.get_control_point(i,j)-pt_ref[i][j]).norm()<1e-6);
1637  }
1638  }
1639  }
1640 
1641  // test whether can promote and demote and get the same control points
1642  {
1643  index_type n(3), m(3);
1644  point_type pt[3+1][3+1];
1645  bool demoted;
1646 
1647  // create surface with specified control points
1648  pt[0][0] << -15, 0, 15;
1649  pt[1][0] << -5, 5, 15;
1650  pt[2][0] << 5, 5, 15;
1651  pt[3][0] << 15, 0, 15;
1652  pt[0][1] << -15, 5, 5;
1653  pt[1][1] << -5, 5, 5;
1654  pt[2][1] << 5, 5, 5;
1655  pt[3][1] << 15, 5, 5;
1656  pt[0][2] << -15, 5, -5;
1657  pt[1][2] << -5, 5, -5;
1658  pt[2][2] << 5, 5, -5;
1659  pt[3][2] << 15, 5, -5;
1660  pt[0][3] << -15, 0, -15;
1661  pt[1][3] << -5, 5, -15;
1662  pt[2][3] << 5, 5, -15;
1663  pt[3][3] << 15, 0, -15;
1664 
1665  // create surface with specified dimensions and set control points
1666  bezier_type bez(n, m), bez2;
1667 
1668  for (index_type i=0; i<=n; ++i)
1669  {
1670  for (index_type j=0; j<=m; ++j)
1671  {
1672  bez.set_control_point(pt[i][j], i, j);
1673  }
1674  }
1675  bez2=bez;
1676 
1677  // promote -> demote u-direction
1678  bez2.promote_u();
1680 
1681  TEST_ASSERT(demoted);
1682 
1683  // test to see if degree has changed
1684  TEST_ASSERT(bez.degree_u()==bez2.degree_u());
1685  TEST_ASSERT(bez.degree_v()==bez2.degree_v());
1686 
1687  // test to see if get expected polygon
1688  for (index_type i=0; i<=n; ++i)
1689  {
1690  for (index_type j=0; j<=m; ++j)
1691  {
1692  TEST_ASSERT(bez2.get_control_point(i, j)==bez.get_control_point(i, j));
1693  }
1694  }
1695 
1696  // promote -> demote v-direction
1697  bez2=bez;
1698  bez2.promote_v();
1700 
1701  TEST_ASSERT(demoted);
1702 
1703  // test to see if degree has changed
1704  TEST_ASSERT(bez.degree_u()==bez2.degree_u());
1705  TEST_ASSERT(bez.degree_v()==bez2.degree_v());
1706 
1707  // test to see if get expected polygon
1708  for (index_type i=0; i<=n; ++i)
1709  {
1710  for (index_type j=0; j<=m; ++j)
1711  {
1712  TEST_ASSERT(bez2.get_control_point(i, j)==bez.get_control_point(i, j));
1713  }
1714  }
1715  }
1716 
1717  // test whether can demote and maintain continuity at end points
1718  {
1719  index_type n(6), m(4);
1720  point_type pt[6+1][4+1];
1721 
1722  // create surface with specified control points
1723  pt[0][0] << 0, 0, 15;
1724  pt[1][0] << 2, 6, 15;
1725  pt[2][0] << 3, 0, 15;
1726  pt[3][0] << 5, 4, 15;
1727  pt[4][0] << 7, 1, 15;
1728  pt[5][0] << 5, 5, 15;
1729  pt[6][0] << 10, 6, 15;
1730  pt[0][1] << 0, 0, 11;
1731  pt[1][1] << 2, 6, 11;
1732  pt[2][1] << 3, 0, 11;
1733  pt[3][1] << 5, 4, 11;
1734  pt[4][1] << 7, 1, 11;
1735  pt[5][1] << 5, 5, 11;
1736  pt[6][1] << 10, 6, 11;
1737  pt[0][2] << 0, 0, 3;
1738  pt[1][2] << 2, 6, 3;
1739  pt[2][2] << 3, 0, 3;
1740  pt[3][2] << 5, 4, 3;
1741  pt[4][2] << 7, 1, 3;
1742  pt[5][2] << 5, 5, 3;
1743  pt[6][2] << 10, 6, 3;
1744  pt[0][3] << 0, 0, 0;
1745  pt[1][3] << 2, 6, 0;
1746  pt[2][3] << 3, 0, 0;
1747  pt[3][3] << 5, 4, 0;
1748  pt[4][3] << 7, 1, 0;
1749  pt[5][3] << 5, 5, 0;
1750  pt[6][3] << 10, 6, 0;
1751  pt[0][4] << 0, 0, -5;
1752  pt[1][4] << 2, 6, -5;
1753  pt[2][4] << 3, 0, -5;
1754  pt[3][4] << 5, 4, -5;
1755  pt[4][4] << 7, 1, -5;
1756  pt[5][4] << 5, 5, -5;
1757  pt[6][4] << 10, 6, -5;
1758 
1759  // create surface with specified dimensions and set control points
1760  bezier_type bez(n, m);
1761 
1762  for (index_type i=0; i<=n; ++i)
1763  {
1764  for (index_type j=0; j<=m; ++j)
1765  {
1766  bez.set_control_point(pt[i][j], i, j);
1767  }
1768  }
1769 
1770  // No End Constraint
1771  {
1772  bezier_type bez2(bez);
1773  point_type eval_out[6], eval_ref[6];
1774  data_type d[6], u[6]={0, 1, 0, 1, 0.5, 0.25}, v[6]={0, 0, 1, 1, 0.5, 0.75};
1775  bool demoted;
1776 
1777  // demote u-direction
1779  TEST_ASSERT(demoted);
1780  TEST_ASSERT(bez.degree_u()==bez2.degree_u()+1);
1781  TEST_ASSERT(bez.degree_v()==bez2.degree_v());
1782 
1783  for (size_t k=0; k<6; ++k)
1784  {
1785  eval_out[0]=bez2.f(u[k], v[k]);
1786  eval_out[1]=bez2.f_u(u[k], v[k]);
1787  eval_out[2]=bez2.f_v(u[k], v[k]);
1788  eval_out[3]=bez2.f_uu(u[k], v[k]);
1789  eval_out[4]=bez2.f_uv(u[k], v[k]);
1790  eval_out[5]=bez2.f_vv(u[k], v[k]);
1791  eval_ref[0]=bez.f(u[k], v[k]);
1792  eval_ref[1]=bez.f_u(u[k], v[k]);
1793  eval_ref[2]=bez.f_v(u[k], v[k]);
1794  eval_ref[3]=bez.f_uu(u[k], v[k]);
1795  eval_ref[4]=bez.f_uv(u[k], v[k]);
1796  eval_ref[5]=bez.f_vv(u[k], v[k]);
1797  for (index_type i=0; i<6; ++i)
1798  {
1799  d[i]=eli::geom::point::distance(eval_out[i], eval_ref[i]);
1800  }
1801 
1802  // boundaries have different error than interior points
1803  if (k<4)
1804  {
1805  TEST_ASSERT_DELTA(d[0], 0.061665, 5e-6);
1806  TEST_ASSERT_DELTA(d[1], 4.43986, 5e-5);
1807  TEST_ASSERT_DELTA(d[2], 0, std::numeric_limits<data_type>::epsilon());
1808  TEST_ASSERT_DELTA(d[3], 103.597, 5e-3);
1809  TEST_ASSERT_DELTA(d[4], 0, std::numeric_limits<data_type>::epsilon());
1810  TEST_ASSERT_DELTA(d[5], 0, std::numeric_limits<data_type>::epsilon());
1811  }
1812  else if (k==4)
1813  {
1814  TEST_ASSERT_DELTA(d[0], 0.061665, 5e-6);
1815  TEST_ASSERT_DELTA(d[1], 0, 20*std::numeric_limits<data_type>::epsilon());
1816  TEST_ASSERT_DELTA(d[2], 0, std::numeric_limits<data_type>::epsilon());
1817  TEST_ASSERT_DELTA(d[3], 8.8797199432, 5e-6);
1818  TEST_ASSERT_DELTA(d[4], 0, std::numeric_limits<data_type>::epsilon());
1819  TEST_ASSERT_DELTA(d[5], 0, std::numeric_limits<data_type>::epsilon());
1820  }
1821  else
1822  {
1823  TEST_ASSERT_DELTA(d[0], 0.061665, 5e-6);
1824  TEST_ASSERT_DELTA(d[1], 0, 20*std::numeric_limits<data_type>::epsilon());
1825  TEST_ASSERT_DELTA(d[2], 0, std::numeric_limits<data_type>::epsilon());
1826  TEST_ASSERT_DELTA(d[3], 11.83962659, 5e-6);
1827  TEST_ASSERT_DELTA(d[4], 0, std::numeric_limits<data_type>::epsilon());
1828  TEST_ASSERT_DELTA(d[5], 0, std::numeric_limits<data_type>::epsilon());
1829  }
1830  }
1831 
1832  // demote v-direction
1833  bez2=bez;
1835  TEST_ASSERT(demoted);
1836  TEST_ASSERT(bez.degree_u()==bez2.degree_u());
1837  TEST_ASSERT(bez.degree_v()==bez2.degree_v()+1);
1838 
1839  for (size_t k=0; k<6; ++k)
1840  {
1841  eval_out[0]=bez2.f(u[k], v[k]);
1842  eval_out[1]=bez2.f_u(u[k], v[k]);
1843  eval_out[2]=bez2.f_v(u[k], v[k]);
1844  eval_out[3]=bez2.f_uu(u[k], v[k]);
1845  eval_out[4]=bez2.f_uv(u[k], v[k]);
1846  eval_out[5]=bez2.f_vv(u[k], v[k]);
1847  eval_ref[0]=bez.f(u[k], v[k]);
1848  eval_ref[1]=bez.f_u(u[k], v[k]);
1849  eval_ref[2]=bez.f_v(u[k], v[k]);
1850  eval_ref[3]=bez.f_uu(u[k], v[k]);
1851  eval_ref[4]=bez.f_uv(u[k], v[k]);
1852  eval_ref[5]=bez.f_vv(u[k], v[k]);
1853  for (index_type i=0; i<6; ++i)
1854  {
1855  d[i]=eli::geom::point::distance(eval_out[i], eval_ref[i]);
1856  }
1857 
1858  // boundaries have different error than interior points
1859  if (k<4)
1860  {
1861  TEST_ASSERT_DELTA(d[0], 0.125, 1e-6);
1862  TEST_ASSERT_DELTA(d[1], 0, std::numeric_limits<data_type>::epsilon());
1863  TEST_ASSERT_DELTA(d[2], 4, 1e-5);
1864  TEST_ASSERT_DELTA(d[3], 0, std::numeric_limits<data_type>::epsilon());
1865  TEST_ASSERT_DELTA(d[4], 0, std::numeric_limits<data_type>::epsilon());
1866  TEST_ASSERT_DELTA(d[5], 40, 1e-5);
1867  }
1868  else if (k==4)
1869  {
1870  TEST_ASSERT_DELTA(d[0], 0.125, 1e-6);
1871  TEST_ASSERT_DELTA(d[1], 0, std::numeric_limits<data_type>::epsilon());
1872  TEST_ASSERT_DELTA(d[2], 0, std::numeric_limits<data_type>::epsilon());
1873  TEST_ASSERT_DELTA(d[3], 0, std::numeric_limits<data_type>::epsilon());
1874  TEST_ASSERT_DELTA(d[4], 0, std::numeric_limits<data_type>::epsilon());
1875  TEST_ASSERT_DELTA(d[5], 8, 1e-6);
1876  }
1877  else
1878  {
1879  TEST_ASSERT_DELTA(d[0], 0.0625, 1e-6);
1880  TEST_ASSERT_DELTA(d[1], 0, std::numeric_limits<data_type>::epsilon());
1881  TEST_ASSERT_DELTA(d[2], 1, 1e-6);
1882  TEST_ASSERT_DELTA(d[3], 0, std::numeric_limits<data_type>::epsilon());
1883  TEST_ASSERT_DELTA(d[4], 0, std::numeric_limits<data_type>::epsilon());
1884  TEST_ASSERT_DELTA(d[5], 4, 2e-6);
1885  }
1886  }
1887  }
1888 
1889  // Point End Constraint
1890  {
1891  bezier_type bez2(bez);
1892  point_type eval_out[6], eval_ref[6];
1893  data_type d[6], u[6]={0, 1, 0, 1, 0.5, 0.25}, v[6]={0, 0, 1, 1, 0.5, 0.75};
1894  bool demoted;
1895 
1896  // demote u-direction
1897  demoted=bez2.demote_u(eli::geom::general::C0);
1898  TEST_ASSERT(demoted);
1899  TEST_ASSERT(bez.degree_u()==bez2.degree_u()+1);
1900  TEST_ASSERT(bez.degree_v()==bez2.degree_v());
1901 
1902  for (size_t k=0; k<6; ++k)
1903  {
1904  eval_out[0]=bez2.f(u[k], v[k]);
1905  eval_out[1]=bez2.f_u(u[k], v[k]);
1906  eval_out[2]=bez2.f_v(u[k], v[k]);
1907  eval_out[3]=bez2.f_uu(u[k], v[k]);
1908  eval_out[4]=bez2.f_uv(u[k], v[k]);
1909  eval_out[5]=bez2.f_vv(u[k], v[k]);
1910  eval_ref[0]=bez.f(u[k], v[k]);
1911  eval_ref[1]=bez.f_u(u[k], v[k]);
1912  eval_ref[2]=bez.f_v(u[k], v[k]);
1913  eval_ref[3]=bez.f_uu(u[k], v[k]);
1914  eval_ref[4]=bez.f_uv(u[k], v[k]);
1915  eval_ref[5]=bez.f_vv(u[k], v[k]);
1916  for (index_type i=0; i<6; ++i)
1917  {
1918  d[i]=eli::geom::point::distance(eval_out[i], eval_ref[i]);
1919  }
1920 
1921  // boundaries have different error than interior points
1922  if (k<4)
1923  {
1924  TEST_ASSERT_DELTA(d[0], 0, 61*std::numeric_limits<data_type>::epsilon());
1925  TEST_ASSERT_DELTA(d[1], 3.82695, 1e-5);
1926  TEST_ASSERT_DELTA(d[2], 0, std::numeric_limits<data_type>::epsilon());
1927  TEST_ASSERT_DELTA(d[3], 99.5007, 1e-4);
1928  TEST_ASSERT_DELTA(d[4], 0, 41*std::numeric_limits<data_type>::epsilon());
1929  TEST_ASSERT_DELTA(d[5], 0, std::numeric_limits<data_type>::epsilon());
1930  }
1931  else if (k==4)
1932  {
1933  TEST_ASSERT_DELTA(d[0], 0.0597961, 1e-6);
1934  TEST_ASSERT_DELTA(d[1], 0, 20*std::numeric_limits<data_type>::epsilon());
1935  TEST_ASSERT_DELTA(d[2], 0, std::numeric_limits<data_type>::epsilon());
1936  TEST_ASSERT_DELTA(d[3], 9.08901, 1e-5);
1937  TEST_ASSERT_DELTA(d[4], 0, std::numeric_limits<data_type>::epsilon());
1938  TEST_ASSERT_DELTA(d[5], 0, std::numeric_limits<data_type>::epsilon());
1939  }
1940  else
1941  {
1942  TEST_ASSERT_DELTA(d[0], 0.0644677, 1e-6);
1943  TEST_ASSERT_DELTA(d[1], 0.0373726, 1e-6);
1944  TEST_ASSERT_DELTA(d[2], 0, std::numeric_limits<data_type>::epsilon());
1945  TEST_ASSERT_DELTA(d[3], 12.7067, 1e-4);
1946  TEST_ASSERT_DELTA(d[4], 0, std::numeric_limits<data_type>::epsilon());
1947  TEST_ASSERT_DELTA(d[5], 0, 9*std::numeric_limits<data_type>::epsilon());
1948  }
1949  }
1950 
1951  // demote v-direction
1952  bez2=bez;
1953  demoted=bez2.demote_v(eli::geom::general::C0);
1954  TEST_ASSERT(demoted);
1955  TEST_ASSERT(bez.degree_u()==bez2.degree_u());
1956  TEST_ASSERT(bez.degree_v()==bez2.degree_v()+1);
1957 
1958  for (size_t k=0; k<6; ++k)
1959  {
1960  eval_out[0]=bez2.f(u[k], v[k]);
1961  eval_out[1]=bez2.f_u(u[k], v[k]);
1962  eval_out[2]=bez2.f_v(u[k], v[k]);
1963  eval_out[3]=bez2.f_uu(u[k], v[k]);
1964  eval_out[4]=bez2.f_uv(u[k], v[k]);
1965  eval_out[5]=bez2.f_vv(u[k], v[k]);
1966  eval_ref[0]=bez.f(u[k], v[k]);
1967  eval_ref[1]=bez.f_u(u[k], v[k]);
1968  eval_ref[2]=bez.f_v(u[k], v[k]);
1969  eval_ref[3]=bez.f_uu(u[k], v[k]);
1970  eval_ref[4]=bez.f_uv(u[k], v[k]);
1971  eval_ref[5]=bez.f_vv(u[k], v[k]);
1972  for (index_type i=0; i<6; ++i)
1973  {
1974  d[i]=eli::geom::point::distance(eval_out[i], eval_ref[i]);
1975  }
1976 
1977  // boundaries have different error than interior points
1978  if (k<4)
1979  {
1980  TEST_ASSERT_DELTA(d[0], 0, std::numeric_limits<data_type>::epsilon());
1981  TEST_ASSERT_DELTA(d[1], 0, std::numeric_limits<data_type>::epsilon());
1982  TEST_ASSERT_DELTA(d[2], 3.428571, 1e-5);
1983  TEST_ASSERT_DELTA(d[3], 0, std::numeric_limits<data_type>::epsilon());
1984  TEST_ASSERT_DELTA(d[4], 0, std::numeric_limits<data_type>::epsilon());
1985  TEST_ASSERT_DELTA(d[5], 38.8571, 1e-4);
1986  }
1987  else if (k==4)
1988  {
1989  TEST_ASSERT_DELTA(d[0], 0.1428571, 1e-6);
1990  TEST_ASSERT_DELTA(d[1], 0, std::numeric_limits<data_type>::epsilon());
1991  TEST_ASSERT_DELTA(d[2], 0, 17*std::numeric_limits<data_type>::epsilon());
1992  TEST_ASSERT_DELTA(d[3], 0, std::numeric_limits<data_type>::epsilon());
1993  TEST_ASSERT_DELTA(d[4], 0, std::numeric_limits<data_type>::epsilon());
1994  TEST_ASSERT_DELTA(d[5], 9.142857, 1e-5);
1995  }
1996  else
1997  {
1998  TEST_ASSERT_DELTA(d[0], 0.0803571, 1e-6);
1999  TEST_ASSERT_DELTA(d[1], 0, std::numeric_limits<data_type>::epsilon());
2000  TEST_ASSERT_DELTA(d[2], 1.285714, 1e-6);
2001  TEST_ASSERT_DELTA(d[3], 0, std::numeric_limits<data_type>::epsilon());
2002  TEST_ASSERT_DELTA(d[4], 0, std::numeric_limits<data_type>::epsilon());
2003  TEST_ASSERT_DELTA(d[5], 2.857143, 1e-5);
2004  }
2005  }
2006  }
2007 
2008  // Slope End Constraint
2009  {
2010  bezier_type bez2(bez);
2011  point_type eval_out[6], eval_ref[6];
2012  data_type d[6], u[6]={0, 1, 0, 1, 0.5, 0.25}, v[6]={0, 0, 1, 1, 0.5, 0.75};
2013  bool demoted;
2014 
2015  // demote u-direction
2016  demoted=bez2.demote_u(eli::geom::general::C1);
2017  TEST_ASSERT(demoted);
2018  TEST_ASSERT(bez.degree_u()==bez2.degree_u()+1);
2019  TEST_ASSERT(bez.degree_v()==bez2.degree_v());
2020 
2021  for (size_t k=0; k<6; ++k)
2022  {
2023  eval_out[0]=bez2.f(u[k], v[k]);
2024  eval_out[1]=bez2.f_u(u[k], v[k]);
2025  eval_out[2]=bez2.f_v(u[k], v[k]);
2026  eval_out[3]=bez2.f_uu(u[k], v[k]);
2027  eval_out[4]=bez2.f_uv(u[k], v[k]);
2028  eval_out[5]=bez2.f_vv(u[k], v[k]);
2029  eval_ref[0]=bez.f(u[k], v[k]);
2030  eval_ref[1]=bez.f_u(u[k], v[k]);
2031  eval_ref[2]=bez.f_v(u[k], v[k]);
2032  eval_ref[3]=bez.f_uu(u[k], v[k]);
2033  eval_ref[4]=bez.f_uv(u[k], v[k]);
2034  eval_ref[5]=bez.f_vv(u[k], v[k]);
2035  for (index_type i=0; i<6; ++i)
2036  {
2037  d[i]=eli::geom::point::distance(eval_out[i], eval_ref[i]);
2038  }
2039 
2040  // boundaries have different error than interior points
2041  if (k<4)
2042  {
2043  TEST_ASSERT_DELTA(d[0], 0, std::numeric_limits<data_type>::epsilon());
2044  TEST_ASSERT_DELTA(d[1], 0, 20*std::numeric_limits<data_type>::epsilon());
2045  TEST_ASSERT_DELTA(d[2], 0, std::numeric_limits<data_type>::epsilon());
2046  TEST_ASSERT_DELTA(d[3], 57.40425, 1e-4);
2047  TEST_ASSERT_DELTA(d[4], 0, std::numeric_limits<data_type>::epsilon());
2048  TEST_ASSERT_DELTA(d[5], 0, std::numeric_limits<data_type>::epsilon());
2049  }
2050  else if (k==4)
2051  {
2052  TEST_ASSERT_DELTA(d[0], 0.179388, 1e-6);
2053  TEST_ASSERT_DELTA(d[1], 0, 20*std::numeric_limits<data_type>::epsilon());
2054  TEST_ASSERT_DELTA(d[2], 0, std::numeric_limits<data_type>::epsilon());
2055  TEST_ASSERT_DELTA(d[3], 18.65638, 1e-5);
2056  TEST_ASSERT_DELTA(d[4], 0, std::numeric_limits<data_type>::epsilon());
2057  TEST_ASSERT_DELTA(d[5], 0, std::numeric_limits<data_type>::epsilon());
2058  }
2059  else
2060  {
2061  TEST_ASSERT_DELTA(d[0], 0.1765853, 1e-6);
2062  TEST_ASSERT_DELTA(d[1], 1.278142, 1e-5);
2063  TEST_ASSERT_DELTA(d[2], 0, std::numeric_limits<data_type>::epsilon());
2064  TEST_ASSERT_DELTA(d[3], 16.05525, 1e-5);
2065  TEST_ASSERT_DELTA(d[4], 0, std::numeric_limits<data_type>::epsilon());
2066  TEST_ASSERT_DELTA(d[5], 0, std::numeric_limits<data_type>::epsilon());
2067  }
2068  }
2069 
2070  // demote v-direction (cannot do, degree is to low)
2071  bez2=bez;
2072  demoted=bez2.demote_v(eli::geom::general::C1);
2073  TEST_ASSERT(!demoted);
2074  }
2075 
2076  // Second Derivative Constraint
2077  {
2078  bezier_type bez2(bez);
2079  bool demoted;
2080 
2081  // demote u-direction (cannot do, degree is to low)
2082  demoted=bez2.demote_u(eli::geom::general::C2);
2083  TEST_ASSERT(!demoted);
2084 
2085  // demote v-direction (cannot do, degree is to low)
2086  bez2=bez;
2087  demoted=bez2.demote_v(eli::geom::general::C2);
2088  TEST_ASSERT(!demoted);
2089  }
2090  }
2091  }
2092 
2093 
2095  {
2096  data_type eps(std::numeric_limits<data__>::epsilon());
2097  index_type n(6), m(4);
2098  point_type pt[6+1][4+1];
2099 
2100  // create surface with specified control points
2101  pt[0][0] << 0, 0, 15;
2102  pt[1][0] << 2, 6, 15;
2103  pt[2][0] << 3, 0, 15;
2104  pt[3][0] << 5, 4, 15;
2105  pt[4][0] << 7, 1, 15;
2106  pt[5][0] << 5, 5, 15;
2107  pt[6][0] << 10, 6, 15;
2108  pt[0][1] << 0, 0, 11;
2109  pt[1][1] << 2, 6, 11;
2110  pt[2][1] << 3, 0, 11;
2111  pt[3][1] << 5, 4, 11;
2112  pt[4][1] << 7, 1, 11;
2113  pt[5][1] << 5, 5, 11;
2114  pt[6][1] << 10, 6, 11;
2115  pt[0][2] << 0, 0, 3;
2116  pt[1][2] << 2, 6, 3;
2117  pt[2][2] << 3, 0, 3;
2118  pt[3][2] << 5, 4, 3;
2119  pt[4][2] << 7, 1, 3;
2120  pt[5][2] << 5, 5, 3;
2121  pt[6][2] << 10, 6, 3;
2122  pt[0][3] << 0, 0, 0;
2123  pt[1][3] << 2, 6, 0;
2124  pt[2][3] << 3, 0, 0;
2125  pt[3][3] << 5, 4, 0;
2126  pt[4][3] << 7, 1, 0;
2127  pt[5][3] << 5, 5, 0;
2128  pt[6][3] << 10, 6, 0;
2129  pt[0][4] << 0, 0, -5;
2130  pt[1][4] << 2, 6, -5;
2131  pt[2][4] << 3, 0, -5;
2132  pt[3][4] << 5, 4, -5;
2133  pt[4][4] << 7, 1, -5;
2134  pt[5][4] << 5, 5, -5;
2135  pt[6][4] << 10, 6, -5;
2136 
2137  // create surface with specified dimensions and set control points
2138  bezier_type bez(n, m);
2139 
2140  for (index_type i=0; i<=n; ++i)
2141  {
2142  for (index_type j=0; j<=m; ++j)
2143  {
2144  bez.set_control_point(pt[i][j], i, j);
2145  }
2146  }
2147 
2148  {
2149  bezier_type bez2(bez);
2150 
2151  bez2.to_cubic_u();
2152 
2153  TEST_ASSERT(bez2.degree_u()==3);
2154 
2155  data_type u, v;
2156 
2157  // test evaluation and first derivative at u endpoints
2158  u=0.0;
2159  v=0.0;
2160  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2161  TEST_ASSERT(bez.f_u(u, v)==bez2.f_u(u, v));
2162 
2163  v=static_cast<data_type>(0.2321);
2164  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2165  TEST_ASSERT(bez.f_u(u, v)==bez2.f_u(u, v));
2166 
2167  v=0.5;
2168  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2169  TEST_ASSERT(bez.f_u(u, v)==bez2.f_u(u, v));
2170 
2171  v=1.0;
2172  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2173  TEST_ASSERT(bez.f_u(u, v)==bez2.f_u(u, v));
2174 
2175  u=1.0;
2176  v=0.0;
2177  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2178  TEST_ASSERT(bez.f_u(u, v)==bez2.f_u(u, v));
2179 
2180  v=static_cast<data_type>(0.2321);
2181  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2182  TEST_ASSERT((bez.f_u(u, v)-bez2.f_u(u, v)).norm() < 50*eps);
2183 
2184  v=0.5;
2185  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2186  TEST_ASSERT(bez.f_u(u, v)==bez2.f_u(u, v));
2187 
2188  v=1.0;
2189  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2190  TEST_ASSERT(bez.f_u(u, v)==bez2.f_u(u, v));
2191  }
2192 
2193  {
2194  bezier_type bez2(bez);
2195 
2196  bez2.to_cubic_v();
2197 
2198  TEST_ASSERT(bez2.degree_v()==3);
2199 
2200  data_type u, v;
2201 
2202  // test evaluation and first derivative at v endpoints
2203  v=0.0;
2204  u=0.0;
2205  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2206  TEST_ASSERT((bez.f_v(u, v)-bez2.f_v(u, v)).norm() < 50*eps);
2207 
2208  u=static_cast<data_type>(0.2321);
2209  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2210  TEST_ASSERT((bez.f_v(u, v)-bez2.f_v(u, v)).norm() < 50*eps);
2211 
2212  u=0.5;
2213  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2214  TEST_ASSERT((bez.f_v(u, v)-bez2.f_v(u, v)).norm() < 50*eps);
2215 
2216  u=1.0;
2217  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2218  TEST_ASSERT((bez.f_v(u, v)-bez2.f_v(u, v)).norm() < 50*eps);
2219 
2220  v=1.0;
2221  u=0.0;
2222  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2223  TEST_ASSERT(bez.f_v(u, v)==bez2.f_v(u, v));
2224 
2225  u=static_cast<data_type>(0.2321);
2226  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2227  TEST_ASSERT(bez.f_v(u, v)==bez2.f_v(u, v));
2228 
2229  u=0.5;
2230  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2231  TEST_ASSERT(bez.f_v(u, v)==bez2.f_v(u, v));
2232 
2233  u=1.0;
2234  TEST_ASSERT(bez.f(u, v)==bez2.f(u, v));
2235  TEST_ASSERT(bez.f_v(u, v)==bez2.f_v(u, v));
2236  }
2237 
2238  }
2239 
2241  {
2242  index_type n(6), m(4);
2243  point_type pt[6+1][4+1];
2244 
2245  // create surface with specified control points
2246  pt[0][0] << 0, 0, 15;
2247  pt[1][0] << 2, 6, 15;
2248  pt[2][0] << 3, 0, 15;
2249  pt[3][0] << 5, 4, 15;
2250  pt[4][0] << 7, 1, 15;
2251  pt[5][0] << 5, 5, 15;
2252  pt[6][0] << 10, 6, 15;
2253  pt[0][1] << 0, 0, 11;
2254  pt[1][1] << 2, 6, 11;
2255  pt[2][1] << 3, 0, 11;
2256  pt[3][1] << 5, 4, 11;
2257  pt[4][1] << 7, 1, 11;
2258  pt[5][1] << 5, 5, 11;
2259  pt[6][1] << 10, 6, 11;
2260  pt[0][2] << 0, 0, 3;
2261  pt[1][2] << 2, 6, 3;
2262  pt[2][2] << 3, 0, 3;
2263  pt[3][2] << 5, 4, 3;
2264  pt[4][2] << 7, 1, 3;
2265  pt[5][2] << 5, 5, 3;
2266  pt[6][2] << 10, 6, 3;
2267  pt[0][3] << 0, 0, 0;
2268  pt[1][3] << 2, 6, 0;
2269  pt[2][3] << 3, 0, 0;
2270  pt[3][3] << 5, 4, 0;
2271  pt[4][3] << 7, 1, 0;
2272  pt[5][3] << 5, 5, 0;
2273  pt[6][3] << 10, 6, 0;
2274  pt[0][4] << 0, 0, -5;
2275  pt[1][4] << 2, 6, -5;
2276  pt[2][4] << 3, 0, -5;
2277  pt[3][4] << 5, 4, -5;
2278  pt[4][4] << 7, 1, -5;
2279  pt[5][4] << 5, 5, -5;
2280  pt[6][4] << 10, 6, -5;
2281 
2282  // create surface with specified dimensions and set control points
2283  bezier_type bez(n, m);
2284 
2285  for (index_type i=0; i<=n; ++i)
2286  {
2287  for (index_type j=0; j<=m; ++j)
2288  {
2289  bez.set_control_point(pt[i][j], i, j);
2290  }
2291  }
2292 
2293  {
2294  data_type d;
2295 
2296  d = bez.eqp_distance_bound(bez);
2297 
2298  TEST_ASSERT(d==0);
2299  }
2300  {
2301  bezier_type bez2(bez);
2302 
2303  data_type d;
2304 
2305  bez2.promote_u_to(bez.degree_u()+3);
2306  bez2.promote_v_to(bez.degree_v()+3);
2307 
2308  d = bez.eqp_distance_bound(bez2);
2309 
2310  TEST_ASSERT(d==0);
2311  }
2312  {
2313  bezier_type bez2(bez);
2314 
2315  data_type d;
2316  data_type offset(3);
2317 
2318  point_type trans;
2319 
2320  trans << 0, 0, offset;
2321 
2322  bez2.translate(trans);
2323 
2324  d = bez.eqp_distance_bound(bez2);
2325 
2326  TEST_ASSERT(d==offset);
2327  }
2328  }
2329 
2330  void split_test()
2331  {
2332  // split with known control points
2333  {
2334  index_type n(3), m(3);
2335  point_type pt[3+1][3+1], pt_lo[3+1][3+1], pt_hi[3+1][3+1];
2336 
2337  // create surface with specified control points
2338  pt[0][0] << 0, 0, 0;
2339  pt[1][0] << 0, 2, 0;
2340  pt[2][0] << 8, 2, 0;
2341  pt[3][0] << 4, 0, 0;
2342  pt[0][1] << 0, 0, 1;
2343  pt[1][1] << 0, 2, 1;
2344  pt[2][1] << 8, 2, 1;
2345  pt[3][1] << 4, 0, 1;
2346  pt[0][2] << 0, 0, 3;
2347  pt[1][2] << 0, 2, 3;
2348  pt[2][2] << 8, 2, 3;
2349  pt[3][2] << 4, 0, 3;
2350  pt[0][3] << 0, 0, 4;
2351  pt[1][3] << 0, 2, 4;
2352  pt[2][3] << 8, 2, 4;
2353  pt[3][3] << 4, 0, 4;
2354 
2355  // create surface with specified dimensions and set control points
2356  bezier_type bez(n, m), bez_lo, bez_hi;
2357 
2358  for (index_type i=0; i<=n; ++i)
2359  {
2360  for (index_type j=0; j<=m; ++j)
2361  {
2362  bez.set_control_point(pt[i][j], i, j);
2363  }
2364  }
2365 
2366  // split the surface
2367  bez.split_u(bez_lo, bez_hi, 0.5);
2368 
2369  // set the reference control points
2370  pt_lo[0][0] << 0.0, 0.0, 0;
2371  pt_lo[1][0] << 0.0, 1.0, 0;
2372  pt_lo[2][0] << 2.0, 1.5, 0;
2373  pt_lo[3][0] << 3.5, 1.5, 0;
2374  pt_lo[0][1] << 0.0, 0.0, 1;
2375  pt_lo[1][1] << 0.0, 1.0, 1;
2376  pt_lo[2][1] << 2.0, 1.5, 1;
2377  pt_lo[3][1] << 3.5, 1.5, 1;
2378  pt_lo[0][2] << 0.0, 0.0, 3;
2379  pt_lo[1][2] << 0.0, 1.0, 3;
2380  pt_lo[2][2] << 2.0, 1.5, 3;
2381  pt_lo[3][2] << 3.5, 1.5, 3;
2382  pt_lo[0][3] << 0.0, 0.0, 4;
2383  pt_lo[1][3] << 0.0, 1.0, 4;
2384  pt_lo[2][3] << 2.0, 1.5, 4;
2385  pt_lo[3][3] << 3.5, 1.5, 4;
2386 
2387  pt_hi[0][0] << 3.5, 1.5, 0;
2388  pt_hi[1][0] << 5.0, 1.5, 0;
2389  pt_hi[2][0] << 6.0, 1.0, 0;
2390  pt_hi[3][0] << 4.0, 0.0, 0;
2391  pt_hi[0][1] << 3.5, 1.5, 1;
2392  pt_hi[1][1] << 5.0, 1.5, 1;
2393  pt_hi[2][1] << 6.0, 1.0, 1;
2394  pt_hi[3][1] << 4.0, 0.0, 1;
2395  pt_hi[0][2] << 3.5, 1.5, 3;
2396  pt_hi[1][2] << 5.0, 1.5, 3;
2397  pt_hi[2][2] << 6.0, 1.0, 3;
2398  pt_hi[3][2] << 4.0, 0.0, 3;
2399  pt_hi[0][3] << 3.5, 1.5, 4;
2400  pt_hi[1][3] << 5.0, 1.5, 4;
2401  pt_hi[2][3] << 6.0, 1.0, 4;
2402  pt_hi[3][3] << 4.0, 0.0, 4;
2403 
2404  // compare control points
2405  for (index_type i=0; i<=n; ++i)
2406  {
2407  for (index_type j=0; j<=m; ++j)
2408  {
2409  TEST_ASSERT(bez_lo.get_control_point(i,j)==pt_lo[i][j]);
2410  TEST_ASSERT(bez_hi.get_control_point(i,j)==pt_hi[i][j]);
2411  }
2412  }
2413  }
2414 
2415  // split in u-direction
2416  {
2417  index_type n(3), m(3);
2418  point_type pt[3+1][3+1];
2419  data_type u0(0.25), ul, uh, vl, vh;
2420 
2421  // create surface with specified control points
2422  pt[0][0] << -15, 0, 15;
2423  pt[1][0] << -5, 5, 15;
2424  pt[2][0] << 5, 5, 15;
2425  pt[3][0] << 15, 0, 15;
2426  pt[0][1] << -15, 5, 5;
2427  pt[1][1] << -5, 5, 5;
2428  pt[2][1] << 5, 5, 5;
2429  pt[3][1] << 15, 5, 5;
2430  pt[0][2] << -15, 5, -5;
2431  pt[1][2] << -5, 5, -5;
2432  pt[2][2] << 5, 5, -5;
2433  pt[3][2] << 15, 5, -5;
2434  pt[0][3] << -15, 0, -15;
2435  pt[1][3] << -5, 5, -15;
2436  pt[2][3] << 5, 5, -15;
2437  pt[3][3] << 15, 0, -15;
2438 
2439  // create surface with specified dimensions and set control points
2440  bezier_type bez(n, m), bez_lo, bez_hi;
2441  point_type pt_out, pt_ref;
2442 
2443  for (index_type i=0; i<=n; ++i)
2444  {
2445  for (index_type j=0; j<=m; ++j)
2446  {
2447  bez.set_control_point(pt[i][j], i, j);
2448  }
2449  }
2450 
2451  // split surface
2452  bez.split_u(bez_lo, bez_hi, u0);
2453 
2454  // test lower side matches
2455  ul=0.125;
2456  vl=0.25;
2457  pt_ref=bez.f(ul, vl);
2458  pt_out=bez_lo.f(ul/u0, vl);
2459  TEST_ASSERT(pt_ref==pt_out);
2460  pt_ref=bez.f_u(ul, vl);
2461  pt_out=bez_lo.f_u(ul/u0, vl)/u0;
2462  TEST_ASSERT(pt_ref==pt_out);
2463  pt_ref=bez.f_v(ul, vl);
2464  pt_out=bez_lo.f_v(ul/u0, vl);
2465  TEST_ASSERT(pt_ref==pt_out);
2466  pt_ref=bez.f_uu(ul, vl);
2467  pt_out=bez_lo.f_uu(ul/u0, vl)/u0/u0;
2468  TEST_ASSERT(pt_ref==pt_out);
2469  pt_ref=bez.f_uv(ul, vl);
2470  pt_out=bez_lo.f_uv(ul/u0, vl)/u0;
2471  TEST_ASSERT(pt_ref==pt_out);
2472  pt_ref=bez.f_vv(ul, vl);
2473  pt_out=bez_lo.f_vv(ul/u0, vl);
2474  TEST_ASSERT(pt_ref==pt_out);
2475  pt_ref=bez.f_uuu(ul, vl);
2476  pt_out=bez_lo.f_uuu(ul/u0, vl)/u0/u0/u0;
2477  TEST_ASSERT(pt_ref==pt_out);
2478  pt_ref=bez.f_uuv(ul, vl);
2479  pt_out=bez_lo.f_uuv(ul/u0, vl)/u0/u0;
2480  TEST_ASSERT(pt_ref==pt_out);
2481  pt_ref=bez.f_uvv(ul, vl);
2482  pt_out=bez_lo.f_uvv(ul/u0, vl)/u0;
2483  TEST_ASSERT(pt_ref==pt_out);
2484  pt_ref=bez.f_vvv(ul, vl);
2485  pt_out=bez_lo.f_vvv(ul/u0, vl);
2486  TEST_ASSERT(pt_ref==pt_out);
2487 
2488  // test upper side matches
2489  uh=0.75;
2490  vh=0.75;
2491  pt_ref=bez.f(uh, vh);
2492  pt_out=bez_hi.f((uh-u0)/(1-u0), vh);
2493  TEST_ASSERT((pt_ref - pt_out).norm()<5*std::numeric_limits<data_type>::epsilon());
2494  pt_ref=bez.f_u(uh, vh);
2495  pt_out=bez_hi.f_u((uh-u0)/(1-u0), vh)/(1-u0);
2496  TEST_ASSERT((pt_ref - pt_out).norm()<2*std::numeric_limits<data_type>::epsilon());
2497  pt_ref=bez.f_v(uh, vh);
2498  pt_out=bez_hi.f_v((uh-u0)/(1-u0), vh);
2499  TEST_ASSERT((pt_ref - pt_out).norm()<3*std::numeric_limits<data_type>::epsilon());
2500  pt_ref=bez.f_uu(uh, vh);
2501  pt_out=bez_hi.f_uu((uh-u0)/(1-u0), vh)/(1-u0)/(1-u0);
2502  TEST_ASSERT((pt_ref - pt_out).norm()<2*std::numeric_limits<data_type>::epsilon());
2503  pt_ref=bez.f_uv(uh, vh);
2504  pt_out=bez_hi.f_uv((uh-u0)/(1-u0), vh)/(1-u0);
2505  TEST_ASSERT((pt_ref - pt_out).norm()<2*std::numeric_limits<data_type>::epsilon());
2506  pt_ref=bez.f_vv(uh, vh);
2507  pt_out=bez_hi.f_vv((uh-u0)/(1-u0), vh);
2508  TEST_ASSERT((pt_ref - pt_out).norm()<2*std::numeric_limits<data_type>::epsilon());
2509  pt_ref=bez.f_uuu(uh, vh);
2510  pt_out=bez_hi.f_uuu((uh-u0)/(1-u0), vh)/(1-u0)/(1-u0)/(1-u0);
2511  TEST_ASSERT((pt_ref - pt_out).norm()<2*std::numeric_limits<data_type>::epsilon());
2512  pt_ref=bez.f_uuv(uh, vh);
2513  pt_out=bez_hi.f_uuv((uh-u0)/(1-u0), vh)/(1-u0)/(1-u0);
2514  TEST_ASSERT((pt_ref - pt_out).norm()<2*std::numeric_limits<data_type>::epsilon());
2515  pt_ref=bez.f_uvv(uh, vh);
2516  pt_out=bez_hi.f_uvv((uh-u0)/(1-u0), vh)/(1-u0);
2517  TEST_ASSERT((pt_ref - pt_out).norm()<2*std::numeric_limits<data_type>::epsilon());
2518  pt_ref=bez.f_vvv(uh, vh);
2519  pt_out=bez_hi.f_vvv((uh-u0)/(1-u0), vh);
2520  TEST_ASSERT((pt_ref - pt_out).norm()<2*std::numeric_limits<data_type>::epsilon());
2521  }
2522 
2523  // split in v-direction
2524  {
2525  index_type n(3), m(3);
2526  point_type pt[3+1][3+1];
2527  data_type v0(0.25), ul, uh, vl, vh;
2528 
2529  // create surface with specified control points
2530  pt[0][0] << -15, 0, 15;
2531  pt[1][0] << -5, 5, 15;
2532  pt[2][0] << 5, 5, 15;
2533  pt[3][0] << 15, 0, 15;
2534  pt[0][1] << -15, 5, 5;
2535  pt[1][1] << -5, 5, 5;
2536  pt[2][1] << 5, 5, 5;
2537  pt[3][1] << 15, 5, 5;
2538  pt[0][2] << -15, 5, -5;
2539  pt[1][2] << -5, 5, -5;
2540  pt[2][2] << 5, 5, -5;
2541  pt[3][2] << 15, 5, -5;
2542  pt[0][3] << -15, 0, -15;
2543  pt[1][3] << -5, 5, -15;
2544  pt[2][3] << 5, 5, -15;
2545  pt[3][3] << 15, 0, -15;
2546 
2547  // create surface with specified dimensions and set control points
2548  bezier_type bez(n, m), bez_lo, bez_hi;
2549  point_type pt_out, pt_ref;
2550 
2551  for (index_type i=0; i<=n; ++i)
2552  {
2553  for (index_type j=0; j<=m; ++j)
2554  {
2555  bez.set_control_point(pt[i][j], i, j);
2556  }
2557  }
2558 
2559  // split surface
2560  bez.split_v(bez_lo, bez_hi, v0);
2561 
2562  // test lower side matches
2563  ul=0.25;
2564  vl=0.125;
2565  pt_ref=bez.f(ul, vl);
2566  pt_out=bez_lo.f(ul, vl/v0);
2567  TEST_ASSERT(pt_ref==pt_out);
2568  pt_ref=bez.f_u(ul, vl);
2569  pt_out=bez_lo.f_u(ul, vl/v0);
2570  TEST_ASSERT(pt_ref==pt_out);
2571  pt_ref=bez.f_v(ul, vl);
2572  pt_out=bez_lo.f_v(ul, vl/v0)/v0;
2573  TEST_ASSERT(pt_ref==pt_out);
2574  pt_ref=bez.f_uu(ul, vl);
2575  pt_out=bez_lo.f_uu(ul, vl/v0);
2576  TEST_ASSERT(pt_ref==pt_out);
2577  pt_ref=bez.f_uv(ul, vl);
2578  pt_out=bez_lo.f_uv(ul, vl/v0)/v0;
2579  TEST_ASSERT(pt_ref==pt_out);
2580  pt_ref=bez.f_vv(ul, vl);
2581  pt_out=bez_lo.f_vv(ul, vl/v0)/v0/v0;
2582  TEST_ASSERT(pt_ref==pt_out);
2583  pt_ref=bez.f_uuu(ul, vl);
2584  pt_out=bez_lo.f_uuu(ul, vl/v0);
2585  TEST_ASSERT(pt_ref==pt_out);
2586  pt_ref=bez.f_uuv(ul, vl);
2587  pt_out=bez_lo.f_uuv(ul, vl/v0)/v0;
2588  TEST_ASSERT(pt_ref==pt_out);
2589  pt_ref=bez.f_uvv(ul, vl);
2590  pt_out=bez_lo.f_uvv(ul, vl/v0)/v0/v0;
2591  TEST_ASSERT(pt_ref==pt_out);
2592  pt_ref=bez.f_vvv(ul, vl);
2593  pt_out=bez_lo.f_vvv(ul, vl/v0)/v0/v0/v0;
2594  TEST_ASSERT(pt_ref==pt_out);
2595 
2596  // test upper side matches
2597  uh=0.75;
2598  vh=0.75;
2599  pt_ref=bez.f(uh, vh);
2600  pt_out=bez_hi.f(uh, (vh-v0)/(1-v0));
2601  TEST_ASSERT((pt_ref - pt_out).norm()<5*std::numeric_limits<data_type>::epsilon());
2602  pt_ref=bez.f_u(uh, vh);
2603  pt_out=bez_hi.f_u(uh, (vh-v0)/(1-v0));
2604  TEST_ASSERT((pt_ref - pt_out).norm()<3*std::numeric_limits<data_type>::epsilon());
2605  pt_ref=bez.f_v(uh, vh);
2606  pt_out=bez_hi.f_v(uh, (vh-v0)/(1-v0))/(1-v0);
2607  TEST_ASSERT((pt_ref - pt_out).norm()<3*std::numeric_limits<data_type>::epsilon());
2608  pt_ref=bez.f_uu(uh, vh);
2609  pt_out=bez_hi.f_uu(uh, (vh-v0)/(1-v0));
2610  TEST_ASSERT((pt_ref - pt_out).norm()<2*std::numeric_limits<data_type>::epsilon());
2611  pt_ref=bez.f_uv(uh, vh);
2612  pt_out=bez_hi.f_uv(uh, (vh-v0)/(1-v0))/(1-v0);
2613  TEST_ASSERT((pt_ref - pt_out).norm()<2*std::numeric_limits<data_type>::epsilon());
2614  pt_ref=bez.f_vv(uh, vh);
2615  pt_out=bez_hi.f_vv(uh, (vh-v0)/(1-v0))/(1-v0)/(1-v0);
2616  TEST_ASSERT((pt_ref - pt_out).norm()<2*std::numeric_limits<data_type>::epsilon());
2617  pt_ref=bez.f_uuu(uh, vh);
2618  pt_out=bez_hi.f_uuu(uh, (vh-v0)/(1-v0));
2619  TEST_ASSERT((pt_ref - pt_out).norm()<2*std::numeric_limits<data_type>::epsilon());
2620  pt_ref=bez.f_uuv(uh, vh);
2621  pt_out=bez_hi.f_uuv(uh, (vh-v0)/(1-v0))/(1-v0);
2622  TEST_ASSERT((pt_ref - pt_out).norm()<2*std::numeric_limits<data_type>::epsilon());
2623  pt_ref=bez.f_uvv(uh, vh);
2624  pt_out=bez_hi.f_uvv(uh, (vh-v0)/(1-v0))/(1-v0)/(1-v0);
2625  TEST_ASSERT((pt_ref - pt_out).norm()<2*std::numeric_limits<data_type>::epsilon());
2626  pt_ref=bez.f_vvv(uh, vh);
2627  pt_out=bez_hi.f_vvv(uh, (vh-v0)/(1-v0))/(1-v0)/(1-v0)/(1-v0);
2628  TEST_ASSERT((pt_ref - pt_out).norm()<2*std::numeric_limits<data_type>::epsilon());
2629  }
2630  }
2631 
2633  {
2634  // simple surface test
2635  {
2636  index_type n(3), m(3);
2637  point_type pt[3+1][3+1], pt_out, pt_ref;
2638  data_type u, v;
2639 
2640  // create surface with specified control points
2641  pt[0][0] << -15, 0, 15;
2642  pt[1][0] << -5, 5, 15;
2643  pt[2][0] << 5, 5, 15;
2644  pt[3][0] << 15, 0, 15;
2645  pt[0][1] << -15, 5, 5;
2646  pt[1][1] << -5, 5, 5;
2647  pt[2][1] << 5, 5, 5;
2648  pt[3][1] << 15, 5, 5;
2649  pt[0][2] << -15, 5, -5;
2650  pt[1][2] << -5, 5, -5;
2651  pt[2][2] << 5, 5, -5;
2652  pt[3][2] << 15, 5, -5;
2653  pt[0][3] << -15, 0, -15;
2654  pt[1][3] << -5, 5, -15;
2655  pt[2][3] << 5, 5, -15;
2656  pt[3][3] << 15, 0, -15;
2657 
2658  // create surface with specified dimensions and set control points
2659  bezier_type bez(n, m);
2660 
2661  for (index_type i=0; i<=n; ++i)
2662  {
2663  for (index_type j=0; j<=m; ++j)
2664  {
2665  bez.set_control_point(pt[i][j], i, j);
2666  }
2667  }
2668 
2669  // get the normal and the tangent vectors
2670  point_type normal, S_u, S_v, tmp;
2671 
2672  u=0.5;
2673  v=0.5;
2674  normal=bez.normal(u, v);
2675  S_u=bez.f_u(u, v);
2676  S_v=bez.f_v(u, v);
2677 
2678  // test the normal vector
2679  TEST_ASSERT(tol.approximately_equal(normal.norm(), 1));
2680  tmp=S_u.cross(S_v);
2681  tmp.normalize();
2682  TEST_ASSERT(tol.approximately_equal(tmp, normal));
2683  }
2684 
2685  // flat degenerate surface test
2686  {
2687  bezier_type bez(3, 3);
2688  point_type cp[4], x, y, origin;
2689  data_type k, xr, yr;
2690  index_type i;
2691 
2693  x << 1, 0, 0;
2694  y << 0, 1, 0;
2695  origin << 0, 0, 0;
2696 
2697  // set the first section
2698  xr=0;
2699  yr=0;
2700  cp[0]=xr*x+origin;
2701  cp[1]=xr*x+yr*k*y+origin;
2702  cp[2]=xr*k*x+yr*y+origin;
2703  cp[3]=yr*y+origin;
2704  for (i=0; i<4; ++i)
2705  {
2706  bez.set_control_point(cp[i], i, 0);
2707  }
2708 
2709  // set the second section
2710  xr=0.5;
2711  yr=0.5;
2712  cp[0]=xr*x+origin;
2713  cp[1]=xr*x+yr*k*y+origin;
2714  cp[2]=xr*k*x+yr*y+origin;
2715  cp[3]=yr*y+origin;
2716  for (i=0; i<4; ++i)
2717  {
2718  bez.set_control_point(cp[i], i, 1);
2719  }
2720 
2721  // set the third section
2722  xr=1;
2723  yr=1;
2724  cp[0]=xr*x+origin;
2725  cp[1]=xr*x+yr*k*y+origin;
2726  cp[2]=xr*k*x+yr*y+origin;
2727  cp[3]=yr*y+origin;
2728  for (i=0; i<4; ++i)
2729  {
2730  bez.set_control_point(cp[i], i, 2);
2731  }
2732 
2733  // set the last section
2734  xr=1.5;
2735  yr=1.5;
2736  cp[0]=xr*x+origin;
2737  cp[1]=xr*x+yr*k*y+origin;
2738  cp[2]=xr*k*x+yr*y+origin;
2739  cp[3]=yr*y+origin;
2740  for (i=0; i<4; ++i)
2741  {
2742  bez.set_control_point(cp[i], i, 3);
2743  }
2744 
2745  // get the normal and the tangent vectors
2746  point_type normal, ref_normal, S_u, S_v;
2747  data_type u, v;
2748 
2749  u=0.5;
2750  v=0.5;
2751  normal=bez.normal(u, v);
2752  S_u=bez.f_u(u, v);
2753  S_v=bez.f_v(u, v);
2754 
2755  // test the normal vector
2756  TEST_ASSERT(tol.approximately_equal(normal.norm(), 1));
2757  ref_normal=S_u.cross(S_v);
2758  ref_normal.normalize();
2759  TEST_ASSERT(tol.approximately_equal(ref_normal, normal));
2760 
2761  u=0.5;
2762  v=0;
2763  normal=bez.normal(u, v);
2764  ref_normal=bez.normal(u, v+std::numeric_limits<data_type>::epsilon());
2765  TEST_ASSERT(tol.approximately_equal(ref_normal, normal));
2766  }
2767 
2768  // 3D degenerate surface test
2769  {
2770  bezier_type bez(3, 3);
2771  point_type cp[4], x, y, origin;
2772  data_type k, xr, yr;
2773  index_type i;
2774 
2776  x << 1, 0, 0;
2777  y << 0, 1, 0;
2778 
2779  // set the first section
2780  xr=0;
2781  yr=0;
2782  origin << 0, 0, 0;
2783  cp[0]=xr*x+origin;
2784  cp[1]=xr*x+yr*k*y+origin;
2785  cp[2]=xr*k*x+yr*y+origin;
2786  cp[3]=yr*y+origin;
2787  for (i=0; i<4; ++i)
2788  {
2789  bez.set_control_point(cp[i], i, 0);
2790  }
2791 
2792  // set the second section
2793  xr=0.5;
2794  yr=0.5;
2795  origin << 0, 0, 1;
2796  cp[0]=xr*x+origin;
2797  cp[1]=xr*x+yr*k*y+origin;
2798  cp[2]=xr*k*x+yr*y+origin;
2799  cp[3]=yr*y+origin;
2800  for (i=0; i<4; ++i)
2801  {
2802  bez.set_control_point(cp[i], i, 1);
2803  }
2804 
2805  // set the third section
2806  xr=1;
2807  yr=1;
2808  origin << 0, 0, 2;
2809  cp[0]=xr*x+origin;
2810  cp[1]=xr*x+yr*k*y+origin;
2811  cp[2]=xr*k*x+yr*y+origin;
2812  cp[3]=yr*y+origin;
2813  for (i=0; i<4; ++i)
2814  {
2815  bez.set_control_point(cp[i], i, 2);
2816  }
2817 
2818  // set the last section
2819  xr=1.5;
2820  yr=1.5;
2821  origin << 0, 0, 3;
2822  cp[0]=xr*x+origin;
2823  cp[1]=xr*x+yr*k*y+origin;
2824  cp[2]=xr*k*x+yr*y+origin;
2825  cp[3]=yr*y+origin;
2826  for (i=0; i<4; ++i)
2827  {
2828  bez.set_control_point(cp[i], i, 3);
2829  }
2830 
2831  // get the normal and the tangent vectors
2832  point_type normal, ref_normal, S_u, S_v;
2833  data_type u, v;
2834 
2835  u=0.5;
2836  v=0.5;
2837  normal=bez.normal(u, v);
2838  S_u=bez.f_u(u, v);
2839  S_v=bez.f_v(u, v);
2840 
2841  // test the normal vector
2842  TEST_ASSERT(tol.approximately_equal(normal.norm(), 1));
2843  ref_normal=S_u.cross(S_v);
2844  ref_normal.normalize();
2845  TEST_ASSERT(tol.approximately_equal(ref_normal, normal));
2846 
2847  u=0.5;
2848  v=0;
2849  normal=bez.normal(u, v);
2850  ref_normal=bez.normal(u, v+std::numeric_limits<data_type>::epsilon());
2851  TEST_ASSERT(tol.approximately_equal(ref_normal, normal));
2852  }
2853  }
2854 
2856  {
2857  index_type n(3), m(3);
2858  point_type pt[3+1][3+1], pt_out, pt_ref;
2859  data_type u, v, d;
2860  curve_type c;
2861 
2862  // create surface with specified control points
2863  pt[0][0] << -15, 0, 15;
2864  pt[1][0] << -5, 5, 15;
2865  pt[2][0] << 5, 5, 15;
2866  pt[3][0] << 15, 0, 15;
2867  pt[0][1] << -15, 5, 5;
2868  pt[1][1] << -5, 5, 5;
2869  pt[2][1] << 5, 5, 5;
2870  pt[3][1] << 15, 5, 5;
2871  pt[0][2] << -15, 5, -5;
2872  pt[1][2] << -5, 5, -5;
2873  pt[2][2] << 5, 5, -5;
2874  pt[3][2] << 15, 5, -5;
2875  pt[0][3] << -15, 0, -15;
2876  pt[1][3] << -5, 5, -15;
2877  pt[2][3] << 5, 5, -15;
2878  pt[3][3] << 15, 0, -15;
2879 
2880  // create surface with specified dimensions and set control points
2881  bezier_type bez(n, m);
2882 
2883  for (index_type i=0; i<=n; ++i)
2884  {
2885  for (index_type j=0; j<=m; ++j)
2886  {
2887  bez.set_control_point(pt[i][j], i, j);
2888  }
2889  }
2890 
2891  // Extract u curve
2892  u=0;
2893  bez.get_uconst_curve(c, u);
2894 
2895  v=0;
2896  pt_ref=bez.f(u, v);
2897  pt_out=c.f(v);
2898  TEST_ASSERT(pt_out==pt_ref);
2899 
2900  v=static_cast<data_type>(0.1);
2901  pt_ref=bez.f(u, v);
2902  pt_out=c.f(v);
2903  TEST_ASSERT(pt_out==pt_ref);
2904 
2905  v=static_cast<data_type>(0.2);
2906  pt_ref=bez.f(u, v);
2907  pt_out=c.f(v);
2908  TEST_ASSERT(pt_out==pt_ref);
2909 
2910  v=static_cast<data_type>(0.4);
2911  pt_ref=bez.f(u, v);
2912  pt_out=c.f(v);
2913  TEST_ASSERT(pt_out==pt_ref);
2914 
2915  v=static_cast<data_type>(0.6);
2916  pt_ref=bez.f(u, v);
2917  pt_out=c.f(v);
2918  TEST_ASSERT(pt_out==pt_ref);
2919 
2920  v=static_cast<data_type>(0.8);
2921  pt_ref=bez.f(u, v);
2922  pt_out=c.f(v);
2923  TEST_ASSERT(pt_out==pt_ref);
2924 
2925  v=1;
2926  pt_ref=bez.f(u, v);
2927  pt_out=c.f(v);
2928  TEST_ASSERT(pt_out==pt_ref);
2929 
2930  // Extract u curve
2931  u=static_cast<data_type>(0.1);
2932  bez.get_uconst_curve(c, u);
2933 
2934  v=0;
2935  pt_ref=bez.f(u, v);
2936  pt_out=c.f(v);
2937  TEST_ASSERT(pt_out==pt_ref);
2938 
2939  v=static_cast<data_type>(0.1);
2940  pt_ref=bez.f(u, v);
2941  pt_out=c.f(v);
2942  TEST_ASSERT(pt_out==pt_ref);
2943 
2944  v=static_cast<data_type>(0.2);
2945  pt_ref=bez.f(u, v);
2946  pt_out=c.f(v);
2947  TEST_ASSERT(pt_out==pt_ref);
2948 
2949  v=static_cast<data_type>(0.4);
2950  pt_ref=bez.f(u, v);
2951  pt_out=c.f(v);
2952  TEST_ASSERT(pt_out==pt_ref);
2953 
2954  v=static_cast<data_type>(0.6);
2955  pt_ref=bez.f(u, v);
2956  pt_out=c.f(v);
2957  TEST_ASSERT(pt_out==pt_ref);
2958 
2959  v=static_cast<data_type>(0.8);
2960  pt_ref=bez.f(u, v);
2961  pt_out=c.f(v);
2962  TEST_ASSERT(pt_out==pt_ref);
2963 
2964  v=1;
2965  pt_ref=bez.f(u, v);
2966  pt_out=c.f(v);
2967  TEST_ASSERT(pt_out==pt_ref);
2968 
2969  // Extract u curve
2970  u=static_cast<data_type>(0.2);
2971  bez.get_uconst_curve(c, u);
2972 
2973  v=0;
2974  pt_ref=bez.f(u, v);
2975  pt_out=c.f(v);
2976  TEST_ASSERT(pt_out==pt_ref);
2977 
2978  v=static_cast<data_type>(0.1);
2979  pt_ref=bez.f(u, v);
2980  pt_out=c.f(v);
2981  TEST_ASSERT(pt_out==pt_ref);
2982 
2983  v=static_cast<data_type>(0.2);
2984  pt_ref=bez.f(u, v);
2985  pt_out=c.f(v);
2986  TEST_ASSERT(pt_out==pt_ref);
2987 
2988  v=static_cast<data_type>(0.4);
2989  pt_ref=bez.f(u, v);
2990  pt_out=c.f(v);
2991  TEST_ASSERT(pt_out==pt_ref);
2992 
2993  v=static_cast<data_type>(0.6);
2994  pt_ref=bez.f(u, v);
2995  pt_out=c.f(v);
2996  TEST_ASSERT(pt_out==pt_ref);
2997 
2998  v=static_cast<data_type>(0.8);
2999  pt_ref=bez.f(u, v);
3000  pt_out=c.f(v);
3001  TEST_ASSERT(pt_out==pt_ref);
3002 
3003  v=1;
3004  pt_ref=bez.f(u, v);
3005  pt_out=c.f(v);
3006  TEST_ASSERT(pt_out==pt_ref);
3007 
3008  // Extract u curve
3009  u=0.5;
3010  bez.get_uconst_curve(c, u);
3011 
3012  v=0;
3013  pt_ref=bez.f(u, v);
3014  pt_out=c.f(v);
3015  TEST_ASSERT(pt_out==pt_ref);
3016 
3017  v=static_cast<data_type>(0.1);
3018  pt_ref=bez.f(u, v);
3019  pt_out=c.f(v);
3020  TEST_ASSERT(pt_out==pt_ref);
3021 
3022  v=static_cast<data_type>(0.2);
3023  pt_ref=bez.f(u, v);
3024  pt_out=c.f(v);
3025  TEST_ASSERT(pt_out==pt_ref);
3026 
3027  v=static_cast<data_type>(0.4);
3028  pt_ref=bez.f(u, v);
3029  pt_out=c.f(v);
3030  TEST_ASSERT(pt_out==pt_ref);
3031 
3032  v=static_cast<data_type>(0.6);
3033  pt_ref=bez.f(u, v);
3034  pt_out=c.f(v);
3035  TEST_ASSERT(pt_out==pt_ref);
3036 
3037  v=static_cast<data_type>(0.8);
3038  pt_ref=bez.f(u, v);
3039  pt_out=c.f(v);
3040  TEST_ASSERT(pt_out==pt_ref);
3041 
3042  v=1;
3043  pt_ref=bez.f(u, v);
3044  pt_out=c.f(v);
3045  TEST_ASSERT(pt_out==pt_ref);
3046 
3047  // Extract u curve
3048  u=1;
3049  bez.get_uconst_curve(c, u);
3050 
3051  v=0;
3052  pt_ref=bez.f(u, v);
3053  pt_out=c.f(v);
3054  TEST_ASSERT(pt_out==pt_ref);
3055 
3056  v=static_cast<data_type>(0.1);
3057  pt_ref=bez.f(u, v);
3058  pt_out=c.f(v);
3059  TEST_ASSERT(pt_out==pt_ref);
3060 
3061  v=static_cast<data_type>(0.2);
3062  pt_ref=bez.f(u, v);
3063  pt_out=c.f(v);
3064  TEST_ASSERT(pt_out==pt_ref);
3065 
3066  v=static_cast<data_type>(0.4);
3067  pt_ref=bez.f(u, v);
3068  pt_out=c.f(v);
3069  TEST_ASSERT(pt_out==pt_ref);
3070 
3071  v=static_cast<data_type>(0.6);
3072  pt_ref=bez.f(u, v);
3073  pt_out=c.f(v);
3074  TEST_ASSERT(pt_out==pt_ref);
3075 
3076  v=static_cast<data_type>(0.8);
3077  pt_ref=bez.f(u, v);
3078  pt_out=c.f(v);
3079  TEST_ASSERT(pt_out==pt_ref);
3080 
3081  v=1;
3082  pt_ref=bez.f(u, v);
3083  pt_out=c.f(v);
3084  TEST_ASSERT(pt_out==pt_ref);
3085 
3086 
3087 
3088  // Extract v curve
3089  v=0;
3090  bez.get_vconst_curve(c, v);
3091 
3092  u=0;
3093  pt_ref=bez.f(u, v);
3094  pt_out=c.f(u);
3095  TEST_ASSERT(pt_out==pt_ref);
3096 
3097  u=static_cast<data_type>(0.1);
3098  pt_ref=bez.f(u, v);
3099  pt_out=c.f(u);
3100  TEST_ASSERT(pt_out==pt_ref);
3101 
3102  u=static_cast<data_type>(0.2);
3103  pt_ref=bez.f(u, v);
3104  pt_out=c.f(u);
3105  TEST_ASSERT(pt_out==pt_ref);
3106 
3107  u=static_cast<data_type>(0.4);
3108  pt_ref=bez.f(u, v);
3109  pt_out=c.f(u);
3110  TEST_ASSERT(pt_out==pt_ref);
3111 
3112  u=static_cast<data_type>(0.6);
3113  pt_ref=bez.f(u, v);
3114  pt_out=c.f(u);
3115  TEST_ASSERT(pt_out==pt_ref);
3116 
3117  u=static_cast<data_type>(0.8);
3118  pt_ref=bez.f(u, v);
3119  pt_out=c.f(u);
3120  TEST_ASSERT(pt_out==pt_ref);
3121 
3122  u=1;
3123  pt_ref=bez.f(u, v);
3124  pt_out=c.f(u);
3125  TEST_ASSERT(pt_out==pt_ref);
3126 
3127  // Extract v curve
3128  v=static_cast<data_type>(0.1);
3129  bez.get_vconst_curve(c, v);
3130 
3131  u=0;
3132  pt_ref=bez.f(u, v);
3133  pt_out=c.f(u);
3134  TEST_ASSERT(pt_out==pt_ref);
3135 
3136  u=static_cast<data_type>(0.1);
3137  pt_ref=bez.f(u, v);
3138  pt_out=c.f(u);
3139  d=(pt_ref - pt_out).norm();
3140  TEST_ASSERT(d<std::numeric_limits<data_type>::epsilon()*30);
3141 
3142  u=static_cast<data_type>(0.2);
3143  pt_ref=bez.f(u, v);
3144  pt_out=c.f(u);
3145  if (typeid(data_type)==typeid(float))
3146  {
3147  d=(pt_ref - pt_out).norm();
3148  TEST_ASSERT(d<std::numeric_limits<data_type>::epsilon()*30);
3149  }
3150  else
3151  {
3152  TEST_ASSERT(pt_out==pt_ref);
3153  }
3154 
3155  u=static_cast<data_type>(0.4);
3156  pt_ref=bez.f(u, v);
3157  pt_out=c.f(u);
3158  d=(pt_ref - pt_out).norm();
3159  TEST_ASSERT(d<std::numeric_limits<data_type>::epsilon()*30);
3160 
3161  u=static_cast<data_type>(0.6);
3162  pt_ref=bez.f(u, v);
3163  pt_out=c.f(u);
3164  d=(pt_ref - pt_out).norm();
3165  TEST_ASSERT(d<std::numeric_limits<data_type>::epsilon()*30);
3166 
3167  u=static_cast<data_type>(0.8);
3168  pt_ref=bez.f(u, v);
3169  pt_out=c.f(u);
3170  if (typeid(data_type)==typeid(float))
3171  {
3172  d=(pt_ref - pt_out).norm();
3173  TEST_ASSERT(d<std::numeric_limits<data_type>::epsilon()*30);
3174  }
3175  else
3176  {
3177  TEST_ASSERT(pt_out==pt_ref);
3178  }
3179 
3180  u=1;
3181  pt_ref=bez.f(u, v);
3182  pt_out=c.f(u);
3183  TEST_ASSERT(pt_out==pt_ref);
3184 
3185  // Extract v curve
3186  v=static_cast<data_type>(0.2);
3187  bez.get_vconst_curve(c, v);
3188 
3189  u=0;
3190  pt_ref=bez.f(u, v);
3191  pt_out=c.f(u);
3192  TEST_ASSERT(pt_out==pt_ref);
3193 
3194  u=static_cast<data_type>(0.1);
3195  pt_ref=bez.f(u, v);
3196  pt_out=c.f(u);
3197  if (typeid(data_type)==typeid(float))
3198  {
3199  d=(pt_ref - pt_out).norm();
3200  TEST_ASSERT(d<std::numeric_limits<data_type>::epsilon()*30);
3201  }
3202  else
3203  {
3204  TEST_ASSERT(pt_out==pt_ref);
3205  }
3206 
3207  u=static_cast<data_type>(0.2);
3208  pt_ref=bez.f(u, v);
3209  pt_out=c.f(u);
3210  TEST_ASSERT(pt_out==pt_ref);
3211 
3212  u=static_cast<data_type>(0.4);
3213  pt_ref=bez.f(u, v);
3214  pt_out=c.f(u);
3215  d=(pt_ref - pt_out).norm();
3216  TEST_ASSERT(d<std::numeric_limits<data_type>::epsilon()*30);
3217 
3218  u=static_cast<data_type>(0.6);
3219  pt_ref=bez.f(u, v);
3220  pt_out=c.f(u);
3221  d=(pt_ref - pt_out).norm();
3222  TEST_ASSERT(d<std::numeric_limits<data_type>::epsilon()*30);
3223 
3224  u=static_cast<data_type>(0.8);
3225  pt_ref=bez.f(u, v);
3226  pt_out=c.f(u);
3227  d=(pt_ref - pt_out).norm();
3228  TEST_ASSERT(d<std::numeric_limits<data_type>::epsilon()*30);
3229 
3230  u=1;
3231  pt_ref=bez.f(u, v);
3232  pt_out=c.f(u);
3233  TEST_ASSERT(pt_out==pt_ref);
3234 
3235  // Extract v curve
3236  v=static_cast<data_type>(0.4);
3237  bez.get_vconst_curve(c, v);
3238 
3239  u=0;
3240  pt_ref=bez.f(u, v);
3241  pt_out=c.f(u);
3242  TEST_ASSERT(pt_out==pt_ref);
3243 
3244  u=static_cast<data_type>(0.1);
3245  pt_ref=bez.f(u, v);
3246  pt_out=c.f(u);
3247  d=(pt_ref - pt_out).norm();
3248  TEST_ASSERT(d<std::numeric_limits<data_type>::epsilon()*30);
3249 
3250  u=static_cast<data_type>(0.2);
3251  pt_ref=bez.f(u, v);
3252  pt_out=c.f(u);
3253  d=(pt_ref - pt_out).norm();
3254  TEST_ASSERT(d<std::numeric_limits<data_type>::epsilon()*30);
3255 
3256  u=static_cast<data_type>(0.4);
3257  pt_ref=bez.f(u, v);
3258  pt_out=c.f(u);
3259  if (typeid(data_type)==typeid(float))
3260  {
3261  d=(pt_ref - pt_out).norm();
3262  TEST_ASSERT(d<std::numeric_limits<data_type>::epsilon()*30);
3263  }
3264  else
3265  {
3266  TEST_ASSERT(pt_out==pt_ref);
3267  }
3268 
3269  u=static_cast<data_type>(0.6);
3270  pt_ref=bez.f(u, v);
3271  pt_out=c.f(u);
3272  if (typeid(data_type)==typeid(float))
3273  {
3274  d=(pt_ref - pt_out).norm();
3275  TEST_ASSERT(d<std::numeric_limits<data_type>::epsilon()*30);
3276  }
3277  else
3278  {
3279  TEST_ASSERT(pt_out==pt_ref);
3280  }
3281 
3282  u=static_cast<data_type>(0.8);
3283  pt_ref=bez.f(u, v);
3284  pt_out=c.f(u);
3285  if (typeid(data_type)==typeid(float))
3286  {
3287  d=(pt_ref - pt_out).norm();
3288  TEST_ASSERT(d<std::numeric_limits<data_type>::epsilon()*30);
3289  }
3290  else
3291  {
3292  TEST_ASSERT(pt_out==pt_ref);
3293  }
3294 
3295  u=1;
3296  pt_ref=bez.f(u, v);
3297  pt_out=c.f(u);
3298  if (typeid(data_type)==typeid(float))
3299  {
3300  d=(pt_ref - pt_out).norm();
3301  TEST_ASSERT(d<std::numeric_limits<data_type>::epsilon()*30);
3302  }
3303  else
3304  {
3305  TEST_ASSERT(pt_out==pt_ref);
3306  }
3307 
3308  // Extract v curve
3309  v=1;
3310  bez.get_vconst_curve(c, v);
3311 
3312  u=0;
3313  pt_ref=bez.f(u, v);
3314  pt_out=c.f(u);
3315  TEST_ASSERT(pt_out==pt_ref);
3316 
3317  u=static_cast<data_type>(0.1);
3318  pt_ref=bez.f(u, v);
3319  pt_out=c.f(u);
3320  TEST_ASSERT(pt_out==pt_ref);
3321 
3322  u=static_cast<data_type>(0.2);
3323  pt_ref=bez.f(u, v);
3324  pt_out=c.f(u);
3325  TEST_ASSERT(pt_out==pt_ref);
3326 
3327  u=static_cast<data_type>(0.4);
3328  pt_ref=bez.f(u, v);
3329  pt_out=c.f(u);
3330  TEST_ASSERT(pt_out==pt_ref);
3331 
3332  u=static_cast<data_type>(0.6);
3333  pt_ref=bez.f(u, v);
3334  pt_out=c.f(u);
3335  TEST_ASSERT(pt_out==pt_ref);
3336 
3337  u=static_cast<data_type>(0.8);
3338  pt_ref=bez.f(u, v);
3339  pt_out=c.f(u);
3340  TEST_ASSERT(pt_out==pt_ref);
3341 
3342  u=1;
3343  pt_ref=bez.f(u, v);
3344  pt_out=c.f(u);
3345  TEST_ASSERT(pt_out==pt_ref);
3346  }
3347 };
3348 
3349 #endif
3350 
void curvature_test()
Definition: bezier_surface_test_suite.hpp:1009
void set_control_point(const point_type &cp, const index_type &i, const index_type &j)
Definition: bezier.hpp:265
point_type f_uuv(const data_type &u, const data_type &v) const
Definition: bezier.hpp:735
void derivative_3_test()
Definition: bezier_surface_test_suite.hpp:905
point_type f_uuu(const data_type &u, const data_type &v) const
Definition: bezier.hpp:682
void create_circle(std::vector< point_type > &)
Definition: bezier_surface_test_suite.hpp:236
Definition: continuity.hpp:26
point_type f_vvv(const data_type &u, const data_type &v) const
Definition: bezier.hpp:845
void AddTests(const long double &)
Definition: bezier_surface_test_suite.hpp:88
Definition: continuity.hpp:28
Derived1__::Scalar distance(const Eigen::MatrixBase< Derived1__ > &p1, const Eigen::MatrixBase< Derived2__ > &p2)
Definition: distance.hpp:33
void derivative_2_test()
Definition: bezier_surface_test_suite.hpp:815
Definition: bounding_box.hpp:27
void promote_v_to(index_type target_degree)
Definition: bezier.hpp:1027
point_type f_uu(const data_type &u, const data_type &v) const
Definition: bezier.hpp:524
void distance_bound_test()
Definition: bezier_surface_test_suite.hpp:2240
data_type eqp_distance_bound(const bezier< data_type, dim__, tol__ > &bs) const
Definition: bezier.hpp:1257
point_type get_max() const
Definition: bounding_box.hpp:104
static dimension_type dimension()
Definition: bezier.hpp:144
point_type normal(const data_type &u, const data_type &v) const
Definition: bezier.hpp:896
point_type f(const data_type &t) const
Definition: bezier.hpp:324
Definition: bezier.hpp:38
bool demote_u(const geom::general::continuity &u_continuity_degree=geom::general::C0)
Definition: bezier.hpp:1053
void to_cubic_u()
Definition: bezier.hpp:1153
void get_bounding_box(bounding_box_type &bb) const
Definition: bezier.hpp:190
void promotion_to_test()
Definition: bezier_surface_test_suite.hpp:1374
void AddTests(const float &)
Definition: bezier_surface_test_suite.hpp:44
void resize(const index_type &u_dim, const index_type &v_dim)
Definition: bezier.hpp:166
void to_cubic_test()
Definition: bezier_surface_test_suite.hpp:2094
Eigen::Matrix< data_type, dim__, dim__ > rotation_matrix_type
Definition: bezier.hpp:47
point_type f_uvv(const data_type &u, const data_type &v) const
Definition: bezier.hpp:790
bezier_type::index_type index_type
Definition: bezier_surface_test_suite.hpp:37
Definition: math.hpp:25
index_type degree_u() const
Definition: bezier.hpp:146
tolerance_type tol
Definition: bezier_surface_test_suite.hpp:41
void gaussian_curvature(typename surface__::data_type &k, const surface__ &s, const typename surface__::data_type &u, const typename surface__::data_type &v)
Definition: curvature.hpp:206
void get_uconst_curve(curve_type &bc, const data_type &u) const
Definition: bezier.hpp:340
void reverse_v()
Definition: bezier.hpp:296
Definition: continuity.hpp:27
void split_test()
Definition: bezier_surface_test_suite.hpp:2330
point_type f_u(const data_type &u, const data_type &v) const
Definition: bezier.hpp:422
bezier_type::control_point_type control_point_type
Definition: bezier_surface_test_suite.hpp:34
bezier_type::data_type data_type
Definition: bezier_surface_test_suite.hpp:36
void assignment_test()
Definition: bezier_surface_test_suite.hpp:252
void get_curve_test()
Definition: bezier_surface_test_suite.hpp:2855
bezier_type::point_type point_type
Definition: bezier_surface_test_suite.hpp:35
bezier_surface_test_suite()
Definition: bezier_surface_test_suite.hpp:112
point_type get_control_point(const index_type &i, const index_type &j) const
Definition: bezier.hpp:175
void rotate(const rotation_matrix_type &rmat)
Definition: bezier.hpp:204
void principal_curvature(typename surface__::data_type &kmax, typename surface__::data_type &kmin, const surface__ &s, const typename surface__::data_type &u, const typename surface__::data_type &v)
Definition: curvature.hpp:234
Eigen::Matrix< data_type, 1, dim__ > point_type
Definition: bezier.hpp:43
void demotion_test()
Definition: bezier_surface_test_suite.hpp:1534
~bezier_surface_test_suite()
Definition: bezier_surface_test_suite.hpp:116
void mean_curvature(typename surface__::data_type &k, const surface__ &s, const typename surface__::data_type &u, const typename surface__::data_type &v)
Definition: curvature.hpp:178
point_type get_min() const
Definition: bounding_box.hpp:93
void octave_print(int figno, const bezier_type &bez) const
Definition: bezier_surface_test_suite.hpp:121
eli::geom::surface::bezier< data__, 3 > bezier_type
Definition: bezier_surface_test_suite.hpp:33
data__ data_type
Definition: bezier.hpp:42
point_type f(const data_type &u, const data_type &v) const
Definition: bezier.hpp:382
void swap_uv()
Definition: bezier.hpp:315
void transformation_test()
Definition: bezier_surface_test_suite.hpp:500
tol__ tolerance_type
Definition: bezier.hpp:46
void reverse_u()
Definition: bezier.hpp:277
bezier_type::curve_type curve_type
Definition: bezier_surface_test_suite.hpp:39
control_point_type::Index index_type
Definition: bezier.hpp:45
void get_vconst_curve(curve_type &bc, const data_type &v) const
Definition: bezier.hpp:361
void normal_test()
Definition: bezier_surface_test_suite.hpp:2632
void to_cubic_v()
Definition: bezier.hpp:1179
void promote_u_to(index_type target_degree)
Definition: bezier.hpp:975
bool demote_v(const geom::general::continuity &v_continuity_degree=geom::general::C0)
Definition: bezier.hpp:1103
void bounding_box_test()
Definition: bezier_surface_test_suite.hpp:343
void promote_u()
Definition: bezier.hpp:951
void reverse_test()
Definition: bezier_surface_test_suite.hpp:395
bezier_type::tolerance_type tolerance_type
Definition: bezier_surface_test_suite.hpp:38
void translate(const point_type &trans)
Definition: bezier.hpp:220
index_type degree_v() const
Definition: bezier.hpp:147
void evaluation_test()
Definition: bezier_surface_test_suite.hpp:677
Definition: bezier_surface_test_suite.hpp:30
void swap_test()
Definition: bezier_surface_test_suite.hpp:452
point_type f_vv(const data_type &u, const data_type &v) const
Definition: bezier.hpp:631
Definition: bezier.hpp:109
point_type f_v(const data_type &u, const data_type &v) const
Definition: bezier.hpp:473
void promotion_test()
Definition: bezier_surface_test_suite.hpp:1264
void AddTests(const double &)
Definition: bezier_surface_test_suite.hpp:66
point_type control_point_type
Definition: bezier.hpp:44
Definition: continuity.hpp:29
void promote_v()
Definition: bezier.hpp:1001
void derivative_1_test()
Definition: bezier_surface_test_suite.hpp:739
point_type f_uv(const data_type &u, const data_type &v) const
Definition: bezier.hpp:576