curve.tpp

Go to the documentation of this file.

/*
  CLAW - a C++ Library Absolutely Wonderful

  CLAW is a free library without any particular aim but being useful to 
  anyone.

  Copyright (C) 2005-2011 Julien Jorge

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  contact: julien.jorge@gamned.org
*/
#include <boost/math/special_functions/cbrt.hpp>
#include <boost/math/constants/constants.hpp>

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
claw::math::curve<C, Traits>::control_point::control_point()
{

} // curve::control_point::control_point()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
claw::math::curve<C, Traits>::control_point::control_point
( const coordinate_type& p )
  : m_position(p), m_input_direction(p), m_output_direction(p)
{

} // curve::control_point::control_point()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
claw::math::curve<C, Traits>::control_point::control_point
( const coordinate_type& p, const coordinate_type& input_direction,
  const coordinate_type& output_direction )
  : m_position(p), m_input_direction(input_direction),
    m_output_direction(output_direction)
{

} // curve::control_point::control_point()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
const typename claw::math::curve<C, Traits>::control_point::coordinate_type&
claw::math::curve<C, Traits>::control_point::get_position() const
{
  return m_position;
} // curve::control_point::get_position()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
const typename claw::math::curve<C, Traits>::control_point::coordinate_type&
claw::math::curve<C, Traits>::control_point::get_input_direction() const
{
  return m_input_direction;
} // curve::control_point::get_input_direction()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
const typename claw::math::curve<C, Traits>::control_point::coordinate_type&
claw::math::curve<C, Traits>::control_point::get_output_direction() const
{
  return m_output_direction;
} // curve::control_point::get_output_direction()





/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
claw::math::curve<C, Traits>::section::resolved_point::resolved_point
( const coordinate_type& position, const section& s, const double t )
  : m_position(position), m_section(s), m_date(t)
{

} // curve::section::resolved_point::resolved_point()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
const
typename claw::math::curve<C, Traits>::section::resolved_point::coordinate_type&
claw::math::curve<C, Traits>::section::resolved_point::get_position() const
{
  return m_position;
} // curve::section::resolved_point::get_position()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
const typename claw::math::curve<C, Traits>::section&
claw::math::curve<C, Traits>::section::resolved_point::get_section() const
{
  return m_section;
} // curve::section::::resolved_point::get_section()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
double
claw::math::curve<C, Traits>::section::resolved_point::get_date() const
{
  return m_date;
} // curve::section::resolved_point::get_date()






