minor gitignore additions
[summerhack] / 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 }