2 Copyright (C) 2001-2006, William Joseph.
5 This file is part of GtkRadiant.
7 GtkRadiant is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
12 GtkRadiant is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
17 You should have received a copy of the GNU General Public License
18 along with GtkRadiant; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
22 #if !defined( INCLUDED_MATH_CURVE_H )
23 #define INCLUDED_MATH_CURVE_H
26 /// \brief Curve data types and related operations.
28 #include "debugging/debugging.h"
29 #include "container/array.h"
30 #include <math/matrix.h>
33 template<typename I, typename Degree>
34 struct BernsteinPolynomial
36 static double apply( double t ){
37 return 1; // general case not implemented
41 typedef IntegralConstant<0> Zero;
42 typedef IntegralConstant<1> One;
43 typedef IntegralConstant<2> Two;
44 typedef IntegralConstant<3> Three;
45 typedef IntegralConstant<4> Four;
48 struct BernsteinPolynomial<Zero, Zero>
50 static double apply( double t ){
56 struct BernsteinPolynomial<Zero, One>
58 static double apply( double t ){
64 struct BernsteinPolynomial<One, One>
66 static double apply( double t ){
72 struct BernsteinPolynomial<Zero, Two>
74 static double apply( double t ){
75 return ( 1 - t ) * ( 1 - t );
80 struct BernsteinPolynomial<One, Two>
82 static double apply( double t ){
83 return 2 * ( 1 - t ) * t;
88 struct BernsteinPolynomial<Two, Two>
90 static double apply( double t ){
96 struct BernsteinPolynomial<Zero, Three>
98 static double apply( double t ){
99 return ( 1 - t ) * ( 1 - t ) * ( 1 - t );
104 struct BernsteinPolynomial<One, Three>
106 static double apply( double t ){
107 return 3 * ( 1 - t ) * ( 1 - t ) * t;
112 struct BernsteinPolynomial<Two, Three>
114 static double apply( double t ){
115 return 3 * ( 1 - t ) * t * t;
120 struct BernsteinPolynomial<Three, Three>
122 static double apply( double t ){
127 typedef Array<Vector3> ControlPoints;
129 inline Vector3 CubicBezier_evaluate( const Vector3* firstPoint, double t ){
130 Vector3 result( 0, 0, 0 );
131 double denominator = 0;
134 double weight = BernsteinPolynomial<Zero, Three>::apply( t );
135 result += vector3_scaled( *firstPoint++, weight );
136 denominator += weight;
139 double weight = BernsteinPolynomial<One, Three>::apply( t );
140 result += vector3_scaled( *firstPoint++, weight );
141 denominator += weight;
144 double weight = BernsteinPolynomial<Two, Three>::apply( t );
145 result += vector3_scaled( *firstPoint++, weight );
146 denominator += weight;
149 double weight = BernsteinPolynomial<Three, Three>::apply( t );
150 result += vector3_scaled( *firstPoint++, weight );
151 denominator += weight;
154 return result / denominator;
157 inline Vector3 CubicBezier_evaluateMid( const Vector3* firstPoint ){
158 return vector3_scaled( firstPoint[0], 0.125 )
159 + vector3_scaled( firstPoint[1], 0.375 )
160 + vector3_scaled( firstPoint[2], 0.375 )
161 + vector3_scaled( firstPoint[3], 0.125 );
164 inline Vector3 CatmullRom_evaluate( const ControlPoints& controlPoints, double t ){
165 // scale t to be segment-relative
166 t *= double(controlPoints.size() - 1);
168 // subtract segment index;
169 std::size_t segment = 0;
170 for ( std::size_t i = 0; i < controlPoints.size() - 1; ++i )
172 if ( t <= double(i + 1) ) {
179 const double reciprocal_alpha = 1.0 / 3.0;
181 Vector3 bezierPoints[4];
182 bezierPoints[0] = controlPoints[segment];
183 bezierPoints[1] = ( segment > 0 )
184 ? controlPoints[segment] + vector3_scaled( controlPoints[segment + 1] - controlPoints[segment - 1], reciprocal_alpha * 0.5 )
185 : controlPoints[segment] + vector3_scaled( controlPoints[segment + 1] - controlPoints[segment], reciprocal_alpha );
186 bezierPoints[2] = ( segment < controlPoints.size() - 2 )
187 ? controlPoints[segment + 1] + vector3_scaled( controlPoints[segment] - controlPoints[segment + 2], reciprocal_alpha * 0.5 )
188 : controlPoints[segment + 1] + vector3_scaled( controlPoints[segment] - controlPoints[segment + 1], reciprocal_alpha );
189 bezierPoints[3] = controlPoints[segment + 1];
190 return CubicBezier_evaluate( bezierPoints, t );
193 typedef Array<float> Knots;
195 inline double BSpline_basis( const Knots& knots, std::size_t i, std::size_t degree, double t ){
199 && knots[i] < knots[i + 1] ) {
204 double leftDenom = knots[i + degree] - knots[i];
205 double left = ( leftDenom == 0 ) ? 0 : ( ( t - knots[i] ) / leftDenom ) * BSpline_basis( knots, i, degree - 1, t );
206 double rightDenom = knots[i + degree + 1] - knots[i + 1];
207 double right = ( rightDenom == 0 ) ? 0 : ( ( knots[i + degree + 1] - t ) / rightDenom ) * BSpline_basis( knots, i + 1, degree - 1, t );
211 inline Vector3 BSpline_evaluate( const ControlPoints& controlPoints, const Knots& knots, std::size_t degree, double t ){
212 Vector3 result( 0, 0, 0 );
213 for ( std::size_t i = 0; i < controlPoints.size(); ++i )
215 result += vector3_scaled( controlPoints[i], BSpline_basis( knots, i, degree, t ) );
220 typedef Array<float> NURBSWeights;
222 inline Vector3 NURBS_evaluate( const ControlPoints& controlPoints, const NURBSWeights& weights, const Knots& knots, std::size_t degree, double t ){
223 Vector3 result( 0, 0, 0 );
224 double denominator = 0;
225 for ( std::size_t i = 0; i < controlPoints.size(); ++i )
227 double weight = weights[i] * BSpline_basis( knots, i, degree, t );
228 result += vector3_scaled( controlPoints[i], weight );
229 denominator += weight;
231 return result / denominator;
234 inline void KnotVector_openUniform( Knots& knots, std::size_t count, std::size_t degree ){
235 knots.resize( count + degree + 1 );
237 std::size_t equalKnots = 1;
239 for ( std::size_t i = 0; i < equalKnots; ++i )
242 knots[knots.size() - ( i + 1 )] = 1;
245 std::size_t difference = knots.size() - 2 * ( equalKnots );
246 for ( std::size_t i = 0; i < difference; ++i )
248 knots[i + equalKnots] = Knots::value_type( double(i + 1) * 1.0 / double(difference + 1) );