Claw 1.7.3
curve.tpp
1/*
2 CLAW - a C++ Library Absolutely Wonderful
3
4 CLAW is a free library without any particular aim but being useful to
5 anyone.
6
7 Copyright (C) 2005-2011 Julien Jorge
8
9 This library is free software; you can redistribute it and/or
10 modify it under the terms of the GNU Lesser General Public
11 License as published by the Free Software Foundation; either
12 version 2.1 of the License, or (at your option) any later version.
13
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public
20 License along with this library; if not, write to the Free Software
21 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
22 contact: julien.jorge@gamned.org
23*/
24/**
25 * \file claw/impl/curve.tpp
26 * \brief Implementation of claw::math::curve.
27 * \author Julien Jorge
28 */
29#include <boost/math/special_functions/cbrt.hpp>
30#include <boost/math/constants/constants.hpp>
31
32/*----------------------------------------------------------------------------*/
33/**
34 * \brief Constructor.
35 */
36template<typename C, typename Traits>
37claw::math::curve<C, Traits>::control_point::control_point()
38{
39
40} // curve::control_point::control_point()
41
42/*----------------------------------------------------------------------------*/
43/**
44 * \brief Constructor.
45 * \param p The position of this control point. It will be assigned to the input
46 * and output directions too.
47 */
48template<typename C, typename Traits>
49claw::math::curve<C, Traits>::control_point::control_point
50( const coordinate_type& p )
51 : m_position(p), m_input_direction(p), m_output_direction(p)
52{
53
54} // curve::control_point::control_point()
55
56/*----------------------------------------------------------------------------*/
57/**
58 * \brief Constructor.
59 * \param p The position of this control point.
60 * \param input_direction The point in the direction of which the curve enters
61 * \a p.
62 * \param output_direction The point in the direction of which the curve leaves
63 * \a p.
64 */
65template<typename C, typename Traits>
66claw::math::curve<C, Traits>::control_point::control_point
67( const coordinate_type& p, const coordinate_type& input_direction,
68 const coordinate_type& output_direction )
69 : m_position(p), m_input_direction(input_direction),
70 m_output_direction(output_direction)
71{
72
73} // curve::control_point::control_point()
74
75/*----------------------------------------------------------------------------*/
76/**
77 * \brief Get the position of this control point.
78 */
79template<typename C, typename Traits>
80const typename claw::math::curve<C, Traits>::control_point::coordinate_type&
81claw::math::curve<C, Traits>::control_point::get_position() const
82{
83 return m_position;
84} // curve::control_point::get_position()
85
86/*----------------------------------------------------------------------------*/
87/**
88 * \brief Get the point in the direction of which the curve enters this
89 * position.
90 */
91template<typename C, typename Traits>
92const typename claw::math::curve<C, Traits>::control_point::coordinate_type&
93claw::math::curve<C, Traits>::control_point::get_input_direction() const
94{
95 return m_input_direction;
96} // curve::control_point::get_input_direction()
97
98/*----------------------------------------------------------------------------*/
99/**
100 * \brief Get the point in the direction of which the curve leaves this
101 * position.
102 */
103template<typename C, typename Traits>
104const typename claw::math::curve<C, Traits>::control_point::coordinate_type&
105claw::math::curve<C, Traits>::control_point::get_output_direction() const
106{
107 return m_output_direction;
108} // curve::control_point::get_output_direction()
109
110
111
112
113
114/*----------------------------------------------------------------------------*/
115/**
116 * \brief Constructor.
117 * \param position The position of the point.
118 * \param s The section on which the point has been found.
119 * \param t The date of the point on the section.
120 */
121template<typename C, typename Traits>
122claw::math::curve<C, Traits>::section::resolved_point::resolved_point
123( const coordinate_type& position, const section& s, const double t )
124 : m_position(position), m_section(s), m_date(t)
125{
126
127} // curve::section::resolved_point::resolved_point()
128
129/*----------------------------------------------------------------------------*/
130/**
131 * \brief Get The position of the point.
132 */
133template<typename C, typename Traits>
134const
135typename claw::math::curve<C, Traits>::section::resolved_point::coordinate_type&
136claw::math::curve<C, Traits>::section::resolved_point::get_position() const
137{
138 return m_position;
139} // curve::section::resolved_point::get_position()
140
141/*----------------------------------------------------------------------------*/
142/**
143 * \brief Get the section on which the point has been found.
144 */
145template<typename C, typename Traits>
146const typename claw::math::curve<C, Traits>::section&
147claw::math::curve<C, Traits>::section::resolved_point::get_section() const
148{
149 return m_section;
150} // curve::section::::resolved_point::get_section()
151
152/*----------------------------------------------------------------------------*/
153/**
154 * \brief Get the date of the point on the section.
155 */
156template<typename C, typename Traits>
157double
158claw::math::curve<C, Traits>::section::resolved_point::get_date() const
159{
160 return m_date;
161} // curve::section::resolved_point::get_date()
162
163
164
165
166
167
168/*----------------------------------------------------------------------------*/
169/**
170 * \brief Constructor.
171 * \param origin The point at the beginning of the section.
172 * \param end The point at the end of the section.
173 */
174template<typename C, typename Traits>
175claw::math::curve<C, Traits>::section::section
176( const iterator_type& origin, const iterator_type& end )
177 : m_origin(origin), m_end(end)
178{
179
180} // curve::section::section()
181
182/*----------------------------------------------------------------------------*/
183/**
184 * \brief Get the point of this section at a given date.
185 * \param t The date at which the point is computed.
186 */
187template<typename C, typename Traits>
188typename claw::math::curve<C, Traits>::section::coordinate_type
189claw::math::curve<C, Traits>::section::get_point_at( double t ) const
190{
191 if ( m_origin == m_end )
192 return m_origin->get_position();
193
194 const value_type x
195 ( evaluate
196 ( t, traits_type::get_x(m_origin->get_position()),
197 traits_type::get_x(m_origin->get_output_direction()),
198 traits_type::get_x(m_end->get_input_direction()),
199 traits_type::get_x(m_end->get_position()) ) );
200 const value_type y
201 ( evaluate
202 ( t, traits_type::get_y(m_origin->get_position()),
203 traits_type::get_y(m_origin->get_output_direction()),
204 traits_type::get_y(m_end->get_input_direction()),
205 traits_type::get_y(m_end->get_position()) ) );
206
207 return traits_type::make_coordinate( x, y );
208} // curve::section::get_point_at()
209
210/*----------------------------------------------------------------------------*/
211/**
212 * \brief Get the direction of the tangent at a given date.
213 * \param t The date for which we want the tangent.
214 */
215template<typename C, typename Traits>
216typename claw::math::curve<C, Traits>::section::coordinate_type
217claw::math::curve<C, Traits>::section::get_tangent_at( double t ) const
218{
219 const value_type dx = evaluate_derived
220 ( t, traits_type::get_x(m_origin->get_position()),
221 traits_type::get_x(m_origin->get_output_direction()),
222 traits_type::get_x(m_end->get_input_direction()),
223 traits_type::get_x(m_end->get_position()) );
224
225 const value_type dy = evaluate_derived
226 ( t, traits_type::get_y(m_origin->get_position()),
227 traits_type::get_y(m_origin->get_output_direction()),
228 traits_type::get_y(m_end->get_input_direction()),
229 traits_type::get_y(m_end->get_position()) );
230
231 if ( dx == 0 )
232 return traits_type::make_coordinate( 0, (dy < 0) ? -1 : 1 );
233 else
234 return traits_type::make_coordinate( 1, dy / dx );
235} // curve::section::get_tangent_at()
236
237/*----------------------------------------------------------------------------*/
238/**
239 * \brief Get the points having the given x-coordinate on this section.
240 * \param x The coordinate for which we want the points.
241 * \param off_domain Tell the method to keep the points found at a date outside
242 * [0, 1].
243 */
244template<typename C, typename Traits>
245std::vector<typename claw::math::curve<C, Traits>::section::resolved_point>
246claw::math::curve<C, Traits>::section::get_point_at_x
247( value_type x, bool off_domain ) const
248{
249 std::vector<resolved_point> result;
250
251 if ( empty() )
252 return result;
253
254 const std::vector<double> roots
255 ( get_roots
256 ( x, traits_type::get_x(m_origin->get_position()),
257 traits_type::get_x(m_origin->get_output_direction()),
258 traits_type::get_x(m_end->get_input_direction()),
259 traits_type::get_x(m_end->get_position() ) ) );
260
261 for ( std::size_t i=0; i!=roots.size(); ++i )
262 result.push_back
263 ( resolved_point( get_point_at( roots[i] ), *this, roots[i] ) );
264
265 ensure_ends_in_points
266 ( result,
267 (x == m_origin->get_position().x), (x == m_end->get_position().x) );
268
269 if ( !off_domain )
270 return extract_domain_points( result );
271 else
272 return result;
273} // curve::section::get_point_at_x()
274
275/*----------------------------------------------------------------------------*/
276/**
277 * \brief Get an iterator on the control point at the origin of the section in
278 * the curve from which it was created.
279 */
280template<typename C, typename Traits>
281const typename claw::math::curve<C, Traits>::section::iterator_type&
282claw::math::curve<C, Traits>::section::get_origin() const
283{
284 return m_origin;
285} // curve::section::get_origin()
286
287/*----------------------------------------------------------------------------*/
288/**
289 * \brief Tell if there is no points on this section.
290 */
291template<typename C, typename Traits>
292bool claw::math::curve<C, Traits>::section::empty() const
293{
294 return m_origin == m_end;
295} // curve::section::empty()
296
297/*----------------------------------------------------------------------------*/
298/**
299 * \brief Get the value of the curve's equation on one dimension, at a given
300 * date.
301 * \param t The date at which the value us computed.
302 * \param origin The value on the computed dimension of the first point of the
303 * section of the curve.
304 * \param output_direction The value on the computed dimension of the point in
305 * the direction of which the curve leaves \a origin.
306 * \param input_direction The value on the computed dimension of the point in
307 * the direction of which the curve enters \a end.
308 * \param origin The value on the computed dimension of the last point of the
309 * section of the curve.
310 */
311template<typename C, typename Traits>
312typename claw::math::curve<C, Traits>::section::value_type
313claw::math::curve<C, Traits>::section::evaluate
314( double t, value_type origin, value_type output_direction,
315 value_type input_direction, value_type end ) const
316{
317 const double dt(1 - t);
318
319 return origin * dt * dt * dt
320 + 3 * output_direction * t * dt * dt
321 + 3 * input_direction * t * t * dt
322 + end * t * t * t;
323} // curve::section::evaluate()
324
325/*----------------------------------------------------------------------------*/
326/**
327 * \brief Get the value at a given date of the curve's derived equation on one
328 * dimension.
329 * \param t The date at which the value us computed.
330 * \param origin The value on the computed dimension of the first point of the
331 * section of the curve.
332 * \param output_direction The value on the computed dimension of the point in
333 * the direction of which the curve leaves \a origin.
334 * \param input_direction The value on the computed dimension of the point in
335 * the direction of which the curve enters \a end.
336 * \param origin The value on the computed dimension of the last point of the
337 * section of the curve.
338 */
339template<typename C, typename Traits>
340typename claw::math::curve<C, Traits>::section::value_type
341claw::math::curve<C, Traits>::section::evaluate_derived
342( double t, value_type origin, value_type output_direction,
343 value_type input_direction, value_type end ) const
344{
345 return - 3 * origin + 3 * output_direction
346 + (6 * origin - 12 * output_direction + 6 * input_direction) * t
347 + (-3 * origin + 9 * output_direction - 9 * input_direction + 3 * end)
348 * t * t;
349} // curve::section::evaluate_derived()
350
351/*----------------------------------------------------------------------------*/
352/**
353 * \brief Ensure that a vector of resolved_point contains some ends of the
354 * curve.
355 *
356 * The computation on the doubles values may produce approximated resolved
357 * points. If one end of the curve is expected to be in the resolved point, then
358 * this procedure replaces the nearest resolved point by a resolved point on
359 * the adequate end.
360 *
361 * \param p The points to filter.
362 * \param ensure_origin Tell to guarantee that the origin is present in \a p.
363 * \param ensure_end Tell to guarantee that the end is present in \a p.
364 */
365template<typename C, typename Traits>
366void claw::math::curve<C, Traits>::section::ensure_ends_in_points
367( std::vector<resolved_point>& p, bool ensure_origin, bool ensure_end ) const
368{
369 double min_distance_origin( std::numeric_limits<double>::max() );
370 double min_distance_end( std::numeric_limits<double>::max() );
371 std::size_t origin_index(p.size());
372 std::size_t end_index(p.size());
373
374 for ( std::size_t i=0; i!=p.size(); ++i )
375 {
376 const double distance_origin( std::abs( p[i].get_date() ) );
377
378 if ( distance_origin <= min_distance_origin )
379 {
380 min_distance_origin = distance_origin;
381 origin_index = i;
382 }
383
384 const double distance_end( std::abs( 1 - p[i].get_date() ) );
385
386 if ( distance_end <= min_distance_end )
387 {
388 min_distance_end = distance_end;
389 end_index = i;
390 }
391 }
392
393 if ( ensure_origin )
394 p[origin_index] = resolved_point( m_origin->get_position(), *this, 0.0 );
395
396 if ( ensure_end )
397 p[end_index] = resolved_point( m_end->get_position(), *this, 1.0 );
398} // curve::section::ensure_ends_in_points()
399
400/*----------------------------------------------------------------------------*/
401/**
402 * \brief Extract the points in the domain of the curve.
403 * \param p The points from which we extract the ones whose date is in [0, 1].
404 */
405template<typename C, typename Traits>
406std::vector<typename claw::math::curve<C, Traits>::section::resolved_point>
407claw::math::curve<C, Traits>::section::extract_domain_points
408( const std::vector<resolved_point>& p ) const
409{
410 std::vector<resolved_point> clean_result;
411
412 for ( std::size_t i=0; i!=p.size(); ++i )
413 if ( (p[i].get_date() >= 0) && (p[i].get_date() <= 1) )
414 clean_result.push_back( p[i] );
415
416 return clean_result;
417} // curve::section::extract_domain_points()
418
419/*----------------------------------------------------------------------------*/
420/**
421 * \brief Get the dates at which the curve passes at a given coordinate, on a
422 given dimension.
423 * \param x The coordinate for which we want the dates.
424 * \param origin The value on the computed dimension of the first point of the
425 * section of the curve.
426 * \param output_direction The value on the computed dimension of the point in
427 * the direction of which the curve leaves \a origin.
428 * \param input_direction The value on the computed dimension of the point in
429 * the direction of which the curve enters \a end.
430 * \param end The value on the computed dimension of the last point of the
431 * section of the curve.
432 */
433template<typename C, typename Traits>
434std::vector<double>
435claw::math::curve<C, Traits>::section::get_roots
436( value_type x, value_type origin, value_type output_direction,
437 value_type input_direction, value_type end ) const
438{
439 const value_type a
440 (-origin + 3 * output_direction - 3 * input_direction + end );
441 const value_type b( 3 * origin - 6 * output_direction + 3 * input_direction );
442 const value_type c( -3 * origin + 3 * output_direction );
443 const value_type d( origin - x );
444
445 if ( a == 0 )
446 return get_roots_degree_2(b, c, d);
447 else
448 return get_roots_degree_3(a, b, c, d);
449} // curve::section::get_roots()
450
451/*----------------------------------------------------------------------------*/
452/**
453 * \brief Get the dates at which the curve passes at a given coordinate, in the
454 * case where the equation is a reduced to a polynom of degree 2.
455 * \param a The coefficient of the square part of the equation.
456 * \param b The coefficient of the linear part of the equation.
457 * \param c The constant of the equation.
458 */
459template<typename C, typename Traits>
460std::vector<double>
461claw::math::curve<C, Traits>::section::get_roots_degree_2
462( value_type a, value_type b, value_type c ) const
463{
464 const value_type delta( b * b - 4 * a * c );
465
466 std::vector<double> result;
467
468 if ( delta == 0 )
469 result.push_back( - b / ( 2 * a ) );
470 else if ( delta > 0 )
471 {
472 result.push_back( (-b - std::sqrt(delta)) / (2 * a) );
473 result.push_back( (-b + std::sqrt(delta)) / (2 * a) );
474 }
475
476 return result;
477} // curve::section::get_roots_degree_2()
478
479/*----------------------------------------------------------------------------*/
480/**
481 * \brief Get the dates at which the curve passes at a given coordinate, in the
482 * case where the equation is a reduced to a polynom of degree 3.
483 * \param a The coefficient of the cubic part of the equation.
484 * \param b The coefficient of the square part of the equation.
485 * \param c The coefficient of the linear part of the equation.
486 * \param d The constant of the equation.
487 */
488template<typename C, typename Traits>
489std::vector<double>
490claw::math::curve<C, Traits>::section::get_roots_degree_3
491( value_type a, value_type b, value_type c, value_type d ) const
492{
493 // The following is the application of the method of Cardan
494
495 const value_type p( -(b * b) / (3.0 * a * a) + c / a );
496 const value_type q
497 ( ( b / (27.0 * a) )
498 * ( (2.0 * b * b) / (a * a)
499 - 9.0 * c / a )
500 + d / a );
501
502 const value_type delta( q * q + 4.0 * p * p * p / 27.0 );
503
504 std::vector<double> result;
505
506 if ( delta == 0 )
507 {
508 if ( p == 0 )
509 result.push_back(0);
510 else
511 {
512 result.push_back( 3.0 * q / p - b / (3.0 * a) );
513 result.push_back( - 3.0 * q / (2.0 * p) - b / (3.0 * a) );
514 }
515 }
516 else if ( delta > 0 )
517 {
518 result.push_back
519 ( boost::math::cbrt
520 ( (-q + std::sqrt(delta)) / 2.0 )
521 + boost::math::cbrt
522 ( (-q - std::sqrt(delta)) / 2.0 ) - b / (3.0 * a));
523 }
524 else
525 for ( std::size_t i=0; i!=3; ++i )
526 result.push_back
527 ( 2.0 * std::sqrt( -p / 3.0 )
528 * std::cos
529 ( std::acos( std::sqrt(27.0 / (- p * p * p)) * - q / 2.0 ) / 3.0
530 + 2.0 * i * boost::math::constants::pi<double>() / 3.0 )
531 - b / (3.0 * a));
532
533 return result;
534} // curve::section::get_roots_degree_3()
535
536
537
538
539
540
541/*----------------------------------------------------------------------------*/
542/**
543 * \brief Add a point at the end of the curve.
544 * \param p The point to add.
545 */
546template<typename C, typename Traits>
547void claw::math::curve<C, Traits>::push_back( const control_point& p )
548{
549 m_points.push_back(p);
550} // curve::push_back()
551
552/*----------------------------------------------------------------------------*/
553/**
554 * \brief Add a point at the beginning of the curve.
555 * \param p The point to add.
556 */
557template<typename C, typename Traits>
558void claw::math::curve<C, Traits>::push_front( const control_point& p )
559{
560 m_points.push_front(p);
561} // curve::push_front()
562
563/*----------------------------------------------------------------------------*/
564/**
565 * \brief Add a point before an other point of the curve.
566 * \param pos An iterator on the point before which the control point is added.
567 * \param p The point to add.
568 */
569template<typename C, typename Traits>
570void claw::math::curve<C, Traits>::insert
571( const iterator& pos, const control_point& p )
572{
573 m_points.insert( pos, p );
574} // curve::insert()
575
576/*----------------------------------------------------------------------------*/
577/**
578 * \brief Get the section of the curve starting at a given control point.
579 * \param pos An iterator of the control point at which the returned section
580 * begins.
581 */
582template<typename C, typename Traits>
583typename claw::math::curve<C, Traits>::section
584claw::math::curve<C, Traits>::get_section( const const_iterator& pos ) const
585{
586 const_iterator it(pos);
587 ++it;
588
589 if ( it == end() )
590 return section( pos, pos );
591 else
592 return section( pos, it );
593} // curve::get_section()
594
595/*----------------------------------------------------------------------------*/
596/**
597 * \brief Get the points having the given x-coordinate on the curve.
598 * \param x The coordinate for which we want the points.
599 * \param off_domain Tell the method to keep the points found at a date outside
600 * [0, 1].
601 *
602 * \todo Remove the duplicates in the result.
603 */
604template<typename C, typename Traits>
605std::vector< typename claw::math::curve<C, Traits>::section::resolved_point >
606claw::math::curve<C, Traits>::get_point_at_x
607( value_type x, bool off_domain ) const
608{
609 typedef std::vector<typename section::resolved_point> result_type;
610 result_type result;
611
612 for ( const_iterator it=begin(); it!=end(); ++it )
613 {
614 const section s( get_section(it) );
615
616 if ( !s.empty() )
617 {
618 const result_type new_points( s.get_point_at_x(x) );
619 result.insert( result.end(), new_points.begin(), new_points.end() );
620 }
621 }
622
623 if ( off_domain )
624 {
625 const result_type before( get_point_at_x_before_origin(x) );
626 result.insert( result.begin(), before.begin(), before.end() );
627
628 const result_type after( get_point_at_x_after_end(x) );
629 result.insert( result.end(), after.begin(), after.end() );
630 }
631
632 return result;
633} // curve::get_point_at_x()
634
635/*----------------------------------------------------------------------------*/
636/**
637 * \brief Get an iterator on the first control point.
638 */
639template<typename C, typename Traits>
640typename claw::math::curve<C, Traits>::iterator
641claw::math::curve<C, Traits>::begin()
642{
643 return m_points.begin();
644} // curve::begin()
645
646/*----------------------------------------------------------------------------*/
647/**
648 * \brief Get an iterator past the last control point.
649 */
650template<typename C, typename Traits>
651typename claw::math::curve<C, Traits>::iterator
652claw::math::curve<C, Traits>::end()
653{
654 return m_points.end();
655} // curve::end()
656
657/*----------------------------------------------------------------------------*/
658/**
659 * \brief Get an iterator on the first control point.
660 */
661template<typename C, typename Traits>
662typename claw::math::curve<C, Traits>::const_iterator
663claw::math::curve<C, Traits>::begin() const
664{
665 return m_points.begin();
666} // curve::begin()
667
668/*----------------------------------------------------------------------------*/
669/**
670 * \brief Get an iterator past the last control point.
671 */
672template<typename C, typename Traits>
673typename claw::math::curve<C, Traits>::const_iterator
674claw::math::curve<C, Traits>::end() const
675{
676 return m_points.end();
677} // curve::end()
678
679/*----------------------------------------------------------------------------*/
680/**
681 * \brief Get the points having the given x-coordinate before the origin of the
682 * curve.
683 * \param x The coordinate for which we want the points.
684 */
685template<typename C, typename Traits>
686std::vector< typename claw::math::curve<C, Traits>::section::resolved_point >
687claw::math::curve<C, Traits>::get_point_at_x_before_origin( value_type x ) const
688{
689 typedef std::vector<typename section::resolved_point> result_type;
690 result_type result;
691
692 const section s( get_section(begin()) );
693
694 if ( !s.empty() )
695 {
696 const result_type points( s.get_point_at_x(x, true) );
697
698 for ( std::size_t i(0); i!=points.size(); ++i )
699 if ( points[i].get_date() < 0 )
700 result.push_back( points[i] );
701 }
702
703 return result;
704} // curve::get_point_at_x_before_origin()
705
706/*----------------------------------------------------------------------------*/
707/**
708 * \brief Get the points having the given x-coordinate after the end of the
709 * curve.
710 * \param x The coordinate for which we want the points.
711 */
712template<typename C, typename Traits>
713std::vector< typename claw::math::curve<C, Traits>::section::resolved_point >
714claw::math::curve<C, Traits>::get_point_at_x_after_end( value_type x ) const
715{
716 typedef std::vector<typename section::resolved_point> result_type;
717 result_type result;
718
719 if ( m_points.size() < 2 )
720 return result;
721
722 const_iterator it(end());
723 std::advance(it, -2);
724
725 const section s( get_section( it ) );
726
727 if ( !s.empty() )
728 {
729 const result_type points( s.get_point_at_x(x, true) );
730
731 for ( std::size_t i(0); i!=points.size(); ++i )
732 if ( points[i].get_date() > 1 )
733 result.push_back( points[i] );
734 }
735
736 return result;
737} // curve::get_point_at_x_after_end()