Line data Source code
1 : /**
2 : * @file geomexts.h
3 : * @brief Internal engine geometry structures not needed in the API
4 : */
5 : #ifndef GEOMEXTS_H_
6 : #define GEOMEXTS_H_
7 :
8 : /**
9 : * @brief half precision floating point object, similar to `float` but with less data bits
10 : */
11 : struct half
12 : {
13 : ushort val;
14 :
15 48 : half() {}
16 22 : half(float f) //truncates a 32 bit float to a 16 bit half
17 22 : {
18 : union
19 : {
20 : int i;
21 : float f;
22 : } conv;
23 22 : conv.f = f;
24 22 : ushort signbit = (conv.i>>(31-15)) & (1<<15),
25 22 : mantissa = (conv.i>>(23-10)) & 0x3FF;
26 22 : int exponent = ((conv.i>>23)&0xFF) - 127 + 15;
27 22 : if(exponent <= 0)
28 : {
29 8 : mantissa |= 0x400;
30 8 : mantissa >>= min(1-exponent, 10+1);
31 8 : exponent = 0;
32 : }
33 14 : else if(exponent >= 0x1F)
34 : {
35 0 : mantissa = 0;
36 0 : exponent = 0x1F;
37 : }
38 22 : val = signbit | (static_cast<ushort>(exponent)<<10) | mantissa;
39 22 : }
40 :
41 2 : bool operator==(const half &h) const { return val == h.val; }
42 1 : bool operator!=(const half &h) const { return val != h.val; }
43 : };
44 :
45 : struct triangle
46 : {
47 : vec a, b, c;
48 :
49 7 : triangle(const vec &a, const vec &b, const vec &c) : a(a), b(b), c(c) {}
50 1 : triangle() {}
51 :
52 1 : triangle &add(const vec &o)
53 : {
54 1 : a.add(o);
55 1 : b.add(o);
56 1 : c.add(o);
57 1 : return *this;
58 : }
59 :
60 1 : triangle &sub(const vec &o)
61 : {
62 1 : a.sub(o);
63 1 : b.sub(o);
64 1 : c.sub(o);
65 1 : return *this;
66 : }
67 :
68 6 : bool operator==(const triangle &t) const
69 : {
70 6 : return a == t.a && b == t.b && c == t.c;
71 : }
72 : };
73 :
74 : /**
75 : * @brief An object in four-dimensional complex space.
76 : *
77 : * This object is a four component "vector" with three imaginary components and
78 : * one scalar component used for object rotations: quats have 6 DoF among their
79 : * four terms, avoiding gimbal lock inherent to Euler angles.
80 : *
81 : * Follows the Hamiltonian convention (i*j = k) for dimensional relationships.
82 : *
83 : * x represents the i component
84 : * y represents the j component
85 : * z represents the k component
86 : * w represents the real (scalar) component
87 : */
88 : struct quat : vec4<float>
89 : {
90 : public:
91 558 : quat() {}
92 249 : quat(float x, float y, float z, float w) : vec4<float>(x, y, z, w) {}
93 4 : quat(const vec &axis, float angle)
94 4 : {
95 4 : w = cosf(angle/2);
96 4 : float s = std::sin(angle/2);
97 4 : x = s*axis.x;
98 4 : y = s*axis.y;
99 4 : z = s*axis.z;
100 4 : }
101 3 : quat(const vec &u, const vec &v)
102 3 : {
103 3 : w = sqrtf(u.squaredlen() * v.squaredlen()) + u.dot(v);
104 3 : cross(u, v);
105 3 : normalize();
106 3 : }
107 3 : explicit quat(const vec &v)
108 3 : {
109 3 : x = v.x;
110 3 : y = v.y;
111 3 : z = v.z;
112 3 : restorew();
113 3 : }
114 442 : explicit quat(const matrix3 &m) { convertmatrix(m); }
115 3 : explicit quat(const matrix4x3 &m) { convertmatrix(m); }
116 1 : explicit quat(const matrix4 &m) { convertmatrix(m); }
117 :
118 15 : void restorew() { w = 1.0f-x*x-y*y-z*z; w = w<0 ? 0 : -sqrtf(w); }
119 :
120 42 : quat &add(const vec4<float> &o) { vec4<float>::add(o); return *this; }
121 106 : quat &sub(const vec4<float> &o) { vec4<float>::sub(o); return *this; }
122 14 : quat &mul(float k) { vec4<float>::mul(k); return *this; }
123 16 : quat &madd(const vec4<float> &a, const float &b) { return add(vec4<float>(a).mul(b)); }
124 11 : quat &msub(const vec4<float> &a, const float &b) { return sub(vec4<float>(a).mul(b)); }
125 :
126 90 : quat &mul(const quat &p, const quat &o)
127 : {
128 90 : x = p.w*o.x + p.x*o.w + p.y*o.z - p.z*o.y;
129 90 : y = p.w*o.y - p.x*o.z + p.y*o.w + p.z*o.x;
130 90 : z = p.w*o.z + p.x*o.y - p.y*o.x + p.z*o.w;
131 90 : w = p.w*o.w - p.x*o.x - p.y*o.y - p.z*o.z;
132 90 : return *this;
133 : }
134 9 : quat &mul(const quat &o) { return mul(quat(*this), o); }
135 :
136 : /**
137 : * @brief Turns this quaternion into its inverse.
138 : *
139 : * The inverse q⁻¹ of the quaternion (used in the rotation formula q*a*q⁻¹)
140 : * is one with the imaginary components equal to their additive inverse, with
141 : * no change to the real (scalar) component.
142 : *
143 : * @return a reference to `this` object
144 : */
145 26 : quat &invert() { neg3(); return *this; }
146 :
147 6 : quat &normalize() { vec4<float>::safenormalize(); return *this; }
148 :
149 503 : vec rotate(const vec &v) const
150 : {
151 503 : return vec().cross(*this, vec().cross(*this, v).madd(v, w)).mul(2).add(v);
152 : }
153 :
154 3 : vec invertedrotate(const vec &v) const
155 : {
156 3 : return vec().cross(*this, vec().cross(*this, v).msub(v, w)).mul(2).add(v);
157 : }
158 :
159 : private:
160 : template<class M>
161 446 : void convertmatrix(const M &m)
162 : {
163 446 : float trace = m.a.x + m.b.y + m.c.z;
164 446 : if(trace>0)
165 : {
166 294 : float r = sqrtf(1 + trace), inv = 0.5f/r;
167 294 : w = 0.5f*r;
168 294 : x = (m.b.z - m.c.y)*inv;
169 294 : y = (m.c.x - m.a.z)*inv;
170 294 : z = (m.a.y - m.b.x)*inv;
171 : }
172 152 : else if(m.a.x > m.b.y && m.a.x > m.c.z)
173 : {
174 25 : float r = sqrtf(1 + m.a.x - m.b.y - m.c.z), inv = 0.5f/r;
175 25 : x = 0.5f*r;
176 25 : y = (m.a.y + m.b.x)*inv;
177 25 : z = (m.c.x + m.a.z)*inv;
178 25 : w = (m.b.z - m.c.y)*inv;
179 25 : }
180 127 : else if(m.b.y > m.c.z)
181 : {
182 60 : float r = sqrtf(1 + m.b.y - m.a.x - m.c.z), inv = 0.5f/r;
183 60 : x = (m.a.y + m.b.x)*inv;
184 60 : y = 0.5f*r;
185 60 : z = (m.b.z + m.c.y)*inv;
186 60 : w = (m.c.x - m.a.z)*inv;
187 : }
188 : else
189 : {
190 67 : float r = sqrtf(1 + m.c.z - m.a.x - m.b.y), inv = 0.5f/r;
191 67 : x = (m.c.x + m.a.z)*inv;
192 67 : y = (m.b.z + m.c.y)*inv;
193 67 : z = 0.5f*r;
194 67 : w = (m.a.y - m.b.x)*inv;
195 : }
196 446 : }
197 : };
198 :
199 : /**
200 : * @brief dual quaternion object
201 : *
202 : * 8 dimensional numbers that are the product of dual numbers and quaternions
203 : * used for rigid body transformations (like animations) (dualquats have 6 DoF)
204 : */
205 : struct dualquat
206 : {
207 : quat real, dual;
208 :
209 35 : dualquat() {}
210 12 : dualquat(const quat &q, const vec &p)
211 12 : : real(q),
212 12 : dual(0.5f*( p.x*q.w + p.y*q.z - p.z*q.y),
213 12 : 0.5f*(-p.x*q.z + p.y*q.w + p.z*q.x),
214 12 : 0.5f*( p.x*q.y - p.y*q.x + p.z*q.w),
215 12 : -0.5f*( p.x*q.x + p.y*q.y + p.z*q.z))
216 : {
217 12 : }
218 31 : explicit dualquat(const quat &q) : real(q), dual(0, 0, 0, 0) {}
219 : explicit dualquat(const matrix4x3 &m);
220 :
221 4 : dualquat &mul(float k) { real.mul(k); dual.mul(k); return *this; }
222 :
223 7 : dualquat &invert()
224 : {
225 7 : real.invert();
226 7 : dual.invert();
227 7 : dual.msub(real, 2*real.dot(dual));
228 7 : return *this;
229 : }
230 :
231 20 : void mul(const dualquat &p, const dualquat &o)
232 : {
233 20 : real.mul(p.real, o.real);
234 20 : dual.mul(p.real, o.dual).add(quat().mul(p.dual, o.real));
235 20 : }
236 6 : void mul(const dualquat &o) { mul(dualquat(*this), o); }
237 :
238 3 : void mulorient(const quat &q)
239 : {
240 3 : real.mul(q, quat(real));
241 3 : dual.mul(quat(q).invert(), quat(dual));
242 3 : }
243 :
244 2 : void mulorient(const quat &q, const dualquat &base)
245 : {
246 2 : quat trans;
247 2 : trans.mul(base.dual, quat(base.real).invert());
248 2 : dual.mul(quat(q).invert(), quat().mul(real, trans).add(dual));
249 :
250 2 : real.mul(q, quat(real));
251 2 : dual.add(quat().mul(real, trans.invert())).msub(real, 2*base.real.dot(base.dual));
252 2 : }
253 :
254 2 : void normalize()
255 : {
256 2 : float invlen = 1/real.magnitude();
257 2 : real.mul(invlen);
258 2 : dual.mul(invlen);
259 2 : }
260 :
261 2 : void translate(const vec &p)
262 : {
263 2 : dual.x += 0.5f*( p.x*real.w + p.y*real.z - p.z*real.y);
264 2 : dual.y += 0.5f*(-p.x*real.z + p.y*real.w + p.z*real.x);
265 2 : dual.z += 0.5f*( p.x*real.y - p.y*real.x + p.z*real.w);
266 2 : dual.w += -0.5f*( p.x*real.x + p.y*real.y + p.z*real.z);
267 2 : }
268 :
269 6 : void fixantipodal(const dualquat &d)
270 : {
271 6 : if(real.dot(d.real) < 0)
272 : {
273 1 : real.neg();
274 1 : dual.neg();
275 : }
276 6 : }
277 :
278 7 : void accumulate(const dualquat &d, float k)
279 : {
280 7 : if(real.dot(d.real) < 0)
281 : {
282 1 : k = -k;
283 : }
284 7 : real.madd(d.real, k);
285 7 : dual.madd(d.dual, k);
286 7 : }
287 :
288 2 : vec transform(const vec &v) const
289 : {
290 2 : return vec().cross(real, vec().cross(real, v).madd(v, real.w).add(vec(dual))).madd(vec(dual), real.w).msub(vec(real), dual.w).mul(2).add(v);
291 : }
292 :
293 2 : quat transform(const quat &q) const
294 : {
295 2 : return quat().mul(real, q);
296 : }
297 :
298 2 : vec transformnormal(const vec &v) const
299 : {
300 2 : return real.rotate(v);
301 : }
302 :
303 : };
304 :
305 2 : inline dualquat::dualquat(const matrix4x3 &m) : real(m)
306 : {
307 2 : dual.x = 0.5f*( m.d.x*real.w + m.d.y*real.z - m.d.z*real.y);
308 2 : dual.y = 0.5f*(-m.d.x*real.z + m.d.y*real.w + m.d.z*real.x);
309 2 : dual.z = 0.5f*( m.d.x*real.y - m.d.y*real.x + m.d.z*real.w);
310 2 : dual.w = -0.5f*( m.d.x*real.x + m.d.y*real.y + m.d.z*real.z);
311 2 : }
312 :
313 : struct plane : vec
314 : {
315 : float offset;
316 :
317 2 : float dist(const vec &p) const { return dot(p)+offset; }
318 3 : float dist(const vec4<float> &p) const { return p.dot3(*this) + p.w*offset; }
319 15 : bool operator==(const plane &p) const { return x==p.x && y==p.y && z==p.z && offset==p.offset; }
320 : bool operator!=(const plane &p) const { return x!=p.x || y!=p.y || z!=p.z || offset!=p.offset; }
321 :
322 12332 : plane() : vec(0,0,0), offset(0)
323 : {
324 12332 : }
325 :
326 27 : plane(const vec &c, float off) : vec(c), offset(off)
327 : {
328 27 : if(x == 0 && y == 0 && z == 0)
329 : {
330 1 : throw std::invalid_argument("cannot create plane with no normal vector");
331 : }
332 26 : }
333 2 : plane(const vec4<float> &p) : vec(p), offset(p.w)
334 : {
335 2 : if(x == 0 && y == 0 && z == 0)
336 : {
337 1 : throw std::invalid_argument("cannot create plane with no normal vector");
338 : }
339 1 : }
340 5 : plane(int d, float off)
341 5 : {
342 5 : if(d < 0 || d > 2)
343 : {
344 2 : throw std::invalid_argument("cannot specify plane index outside 0..2");
345 : }
346 3 : x = y = z = 0.0f;
347 3 : (*this)[d] = 1.0f;
348 3 : offset = -off;
349 3 : }
350 22 : plane(float a, float b, float c, float d) : vec(a, b, c), offset(d)
351 : {
352 22 : if(x == 0 && y == 0 && z == 0)
353 : {
354 1 : throw std::invalid_argument("cannot create plane with no normal vector");
355 : }
356 21 : }
357 :
358 1 : void toplane(const vec &n, const vec &p)
359 : {
360 1 : x = n.x;
361 1 : y = n.y;
362 1 : z = n.z;
363 1 : offset = -dot(p);
364 1 : }
365 :
366 5 : bool toplane(const vec &a, const vec &b, const vec &c)
367 : {
368 5 : cross(vec(b).sub(a), vec(c).sub(a));
369 5 : float mag = magnitude();
370 5 : if(!mag)
371 : {
372 4 : return false;
373 : }
374 1 : div(mag);
375 1 : offset = -dot(a);
376 1 : return true;
377 : }
378 :
379 4 : bool rayintersect(const vec &o, const vec &ray, float &dist) const
380 : {
381 4 : float cosalpha = dot(ray);
382 4 : if(cosalpha==0)
383 : {
384 2 : return false;
385 : }
386 2 : float deltac = offset+dot(o);
387 2 : dist -= deltac/cosalpha;
388 2 : return true;
389 : }
390 :
391 2 : plane &reflectz(float rz)
392 : {
393 2 : offset += 2*rz*z;
394 2 : z = -z;
395 2 : return *this;
396 : }
397 :
398 2 : plane &invert()
399 : {
400 2 : neg();
401 2 : offset = -offset;
402 2 : return *this;
403 : }
404 :
405 2 : plane &scale(float k)
406 : {
407 2 : mul(k);
408 2 : return *this;
409 : }
410 :
411 : plane &translate(const vec &p)
412 : {
413 : offset += dot(p);
414 : return *this;
415 : }
416 :
417 2 : plane &normalize()
418 : {
419 2 : float mag = magnitude();
420 2 : div(mag);
421 2 : offset /= mag;
422 2 : return *this;
423 : }
424 :
425 5 : float zdelta(const vec &p) const
426 : {
427 5 : return -(x*p.x+y*p.y)/z;
428 : }
429 :
430 : };
431 :
432 : //short integer quaternion
433 : class squat
434 : {
435 : public:
436 : short x, y, z, w;
437 :
438 11 : squat() {}
439 : //all dimensions of `q` should be <= 1 (normalized)
440 6 : squat(const vec4<float> &q)
441 6 : {
442 6 : convert(q);
443 6 : }
444 :
445 3 : void lerp(const vec4<float> &a, const vec4<float> &b, float t)
446 : {
447 3 : vec4<float> q;
448 3 : q.lerp(a, b, t);
449 3 : convert(q);
450 3 : }
451 : private:
452 9 : void convert(const vec4<float> &q)
453 : {
454 9 : x = static_cast<short>(q.x*32767.5f-0.5f);
455 9 : y = static_cast<short>(q.y*32767.5f-0.5f);
456 9 : z = static_cast<short>(q.z*32767.5f-0.5f);
457 9 : w = static_cast<short>(q.w*32767.5f-0.5f);
458 9 : }
459 : };
460 :
461 : struct matrix2
462 : {
463 : vec2 a, b;
464 :
465 : matrix2() {}
466 : matrix2(const vec2 &a, const vec2 &b) : a(a), b(b) {}
467 : explicit matrix2(const matrix4 &m) : a(m.a), b(m.b) {}
468 : explicit matrix2(const matrix3 &m) : a(m.a), b(m.b) {}
469 : };
470 :
471 : #endif
|
