1 /* This file is part of Shapes.
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
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.
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/>.
16 * Copyright 2008 Henrik Tidefelt
24 Bezier::bezierRootsOfPolynomial( double t
[2], double k0
, double k1
)
26 double * dst
= & t
[0];
30 double tmp
= - k0
/ k1
;
31 if( 0 <= tmp
&& tmp
<= 1 )
41 Bezier::bezierRootsOfPolynomial( double t
[3], double k0
, double k1
, double k2
)
43 double * dst
= & t
[0];
49 double tmp
= - k0
/ k1
;
50 if( 0 <= tmp
&& tmp
<= 1 )
61 double r2
= k1
* k1
* 0.25 - k0
;
64 double r
= sqrt( r2
);
66 double tmp
= - 0.5 * k1
- r
;
67 if( 0 <= tmp
&& tmp
<= 1 )
74 double tmp
= - 0.5 * k1
+ r
;
75 if( 0 <= tmp
&& tmp
<= 1 )
87 Bezier::bezierRootsOfPolynomial( double t
[4], double k0
, double k1
, double k2
, double k3
)
89 double * dst
= & t
[0];
96 // If we reach here, the equation is identically solved for all t.
106 double tmp
= - k0
/ k1
;
107 if( 0 <= tmp
&& tmp
<= 1 )
119 double tmp_c
= - 0.5 * k1
;
120 double r2
= tmp_c
* tmp_c
- k0
;
123 double r
= sqrt( r2
);
125 double tmp
= tmp_c
- r
;
126 if( 0 <= tmp
&& tmp
<= 1 )
133 double tmp
= tmp_c
+ r
;
134 if( 0 <= tmp
&& tmp
<= 1 )
150 /* The solution that now follows is an implementation of the first cubic formula
151 * presented on MathWorld.
154 typedef std::complex< double > Complex
;
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
);
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 ) );
166 Complex t1
= w1
- Q
/ w1
- k2
/ 3;
167 Complex t2
= w2
- Q
/ w2
- k2
/ 3;
168 Complex t3
= w3
- Q
/ w3
- k2
/ 3;
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
)
173 double tmp
= t1
.real( );
174 if( 0 <= tmp
&& tmp
<= 1 )
180 if( fabs( t2
.imag( ) ) < IMAG_TOL
)
182 double tmp
= t2
.real( );
183 if( 0 <= tmp
&& tmp
<= 1 )
189 if( fabs( t3
.imag( ) ) < IMAG_TOL
)
191 double tmp
= t3
.real( );
192 if( 0 <= tmp
&& tmp
<= 1 )