Code-Eli  0.3.6
piecewise_superellipse_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_superellipse_creator_hpp
14 #define eli_geom_curve_piecewise_superellipse_creator_hpp
15 
16 #include "eli/code_eli.hpp"
17 
21 
22 namespace eli
23 {
24  namespace geom
25  {
26  namespace curve
27  {
28  template<typename data__, unsigned short dim__, typename tol__>
29  class piecewise_2d_curve_creator : public piecewise_creator_base<data__, dim__, tol__>
30  {
31  public:
37 
38  piecewise_2d_curve_creator() : piecewise_creator_base<data_type, dim__, tolerance_type>(4, 0)
39  {
40  origin.setZero();
41  }
42  piecewise_2d_curve_creator(const index_type &ns, const data_type &tt0) : piecewise_creator_base<data_type, dim__, tolerance_type>(ns, tt0)
43  {
44  origin.setZero();
45  }
47  : piecewise_creator_base<data_type, dim__, tolerance_type>(p2dc), origin(p2dc.origin) {}
49 
50  void set_origin(const point_type &orig) {origin=orig;}
51  point_type get_origin() const {return origin;}
52 
53  protected:
54  virtual void fun(point_type &f, const data_type &t) const = 0;
55  virtual void fun(point_type &f, point_type &fp, const data_type &t) const = 0;
56 
57  private:
58  point_type origin;
59  };
60 
61 
62  template<typename data__, unsigned short dim__, typename tol__>
63  class piecewise_superellipse_creator : public piecewise_2d_curve_creator<data__, dim__, tol__>
64  {
65  public:
70 
71  piecewise_superellipse_creator() : piecewise_2d_curve_creator<data_type, dim__, tolerance_type>(4, 0), a(1), b(1), m(2), n(2), max_degree(3) {}
72  piecewise_superellipse_creator(const index_type &ns) : piecewise_2d_curve_creator<data_type, dim__, tolerance_type>(ns, 0), a(1), b(1), m(2), n(2), max_degree(3) {}
74  : piecewise_2d_curve_creator<data_type, dim__, tolerance_type>(ppc), a(ppc.a), b(ppc.b), m(ppc.m), n(ppc.n), max_degree(ppc.max_degree) {}
76 
77  void set_axis(const data_type &aa, const data_type &bb)
78  {
79  set_a_axis(aa);
80  set_b_axis(bb);
81  }
82  void get_axis(data_type &aa, data_type &bb) const
83  {
84  aa = a;
85  bb = b;
86  }
87 
88  void set_a_axis(const data_type &aa)
89  {
90  if (aa>=0)
91  {
92  a=aa;
93  }
94  }
95  const data_type & get_a_axis() const {return a;}
96 
97  void set_b_axis(const data_type &bb)
98  {
99  if (bb>=0)
100  {
101  b=bb;
102  }
103  }
104  const data_type & get_b_axis() const {return b;}
105 
106  void set_exponents(const data_type &mm, const data_type &nn)
107  {
108  set_m_exponent(mm);
109  set_n_exponent(nn);
110  }
111  void get_exponents(data_type &mm, data_type &nn) const
112  {
113  mm = m;
114  nn = n;
115  }
116 
117  void set_m_exponent(const data_type &mm)
118  {
119  if (mm>0)
120  {
121  m=mm;
122  }
123  }
124  const data_type & get_m_exponent() const {return m;}
125 
126  void set_n_exponent(const data_type &nn)
127  {
128  if (nn>0)
129  {
130  n=nn;
131  }
132  }
133  const data_type & get_n_exponent() const {return n;}
134 
135  void set_max_degree(const index_type &md)
136  {
137  if (md>2)
138  {
139  max_degree=md;
140  }
141  }
142  index_type get_max_degree() const {return max_degree;}
143 
145  {
146  typedef piecewise<bezier, data_type, dim__, tolerance_type> piecewise_curve_type;
147  typedef typename piecewise_curve_type::curve_type curve_type;
148  typedef typename piecewise_curve_type::control_point_type control_point_type;
149  typedef typename piecewise_curve_type::error_code error_code;
150  typedef typename curve_type::fit_container_type fit_container_type;
151 
152  pc.clear();
153 
154  // TODO: should this be an adaptive thing instead of just start with max?
155 
156  curve_type crv(max_degree);
157  index_type nsegs(this->get_number_segments());
158  error_code err;
159 
160  // set the start parameter
161  pc.set_t0(this->get_t0());
162 
163  // TODO: REMOVE THIS AFTER TESTING
164  assert(nsegs%2==0);
165 
166  // if even number of segments then can create top half
167  if (nsegs%2==0)
168  {
169  piecewise_curve_type pc_bottom;
170  index_type iseg, nseg_top(nsegs/2);
171 
172  // if have multiple of 4 segments then can create first quarter
173  if (nsegs%4==0)
174  {
175  piecewise_curve_type pc_left;
176  index_type nseg_first(nseg_top/2), nsample_pts(nseg_first*(max_degree)+1);
177  std::vector<point_type> f(nsample_pts);
178 
179  // ensure that the first and last points are set appropriately
180  f[0].setZero();
181  f[0].x()=a;
182  f[nsample_pts-1].setZero();
183  f[nsample_pts-1].y()=b;
184 
185  // create rest of sample points
186  for (index_type i=1; i<(nsample_pts-1); ++i)
187  {
188  data_type t, argx, argy, tmp;
189 
190  // to get the rapid parameterization variation, find the parameters that
191  // produce uniform spacing in x and y and then average the two
192  tmp=static_cast<data_type>(i)/(nsample_pts-1);
193  argx=std::pow(1-tmp, m/2);
194  argx=std::min(argx, static_cast<data_type>(1));
195  argx=std::max(argx, static_cast<data_type>(0));
196  argy=std::pow(tmp, n/2);
197  argy=std::min(argy, static_cast<data_type>(1));
198  argy=std::max(argy, static_cast<data_type>(0));
199 
200  // TODO: Fix this to get the parameter that better captures extreme n and m cases
201  t=(std::acos(argx)+std::asin(argy))/(2*eli::constants::math<data__>::two_pi());
202 
203  fun(f[i], t);
204  }
205 
206  // create first quarter
207  for (iseg=0; iseg<nseg_first; ++iseg)
208  {
209  control_point_type cp;
210  fit_container_type fcon;
211  crv.clear();
212 
213  fcon.set_points(f.begin()+(iseg*max_degree), f.begin()+((iseg+1)*max_degree+1));
214 
215  crv.interpolate(fcon);
216 
217  // check the first slope to make sure it is reasonable
218  if (iseg==0)
219  {
220  bool need_set(false);
221  control_point_type cp(crv.get_control_point(1));
222 
223  if (cp.y()<0)
224  {
225  cp.y()=0;
226  need_set=true;
227  }
228  if (cp.x()>a)
229  {
230  cp.x()=a;
231  need_set=true;
232  }
233  if (need_set)
234  {
235  crv.set_control_point(cp, 1);
236  }
237  }
238 
239  // check the last slope to make sure it is reasonable
240  if (iseg==(nseg_first-1))
241  {
242  bool need_set(false);
243  control_point_type cp(crv.get_control_point(max_degree-1));
244 
245  if (cp.x()<0)
246  {
247  cp.x()=0;
248  need_set=true;
249  }
250  if (cp.y()>b)
251  {
252  cp.y()=b;
253  need_set=true;
254  }
255  if (need_set)
256  {
257  crv.set_control_point(cp, max_degree-1);
258  }
259  }
260 
261  err=pc.push_back(crv, this->get_segment_dt(iseg));
262  if (err!=piecewise_curve_type::NO_ERRORS)
263  {
264  std::cout << "error number: " << err << std::endl;
265  assert(false);
266  pc.clear();
267  pc.set_t0(0);
268  return false;
269  }
270  }
271 
272  // mirror to create second quarter
273  pc_left=pc;
274  pc_left.reflect_yz();
275  pc_left.reverse();
276 
277 // for (int ii=0; ii<=max_degree; ++ii)
278 // {
279 // std::cout << "ii = " << ii << "\tf_r=" << pc.f(ii/(1.0*max_degree))
280 // << "\tf_l=" << pc_left.f(1-ii/(1.0*max_degree)) << std::endl;
281 // }
282 // assert(false);
283 
284  // push back left side
285  curve_type c;
286  index_type rnseg(pc.number_segments());
287  for (iseg=0; iseg<pc_left.number_segments(); ++iseg)
288  {
289  pc_left.get(c, iseg);
290  err=pc.push_back(c, this->get_segment_dt(rnseg+iseg));
291  if (err!=piecewise_curve_type::NO_ERRORS)
292  {
293  std::cout << "error number: " << err << std::endl;
294  assert(false);
295  pc.clear();
296  pc.set_t0(0);
297  return false;
298  }
299  }
300  }
301  // else create entire top half
302  else
303  {
304  index_type nsample_pts(nseg_top*(max_degree)+1), nhalf(nsample_pts/2);
305  std::vector<point_type> f(nsample_pts);
306 
307  // ensure that the first and last points are set appropriately
308  f[0].setZero();
309  f[0].x()=a;
310  f[nsample_pts-1].setZero();
311  f[nsample_pts-1].x()=-a;
312  if (nsample_pts%2==1)
313  {
314  f[nhalf].setZero();
315  f[nhalf].y()=b;
316  }
317 
318  // create rest of sample points
319  for (index_type i=1; i<nhalf; ++i)
320  {
321  data_type t, argx, argy, tmp;
322 
323  // to get the rapid parameterization variation, find the parameters that
324  // produce uniform spacing in x and y and then average the two
325  tmp=static_cast<data_type>(i)/(nhalf);
326  argx=std::pow(1-tmp, m/2);
327  argx=std::min(argx, static_cast<data_type>(1));
328  argx=std::max(argx, static_cast<data_type>(0));
329  argy=std::pow(tmp, n/2);
330  argy=std::min(argy, static_cast<data_type>(1));
331  argy=std::max(argy, static_cast<data_type>(0));
332 
333  // TODO: Fix this to get the parameter that better captures extreme n and m cases
334  t=(std::acos(argx)+std::asin(argy))/(2*eli::constants::math<data__>::two_pi());
335 // std::cout << "argx=" << argx << "\targy=" << argy << "\tt=" << t << std::endl;
336 
337  fun(f[i], t);
338  fun(f[nsample_pts-i-1], static_cast<data_type>(0.5)-t);
339  }
340 
341 // if (max_degree==6)
342 // {
343 // for (index_type i=0; i<nsample_pts; ++i)
344 // {
345 // std::cout << "f[" << i << "] = " << f[i] << std::endl;
346 // }
347 // }
348 
349  // create first quarter
350  for (iseg=0; iseg<nseg_top; ++iseg)
351  {
352  control_point_type cp;
353  fit_container_type fcon;
354  crv.clear();
355 
356  fcon.set_points(f.begin()+(iseg*max_degree), f.begin()+((iseg+1)*max_degree+1));
357 
358  crv.interpolate(fcon);
359 
360  // check the first slope to make sure it is reasonable
361  if (iseg==0)
362  {
363  bool need_set(false);
364  control_point_type cp(crv.get_control_point(1));
365 
366  if (cp.y()<0)
367  {
368  cp.y()=0;
369  need_set=true;
370  }
371  if (cp.x()>a)
372  {
373  cp.x()=a;
374  need_set=true;
375  }
376  if (need_set)
377  {
378  crv.set_control_point(cp, 1);
379  }
380  }
381 
382  // check the last slope to make sure it is reasonable
383  if (iseg==(nseg_top-1))
384  {
385  bool need_set(false);
386  control_point_type cp(crv.get_control_point(max_degree-1));
387 
388  if (cp.y()<0)
389  {
390  cp.y()=0;
391  need_set=true;
392  }
393  if (cp.x()<-a)
394  {
395  cp.x()=-a;
396  need_set=true;
397  }
398  if (need_set)
399  {
400  crv.set_control_point(cp, max_degree-1);
401  }
402  }
403 
404  err=pc.push_back(crv, this->get_segment_dt(iseg));
405  if (err!=piecewise_curve_type::NO_ERRORS)
406  {
407  std::cout << "error number: " << err << std::endl;
408  assert(false);
409  pc.clear();
410  pc.set_t0(0);
411  return false;
412  }
413  }
414  }
415 
416  // mirror for bottom half
417  pc_bottom=pc;
418  pc_bottom.reflect_xz();
419  pc_bottom.reverse();
420 
421  // push back the bottom curve
422  curve_type c;
423  index_type rtseg(pc.number_segments());
424  for (iseg=0; iseg<pc_bottom.number_segments(); ++iseg)
425  {
426  pc_bottom.get(c, iseg);
427  err=pc.push_back(c, this->get_segment_dt(rtseg+iseg));
428  if (err!=piecewise_curve_type::NO_ERRORS)
429  {
430  std::cout << "error number: " << err << std::endl;
431  assert(false);
432  pc.clear();
433  pc.set_t0(0);
434  return false;
435  }
436  }
437  }
438  // else odd number of segments
439  else
440  {
441  assert(false);
442  return false;
443  }
444 
445  // translate to center on origin
446  pc.translate(this->get_origin());
447 
448  return true;
449  }
450 
451  protected:
452  // TODO: These should be in the base class that is shared with circle, ellipse, and other 2D curves
453  void set_x_radius(const data_type &xr)
454  {
455  if(xr>=0)
456  {
457  xradius=xr;
458  }
459  else
460  {
461  assert(false);
462  }
463  }
464  const data_type & get_x_radius() const {return xradius;}
465  void set_y_radius(const data_type &yr)
466  {
467  if (yr>=0)
468  {
469  yradius=yr;
470  }
471  else
472  {
473  assert(false);
474  }
475  }
476  const data_type & get_y_radius() const {return yradius;}
477  // TODO: end note
478 
479  protected:
480  virtual void fun(point_type &f, const data_type &t) const
481  {
482  // short circuit if given bad parameter
483  if ((t<0) || (t>1))
484  return;
485 
486  // zero out point
487  f.setZero();
488 
489  // determine the sign terms
490  data_type sign_cos(1), sign_sin(1);
491 
492  if ((t>0.25) && (t<0.5))
493  {
494  sign_cos=-1;
495  }
496  else if ((t>0.5) && (t>0.75))
497  {
498  sign_cos=-1;
499  sign_sin=-1;
500  }
501  else if (t>0.75)
502  {
503  sign_sin=-1;
504  }
505 
506  data_type theta(eli::constants::math<data__>::two_pi()*t);
507  data_type abs_cos(std::abs(std::cos(theta)));
508  data_type abs_sin(std::abs(std::sin(theta)));
509 
510  // calculate the function
511  f.x()=a*sign_cos*std::pow(abs_cos, 2/m);
512  f.y()=b*sign_sin*std::pow(abs_sin, 2/n);
513  }
514 
515  virtual void fun(point_type &f, point_type &fp, const data_type &t) const
516  {
517  // short circuit if given bad parameter
518  if ((t<0) || (t>1))
519  return;
520 
521  // zero out point
522  f.setZero();
523  fp.setZero();
524 
525  // determine the sign terms
526  data_type sign_cos(1), sign_sin(1);
527 
528  if ((t>0.25) && (t<0.5))
529  {
530  sign_cos=-1;
531  }
532  else if ((t>0.5) && (t>0.75))
533  {
534  sign_cos=-1;
535  sign_sin=-1;
536  }
537  else if (t>0.75)
538  {
539  sign_sin=-1;
540  }
541 
542  data_type theta(eli::constants::math<data__>::two_pi()*t);
543  data_type abs_cos(std::abs(std::cos(theta)));
544  data_type abs_sin(std::abs(std::sin(theta)));
545 
546  // calculate the function
547  f.x()=a*sign_cos*std::pow(abs_cos, 2/m);
548  f.y()=b*sign_sin*std::pow(abs_sin, 2/n);
549 
550  // calculate the 1st derivatives
551  fp.x()=(-2*eli::constants::math<data__>::two_pi()*a/m)*std::sin(theta)*std::pow(abs_cos, 2/m-1);
552  fp.y()=(2*eli::constants::math<data__>::two_pi()*b/n)*std::cos(theta)*std::pow(abs_sin, 2/n-1);
553  }
554 
555  private:
556  data_type xradius, yradius;
557 
558  data_type a, b, m, n;
559  index_type max_degree;
560  };
561  }
562  }
563 }
564 #endif
void get_exponents(data_type &mm, data_type &nn) const
Definition: piecewise_superellipse_creator.hpp:111
void set_b_axis(const data_type &bb)
Definition: piecewise_superellipse_creator.hpp:97
virtual void fun(point_type &f, point_type &fp, const data_type &t) const
Definition: piecewise_superellipse_creator.hpp:515
piecewise_2d_curve_creator(const piecewise_2d_curve_creator< data_type, dim__, tolerance_type > &p2dc)
Definition: piecewise_superellipse_creator.hpp:46
virtual bool create(piecewise< bezier, data_type, dim__, tolerance_type > &pc) const
Definition: piecewise_superellipse_creator.hpp:144
void get_axis(data_type &aa, data_type &bb) const
Definition: piecewise_superellipse_creator.hpp:82
data__ data_type
Definition: piecewise_creator_base.hpp:33
Definition: math.hpp:20
const data_type & get_n_exponent() const
Definition: piecewise_superellipse_creator.hpp:133
data_type a
Definition: piecewise_superellipse_creator.hpp:558
piecewise_2d_curve_creator< data__, dim__, tol__ >::tolerance_type tolerance_type
Definition: piecewise_superellipse_creator.hpp:69
void set_max_degree(const index_type &md)
Definition: piecewise_superellipse_creator.hpp:135
point_type get_origin() const
Definition: piecewise_superellipse_creator.hpp:51
error_code push_back(const curve_type &curve, const data_type &dt=1.0)
Definition: piecewise.hpp:688
const data_type & get_b_axis() const
Definition: piecewise_superellipse_creator.hpp:104
const data_type & get_y_radius() const
Definition: piecewise_superellipse_creator.hpp:476
piecewise_2d_curve_creator(const index_type &ns, const data_type &tt0)
Definition: piecewise_superellipse_creator.hpp:42
index_type number_segments() const
Definition: piecewise.hpp:419
void set_exponents(const data_type &mm, const data_type &nn)
Definition: piecewise_superellipse_creator.hpp:106
void reflect_xz()
Definition: piecewise.hpp:559
void set_n_exponent(const data_type &nn)
Definition: piecewise_superellipse_creator.hpp:126
piecewise_superellipse_creator(const index_type &ns)
Definition: piecewise_superellipse_creator.hpp:72
Definition: math.hpp:25
Definition: piecewise.hpp:244
Eigen::Matrix< data_type, 1, dim__ > point_type
Definition: piecewise_creator_base.hpp:34
piecewise_creator_base< data__, dim__, tol__ > base_class_type
Definition: piecewise_superellipse_creator.hpp:32
piecewise_superellipse_creator(const piecewise_superellipse_creator< data_type, dim__, tolerance_type > &ppc)
Definition: piecewise_superellipse_creator.hpp:73
index_type max_degree
Definition: piecewise_superellipse_creator.hpp:559
data_type yradius
Definition: piecewise_superellipse_creator.hpp:556
index_type get_number_segments() const
Definition: piecewise_creator_base.hpp:47
point_type origin
Definition: piecewise_superellipse_creator.hpp:58
piecewise_superellipse_creator()
Definition: piecewise_superellipse_creator.hpp:71
void translate(const point_type &trans)
Definition: piecewise.hpp:491
void set_m_exponent(const data_type &mm)
Definition: piecewise_superellipse_creator.hpp:117
data_type m
Definition: piecewise_superellipse_creator.hpp:558
void set_a_axis(const data_type &aa)
Definition: piecewise_superellipse_creator.hpp:88
const data_type & get_m_exponent() const
Definition: piecewise_superellipse_creator.hpp:124
const data_type & get_x_radius() const
Definition: piecewise_superellipse_creator.hpp:464
~piecewise_2d_curve_creator()
Definition: piecewise_superellipse_creator.hpp:48
piecewise_2d_curve_creator< data__, dim__, tol__ >::index_type index_type
Definition: piecewise_superellipse_creator.hpp:67
base_class_type::point_type point_type
Definition: piecewise_superellipse_creator.hpp:34
data_type xradius
Definition: piecewise_superellipse_creator.hpp:556
piecewise_2d_curve_creator< data__, dim__, tol__ >::point_type point_type
Definition: piecewise_superellipse_creator.hpp:68
data_type get_segment_dt(const index_type &i) const
Definition: piecewise_creator_base.hpp:78
void clear()
Definition: piecewise.hpp:599
piecewise_2d_curve_creator()
Definition: piecewise_superellipse_creator.hpp:38
virtual void fun(point_type &f, const data_type &t) const
Definition: piecewise_superellipse_creator.hpp:480
void set_x_radius(const data_type &xr)
Definition: piecewise_superellipse_creator.hpp:453
tol__ tolerance_type
Definition: piecewise_creator_base.hpp:36
void set_origin(const point_type &orig)
Definition: piecewise_superellipse_creator.hpp:50
void set_axis(const data_type &aa, const data_type &bb)
Definition: piecewise_superellipse_creator.hpp:77
virtual void fun(point_type &f, const data_type &t) const =0
void set_t0(const data_type &t0_in)
Definition: piecewise.hpp:340
void set_y_radius(const data_type &yr)
Definition: piecewise_superellipse_creator.hpp:465
base_class_type::data_type data_type
Definition: piecewise_superellipse_creator.hpp:33
index_type get_max_degree() const
Definition: piecewise_superellipse_creator.hpp:142
data_type get_t0() const
Definition: piecewise_creator_base.hpp:62
base_class_type::index_type index_type
Definition: piecewise_superellipse_creator.hpp:35
point_type::Index index_type
Definition: piecewise_creator_base.hpp:35
void reflect_yz()
Definition: piecewise.hpp:569
piecewise_2d_curve_creator< data__, dim__, tol__ >::data_type data_type
Definition: piecewise_superellipse_creator.hpp:66
data_type b
Definition: piecewise_superellipse_creator.hpp:558
Definition: piecewise_superellipse_creator.hpp:63
~piecewise_superellipse_creator()
Definition: piecewise_superellipse_creator.hpp:75
data_type n
Definition: piecewise_superellipse_creator.hpp:558
base_class_type::tolerance_type tolerance_type
Definition: piecewise_superellipse_creator.hpp:36
Definition: piecewise_superellipse_creator.hpp:29
Definition: piecewise_creator_base.hpp:30
const data_type & get_a_axis() const
Definition: piecewise_superellipse_creator.hpp:95