| 1 | // Copyright 2009-2021 Intel Corporation | 
| 2 | // SPDX-License-Identifier: Apache-2.0 | 
| 3 |  | 
| 4 | #pragma once | 
| 5 |  | 
| 6 | #include "vec3.h" | 
| 7 | #include "quaternion.h" | 
| 8 |  | 
| 9 | namespace embree | 
| 10 | { | 
| 11 |   //////////////////////////////////////////////////////////////////////////////// | 
| 12 |   /// 3D Linear Transform (3x3 Matrix) | 
| 13 |   //////////////////////////////////////////////////////////////////////////////// | 
| 14 |  | 
| 15 |   template<typename T> struct LinearSpace3 | 
| 16 |   { | 
| 17 |     typedef T Vector; | 
| 18 |     typedef typename T::Scalar Scalar; | 
| 19 |  | 
| 20 |     /*! default matrix constructor */ | 
| 21 |     __forceinline LinearSpace3           ( ) {} | 
| 22 |     __forceinline LinearSpace3           ( const LinearSpace3& other ) { vx = other.vx; vy = other.vy; vz = other.vz; } | 
| 23 |     __forceinline LinearSpace3& operator=( const LinearSpace3& other ) { vx = other.vx; vy = other.vy; vz = other.vz; return *this; } | 
| 24 |  | 
| 25 |     template<typename L1> __forceinline LinearSpace3( const LinearSpace3<L1>& s ) : vx(s.vx), vy(s.vy), vz(s.vz) {} | 
| 26 |  | 
| 27 |     /*! matrix construction from column vectors */ | 
| 28 |     __forceinline LinearSpace3(const Vector& vx, const Vector& vy, const Vector& vz) | 
| 29 |       : vx(vx), vy(vy), vz(vz) {} | 
| 30 |  | 
| 31 |     /*! construction from quaternion */ | 
| 32 |     __forceinline LinearSpace3( const QuaternionT<Scalar>& q ) | 
| 33 |       : vx((q.r*q.r + q.i*q.i - q.j*q.j - q.k*q.k), 2.0f*(q.i*q.j + q.r*q.k), 2.0f*(q.i*q.k - q.r*q.j)) | 
| 34 |       , vy(2.0f*(q.i*q.j - q.r*q.k), (q.r*q.r - q.i*q.i + q.j*q.j - q.k*q.k), 2.0f*(q.j*q.k + q.r*q.i)) | 
| 35 |       , vz(2.0f*(q.i*q.k + q.r*q.j), 2.0f*(q.j*q.k - q.r*q.i), (q.r*q.r - q.i*q.i - q.j*q.j + q.k*q.k)) {} | 
| 36 |  | 
| 37 |     /*! matrix construction from row mayor data */ | 
| 38 |     __forceinline LinearSpace3(const Scalar& m00, const Scalar& m01, const Scalar& m02, | 
| 39 |                                const Scalar& m10, const Scalar& m11, const Scalar& m12, | 
| 40 |                                const Scalar& m20, const Scalar& m21, const Scalar& m22) | 
| 41 |       : vx(m00,m10,m20), vy(m01,m11,m21), vz(m02,m12,m22) {} | 
| 42 |  | 
| 43 |     /*! compute the determinant of the matrix */ | 
| 44 |     __forceinline const Scalar det() const { return dot(vx,cross(vy,vz)); } | 
| 45 |  | 
| 46 |     /*! compute adjoint matrix */ | 
| 47 |     __forceinline const LinearSpace3 adjoint() const { return LinearSpace3(cross(vy,vz),cross(vz,vx),cross(vx,vy)).transposed(); } | 
| 48 |  | 
| 49 |     /*! compute inverse matrix */ | 
| 50 |     __forceinline const LinearSpace3 inverse() const { return adjoint()/det(); } | 
| 51 |  | 
| 52 |     /*! compute transposed matrix */ | 
| 53 |     __forceinline const LinearSpace3 transposed() const { return LinearSpace3(vx.x,vx.y,vx.z,vy.x,vy.y,vy.z,vz.x,vz.y,vz.z); } | 
| 54 |  | 
| 55 |     /*! returns first row of matrix */ | 
| 56 |     __forceinline Vector row0() const { return Vector(vx.x,vy.x,vz.x); } | 
| 57 |  | 
| 58 |     /*! returns second row of matrix */ | 
| 59 |     __forceinline Vector row1() const { return Vector(vx.y,vy.y,vz.y); } | 
| 60 |  | 
| 61 |     /*! returns third row of matrix */ | 
| 62 |     __forceinline Vector row2() const { return Vector(vx.z,vy.z,vz.z); } | 
| 63 |  | 
| 64 |     //////////////////////////////////////////////////////////////////////////////// | 
| 65 |     /// Constants | 
| 66 |     //////////////////////////////////////////////////////////////////////////////// | 
| 67 |  | 
| 68 |     __forceinline LinearSpace3( ZeroTy ) : vx(zero), vy(zero), vz(zero) {} | 
| 69 |     __forceinline LinearSpace3( OneTy ) : vx(one, zero, zero), vy(zero, one, zero), vz(zero, zero, one) {} | 
| 70 |  | 
| 71 |     /*! return matrix for scaling */ | 
| 72 |     static __forceinline LinearSpace3 scale(const Vector& s) { | 
| 73 |       return LinearSpace3(s.x,   0,   0, | 
| 74 |                           0  , s.y,   0, | 
| 75 |                           0  ,   0, s.z); | 
| 76 |     } | 
| 77 |  | 
| 78 |     /*! return matrix for rotation around arbitrary axis */ | 
| 79 |     static __forceinline LinearSpace3 rotate(const Vector& _u, const Scalar& r) { | 
| 80 |       Vector u = normalize(_u); | 
| 81 |       Scalar s = sin(r), c = cos(r); | 
| 82 |       return LinearSpace3(u.x*u.x+(1-u.x*u.x)*c,  u.x*u.y*(1-c)-u.z*s,    u.x*u.z*(1-c)+u.y*s, | 
| 83 |                           u.x*u.y*(1-c)+u.z*s,    u.y*u.y+(1-u.y*u.y)*c,  u.y*u.z*(1-c)-u.x*s, | 
| 84 |                           u.x*u.z*(1-c)-u.y*s,    u.y*u.z*(1-c)+u.x*s,    u.z*u.z+(1-u.z*u.z)*c); | 
| 85 |     } | 
| 86 |  | 
| 87 |   public: | 
| 88 |  | 
| 89 |     /*! the column vectors of the matrix */ | 
| 90 |     Vector vx,vy,vz; | 
| 91 |   }; | 
| 92 |  | 
| 93 |   /*! compute transposed matrix */ | 
| 94 |   template<> __forceinline const LinearSpace3<Vec3fa> LinearSpace3<Vec3fa>::transposed() const {  | 
| 95 |     vfloat4 rx,ry,rz; transpose(r0: (vfloat4&)vx,r1: (vfloat4&)vy,r2: (vfloat4&)vz,r3: vfloat4(zero),c0&: rx,c1&: ry,c2&: rz); | 
| 96 |     return LinearSpace3<Vec3fa>(Vec3fa(rx),Vec3fa(ry),Vec3fa(rz));  | 
| 97 |   } | 
| 98 |  | 
| 99 |   template<typename T> | 
| 100 |     __forceinline const LinearSpace3<T> transposed(const LinearSpace3<T>& xfm) {  | 
| 101 |     return xfm.transposed(); | 
| 102 |   } | 
| 103 |  | 
| 104 |   //////////////////////////////////////////////////////////////////////////////// | 
| 105 |   // Unary Operators | 
| 106 |   //////////////////////////////////////////////////////////////////////////////// | 
| 107 |  | 
| 108 |   template<typename T> __forceinline LinearSpace3<T> operator -( const LinearSpace3<T>& a ) { return LinearSpace3<T>(-a.vx,-a.vy,-a.vz); } | 
| 109 |   template<typename T> __forceinline LinearSpace3<T> operator +( const LinearSpace3<T>& a ) { return LinearSpace3<T>(+a.vx,+a.vy,+a.vz); } | 
| 110 |   template<typename T> __forceinline LinearSpace3<T> rcp       ( const LinearSpace3<T>& a ) { return a.inverse(); } | 
| 111 |  | 
| 112 |   /* constructs a coordinate frame form a normalized normal */ | 
| 113 |   template<typename T> __forceinline LinearSpace3<T> frame(const T& N)  | 
| 114 |   { | 
| 115 |     const T dx0(0,N.z,-N.y); | 
| 116 |     const T dx1(-N.z,0,N.x); | 
| 117 |     const T dx = normalize(select(dot(dx0,dx0) > dot(dx1,dx1),dx0,dx1)); | 
| 118 |     const T dy = normalize(cross(N,dx)); | 
| 119 |     return LinearSpace3<T>(dx,dy,N); | 
| 120 |   } | 
| 121 |  | 
| 122 |   /* constructs a coordinate frame from a normal and approximate x-direction */ | 
| 123 |   template<typename T> __forceinline LinearSpace3<T> frame(const T& N, const T& dxi) | 
| 124 |   { | 
| 125 |     if (abs(dot(dxi,N)) > 0.99f) return frame(N); // fallback in case N and dxi are very parallel | 
| 126 |     const T dx = normalize(cross(dxi,N)); | 
| 127 |     const T dy = normalize(cross(N,dx)); | 
| 128 |     return LinearSpace3<T>(dx,dy,N); | 
| 129 |   } | 
| 130 |    | 
| 131 |   /* clamps linear space to range -1 to +1 */ | 
| 132 |   template<typename T> __forceinline LinearSpace3<T> clamp(const LinearSpace3<T>& space) { | 
| 133 |     return LinearSpace3<T>(clamp(space.vx,T(-1.0f),T(1.0f)), | 
| 134 |                            clamp(space.vy,T(-1.0f),T(1.0f)), | 
| 135 |                            clamp(space.vz,T(-1.0f),T(1.0f))); | 
| 136 |   } | 
| 137 |  | 
| 138 |   //////////////////////////////////////////////////////////////////////////////// | 
| 139 |   // Binary Operators | 
| 140 |   //////////////////////////////////////////////////////////////////////////////// | 
| 141 |  | 
| 142 |   template<typename T> __forceinline LinearSpace3<T> operator +( const LinearSpace3<T>& a, const LinearSpace3<T>& b ) { return LinearSpace3<T>(a.vx+b.vx,a.vy+b.vy,a.vz+b.vz); } | 
| 143 |   template<typename T> __forceinline LinearSpace3<T> operator -( const LinearSpace3<T>& a, const LinearSpace3<T>& b ) { return LinearSpace3<T>(a.vx-b.vx,a.vy-b.vy,a.vz-b.vz); } | 
| 144 |  | 
| 145 |   template<typename T> __forceinline LinearSpace3<T> operator*(const typename T::Scalar & a, const LinearSpace3<T>& b) { return LinearSpace3<T>(a*b.vx, a*b.vy, a*b.vz); } | 
| 146 |   template<typename T> __forceinline T               operator*(const LinearSpace3<T>& a, const T              & b) { return madd(T(b.x),a.vx,madd(T(b.y),a.vy,T(b.z)*a.vz)); } | 
| 147 |   template<typename T> __forceinline LinearSpace3<T> operator*(const LinearSpace3<T>& a, const LinearSpace3<T>& b) { return LinearSpace3<T>(a*b.vx, a*b.vy, a*b.vz); } | 
| 148 |  | 
| 149 |   template<typename T> __forceinline LinearSpace3<T> operator/(const LinearSpace3<T>& a, const typename T::Scalar & b) { return LinearSpace3<T>(a.vx/b, a.vy/b, a.vz/b); } | 
| 150 |   template<typename T> __forceinline LinearSpace3<T> operator/(const LinearSpace3<T>& a, const LinearSpace3<T>& b) { return a * rcp(b); } | 
| 151 |  | 
| 152 |   template<typename T> __forceinline LinearSpace3<T>& operator *=( LinearSpace3<T>& a, const LinearSpace3<T>& b ) { return a = a * b; } | 
| 153 |   template<typename T> __forceinline LinearSpace3<T>& operator /=( LinearSpace3<T>& a, const LinearSpace3<T>& b ) { return a = a / b; } | 
| 154 |  | 
| 155 |   template<typename T> __forceinline T       xfmPoint (const LinearSpace3<T>& s, const T      & a) { return madd(T(a.x),s.vx,madd(T(a.y),s.vy,T(a.z)*s.vz)); } | 
| 156 |   template<typename T> __forceinline T       xfmVector(const LinearSpace3<T>& s, const T      & a) { return madd(T(a.x),s.vx,madd(T(a.y),s.vy,T(a.z)*s.vz)); } | 
| 157 |   template<typename T> __forceinline T       xfmNormal(const LinearSpace3<T>& s, const T      & a) { return xfmVector(s.inverse().transposed(),a); } | 
| 158 |  | 
| 159 |   //////////////////////////////////////////////////////////////////////////////// | 
| 160 |   /// Comparison Operators | 
| 161 |   //////////////////////////////////////////////////////////////////////////////// | 
| 162 |  | 
| 163 |   template<typename T> __forceinline bool operator ==( const LinearSpace3<T>& a, const LinearSpace3<T>& b ) { return a.vx == b.vx && a.vy == b.vy && a.vz == b.vz; } | 
| 164 |   template<typename T> __forceinline bool operator !=( const LinearSpace3<T>& a, const LinearSpace3<T>& b ) { return a.vx != b.vx || a.vy != b.vy || a.vz != b.vz; } | 
| 165 |  | 
| 166 |   //////////////////////////////////////////////////////////////////////////////// | 
| 167 |   /// Select | 
| 168 |   //////////////////////////////////////////////////////////////////////////////// | 
| 169 |  | 
| 170 |   template<typename T> __forceinline LinearSpace3<T> select ( const typename T::Scalar::Bool& s, const LinearSpace3<T>& t, const LinearSpace3<T>& f ) { | 
| 171 |     return LinearSpace3<T>(select(s,t.vx,f.vx),select(s,t.vy,f.vy),select(s,t.vz,f.vz)); | 
| 172 |   } | 
| 173 |  | 
| 174 |   /*! blending */ | 
| 175 |   template<typename T> | 
| 176 |     __forceinline LinearSpace3<T> lerp(const LinearSpace3<T>& l0, const LinearSpace3<T>& l1, const float t)  | 
| 177 |   { | 
| 178 |     return LinearSpace3<T>(lerp(l0.vx,l1.vx,t), | 
| 179 |                            lerp(l0.vy,l1.vy,t), | 
| 180 |                            lerp(l0.vz,l1.vz,t)); | 
| 181 |   } | 
| 182 |  | 
| 183 |   //////////////////////////////////////////////////////////////////////////////// | 
| 184 |   /// Output Operators | 
| 185 |   //////////////////////////////////////////////////////////////////////////////// | 
| 186 |  | 
| 187 |   template<typename T> static embree_ostream operator<<(embree_ostream cout, const LinearSpace3<T>& m) { | 
| 188 |     return cout << "{ vx = "  << m.vx << ", vy = "  << m.vy << ", vz = "  << m.vz << "}" ; | 
| 189 |   } | 
| 190 |  | 
| 191 |   /*! Shortcuts for common linear spaces. */ | 
| 192 |   typedef LinearSpace3<Vec3f> LinearSpace3f; | 
| 193 |   typedef LinearSpace3<Vec3fa> LinearSpace3fa; | 
| 194 |   typedef LinearSpace3<Vec3fx> LinearSpace3fx; | 
| 195 |   typedef LinearSpace3<Vec3ff> LinearSpace3ff; | 
| 196 |  | 
| 197 |   template<int N> using LinearSpace3vf = LinearSpace3<Vec3<vfloat<N>>>; | 
| 198 |   typedef LinearSpace3<Vec3<vfloat<4>>>  LinearSpace3vf4; | 
| 199 |   typedef LinearSpace3<Vec3<vfloat<8>>>  LinearSpace3vf8; | 
| 200 |   typedef LinearSpace3<Vec3<vfloat<16>>> LinearSpace3vf16; | 
| 201 |  | 
| 202 |   /*! blending */ | 
| 203 |   template<typename T, typename S> | 
| 204 |     __forceinline LinearSpace3<T> lerp(const LinearSpace3<T>& l0, | 
| 205 |                                        const LinearSpace3<T>& l1, | 
| 206 |                                        const S& t) | 
| 207 |   { | 
| 208 |     return LinearSpace3<T>(lerp(l0.vx,l1.vx,t), | 
| 209 |                            lerp(l0.vy,l1.vy,t), | 
| 210 |                            lerp(l0.vz,l1.vz,t)); | 
| 211 |   } | 
| 212 |  | 
| 213 | } | 
| 214 |  |