1 // Copyright 2009-2021 Intel Corporation
2 // SPDX-License-Identifier: Apache-2.0
3 
4 #pragma once
5 
6 #include "vec2.h"
7 
8 namespace embree
9 {
10   ////////////////////////////////////////////////////////////////////////////////
11   /// 2D Linear Transform (2x2 Matrix)
12   ////////////////////////////////////////////////////////////////////////////////
13 
14   template<typename T> struct LinearSpace2
15   {
16     typedef T Vector;
17     typedef typename T::Scalar Scalar;
18 
19     /*! default matrix constructor */
LinearSpace2LinearSpace220     __forceinline LinearSpace2           ( ) {}
LinearSpace2LinearSpace221     __forceinline LinearSpace2           ( const LinearSpace2& other ) { vx = other.vx; vy = other.vy; }
22     __forceinline LinearSpace2& operator=( const LinearSpace2& other ) { vx = other.vx; vy = other.vy; return *this; }
23 
LinearSpace2LinearSpace224     template<typename L1> __forceinline LinearSpace2( const LinearSpace2<L1>& s ) : vx(s.vx), vy(s.vy) {}
25 
26     /*! matrix construction from column vectors */
LinearSpace2LinearSpace227     __forceinline LinearSpace2(const Vector& vx, const Vector& vy)
28       : vx(vx), vy(vy) {}
29 
30     /*! matrix construction from row mayor data */
LinearSpace2LinearSpace231     __forceinline LinearSpace2(const Scalar& m00, const Scalar& m01,
32                                const Scalar& m10, const Scalar& m11)
33       : vx(m00,m10), vy(m01,m11) {}
34 
35     /*! compute the determinant of the matrix */
detLinearSpace236     __forceinline const Scalar det() const { return vx.x*vy.y - vx.y*vy.x; }
37 
38     /*! compute adjoint matrix */
adjointLinearSpace239     __forceinline const LinearSpace2 adjoint() const { return LinearSpace2(vy.y,-vy.x,-vx.y,vx.x); }
40 
41     /*! compute inverse matrix */
inverseLinearSpace242     __forceinline const LinearSpace2 inverse() const { return adjoint()/det(); }
43 
44     /*! compute transposed matrix */
transposedLinearSpace245     __forceinline const LinearSpace2 transposed() const { return LinearSpace2(vx.x,vx.y,vy.x,vy.y); }
46 
47     /*! returns first row of matrix */
row0LinearSpace248     __forceinline Vector row0() const { return Vector(vx.x,vy.x); }
49 
50     /*! returns second row of matrix */
row1LinearSpace251     __forceinline Vector row1() const { return Vector(vx.y,vy.y); }
52 
53     ////////////////////////////////////////////////////////////////////////////////
54     /// Constants
55     ////////////////////////////////////////////////////////////////////////////////
56 
LinearSpace2LinearSpace257     __forceinline LinearSpace2( ZeroTy ) : vx(zero), vy(zero) {}
LinearSpace2LinearSpace258     __forceinline LinearSpace2( OneTy ) : vx(one, zero), vy(zero, one) {}
59 
60     /*! return matrix for scaling */
scaleLinearSpace261     static __forceinline LinearSpace2 scale(const Vector& s) {
62       return LinearSpace2(s.x,   0,
63                           0  , s.y);
64     }
65 
66     /*! return matrix for rotation */
rotateLinearSpace267     static __forceinline LinearSpace2 rotate(const Scalar& r) {
68       Scalar s = sin(r), c = cos(r);
69       return LinearSpace2(c, -s,
70                           s,  c);
71     }
72 
73     /*! return closest orthogonal matrix (i.e. a general rotation including reflection) */
orthogonalLinearSpace274     LinearSpace2 orthogonal() const
75     {
76       LinearSpace2 m = *this;
77 
78       // mirrored?
79       Scalar mirror(one);
80       if (m.det() < Scalar(zero)) {
81         m.vx = -m.vx;
82         mirror = -mirror;
83       }
84 
85       // rotation
86       for (int i = 0; i < 99; i++) {
87         const LinearSpace2 m_next = 0.5 * (m + m.transposed().inverse());
88         const LinearSpace2 d = m_next - m;
89         m = m_next;
90         // norm^2 of difference small enough?
91         if (max(dot(d.vx, d.vx), dot(d.vy, d.vy)) < 1e-8)
92           break;
93       }
94 
95       // rotation * mirror_x
96       return LinearSpace2(mirror*m.vx, m.vy);
97     }
98 
99   public:
100 
101     /*! the column vectors of the matrix */
102     Vector vx,vy;
103   };
104 
105   ////////////////////////////////////////////////////////////////////////////////
106   // Unary Operators
107   ////////////////////////////////////////////////////////////////////////////////
108 
109   template<typename T> __forceinline LinearSpace2<T> operator -( const LinearSpace2<T>& a ) { return LinearSpace2<T>(-a.vx,-a.vy); }
110   template<typename T> __forceinline LinearSpace2<T> operator +( const LinearSpace2<T>& a ) { return LinearSpace2<T>(+a.vx,+a.vy); }
rcp(const LinearSpace2<T> & a)111   template<typename T> __forceinline LinearSpace2<T> rcp       ( const LinearSpace2<T>& a ) { return a.inverse(); }
112 
113   ////////////////////////////////////////////////////////////////////////////////
114   // Binary Operators
115   ////////////////////////////////////////////////////////////////////////////////
116 
117   template<typename T> __forceinline LinearSpace2<T> operator +( const LinearSpace2<T>& a, const LinearSpace2<T>& b ) { return LinearSpace2<T>(a.vx+b.vx,a.vy+b.vy); }
118   template<typename T> __forceinline LinearSpace2<T> operator -( const LinearSpace2<T>& a, const LinearSpace2<T>& b ) { return LinearSpace2<T>(a.vx-b.vx,a.vy-b.vy); }
119 
120   template<typename T> __forceinline LinearSpace2<T> operator*(const typename T::Scalar & a, const LinearSpace2<T>& b) { return LinearSpace2<T>(a*b.vx, a*b.vy); }
121   template<typename T> __forceinline T               operator*(const LinearSpace2<T>& a, const T              & b) { return b.x*a.vx + b.y*a.vy; }
122   template<typename T> __forceinline LinearSpace2<T> operator*(const LinearSpace2<T>& a, const LinearSpace2<T>& b) { return LinearSpace2<T>(a*b.vx, a*b.vy); }
123 
124   template<typename T> __forceinline LinearSpace2<T> operator/(const LinearSpace2<T>& a, const typename T::Scalar & b) { return LinearSpace2<T>(a.vx/b, a.vy/b); }
125   template<typename T> __forceinline LinearSpace2<T> operator/(const LinearSpace2<T>& a, const LinearSpace2<T>& b) { return a * rcp(b); }
126 
127   template<typename T> __forceinline LinearSpace2<T>& operator *=( LinearSpace2<T>& a, const LinearSpace2<T>& b ) { return a = a * b; }
128   template<typename T> __forceinline LinearSpace2<T>& operator /=( LinearSpace2<T>& a, const LinearSpace2<T>& b ) { return a = a / b; }
129 
130   ////////////////////////////////////////////////////////////////////////////////
131   /// Comparison Operators
132   ////////////////////////////////////////////////////////////////////////////////
133 
134   template<typename T> __forceinline bool operator ==( const LinearSpace2<T>& a, const LinearSpace2<T>& b ) { return a.vx == b.vx && a.vy == b.vy; }
135   template<typename T> __forceinline bool operator !=( const LinearSpace2<T>& a, const LinearSpace2<T>& b ) { return a.vx != b.vx || a.vy != b.vy; }
136 
137   ////////////////////////////////////////////////////////////////////////////////
138   /// Output Operators
139   ////////////////////////////////////////////////////////////////////////////////
140 
141   template<typename T> static embree_ostream operator<<(embree_ostream cout, const LinearSpace2<T>& m) {
142     return cout << "{ vx = " << m.vx << ", vy = " << m.vy << "}";
143   }
144 
145   /*! Shortcuts for common linear spaces. */
146   typedef LinearSpace2<Vec2f> LinearSpace2f;
147   typedef LinearSpace2<Vec2fa> LinearSpace2fa;
148 }
149