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                 return 1; // general case not implemented
  38         }
  39 };
  40
  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;
  46
  47 template<>
  48 struct BernsteinPolynomial<Zero, Zero>
  49 {
  50         static double apply( double t ){
  51                 return 1;
  52         }
  53 };
  54
  55 template<>
  56 struct BernsteinPolynomial<Zero, One>
  57 {
  58         static double apply( double t ){
  59                 return 1 - t;
  60         }
  61 };
  62
  63 template<>
  64 struct BernsteinPolynomial<One, One>
  65 {
  66         static double apply( double t ){
  67                 return t;
  68         }
  69 };
  70
  71 template<>
  72 struct BernsteinPolynomial<Zero, Two>
  73 {
  74         static double apply( double t ){
  75                 return ( 1 - t ) * ( 1 - t );
  76         }
  77 };
  78
  79 template<>
  80 struct BernsteinPolynomial<One, Two>
  81 {
  82         static double apply( double t ){
  83                 return 2 * ( 1 - t ) * t;
  84         }
  85 };
  86
  87 template<>
  88 struct BernsteinPolynomial<Two, Two>
  89 {
  90         static double apply( double t ){
  91                 return t * t;
  92         }
  93 };
  94
  95 template<>
  96 struct BernsteinPolynomial<Zero, Three>
  97 {
  98         static double apply( double t ){
  99                 return ( 1 - t ) * ( 1 - t ) * ( 1 - t );
 100         }
 101 };
 102
 103 template<>
 104 struct BernsteinPolynomial<One, Three>
 105 {
 106         static double apply( double t ){
 107                 return 3 * ( 1 - t ) * ( 1 - t ) * t;
 108         }
 109 };
 110
 111 template<>
 112 struct BernsteinPolynomial<Two, Three>
 113 {
 114         static double apply( double t ){
 115                 return 3 * ( 1 - t ) * t * t;
 116         }
 117 };
 118
 119 template<>
 120 struct BernsteinPolynomial<Three, Three>
 121 {
 122         static double apply( double t ){
 123                 return t * t * t;
 124         }
 125 };
 126
 127 typedef Array<Vector3> ControlPoints;
 128
 129 inline Vector3 CubicBezier_evaluate( const Vector3* firstPoint, double t ){
 130         Vector3 result( 0, 0, 0 );
 131         double denominator = 0;
 132
 133         {
 134                 double weight = BernsteinPolynomial<Zero, Three>::apply( t );
 135                 result += vector3_scaled( *firstPoint++, weight );
 136                 denominator += weight;
 137         }
 138         {
 139                 double weight = BernsteinPolynomial<One, Three>::apply( t );
 140                 result += vector3_scaled( *firstPoint++, weight );
 141                 denominator += weight;
 142         }
 143         {
 144                 double weight = BernsteinPolynomial<Two, Three>::apply( t );
 145                 result += vector3_scaled( *firstPoint++, weight );
 146                 denominator += weight;
 147         }
 148         {
 149                 double weight = BernsteinPolynomial<Three, Three>::apply( t );
 150                 result += vector3_scaled( *firstPoint++, weight );
 151                 denominator += weight;
 152         }
 153
 154         return result / denominator;
 155 }
 156
 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 );
 162 }
 163
 164 inline Vector3 CatmullRom_evaluate( const ControlPoints& controlPoints, double t ){
 165         // scale t to be segment-relative
 166         t *= double(controlPoints.size() - 1);
 167
 168         // subtract segment index;
 169         std::size_t segment = 0;
 170         for ( std::size_t i = 0; i < controlPoints.size() - 1; ++i )
 171         {
 172                 if ( t <= double(i + 1) ) {
 173                         segment = i;
 174                         break;
 175                 }
 176         }
 177         t -= segment;
 178
 179         const double reciprocal_alpha = 1.0 / 3.0;
 180
 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 );
 191 }
 192
 193 typedef Array<float> Knots;
 194
 195 inline double BSpline_basis( const Knots& knots, std::size_t i, std::size_t degree, double t ){
 196         if ( degree == 0 ) {
 197                 if ( knots[i] <= t
 198                          && t < knots[i + 1]
 199                          && knots[i] < knots[i + 1] ) {
 200                         return 1;
 201                 }
 202                 return 0;
 203         }
 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 );
 208         return left + right;
 209 }
 210
 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 )
 214         {
 215                 result += vector3_scaled( controlPoints[i], BSpline_basis( knots, i, degree, t ) );
 216         }
 217         return result;
 218 }
 219
 220 typedef Array<float> NURBSWeights;
 221
 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 )
 226         {
 227                 double weight = weights[i] * BSpline_basis( knots, i, degree, t );
 228                 result += vector3_scaled( controlPoints[i], weight );
 229                 denominator += weight;
 230         }
 231         return result / denominator;
 232 }
 233
 234 inline void KnotVector_openUniform( Knots& knots, std::size_t count, std::size_t degree ){
 235         knots.resize( count + degree + 1 );
 236
 237         std::size_t equalKnots = 1;
 238
 239         for ( std::size_t i = 0; i < equalKnots; ++i )
 240         {
 241                 knots[i] = 0;
 242                 knots[knots.size() - ( i + 1 )] = 1;
 243         }
 244
 245         std::size_t difference = knots.size() - 2 * ( equalKnots );
 246         for ( std::size_t i = 0; i < difference; ++i )
 247         {
 248                 knots[i + equalKnots] = Knots::value_type( double(i + 1) * 1.0 / double(difference + 1) );
 249         }
 250 }
 251
 252 #endif