LCOV - code coverage report
Current view: top level - shared - geomexts.h (source / functions) Coverage Total Hit
Test: Libprimis Test Coverage Lines: 99.2 % 249 247
Test Date: 2025-02-18 06:21:28 Functions: 100.0 % 69 69

            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
        

Generated by: LCOV version 2.0-1