1 // This is an incomplete port of http://glmatrix.net/
2 //
3 // Copyright (c) 2013 Brandon Jones, Colin MacKenzie IV
4 //
5 // This software is provided 'as-is', without any express or implied warranty.
6 // In no event will the authors be held liable for any damages arising from the
7 // use of this software.
8 //
9 // Permission is granted to anyone to use this software for any purpose,
10 // including commercial applications, and to alter it and redistribute it
11 // freely, subject to the following restrictions:
12 //
13 // 1. The origin of this software must not be misrepresented; you must not claim
14 //    that you wrote the original software. If you use this software in a
15 //    product, an acknowledgment in the product documentation would be
16 //    appreciated but is not required.
17 //
18 // 2. Altered source versions must be plainly marked as such, and must not be
19 //    misrepresented as being the original software.
20 //
21 // 3. This notice may not be removed or altered from any source distribution.
22 
23 #include <mbgl/util/mat4.hpp>
24 
25 #include <cmath>
26 
27 namespace mbgl {
28 
29 namespace matrix {
30 
identity(mat4 & out)31 void identity(mat4& out) {
32     out[0] = 1.0f;
33     out[1] = 0.0f;
34     out[2] = 0.0f;
35     out[3] = 0.0f;
36     out[4] = 0.0f;
37     out[5] = 1.0f;
38     out[6] = 0.0f;
39     out[7] = 0.0f;
40     out[8] = 0.0f;
41     out[9] = 0.0f;
42     out[10] = 1.0f;
43     out[11] = 0.0f;
44     out[12] = 0.0f;
45     out[13] = 0.0f;
46     out[14] = 0.0f;
47     out[15] = 1.0f;
48 }
49 
invert(mat4 & out,mat4 & a)50 bool invert(mat4& out, mat4& a) {
51     double a00 = a[0], a01 = a[1], a02 = a[2], a03 = a[3],
52           a10 = a[4], a11 = a[5], a12 = a[6], a13 = a[7],
53           a20 = a[8], a21 = a[9], a22 = a[10], a23 = a[11],
54           a30 = a[12], a31 = a[13], a32 = a[14], a33 = a[15],
55 
56           b00 = a00 * a11 - a01 * a10,
57           b01 = a00 * a12 - a02 * a10,
58           b02 = a00 * a13 - a03 * a10,
59           b03 = a01 * a12 - a02 * a11,
60           b04 = a01 * a13 - a03 * a11,
61           b05 = a02 * a13 - a03 * a12,
62           b06 = a20 * a31 - a21 * a30,
63           b07 = a20 * a32 - a22 * a30,
64           b08 = a20 * a33 - a23 * a30,
65           b09 = a21 * a32 - a22 * a31,
66           b10 = a21 * a33 - a23 * a31,
67           b11 = a22 * a33 - a23 * a32,
68 
69           // Calculate the determinant
70           det = b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06;
71 
72     if (!det) {
73         return true;
74     }
75     det = 1.0 / det;
76 
77     out[0] = (a11 * b11 - a12 * b10 + a13 * b09) * det;
78     out[1] = (a02 * b10 - a01 * b11 - a03 * b09) * det;
79     out[2] = (a31 * b05 - a32 * b04 + a33 * b03) * det;
80     out[3] = (a22 * b04 - a21 * b05 - a23 * b03) * det;
81     out[4] = (a12 * b08 - a10 * b11 - a13 * b07) * det;
82     out[5] = (a00 * b11 - a02 * b08 + a03 * b07) * det;
83     out[6] = (a32 * b02 - a30 * b05 - a33 * b01) * det;
84     out[7] = (a20 * b05 - a22 * b02 + a23 * b01) * det;
85     out[8] = (a10 * b10 - a11 * b08 + a13 * b06) * det;
86     out[9] = (a01 * b08 - a00 * b10 - a03 * b06) * det;
87     out[10] = (a30 * b04 - a31 * b02 + a33 * b00) * det;
88     out[11] = (a21 * b02 - a20 * b04 - a23 * b00) * det;
89     out[12] = (a11 * b07 - a10 * b09 - a12 * b06) * det;
90     out[13] = (a00 * b09 - a01 * b07 + a02 * b06) * det;
91     out[14] = (a31 * b01 - a30 * b03 - a32 * b00) * det;
92     out[15] = (a20 * b03 - a21 * b01 + a22 * b00) * det;
93 
94     return false;
95 }
96 
ortho(mat4 & out,double left,double right,double bottom,double top,double near,double far)97 void ortho(mat4& out, double left, double right, double bottom, double top, double near, double far) {
98     double lr = 1.0f / (left - right),
99           bt = 1.0f / (bottom - top),
100           nf = 1.0f / (near - far);
101     out[0] = -2.0f * lr;
102     out[1] = 0.0f;
103     out[2] = 0.0f;
104     out[3] = 0.0f;
105     out[4] = 0.0f;
106     out[5] = -2.0f * bt;
107     out[6] = 0.0f;
108     out[7] = 0.0f;
109     out[8] = 0.0f;
110     out[9] = 0.0f;
111     out[10] = 2.0f * nf;
112     out[11] = 0.0f;
113     out[12] = (left + right) * lr;
114     out[13] = (top + bottom) * bt;
115     out[14] = (far + near) * nf;
116     out[15] = 1.0f;
117 }
118 
perspective(mat4 & out,double fovy,double aspect,double near,double far)119 void perspective(mat4& out, double fovy, double aspect, double near, double far) {
120     double f = 1.0f / std::tan(fovy / 2.0f),
121           nf = 1.0f / (near - far);
122     out[0] = f / aspect;
123     out[1] = 0.0f;
124     out[2] = 0.0f;
125     out[3] = 0.0f;
126     out[4] = 0.0f;
127     out[5] = f;
128     out[6] = 0.0f;
129     out[7] = 0.0f;
130     out[8] = 0.0f;
131     out[9] = 0.0f;
132     out[10] = (far + near) * nf;
133     out[11] = -1.0f;
134     out[12] = 0.0f;
135     out[13] = 0.0f;
136     out[14] = (2.0f * far * near) * nf;
137     out[15] = 0.0f;
138 }
139 
copy(mat4 & out,const mat4 & a)140 void copy(mat4& out, const mat4& a) {
141     out[0] = a[0];
142     out[1] = a[1];
143     out[2] = a[2];
144     out[3] = a[3];
145     out[4] = a[4];
146     out[5] = a[5];
147     out[6] = a[6];
148     out[7] = a[7];
149     out[8] = a[8];
150     out[9] = a[9];
151     out[10] = a[10];
152     out[11] = a[11];
153     out[12] = a[12];
154     out[13] = a[13];
155     out[14] = a[14];
156     out[15] = a[15];
157 }
158 
translate(mat4 & out,const mat4 & a,double x,double y,double z)159 void translate(mat4& out, const mat4& a, double x, double y, double z) {
160     if (&a == &out) {
161         out[12] = a[0] * x + a[4] * y + a[8] * z + a[12];
162         out[13] = a[1] * x + a[5] * y + a[9] * z + a[13];
163         out[14] = a[2] * x + a[6] * y + a[10] * z + a[14];
164         out[15] = a[3] * x + a[7] * y + a[11] * z + a[15];
165     } else {
166         double a00, a01, a02, a03,
167               a10, a11, a12, a13,
168               a20, a21, a22, a23;
169 
170         a00 = a[0]; a01 = a[1]; a02 = a[2]; a03 = a[3];
171         a10 = a[4]; a11 = a[5]; a12 = a[6]; a13 = a[7];
172         a20 = a[8]; a21 = a[9]; a22 = a[10]; a23 = a[11];
173 
174         out[0] = a00; out[1] = a01; out[2] = a02; out[3] = a03;
175         out[4] = a10; out[5] = a11; out[6] = a12; out[7] = a13;
176         out[8] = a20; out[9] = a21; out[10] = a22; out[11] = a23;
177 
178         out[12] = a00 * x + a10 * y + a20 * z + a[12];
179         out[13] = a01 * x + a11 * y + a21 * z + a[13];
180         out[14] = a02 * x + a12 * y + a22 * z + a[14];
181         out[15] = a03 * x + a13 * y + a23 * z + a[15];
182     }
183 }
184 
rotate_x(mat4 & out,const mat4 & a,double rad)185 void rotate_x(mat4& out, const mat4& a, double rad) {
186     double s = std::sin(rad),
187         c = std::cos(rad),
188         a10 = a[4],
189         a11 = a[5],
190         a12 = a[6],
191         a13 = a[7],
192         a20 = a[8],
193         a21 = a[9],
194         a22 = a[10],
195         a23 = a[11];
196 
197     if (&a != &out) { // If the source and destination differ, copy the unchanged rows
198         out[0]  = a[0];
199         out[1]  = a[1];
200         out[2]  = a[2];
201         out[3]  = a[3];
202         out[12] = a[12];
203         out[13] = a[13];
204         out[14] = a[14];
205         out[15] = a[15];
206     }
207 
208     // Perform axis-specific matrix multiplication
209     out[4] = a10 * c + a20 * s;
210     out[5] = a11 * c + a21 * s;
211     out[6] = a12 * c + a22 * s;
212     out[7] = a13 * c + a23 * s;
213     out[8] = a20 * c - a10 * s;
214     out[9] = a21 * c - a11 * s;
215     out[10] = a22 * c - a12 * s;
216     out[11] = a23 * c - a13 * s;
217 }
218 
rotate_y(mat4 & out,const mat4 & a,double rad)219 void rotate_y(mat4& out, const mat4& a, double rad) {
220     double s = std::sin(rad),
221         c = std::cos(rad),
222         a00 = a[0],
223         a01 = a[1],
224         a02 = a[2],
225         a03 = a[3],
226         a20 = a[8],
227         a21 = a[9],
228         a22 = a[10],
229         a23 = a[11];
230 
231     if (&a != &out) { // If the source and destination differ, copy the unchanged rows
232         out[4]  = a[4];
233         out[5]  = a[5];
234         out[6]  = a[6];
235         out[7]  = a[7];
236         out[12] = a[12];
237         out[13] = a[13];
238         out[14] = a[14];
239         out[15] = a[15];
240     }
241 
242     // Perform axis-specific matrix multiplication
243     out[0] = a00 * c - a20 * s;
244     out[1] = a01 * c - a21 * s;
245     out[2] = a02 * c - a22 * s;
246     out[3] = a03 * c - a23 * s;
247     out[8] = a00 * s + a20 * c;
248     out[9] = a01 * s + a21 * c;
249     out[10] = a02 * s + a22 * c;
250     out[11] = a03 * s + a23 * c;
251 }
252 
rotate_z(mat4 & out,const mat4 & a,double rad)253 void rotate_z(mat4& out, const mat4& a, double rad) {
254     double s = std::sin(rad),
255           c = std::cos(rad),
256           a00 = a[0],
257           a01 = a[1],
258           a02 = a[2],
259           a03 = a[3],
260           a10 = a[4],
261           a11 = a[5],
262           a12 = a[6],
263           a13 = a[7];
264 
265     if (&a != &out) { // If the source and destination differ, copy the unchanged last row
266         out[8]  = a[8];
267         out[9]  = a[9];
268         out[10] = a[10];
269         out[11] = a[11];
270         out[12] = a[12];
271         out[13] = a[13];
272         out[14] = a[14];
273         out[15] = a[15];
274     }
275 
276     // Perform axis-specific matrix multiplication
277     out[0] = a00 * c + a10 * s;
278     out[1] = a01 * c + a11 * s;
279     out[2] = a02 * c + a12 * s;
280     out[3] = a03 * c + a13 * s;
281     out[4] = a10 * c - a00 * s;
282     out[5] = a11 * c - a01 * s;
283     out[6] = a12 * c - a02 * s;
284     out[7] = a13 * c - a03 * s;
285 }
286 
scale(mat4 & out,const mat4 & a,double x,double y,double z)287 void scale(mat4& out, const mat4& a, double x, double y, double z) {
288     out[0] = a[0] * x;
289     out[1] = a[1] * x;
290     out[2] = a[2] * x;
291     out[3] = a[3] * x;
292     out[4] = a[4] * y;
293     out[5] = a[5] * y;
294     out[6] = a[6] * y;
295     out[7] = a[7] * y;
296     out[8] = a[8] * z;
297     out[9] = a[9] * z;
298     out[10] = a[10] * z;
299     out[11] = a[11] * z;
300     out[12] = a[12];
301     out[13] = a[13];
302     out[14] = a[14];
303     out[15] = a[15];
304 }
305 
multiply(mat4 & out,const mat4 & a,const mat4 & b)306 void multiply(mat4& out, const mat4& a, const mat4& b) {
307     double a00 = a[0], a01 = a[1], a02 = a[2], a03 = a[3],
308           a10 = a[4], a11 = a[5], a12 = a[6], a13 = a[7],
309           a20 = a[8], a21 = a[9], a22 = a[10], a23 = a[11],
310           a30 = a[12], a31 = a[13], a32 = a[14], a33 = a[15];
311 
312     // Cache only the current line of the second matrix
313     double b0 = b[0], b1 = b[1], b2 = b[2], b3 = b[3];
314     out[0] = b0 * a00 + b1 * a10 + b2 * a20 + b3 * a30;
315     out[1] = b0 * a01 + b1 * a11 + b2 * a21 + b3 * a31;
316     out[2] = b0 * a02 + b1 * a12 + b2 * a22 + b3 * a32;
317     out[3] = b0 * a03 + b1 * a13 + b2 * a23 + b3 * a33;
318 
319     b0 = b[4]; b1 = b[5]; b2 = b[6]; b3 = b[7];
320     out[4] = b0 * a00 + b1 * a10 + b2 * a20 + b3 * a30;
321     out[5] = b0 * a01 + b1 * a11 + b2 * a21 + b3 * a31;
322     out[6] = b0 * a02 + b1 * a12 + b2 * a22 + b3 * a32;
323     out[7] = b0 * a03 + b1 * a13 + b2 * a23 + b3 * a33;
324 
325     b0 = b[8]; b1 = b[9]; b2 = b[10]; b3 = b[11];
326     out[8] = b0 * a00 + b1 * a10 + b2 * a20 + b3 * a30;
327     out[9] = b0 * a01 + b1 * a11 + b2 * a21 + b3 * a31;
328     out[10] = b0 * a02 + b1 * a12 + b2 * a22 + b3 * a32;
329     out[11] = b0 * a03 + b1 * a13 + b2 * a23 + b3 * a33;
330 
331     b0 = b[12]; b1 = b[13]; b2 = b[14]; b3 = b[15];
332     out[12] = b0 * a00 + b1 * a10 + b2 * a20 + b3 * a30;
333     out[13] = b0 * a01 + b1 * a11 + b2 * a21 + b3 * a31;
334     out[14] = b0 * a02 + b1 * a12 + b2 * a22 + b3 * a32;
335     out[15] = b0 * a03 + b1 * a13 + b2 * a23 + b3 * a33;
336 }
337 
transformMat4(vec4 & out,const vec4 & a,const mat4 & m)338 void transformMat4(vec4& out, const vec4& a, const mat4& m) {
339     double x = a[0], y = a[1], z = a[2], w = a[3];
340     out[0] = m[0] * x + m[4] * y + m[8] * z + m[12] * w;
341     out[1] = m[1] * x + m[5] * y + m[9] * z + m[13] * w;
342     out[2] = m[2] * x + m[6] * y + m[10] * z + m[14] * w;
343     out[3] = m[3] * x + m[7] * y + m[11] * z + m[15] * w;
344 }
345 
346 } // namespace matrix
347 
348 } // namespace mbgl
349