remove RSA's md4.c, replace by DP's
[xonotic/netradiant.git] / libs / math / curve.h
1 /*
2 Copyright (C) 2001-2006, William Joseph.
3 All Rights Reserved.
4
5 This file is part of GtkRadiant.
6
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.
11
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.
16
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
20 */
21
22 #if !defined(INCLUDED_MATH_CURVE_H)
23 #define INCLUDED_MATH_CURVE_H
24
25 /// \file
26 /// \brief Curve data types and related operations.
27
28 #include "debugging/debugging.h"
29 #include "container/array.h"
30 #include <math/matrix.h>
31
32
33 template<typename I, typename Degree>
34 struct BernsteinPolynomial
35 {
36   static double apply(double t)
37   {
38     return 1; // general case not implemented
39   }
40 };
41
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;
47
48 template<>
49 struct BernsteinPolynomial<Zero, Zero>
50 {
51   static double apply(double t)
52   {
53     return 1;
54   }
55 };
56
57 template<>
58 struct BernsteinPolynomial<Zero, One>
59 {
60   static double apply(double t)
61   {
62     return 1 - t;
63   }
64 };
65
66 template<>
67 struct BernsteinPolynomial<One, One>
68 {
69   static double apply(double t)
70   {
71     return t;
72   }
73 };
74
75 template<>
76 struct BernsteinPolynomial<Zero, Two>
77 {
78   static double apply(double t)
79   {
80     return (1 - t) * (1 - t);
81   }
82 };
83
84 template<>
85 struct BernsteinPolynomial<One, Two>
86 {
87   static double apply(double t)
88   {
89     return 2 * (1 - t) * t;
90   }
91 };
92
93 template<>
94 struct BernsteinPolynomial<Two, Two>
95 {
96   static double apply(double t)
97   {
98     return t * t;
99   }
100 };
101
102 template<>
103 struct BernsteinPolynomial<Zero, Three>
104 {
105   static double apply(double t)
106   {
107     return (1 - t) * (1 - t) * (1 - t);
108   }
109 };
110
111 template<>
112 struct BernsteinPolynomial<One, Three>
113 {
114   static double apply(double t)
115   {
116     return 3 * (1 - t) * (1 - t) * t;
117   }
118 };
119
120 template<>
121 struct BernsteinPolynomial<Two, Three>
122 {
123   static double apply(double t)
124   {
125     return 3 * (1 - t) * t * t;
126   }
127 };
128
129 template<>
130 struct BernsteinPolynomial<Three, Three>
131 {
132   static double apply(double t)
133   {
134     return t * t * t;
135   }
136 };
137
138 typedef Array<Vector3> ControlPoints;
139
140 inline Vector3 CubicBezier_evaluate(const Vector3* firstPoint, double t)
141 {
142   Vector3 result(0, 0, 0);
143   double denominator = 0;
144
145   {
146     double weight = BernsteinPolynomial<Zero, Three>::apply(t);
147     result += vector3_scaled(*firstPoint++, weight);
148     denominator += weight;
149   }
150   {
151     double weight = BernsteinPolynomial<One, Three>::apply(t);
152     result += vector3_scaled(*firstPoint++, weight);
153     denominator += weight;
154   }
155   {
156     double weight = BernsteinPolynomial<Two, Three>::apply(t);
157     result += vector3_scaled(*firstPoint++, weight);
158     denominator += weight;
159   }
160   {
161     double weight = BernsteinPolynomial<Three, Three>::apply(t);
162     result += vector3_scaled(*firstPoint++, weight);
163     denominator += weight;
164   }
165
166   return result / denominator;
167 }
168
169 inline Vector3 CubicBezier_evaluateMid(const Vector3* firstPoint)
170 {
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);
175 }
176
177 inline Vector3 CatmullRom_evaluate(const ControlPoints& controlPoints, double t)
178 {
179   // scale t to be segment-relative
180   t *= double(controlPoints.size() - 1);
181
182   // subtract segment index;
183   std::size_t segment = 0;
184   for(std::size_t i = 0; i < controlPoints.size() - 1; ++i)
185   {
186     if(t <= double(i+1))
187     {
188       segment = i;
189       break;
190     }
191   }
192   t -= segment;
193
194   const double reciprocal_alpha = 1.0 / 3.0;
195
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);
206 }
207
208 typedef Array<float> Knots;
209
210 inline double BSpline_basis(const Knots& knots, std::size_t i, std::size_t degree, double t)
211 {
212   if(degree == 0)
213   {
214     if(knots[i] <= t
215       && t < knots[i + 1]
216       && knots[i] < knots[i + 1])
217     {
218       return 1;
219     }
220     return 0;
221   }
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);
226   return left + right;
227 }
228
229 inline Vector3 BSpline_evaluate(const ControlPoints& controlPoints, const Knots& knots, std::size_t degree, double t)
230 {
231   Vector3 result(0, 0, 0);
232   for(std::size_t i = 0; i < controlPoints.size(); ++i)
233   {
234     result += vector3_scaled(controlPoints[i], BSpline_basis(knots, i, degree, t));
235   }
236   return result;
237 }
238
239 typedef Array<float> NURBSWeights;
240
241 inline Vector3 NURBS_evaluate(const ControlPoints& controlPoints, const NURBSWeights& weights, const Knots& knots, std::size_t degree, double t)
242 {
243   Vector3 result(0, 0, 0);
244   double denominator = 0;
245   for(std::size_t i = 0; i < controlPoints.size(); ++i)
246   {
247     double weight = weights[i] * BSpline_basis(knots, i, degree, t);
248     result += vector3_scaled(controlPoints[i], weight);
249     denominator += weight;
250   }
251   return result / denominator;
252 }
253
254 inline void KnotVector_openUniform(Knots& knots, std::size_t count, std::size_t degree)
255 {
256   knots.resize(count + degree + 1);
257
258   std::size_t equalKnots = 1;
259
260   for(std::size_t i = 0; i < equalKnots; ++i)
261   {
262     knots[i] = 0;
263     knots[knots.size() - (i + 1)] = 1;
264   }
265
266   std::size_t difference = knots.size() - 2 * (equalKnots);
267   for(std::size_t i = 0; i < difference; ++i)
268   {
269           knots[i + equalKnots] = Knots::value_type(double(i + 1) * 1.0 / double(difference + 1));
270   }
271 }
272
273 #endif