2 CLAW - a C++ Library Absolutely Wonderful
4 CLAW is a free library without any particular aim but being useful to
7 Copyright (C) 2005-2011 Julien Jorge
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.
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.
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
25 * \file claw/impl/curve.tpp
26 * \brief Implementation of claw::math::curve.
27 * \author Julien Jorge
29#include <boost/math/special_functions/cbrt.hpp>
30#include <boost/math/constants/constants.hpp>
32/*----------------------------------------------------------------------------*/
36template<typename C, typename Traits>
37claw::math::curve<C, Traits>::control_point::control_point()
40} // curve::control_point::control_point()
42/*----------------------------------------------------------------------------*/
45 * \param p The position of this control point. It will be assigned to the input
46 * and output directions too.
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)
54} // curve::control_point::control_point()
56/*----------------------------------------------------------------------------*/
59 * \param p The position of this control point.
60 * \param input_direction The point in the direction of which the curve enters
62 * \param output_direction The point in the direction of which the curve leaves
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)
73} // curve::control_point::control_point()
75/*----------------------------------------------------------------------------*/
77 * \brief Get the position of this control point.
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
84} // curve::control_point::get_position()
86/*----------------------------------------------------------------------------*/
88 * \brief Get the point in the direction of which the curve enters this
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
95 return m_input_direction;
96} // curve::control_point::get_input_direction()
98/*----------------------------------------------------------------------------*/
100 * \brief Get the point in the direction of which the curve leaves this
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
107 return m_output_direction;
108} // curve::control_point::get_output_direction()
114/*----------------------------------------------------------------------------*/
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.
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)
127} // curve::section::resolved_point::resolved_point()
129/*----------------------------------------------------------------------------*/
131 * \brief Get The position of the point.
133template<typename C, typename Traits>
135typename claw::math::curve<C, Traits>::section::resolved_point::coordinate_type&
136claw::math::curve<C, Traits>::section::resolved_point::get_position() const
139} // curve::section::resolved_point::get_position()
141/*----------------------------------------------------------------------------*/
143 * \brief Get the section on which the point has been found.
145template<typename C, typename Traits>
146const typename claw::math::curve<C, Traits>::section&
147claw::math::curve<C, Traits>::section::resolved_point::get_section() const
150} // curve::section::::resolved_point::get_section()
152/*----------------------------------------------------------------------------*/
154 * \brief Get the date of the point on the section.
156template<typename C, typename Traits>
158claw::math::curve<C, Traits>::section::resolved_point::get_date() const
161} // curve::section::resolved_point::get_date()
168/*----------------------------------------------------------------------------*/
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.
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)
180} // curve::section::section()
182/*----------------------------------------------------------------------------*/
184 * \brief Get the point of this section at a given date.
185 * \param t The date at which the point is computed.
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
191 if ( m_origin == m_end )
192 return m_origin->get_position();
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()) ) );
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()) ) );
207 return traits_type::make_coordinate( x, y );
208} // curve::section::get_point_at()
210/*----------------------------------------------------------------------------*/
212 * \brief Get the direction of the tangent at a given date.
213 * \param t The date for which we want the tangent.
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
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()) );
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()) );
232 return traits_type::make_coordinate( 0, (dy < 0) ? -1 : 1 );
234 return traits_type::make_coordinate( 1, dy / dx );
235} // curve::section::get_tangent_at()
237/*----------------------------------------------------------------------------*/
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
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
249 std::vector<resolved_point> result;
254 const std::vector<double> 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() ) ) );
261 for ( std::size_t i=0; i!=roots.size(); ++i )
263 ( resolved_point( get_point_at( roots[i] ), *this, roots[i] ) );
265 ensure_ends_in_points
267 (x == m_origin->get_position().x), (x == m_end->get_position().x) );
270 return extract_domain_points( result );
273} // curve::section::get_point_at_x()
275/*----------------------------------------------------------------------------*/
277 * \brief Get an iterator on the control point at the origin of the section in
278 * the curve from which it was created.
280template<typename C, typename Traits>
281const typename claw::math::curve<C, Traits>::section::iterator_type&
282claw::math::curve<C, Traits>::section::get_origin() const
285} // curve::section::get_origin()
287/*----------------------------------------------------------------------------*/
289 * \brief Tell if there is no points on this section.
291template<typename C, typename Traits>
292bool claw::math::curve<C, Traits>::section::empty() const
294 return m_origin == m_end;
295} // curve::section::empty()
297/*----------------------------------------------------------------------------*/
299 * \brief Get the value of the curve's equation on one dimension, at a given
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.
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
317 const double dt(1 - t);
319 return origin * dt * dt * dt
320 + 3 * output_direction * t * dt * dt
321 + 3 * input_direction * t * t * dt
323} // curve::section::evaluate()
325/*----------------------------------------------------------------------------*/
327 * \brief Get the value at a given date of the curve's derived equation on one
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.
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
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)
349} // curve::section::evaluate_derived()
351/*----------------------------------------------------------------------------*/
353 * \brief Ensure that a vector of resolved_point contains some ends of the
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
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.
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
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());
374 for ( std::size_t i=0; i!=p.size(); ++i )
376 const double distance_origin( std::abs( p[i].get_date() ) );
378 if ( distance_origin <= min_distance_origin )
380 min_distance_origin = distance_origin;
384 const double distance_end( std::abs( 1 - p[i].get_date() ) );
386 if ( distance_end <= min_distance_end )
388 min_distance_end = distance_end;
394 p[origin_index] = resolved_point( m_origin->get_position(), *this, 0.0 );
397 p[end_index] = resolved_point( m_end->get_position(), *this, 1.0 );
398} // curve::section::ensure_ends_in_points()
400/*----------------------------------------------------------------------------*/
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].
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
410 std::vector<resolved_point> clean_result;
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] );
417} // curve::section::extract_domain_points()
419/*----------------------------------------------------------------------------*/
421 * \brief Get the dates at which the curve passes at a given coordinate, on a
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.
433template<typename C, typename Traits>
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
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 );
446 return get_roots_degree_2(b, c, d);
448 return get_roots_degree_3(a, b, c, d);
449} // curve::section::get_roots()
451/*----------------------------------------------------------------------------*/
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.
459template<typename C, typename Traits>
461claw::math::curve<C, Traits>::section::get_roots_degree_2
462( value_type a, value_type b, value_type c ) const
464 const value_type delta( b * b - 4 * a * c );
466 std::vector<double> result;
469 result.push_back( - b / ( 2 * a ) );
470 else if ( delta > 0 )
472 result.push_back( (-b - std::sqrt(delta)) / (2 * a) );
473 result.push_back( (-b + std::sqrt(delta)) / (2 * a) );
477} // curve::section::get_roots_degree_2()
479/*----------------------------------------------------------------------------*/
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.
488template<typename C, typename Traits>
490claw::math::curve<C, Traits>::section::get_roots_degree_3
491( value_type a, value_type b, value_type c, value_type d ) const
493 // The following is the application of the method of Cardan
495 const value_type p( -(b * b) / (3.0 * a * a) + c / a );
498 * ( (2.0 * b * b) / (a * a)
502 const value_type delta( q * q + 4.0 * p * p * p / 27.0 );
504 std::vector<double> result;
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) );
516 else if ( delta > 0 )
520 ( (-q + std::sqrt(delta)) / 2.0 )
522 ( (-q - std::sqrt(delta)) / 2.0 ) - b / (3.0 * a));
525 for ( std::size_t i=0; i!=3; ++i )
527 ( 2.0 * std::sqrt( -p / 3.0 )
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 )
534} // curve::section::get_roots_degree_3()
541/*----------------------------------------------------------------------------*/
543 * \brief Add a point at the end of the curve.
544 * \param p The point to add.
546template<typename C, typename Traits>
547void claw::math::curve<C, Traits>::push_back( const control_point& p )
549 m_points.push_back(p);
550} // curve::push_back()
552/*----------------------------------------------------------------------------*/
554 * \brief Add a point at the beginning of the curve.
555 * \param p The point to add.
557template<typename C, typename Traits>
558void claw::math::curve<C, Traits>::push_front( const control_point& p )
560 m_points.push_front(p);
561} // curve::push_front()
563/*----------------------------------------------------------------------------*/
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.
569template<typename C, typename Traits>
570void claw::math::curve<C, Traits>::insert
571( const iterator& pos, const control_point& p )
573 m_points.insert( pos, p );
576/*----------------------------------------------------------------------------*/
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
582template<typename C, typename Traits>
583typename claw::math::curve<C, Traits>::section
584claw::math::curve<C, Traits>::get_section( const const_iterator& pos ) const
586 const_iterator it(pos);
590 return section( pos, pos );
592 return section( pos, it );
593} // curve::get_section()
595/*----------------------------------------------------------------------------*/
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
602 * \todo Remove the duplicates in the result.
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
609 typedef std::vector<typename section::resolved_point> result_type;
612 for ( const_iterator it=begin(); it!=end(); ++it )
614 const section s( get_section(it) );
618 const result_type new_points( s.get_point_at_x(x) );
619 result.insert( result.end(), new_points.begin(), new_points.end() );
625 const result_type before( get_point_at_x_before_origin(x) );
626 result.insert( result.begin(), before.begin(), before.end() );
628 const result_type after( get_point_at_x_after_end(x) );
629 result.insert( result.end(), after.begin(), after.end() );
633} // curve::get_point_at_x()
635/*----------------------------------------------------------------------------*/
637 * \brief Get an iterator on the first control point.
639template<typename C, typename Traits>
640typename claw::math::curve<C, Traits>::iterator
641claw::math::curve<C, Traits>::begin()
643 return m_points.begin();
646/*----------------------------------------------------------------------------*/
648 * \brief Get an iterator past the last control point.
650template<typename C, typename Traits>
651typename claw::math::curve<C, Traits>::iterator
652claw::math::curve<C, Traits>::end()
654 return m_points.end();
657/*----------------------------------------------------------------------------*/
659 * \brief Get an iterator on the first control point.
661template<typename C, typename Traits>
662typename claw::math::curve<C, Traits>::const_iterator
663claw::math::curve<C, Traits>::begin() const
665 return m_points.begin();
668/*----------------------------------------------------------------------------*/
670 * \brief Get an iterator past the last control point.
672template<typename C, typename Traits>
673typename claw::math::curve<C, Traits>::const_iterator
674claw::math::curve<C, Traits>::end() const
676 return m_points.end();
679/*----------------------------------------------------------------------------*/
681 * \brief Get the points having the given x-coordinate before the origin of the
683 * \param x The coordinate for which we want the points.
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
689 typedef std::vector<typename section::resolved_point> result_type;
692 const section s( get_section(begin()) );
696 const result_type points( s.get_point_at_x(x, true) );
698 for ( std::size_t i(0); i!=points.size(); ++i )
699 if ( points[i].get_date() < 0 )
700 result.push_back( points[i] );
704} // curve::get_point_at_x_before_origin()
706/*----------------------------------------------------------------------------*/
708 * \brief Get the points having the given x-coordinate after the end of the
710 * \param x The coordinate for which we want the points.
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
716 typedef std::vector<typename section::resolved_point> result_type;
719 if ( m_points.size() < 2 )
722 const_iterator it(end());
723 std::advance(it, -2);
725 const section s( get_section( it ) );
729 const result_type points( s.get_point_at_x(x, true) );
731 for ( std::size_t i(0); i!=points.size(); ++i )
732 if ( points[i].get_date() > 1 )
733 result.push_back( points[i] );
737} // curve::get_point_at_x_after_end()