1/*
2 * Copyright 2020 Google Inc.
3 *
4 * Use of this source code is governed by a BSD-style license that can be
5 * found in the LICENSE file.
6 */
7
8#include "include/core/SkM44.h"
9#include "include/core/SkMatrix.h"
10#include "include/private/SkVx.h"
11
12typedef skvx::Vec<4, float> sk4f;
13
14bool SkM44::operator==(const SkM44& other) const {
15 if (this == &other) {
16 return true;
17 }
18
19 sk4f a0 = sk4f::Load(fMat + 0);
20 sk4f a1 = sk4f::Load(fMat + 4);
21 sk4f a2 = sk4f::Load(fMat + 8);
22 sk4f a3 = sk4f::Load(fMat + 12);
23
24 sk4f b0 = sk4f::Load(other.fMat + 0);
25 sk4f b1 = sk4f::Load(other.fMat + 4);
26 sk4f b2 = sk4f::Load(other.fMat + 8);
27 sk4f b3 = sk4f::Load(other.fMat + 12);
28
29 auto eq = (a0 == b0) & (a1 == b1) & (a2 == b2) & (a3 == b3);
30 return (eq[0] & eq[1] & eq[2] & eq[3]) == ~0;
31}
32
33static void transpose_arrays(SkScalar dst[], const SkScalar src[]) {
34 dst[0] = src[0]; dst[1] = src[4]; dst[2] = src[8]; dst[3] = src[12];
35 dst[4] = src[1]; dst[5] = src[5]; dst[6] = src[9]; dst[7] = src[13];
36 dst[8] = src[2]; dst[9] = src[6]; dst[10] = src[10]; dst[11] = src[14];
37 dst[12] = src[3]; dst[13] = src[7]; dst[14] = src[11]; dst[15] = src[15];
38}
39
40void SkM44::getRowMajor(SkScalar v[]) const {
41 transpose_arrays(v, fMat);
42}
43
44SkM44& SkM44::setConcat(const SkM44& a, const SkM44& b) {
45 sk4f c0 = sk4f::Load(a.fMat + 0);
46 sk4f c1 = sk4f::Load(a.fMat + 4);
47 sk4f c2 = sk4f::Load(a.fMat + 8);
48 sk4f c3 = sk4f::Load(a.fMat + 12);
49
50 auto compute = [&](sk4f r) {
51 return c0*r[0] + (c1*r[1] + (c2*r[2] + c3*r[3]));
52 };
53
54 sk4f m0 = compute(sk4f::Load(b.fMat + 0));
55 sk4f m1 = compute(sk4f::Load(b.fMat + 4));
56 sk4f m2 = compute(sk4f::Load(b.fMat + 8));
57 sk4f m3 = compute(sk4f::Load(b.fMat + 12));
58
59 m0.store(fMat + 0);
60 m1.store(fMat + 4);
61 m2.store(fMat + 8);
62 m3.store(fMat + 12);
63 return *this;
64}
65
66SkM44& SkM44::preConcat(const SkMatrix& b) {
67 sk4f c0 = sk4f::Load(fMat + 0);
68 sk4f c1 = sk4f::Load(fMat + 4);
69 sk4f c3 = sk4f::Load(fMat + 12);
70
71 auto compute = [&](float r0, float r1, float r3) {
72 return (c0*r0 + (c1*r1 + c3*r3));
73 };
74
75 sk4f m0 = compute(b[0], b[3], b[6]);
76 sk4f m1 = compute(b[1], b[4], b[7]);
77 sk4f m3 = compute(b[2], b[5], b[8]);
78
79 m0.store(fMat + 0);
80 m1.store(fMat + 4);
81 m3.store(fMat + 12);
82 return *this;
83}
84
85SkM44& SkM44::preTranslate(SkScalar x, SkScalar y, SkScalar z) {
86 sk4f c0 = sk4f::Load(fMat + 0);
87 sk4f c1 = sk4f::Load(fMat + 4);
88 sk4f c2 = sk4f::Load(fMat + 8);
89 sk4f c3 = sk4f::Load(fMat + 12);
90
91 // only need to update the last column
92 (c0*x + (c1*y + (c2*z + c3))).store(fMat + 12);
93 return *this;
94}
95
96SkM44& SkM44::postTranslate(SkScalar x, SkScalar y, SkScalar z) {
97 sk4f t = { x, y, z, 0 };
98 (t * fMat[ 3] + sk4f::Load(fMat + 0)).store(fMat + 0);
99 (t * fMat[ 7] + sk4f::Load(fMat + 4)).store(fMat + 4);
100 (t * fMat[11] + sk4f::Load(fMat + 8)).store(fMat + 8);
101 (t * fMat[15] + sk4f::Load(fMat + 12)).store(fMat + 12);
102 return *this;
103}
104
105SkM44& SkM44::preScale(SkScalar x, SkScalar y) {
106 sk4f c0 = sk4f::Load(fMat + 0);
107 sk4f c1 = sk4f::Load(fMat + 4);
108
109 (c0 * x).store(fMat + 0);
110 (c1 * y).store(fMat + 4);
111 return *this;
112}
113
114SkV4 SkM44::map(float x, float y, float z, float w) const {
115 sk4f c0 = sk4f::Load(fMat + 0);
116 sk4f c1 = sk4f::Load(fMat + 4);
117 sk4f c2 = sk4f::Load(fMat + 8);
118 sk4f c3 = sk4f::Load(fMat + 12);
119
120 SkV4 v;
121 (c0*x + (c1*y + (c2*z + c3*w))).store(&v.x);
122 return v;
123}
124
125void SkM44::normalizePerspective() {
126 // If the bottom row of the matrix is [0, 0, 0, not_one], we will treat the matrix as if it
127 // is in perspective, even though it stills behaves like its affine. If we divide everything
128 // by the not_one value, then it will behave the same, but will be treated as affine,
129 // and therefore faster (e.g. clients can forward-difference calculations).
130 if (fMat[15] != 1 && fMat[15] != 0 && fMat[3] == 0 && fMat[7] == 0 && fMat[11] == 0) {
131 double inv = 1.0 / fMat[15];
132 (sk4f::Load(fMat + 0) * inv).store(fMat + 0);
133 (sk4f::Load(fMat + 4) * inv).store(fMat + 4);
134 (sk4f::Load(fMat + 8) * inv).store(fMat + 8);
135 (sk4f::Load(fMat + 12) * inv).store(fMat + 12);
136 fMat[15] = 1.0f;
137 }
138}
139
140///////////////////////////////////////////////////////////////////////////////
141
142/** We always perform the calculation in doubles, to avoid prematurely losing
143 precision along the way. This relies on the compiler automatically
144 promoting our SkScalar values to double (if needed).
145 */
146bool SkM44::invert(SkM44* inverse) const {
147 double a00 = fMat[0];
148 double a01 = fMat[1];
149 double a02 = fMat[2];
150 double a03 = fMat[3];
151 double a10 = fMat[4];
152 double a11 = fMat[5];
153 double a12 = fMat[6];
154 double a13 = fMat[7];
155 double a20 = fMat[8];
156 double a21 = fMat[9];
157 double a22 = fMat[10];
158 double a23 = fMat[11];
159 double a30 = fMat[12];
160 double a31 = fMat[13];
161 double a32 = fMat[14];
162 double a33 = fMat[15];
163
164 double b00 = a00 * a11 - a01 * a10;
165 double b01 = a00 * a12 - a02 * a10;
166 double b02 = a00 * a13 - a03 * a10;
167 double b03 = a01 * a12 - a02 * a11;
168 double b04 = a01 * a13 - a03 * a11;
169 double b05 = a02 * a13 - a03 * a12;
170 double b06 = a20 * a31 - a21 * a30;
171 double b07 = a20 * a32 - a22 * a30;
172 double b08 = a20 * a33 - a23 * a30;
173 double b09 = a21 * a32 - a22 * a31;
174 double b10 = a21 * a33 - a23 * a31;
175 double b11 = a22 * a33 - a23 * a32;
176
177 // Calculate the determinant
178 double det = b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06;
179
180 double invdet = sk_ieee_double_divide(1.0, det);
181 // If det is zero, we want to return false. However, we also want to return false if 1/det
182 // overflows to infinity (i.e. det is denormalized). All of this is subsumed by our final check
183 // at the bottom (that all 16 scalar matrix entries are finite).
184
185 b00 *= invdet;
186 b01 *= invdet;
187 b02 *= invdet;
188 b03 *= invdet;
189 b04 *= invdet;
190 b05 *= invdet;
191 b06 *= invdet;
192 b07 *= invdet;
193 b08 *= invdet;
194 b09 *= invdet;
195 b10 *= invdet;
196 b11 *= invdet;
197
198 SkScalar tmp[16] = {
199 SkDoubleToScalar(a11 * b11 - a12 * b10 + a13 * b09),
200 SkDoubleToScalar(a02 * b10 - a01 * b11 - a03 * b09),
201 SkDoubleToScalar(a31 * b05 - a32 * b04 + a33 * b03),
202 SkDoubleToScalar(a22 * b04 - a21 * b05 - a23 * b03),
203 SkDoubleToScalar(a12 * b08 - a10 * b11 - a13 * b07),
204 SkDoubleToScalar(a00 * b11 - a02 * b08 + a03 * b07),
205 SkDoubleToScalar(a32 * b02 - a30 * b05 - a33 * b01),
206 SkDoubleToScalar(a20 * b05 - a22 * b02 + a23 * b01),
207 SkDoubleToScalar(a10 * b10 - a11 * b08 + a13 * b06),
208 SkDoubleToScalar(a01 * b08 - a00 * b10 - a03 * b06),
209 SkDoubleToScalar(a30 * b04 - a31 * b02 + a33 * b00),
210 SkDoubleToScalar(a21 * b02 - a20 * b04 - a23 * b00),
211 SkDoubleToScalar(a11 * b07 - a10 * b09 - a12 * b06),
212 SkDoubleToScalar(a00 * b09 - a01 * b07 + a02 * b06),
213 SkDoubleToScalar(a31 * b01 - a30 * b03 - a32 * b00),
214 SkDoubleToScalar(a20 * b03 - a21 * b01 + a22 * b00),
215 };
216 if (!SkScalarsAreFinite(tmp, 16)) {
217 return false;
218 }
219 memcpy(inverse->fMat, tmp, sizeof(tmp));
220 return true;
221}
222
223SkM44 SkM44::transpose() const {
224 SkM44 trans(SkM44::kUninitialized_Constructor);
225 transpose_arrays(trans.fMat, fMat);
226 return trans;
227}
228
229SkM44& SkM44::setRotateUnitSinCos(SkV3 axis, SkScalar sinAngle, SkScalar cosAngle) {
230 // Taken from "Essential Mathematics for Games and Interactive Applications"
231 // James M. Van Verth and Lars M. Bishop -- third edition
232 SkScalar x = axis.x;
233 SkScalar y = axis.y;
234 SkScalar z = axis.z;
235 SkScalar c = cosAngle;
236 SkScalar s = sinAngle;
237 SkScalar t = 1 - c;
238
239 *this = { t*x*x + c, t*x*y - s*z, t*x*z + s*y, 0,
240 t*x*y + s*z, t*y*y + c, t*y*z - s*x, 0,
241 t*x*z - s*y, t*y*z + s*x, t*z*z + c, 0,
242 0, 0, 0, 1 };
243 return *this;
244}
245
246SkM44& SkM44::setRotate(SkV3 axis, SkScalar radians) {
247 SkScalar len = axis.length();
248 if (len > 0 && SkScalarIsFinite(len)) {
249 this->setRotateUnit(axis * (SK_Scalar1 / len), radians);
250 } else {
251 this->setIdentity();
252 }
253 return *this;
254}
255
256///////////////////////////////////////////////////////////////////////////////
257
258void SkM44::dump() const {
259 static const char* format = "|%g %g %g %g|\n"
260 "|%g %g %g %g|\n"
261 "|%g %g %g %g|\n"
262 "|%g %g %g %g|\n";
263 SkDebugf(format,
264 fMat[0], fMat[4], fMat[8], fMat[12],
265 fMat[1], fMat[5], fMat[9], fMat[13],
266 fMat[2], fMat[6], fMat[10], fMat[14],
267 fMat[3], fMat[7], fMat[11], fMat[15]);
268}
269
270static SkV3 normalize(SkV3 v) { return v * (1.0f / v.length()); }
271
272static SkV4 v4(SkV3 v, SkScalar w) { return {v.x, v.y, v.z, w}; }
273
274SkM44 Sk3LookAt(const SkV3& eye, const SkV3& center, const SkV3& up) {
275 SkV3 f = normalize(center - eye);
276 SkV3 u = normalize(up);
277 SkV3 s = normalize(f.cross(u));
278
279 SkM44 m(SkM44::kUninitialized_Constructor);
280 if (!SkM44::Cols(v4(s, 0), v4(s.cross(f), 0), v4(-f, 0), v4(eye, 1)).invert(&m)) {
281 m.setIdentity();
282 }
283 return m;
284}
285
286SkM44 Sk3Perspective(float near, float far, float angle) {
287 SkASSERT(far > near);
288
289 float denomInv = sk_ieee_float_divide(1, far - near);
290 float halfAngle = angle * 0.5f;
291 float cot = sk_float_cos(halfAngle) / sk_float_sin(halfAngle);
292
293 SkM44 m;
294 m.setRC(0, 0, cot);
295 m.setRC(1, 1, cot);
296 m.setRC(2, 2, (far + near) * denomInv);
297 m.setRC(2, 3, 2 * far * near * denomInv);
298 m.setRC(3, 2, -1);
299 return m;
300}
301