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 &a0, const vec &b0, const vec &c0) : a(a0), b(b0), c(c0) {}
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 x0, float y0, float z0, float w0) : vec4<float>(x0, y0, z0, w0) {}
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)
222 : {
223 4 : real.mul(k);
224 4 : dual.mul(k);
225 4 : return *this;
226 : }
227 :
228 7 : dualquat &invert()
229 : {
230 7 : real.invert();
231 7 : dual.invert();
232 7 : dual.msub(real, 2*real.dot(dual));
233 7 : return *this;
234 : }
235 :
236 20 : void mul(const dualquat &p, const dualquat &o)
237 : {
238 20 : real.mul(p.real, o.real);
239 20 : dual.mul(p.real, o.dual).add(quat().mul(p.dual, o.real));
240 20 : }
241 :
242 6 : void mul(const dualquat &o)
243 : {
244 6 : mul(dualquat(*this), o);
245 6 : }
246 :
247 3 : void mulorient(const quat &q)
248 : {
249 3 : real.mul(q, quat(real));
250 3 : dual.mul(quat(q).invert(), quat(dual));
251 3 : }
252 :
253 2 : void mulorient(const quat &q, const dualquat &base)
254 : {
255 2 : quat trans;
256 2 : trans.mul(base.dual, quat(base.real).invert());
257 2 : dual.mul(quat(q).invert(), quat().mul(real, trans).add(dual));
258 :
259 2 : real.mul(q, quat(real));
260 2 : dual.add(quat().mul(real, trans.invert())).msub(real, 2*base.real.dot(base.dual));
261 2 : }
262 :
263 2 : void normalize()
264 : {
265 2 : float invlen = 1/real.magnitude();
266 2 : real.mul(invlen);
267 2 : dual.mul(invlen);
268 2 : }
269 :
270 2 : void translate(const vec &p)
271 : {
272 2 : dual.x += 0.5f*( p.x*real.w + p.y*real.z - p.z*real.y);
273 2 : dual.y += 0.5f*(-p.x*real.z + p.y*real.w + p.z*real.x);
274 2 : dual.z += 0.5f*( p.x*real.y - p.y*real.x + p.z*real.w);
275 2 : dual.w += -0.5f*( p.x*real.x + p.y*real.y + p.z*real.z);
276 2 : }
277 :
278 6 : void fixantipodal(const dualquat &d)
279 : {
280 6 : if(real.dot(d.real) < 0)
281 : {
282 1 : real.neg();
283 1 : dual.neg();
284 : }
285 6 : }
286 :
287 7 : void accumulate(const dualquat &d, float k)
288 : {
289 7 : if(real.dot(d.real) < 0)
290 : {
291 1 : k = -k;
292 : }
293 7 : real.madd(d.real, k);
294 7 : dual.madd(d.dual, k);
295 7 : }
296 :
297 2 : vec transform(const vec &v) const
298 : {
299 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);
300 : }
301 :
302 2 : quat transform(const quat &q) const
303 : {
304 2 : return quat().mul(real, q);
305 : }
306 :
307 2 : vec transformnormal(const vec &v) const
308 : {
309 2 : return real.rotate(v);
310 : }
311 :
312 : };
313 :
314 2 : inline dualquat::dualquat(const matrix4x3 &m) : real(m)
315 : {
316 2 : dual.x = 0.5f*( m.d.x*real.w + m.d.y*real.z - m.d.z*real.y);
317 2 : dual.y = 0.5f*(-m.d.x*real.z + m.d.y*real.w + m.d.z*real.x);
318 2 : dual.z = 0.5f*( m.d.x*real.y - m.d.y*real.x + m.d.z*real.w);
319 2 : dual.w = -0.5f*( m.d.x*real.x + m.d.y*real.y + m.d.z*real.z);
320 2 : }
321 :
322 : struct plane : vec
323 : {
324 : float offset;
325 :
326 2 : float dist(const vec &p) const { return dot(p)+offset; }
327 3 : float dist(const vec4<float> &p) const { return p.dot3(*this) + p.w*offset; }
328 15 : bool operator==(const plane &p) const { return x==p.x && y==p.y && z==p.z && offset==p.offset; }
329 : bool operator!=(const plane &p) const { return x!=p.x || y!=p.y || z!=p.z || offset!=p.offset; }
330 :
331 12440 : plane() : vec(0,0,0), offset(0)
332 : {
333 12440 : }
334 :
335 27 : plane(const vec &c, float off) : vec(c), offset(off)
336 : {
337 27 : if(x == 0 && y == 0 && z == 0)
338 : {
339 1 : throw std::invalid_argument("cannot create plane with no normal vector");
340 : }
341 26 : }
342 :
343 2 : plane(const vec4<float> &p) : vec(p), offset(p.w)
344 : {
345 2 : if(x == 0 && y == 0 && z == 0)
346 : {
347 1 : throw std::invalid_argument("cannot create plane with no normal vector");
348 : }
349 1 : }
350 :
351 5 : plane(int d, float off)
352 5 : {
353 5 : if(d < 0 || d > 2)
354 : {
355 2 : throw std::invalid_argument("cannot specify plane index outside 0..2");
356 : }
357 3 : x = y = z = 0.0f;
358 3 : (*this)[d] = 1.0f;
359 3 : offset = -off;
360 3 : }
361 :
362 22 : plane(float a, float b, float c, float d) : vec(a, b, c), offset(d)
363 : {
364 22 : if(x == 0 && y == 0 && z == 0)
365 : {
366 1 : throw std::invalid_argument("cannot create plane with no normal vector");
367 : }
368 21 : }
369 :
370 1 : void toplane(const vec &n, const vec &p)
371 : {
372 1 : x = n.x;
373 1 : y = n.y;
374 1 : z = n.z;
375 1 : offset = -dot(p);
376 1 : }
377 :
378 5 : bool toplane(const vec &a, const vec &b, const vec &c)
379 : {
380 5 : cross(vec(b).sub(a), vec(c).sub(a));
381 5 : float mag = magnitude();
382 5 : if(!mag)
383 : {
384 4 : return false;
385 : }
386 1 : div(mag);
387 1 : offset = -dot(a);
388 1 : return true;
389 : }
390 :
391 4 : bool rayintersect(const vec &o, const vec &ray, float &dist) const
392 : {
393 4 : float cosalpha = dot(ray);
394 4 : if(cosalpha==0)
395 : {
396 2 : return false;
397 : }
398 2 : float deltac = offset+dot(o);
399 2 : dist -= deltac/cosalpha;
400 2 : return true;
401 : }
402 :
403 2 : plane &reflectz(float rz)
404 : {
405 2 : offset += 2*rz*z;
406 2 : z = -z;
407 2 : return *this;
408 : }
409 :
410 2 : plane &invert()
411 : {
412 2 : neg();
413 2 : offset = -offset;
414 2 : return *this;
415 : }
416 :
417 2 : plane &scale(float k)
418 : {
419 2 : mul(k);
420 2 : return *this;
421 : }
422 :
423 : plane &translate(const vec &p)
424 : {
425 : offset += dot(p);
426 : return *this;
427 : }
428 :
429 2 : plane &normalize()
430 : {
431 2 : float mag = magnitude();
432 2 : div(mag);
433 2 : offset /= mag;
434 2 : return *this;
435 : }
436 :
437 5 : float zdelta(const vec &p) const
438 : {
439 5 : return -(x*p.x+y*p.y)/z;
440 : }
441 :
442 : };
443 :
444 : //short integer quaternion
445 : class squat
446 : {
447 : public:
448 : short x, y, z, w;
449 :
450 11 : squat() {}
451 : //all dimensions of `q` should be <= 1 (normalized)
452 6 : squat(const vec4<float> &q)
453 6 : {
454 6 : convert(q);
455 6 : }
456 :
457 3 : void lerp(const vec4<float> &a, const vec4<float> &b, float t)
458 : {
459 3 : vec4<float> q;
460 3 : q.lerp(a, b, t);
461 3 : convert(q);
462 3 : }
463 : private:
464 9 : void convert(const vec4<float> &q)
465 : {
466 9 : x = static_cast<short>(q.x*32767.5f-0.5f);
467 9 : y = static_cast<short>(q.y*32767.5f-0.5f);
468 9 : z = static_cast<short>(q.z*32767.5f-0.5f);
469 9 : w = static_cast<short>(q.w*32767.5f-0.5f);
470 9 : }
471 : };
472 :
473 : struct matrix2
474 : {
475 : vec2 a, b;
476 :
477 : matrix2() {}
478 : matrix2(const vec2 &a0, const vec2 &b0) : a(a0), b(b0) {}
479 : explicit matrix2(const matrix4 &m) : a(m.a), b(m.b) {}
480 : explicit matrix2(const matrix3 &m) : a(m.a), b(m.b) {}
481 : };
482 :
483 : #endif
|