/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
claw::math::curve<C, Traits>::section::section
( const iterator_type& origin, const iterator_type& end )
  : m_origin(origin), m_end(end)
{

} // curve::section::section()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
typename claw::math::curve<C, Traits>::section::coordinate_type
claw::math::curve<C, Traits>::section::get_point_at( double t ) const
{
  if ( m_origin == m_end )
    return m_origin->get_position();

  const value_type x
    ( evaluate
      ( t, traits_type::get_x(m_origin->get_position()),
        traits_type::get_x(m_origin->get_output_direction()),
        traits_type::get_x(m_end->get_input_direction()),
        traits_type::get_x(m_end->get_position()) ) );
  const value_type y
    ( evaluate
      ( t, traits_type::get_y(m_origin->get_position()),
        traits_type::get_y(m_origin->get_output_direction()),
        traits_type::get_y(m_end->get_input_direction()),
        traits_type::get_y(m_end->get_position()) ) );

  return traits_type::make_coordinate( x, y );
} // curve::section::get_point_at()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
typename claw::math::curve<C, Traits>::section::coordinate_type
claw::math::curve<C, Traits>::section::get_tangent_at( double t ) const
{
  const value_type dx = evaluate_derived
    ( t, traits_type::get_x(m_origin->get_position()),
      traits_type::get_x(m_origin->get_output_direction()),
      traits_type::get_x(m_end->get_input_direction()),
      traits_type::get_x(m_end->get_position()) );

  const value_type dy = evaluate_derived
    ( t, traits_type::get_y(m_origin->get_position()),
      traits_type::get_y(m_origin->get_output_direction()),
      traits_type::get_y(m_end->get_input_direction()),
      traits_type::get_y(m_end->get_position()) );

  if ( dx == 0 )
    return traits_type::make_coordinate( 0, (dy < 0) ? -1 : 1 );
  else
    return traits_type::make_coordinate( 1, dy / dx );
} // curve::section::get_tangent_at()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
std::vector<typename claw::math::curve<C, Traits>::section::resolved_point>
claw::math::curve<C, Traits>::section::get_point_at_x
( value_type x, bool off_domain ) const
{
  std::vector<resolved_point> result;

  if ( empty() )
    return result;
  
  const std::vector<double> roots
    ( get_roots
      ( x, traits_type::get_x(m_origin->get_position()),
        traits_type::get_x(m_origin->get_output_direction()),
        traits_type::get_x(m_end->get_input_direction()),
        traits_type::get_x(m_end->get_position() ) ) );
  
  for ( std::size_t i=0; i!=roots.size(); ++i )
    result.push_back
      ( resolved_point( get_point_at( roots[i] ), *this, roots[i] ) );

  ensure_ends_in_points
    ( result,
      (x == m_origin->get_position().x), (x == m_end->get_position().x) );

  if ( !off_domain )
    return extract_domain_points( result );
  else
    return result;
} // curve::section::get_point_at_x()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
const typename claw::math::curve<C, Traits>::section::iterator_type&
claw::math::curve<C, Traits>::section::get_origin() const
{
  return m_origin;
} // curve::section::get_origin()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
bool claw::math::curve<C, Traits>::section::empty() const
{
  return m_origin == m_end;
} // curve::section::empty()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
typename claw::math::curve<C, Traits>::section::value_type
claw::math::curve<C, Traits>::section::evaluate
( double t, value_type origin, value_type output_direction,
  value_type input_direction, value_type end ) const
{
  const double dt(1 - t);
 
  return origin * dt * dt * dt
    + 3 * output_direction * t * dt * dt
    + 3 * input_direction * t * t * dt
    + end * t * t * t;
} // curve::section::evaluate()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
typename claw::math::curve<C, Traits>::section::value_type
claw::math::curve<C, Traits>::section::evaluate_derived
( double t, value_type origin, value_type output_direction,
  value_type input_direction, value_type end ) const
{
  return - 3 * origin + 3 * output_direction
    + (6 * origin - 12 * output_direction + 6 * input_direction) * t
    + (-3 * origin + 9 * output_direction - 9 * input_direction + 3 * end)
    * t * t;
} // curve::section::evaluate_derived()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
void claw::math::curve<C, Traits>::section::ensure_ends_in_points
( std::vector<resolved_point>& p, bool ensure_origin, bool ensure_end ) const
{
  double min_distance_origin( std::numeric_limits<double>::max() );
  double min_distance_end( std::numeric_limits<double>::max() );
  std::size_t origin_index(p.size());
  std::size_t end_index(p.size());

  for ( std::size_t i=0; i!=p.size(); ++i )
    {
      const double distance_origin( std::abs( p[i].get_date() ) );

      if ( distance_origin <= min_distance_origin )
        {
          min_distance_origin = distance_origin;
          origin_index = i;
        }

      const double distance_end( std::abs( 1 - p[i].get_date() ) );

      if ( distance_end <= min_distance_end )
        {
          min_distance_end = distance_end;
          end_index = i;
        }
    }

  if ( ensure_origin )
    p[origin_index] = resolved_point( m_origin->get_position(), *this, 0.0 );

  if ( ensure_end )
    p[end_index] = resolved_point( m_end->get_position(), *this, 1.0 );
} // curve::section::ensure_ends_in_points()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
std::vector<typename claw::math::curve<C, Traits>::section::resolved_point>
claw::math::curve<C, Traits>::section::extract_domain_points
( const std::vector<resolved_point>& p ) const
{
  std::vector<resolved_point> clean_result;

  for ( std::size_t i=0; i!=p.size(); ++i )
    if ( (p[i].get_date() >= 0) && (p[i].get_date() <= 1) )
      clean_result.push_back( p[i] );

  return clean_result;
} // curve::section::extract_domain_points()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
std::vector<double>
claw::math::curve<C, Traits>::section::get_roots
( value_type x, value_type origin, value_type output_direction,
  value_type input_direction, value_type end ) const
{
  const value_type a
    (-origin + 3 * output_direction - 3 * input_direction + end );
  const value_type b( 3 * origin - 6 * output_direction + 3 * input_direction );
  const value_type c( -3 * origin + 3 * output_direction );
  const value_type d( origin - x );

  if ( a == 0 )
    return get_roots_degree_2(b, c, d);
  else
    return get_roots_degree_3(a, b, c, d);
} // curve::section::get_roots()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
std::vector<double>
claw::math::curve<C, Traits>::section::get_roots_degree_2
( value_type a, value_type b, value_type c ) const
{
  const value_type delta( b * b - 4 * a * c );

  std::vector<double> result;

  if ( delta == 0 )
    result.push_back( - b / ( 2 * a ) );
  else if ( delta > 0 )
    {
      result.push_back( (-b - std::sqrt(delta)) / (2 * a) );
      result.push_back( (-b + std::sqrt(delta)) / (2 * a) );
    }

  return result;
} // curve::section::get_roots_degree_2()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
std::vector<double>
claw::math::curve<C, Traits>::section::get_roots_degree_3
( value_type a, value_type b, value_type c, value_type d ) const
{
  // The following is the application of the method of Cardan

  const value_type p( -(b * b) / (3.0 * a * a) + c / a );
  const value_type q
    ( ( b / (27.0 * a) )
      * ( (2.0 * b * b) / (a * a)
          - 9.0 * c / a )
      + d / a );

  const value_type delta( q * q + 4.0 * p * p * p / 27.0 );

  std::vector<double> result;

  if ( delta == 0 )
    {
      if ( p == 0 )
        result.push_back(0);
      else
        {
          result.push_back( 3.0 * q / p - b / (3.0 * a) );
          result.push_back( - 3.0 * q / (2.0 * p)  - b / (3.0 * a) );
        }
    }
  else if ( delta > 0 )
    {
    result.push_back
      ( boost::math::cbrt
        ( (-q + std::sqrt(delta)) / 2.0 )
        + boost::math::cbrt
        ( (-q - std::sqrt(delta)) / 2.0 ) - b / (3.0 * a));
    }
  else
    for ( std::size_t i=0; i!=3; ++i )
      result.push_back
        ( 2.0 * std::sqrt( -p / 3.0 )
          * std::cos
          ( std::acos( std::sqrt(27.0 / (- p * p * p)) * - q / 2.0 ) / 3.0
            + 2.0 * i * boost::math::constants::pi<double>() / 3.0 ) 
          - b / (3.0 * a));

  return result;
} // curve::section::get_roots_degree_3()






