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
12 #include "src/core/SkMatrixPriv.h"
13 #include "src/core/SkPathPriv.h"
14
15 using sk4f = skvx::Vec<4, float>;
16 using sk2f = skvx::Vec<2, float>;
17
operator ==(const SkM44 & other) const18 bool SkM44::operator==(const SkM44& other) const {
19 if (this == &other) {
20 return true;
21 }
22
23 sk4f a0 = sk4f::Load(fMat + 0);
24 sk4f a1 = sk4f::Load(fMat + 4);
25 sk4f a2 = sk4f::Load(fMat + 8);
26 sk4f a3 = sk4f::Load(fMat + 12);
27
28 sk4f b0 = sk4f::Load(other.fMat + 0);
29 sk4f b1 = sk4f::Load(other.fMat + 4);
30 sk4f b2 = sk4f::Load(other.fMat + 8);
31 sk4f b3 = sk4f::Load(other.fMat + 12);
32
33 auto eq = (a0 == b0) & (a1 == b1) & (a2 == b2) & (a3 == b3);
34 return (eq[0] & eq[1] & eq[2] & eq[3]) == ~0;
35 }
36
transpose_arrays(SkScalar dst[],const SkScalar src[])37 static void transpose_arrays(SkScalar dst[], const SkScalar src[]) {
38 dst[0] = src[0]; dst[1] = src[4]; dst[2] = src[8]; dst[3] = src[12];
39 dst[4] = src[1]; dst[5] = src[5]; dst[6] = src[9]; dst[7] = src[13];
40 dst[8] = src[2]; dst[9] = src[6]; dst[10] = src[10]; dst[11] = src[14];
41 dst[12] = src[3]; dst[13] = src[7]; dst[14] = src[11]; dst[15] = src[15];
42 }
43
getRowMajor(SkScalar v[]) const44 void SkM44::getRowMajor(SkScalar v[]) const {
45 transpose_arrays(v, fMat);
46 }
47
setConcat(const SkM44 & a,const SkM44 & b)48 SkM44& SkM44::setConcat(const SkM44& a, const SkM44& b) {
49 sk4f c0 = sk4f::Load(a.fMat + 0);
50 sk4f c1 = sk4f::Load(a.fMat + 4);
51 sk4f c2 = sk4f::Load(a.fMat + 8);
52 sk4f c3 = sk4f::Load(a.fMat + 12);
53
54 auto compute = [&](sk4f r) {
55 return c0*r[0] + (c1*r[1] + (c2*r[2] + c3*r[3]));
56 };
57
58 sk4f m0 = compute(sk4f::Load(b.fMat + 0));
59 sk4f m1 = compute(sk4f::Load(b.fMat + 4));
60 sk4f m2 = compute(sk4f::Load(b.fMat + 8));
61 sk4f m3 = compute(sk4f::Load(b.fMat + 12));
62
63 m0.store(fMat + 0);
64 m1.store(fMat + 4);
65 m2.store(fMat + 8);
66 m3.store(fMat + 12);
67 return *this;
68 }
69
preConcat(const SkMatrix & b)70 SkM44& SkM44::preConcat(const SkMatrix& b) {
71 sk4f c0 = sk4f::Load(fMat + 0);
72 sk4f c1 = sk4f::Load(fMat + 4);
73 sk4f c3 = sk4f::Load(fMat + 12);
74
75 auto compute = [&](float r0, float r1, float r3) {
76 return (c0*r0 + (c1*r1 + c3*r3));
77 };
78
79 sk4f m0 = compute(b[0], b[3], b[6]);
80 sk4f m1 = compute(b[1], b[4], b[7]);
81 sk4f m3 = compute(b[2], b[5], b[8]);
82
83 m0.store(fMat + 0);
84 m1.store(fMat + 4);
85 m3.store(fMat + 12);
86 return *this;
87 }
88
preTranslate(SkScalar x,SkScalar y,SkScalar z)89 SkM44& SkM44::preTranslate(SkScalar x, SkScalar y, SkScalar z) {
90 sk4f c0 = sk4f::Load(fMat + 0);
91 sk4f c1 = sk4f::Load(fMat + 4);
92 sk4f c2 = sk4f::Load(fMat + 8);
93 sk4f c3 = sk4f::Load(fMat + 12);
94
95 // only need to update the last column
96 (c0*x + (c1*y + (c2*z + c3))).store(fMat + 12);
97 return *this;
98 }
99
postTranslate(SkScalar x,SkScalar y,SkScalar z)100 SkM44& SkM44::postTranslate(SkScalar x, SkScalar y, SkScalar z) {
101 sk4f t = { x, y, z, 0 };
102 (t * fMat[ 3] + sk4f::Load(fMat + 0)).store(fMat + 0);
103 (t * fMat[ 7] + sk4f::Load(fMat + 4)).store(fMat + 4);
104 (t * fMat[11] + sk4f::Load(fMat + 8)).store(fMat + 8);
105 (t * fMat[15] + sk4f::Load(fMat + 12)).store(fMat + 12);
106 return *this;
107 }
108
preScale(SkScalar x,SkScalar y)109 SkM44& SkM44::preScale(SkScalar x, SkScalar y) {
110 sk4f c0 = sk4f::Load(fMat + 0);
111 sk4f c1 = sk4f::Load(fMat + 4);
112
113 (c0 * x).store(fMat + 0);
114 (c1 * y).store(fMat + 4);
115 return *this;
116 }
117
preScale(SkScalar x,SkScalar y,SkScalar z)118 SkM44& SkM44::preScale(SkScalar x, SkScalar y, SkScalar z) {
119 sk4f c0 = sk4f::Load(fMat + 0);
120 sk4f c1 = sk4f::Load(fMat + 4);
121 sk4f c2 = sk4f::Load(fMat + 8);
122
123 (c0 * x).store(fMat + 0);
124 (c1 * y).store(fMat + 4);
125 (c2 * z).store(fMat + 8);
126 return *this;
127 }
128
map(float x,float y,float z,float w) const129 SkV4 SkM44::map(float x, float y, float z, float w) const {
130 sk4f c0 = sk4f::Load(fMat + 0);
131 sk4f c1 = sk4f::Load(fMat + 4);
132 sk4f c2 = sk4f::Load(fMat + 8);
133 sk4f c3 = sk4f::Load(fMat + 12);
134
135 SkV4 v;
136 (c0*x + (c1*y + (c2*z + c3*w))).store(&v.x);
137 return v;
138 }
139
map_rect_affine(const SkRect & src,const float mat[16])140 static SkRect map_rect_affine(const SkRect& src, const float mat[16]) {
141 // When multiplied against vectors of the form <x,y,x,y>, 'flip' allows a single min(sk4f, sk4f)
142 // to compute both the min and "negated" max between the xy coordinates. Once finished, another
143 // multiplication produces the original max.
144 const sk4f flip{1.f, 1.f, -1.f, -1.f};
145
146 // Since z = 0 and it's assumed ther's no perspective, only load the upper 2x2 and (tx,ty) in c3
147 sk4f c0 = skvx::shuffle<0,1,0,1>(sk2f::Load(mat + 0)) * flip;
148 sk4f c1 = skvx::shuffle<0,1,0,1>(sk2f::Load(mat + 4)) * flip;
149 sk4f c3 = skvx::shuffle<0,1,0,1>(sk2f::Load(mat + 12));
150
151 // Compute the min and max of the four transformed corners pre-translation; then translate once
152 // at the end.
153 sk4f minMax = c3 + flip * min(min(c0 * src.fLeft + c1 * src.fTop,
154 c0 * src.fRight + c1 * src.fTop),
155 min(c0 * src.fLeft + c1 * src.fBottom,
156 c0 * src.fRight + c1 * src.fBottom));
157
158 // minMax holds (min x, min y, max x, max y) so can be copied into an SkRect expecting l,t,r,b
159 SkRect r;
160 minMax.store(&r);
161 return r;
162 }
163
map_rect_perspective(const SkRect & src,const float mat[16])164 static SkRect map_rect_perspective(const SkRect& src, const float mat[16]) {
165 // Like map_rect_affine, z = 0 so we can skip the 3rd column, but we do need to compute w's
166 // for each corner of the src rect.
167 sk4f c0 = sk4f::Load(mat + 0);
168 sk4f c1 = sk4f::Load(mat + 4);
169 sk4f c3 = sk4f::Load(mat + 12);
170
171 // Unlike map_rect_affine, we do not defer the 4th column since we may need to homogeneous
172 // coordinates to clip against the w=0 plane
173 sk4f tl = c0 * src.fLeft + c1 * src.fTop + c3;
174 sk4f tr = c0 * src.fRight + c1 * src.fTop + c3;
175 sk4f bl = c0 * src.fLeft + c1 * src.fBottom + c3;
176 sk4f br = c0 * src.fRight + c1 * src.fBottom + c3;
177
178 // After clipping to w>0 and projecting to 2d, 'project' employs the same negation trick to
179 // compute min and max at the same time.
180 const sk4f flip{1.f, 1.f, -1.f, -1.f};
181 auto project = [&flip](const sk4f& p0, const sk4f& p1, const sk4f& p2) {
182 float w0 = p0[3];
183 if (w0 >= SkPathPriv::kW0PlaneDistance) {
184 // Unclipped, just divide by w
185 return flip * skvx::shuffle<0,1,0,1>(p0) / w0;
186 } else {
187 auto clip = [&](const sk4f& p) {
188 float w = p[3];
189 if (w >= SkPathPriv::kW0PlaneDistance) {
190 float t = (SkPathPriv::kW0PlaneDistance - w0) / (w - w0);
191 sk2f c = (t * skvx::shuffle<0,1>(p) + (1.f - t) * skvx::shuffle<0,1>(p0)) /
192 SkPathPriv::kW0PlaneDistance;
193
194 return flip * skvx::shuffle<0,1,0,1>(c);
195 } else {
196 return sk4f(SK_ScalarInfinity);
197 }
198 };
199 // Clip both edges leaving p0, and return the min/max of the two clipped points
200 // (since clip returns infinity when both p0 and 2nd vertex have w<0, it'll
201 // automatically be ignored).
202 return min(clip(p1), clip(p2));
203 }
204 };
205
206 // Project all 4 corners, and pass in their adjacent vertices for clipping if it has w < 0,
207 // then accumulate the min and max xy's.
208 sk4f minMax = flip * min(min(project(tl, tr, bl), project(tr, br, tl)),
209 min(project(br, bl, tr), project(bl, tl, br)));
210
211 SkRect r;
212 minMax.store(&r);
213 return r;
214 }
215
MapRect(const SkM44 & m,const SkRect & src)216 SkRect SkMatrixPriv::MapRect(const SkM44& m, const SkRect& src) {
217 const bool hasPerspective =
218 m.fMat[3] != 0 || m.fMat[7] != 0 || m.fMat[11] != 0 || m.fMat[15] != 1;
219 if (hasPerspective) {
220 return map_rect_perspective(src, m.fMat);
221 } else {
222 return map_rect_affine(src, m.fMat);
223 }
224 }
225
normalizePerspective()226 void SkM44::normalizePerspective() {
227 // If the bottom row of the matrix is [0, 0, 0, not_one], we will treat the matrix as if it
228 // is in perspective, even though it stills behaves like its affine. If we divide everything
229 // by the not_one value, then it will behave the same, but will be treated as affine,
230 // and therefore faster (e.g. clients can forward-difference calculations).
231 if (fMat[15] != 1 && fMat[15] != 0 && fMat[3] == 0 && fMat[7] == 0 && fMat[11] == 0) {
232 double inv = 1.0 / fMat[15];
233 (sk4f::Load(fMat + 0) * inv).store(fMat + 0);
234 (sk4f::Load(fMat + 4) * inv).store(fMat + 4);
235 (sk4f::Load(fMat + 8) * inv).store(fMat + 8);
236 (sk4f::Load(fMat + 12) * inv).store(fMat + 12);
237 fMat[15] = 1.0f;
238 }
239 }
240
241 ///////////////////////////////////////////////////////////////////////////////
242
243 /** We always perform the calculation in doubles, to avoid prematurely losing
244 precision along the way. This relies on the compiler automatically
245 promoting our SkScalar values to double (if needed).
246 */
invert(SkM44 * inverse) const247 bool SkM44::invert(SkM44* inverse) const {
248 double a00 = fMat[0];
249 double a01 = fMat[1];
250 double a02 = fMat[2];
251 double a03 = fMat[3];
252 double a10 = fMat[4];
253 double a11 = fMat[5];
254 double a12 = fMat[6];
255 double a13 = fMat[7];
256 double a20 = fMat[8];
257 double a21 = fMat[9];
258 double a22 = fMat[10];
259 double a23 = fMat[11];
260 double a30 = fMat[12];
261 double a31 = fMat[13];
262 double a32 = fMat[14];
263 double a33 = fMat[15];
264
265 double b00 = a00 * a11 - a01 * a10;
266 double b01 = a00 * a12 - a02 * a10;
267 double b02 = a00 * a13 - a03 * a10;
268 double b03 = a01 * a12 - a02 * a11;
269 double b04 = a01 * a13 - a03 * a11;
270 double b05 = a02 * a13 - a03 * a12;
271 double b06 = a20 * a31 - a21 * a30;
272 double b07 = a20 * a32 - a22 * a30;
273 double b08 = a20 * a33 - a23 * a30;
274 double b09 = a21 * a32 - a22 * a31;
275 double b10 = a21 * a33 - a23 * a31;
276 double b11 = a22 * a33 - a23 * a32;
277
278 // Calculate the determinant
279 double det = b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06;
280
281 double invdet = sk_ieee_double_divide(1.0, det);
282 // If det is zero, we want to return false. However, we also want to return false if 1/det
283 // overflows to infinity (i.e. det is denormalized). All of this is subsumed by our final check
284 // at the bottom (that all 16 scalar matrix entries are finite).
285
286 b00 *= invdet;
287 b01 *= invdet;
288 b02 *= invdet;
289 b03 *= invdet;
290 b04 *= invdet;
291 b05 *= invdet;
292 b06 *= invdet;
293 b07 *= invdet;
294 b08 *= invdet;
295 b09 *= invdet;
296 b10 *= invdet;
297 b11 *= invdet;
298
299 SkScalar tmp[16] = {
300 SkDoubleToScalar(a11 * b11 - a12 * b10 + a13 * b09),
301 SkDoubleToScalar(a02 * b10 - a01 * b11 - a03 * b09),
302 SkDoubleToScalar(a31 * b05 - a32 * b04 + a33 * b03),
303 SkDoubleToScalar(a22 * b04 - a21 * b05 - a23 * b03),
304 SkDoubleToScalar(a12 * b08 - a10 * b11 - a13 * b07),
305 SkDoubleToScalar(a00 * b11 - a02 * b08 + a03 * b07),
306 SkDoubleToScalar(a32 * b02 - a30 * b05 - a33 * b01),
307 SkDoubleToScalar(a20 * b05 - a22 * b02 + a23 * b01),
308 SkDoubleToScalar(a10 * b10 - a11 * b08 + a13 * b06),
309 SkDoubleToScalar(a01 * b08 - a00 * b10 - a03 * b06),
310 SkDoubleToScalar(a30 * b04 - a31 * b02 + a33 * b00),
311 SkDoubleToScalar(a21 * b02 - a20 * b04 - a23 * b00),
312 SkDoubleToScalar(a11 * b07 - a10 * b09 - a12 * b06),
313 SkDoubleToScalar(a00 * b09 - a01 * b07 + a02 * b06),
314 SkDoubleToScalar(a31 * b01 - a30 * b03 - a32 * b00),
315 SkDoubleToScalar(a20 * b03 - a21 * b01 + a22 * b00),
316 };
317 if (!SkScalarsAreFinite(tmp, 16)) {
318 return false;
319 }
320 memcpy(inverse->fMat, tmp, sizeof(tmp));
321 return true;
322 }
323
transpose() const324 SkM44 SkM44::transpose() const {
325 SkM44 trans(SkM44::kUninitialized_Constructor);
326 transpose_arrays(trans.fMat, fMat);
327 return trans;
328 }
329
setRotateUnitSinCos(SkV3 axis,SkScalar sinAngle,SkScalar cosAngle)330 SkM44& SkM44::setRotateUnitSinCos(SkV3 axis, SkScalar sinAngle, SkScalar cosAngle) {
331 // Taken from "Essential Mathematics for Games and Interactive Applications"
332 // James M. Van Verth and Lars M. Bishop -- third edition
333 SkScalar x = axis.x;
334 SkScalar y = axis.y;
335 SkScalar z = axis.z;
336 SkScalar c = cosAngle;
337 SkScalar s = sinAngle;
338 SkScalar t = 1 - c;
339
340 *this = { t*x*x + c, t*x*y - s*z, t*x*z + s*y, 0,
341 t*x*y + s*z, t*y*y + c, t*y*z - s*x, 0,
342 t*x*z - s*y, t*y*z + s*x, t*z*z + c, 0,
343 0, 0, 0, 1 };
344 return *this;
345 }
346
setRotate(SkV3 axis,SkScalar radians)347 SkM44& SkM44::setRotate(SkV3 axis, SkScalar radians) {
348 SkScalar len = axis.length();
349 if (len > 0 && SkScalarIsFinite(len)) {
350 this->setRotateUnit(axis * (SK_Scalar1 / len), radians);
351 } else {
352 this->setIdentity();
353 }
354 return *this;
355 }
356
357 ///////////////////////////////////////////////////////////////////////////////
358
dump() const359 void SkM44::dump() const {
360 static const char* format = "|%g %g %g %g|\n"
361 "|%g %g %g %g|\n"
362 "|%g %g %g %g|\n"
363 "|%g %g %g %g|\n";
364 SkDebugf(format,
365 fMat[0], fMat[4], fMat[8], fMat[12],
366 fMat[1], fMat[5], fMat[9], fMat[13],
367 fMat[2], fMat[6], fMat[10], fMat[14],
368 fMat[3], fMat[7], fMat[11], fMat[15]);
369 }
370
371 ///////////////////////////////////////////////////////////////////////////////
372
RectToRect(const SkRect & src,const SkRect & dst)373 SkM44 SkM44::RectToRect(const SkRect& src, const SkRect& dst) {
374 if (src.isEmpty()) {
375 return SkM44();
376 } else if (dst.isEmpty()) {
377 return SkM44::Scale(0.f, 0.f, 0.f);
378 }
379
380 float sx = dst.width() / src.width();
381 float sy = dst.height() / src.height();
382
383 float tx = dst.fLeft - sx * src.fLeft;
384 float ty = dst.fTop - sy * src.fTop;
385
386 return SkM44{sx, 0.f, 0.f, tx,
387 0.f, sy, 0.f, ty,
388 0.f, 0.f, 1.f, 0.f,
389 0.f, 0.f, 0.f, 1.f};
390 }
391
normalize(SkV3 v)392 static SkV3 normalize(SkV3 v) { return v * (1.0f / v.length()); }
393
v4(SkV3 v,SkScalar w)394 static SkV4 v4(SkV3 v, SkScalar w) { return {v.x, v.y, v.z, w}; }
395
LookAt(const SkV3 & eye,const SkV3 & center,const SkV3 & up)396 SkM44 SkM44::LookAt(const SkV3& eye, const SkV3& center, const SkV3& up) {
397 SkV3 f = normalize(center - eye);
398 SkV3 u = normalize(up);
399 SkV3 s = normalize(f.cross(u));
400
401 SkM44 m(SkM44::kUninitialized_Constructor);
402 if (!SkM44::Cols(v4(s, 0), v4(s.cross(f), 0), v4(-f, 0), v4(eye, 1)).invert(&m)) {
403 m.setIdentity();
404 }
405 return m;
406 }
407
Perspective(float near,float far,float angle)408 SkM44 SkM44::Perspective(float near, float far, float angle) {
409 SkASSERT(far > near);
410
411 float denomInv = sk_ieee_float_divide(1, far - near);
412 float halfAngle = angle * 0.5f;
413 float cot = sk_float_cos(halfAngle) / sk_float_sin(halfAngle);
414
415 SkM44 m;
416 m.setRC(0, 0, cot);
417 m.setRC(1, 1, cot);
418 m.setRC(2, 2, (far + near) * denomInv);
419 m.setRC(2, 3, 2 * far * near * denomInv);
420 m.setRC(3, 2, -1);
421 return m;
422 }
423