Code-Eli  0.3.6
piecewise_cubic_spline_creator.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 eli_geom_curve_piecewise_spline_creator_hpp
14 #define eli_geom_curve_piecewise_spline_creator_hpp
15 
16 #include <iterator>
17 #include <vector>
18 
19 #include "eli/code_eli.hpp"
20 
21 #include "eli/mutil/fd/d1o2.hpp"
22 
27 
28 namespace eli
29 {
30  namespace geom
31  {
32  namespace curve
33  {
34  template<typename data__, unsigned short dim__, typename tol__>
35  class piecewise_cubic_spline_creator : public piecewise_creator_base<data__, dim__, tol__>
36  {
37  public:
43 
44  piecewise_cubic_spline_creator() : piecewise_creator_base<data_type, dim__, tolerance_type>(0, 0), control_point(0) {}
45  piecewise_cubic_spline_creator(const index_type &ns)
46  : piecewise_creator_base<data_type, dim__, tolerance_type>(ns, 0), control_point(3*ns+1) {}
48  : piecewise_creator_base<data_type, dim__, tolerance_type>(pcc), control_point(pcc.control_point) {}
49 
50  void get_segment_control_points(point_type &cp0, point_type &cp1,
51  point_type &cp2, point_type &cp3, const index_type &i) const
52  {
53  if ((3*i+1)<static_cast<index_type>(control_point.size()))
54  {
55  cp0=control_point[3*i];
56  cp1=control_point[3*i+1];
57  cp2=control_point[3*i+2];
58  cp3=control_point[3*i+3];
59  }
60  }
61 
62  void set_segment_control_points(const point_type &cp0, const point_type &cp1,
63  const point_type &cp2, const point_type &cp3, const index_type &i)
64  {
65  if ((3*i+1)<static_cast<index_type>(control_point.size()))
66  {
67  control_point[3*i] =cp0;
68  control_point[3*i+1]=cp1;
69  control_point[3*i+2]=cp2;
70  control_point[3*i+3]=cp3;
71  }
72  }
73 
74  void set_segment_point_slope(const point_type &p0, const point_type &m0,
75  const point_type &p1, const point_type &m1, const index_type &i)
76  {
77  if ((3*i+1)<static_cast<index_type>(control_point.size()))
78  {
79  point_type cp[4];
80  data_type dt(this->get_segment_dt(i));
81 
82  cp[0]=p0;
83  cp[1]=p0+(dt*m0/3);
84  cp[2]=p1-(dt*m1/3);
85  cp[3]=p1;
86  set_segment_control_points(cp[0], cp[1], cp[2], cp[3], i);
87  }
88  }
89 
91  {
92  typedef piecewise<bezier, data_type, dim__, tolerance_type> piecewise_curve_type;
93  typedef typename piecewise_curve_type::curve_type curve_type;
94  typedef typename piecewise_curve_type::error_code error_code;
95 
96  pc.clear();
97 
98  curve_type c(3);
99  error_code err;
100  index_type nsegs(this->get_number_segments());
101 
102  // do sanity check
103  if (control_point.size()!=(3*static_cast<size_t>(nsegs)+1))
104  {
105  assert(false);
106  return false;
107  }
108 
109  // set the start parameter
110  pc.set_t0(this->get_t0());
111 
112  // set each segment
113  for (index_type i=0; i<nsegs; ++i)
114  {
115  c.set_control_point(control_point[3*i ], 0);
116  c.set_control_point(control_point[3*i+1], 1);
117  c.set_control_point(control_point[3*i+2], 2);
118  c.set_control_point(control_point[3*i+3], 3);
119  err=pc.push_back(c, this->get_segment_dt(i));
120  if (err!=piecewise_curve_type::NO_ERRORS)
121  {
122  pc.clear();
123  pc.set_t0(0);
124  assert(false);
125  return false;
126  }
127  }
128 
129  return true;
130  }
131 
140  template<typename point_it__>
141  void set_chip(point_it__ itb, const eli::geom::general::continuity &end_cont)
142  {
143  index_type i, j, nsegs=this->get_number_segments(), npts;
144 
145  npts=nsegs;
146  if (end_cont==eli::geom::general::NOT_CONNECTED)
147  {
148  ++npts;
149  }
150 
151  // can't work with less than three points
152  if (npts<3)
153  {
154  assert(false);
155  return;
156  }
157 
158  point_it__ it, itm1, itp1, ite;
159  data_type tmp[3], t[3], dt;
160  point_type m[2];
162 
163  // set the end iterator if needed
164  ite=itb;
165  switch(end_cont)
166  {
168  {
169  break;
170  }
174  {
175  std::advance(ite, npts);
176  break;
177  }
178  default:
179  {
180  assert(false);
181  return;
182  break;
183  }
184  }
185 
186  // need to do first segment separately
187  // this is the mapping between iterators and indexes
188  // itm1 -> 0
189  // it ---> 1
190  // itp1 -> 2
191  itm1=itb;
192  it=itm1; ++it;
193  itp1=it; ++itp1;
194  i=0;
195 
196  // calculate the slope at the start of curve
197  switch (end_cont)
198  {
201  {
202  // set the parameter values
203  t[0]=this->get_t0();
204  t[1]=t[0]+this->get_segment_dt(0);
205  t[2]=t[1]+this->get_segment_dt(1);
206 
208  for (j=0; j<dim__; ++j)
209  {
210  tmp[0]=(*itm1)(j);
211  tmp[1]=(*it)(j);
212  tmp[2]=(*itp1)(j);
213  d1approx.evaluate(m[0](j), tmp, t);
214  }
215 
216  break;
217  }
220  {
221  point_it__ item1(ite);
222  --item1;
223 
224  // set the parameter values
225  t[0]=this->get_t0()-this->get_segment_dt(nsegs-1);
226  t[1]=this->get_t0();
227  t[2]=t[1]+this->get_segment_dt(0);
228 
230  for (j=0; j<dim__; ++j)
231  {
232  tmp[0]=(*item1)(j);
233  tmp[1]=(*itm1)(j);
234  tmp[2]=(*it)(j);
235  d1approx.evaluate(m[0](j), tmp, t);
236  }
237  break;
238  }
239  default:
240  {
241  return;
242  break;
243  }
244  }
245 
246  // calculate the slope at end of first segment
247  t[0]=this->get_t0();
248  t[1]=t[0]+this->get_segment_dt(0);
249  t[2]=t[1]+this->get_segment_dt(1);
251  for (j=0; j<dim__; ++j)
252  {
253  tmp[0]=(*itm1)(j);
254  tmp[1]=(*it)(j);
255  tmp[2]=(*itp1)(j);
256  d1approx.evaluate(m[1](j), tmp, t);
257  }
258 
259  // calculate the first set of control points
260  dt=this->get_segment_dt(i);
261 
262  set_segment_control_points(*itm1, (*itm1)+dt*m[0]/3, (*it)-dt*m[1]/3, *it, i);
263  m[0]=m[1];
264 
265  // do all interior segments
266  // this is the mapping between iterators and indexes
267  // itm1 -> i-1
268  // it ---> i
269  // itp1 -> i+1
270  ++itm1; ++it; ++itp1;
271  for (i=1; i<npts-2; ++i, ++itm1, ++it, ++itp1)
272  {
273  t[0]=t[1];
274  t[1]=t[2];
275  t[2]+=this->get_segment_dt(i+1);
276 
277  for (j=0; j<dim__; ++j)
278  {
279  tmp[0]=(*itm1)(j);
280  tmp[1]=(*it)(j);
281  tmp[2]=(*itp1)(j);
282  d1approx.evaluate(m[1](j), tmp, t);
283  }
284  dt=this->get_segment_dt(i);
285 
286  set_segment_control_points(*itm1, (*itm1)+dt*m[0]/3, (*it)-dt*m[1]/3, *it, i);
287  m[0]=m[1];
288  }
289 
290  // need to do the remaining segments separately
291  // this is the mapping between iterators and indexes
292  // itm1 -> npts-3
293  // it ---> npts-2
294  // itp1 -> npts-1
295  itp1=it;
296  it=itm1;
297  itm1=it; --itm1;
298  dt=this->get_segment_dt(npts-2);
299  switch (end_cont)
300  {
302  {
303  // last regular segment
305  for (j=0; j<dim__; ++j)
306  {
307  tmp[0]=(*itm1)(j);
308  tmp[1]=(*it)(j);
309  tmp[2]=(*itp1)(j);
310  d1approx.evaluate(m[1](j), tmp, t);
311  }
312 
313  set_segment_control_points(*it, (*it)+dt*m[0]/3, (*itp1)-dt*m[1]/3, *itp1, i);
314  break;
315  }
319  {
320  t[0]=t[1];
321  t[1]=t[2];
322  t[2]+=this->get_segment_dt(i+1);
323 
324  // need to do last regular segment separately
326  for (j=0; j<dim__; ++j)
327  {
328  tmp[0]=(*it)(j);
329  tmp[1]=(*itp1)(j);
330  tmp[2]=(*itb)(j);
331  d1approx.evaluate(m[1](j), tmp, t);
332  }
333 
334  set_segment_control_points(*it, (*it)+dt*m[0]/3, (*itp1)-dt*m[1]/3, *itp1, i);
335  m[0]=m[1];
336 
337  // need to do closing segment separately
338  dt=this->get_segment_dt(i+1);
339  if (end_cont==eli::geom::general::C0)
340  {
342  for (j=0; j<dim__; ++j)
343  {
344  tmp[0]=(*it)(j);
345  tmp[1]=(*itp1)(j);
346  tmp[2]=(*itb)(j);
347  d1approx.evaluate(m[1](j), tmp, t);
348  }
349  }
350  else
351  {
352  point_it__ it2(itb);
353  ++it2;
354 
355  t[0]=t[1];
356  t[1]=t[2];
357  t[2]+=this->get_segment_dt(0);
358 
360  for (j=0; j<dim__; ++j)
361  {
362  tmp[0]=(*itp1)(j);
363  tmp[1]=(*itb)(j);
364  tmp[2]=(*it2)(j);
365  d1approx.evaluate(m[1](j), tmp, t);
366  }
367  }
368 
369  set_segment_control_points(*itp1, (*itp1)+dt*m[0]/3, (*itb)-dt*m[1]/3, *itb, i+1);
370  break;
371  }
372  default:
373  {
374  assert(false);
375  return;
376  break;
377  }
378  }
379  }
380 
389  template<typename point_it__>
390  void set_cardinal(point_it__ itb, const data__ &c, const eli::geom::general::continuity &end_cont)
391  {
392  index_type i, nsegs=this->get_number_segments(), npts;
393 
394  npts=nsegs;
395  if (end_cont==eli::geom::general::NOT_CONNECTED)
396  {
397  ++npts;
398  }
399 
400  // can't work with less than three points
401  if (npts<3)
402  {
403  assert(false);
404  return;
405  }
406 
407  // check the parameter
408  if ((c<0) || (c>=1))
409  {
410  assert(false);
411  return;
412  }
413 
414  point_it__ it, itm1, itp1, ite;
415  data_type dt;
416  point_type m[2];
417 
418  // set the end iterator if needed
419  ite=itb;
420  switch(end_cont)
421  {
423  {
424  break;
425  }
429  {
430  std::advance(ite, npts);
431  break;
432  }
433  default:
434  {
435  assert(false);
436  return;
437  break;
438  }
439  }
440 
441  // need to do first segment separately
442  // this is the mapping between iterators and indexes
443  // itm1 -> 0
444  // it ---> 1
445  // itp1 -> 2
446  itm1=itb;
447  it=itm1; ++it;
448  itp1=it; ++itp1;
449  i=0;
450 
451  // calculate the slope at the start of curve
452  switch (end_cont)
453  {
456  {
457  // set the parameter values
458  m[0]=(1-c)*((*it)-(*itm1))/this->get_segment_dt(0);
459 
460  break;
461  }
464  {
465  point_it__ item1(ite);
466  --item1;
467 
468  // set the parameter values
469  m[0]=(1-c)*((*it)-(*item1))/(this->get_segment_dt(nsegs-1)+this->get_segment_dt(0));
470 
471  break;
472  }
473  default:
474  {
475  return;
476  break;
477  }
478  }
479 
480  // calculate the slope at end of first segment
481  m[1]=(1-c)*((*itp1)-(*itm1))/(this->get_segment_dt(0)+this->get_segment_dt(1));
482 
483  // calculate the first set of control points
484  dt=this->get_segment_dt(i);
485  set_segment_control_points(*itm1, (*itm1)+dt*m[0]/3, (*it)-dt*m[1]/3, *it, i);
486  m[0]=m[1];
487 
488  // do all interior segments
489  // this is the mapping between iterators and indexes
490  // itm1 -> i-1
491  // it ---> i
492  // itp1 -> i+1
493  ++itm1; ++it; ++itp1;
494  for (i=1; i<npts-2; ++i, ++itm1, ++it, ++itp1)
495  {
496  m[1]=(1-c)*((*itp1)-(*itm1))/(this->get_segment_dt(i)+this->get_segment_dt(i+1));
497 
498  dt=this->get_segment_dt(i);
499  set_segment_control_points(*itm1, (*itm1)+dt*m[0]/3, (*it)-dt*m[1]/3, *it, i);
500  m[0]=m[1];
501  }
502 
503  // need to do the remaining segments separately
504  // this is the mapping between iterators and indexes
505  // itm1 -> npts-3
506  // it ---> npts-2
507  // itp1 -> npts-1
508  itp1=it;
509  it=itm1;
510  itm1=it; --itm1;
511  dt=this->get_segment_dt(npts-2);
512  switch (end_cont)
513  {
515  {
516  // last regular segment
517  m[1]=(1-c)*((*itp1)-(*it))/this->get_segment_dt(nsegs-1);
518 
519  set_segment_control_points(*it, (*it)+dt*m[0]/3, (*itp1)-dt*m[1]/3, *itp1, i);
520  break;
521  }
525  {
526  // need to do last regular segment separately
527  m[1]=(1-c)*((*itb)-(*it))/(this->get_segment_dt(nsegs-2)+this->get_segment_dt(nsegs-1));
528 
529  set_segment_control_points(*it, (*it)+dt*m[0]/3, (*itp1)-dt*m[1]/3, *itp1, i);
530  m[0]=m[1];
531 
532  // need to do closing segment separately
533  dt=this->get_segment_dt(i+1);
534  if (end_cont==eli::geom::general::C0)
535  {
536  m[1]=(1-c)*((*itb)-(*itp1))/this->get_segment_dt(nsegs-1);
537  }
538  else
539  {
540  point_it__ it2(itb);
541  ++it2;
542 
543  m[1]=(1-c)*((*it2)-(*itp1))/(this->get_segment_dt(nsegs-1)+this->get_segment_dt(0));
544  }
545 
546  set_segment_control_points(*itp1, (*itp1)+dt*m[0]/3, (*itb)-dt*m[1]/3, *itb, i+1);
547  break;
548  }
549  default:
550  {
551  assert(false);
552  return;
553  break;
554  }
555  }
556  }
557 
563  template<typename point_it__>
564  void set_catmull_rom(point_it__ itb, const eli::geom::general::continuity &end_cont)
565  {
566  set_cardinal(itb, 0, end_cont);
567  }
568 
578  template<typename point_it__>
579  void set_kochanek_bartels(point_it__ itb, const data__ &tension, const data__ &bias,
580  const data__ &continuity, const eli::geom::general::continuity &end_cont)
581  {
582  index_type i, nsegs=this->get_number_segments(), npts;
583 
584  npts=nsegs;
585  if (end_cont==eli::geom::general::NOT_CONNECTED)
586  {
587  ++npts;
588  }
589 
590  // can't work with less than three points
591  if (npts<4)
592  {
593  assert(false);
594  return;
595  }
596 
597  // check some parameters
598  if ( (tension<-1) || (tension>1) )
599  {
600  assert(false);
601  return;
602  }
603  if ( (bias<-1) || (bias>1) )
604  {
605  assert(false);
606  return;
607  }
608  if ( (continuity<-1) || (continuity>1) )
609  {
610  assert(false);
611  return;
612  }
613 
614  point_it__ it, itm1, itp1, itp2, ite;
615  data_type dt;
616  point_type m[2];
617 
618  // set the end iterator if needed
619  ite=itb;
620  switch(end_cont)
621  {
623  {
624  break;
625  }
629  {
630  std::advance(ite, npts);
631  break;
632  }
633  default:
634  {
635  assert(false);
636  return;
637  break;
638  }
639  }
640 
641  // need to do first segment separately
642  // this is the mapping between iterators and indexes
643  // itm1 -> 0
644  // it ---> 1
645  // itp1 -> 2
646  // itp1 -> 3
647  itm1=itb;
648  it=itm1; ++it;
649  itp1=it; ++itp1;
650  itp2=itp1; ++itp2;
651  i=0;
652 
653  // calculate the slope at the start of curve
654  switch (end_cont)
655  {
658  {
659  // set the parameter values
660  m[0]=(1-tension)*(1-bias)*(1-continuity)*((*it)-(*itm1))/this->get_segment_dt(0);
661 
662  break;
663  }
666  {
667  point_it__ item1(ite);
668  --item1;
669 
670  // set the parameter values
671  m[0]=static_cast<data_type>(0.5)*(1-tension)*(1+bias)*(1+continuity)*((*itm1)-(*item1))/this->get_segment_dt(nsegs-1)
672  +static_cast<data_type>(0.5)*(1-tension)*(1-bias)*(1-continuity)*((*it)-(*itm1))/this->get_segment_dt(0);
673 
674  break;
675  }
676  default:
677  {
678  return;
679  break;
680  }
681  }
682 
683  // calculate the slope at end of first segment
684  m[1]=static_cast<data_type>(0.5)*(1-tension)*(1+bias)*(1-continuity)*((*it)-(*itm1))/this->get_segment_dt(i)
685  +static_cast<data_type>(0.5)*(1-tension)*(1-bias)*(1+continuity)*((*itp1)-(*it))/this->get_segment_dt(i+1);
686 
687  // calculate the first set of control points
688  dt=this->get_segment_dt(i);
689  set_segment_control_points(*itm1, (*itm1)+dt*m[0]/3, (*it)-dt*m[1]/3, *it, i);
690 
691  // do all interior segments
692  // this is the mapping between iterators and indexes
693  // itm1 -> i-1
694  // it ---> i
695  // itp1 -> i+1
696  for (i=1; i<npts-2; ++i, ++itm1, ++it, ++itp1, ++itp2)
697  {
698  dt=this->get_segment_dt(i);
699  m[0]=static_cast<data_type>(0.5)*(1-tension)*(1+bias)*(1+continuity)*((*it)-(*itm1))/this->get_segment_dt(i-1)
700  +static_cast<data_type>(0.5)*(1-tension)*(1-bias)*(1-continuity)*((*itp1)-(*it))/dt;
701  m[1]=static_cast<data_type>(0.5)*(1-tension)*(1+bias)*(1-continuity)*((*itp1)-(*it))/dt
702  +static_cast<data_type>(0.5)*(1-tension)*(1-bias)*(1+continuity)*((*itp2)-(*itp1))/this->get_segment_dt(i+1);
703 
704  set_segment_control_points(*it, (*it)+dt*m[0]/3, (*itp1)-dt*m[1]/3, *itp1, i);
705  }
706 
707  // need to do the remaining segments separately
708  // this is the mapping between iterators and indexes
709  // itm1 -> npts-3
710  // it ---> npts-2
711  // itp1 -> npts-1
712  dt=this->get_segment_dt(npts-2);
713  switch (end_cont)
714  {
716  {
717  // last regular segment
718  m[0]=static_cast<data_type>(0.5)*(1-tension)*(1+bias)*(1+continuity)*((*it)-(*itm1))/this->get_segment_dt(nsegs-2)
719  +static_cast<data_type>(0.5)*(1-tension)*(1-bias)*(1-continuity)*((*itp1)-(*it))/this->get_segment_dt(nsegs-1);
720  m[1]=(1-tension)*(1+bias)*(1-continuity)*((*itp1)-(*it))/this->get_segment_dt(nsegs-1);
721 
722  set_segment_control_points(*it, (*it)+dt*m[0]/3, (*itp1)-dt*m[1]/3, *itp1, i);
723  break;
724  }
728  {
729  // need to do last regular segment separately
730  m[0]=static_cast<data_type>(0.5)*(1-tension)*(1+bias)*(1+continuity)*((*it)-(*itm1))/this->get_segment_dt(nsegs-3)
731  +static_cast<data_type>(0.5)*(1-tension)*(1-bias)*(1-continuity)*((*itp1)-(*it))/this->get_segment_dt(nsegs-2);
732  m[1]=static_cast<data_type>(0.5)*(1-tension)*(1+bias)*(1-continuity)*((*itp1)-(*it))/this->get_segment_dt(nsegs-2)
733  +static_cast<data_type>(0.5)*(1-tension)*(1-bias)*(1+continuity)*((*itb)-(*itp1))/this->get_segment_dt(nsegs-1);
734 
735  set_segment_control_points(*it, (*it)+dt*m[0]/3, (*itp1)-dt*m[1]/3, *itp1, i);
736 
737  // need to do closing segment separately
738  dt=this->get_segment_dt(i+1);
739  m[0]=static_cast<data_type>(0.5)*(1-tension)*(1+bias)*(1+continuity)*((*itp1)-(*it))/this->get_segment_dt(nsegs-2)
740  +static_cast<data_type>(0.5)*(1-tension)*(1-bias)*(1-continuity)*((*itb)-(*itp1))/this->get_segment_dt(nsegs-1);
741  if (end_cont==eli::geom::general::C0)
742  {
743  m[1]=(1-tension)*(1+bias)*(1-continuity)*((*itb)-(*itp1))/this->get_segment_dt(nsegs-1);
744  }
745  else
746  {
747  point_it__ it2(itb);
748  ++it2;
749 
750  m[1]=static_cast<data_type>(0.5)*(1-tension)*(1+bias)*(1-continuity)*((*itb)-(*itp1))/this->get_segment_dt(nsegs-1)
751  +static_cast<data_type>(0.5)*(1-tension)*(1-bias)*(1+continuity)*((*it2)-(*itb))/this->get_segment_dt(0);
752  }
753 
754  set_segment_control_points(*itp1, (*itp1)+dt*m[0]/3, (*itb)-dt*m[1]/3, *itb, i+1);
755  break;
756  }
757  default:
758  {
759  assert(false);
760  return;
761  break;
762  }
763  }
764  }
765 
772  template<typename point_it__>
773  void set_cubic_spline(point_it__ itb)
774  {
775  index_type i, nseg(this->get_number_segments()), nunk(3*nseg+1);
776  Eigen::Matrix<data_type, Eigen::Dynamic, Eigen::Dynamic> M(nunk, nunk);
777  Eigen::Matrix<data_type, Eigen::Dynamic, dim__> b(nunk, dim__);
778 
779  // cannot apply this condition for less than 3 segments
780  if (nseg<3)
781  {
782  assert(false);
783  return;
784  }
785 
786  // create all of the common rows of system of equations
788 
789  // set the not-a-knot conditions (f''' is continuous between first two and last two segments)
790  data_type dt1_3, dt2_3;
791  dt1_3=this->get_segment_dt(0);
792  dt1_3*=dt1_3*dt1_3;
793  dt2_3=this->get_segment_dt(1);
794  dt2_3*=dt2_3*dt2_3;
795  M(1, 0)=-1/dt1_3;
796  M(1, 1)=3/dt1_3;
797  M(1, 2)=-3/dt1_3;
798  M(1, 3)=1/dt1_3+1/dt2_3;
799  M(1, 4)=-3/dt2_3;
800  M(1, 5)=3/dt2_3;
801  M(1, 6)=-1/dt2_3;
802  b.row(1).setZero();
803  dt1_3=this->get_segment_dt(nseg-2);
804  dt1_3*=dt1_3*dt1_3;
805  dt2_3=this->get_segment_dt(nseg-1);
806  dt2_3*=dt2_3*dt2_3;
807  M(nunk-2, nunk-7)=-1/dt1_3;
808  M(nunk-2, nunk-6)=3/dt1_3;
809  M(nunk-2, nunk-5)=-3/dt1_3;
810  M(nunk-2, nunk-4)=1/dt1_3+1/dt2_3;
811  M(nunk-2, nunk-3)=-3/dt2_3;
812  M(nunk-2, nunk-2)=3/dt2_3;
813  M(nunk-2, nunk-1)=-1/dt2_3;
814  b.row(nunk-2).setZero();
815 
816  // find the control points and set them
817  b=M.lu().solve(b);
818  for (i=0; i<nunk; ++i)
819  {
820  control_point[i]=b.row(i);
821  }
822  }
823 
830  template<typename point_it__>
831  void set_clamped_cubic_spline(point_it__ itb, const point_type &start_slope, const point_type &end_slope)
832  {
833  index_type i, nseg(this->get_number_segments()), nunk(3*nseg+1);
834  Eigen::Matrix<data_type, Eigen::Dynamic, Eigen::Dynamic> M(nunk, nunk);
835  Eigen::Matrix<data_type, Eigen::Dynamic, dim__> b(nunk, dim__);
836 
837  // create all of the common rows of system of equations
839 
840  // set the clamped conditions
841  data_type tmp;
842  tmp=this->get_segment_dt(0);
843  M(1, 0)=-3/tmp;
844  M(1, 1)=3/tmp;
845  b.row(1)=start_slope;
846  tmp=this->get_segment_dt(nseg-1);
847  M(nunk-2, nunk-2)=-3/tmp;
848  M(nunk-2, nunk-1)=3/tmp;
849  b.row(nunk-2)=end_slope;
850 
851  // find the control points and set them
852  b=M.lu().solve(b);
853  for (i=0; i<nunk; ++i)
854  {
855  control_point[i]=b.row(i);
856  }
857  }
858 
865  template<typename point_it__>
866  void set_natural_cubic_spline(point_it__ itb)
867  {
868  index_type i, nseg(this->get_number_segments()), nunk(3*nseg+1);
869  Eigen::Matrix<data_type, Eigen::Dynamic, Eigen::Dynamic> M(nunk, nunk);
870  Eigen::Matrix<data_type, Eigen::Dynamic, dim__> b(nunk, dim__);
871 
872  // create all of the common rows of system of equations
874 
875  // set the natural conditions (f''=0 at both ends)
876  M(1, 0)=1;
877  M(1, 1)=-2;
878  M(1, 2)=1;
879  b.row(1).setZero();
880  M(nunk-2, nunk-3)=1;
881  M(nunk-2, nunk-2)=-2;
882  M(nunk-2, nunk-1)=1;
883  b.row(nunk-2).setZero();
884 
885  // find the control points and set them
886  b=M.lu().solve(b);
887  for (i=0; i<nunk; ++i)
888  {
889  control_point[i]=b.row(i);
890  }
891  }
892 
899  template<typename point_it__>
900  void set_closed_cubic_spline(point_it__ itb)
901  {
902  index_type i, nseg(this->get_number_segments());
903  std::vector<point_type, Eigen::aligned_allocator<point_type>> pt(nseg+1);
904  point_it__ it;
905 
906  // copy over points
907  for (i=0, it=itb; i<nseg; ++i, ++it)
908  {
909  pt[i]=(*it);
910  }
911  pt[nseg]=pt[0];
912 
913  set_periodic_cubic_spline(pt.begin());
914  }
915 
922  template<typename point_it__>
923  void set_periodic_cubic_spline(point_it__ itb)
924  {
925  index_type i, nseg(this->get_number_segments()), nunk(3*nseg+1);
926  Eigen::Matrix<data_type, Eigen::Dynamic, Eigen::Dynamic> M(nunk, nunk);
927  Eigen::Matrix<data_type, Eigen::Dynamic, dim__> b(nunk, dim__);
928 
929  // create all of the common rows of system of equations
931 
932  // set the periodic conditions (f' and f'' are same at start and end)
933  data_type dt0(this->get_segment_dt(0)), dtn(this->get_segment_dt(nseg-1));
934  M(1, 0)=-1/dt0;
935  M(1, 1)=1/dt0;
936  M(1, nunk-2)=1/dtn;
937  M(1, nunk-1)=-1/dtn;
938  b.row(1).setZero();
939  M(nunk-2, 0)=1/dt0/dt0;
940  M(nunk-2, 1)=-2/dt0/dt0;
941  M(nunk-2, 2)=1/dt0/dt0;
942  M(nunk-2, nunk-3)=-1/dtn/dtn;
943  M(nunk-2, nunk-2)=2/dtn/dtn;
944  M(nunk-2, nunk-1)=-1/dtn/dtn;
945  b.row(nunk-2).setZero();
946 
947  // find the control points and set them
948  b=M.lu().solve(b);
949  for (i=0; i<nunk; ++i)
950  {
951  control_point[i]=b.row(i);
952  }
953  }
954 
955  private:
956  typedef std::vector<point_type, Eigen::aligned_allocator<point_type>> point_collection_type;
957 
959 
960  template<typename Derived1__, typename Derived2__, typename point_it__>
961  void create_cubic_spline_base_matrix(Eigen::MatrixBase<Derived1__> &M, Eigen::MatrixBase<Derived2__> &b, point_it__ itb)
962  {
963  index_type i, nseg(this->get_number_segments());
964  point_it__ it;
965 
966  // cycle through each segment and set conditions
967  it=itb;
968  M.setZero();
969  // set the C0 condition at start of segment 0
970  M(0, 0)=1;
971  b.row(0)=(*it);
972  for (++it, i=1; i<nseg; ++i, ++it)
973  {
974  data_type dt(this->get_segment_dt(i)), dtm1(this->get_segment_dt(i-1));
975 
976  // set the C2 condition at start of segment i
977  M(3*i-1, 3*i-2)=1/dtm1/dtm1;
978  M(3*i-1, 3*i-1)=-2/dtm1/dtm1;
979  M(3*i-1, 3*i)=1/dtm1/dtm1-1/dt/dt;
980  M(3*i-1, 3*i+1)=2/dt/dt;
981  M(3*i-1, 3*i+2)=-1/dt/dt;
982  b.row(3*i-1).setZero();
983 
984  // set the C0 condition at start of segment i
985  M(3*i, 3*i)=1;
986  b.row(3*i)=(*it);
987 
988  // set the C1 condition at start of segment i
989  M(3*i+1, 3*i-1)=1/dtm1;
990  M(3*i+1, 3*i)=-(1/dtm1+1/dt);
991  M(3*i+1, 3*i+1)=1/dt;
992  b.row(3*i+1).setZero();
993  }
994  // set the c0 condition at end of last segement
995  M(3*i, 3*i)=1;
996  b.row(3*i)=(*it);
997  }
998 
999  private:
1000  point_collection_type control_point;
1001  };
1002  }
1003  }
1004 }
1005 #endif
void set_cubic_spline(point_it__ itb)
Definition: piecewise_cubic_spline_creator.hpp:773
void set_segment_point_slope(const point_type &p0, const point_type &m0, const point_type &p1, const point_type &m1, const index_type &i)
Definition: piecewise_cubic_spline_creator.hpp:74
Definition: continuity.hpp:26
data__ data_type
Definition: piecewise_creator_base.hpp:33
Definition: math.hpp:20
Definition: continuity.hpp:28
piecewise_cubic_spline_creator()
Definition: piecewise_cubic_spline_creator.hpp:44
void get_segment_control_points(point_type &cp0, point_type &cp1, point_type &cp2, point_type &cp3, const index_type &i) const
Definition: piecewise_cubic_spline_creator.hpp:50
void set_stencil(const stencil &s)
Definition: d1o2.hpp:70
piecewise_cubic_spline_creator(const index_type &ns)
Definition: piecewise_cubic_spline_creator.hpp:45
void number_segments_changed()
Definition: piecewise_cubic_spline_creator.hpp:958
void set_chip(point_it__ itb, const eli::geom::general::continuity &end_cont)
Definition: piecewise_cubic_spline_creator.hpp:141
error_code push_back(const curve_type &curve, const data_type &dt=1.0)
Definition: piecewise.hpp:688
base_class_type::data_type data_type
Definition: piecewise_cubic_spline_creator.hpp:39
void set_cardinal(point_it__ itb, const data__ &c, const eli::geom::general::continuity &end_cont)
Definition: piecewise_cubic_spline_creator.hpp:390
point_collection_type control_point
Definition: piecewise_cubic_spline_creator.hpp:1000
void set_segment_control_points(const point_type &cp0, const point_type &cp1, const point_type &cp2, const point_type &cp3, const index_type &i)
Definition: piecewise_cubic_spline_creator.hpp:62
Definition: piecewise.hpp:244
Eigen::Matrix< data_type, 1, dim__ > point_type
Definition: piecewise_creator_base.hpp:34
std::vector< point_type, Eigen::aligned_allocator< point_type > > point_collection_type
Definition: piecewise_cubic_spline_creator.hpp:956
Definition: continuity.hpp:27
index_type get_number_segments() const
Definition: piecewise_creator_base.hpp:47
void set_kochanek_bartels(point_it__ itb, const data__ &tension, const data__ &bias, const data__ &continuity, const eli::geom::general::continuity &end_cont)
Definition: piecewise_cubic_spline_creator.hpp:579
Definition: continuity.hpp:33
Definition: piecewise_cubic_spline_creator.hpp:35
piecewise_creator_base< data__, dim__, tol__ > base_class_type
Definition: piecewise_cubic_spline_creator.hpp:38
base_class_type::point_type point_type
Definition: piecewise_cubic_spline_creator.hpp:40
virtual bool create(piecewise< bezier, data_type, dim__, tolerance_type > &pc) const
Definition: piecewise_cubic_spline_creator.hpp:90
void set_catmull_rom(point_it__ itb, const eli::geom::general::continuity &end_cont)
Definition: piecewise_cubic_spline_creator.hpp:564
Definition: d1o2.hpp:27
data_type get_segment_dt(const index_type &i) const
Definition: piecewise_creator_base.hpp:78
void clear()
Definition: piecewise.hpp:599
int evaluate(data__ &d, itphi__ itphi, const data__ &dx) const
Definition: d1o2.hpp:131
base_class_type::tolerance_type tolerance_type
Definition: piecewise_cubic_spline_creator.hpp:42
base_class_type::index_type index_type
Definition: piecewise_cubic_spline_creator.hpp:41
tol__ tolerance_type
Definition: piecewise_creator_base.hpp:36
void set_natural_cubic_spline(point_it__ itb)
Definition: piecewise_cubic_spline_creator.hpp:866
void create_cubic_spline_base_matrix(Eigen::MatrixBase< Derived1__ > &M, Eigen::MatrixBase< Derived2__ > &b, point_it__ itb)
Definition: piecewise_cubic_spline_creator.hpp:961
continuity
Definition: continuity.hpp:24
void set_t0(const data_type &t0_in)
Definition: piecewise.hpp:340
data_type get_t0() const
Definition: piecewise_creator_base.hpp:62
point_type::Index index_type
Definition: piecewise_creator_base.hpp:35
std::vector< data_type > dt
Definition: piecewise_creator_base.hpp:97
void set_closed_cubic_spline(point_it__ itb)
Definition: piecewise_cubic_spline_creator.hpp:900
Definition: piecewise_creator_base.hpp:30
void set_periodic_cubic_spline(point_it__ itb)
Definition: piecewise_cubic_spline_creator.hpp:923
void set_clamped_cubic_spline(point_it__ itb, const point_type &start_slope, const point_type &end_slope)
Definition: piecewise_cubic_spline_creator.hpp:831
piecewise_cubic_spline_creator(const piecewise_cubic_spline_creator< data_type, dim__, tolerance_type > &pcc)
Definition: piecewise_cubic_spline_creator.hpp:47