Claw
1.7.3
|
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()