source/bezier.h

   1 /* This file is part of Shapes.
   2  *
   3  * Shapes is free software: you can redistribute it and/or modify
   4  * it under the terms of the GNU General Public License as published by
   5  * the Free Software Foundation, either version 3 of the License, or
   6  * any later version.
   7  *
   8  * Shapes is distributed in the hope that it will be useful,
   9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11  * GNU General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU General Public License
  14  * along with Shapes.  If not, see <http://www.gnu.org/licenses/>.
  15  *
  16  * Copyright 2008, 2009 Henrik Tidefelt
  17  */
  18
  19 #pragma once
  20
  21 #include <cmath>
  22 #include <iostream>
  23
  24 namespace Bezier
  25 {
  26         template< class T >
  27                 class PolyCoeffs;
  28
  29         template< class T >
  30                 class ControlPoints
  31                 {
  32                 public:
  33                         T p0_;
  34                         T p1_;
  35                         T p2_;
  36                         T p3_;
  37
  38                         ControlPoints( const T & p0, const T & p1, const T & p2, const T & p3 );
  39                         ControlPoints( const PolyCoeffs< T > & polyCoeffs );
  40
  41                         T point( double t ) const;
  42                         void get( double * dst ) const;
  43                         void getEndpoints( double * dst ) const;
  44                 };
  45
  46         template< class T >
  47                 class PolyCoeffs
  48                 {
  49                 public:
  50                         T z0_;
  51                         T z1_;
  52                         T z2_;
  53                         T z3_;
  54
  55                         PolyCoeffs( const T & z0, const T & z1, const T & z2, const T & z3 );
  56                         PolyCoeffs( const ControlPoints< T > & controlPoints );
  57                         PolyCoeffs< T > subSection( double t0, double t1 ) const;
  58
  59                         T point( double t ) const;
  60                         T velocity( double t ) const;
  61                         T acceleration( double t ) const;
  62                         // The number of values returned in the destination argument by the following functions may vary.       A HUGE_VAL is used for termination.
  63                         void hyperplaneIntersections( double t[4], const T & n, const double d ) const; // Finds all t such that inner( point( t ), n ) == d.   Note that there can be at most 3 intersections.
  64                         void stationaryPoints( double t[3], const T & d ) const;        // Finds all t such that inner( velocity( t ), d ) == d.        Note that there can be at most 2 stationary points.
  65                         void inflections( double t[3] ) const;  // Finds all t such that acceleration is parallel to velocity.  Only defined for Coords2D, which makes sense since there are generally no such points in higher dimensions.     Note that there can be at most 2 points of inflection (this is not obvious at first glance)!
  66                 };
  67
  68         /* By "Bezier roots" I mean real roots in the interval [ 0, 1 ].
  69          */
  70         void bezierRootsOfPolynomial( double t[2], double k0, double k1 );
  71         void bezierRootsOfPolynomial( double t[3], double k0, double k1, double k2 );
  72         void bezierRootsOfPolynomial( double t[4], double k0, double k1, double k2, double k3 );
  73
  74         template< class T >
  75                 ControlPoints< T >::ControlPoints( const T & p0, const T & p1, const T & p2, const T & p3 )
  76                 : p0_( p0 ), p1_( p1 ), p2_( p2 ), p3_( p3 )
  77                 { }
  78
  79         template< class T >
  80                 ControlPoints< T >::ControlPoints( const PolyCoeffs< T > & pc )
  81                 : p0_( pc.z0_ ),
  82                         p1_( pc.z0_ + (1.0/3) * pc.z1_ ),
  83                         p2_( pc.z0_ + (1.0/3) * ( 2 * pc.z1_ + pc.z2_ ) ),
  84                         p3_( pc.z0_ + pc.z1_ + pc.z2_ + pc.z3_ )
  85                 { }
  86
  87         template< class T >
  88                 T
  89                 ControlPoints< T >::point( double t ) const
  90                 {
  91                         double ct = 1 - t;
  92                         double ct2 = ct * ct;
  93                         double t2 = t * t;
  94                         return ct2 * ct * p0_ + 3 * ct2 * t * p1_ + 3 * ct * t2 * p2_ + t2 * t * p3_;
  95                 }
  96
  97         template< class T >
  98                 std::ostream &
  99                 operator << ( std::ostream & os, const ControlPoints< T > & self )
 100                 {
 101                         os << "t -> (1-t)^{3} * " << self.p0_ << " + 3*(1-t)^{2}*t * " << self.p1_ << " + 3*(1-t)*t^{2} * " << self.p2_ << " + t^{3} * " << self.p3_ ;
 102                         return os;
 103                 }
 104
 105
 106         template< class T >
 107                 PolyCoeffs< T >::PolyCoeffs( const T & z0, const T & z1, const T & z2, const T & z3 )
 108                 : z0_( z0 ), z1_( z1 ), z2_( z2 ), z3_( z3 )
 109                 { }
 110
 111         template< class T >
 112                 PolyCoeffs< T >::PolyCoeffs( const ControlPoints< T > & cp )
 113                 : z0_( cp.p0_ ),
 114                         z1_( 3 * (cp.p1_ - cp.p0_) ),
 115                         z2_( 3 * ( cp.p0_ - 2 * cp.p1_ + cp.p2_) ),
 116                         z3_( 3 * cp.p1_ - 3 * cp.p2_ + cp.p3_ - cp.p0_ ) // put cp_.0 last since we don't rely un unary negation
 117                         { }
 118
 119         template< class T >
 120                 T
 121                 PolyCoeffs< T >::point( double t ) const
 122                 {
 123                         return ( ( z3_ * t + z2_ ) * t + z1_ ) * t + z0_;
 124                 }
 125
 126         template< class T >
 127                 T
 128                 PolyCoeffs< T >::velocity( double t ) const
 129                 {
 130                         return ( z3_ * ( ( 3./2.) * t ) + z2_ ) * ( 2 * t ) + z1_;
 131                 }
 132
 133         template< class T >
 134                 T
 135                 PolyCoeffs< T >::acceleration( double t ) const
 136                 {
 137                         return z3_ * ( 6 * t ) + z2_ * 2;
 138                 }
 139
 140         template< class T >
 141                 PolyCoeffs< T >
 142                 PolyCoeffs< T >::subSection( double t0, double t1 ) const
 143                 {
 144                         double dt = t1 - t0;
 145                         double dt2 = dt * dt;
 146                         return PolyCoeffs<T>( ( ( z3_ * t0 + z2_ ) * t0 + z1_ ) * t0 + z0_,
 147                                                                                                                 dt * ( ( 3 * t0 * z3_ + 2 * z2_ ) * t0 + z1_ ),
 148                                                                                                                 dt2 * ( z2_ + 3 * t0 * z3_ ),
 149                                                                                                                 dt2 * dt * z3_ );
 150                 }
 151
 152         template< class T >
 153                 void
 154                 PolyCoeffs< T >::hyperplaneIntersections( double t[4], const T & n, const double d ) const
 155                 {
 156                         double k3 = innerScalar( n, z3_ );
 157                         double k2 = innerScalar( n, z2_ );
 158                         double k1 = innerScalar( n, z1_ );
 159                         double k0 = innerScalar( n, z0_ ) - d;
 160
 161                         bezierRootsOfPolynomial( t, k0, k1, k2, k3 );
 162                 }
 163
 164         template< class T >
 165                 void
 166                 PolyCoeffs< T >::stationaryPoints( double t[3], const T & d ) const
 167                 {
 168                         double k0 = innerScalar( z1_, d );
 169                         double k1 = 2 * innerScalar( z2_, d );
 170                         double k2 = 3 * innerScalar( z3_, d );
 171
 172                         bezierRootsOfPolynomial( t, k0, k1, k2 );
 173                 }
 174
 175         template< class T >
 176         void
 177         PolyCoeffs< T >::inflections( double t[3] ) const
 178         {
 179                 /* We begin by constructing the polynomial coefficients of the scalar product betwenn the velocity and a normal to the acceleration.
 180                  */
 181
 182                 /* Note that we cannot just take any orthogonal vector for z2_ and z3_ separately; they must be constructed using the same transform.
 183                  * For generality, we avoid mentioning the components, although the method quarterTurnCounterClockwise is generally defined for classes
 184                  * with components named x_ and y_.
 185                  */
 186                 T z2n = z2_.quarterTurnCounterClockwise( );
 187                 T z3n = z3_.quarterTurnCounterClockwise( );
 188
 189                 //              double k3 = 6 * innerScalar( z3_, z3n ); // Note that this vanishes!
 190                 double k2 = -6 * innerScalar( z3_, z2n ); // (After some simplification.)
 191                 double k1 = 6 * innerScalar( z1_, z3n ); // (After some simplification.)
 192                 double k0 = 2 * innerScalar( z1_, z2n );
 193
 194                 /* If the acceleration ever vanishes, or is close to vanish, then z2_ is parallel to z3_, and k2 becomes very small.
 195                  * When this happens, it is a very small difference between having an inflection and not having one, and in this case we prefer
 196                  * to find the inflection.
 197                  */
 198                 if( k2*k2 < 0.00001 * 36 * ( innerScalar( z3_, z3_ ) * innerScalar( z2_, z2_ ) ) )
 199                         {
 200                                 bezierRootsOfPolynomial( t, k0, k1 );
 201                                 /* In this case we also add the point where the acceleration vanishes to the points of inflection.
 202                                  * We do this by projecting the acceleration on z3_ to get a scalar equation.
 203                                  */
 204                                 double * tmp = & t[0];
 205                                 for( ; *tmp != HUGE_VAL; ++tmp )
 206                                         ;
 207                                 bezierRootsOfPolynomial( tmp, innerScalar( z3_, z2_ ), 3 * innerScalar( z3_, z3_ ) ); // A factor of 2 was taken out.
 208                         }
 209                 else
 210                         {
 211                                 bezierRootsOfPolynomial( t, k0, k1, k2 );
 212                         }
 213         }
 214
 215         template< class T >
 216                 std::ostream &
 217                 operator << ( std::ostream & os, const PolyCoeffs< T > & self )
 218                 {
 219                         os << "t -> " << self.z0 << " + t * " << self.z1 << " + t^{2} * " << self.z2 << " + t^{3} * " << self.z3 ;
 220                         return os;
 221                 }
 222
 223
 224 }