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