source/bezier.cc

   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 Henrik Tidefelt
  17  */
  18
  19 #include "bezier.h"
  20
  21 #include <complex>
  22
  23 void
  24 Bezier::bezierRootsOfPolynomial( double t[2], double k0, double k1 )
  25 {
  26         double * dst = & t[0];
  27
  28         if( k1 != 0 )
  29                 {
  30                         double tmp = - k0 / k1;
  31                         if( 0 <= tmp && tmp <= 1 )
  32                                 {
  33                                         *dst = tmp;
  34                                         ++dst;
  35                                 }
  36                 }
  37         *dst = HUGE_VAL;
  38 }
  39
  40 void
  41 Bezier::bezierRootsOfPolynomial( double t[3], double k0, double k1, double k2 )
  42 {
  43         double * dst = & t[0];
  44
  45         if( k2 == 0 )
  46                 {
  47                         if( k1 != 0 )
  48                                 {
  49                                         double tmp = - k0 / k1;
  50                                         if( 0 <= tmp && tmp <= 1 )
  51                                                 {
  52                                                         *dst = tmp;
  53                                                         ++dst;
  54                                                 }
  55                                 }
  56                 }
  57         else
  58                 {
  59                         k0 /= k2;
  60                         k1 /= k2;
  61                         double r2 = k1 * k1 * 0.25 - k0;
  62                         if( r2 >= 0 )
  63                                 {
  64                                         double r = sqrt( r2 );
  65                                         {
  66                                                 double tmp = - 0.5 * k1 - r;
  67                                                 if( 0 <= tmp && tmp <= 1 )
  68                                                         {
  69                                                                 *dst = tmp;
  70                                                                 ++dst;
  71                                                         }
  72                                         }
  73                                         {
  74                                                 double tmp = - 0.5 * k1 + r;
  75                                                 if( 0 <= tmp && tmp <= 1 )
  76                                                         {
  77                                                                 *dst = tmp;
  78                                                                 ++dst;
  79                                                         }
  80                                         }
  81                                 }
  82                 }
  83         *dst = HUGE_VAL;
  84 }
  85
  86 void
  87 Bezier::bezierRootsOfPolynomial( double t[4], double k0, double k1, double k2, double k3 )
  88 {
  89         double * dst = & t[0];
  90         if( k3 == 0 )
  91                 {
  92                         if( k2 == 0 )
  93                                 {
  94                                         if( k1 == 0 )
  95                                                 {
  96                                                         // If we reach here, the equation is identically solved for all t.
  97                                                         if( k0 == 0 )
  98                                                                 {
  99                                                                         *dst = 0;
 100                                                                         ++dst;
 101                                                                 }
 102                                                 }
 103                                         else
 104                                                 {
 105                                                         // Here k1 != 0
 106                                                         double tmp = - k0 / k1;
 107                                                         if( 0 <= tmp && tmp <= 1 )
 108                                                                 {
 109                                                                         *dst = tmp;
 110                                                                         ++dst;
 111                                                                 }
 112                                                 }
 113                                 }
 114                         else
 115                                 {
 116                                         // Here k2 != 0
 117                                         k1 /= k2;
 118                                         k0 /= k2;
 119                                         double tmp_c = - 0.5 * k1;
 120                                         double r2 = tmp_c * tmp_c - k0;
 121                                         if( r2 >= 0 )
 122                                                 {
 123                                                         double r = sqrt( r2 );
 124                                                         {
 125                                                                 double tmp = tmp_c - r;
 126                                                                 if( 0 <= tmp && tmp <= 1 )
 127                                                                         {
 128                                                                                 *dst = tmp;
 129                                                                                 ++dst;
 130                                                                         }
 131                                                         }
 132                                                         {
 133                                                                 double tmp = tmp_c + r;
 134                                                                 if( 0 <= tmp && tmp <= 1 )
 135                                                                         {
 136                                                                                 *dst = tmp;
 137                                                                                 ++dst;
 138                                                                         }
 139                                                         }
 140                                                 }
 141                                 }
 142                         *dst = HUGE_VAL;
 143                         return;
 144                 }
 145
 146         k2 /= k3;
 147         k1 /= k3;
 148         k0 /= k3;
 149
 150         /* The solution that now follows is an implementation of the first cubic formula
 151          * presented on MathWorld.
 152          */
 153
 154         typedef std::complex< double > Complex;
 155
 156         Complex Q( ( 3 * k1 - k2 * k2 ) / 9, 0 );
 157         Complex R( 0.5 * ( ( 9 * k1 * k2 - 2 * k2 * k2 * k2 ) / 27 - k0 ), 0 );
 158         Complex r = sqrt( R * R + Q * Q * Q );
 159
 160         Complex wCube;
 161         wCube = R + r;  // ( R - r ) is also an option
 162         Complex w1 = pow( wCube, 1./3 );
 163         Complex w2 = w1 * exp( Complex( 0, 2 * M_PI / 3 ) );
 164         Complex w3 = w1 * exp( Complex( 0, -2 * M_PI / 3 ) );
 165
 166         Complex t1 = w1 - Q / w1 - k2 / 3;
 167         Complex t2 = w2 - Q / w2 - k2 / 3;
 168         Complex t3 = w3 - Q / w3 - k2 / 3;
 169
 170         const double IMAG_TOL = 1e-10; // This is "relative" to the interesting t-range being 0..1
 171         if( fabs( t1.imag( ) ) < IMAG_TOL )
 172                 {
 173                         double tmp = t1.real( );
 174                         if( 0 <= tmp && tmp <= 1 )
 175                                 {
 176                                         *dst = tmp;
 177                                         ++dst;
 178                                 }
 179                 }
 180         if( fabs( t2.imag( ) ) < IMAG_TOL )
 181                 {
 182                         double tmp = t2.real( );
 183                         if( 0 <= tmp && tmp <= 1 )
 184                                 {
 185                                         *dst = tmp;
 186                                         ++dst;
 187                                 }
 188                 }
 189         if( fabs( t3.imag( ) ) < IMAG_TOL )
 190                 {
 191                         double tmp = t3.real( );
 192                         if( 0 <= tmp && tmp <= 1 )
 193                                 {
 194                                         *dst = tmp;
 195                                         ++dst;
 196                                 }
 197                 }
 198         *dst = HUGE_VAL;
 199 }