git.mutantstargoat.com Git - summerhack/blob - tools/curve_draw/vmath/vmath.cc

   1 /*
   2 This file is part of XRay, a photorealistic 3D rendering library.
   3 Copyright (C) 2005 John Tsiombikas
   4
   5 XRay is free software; you can redistribute it and/or modify
   6 it under the terms of the GNU General Public License as published by
   7 the Free Software Foundation; either version 2 of the License, or
   8 (at your option) any later version.
   9
  10 XRay is distributed in the hope that it will be useful,
  11 but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 GNU General Public License for more details.
  14
  15 You should have received a copy of the GNU General Public License
  16 along with XRay; if not, write to the Free Software
  17 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  18 */
  19
  20 /**
  21  * @file vmath.cxx
  22  * @author John Tsiombikas
  23  * @date 28 June 2005
  24  *
  25  * vmath provides the mathematical foundations for 3D graphics.
  26  * This specific file contains various constants and basic operations
  27  * as well as the Basis class, which represents an orthogonal basis in
  28  * 3 dimensional space.
  29  */
  30
  31 #include "vmath.h"
  32
  33 const scalar_t e = 2.7182818;
  34
  35 const scalar_t pi = 3.1415926535897932;
  36 const scalar_t half_pi = pi / 2.0;
  37 const scalar_t quarter_pi = pi / 4.0;
  38 const scalar_t two_pi = pi * 2.0;
  39 const scalar_t three_half_pi = 3.0 * pi / 2.0;
  40
  41 const scalar_t xsmall_number = 1.e-8;
  42 const scalar_t small_number = 1.e-4;
  43
  44 const scalar_t error_margin = 1.e-6;
  45
  46
  47 /** Generates a random number in [0, 1) */
  48 scalar_t frand(scalar_t range) {
  49         return range * (float)rand() / (float)RAND_MAX;
  50 }
  51
  52 /** Numerical calculation of integrals using simpson's rule */
  53 scalar_t integral(scalar_t (*f)(scalar_t), scalar_t low, scalar_t high, int samples)  {
  54         scalar_t h = (high - low) / (scalar_t)samples;
  55         scalar_t sum = 0.0;
  56
  57         for(int i=0; i<samples+1; i++) {
  58                 scalar_t y = f((scalar_t)i * h + low);
  59                 sum += ((!i || i == samples) ? y : ((i%2) ? 4.0 * y : 2.0 * y)) * (h / 3.0);
  60         }
  61         return sum;
  62 }
  63
  64 /** Gaussuan function */
  65 scalar_t gaussian(scalar_t x, scalar_t mean, scalar_t sdev) {
  66         scalar_t exponent = -SQ(x - mean) / (2.0 * SQ(sdev));
  67
  68         return 1.0 - -pow(e, exponent) / (sdev * sqrt(two_pi));
  69 }
  70
  71
  72 /** b-spline approximation */
  73 scalar_t bspline(const Vector4 &cpvec, scalar_t t) {
  74         Matrix4x4 bspline_mat(  -1,  3, -3,  1,
  75                                                          3, -6,  3,  0,
  76                                                         -3,  0,  3,  0,
  77                                                          1,  4,  1,  0);
  78
  79         scalar_t t_square = t * t;
  80         scalar_t t_cube = t_square * t;
  81         Vector4 params(t_cube, t_square, t, 1.0);
  82
  83         return dot_product(params, cpvec.transformed(bspline_mat) / 6.0);
  84 }
  85
  86 /** Catmull-rom spline interpolation */
  87 scalar_t catmull_rom_spline(const Vector4 &cpvec, scalar_t t) {
  88         Matrix4x4 crspline_mat( -1,  3, -3,  1,
  89                                                          2, -5,  4, -1,
  90                                                         -1,  0,  1,  0,
  91                                                          0,  2,  0,  0);
  92
  93         scalar_t t_square = t * t;
  94         scalar_t t_cube = t_square * t;
  95         Vector4 params(t_cube, t_square, t, 1.0);
  96
  97         return dot_product(params, cpvec.transformed(crspline_mat) / 2.0);
  98 }
  99
 100 /**
 101  * @author Mihalis Georgoulopoulos
 102  *
 103  * Bezier interpolation.
 104  * returns an interpolated scalar value given 4 control values.
 105  */
 106 scalar_t bezier(const Vector4 &cp, scalar_t t)
 107 {
 108         static scalar_t omt, omt3, t3, f;
 109         t3 = t * t * t;
 110         omt = 1.0f - t;
 111         omt3 = omt * omt * omt;
 112         f = 3 * t * omt;
 113
 114         return (cp.x * omt3) + (cp.y * f * omt) + (cp.z * f * t) + (cp.w * t3);
 115 }
 116
 117 /**
 118  * @author Mihalis Georgoulopoulos
 119  *
 120  * Bezier interpolation.
 121  * Vector3 version of the bezier function.
 122  */
 123 Vector3 bezier(const Vector3 &p0, const Vector3 &p1, const Vector3 &p2, const Vector3 &p3, scalar_t t)
 124 {
 125         static scalar_t omt, omt3, t3, f;
 126         t3 = t * t * t;
 127         omt = 1.0f - t;
 128         omt3 = omt * omt * omt;
 129         f = 3 * t * omt;
 130
 131         return (p0 * omt3) + (p1 * f * omt) + (p2 * f * t) + (p3 * t3);
 132 }
 133
 134 /**
 135  * @author Mihalis Georgoulopoulos
 136  *
 137  * Bezier tangent calculation.
 138  * returns a vector tangent to the bezier curve at the specified point.
 139  */
 140 Vector3 bezier_tangent(const Vector3 &p0, const Vector3 &p1, const Vector3 &p2, const Vector3 &p3, scalar_t t)
 141 {
 142         static scalar_t omt;
 143         omt = 1.0f - t;
 144         static Vector3 p4, p5, p6, p7, p8;
 145         p4 = p0 * omt + p1 * t;
 146         p5 = p1 * omt + p2 * t;
 147         p6 = p2 * omt + p3 * t;
 148
 149         p7 = p4 * omt + p5 * t;
 150         p8 = p5 * omt + p6 * t;
 151
 152         return p8 - p7;
 153 }
 154
 155 Basis::Basis() {
 156         i = Vector3(1, 0, 0);
 157         j = Vector3(0, 1, 0);
 158         k = Vector3(0, 0, 1);
 159 }
 160
 161 Basis::Basis(const Vector3 &i, const Vector3 &j, const Vector3 &k) {
 162         this->i = i;
 163         this->j = j;
 164         this->k = k;
 165 }
 166
 167 Basis::Basis(const Vector3 &dir, bool left_handed) {
 168         k = dir;
 169         j = Vector3(0, 1, 0);
 170         i = cross_product(j, k);
 171         j = cross_product(k, i);
 172 }
 173
 174 /** Rotate with euler angles */
 175 void Basis::rotate(scalar_t x, scalar_t y, scalar_t z) {
 176         Matrix4x4 RotMat;
 177         RotMat.set_rotation(Vector3(x, y, z));
 178         i.transform(RotMat);
 179         j.transform(RotMat);
 180         k.transform(RotMat);
 181 }
 182
 183 /** Rotate by axis-angle specification */
 184 void Basis::rotate(const Vector3 &axis, scalar_t angle) {
 185         Quaternion q;
 186         q.set_rotation(axis, angle);
 187         i.transform(q);
 188         j.transform(q);
 189         k.transform(q);
 190 }
 191
 192 /** Rotate with a 4x4 matrix */
 193 void Basis::rotate(const Matrix4x4 &mat) {
 194         i.transform(mat);
 195         j.transform(mat);
 196         k.transform(mat);
 197 }
 198
 199 /** Rotate with a quaternion */
 200 void Basis::rotate(const Quaternion &quat) {
 201         i.transform(quat);
 202         j.transform(quat);
 203         k.transform(quat);
 204 }
 205
 206 /** Extract a rotation matrix from the orthogonal basis */
 207 Matrix3x3 Basis::create_rotation_matrix() const {
 208         return Matrix3x3(       i.x, j.x, k.x,
 209                                                 i.y, j.y, k.y,
 210                                                 i.z, j.z, k.z);
 211 }