1/*
2 * Copyright 2011 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 "SkMatrix44.h"
9
10static inline bool eq4(const SkMScalar* SK_RESTRICT a,
11                      const SkMScalar* SK_RESTRICT b) {
12    return (a[0] == b[0]) & (a[1] == b[1]) & (a[2] == b[2]) & (a[3] == b[3]);
13}
14
15bool SkMatrix44::operator==(const SkMatrix44& other) const {
16    if (this == &other) {
17        return true;
18    }
19
20    if (this->isTriviallyIdentity() && other.isTriviallyIdentity()) {
21        return true;
22    }
23
24    const SkMScalar* SK_RESTRICT a = &fMat[0][0];
25    const SkMScalar* SK_RESTRICT b = &other.fMat[0][0];
26
27#if 0
28    for (int i = 0; i < 16; ++i) {
29        if (a[i] != b[i]) {
30            return false;
31        }
32    }
33    return true;
34#else
35    // to reduce branch instructions, we compare 4 at a time.
36    // see bench/Matrix44Bench.cpp for test.
37    if (!eq4(&a[0], &b[0])) {
38        return false;
39    }
40    if (!eq4(&a[4], &b[4])) {
41        return false;
42    }
43    if (!eq4(&a[8], &b[8])) {
44        return false;
45    }
46    return eq4(&a[12], &b[12]);
47#endif
48}
49
50///////////////////////////////////////////////////////////////////////////////
51
52int SkMatrix44::computeTypeMask() const {
53    unsigned mask = 0;
54
55    if (0 != perspX() || 0 != perspY() || 0 != perspZ() || 1 != fMat[3][3]) {
56        return kTranslate_Mask | kScale_Mask | kAffine_Mask | kPerspective_Mask;
57    }
58
59    if (0 != transX() || 0 != transY() || 0 != transZ()) {
60        mask |= kTranslate_Mask;
61    }
62
63    if (1 != scaleX() || 1 != scaleY() || 1 != scaleZ()) {
64        mask |= kScale_Mask;
65    }
66
67    if (0 != fMat[1][0] || 0 != fMat[0][1] || 0 != fMat[0][2] ||
68        0 != fMat[2][0] || 0 != fMat[1][2] || 0 != fMat[2][1]) {
69            mask |= kAffine_Mask;
70    }
71
72    return mask;
73}
74
75///////////////////////////////////////////////////////////////////////////////
76
77void SkMatrix44::asColMajorf(float dst[]) const {
78    const SkMScalar* src = &fMat[0][0];
79#ifdef SK_MSCALAR_IS_DOUBLE
80    for (int i = 0; i < 16; ++i) {
81        dst[i] = SkMScalarToFloat(src[i]);
82    }
83#elif defined SK_MSCALAR_IS_FLOAT
84    memcpy(dst, src, 16 * sizeof(float));
85#endif
86}
87
88void SkMatrix44::asColMajord(double dst[]) const {
89    const SkMScalar* src = &fMat[0][0];
90#ifdef SK_MSCALAR_IS_DOUBLE
91    memcpy(dst, src, 16 * sizeof(double));
92#elif defined SK_MSCALAR_IS_FLOAT
93    for (int i = 0; i < 16; ++i) {
94        dst[i] = SkMScalarToDouble(src[i]);
95    }
96#endif
97}
98
99void SkMatrix44::asRowMajorf(float dst[]) const {
100    const SkMScalar* src = &fMat[0][0];
101    for (int i = 0; i < 4; ++i) {
102        dst[0] = SkMScalarToFloat(src[0]);
103        dst[4] = SkMScalarToFloat(src[1]);
104        dst[8] = SkMScalarToFloat(src[2]);
105        dst[12] = SkMScalarToFloat(src[3]);
106        src += 4;
107        dst += 1;
108    }
109}
110
111void SkMatrix44::asRowMajord(double dst[]) const {
112    const SkMScalar* src = &fMat[0][0];
113    for (int i = 0; i < 4; ++i) {
114        dst[0] = SkMScalarToDouble(src[0]);
115        dst[4] = SkMScalarToDouble(src[1]);
116        dst[8] = SkMScalarToDouble(src[2]);
117        dst[12] = SkMScalarToDouble(src[3]);
118        src += 4;
119        dst += 1;
120    }
121}
122
123void SkMatrix44::setColMajorf(const float src[]) {
124    SkMScalar* dst = &fMat[0][0];
125#ifdef SK_MSCALAR_IS_DOUBLE
126    for (int i = 0; i < 16; ++i) {
127        dst[i] = SkMScalarToFloat(src[i]);
128    }
129#elif defined SK_MSCALAR_IS_FLOAT
130    memcpy(dst, src, 16 * sizeof(float));
131#endif
132
133    this->dirtyTypeMask();
134}
135
136void SkMatrix44::setColMajord(const double src[]) {
137    SkMScalar* dst = &fMat[0][0];
138#ifdef SK_MSCALAR_IS_DOUBLE
139    memcpy(dst, src, 16 * sizeof(double));
140#elif defined SK_MSCALAR_IS_FLOAT
141    for (int i = 0; i < 16; ++i) {
142        dst[i] = SkDoubleToMScalar(src[i]);
143    }
144#endif
145
146    this->dirtyTypeMask();
147}
148
149void SkMatrix44::setRowMajorf(const float src[]) {
150    SkMScalar* dst = &fMat[0][0];
151    for (int i = 0; i < 4; ++i) {
152        dst[0] = SkMScalarToFloat(src[0]);
153        dst[4] = SkMScalarToFloat(src[1]);
154        dst[8] = SkMScalarToFloat(src[2]);
155        dst[12] = SkMScalarToFloat(src[3]);
156        src += 4;
157        dst += 1;
158    }
159    this->dirtyTypeMask();
160}
161
162void SkMatrix44::setRowMajord(const double src[]) {
163    SkMScalar* dst = &fMat[0][0];
164    for (int i = 0; i < 4; ++i) {
165        dst[0] = SkDoubleToMScalar(src[0]);
166        dst[4] = SkDoubleToMScalar(src[1]);
167        dst[8] = SkDoubleToMScalar(src[2]);
168        dst[12] = SkDoubleToMScalar(src[3]);
169        src += 4;
170        dst += 1;
171    }
172    this->dirtyTypeMask();
173}
174
175///////////////////////////////////////////////////////////////////////////////
176
177const SkMatrix44& SkMatrix44::I() {
178    static const SkMatrix44 gIdentity44(kIdentity_Constructor);
179    return gIdentity44;
180}
181
182void SkMatrix44::setIdentity() {
183    fMat[0][0] = 1;
184    fMat[0][1] = 0;
185    fMat[0][2] = 0;
186    fMat[0][3] = 0;
187    fMat[1][0] = 0;
188    fMat[1][1] = 1;
189    fMat[1][2] = 0;
190    fMat[1][3] = 0;
191    fMat[2][0] = 0;
192    fMat[2][1] = 0;
193    fMat[2][2] = 1;
194    fMat[2][3] = 0;
195    fMat[3][0] = 0;
196    fMat[3][1] = 0;
197    fMat[3][2] = 0;
198    fMat[3][3] = 1;
199    this->setTypeMask(kIdentity_Mask);
200}
201
202void SkMatrix44::set3x3(SkMScalar m00, SkMScalar m01, SkMScalar m02,
203                        SkMScalar m10, SkMScalar m11, SkMScalar m12,
204                        SkMScalar m20, SkMScalar m21, SkMScalar m22) {
205    fMat[0][0] = m00; fMat[0][1] = m01; fMat[0][2] = m02; fMat[0][3] = 0;
206    fMat[1][0] = m10; fMat[1][1] = m11; fMat[1][2] = m12; fMat[1][3] = 0;
207    fMat[2][0] = m20; fMat[2][1] = m21; fMat[2][2] = m22; fMat[2][3] = 0;
208    fMat[3][0] = 0;   fMat[3][1] = 0;   fMat[3][2] = 0;   fMat[3][3] = 1;
209    this->dirtyTypeMask();
210}
211
212///////////////////////////////////////////////////////////////////////////////
213
214void SkMatrix44::setTranslate(SkMScalar dx, SkMScalar dy, SkMScalar dz) {
215    this->setIdentity();
216
217    if (!dx && !dy && !dz) {
218        return;
219    }
220
221    fMat[3][0] = dx;
222    fMat[3][1] = dy;
223    fMat[3][2] = dz;
224    this->setTypeMask(kTranslate_Mask);
225}
226
227void SkMatrix44::preTranslate(SkMScalar dx, SkMScalar dy, SkMScalar dz) {
228    if (!dx && !dy && !dz) {
229        return;
230    }
231
232    for (int i = 0; i < 4; ++i) {
233        fMat[3][i] = fMat[0][i] * dx + fMat[1][i] * dy + fMat[2][i] * dz + fMat[3][i];
234    }
235    this->dirtyTypeMask();
236}
237
238void SkMatrix44::postTranslate(SkMScalar dx, SkMScalar dy, SkMScalar dz) {
239    if (!dx && !dy && !dz) {
240        return;
241    }
242
243    if (this->getType() & kPerspective_Mask) {
244        for (int i = 0; i < 4; ++i) {
245            fMat[i][0] += fMat[i][3] * dx;
246            fMat[i][1] += fMat[i][3] * dy;
247            fMat[i][2] += fMat[i][3] * dz;
248        }
249    } else {
250        fMat[3][0] += dx;
251        fMat[3][1] += dy;
252        fMat[3][2] += dz;
253        this->dirtyTypeMask();
254    }
255}
256
257///////////////////////////////////////////////////////////////////////////////
258
259void SkMatrix44::setScale(SkMScalar sx, SkMScalar sy, SkMScalar sz) {
260    this->setIdentity();
261
262    if (1 == sx && 1 == sy && 1 == sz) {
263        return;
264    }
265
266    fMat[0][0] = sx;
267    fMat[1][1] = sy;
268    fMat[2][2] = sz;
269    this->setTypeMask(kScale_Mask);
270}
271
272void SkMatrix44::preScale(SkMScalar sx, SkMScalar sy, SkMScalar sz) {
273    if (1 == sx && 1 == sy && 1 == sz) {
274        return;
275    }
276
277    // The implementation matrix * pureScale can be shortcut
278    // by knowing that pureScale components effectively scale
279    // the columns of the original matrix.
280    for (int i = 0; i < 4; i++) {
281        fMat[0][i] *= sx;
282        fMat[1][i] *= sy;
283        fMat[2][i] *= sz;
284    }
285    this->dirtyTypeMask();
286}
287
288void SkMatrix44::postScale(SkMScalar sx, SkMScalar sy, SkMScalar sz) {
289    if (1 == sx && 1 == sy && 1 == sz) {
290        return;
291    }
292
293    for (int i = 0; i < 4; i++) {
294        fMat[i][0] *= sx;
295        fMat[i][1] *= sy;
296        fMat[i][2] *= sz;
297    }
298    this->dirtyTypeMask();
299}
300
301///////////////////////////////////////////////////////////////////////////////
302
303void SkMatrix44::setRotateAbout(SkMScalar x, SkMScalar y, SkMScalar z,
304                                SkMScalar radians) {
305    double len2 = (double)x * x + (double)y * y + (double)z * z;
306    if (1 != len2) {
307        if (0 == len2) {
308            this->setIdentity();
309            return;
310        }
311        double scale = 1 / sqrt(len2);
312        x = SkDoubleToMScalar(x * scale);
313        y = SkDoubleToMScalar(y * scale);
314        z = SkDoubleToMScalar(z * scale);
315    }
316    this->setRotateAboutUnit(x, y, z, radians);
317}
318
319void SkMatrix44::setRotateAboutUnit(SkMScalar x, SkMScalar y, SkMScalar z,
320                                    SkMScalar radians) {
321    double c = cos(radians);
322    double s = sin(radians);
323    double C = 1 - c;
324    double xs = x * s;
325    double ys = y * s;
326    double zs = z * s;
327    double xC = x * C;
328    double yC = y * C;
329    double zC = z * C;
330    double xyC = x * yC;
331    double yzC = y * zC;
332    double zxC = z * xC;
333
334    // if you're looking at wikipedia, remember that we're column major.
335    this->set3x3(SkDoubleToMScalar(x * xC + c),     // scale x
336                 SkDoubleToMScalar(xyC + zs),       // skew x
337                 SkDoubleToMScalar(zxC - ys),       // trans x
338
339                 SkDoubleToMScalar(xyC - zs),       // skew y
340                 SkDoubleToMScalar(y * yC + c),     // scale y
341                 SkDoubleToMScalar(yzC + xs),       // trans y
342
343                 SkDoubleToMScalar(zxC + ys),       // persp x
344                 SkDoubleToMScalar(yzC - xs),       // persp y
345                 SkDoubleToMScalar(z * zC + c));    // persp 2
346}
347
348///////////////////////////////////////////////////////////////////////////////
349
350static bool bits_isonly(int value, int mask) {
351    return 0 == (value & ~mask);
352}
353
354void SkMatrix44::setConcat(const SkMatrix44& a, const SkMatrix44& b) {
355    const SkMatrix44::TypeMask a_mask = a.getType();
356    const SkMatrix44::TypeMask b_mask = b.getType();
357
358    if (kIdentity_Mask == a_mask) {
359        *this = b;
360        return;
361    }
362    if (kIdentity_Mask == b_mask) {
363        *this = a;
364        return;
365    }
366
367    bool useStorage = (this == &a || this == &b);
368    SkMScalar storage[16];
369    SkMScalar* result = useStorage ? storage : &fMat[0][0];
370
371    // Both matrices are at most scale+translate
372    if (bits_isonly(a_mask | b_mask, kScale_Mask | kTranslate_Mask)) {
373        result[0] = a.fMat[0][0] * b.fMat[0][0];
374        result[1] = result[2] = result[3] = result[4] = 0;
375        result[5] = a.fMat[1][1] * b.fMat[1][1];
376        result[6] = result[7] = result[8] = result[9] = 0;
377        result[10] = a.fMat[2][2] * b.fMat[2][2];
378        result[11] = 0;
379        result[12] = a.fMat[0][0] * b.fMat[3][0] + a.fMat[3][0];
380        result[13] = a.fMat[1][1] * b.fMat[3][1] + a.fMat[3][1];
381        result[14] = a.fMat[2][2] * b.fMat[3][2] + a.fMat[3][2];
382        result[15] = 1;
383    } else {
384        for (int j = 0; j < 4; j++) {
385            for (int i = 0; i < 4; i++) {
386                double value = 0;
387                for (int k = 0; k < 4; k++) {
388                    value += SkMScalarToDouble(a.fMat[k][i]) * b.fMat[j][k];
389                }
390                *result++ = SkDoubleToMScalar(value);
391            }
392        }
393    }
394
395    if (useStorage) {
396        memcpy(fMat, storage, sizeof(storage));
397    }
398    this->dirtyTypeMask();
399}
400
401///////////////////////////////////////////////////////////////////////////////
402
403/** We always perform the calculation in doubles, to avoid prematurely losing
404    precision along the way. This relies on the compiler automatically
405    promoting our SkMScalar values to double (if needed).
406 */
407double SkMatrix44::determinant() const {
408    if (this->isIdentity()) {
409        return 1;
410    }
411    if (this->isScaleTranslate()) {
412        return fMat[0][0] * fMat[1][1] * fMat[2][2] * fMat[3][3];
413    }
414
415    double a00 = fMat[0][0];
416    double a01 = fMat[0][1];
417    double a02 = fMat[0][2];
418    double a03 = fMat[0][3];
419    double a10 = fMat[1][0];
420    double a11 = fMat[1][1];
421    double a12 = fMat[1][2];
422    double a13 = fMat[1][3];
423    double a20 = fMat[2][0];
424    double a21 = fMat[2][1];
425    double a22 = fMat[2][2];
426    double a23 = fMat[2][3];
427    double a30 = fMat[3][0];
428    double a31 = fMat[3][1];
429    double a32 = fMat[3][2];
430    double a33 = fMat[3][3];
431
432    double b00 = a00 * a11 - a01 * a10;
433    double b01 = a00 * a12 - a02 * a10;
434    double b02 = a00 * a13 - a03 * a10;
435    double b03 = a01 * a12 - a02 * a11;
436    double b04 = a01 * a13 - a03 * a11;
437    double b05 = a02 * a13 - a03 * a12;
438    double b06 = a20 * a31 - a21 * a30;
439    double b07 = a20 * a32 - a22 * a30;
440    double b08 = a20 * a33 - a23 * a30;
441    double b09 = a21 * a32 - a22 * a31;
442    double b10 = a21 * a33 - a23 * a31;
443    double b11 = a22 * a33 - a23 * a32;
444
445    // Calculate the determinant
446    return b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06;
447}
448
449///////////////////////////////////////////////////////////////////////////////
450
451static bool is_matrix_finite(const SkMatrix44& matrix) {
452    SkMScalar accumulator = 0;
453    for (int row = 0; row < 4; ++row) {
454        for (int col = 0; col < 4; ++col) {
455            accumulator *= matrix.get(row, col);
456        }
457    }
458    return accumulator == 0;
459}
460
461bool SkMatrix44::invert(SkMatrix44* storage) const {
462    if (this->isIdentity()) {
463        if (storage) {
464            storage->setIdentity();
465        }
466        return true;
467    }
468
469    if (this->isTranslate()) {
470        if (storage) {
471            storage->setTranslate(-fMat[3][0], -fMat[3][1], -fMat[3][2]);
472        }
473        return true;
474    }
475
476    SkMatrix44 tmp(kUninitialized_Constructor);
477    // Use storage if it's available and distinct from this matrix.
478    SkMatrix44* inverse = (storage && storage != this) ? storage : &tmp;
479    if (this->isScaleTranslate()) {
480        if (0 == fMat[0][0] * fMat[1][1] * fMat[2][2]) {
481            return false;
482        }
483
484        double invXScale = 1 / fMat[0][0];
485        double invYScale = 1 / fMat[1][1];
486        double invZScale = 1 / fMat[2][2];
487
488        inverse->fMat[0][0] = SkDoubleToMScalar(invXScale);
489        inverse->fMat[0][1] = 0;
490        inverse->fMat[0][2] = 0;
491        inverse->fMat[0][3] = 0;
492
493        inverse->fMat[1][0] = 0;
494        inverse->fMat[1][1] = SkDoubleToMScalar(invYScale);
495        inverse->fMat[1][2] = 0;
496        inverse->fMat[1][3] = 0;
497
498        inverse->fMat[2][0] = 0;
499        inverse->fMat[2][1] = 0;
500        inverse->fMat[2][2] = SkDoubleToMScalar(invZScale);
501        inverse->fMat[2][3] = 0;
502
503        inverse->fMat[3][0] = SkDoubleToMScalar(-fMat[3][0] * invXScale);
504        inverse->fMat[3][1] = SkDoubleToMScalar(-fMat[3][1] * invYScale);
505        inverse->fMat[3][2] = SkDoubleToMScalar(-fMat[3][2] * invZScale);
506        inverse->fMat[3][3] = 1;
507
508        inverse->setTypeMask(this->getType());
509
510        if (!is_matrix_finite(*inverse)) {
511            return false;
512        }
513        if (storage && inverse != storage) {
514            *storage = *inverse;
515        }
516        return true;
517    }
518
519    double a00 = fMat[0][0];
520    double a01 = fMat[0][1];
521    double a02 = fMat[0][2];
522    double a03 = fMat[0][3];
523    double a10 = fMat[1][0];
524    double a11 = fMat[1][1];
525    double a12 = fMat[1][2];
526    double a13 = fMat[1][3];
527    double a20 = fMat[2][0];
528    double a21 = fMat[2][1];
529    double a22 = fMat[2][2];
530    double a23 = fMat[2][3];
531    double a30 = fMat[3][0];
532    double a31 = fMat[3][1];
533    double a32 = fMat[3][2];
534    double a33 = fMat[3][3];
535
536    if (!(this->getType() & kPerspective_Mask)) {
537        // If we know the matrix has no perspective, then the perspective
538        // component is (0, 0, 0, 1). We can use this information to save a lot
539        // of arithmetic that would otherwise be spent to compute the inverse
540        // of a general matrix.
541
542        SkASSERT(a03 == 0);
543        SkASSERT(a13 == 0);
544        SkASSERT(a23 == 0);
545        SkASSERT(a33 == 1);
546
547        double b00 = a00 * a11 - a01 * a10;
548        double b01 = a00 * a12 - a02 * a10;
549        double b03 = a01 * a12 - a02 * a11;
550        double b06 = a20 * a31 - a21 * a30;
551        double b07 = a20 * a32 - a22 * a30;
552        double b08 = a20;
553        double b09 = a21 * a32 - a22 * a31;
554        double b10 = a21;
555        double b11 = a22;
556
557        // Calculate the determinant
558        double det = b00 * b11 - b01 * b10 + b03 * b08;
559
560        double invdet = 1.0 / det;
561        // If det is zero, we want to return false. However, we also want to return false
562        // if 1/det overflows to infinity (i.e. det is denormalized). Both of these are
563        // handled by checking that 1/det is finite.
564        if (!sk_float_isfinite(invdet)) {
565            return false;
566        }
567
568        b00 *= invdet;
569        b01 *= invdet;
570        b03 *= invdet;
571        b06 *= invdet;
572        b07 *= invdet;
573        b08 *= invdet;
574        b09 *= invdet;
575        b10 *= invdet;
576        b11 *= invdet;
577
578        inverse->fMat[0][0] = SkDoubleToMScalar(a11 * b11 - a12 * b10);
579        inverse->fMat[0][1] = SkDoubleToMScalar(a02 * b10 - a01 * b11);
580        inverse->fMat[0][2] = SkDoubleToMScalar(b03);
581        inverse->fMat[0][3] = 0;
582        inverse->fMat[1][0] = SkDoubleToMScalar(a12 * b08 - a10 * b11);
583        inverse->fMat[1][1] = SkDoubleToMScalar(a00 * b11 - a02 * b08);
584        inverse->fMat[1][2] = SkDoubleToMScalar(-b01);
585        inverse->fMat[1][3] = 0;
586        inverse->fMat[2][0] = SkDoubleToMScalar(a10 * b10 - a11 * b08);
587        inverse->fMat[2][1] = SkDoubleToMScalar(a01 * b08 - a00 * b10);
588        inverse->fMat[2][2] = SkDoubleToMScalar(b00);
589        inverse->fMat[2][3] = 0;
590        inverse->fMat[3][0] = SkDoubleToMScalar(a11 * b07 - a10 * b09 - a12 * b06);
591        inverse->fMat[3][1] = SkDoubleToMScalar(a00 * b09 - a01 * b07 + a02 * b06);
592        inverse->fMat[3][2] = SkDoubleToMScalar(a31 * b01 - a30 * b03 - a32 * b00);
593        inverse->fMat[3][3] = 1;
594
595        inverse->setTypeMask(this->getType());
596        if (!is_matrix_finite(*inverse)) {
597            return false;
598        }
599        if (storage && inverse != storage) {
600            *storage = *inverse;
601        }
602        return true;
603    }
604
605    double b00 = a00 * a11 - a01 * a10;
606    double b01 = a00 * a12 - a02 * a10;
607    double b02 = a00 * a13 - a03 * a10;
608    double b03 = a01 * a12 - a02 * a11;
609    double b04 = a01 * a13 - a03 * a11;
610    double b05 = a02 * a13 - a03 * a12;
611    double b06 = a20 * a31 - a21 * a30;
612    double b07 = a20 * a32 - a22 * a30;
613    double b08 = a20 * a33 - a23 * a30;
614    double b09 = a21 * a32 - a22 * a31;
615    double b10 = a21 * a33 - a23 * a31;
616    double b11 = a22 * a33 - a23 * a32;
617
618    // Calculate the determinant
619    double det = b00 * b11 - b01 * b10 + b02 * b09 + b03 * b08 - b04 * b07 + b05 * b06;
620
621    double invdet = 1.0 / det;
622    // If det is zero, we want to return false. However, we also want to return false
623    // if 1/det overflows to infinity (i.e. det is denormalized). Both of these are
624    // handled by checking that 1/det is finite.
625    if (!sk_float_isfinite(invdet)) {
626        return false;
627    }
628
629    b00 *= invdet;
630    b01 *= invdet;
631    b02 *= invdet;
632    b03 *= invdet;
633    b04 *= invdet;
634    b05 *= invdet;
635    b06 *= invdet;
636    b07 *= invdet;
637    b08 *= invdet;
638    b09 *= invdet;
639    b10 *= invdet;
640    b11 *= invdet;
641
642    inverse->fMat[0][0] = SkDoubleToMScalar(a11 * b11 - a12 * b10 + a13 * b09);
643    inverse->fMat[0][1] = SkDoubleToMScalar(a02 * b10 - a01 * b11 - a03 * b09);
644    inverse->fMat[0][2] = SkDoubleToMScalar(a31 * b05 - a32 * b04 + a33 * b03);
645    inverse->fMat[0][3] = SkDoubleToMScalar(a22 * b04 - a21 * b05 - a23 * b03);
646    inverse->fMat[1][0] = SkDoubleToMScalar(a12 * b08 - a10 * b11 - a13 * b07);
647    inverse->fMat[1][1] = SkDoubleToMScalar(a00 * b11 - a02 * b08 + a03 * b07);
648    inverse->fMat[1][2] = SkDoubleToMScalar(a32 * b02 - a30 * b05 - a33 * b01);
649    inverse->fMat[1][3] = SkDoubleToMScalar(a20 * b05 - a22 * b02 + a23 * b01);
650    inverse->fMat[2][0] = SkDoubleToMScalar(a10 * b10 - a11 * b08 + a13 * b06);
651    inverse->fMat[2][1] = SkDoubleToMScalar(a01 * b08 - a00 * b10 - a03 * b06);
652    inverse->fMat[2][2] = SkDoubleToMScalar(a30 * b04 - a31 * b02 + a33 * b00);
653    inverse->fMat[2][3] = SkDoubleToMScalar(a21 * b02 - a20 * b04 - a23 * b00);
654    inverse->fMat[3][0] = SkDoubleToMScalar(a11 * b07 - a10 * b09 - a12 * b06);
655    inverse->fMat[3][1] = SkDoubleToMScalar(a00 * b09 - a01 * b07 + a02 * b06);
656    inverse->fMat[3][2] = SkDoubleToMScalar(a31 * b01 - a30 * b03 - a32 * b00);
657    inverse->fMat[3][3] = SkDoubleToMScalar(a20 * b03 - a21 * b01 + a22 * b00);
658    inverse->dirtyTypeMask();
659
660    inverse->setTypeMask(this->getType());
661    if (!is_matrix_finite(*inverse)) {
662        return false;
663    }
664    if (storage && inverse != storage) {
665        *storage = *inverse;
666    }
667    return true;
668}
669
670///////////////////////////////////////////////////////////////////////////////
671
672void SkMatrix44::transpose() {
673    SkTSwap(fMat[0][1], fMat[1][0]);
674    SkTSwap(fMat[0][2], fMat[2][0]);
675    SkTSwap(fMat[0][3], fMat[3][0]);
676    SkTSwap(fMat[1][2], fMat[2][1]);
677    SkTSwap(fMat[1][3], fMat[3][1]);
678    SkTSwap(fMat[2][3], fMat[3][2]);
679
680    if (!this->isTriviallyIdentity()) {
681        this->dirtyTypeMask();
682    }
683}
684
685///////////////////////////////////////////////////////////////////////////////
686
687void SkMatrix44::mapScalars(const SkScalar src[4], SkScalar dst[4]) const {
688    SkScalar storage[4];
689    SkScalar* result = (src == dst) ? storage : dst;
690
691    for (int i = 0; i < 4; i++) {
692        SkMScalar value = 0;
693        for (int j = 0; j < 4; j++) {
694            value += fMat[j][i] * src[j];
695        }
696        result[i] = SkMScalarToScalar(value);
697    }
698
699    if (storage == result) {
700        memcpy(dst, storage, sizeof(storage));
701    }
702}
703
704#ifdef SK_MSCALAR_IS_DOUBLE
705
706void SkMatrix44::mapMScalars(const SkMScalar src[4], SkMScalar dst[4]) const {
707    SkMScalar storage[4];
708    SkMScalar* result = (src == dst) ? storage : dst;
709
710    for (int i = 0; i < 4; i++) {
711        SkMScalar value = 0;
712        for (int j = 0; j < 4; j++) {
713            value += fMat[j][i] * src[j];
714        }
715        result[i] = value;
716    }
717
718    if (storage == result) {
719        memcpy(dst, storage, sizeof(storage));
720    }
721}
722
723#endif
724
725typedef void (*Map2Procf)(const SkMScalar mat[][4], const float src2[], int count, float dst4[]);
726typedef void (*Map2Procd)(const SkMScalar mat[][4], const double src2[], int count, double dst4[]);
727
728static void map2_if(const SkMScalar mat[][4], const float* SK_RESTRICT src2,
729                    int count, float* SK_RESTRICT dst4) {
730    for (int i = 0; i < count; ++i) {
731        dst4[0] = src2[0];
732        dst4[1] = src2[1];
733        dst4[2] = 0;
734        dst4[3] = 1;
735        src2 += 2;
736        dst4 += 4;
737    }
738}
739
740static void map2_id(const SkMScalar mat[][4], const double* SK_RESTRICT src2,
741                    int count, double* SK_RESTRICT dst4) {
742    for (int i = 0; i < count; ++i) {
743        dst4[0] = src2[0];
744        dst4[1] = src2[1];
745        dst4[2] = 0;
746        dst4[3] = 1;
747        src2 += 2;
748        dst4 += 4;
749    }
750}
751
752static void map2_tf(const SkMScalar mat[][4], const float* SK_RESTRICT src2,
753                    int count, float* SK_RESTRICT dst4) {
754    const float mat30 = SkMScalarToFloat(mat[3][0]);
755    const float mat31 = SkMScalarToFloat(mat[3][1]);
756    const float mat32 = SkMScalarToFloat(mat[3][2]);
757    for (int n = 0; n < count; ++n) {
758        dst4[0] = src2[0] + mat30;
759        dst4[1] = src2[1] + mat31;
760        dst4[2] = mat32;
761        dst4[3] = 1;
762        src2 += 2;
763        dst4 += 4;
764    }
765}
766
767static void map2_td(const SkMScalar mat[][4], const double* SK_RESTRICT src2,
768                    int count, double* SK_RESTRICT dst4) {
769    for (int n = 0; n < count; ++n) {
770        dst4[0] = src2[0] + mat[3][0];
771        dst4[1] = src2[1] + mat[3][1];
772        dst4[2] = mat[3][2];
773        dst4[3] = 1;
774        src2 += 2;
775        dst4 += 4;
776    }
777}
778
779static void map2_sf(const SkMScalar mat[][4], const float* SK_RESTRICT src2,
780                    int count, float* SK_RESTRICT dst4) {
781    const float mat32 = SkMScalarToFloat(mat[3][2]);
782    for (int n = 0; n < count; ++n) {
783        dst4[0] = SkMScalarToFloat(mat[0][0] * src2[0] + mat[3][0]);
784        dst4[1] = SkMScalarToFloat(mat[1][1] * src2[1] + mat[3][1]);
785        dst4[2] = mat32;
786        dst4[3] = 1;
787        src2 += 2;
788        dst4 += 4;
789    }
790}
791
792static void map2_sd(const SkMScalar mat[][4], const double* SK_RESTRICT src2,
793                    int count, double* SK_RESTRICT dst4) {
794    for (int n = 0; n < count; ++n) {
795        dst4[0] = mat[0][0] * src2[0] + mat[3][0];
796        dst4[1] = mat[1][1] * src2[1] + mat[3][1];
797        dst4[2] = mat[3][2];
798        dst4[3] = 1;
799        src2 += 2;
800        dst4 += 4;
801    }
802}
803
804static void map2_af(const SkMScalar mat[][4], const float* SK_RESTRICT src2,
805                    int count, float* SK_RESTRICT dst4) {
806    SkMScalar r;
807    for (int n = 0; n < count; ++n) {
808        SkMScalar sx = SkFloatToMScalar(src2[0]);
809        SkMScalar sy = SkFloatToMScalar(src2[1]);
810        r = mat[0][0] * sx + mat[1][0] * sy + mat[3][0];
811        dst4[0] = SkMScalarToFloat(r);
812        r = mat[0][1] * sx + mat[1][1] * sy + mat[3][1];
813        dst4[1] = SkMScalarToFloat(r);
814        r = mat[0][2] * sx + mat[1][2] * sy + mat[3][2];
815        dst4[2] = SkMScalarToFloat(r);
816        dst4[3] = 1;
817        src2 += 2;
818        dst4 += 4;
819    }
820}
821
822static void map2_ad(const SkMScalar mat[][4], const double* SK_RESTRICT src2,
823                    int count, double* SK_RESTRICT dst4) {
824    for (int n = 0; n < count; ++n) {
825        double sx = src2[0];
826        double sy = src2[1];
827        dst4[0] = mat[0][0] * sx + mat[1][0] * sy + mat[3][0];
828        dst4[1] = mat[0][1] * sx + mat[1][1] * sy + mat[3][1];
829        dst4[2] = mat[0][2] * sx + mat[1][2] * sy + mat[3][2];
830        dst4[3] = 1;
831        src2 += 2;
832        dst4 += 4;
833    }
834}
835
836static void map2_pf(const SkMScalar mat[][4], const float* SK_RESTRICT src2,
837                    int count, float* SK_RESTRICT dst4) {
838    SkMScalar r;
839    for (int n = 0; n < count; ++n) {
840        SkMScalar sx = SkFloatToMScalar(src2[0]);
841        SkMScalar sy = SkFloatToMScalar(src2[1]);
842        for (int i = 0; i < 4; i++) {
843            r = mat[0][i] * sx + mat[1][i] * sy + mat[3][i];
844            dst4[i] = SkMScalarToFloat(r);
845        }
846        src2 += 2;
847        dst4 += 4;
848    }
849}
850
851static void map2_pd(const SkMScalar mat[][4], const double* SK_RESTRICT src2,
852                    int count, double* SK_RESTRICT dst4) {
853    for (int n = 0; n < count; ++n) {
854        double sx = src2[0];
855        double sy = src2[1];
856        for (int i = 0; i < 4; i++) {
857            dst4[i] = mat[0][i] * sx + mat[1][i] * sy + mat[3][i];
858        }
859        src2 += 2;
860        dst4 += 4;
861    }
862}
863
864void SkMatrix44::map2(const float src2[], int count, float dst4[]) const {
865    static const Map2Procf gProc[] = {
866        map2_if, map2_tf, map2_sf, map2_sf, map2_af, map2_af, map2_af, map2_af
867    };
868
869    TypeMask mask = this->getType();
870    Map2Procf proc = (mask & kPerspective_Mask) ? map2_pf : gProc[mask];
871    proc(fMat, src2, count, dst4);
872}
873
874void SkMatrix44::map2(const double src2[], int count, double dst4[]) const {
875    static const Map2Procd gProc[] = {
876        map2_id, map2_td, map2_sd, map2_sd, map2_ad, map2_ad, map2_ad, map2_ad
877    };
878
879    TypeMask mask = this->getType();
880    Map2Procd proc = (mask & kPerspective_Mask) ? map2_pd : gProc[mask];
881    proc(fMat, src2, count, dst4);
882}
883
884bool SkMatrix44::preserves2dAxisAlignment (SkMScalar epsilon) const {
885
886    // Can't check (mask & kPerspective_Mask) because Z isn't relevant here.
887    if (0 != perspX() || 0 != perspY()) return false;
888
889    // A matrix with two non-zeroish values in any of the upper right
890    // rows or columns will skew.  If only one value in each row or
891    // column is non-zeroish, we get a scale plus perhaps a 90-degree
892    // rotation.
893    int col0 = 0;
894    int col1 = 0;
895    int row0 = 0;
896    int row1 = 0;
897
898    // Must test against epsilon, not 0, because we can get values
899    // around 6e-17 in the matrix that "should" be 0.
900
901    if (SkMScalarAbs(fMat[0][0]) > epsilon) {
902        col0++;
903        row0++;
904    }
905    if (SkMScalarAbs(fMat[0][1]) > epsilon) {
906        col1++;
907        row0++;
908    }
909    if (SkMScalarAbs(fMat[1][0]) > epsilon) {
910        col0++;
911        row1++;
912    }
913    if (SkMScalarAbs(fMat[1][1]) > epsilon) {
914        col1++;
915        row1++;
916    }
917    if (col0 > 1 || col1 > 1 || row0 > 1 || row1 > 1) {
918        return false;
919    }
920
921    return true;
922}
923
924///////////////////////////////////////////////////////////////////////////////
925
926void SkMatrix44::dump() const {
927    static const char* format =
928        "[%g %g %g %g][%g %g %g %g][%g %g %g %g][%g %g %g %g]\n";
929#if 0
930    SkDebugf(format,
931             fMat[0][0], fMat[1][0], fMat[2][0], fMat[3][0],
932             fMat[0][1], fMat[1][1], fMat[2][1], fMat[3][1],
933             fMat[0][2], fMat[1][2], fMat[2][2], fMat[3][2],
934             fMat[0][3], fMat[1][3], fMat[2][3], fMat[3][3]);
935#else
936    SkDebugf(format,
937             fMat[0][0], fMat[0][1], fMat[0][2], fMat[0][3],
938             fMat[1][0], fMat[1][1], fMat[1][2], fMat[1][3],
939             fMat[2][0], fMat[2][1], fMat[2][2], fMat[2][3],
940             fMat[3][0], fMat[3][1], fMat[3][2], fMat[3][3]);
941#endif
942}
943
944///////////////////////////////////////////////////////////////////////////////
945
946static void initFromMatrix(SkMScalar dst[4][4], const SkMatrix& src) {
947    dst[0][0] = SkScalarToMScalar(src[SkMatrix::kMScaleX]);
948    dst[1][0] = SkScalarToMScalar(src[SkMatrix::kMSkewX]);
949    dst[2][0] = 0;
950    dst[3][0] = SkScalarToMScalar(src[SkMatrix::kMTransX]);
951    dst[0][1] = SkScalarToMScalar(src[SkMatrix::kMSkewY]);
952    dst[1][1] = SkScalarToMScalar(src[SkMatrix::kMScaleY]);
953    dst[2][1] = 0;
954    dst[3][1] = SkScalarToMScalar(src[SkMatrix::kMTransY]);
955    dst[0][2] = 0;
956    dst[1][2] = 0;
957    dst[2][2] = 1;
958    dst[3][2] = 0;
959    dst[0][3] = SkScalarToMScalar(src[SkMatrix::kMPersp0]);
960    dst[1][3] = SkScalarToMScalar(src[SkMatrix::kMPersp1]);
961    dst[2][3] = 0;
962    dst[3][3] = SkScalarToMScalar(src[SkMatrix::kMPersp2]);
963}
964
965SkMatrix44::SkMatrix44(const SkMatrix& src) {
966    this->operator=(src);
967}
968
969SkMatrix44& SkMatrix44::operator=(const SkMatrix& src) {
970    initFromMatrix(fMat, src);
971
972    if (src.isIdentity()) {
973        this->setTypeMask(kIdentity_Mask);
974    } else {
975        this->dirtyTypeMask();
976    }
977    return *this;
978}
979
980SkMatrix44::operator SkMatrix() const {
981    SkMatrix dst;
982
983    dst[SkMatrix::kMScaleX]  = SkMScalarToScalar(fMat[0][0]);
984    dst[SkMatrix::kMSkewX]  = SkMScalarToScalar(fMat[1][0]);
985    dst[SkMatrix::kMTransX] = SkMScalarToScalar(fMat[3][0]);
986
987    dst[SkMatrix::kMSkewY]  = SkMScalarToScalar(fMat[0][1]);
988    dst[SkMatrix::kMScaleY] = SkMScalarToScalar(fMat[1][1]);
989    dst[SkMatrix::kMTransY] = SkMScalarToScalar(fMat[3][1]);
990
991    dst[SkMatrix::kMPersp0] = SkMScalarToScalar(fMat[0][3]);
992    dst[SkMatrix::kMPersp1] = SkMScalarToScalar(fMat[1][3]);
993    dst[SkMatrix::kMPersp2] = SkMScalarToScalar(fMat[3][3]);
994
995    return dst;
996}
997