Claw  1.7.3
impl/curve.tpp
Go to the documentation of this file.
00001 /*
00002   CLAW - a C++ Library Absolutely Wonderful
00003 
00004   CLAW is a free library without any particular aim but being useful to 
00005   anyone.
00006 
00007   Copyright (C) 2005-2011 Julien Jorge
00008 
00009   This library is free software; you can redistribute it and/or
00010   modify it under the terms of the GNU Lesser General Public
00011   License as published by the Free Software Foundation; either
00012   version 2.1 of the License, or (at your option) any later version.
00013 
00014   This library is distributed in the hope that it will be useful,
00015   but WITHOUT ANY WARRANTY; without even the implied warranty of
00016   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017   Lesser General Public License for more details.
00018 
00019   You should have received a copy of the GNU Lesser General Public
00020   License along with this library; if not, write to the Free Software
00021   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00022   contact: julien.jorge@gamned.org
00023 */
00029 #include <boost/math/special_functions/cbrt.hpp>
00030 #include <boost/math/constants/constants.hpp>
00031 
00032 /*----------------------------------------------------------------------------*/
00036 template<typename C, typename Traits>
00037 claw::math::curve<C, Traits>::control_point::control_point()
00038 {
00039 
00040 } // curve::control_point::control_point()
00041 
00042 /*----------------------------------------------------------------------------*/
00048 template<typename C, typename Traits>
00049 claw::math::curve<C, Traits>::control_point::control_point
00050 ( const coordinate_type& p )
00051   : m_position(p), m_input_direction(p), m_output_direction(p)
00052 {
00053 
00054 } // curve::control_point::control_point()
00055 
00056 /*----------------------------------------------------------------------------*/
00065 template<typename C, typename Traits>
00066 claw::math::curve<C, Traits>::control_point::control_point
00067 ( const coordinate_type& p, const coordinate_type& input_direction,
00068   const coordinate_type& output_direction )
00069   : m_position(p), m_input_direction(input_direction),
00070     m_output_direction(output_direction)
00071 {
00072 
00073 } // curve::control_point::control_point()
00074 
00075 /*----------------------------------------------------------------------------*/
00079 template<typename C, typename Traits>
00080 const typename claw::math::curve<C, Traits>::control_point::coordinate_type&
00081 claw::math::curve<C, Traits>::control_point::get_position() const
00082 {
00083   return m_position;
00084 } // curve::control_point::get_position()
00085 
00086 /*----------------------------------------------------------------------------*/
00091 template<typename C, typename Traits>
00092 const typename claw::math::curve<C, Traits>::control_point::coordinate_type&
00093 claw::math::curve<C, Traits>::control_point::get_input_direction() const
00094 {
00095   return m_input_direction;
00096 } // curve::control_point::get_input_direction()
00097 
00098 /*----------------------------------------------------------------------------*/
00103 template<typename C, typename Traits>
00104 const typename claw::math::curve<C, Traits>::control_point::coordinate_type&
00105 claw::math::curve<C, Traits>::control_point::get_output_direction() const
00106 {
00107   return m_output_direction;
00108 } // curve::control_point::get_output_direction()
00109 
00110 
00111 
00112 
00113 
00114 /*----------------------------------------------------------------------------*/
00121 template<typename C, typename Traits>
00122 claw::math::curve<C, Traits>::section::resolved_point::resolved_point
00123 ( const coordinate_type& position, const section& s, const double t )
00124   : m_position(position), m_section(s), m_date(t)
00125 {
00126 
00127 } // curve::section::resolved_point::resolved_point()
00128 
00129 /*----------------------------------------------------------------------------*/
00133 template<typename C, typename Traits>
00134 const
00135 typename claw::math::curve<C, Traits>::section::resolved_point::coordinate_type&
00136 claw::math::curve<C, Traits>::section::resolved_point::get_position() const
00137 {
00138   return m_position;
00139 } // curve::section::resolved_point::get_position()
00140 
00141 /*----------------------------------------------------------------------------*/
00145 template<typename C, typename Traits>
00146 const typename claw::math::curve<C, Traits>::section&
00147 claw::math::curve<C, Traits>::section::resolved_point::get_section() const
00148 {
00149   return m_section;
00150 } // curve::section::::resolved_point::get_section()
00151 
00152 /*----------------------------------------------------------------------------*/
00156 template<typename C, typename Traits>
00157 double
00158 claw::math::curve<C, Traits>::section::resolved_point::get_date() const
00159 {
00160   return m_date;
00161 } // curve::section::resolved_point::get_date()
00162 
00163 
00164 
00165 
00166 
00167 
00168 /*----------------------------------------------------------------------------*/
00174 template<typename C, typename Traits>
00175 claw::math::curve<C, Traits>::section::section
00176 ( const iterator_type& origin, const iterator_type& end )
00177   : m_origin(origin), m_end(end)
00178 {
00179 
00180 } // curve::section::section()
00181 
00182 /*----------------------------------------------------------------------------*/
00187 template<typename C, typename Traits>
00188 typename claw::math::curve<C, Traits>::section::coordinate_type
00189 claw::math::curve<C, Traits>::section::get_point_at( double t ) const
00190 {
00191   if ( m_origin == m_end )
00192     return m_origin->get_position();
00193 
00194   const value_type x
00195     ( evaluate
00196       ( t, traits_type::get_x(m_origin->get_position()),
00197         traits_type::get_x(m_origin->get_output_direction()),
00198         traits_type::get_x(m_end->get_input_direction()),
00199         traits_type::get_x(m_end->get_position()) ) );
00200   const value_type y
00201     ( evaluate
00202       ( t, traits_type::get_y(m_origin->get_position()),
00203         traits_type::get_y(m_origin->get_output_direction()),
00204         traits_type::get_y(m_end->get_input_direction()),
00205         traits_type::get_y(m_end->get_position()) ) );
00206 
00207   return traits_type::make_coordinate( x, y );
00208 } // curve::section::get_point_at()
00209 
00210 /*----------------------------------------------------------------------------*/
00215 template<typename C, typename Traits>
00216 typename claw::math::curve<C, Traits>::section::coordinate_type
00217 claw::math::curve<C, Traits>::section::get_tangent_at( double t ) const
00218 {
00219   const value_type dx = evaluate_derived
00220     ( t, traits_type::get_x(m_origin->get_position()),
00221       traits_type::get_x(m_origin->get_output_direction()),
00222       traits_type::get_x(m_end->get_input_direction()),
00223       traits_type::get_x(m_end->get_position()) );
00224 
00225   const value_type dy = evaluate_derived
00226     ( t, traits_type::get_y(m_origin->get_position()),
00227       traits_type::get_y(m_origin->get_output_direction()),
00228       traits_type::get_y(m_end->get_input_direction()),
00229       traits_type::get_y(m_end->get_position()) );
00230 
00231   if ( dx == 0 )
00232     return traits_type::make_coordinate( 0, (dy < 0) ? -1 : 1 );
00233   else
00234     return traits_type::make_coordinate( 1, dy / dx );
00235 } // curve::section::get_tangent_at()
00236 
00237 /*----------------------------------------------------------------------------*/
00244 template<typename C, typename Traits>
00245 std::vector<typename claw::math::curve<C, Traits>::section::resolved_point>
00246 claw::math::curve<C, Traits>::section::get_point_at_x
00247 ( value_type x, bool off_domain ) const
00248 {
00249   std::vector<resolved_point> result;
00250 
00251   if ( empty() )
00252     return result;
00253   
00254   const std::vector<double> roots
00255     ( get_roots
00256       ( x, traits_type::get_x(m_origin->get_position()),
00257         traits_type::get_x(m_origin->get_output_direction()),
00258         traits_type::get_x(m_end->get_input_direction()),
00259         traits_type::get_x(m_end->get_position() ) ) );
00260   
00261   for ( std::size_t i=0; i!=roots.size(); ++i )
00262     result.push_back
00263       ( resolved_point( get_point_at( roots[i] ), *this, roots[i] ) );
00264 
00265   ensure_ends_in_points
00266     ( result,
00267       (x == m_origin->get_position().x), (x == m_end->get_position().x) );
00268 
00269   if ( !off_domain )
00270     return extract_domain_points( result );
00271   else
00272     return result;
00273 } // curve::section::get_point_at_x()
00274 
00275 /*----------------------------------------------------------------------------*/
00280 template<typename C, typename Traits>
00281 const typename claw::math::curve<C, Traits>::section::iterator_type&
00282 claw::math::curve<C, Traits>::section::get_origin() const
00283 {
00284   return m_origin;
00285 } // curve::section::get_origin()
00286 
00287 /*----------------------------------------------------------------------------*/
00291 template<typename C, typename Traits>
00292 bool claw::math::curve<C, Traits>::section::empty() const
00293 {
00294   return m_origin == m_end;
00295 } // curve::section::empty()
00296 
00297 /*----------------------------------------------------------------------------*/
00311 template<typename C, typename Traits>
00312 typename claw::math::curve<C, Traits>::section::value_type
00313 claw::math::curve<C, Traits>::section::evaluate
00314 ( double t, value_type origin, value_type output_direction,
00315   value_type input_direction, value_type end ) const
00316 {
00317   const double dt(1 - t);
00318  
00319   return origin * dt * dt * dt
00320     + 3 * output_direction * t * dt * dt
00321     + 3 * input_direction * t * t * dt
00322     + end * t * t * t;
00323 } // curve::section::evaluate()
00324 
00325 /*----------------------------------------------------------------------------*/
00339 template<typename C, typename Traits>
00340 typename claw::math::curve<C, Traits>::section::value_type
00341 claw::math::curve<C, Traits>::section::evaluate_derived
00342 ( double t, value_type origin, value_type output_direction,
00343   value_type input_direction, value_type end ) const
00344 {
00345   return - 3 * origin + 3 * output_direction
00346     + (6 * origin - 12 * output_direction + 6 * input_direction) * t
00347     + (-3 * origin + 9 * output_direction - 9 * input_direction + 3 * end)
00348     * t * t;
00349 } // curve::section::evaluate_derived()
00350 
00351 /*----------------------------------------------------------------------------*/
00365 template<typename C, typename Traits>
00366 void claw::math::curve<C, Traits>::section::ensure_ends_in_points
00367 ( std::vector<resolved_point>& p, bool ensure_origin, bool ensure_end ) const
00368 {
00369   double min_distance_origin( std::numeric_limits<double>::max() );
00370   double min_distance_end( std::numeric_limits<double>::max() );
00371   std::size_t origin_index(p.size());
00372   std::size_t end_index(p.size());
00373 
00374   for ( std::size_t i=0; i!=p.size(); ++i )
00375     {
00376       const double distance_origin( std::abs( p[i].get_date() ) );
00377 
00378       if ( distance_origin <= min_distance_origin )
00379         {
00380           min_distance_origin = distance_origin;
00381           origin_index = i;
00382         }
00383 
00384       const double distance_end( std::abs( 1 - p[i].get_date() ) );
00385 
00386       if ( distance_end <= min_distance_end )
00387         {
00388           min_distance_end = distance_end;
00389           end_index = i;
00390         }
00391     }
00392 
00393   if ( ensure_origin )
00394     p[origin_index] = resolved_point( m_origin->get_position(), *this, 0.0 );
00395 
00396   if ( ensure_end )
00397     p[end_index] = resolved_point( m_end->get_position(), *this, 1.0 );
00398 } // curve::section::ensure_ends_in_points()
00399 
00400 /*----------------------------------------------------------------------------*/
00405 template<typename C, typename Traits>
00406 std::vector<typename claw::math::curve<C, Traits>::section::resolved_point>
00407 claw::math::curve<C, Traits>::section::extract_domain_points
00408 ( const std::vector<resolved_point>& p ) const
00409 {
00410   std::vector<resolved_point> clean_result;
00411 
00412   for ( std::size_t i=0; i!=p.size(); ++i )
00413     if ( (p[i].get_date() >= 0) && (p[i].get_date() <= 1) )
00414       clean_result.push_back( p[i] );
00415 
00416   return clean_result;
00417 } // curve::section::extract_domain_points()
00418 
00419 /*----------------------------------------------------------------------------*/
00433 template<typename C, typename Traits>
00434 std::vector<double>
00435 claw::math::curve<C, Traits>::section::get_roots
00436 ( value_type x, value_type origin, value_type output_direction,
00437   value_type input_direction, value_type end ) const
00438 {
00439   const value_type a
00440     (-origin + 3 * output_direction - 3 * input_direction + end );
00441   const value_type b( 3 * origin - 6 * output_direction + 3 * input_direction );
00442   const value_type c( -3 * origin + 3 * output_direction );
00443   const value_type d( origin - x );
00444 
00445   if ( a == 0 )
00446     return get_roots_degree_2(b, c, d);
00447   else
00448     return get_roots_degree_3(a, b, c, d);
00449 } // curve::section::get_roots()
00450 
00451 /*----------------------------------------------------------------------------*/
00459 template<typename C, typename Traits>
00460 std::vector<double>
00461 claw::math::curve<C, Traits>::section::get_roots_degree_2
00462 ( value_type a, value_type b, value_type c ) const
00463 {
00464   const value_type delta( b * b - 4 * a * c );
00465 
00466   std::vector<double> result;
00467 
00468   if ( delta == 0 )
00469     result.push_back( - b / ( 2 * a ) );
00470   else if ( delta > 0 )
00471     {
00472       result.push_back( (-b - std::sqrt(delta)) / (2 * a) );
00473       result.push_back( (-b + std::sqrt(delta)) / (2 * a) );
00474     }
00475 
00476   return result;
00477 } // curve::section::get_roots_degree_2()
00478 
00479 /*----------------------------------------------------------------------------*/
00488 template<typename C, typename Traits>
00489 std::vector<double>
00490 claw::math::curve<C, Traits>::section::get_roots_degree_3
00491 ( value_type a, value_type b, value_type c, value_type d ) const
00492 {
00493   // The following is the application of the method of Cardan
00494 
00495   const value_type p( -(b * b) / (3.0 * a * a) + c / a );
00496   const value_type q
00497     ( ( b / (27.0 * a) )
00498       * ( (2.0 * b * b) / (a * a)
00499           - 9.0 * c / a )
00500       + d / a );
00501 
00502   const value_type delta( q * q + 4.0 * p * p * p / 27.0 );
00503 
00504   std::vector<double> result;
00505 
00506   if ( delta == 0 )
00507     {
00508       if ( p == 0 )
00509         result.push_back(0);
00510       else
00511         {
00512           result.push_back( 3.0 * q / p - b / (3.0 * a) );
00513           result.push_back( - 3.0 * q / (2.0 * p)  - b / (3.0 * a) );
00514         }
00515     }
00516   else if ( delta > 0 )
00517     {
00518     result.push_back
00519       ( boost::math::cbrt
00520         ( (-q + std::sqrt(delta)) / 2.0 )
00521         + boost::math::cbrt
00522         ( (-q - std::sqrt(delta)) / 2.0 ) - b / (3.0 * a));
00523     }
00524   else
00525     for ( std::size_t i=0; i!=3; ++i )
00526       result.push_back
00527         ( 2.0 * std::sqrt( -p / 3.0 )
00528           * std::cos
00529           ( std::acos( std::sqrt(27.0 / (- p * p * p)) * - q / 2.0 ) / 3.0
00530             + 2.0 * i * boost::math::constants::pi<double>() / 3.0 ) 
00531           - b / (3.0 * a));
00532 
00533   return result;
00534 } // curve::section::get_roots_degree_3()
00535 
00536 
00537 
00538 
00539 
00540 
00541 /*----------------------------------------------------------------------------*/
00546 template<typename C, typename Traits>
00547 void claw::math::curve<C, Traits>::push_back( const control_point& p )
00548 {
00549   m_points.push_back(p);
00550 } // curve::push_back()
00551 
00552 /*----------------------------------------------------------------------------*/
00557 template<typename C, typename Traits>
00558 void claw::math::curve<C, Traits>::push_front( const control_point& p )
00559 {
00560   m_points.push_front(p);
00561 } // curve::push_front()
00562 
00563 /*----------------------------------------------------------------------------*/
00569 template<typename C, typename Traits>
00570 void claw::math::curve<C, Traits>::insert
00571 ( const iterator& pos, const control_point& p )
00572 {
00573   m_points.insert( pos, p );
00574 } // curve::insert()
00575 
00576 /*----------------------------------------------------------------------------*/
00582 template<typename C, typename Traits>
00583 typename claw::math::curve<C, Traits>::section
00584 claw::math::curve<C, Traits>::get_section( const const_iterator& pos ) const
00585 {
00586   const_iterator it(pos);
00587   ++it;
00588 
00589   if ( it == end() )
00590     return section( pos, pos );
00591   else
00592     return section( pos, it );
00593 } // curve::get_section()
00594 
00595 /*----------------------------------------------------------------------------*/
00604 template<typename C, typename Traits>
00605 std::vector< typename claw::math::curve<C, Traits>::section::resolved_point >
00606 claw::math::curve<C, Traits>::get_point_at_x
00607 ( value_type x, bool off_domain ) const
00608 {
00609   typedef std::vector<typename section::resolved_point> result_type;
00610   result_type result;
00611   
00612   for ( const_iterator it=begin(); it!=end(); ++it )
00613     {
00614       const section s( get_section(it) );
00615 
00616       if ( !s.empty() )
00617         {
00618           const result_type new_points( s.get_point_at_x(x) );
00619           result.insert( result.end(), new_points.begin(), new_points.end() );
00620         }
00621     }
00622 
00623   if ( off_domain )
00624     {
00625       const result_type before( get_point_at_x_before_origin(x) );
00626       result.insert( result.begin(), before.begin(), before.end() );
00627 
00628       const result_type after( get_point_at_x_after_end(x) );
00629       result.insert( result.end(), after.begin(), after.end() );
00630     }
00631 
00632   return result;
00633 } // curve::get_point_at_x()
00634 
00635 /*----------------------------------------------------------------------------*/
00639 template<typename C, typename Traits>
00640 typename claw::math::curve<C, Traits>::iterator
00641 claw::math::curve<C, Traits>::begin()
00642 {
00643   return m_points.begin();
00644 } // curve::begin()
00645 
00646 /*----------------------------------------------------------------------------*/
00650 template<typename C, typename Traits>
00651 typename claw::math::curve<C, Traits>::iterator
00652 claw::math::curve<C, Traits>::end()
00653 {
00654   return m_points.end();
00655 } // curve::end()
00656 
00657 /*----------------------------------------------------------------------------*/
00661 template<typename C, typename Traits>
00662 typename claw::math::curve<C, Traits>::const_iterator
00663 claw::math::curve<C, Traits>::begin() const
00664 {
00665   return m_points.begin();
00666 } // curve::begin()
00667 
00668 /*----------------------------------------------------------------------------*/
00672 template<typename C, typename Traits>
00673 typename claw::math::curve<C, Traits>::const_iterator
00674 claw::math::curve<C, Traits>::end() const
00675 {
00676   return m_points.end();
00677 } // curve::end()
00678 
00679 /*----------------------------------------------------------------------------*/
00685 template<typename C, typename Traits>
00686 std::vector< typename claw::math::curve<C, Traits>::section::resolved_point >
00687 claw::math::curve<C, Traits>::get_point_at_x_before_origin( value_type x ) const
00688 {
00689   typedef std::vector<typename section::resolved_point> result_type;
00690   result_type result;
00691   
00692   const section s( get_section(begin()) );
00693 
00694   if ( !s.empty() )
00695     {
00696       const result_type points( s.get_point_at_x(x, true) );
00697 
00698       for ( std::size_t i(0); i!=points.size(); ++i )
00699         if ( points[i].get_date() < 0 )
00700           result.push_back( points[i] );
00701     }
00702 
00703   return result;
00704 } // curve::get_point_at_x_before_origin()
00705 
00706 /*----------------------------------------------------------------------------*/
00712 template<typename C, typename Traits>
00713 std::vector< typename claw::math::curve<C, Traits>::section::resolved_point >
00714 claw::math::curve<C, Traits>::get_point_at_x_after_end( value_type x ) const
00715 {
00716   typedef std::vector<typename section::resolved_point> result_type;
00717   result_type result;
00718   
00719   if ( m_points.size() < 2 )
00720     return result;
00721 
00722   const_iterator it(end());
00723   std::advance(it, -2);
00724 
00725   const section s( get_section( it ) );
00726 
00727   if ( !s.empty() )
00728     {
00729       const result_type points( s.get_point_at_x(x, true) );
00730 
00731       for ( std::size_t i(0); i!=points.size(); ++i )
00732         if ( points[i].get_date() > 1 )
00733           result.push_back( points[i] );
00734     }
00735 
00736   return result;
00737 } // curve::get_point_at_x_after_end()