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)
38 return 1; // general case not implemented
42 typedef IntegralConstant<0> Zero;
43 typedef IntegralConstant<1> One;
44 typedef IntegralConstant<2> Two;
45 typedef IntegralConstant<3> Three;
46 typedef IntegralConstant<4> Four;
49 struct BernsteinPolynomial<Zero, Zero>
51 static double apply(double t)
58 struct BernsteinPolynomial<Zero, One>
60 static double apply(double t)
67 struct BernsteinPolynomial<One, One>
69 static double apply(double t)
76 struct BernsteinPolynomial<Zero, Two>
78 static double apply(double t)
80 return (1 - t) * (1 - t);
85 struct BernsteinPolynomial<One, Two>
87 static double apply(double t)
89 return 2 * (1 - t) * t;
94 struct BernsteinPolynomial<Two, Two>
96 static double apply(double t)
103 struct BernsteinPolynomial<Zero, Three>
105 static double apply(double t)
107 return (1 - t) * (1 - t) * (1 - t);
112 struct BernsteinPolynomial<One, Three>
114 static double apply(double t)
116 return 3 * (1 - t) * (1 - t) * t;
121 struct BernsteinPolynomial<Two, Three>
123 static double apply(double t)
125 return 3 * (1 - t) * t * t;
130 struct BernsteinPolynomial<Three, Three>
132 static double apply(double t)
138 typedef Array<Vector3> ControlPoints;
140 inline Vector3 CubicBezier_evaluate(const Vector3* firstPoint, double t)
142 Vector3 result(0, 0, 0);
143 double denominator = 0;
146 double weight = BernsteinPolynomial<Zero, Three>::apply(t);
147 result += vector3_scaled(*firstPoint++, weight);
148 denominator += weight;
151 double weight = BernsteinPolynomial<One, Three>::apply(t);
152 result += vector3_scaled(*firstPoint++, weight);
153 denominator += weight;
156 double weight = BernsteinPolynomial<Two, Three>::apply(t);
157 result += vector3_scaled(*firstPoint++, weight);
158 denominator += weight;
161 double weight = BernsteinPolynomial<Three, Three>::apply(t);
162 result += vector3_scaled(*firstPoint++, weight);
163 denominator += weight;
166 return result / denominator;
169 inline Vector3 CubicBezier_evaluateMid(const Vector3* firstPoint)
171 return vector3_scaled(firstPoint[0], 0.125)
172 + vector3_scaled(firstPoint[1], 0.375)
173 + vector3_scaled(firstPoint[2], 0.375)
174 + vector3_scaled(firstPoint[3], 0.125);
177 inline Vector3 CatmullRom_evaluate(const ControlPoints& controlPoints, double t)
179 // scale t to be segment-relative
180 t *= double(controlPoints.size() - 1);
182 // subtract segment index;
183 std::size_t segment = 0;
184 for(std::size_t i = 0; i < controlPoints.size() - 1; ++i)
194 const double reciprocal_alpha = 1.0 / 3.0;
196 Vector3 bezierPoints[4];
197 bezierPoints[0] = controlPoints[segment];
198 bezierPoints[1] = (segment > 0)
199 ? controlPoints[segment] + vector3_scaled(controlPoints[segment + 1] - controlPoints[segment - 1], reciprocal_alpha * 0.5)
200 : controlPoints[segment] + vector3_scaled(controlPoints[segment + 1] - controlPoints[segment], reciprocal_alpha);
201 bezierPoints[2] = (segment < controlPoints.size() - 2)
202 ? controlPoints[segment + 1] + vector3_scaled(controlPoints[segment] - controlPoints[segment + 2], reciprocal_alpha * 0.5)
203 : controlPoints[segment + 1] + vector3_scaled(controlPoints[segment] - controlPoints[segment + 1], reciprocal_alpha);
204 bezierPoints[3] = controlPoints[segment + 1];
205 return CubicBezier_evaluate(bezierPoints, t);
208 typedef Array<float> Knots;
210 inline double BSpline_basis(const Knots& knots, std::size_t i, std::size_t degree, double t)
216 && knots[i] < knots[i + 1])
222 double leftDenom = knots[i + degree] - knots[i];
223 double left = (leftDenom == 0) ? 0 : ((t - knots[i]) / leftDenom) * BSpline_basis(knots, i, degree - 1, t);
224 double rightDenom = knots[i + degree + 1] - knots[i + 1];
225 double right = (rightDenom == 0) ? 0 : ((knots[i + degree + 1] - t) / rightDenom) * BSpline_basis(knots, i + 1, degree - 1, t);
229 inline Vector3 BSpline_evaluate(const ControlPoints& controlPoints, const Knots& knots, std::size_t degree, double t)
231 Vector3 result(0, 0, 0);
232 for(std::size_t i = 0; i < controlPoints.size(); ++i)
234 result += vector3_scaled(controlPoints[i], BSpline_basis(knots, i, degree, t));
239 typedef Array<float> NURBSWeights;
241 inline Vector3 NURBS_evaluate(const ControlPoints& controlPoints, const NURBSWeights& weights, const Knots& knots, std::size_t degree, double t)
243 Vector3 result(0, 0, 0);
244 double denominator = 0;
245 for(std::size_t i = 0; i < controlPoints.size(); ++i)
247 double weight = weights[i] * BSpline_basis(knots, i, degree, t);
248 result += vector3_scaled(controlPoints[i], weight);
249 denominator += weight;
251 return result / denominator;
254 inline void KnotVector_openUniform(Knots& knots, std::size_t count, std::size_t degree)
256 knots.resize(count + degree + 1);
258 std::size_t equalKnots = 1;
260 for(std::size_t i = 0; i < equalKnots; ++i)
263 knots[knots.size() - (i + 1)] = 1;
266 std::size_t difference = knots.size() - 2 * (equalKnots);
267 for(std::size_t i = 0; i < difference; ++i)
269 knots[i + equalKnots] = Knots::value_type(double(i + 1) * 1.0 / double(difference + 1));