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 */
LinearSpace3LinearSpace321     __forceinline LinearSpace3           ( ) {}
LinearSpace3LinearSpace322     __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 
LinearSpace3LinearSpace325     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 */
LinearSpace3LinearSpace328     __forceinline LinearSpace3(const Vector& vx, const Vector& vy, const Vector& vz)
29       : vx(vx), vy(vy), vz(vz) {}
30 
31     /*! construction from quaternion */
LinearSpace3LinearSpace332     __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 */
LinearSpace3LinearSpace338     __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 */
detLinearSpace344     __forceinline const Scalar det() const { return dot(vx,cross(vy,vz)); }
45 
46     /*! compute adjoint matrix */
adjointLinearSpace347     __forceinline const LinearSpace3 adjoint() const { return LinearSpace3(cross(vy,vz),cross(vz,vx),cross(vx,vy)).transposed(); }
48 
49     /*! compute inverse matrix */
inverseLinearSpace350     __forceinline const LinearSpace3 inverse() const { return adjoint()/det(); }
51 
52     /*! compute transposed matrix */
transposedLinearSpace353     __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 */
row0LinearSpace356     __forceinline Vector row0() const { return Vector(vx.x,vy.x,vz.x); }
57 
58     /*! returns second row of matrix */
row1LinearSpace359     __forceinline Vector row1() const { return Vector(vx.y,vy.y,vz.y); }
60 
61     /*! returns third row of matrix */
row2LinearSpace362     __forceinline Vector row2() const { return Vector(vx.z,vy.z,vz.z); }
63 
64     ////////////////////////////////////////////////////////////////////////////////
65     /// Constants
66     ////////////////////////////////////////////////////////////////////////////////
67 
LinearSpace3LinearSpace368     __forceinline LinearSpace3( ZeroTy ) : vx(zero), vy(zero), vz(zero) {}
LinearSpace3LinearSpace369     __forceinline LinearSpace3( OneTy ) : vx(one, zero, zero), vy(zero, one, zero), vz(zero, zero, one) {}
70 
71     /*! return matrix for scaling */
scaleLinearSpace372     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 */
rotateLinearSpace379     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 */
transposed()94   template<> __forceinline const LinearSpace3<Vec3fa> LinearSpace3<Vec3fa>::transposed() const {
95     vfloat4 rx,ry,rz; transpose((vfloat4&)vx,(vfloat4&)vy,(vfloat4&)vz,vfloat4(zero),rx,ry,rz);
96     return LinearSpace3<Vec3fa>(Vec3fa(rx),Vec3fa(ry),Vec3fa(rz));
97   }
98 
99   template<typename T>
transposed(const LinearSpace3<T> & xfm)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); }
rcp(const LinearSpace3<T> & a)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 */
frame(const T & N)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 */
frame(const T & N,const T & dxi)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 */
clamp(const LinearSpace3<T> & space)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 
xfmPoint(const LinearSpace3<T> & s,const T & a)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)); }
xfmVector(const LinearSpace3<T> & s,const T & a)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)); }
xfmNormal(const LinearSpace3<T> & s,const T & a)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 
select(const typename T::Scalar::Bool & s,const LinearSpace3<T> & t,const LinearSpace3<T> & f)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>
lerp(const LinearSpace3<T> & l0,const LinearSpace3<T> & l1,const float 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>
lerp(const LinearSpace3<T> & l0,const LinearSpace3<T> & l1,const S & t)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