/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
void claw::math::curve<C, Traits>::push_back( const control_point& p )
{
  m_points.push_back(p);
} // curve::push_back()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
void claw::math::curve<C, Traits>::push_front( const control_point& p )
{
  m_points.push_front(p);
} // curve::push_front()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
void claw::math::curve<C, Traits>::insert
( const iterator& pos, const control_point& p )
{
  m_points.insert( pos, p );
} // curve::insert()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
typename claw::math::curve<C, Traits>::section
claw::math::curve<C, Traits>::get_section( const const_iterator& pos ) const
{
  const_iterator it(pos);
  ++it;

  if ( it == end() )
    return section( pos, pos );
  else
    return section( pos, it );
} // curve::get_section()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
std::vector< typename claw::math::curve<C, Traits>::section::resolved_point >
claw::math::curve<C, Traits>::get_point_at_x
( value_type x, bool off_domain ) const
{
  typedef std::vector<typename section::resolved_point> result_type;
  result_type result;
  
  for ( const_iterator it=begin(); it!=end(); ++it )
    {
      const section s( get_section(it) );

      if ( !s.empty() )
        {
          const result_type new_points( s.get_point_at_x(x) );
          result.insert( result.end(), new_points.begin(), new_points.end() );
        }
    }

  if ( off_domain )
    {
      const result_type before( get_point_at_x_before_origin(x) );
      result.insert( result.begin(), before.begin(), before.end() );

      const result_type after( get_point_at_x_after_end(x) );
      result.insert( result.end(), after.begin(), after.end() );
    }

  return result;
} // curve::get_point_at_x()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
typename claw::math::curve<C, Traits>::iterator
claw::math::curve<C, Traits>::begin()
{
  return m_points.begin();
} // curve::begin()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
typename claw::math::curve<C, Traits>::iterator
claw::math::curve<C, Traits>::end()
{
  return m_points.end();
} // curve::end()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
typename claw::math::curve<C, Traits>::const_iterator
claw::math::curve<C, Traits>::begin() const
{
  return m_points.begin();
} // curve::begin()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
typename claw::math::curve<C, Traits>::const_iterator
claw::math::curve<C, Traits>::end() const
{
  return m_points.end();
} // curve::end()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
std::vector< typename claw::math::curve<C, Traits>::section::resolved_point >
claw::math::curve<C, Traits>::get_point_at_x_before_origin( value_type x ) const
{
  typedef std::vector<typename section::resolved_point> result_type;
  result_type result;
  
  const section s( get_section(begin()) );

  if ( !s.empty() )
    {
      const result_type points( s.get_point_at_x(x, true) );

      for ( std::size_t i(0); i!=points.size(); ++i )
        if ( points[i].get_date() < 0 )
          result.push_back( points[i] );
    }

  return result;
} // curve::get_point_at_x_before_origin()

/*----------------------------------------------------------------------------*/
template<typename C, typename Traits>
std::vector< typename claw::math::curve<C, Traits>::section::resolved_point >
claw::math::curve<C, Traits>::get_point_at_x_after_end( value_type x ) const
{
  typedef std::vector<typename section::resolved_point> result_type;
  result_type result;
  
  if ( m_points.size() < 2 )
    return result;

  const_iterator it(end());
  std::advance(it, -2);

  const section s( get_section( it ) );

  if ( !s.empty() )
    {
      const result_type points( s.get_point_at_x(x, true) );

      for ( std::size_t i(0); i!=points.size(); ++i )
        if ( points[i].get_date() > 1 )
          result.push_back( points[i] );
    }

  return result;
} // curve::get_point_at_x_after_end()