1/*
2 * Copyright 2006 The Android Open Source Project
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 "SkMatrix.h"
9#include "SkFloatBits.h"
10#include "SkString.h"
11#include "SkNx.h"
12
13#include <stddef.h>
14
15static void normalize_perspective(SkScalar mat[9]) {
16    // If it was interesting to never store the last element, we could divide all 8 other
17    // elements here by the 9th, making it 1.0...
18    //
19    // When SkScalar was SkFixed, we would sometimes rescale the entire matrix to keep its
20    // component values from getting too large. This is not a concern when using floats/doubles,
21    // so we do nothing now.
22
23    // Disable this for now, but it could be enabled.
24#if 0
25    if (0 == mat[SkMatrix::kMPersp0] && 0 == mat[SkMatrix::kMPersp1]) {
26        SkScalar p2 = mat[SkMatrix::kMPersp2];
27        if (p2 != 0 && p2 != 1) {
28            double inv = 1.0 / p2;
29            for (int i = 0; i < 6; ++i) {
30                mat[i] = SkDoubleToScalar(mat[i] * inv);
31            }
32            mat[SkMatrix::kMPersp2] = 1;
33        }
34    }
35#endif
36}
37
38// In a few places, we performed the following
39//      a * b + c * d + e
40// as
41//      a * b + (c * d + e)
42//
43// sdot and scross are indended to capture these compound operations into a
44// function, with an eye toward considering upscaling the intermediates to
45// doubles for more precision (as we do in concat and invert).
46//
47// However, these few lines that performed the last add before the "dot", cause
48// tiny image differences, so we guard that change until we see the impact on
49// chrome's layouttests.
50//
51#define SK_LEGACY_MATRIX_MATH_ORDER
52
53static inline float SkDoubleToFloat(double x) {
54    return static_cast<float>(x);
55}
56
57/*      [scale-x    skew-x      trans-x]   [X]   [X']
58        [skew-y     scale-y     trans-y] * [Y] = [Y']
59        [persp-0    persp-1     persp-2]   [1]   [1 ]
60*/
61
62void SkMatrix::reset() {
63    fMat[kMScaleX] = fMat[kMScaleY] = fMat[kMPersp2] = 1;
64    fMat[kMSkewX]  = fMat[kMSkewY] =
65    fMat[kMTransX] = fMat[kMTransY] =
66    fMat[kMPersp0] = fMat[kMPersp1] = 0;
67    this->setTypeMask(kIdentity_Mask | kRectStaysRect_Mask);
68}
69
70void SkMatrix::set9(const SkScalar buffer[]) {
71    memcpy(fMat, buffer, 9 * sizeof(SkScalar));
72    normalize_perspective(fMat);
73    this->setTypeMask(kUnknown_Mask);
74}
75
76void SkMatrix::setAffine(const SkScalar buffer[]) {
77    fMat[kMScaleX] = buffer[kAScaleX];
78    fMat[kMSkewX]  = buffer[kASkewX];
79    fMat[kMTransX] = buffer[kATransX];
80    fMat[kMSkewY]  = buffer[kASkewY];
81    fMat[kMScaleY] = buffer[kAScaleY];
82    fMat[kMTransY] = buffer[kATransY];
83    fMat[kMPersp0] = 0;
84    fMat[kMPersp1] = 0;
85    fMat[kMPersp2] = 1;
86    this->setTypeMask(kUnknown_Mask);
87}
88
89// this guy aligns with the masks, so we can compute a mask from a varaible 0/1
90enum {
91    kTranslate_Shift,
92    kScale_Shift,
93    kAffine_Shift,
94    kPerspective_Shift,
95    kRectStaysRect_Shift
96};
97
98static const int32_t kScalar1Int = 0x3f800000;
99
100uint8_t SkMatrix::computePerspectiveTypeMask() const {
101    // Benchmarking suggests that replacing this set of SkScalarAs2sCompliment
102    // is a win, but replacing those below is not. We don't yet understand
103    // that result.
104    if (fMat[kMPersp0] != 0 || fMat[kMPersp1] != 0 || fMat[kMPersp2] != 1) {
105        // If this is a perspective transform, we return true for all other
106        // transform flags - this does not disable any optimizations, respects
107        // the rule that the type mask must be conservative, and speeds up
108        // type mask computation.
109        return SkToU8(kORableMasks);
110    }
111
112    return SkToU8(kOnlyPerspectiveValid_Mask | kUnknown_Mask);
113}
114
115uint8_t SkMatrix::computeTypeMask() const {
116    unsigned mask = 0;
117
118    if (fMat[kMPersp0] != 0 || fMat[kMPersp1] != 0 || fMat[kMPersp2] != 1) {
119        // Once it is determined that that this is a perspective transform,
120        // all other flags are moot as far as optimizations are concerned.
121        return SkToU8(kORableMasks);
122    }
123
124    if (fMat[kMTransX] != 0 || fMat[kMTransY] != 0) {
125        mask |= kTranslate_Mask;
126    }
127
128    int m00 = SkScalarAs2sCompliment(fMat[SkMatrix::kMScaleX]);
129    int m01 = SkScalarAs2sCompliment(fMat[SkMatrix::kMSkewX]);
130    int m10 = SkScalarAs2sCompliment(fMat[SkMatrix::kMSkewY]);
131    int m11 = SkScalarAs2sCompliment(fMat[SkMatrix::kMScaleY]);
132
133    if (m01 | m10) {
134        // The skew components may be scale-inducing, unless we are dealing
135        // with a pure rotation.  Testing for a pure rotation is expensive,
136        // so we opt for being conservative by always setting the scale bit.
137        // along with affine.
138        // By doing this, we are also ensuring that matrices have the same
139        // type masks as their inverses.
140        mask |= kAffine_Mask | kScale_Mask;
141
142        // For rectStaysRect, in the affine case, we only need check that
143        // the primary diagonal is all zeros and that the secondary diagonal
144        // is all non-zero.
145
146        // map non-zero to 1
147        m01 = m01 != 0;
148        m10 = m10 != 0;
149
150        int dp0 = 0 == (m00 | m11) ;  // true if both are 0
151        int ds1 = m01 & m10;        // true if both are 1
152
153        mask |= (dp0 & ds1) << kRectStaysRect_Shift;
154    } else {
155        // Only test for scale explicitly if not affine, since affine sets the
156        // scale bit.
157        if ((m00 - kScalar1Int) | (m11 - kScalar1Int)) {
158            mask |= kScale_Mask;
159        }
160
161        // Not affine, therefore we already know secondary diagonal is
162        // all zeros, so we just need to check that primary diagonal is
163        // all non-zero.
164
165        // map non-zero to 1
166        m00 = m00 != 0;
167        m11 = m11 != 0;
168
169        // record if the (p)rimary diagonal is all non-zero
170        mask |= (m00 & m11) << kRectStaysRect_Shift;
171    }
172
173    return SkToU8(mask);
174}
175
176///////////////////////////////////////////////////////////////////////////////
177
178bool operator==(const SkMatrix& a, const SkMatrix& b) {
179    const SkScalar* SK_RESTRICT ma = a.fMat;
180    const SkScalar* SK_RESTRICT mb = b.fMat;
181
182    return  ma[0] == mb[0] && ma[1] == mb[1] && ma[2] == mb[2] &&
183            ma[3] == mb[3] && ma[4] == mb[4] && ma[5] == mb[5] &&
184            ma[6] == mb[6] && ma[7] == mb[7] && ma[8] == mb[8];
185}
186
187///////////////////////////////////////////////////////////////////////////////
188
189// helper function to determine if upper-left 2x2 of matrix is degenerate
190static inline bool is_degenerate_2x2(SkScalar scaleX, SkScalar skewX,
191                                     SkScalar skewY,  SkScalar scaleY) {
192    SkScalar perp_dot = scaleX*scaleY - skewX*skewY;
193    return SkScalarNearlyZero(perp_dot, SK_ScalarNearlyZero*SK_ScalarNearlyZero);
194}
195
196///////////////////////////////////////////////////////////////////////////////
197
198bool SkMatrix::isSimilarity(SkScalar tol) const {
199    // if identity or translate matrix
200    TypeMask mask = this->getType();
201    if (mask <= kTranslate_Mask) {
202        return true;
203    }
204    if (mask & kPerspective_Mask) {
205        return false;
206    }
207
208    SkScalar mx = fMat[kMScaleX];
209    SkScalar my = fMat[kMScaleY];
210    // if no skew, can just compare scale factors
211    if (!(mask & kAffine_Mask)) {
212        return !SkScalarNearlyZero(mx) && SkScalarNearlyEqual(SkScalarAbs(mx), SkScalarAbs(my));
213    }
214    SkScalar sx = fMat[kMSkewX];
215    SkScalar sy = fMat[kMSkewY];
216
217    if (is_degenerate_2x2(mx, sx, sy, my)) {
218        return false;
219    }
220
221    // upper 2x2 is rotation/reflection + uniform scale if basis vectors
222    // are 90 degree rotations of each other
223    return (SkScalarNearlyEqual(mx, my, tol) && SkScalarNearlyEqual(sx, -sy, tol))
224        || (SkScalarNearlyEqual(mx, -my, tol) && SkScalarNearlyEqual(sx, sy, tol));
225}
226
227bool SkMatrix::preservesRightAngles(SkScalar tol) const {
228    TypeMask mask = this->getType();
229
230    if (mask <= kTranslate_Mask) {
231        // identity, translate and/or scale
232        return true;
233    }
234    if (mask & kPerspective_Mask) {
235        return false;
236    }
237
238    SkASSERT(mask & (kAffine_Mask | kScale_Mask));
239
240    SkScalar mx = fMat[kMScaleX];
241    SkScalar my = fMat[kMScaleY];
242    SkScalar sx = fMat[kMSkewX];
243    SkScalar sy = fMat[kMSkewY];
244
245    if (is_degenerate_2x2(mx, sx, sy, my)) {
246        return false;
247    }
248
249    // upper 2x2 is scale + rotation/reflection if basis vectors are orthogonal
250    SkVector vec[2];
251    vec[0].set(mx, sy);
252    vec[1].set(sx, my);
253
254    return SkScalarNearlyZero(vec[0].dot(vec[1]), SkScalarSquare(tol));
255}
256
257///////////////////////////////////////////////////////////////////////////////
258
259static inline SkScalar sdot(SkScalar a, SkScalar b, SkScalar c, SkScalar d) {
260    return a * b + c * d;
261}
262
263static inline SkScalar sdot(SkScalar a, SkScalar b, SkScalar c, SkScalar d,
264                             SkScalar e, SkScalar f) {
265    return a * b + c * d + e * f;
266}
267
268static inline SkScalar scross(SkScalar a, SkScalar b, SkScalar c, SkScalar d) {
269    return a * b - c * d;
270}
271
272void SkMatrix::setTranslate(SkScalar dx, SkScalar dy) {
273    if (dx || dy) {
274        fMat[kMTransX] = dx;
275        fMat[kMTransY] = dy;
276
277        fMat[kMScaleX] = fMat[kMScaleY] = fMat[kMPersp2] = 1;
278        fMat[kMSkewX]  = fMat[kMSkewY] =
279        fMat[kMPersp0] = fMat[kMPersp1] = 0;
280
281        this->setTypeMask(kTranslate_Mask | kRectStaysRect_Mask);
282    } else {
283        this->reset();
284    }
285}
286
287void SkMatrix::preTranslate(SkScalar dx, SkScalar dy) {
288    if (!dx && !dy) {
289        return;
290    }
291
292    if (this->hasPerspective()) {
293        SkMatrix    m;
294        m.setTranslate(dx, dy);
295        this->preConcat(m);
296    } else {
297        fMat[kMTransX] += sdot(fMat[kMScaleX], dx, fMat[kMSkewX], dy);
298        fMat[kMTransY] += sdot(fMat[kMSkewY], dx, fMat[kMScaleY], dy);
299        this->setTypeMask(kUnknown_Mask | kOnlyPerspectiveValid_Mask);
300    }
301}
302
303void SkMatrix::postTranslate(SkScalar dx, SkScalar dy) {
304    if (!dx && !dy) {
305        return;
306    }
307
308    if (this->hasPerspective()) {
309        SkMatrix    m;
310        m.setTranslate(dx, dy);
311        this->postConcat(m);
312    } else {
313        fMat[kMTransX] += dx;
314        fMat[kMTransY] += dy;
315        this->setTypeMask(kUnknown_Mask | kOnlyPerspectiveValid_Mask);
316    }
317}
318
319///////////////////////////////////////////////////////////////////////////////
320
321void SkMatrix::setScale(SkScalar sx, SkScalar sy, SkScalar px, SkScalar py) {
322    if (1 == sx && 1 == sy) {
323        this->reset();
324    } else {
325        this->setScaleTranslate(sx, sy, px - sx * px, py - sy * py);
326    }
327}
328
329void SkMatrix::setScale(SkScalar sx, SkScalar sy) {
330    if (1 == sx && 1 == sy) {
331        this->reset();
332    } else {
333        fMat[kMScaleX] = sx;
334        fMat[kMScaleY] = sy;
335        fMat[kMPersp2] = 1;
336
337        fMat[kMTransX] = fMat[kMTransY] =
338        fMat[kMSkewX]  = fMat[kMSkewY] =
339        fMat[kMPersp0] = fMat[kMPersp1] = 0;
340
341        this->setTypeMask(kScale_Mask | kRectStaysRect_Mask);
342    }
343}
344
345bool SkMatrix::setIDiv(int divx, int divy) {
346    if (!divx || !divy) {
347        return false;
348    }
349    this->setScale(SkScalarInvert(divx), SkScalarInvert(divy));
350    return true;
351}
352
353void SkMatrix::preScale(SkScalar sx, SkScalar sy, SkScalar px, SkScalar py) {
354    if (1 == sx && 1 == sy) {
355        return;
356    }
357
358    SkMatrix    m;
359    m.setScale(sx, sy, px, py);
360    this->preConcat(m);
361}
362
363void SkMatrix::preScale(SkScalar sx, SkScalar sy) {
364    if (1 == sx && 1 == sy) {
365        return;
366    }
367
368    // the assumption is that these multiplies are very cheap, and that
369    // a full concat and/or just computing the matrix type is more expensive.
370    // Also, the fixed-point case checks for overflow, but the float doesn't,
371    // so we can get away with these blind multiplies.
372
373    fMat[kMScaleX] *= sx;
374    fMat[kMSkewY]  *= sx;
375    fMat[kMPersp0] *= sx;
376
377    fMat[kMSkewX]  *= sy;
378    fMat[kMScaleY] *= sy;
379    fMat[kMPersp1] *= sy;
380
381    this->orTypeMask(kScale_Mask);
382}
383
384void SkMatrix::postScale(SkScalar sx, SkScalar sy, SkScalar px, SkScalar py) {
385    if (1 == sx && 1 == sy) {
386        return;
387    }
388    SkMatrix    m;
389    m.setScale(sx, sy, px, py);
390    this->postConcat(m);
391}
392
393void SkMatrix::postScale(SkScalar sx, SkScalar sy) {
394    if (1 == sx && 1 == sy) {
395        return;
396    }
397    SkMatrix    m;
398    m.setScale(sx, sy);
399    this->postConcat(m);
400}
401
402// this guy perhaps can go away, if we have a fract/high-precision way to
403// scale matrices
404bool SkMatrix::postIDiv(int divx, int divy) {
405    if (divx == 0 || divy == 0) {
406        return false;
407    }
408
409    const float invX = 1.f / divx;
410    const float invY = 1.f / divy;
411
412    fMat[kMScaleX] *= invX;
413    fMat[kMSkewX]  *= invX;
414    fMat[kMTransX] *= invX;
415
416    fMat[kMScaleY] *= invY;
417    fMat[kMSkewY]  *= invY;
418    fMat[kMTransY] *= invY;
419
420    this->setTypeMask(kUnknown_Mask);
421    return true;
422}
423
424////////////////////////////////////////////////////////////////////////////////////
425
426void SkMatrix::setSinCos(SkScalar sinV, SkScalar cosV, SkScalar px, SkScalar py) {
427    const SkScalar oneMinusCosV = 1 - cosV;
428
429    fMat[kMScaleX]  = cosV;
430    fMat[kMSkewX]   = -sinV;
431    fMat[kMTransX]  = sdot(sinV, py, oneMinusCosV, px);
432
433    fMat[kMSkewY]   = sinV;
434    fMat[kMScaleY]  = cosV;
435    fMat[kMTransY]  = sdot(-sinV, px, oneMinusCosV, py);
436
437    fMat[kMPersp0] = fMat[kMPersp1] = 0;
438    fMat[kMPersp2] = 1;
439
440    this->setTypeMask(kUnknown_Mask | kOnlyPerspectiveValid_Mask);
441}
442
443void SkMatrix::setSinCos(SkScalar sinV, SkScalar cosV) {
444    fMat[kMScaleX]  = cosV;
445    fMat[kMSkewX]   = -sinV;
446    fMat[kMTransX]  = 0;
447
448    fMat[kMSkewY]   = sinV;
449    fMat[kMScaleY]  = cosV;
450    fMat[kMTransY]  = 0;
451
452    fMat[kMPersp0] = fMat[kMPersp1] = 0;
453    fMat[kMPersp2] = 1;
454
455    this->setTypeMask(kUnknown_Mask | kOnlyPerspectiveValid_Mask);
456}
457
458void SkMatrix::setRotate(SkScalar degrees, SkScalar px, SkScalar py) {
459    SkScalar sinV, cosV;
460    sinV = SkScalarSinCos(SkDegreesToRadians(degrees), &cosV);
461    this->setSinCos(sinV, cosV, px, py);
462}
463
464void SkMatrix::setRotate(SkScalar degrees) {
465    SkScalar sinV, cosV;
466    sinV = SkScalarSinCos(SkDegreesToRadians(degrees), &cosV);
467    this->setSinCos(sinV, cosV);
468}
469
470void SkMatrix::preRotate(SkScalar degrees, SkScalar px, SkScalar py) {
471    SkMatrix    m;
472    m.setRotate(degrees, px, py);
473    this->preConcat(m);
474}
475
476void SkMatrix::preRotate(SkScalar degrees) {
477    SkMatrix    m;
478    m.setRotate(degrees);
479    this->preConcat(m);
480}
481
482void SkMatrix::postRotate(SkScalar degrees, SkScalar px, SkScalar py) {
483    SkMatrix    m;
484    m.setRotate(degrees, px, py);
485    this->postConcat(m);
486}
487
488void SkMatrix::postRotate(SkScalar degrees) {
489    SkMatrix    m;
490    m.setRotate(degrees);
491    this->postConcat(m);
492}
493
494////////////////////////////////////////////////////////////////////////////////////
495
496void SkMatrix::setSkew(SkScalar sx, SkScalar sy, SkScalar px, SkScalar py) {
497    fMat[kMScaleX]  = 1;
498    fMat[kMSkewX]   = sx;
499    fMat[kMTransX]  = -sx * py;
500
501    fMat[kMSkewY]   = sy;
502    fMat[kMScaleY]  = 1;
503    fMat[kMTransY]  = -sy * px;
504
505    fMat[kMPersp0] = fMat[kMPersp1] = 0;
506    fMat[kMPersp2] = 1;
507
508    this->setTypeMask(kUnknown_Mask | kOnlyPerspectiveValid_Mask);
509}
510
511void SkMatrix::setSkew(SkScalar sx, SkScalar sy) {
512    fMat[kMScaleX]  = 1;
513    fMat[kMSkewX]   = sx;
514    fMat[kMTransX]  = 0;
515
516    fMat[kMSkewY]   = sy;
517    fMat[kMScaleY]  = 1;
518    fMat[kMTransY]  = 0;
519
520    fMat[kMPersp0] = fMat[kMPersp1] = 0;
521    fMat[kMPersp2] = 1;
522
523    this->setTypeMask(kUnknown_Mask | kOnlyPerspectiveValid_Mask);
524}
525
526void SkMatrix::preSkew(SkScalar sx, SkScalar sy, SkScalar px, SkScalar py) {
527    SkMatrix    m;
528    m.setSkew(sx, sy, px, py);
529    this->preConcat(m);
530}
531
532void SkMatrix::preSkew(SkScalar sx, SkScalar sy) {
533    SkMatrix    m;
534    m.setSkew(sx, sy);
535    this->preConcat(m);
536}
537
538void SkMatrix::postSkew(SkScalar sx, SkScalar sy, SkScalar px, SkScalar py) {
539    SkMatrix    m;
540    m.setSkew(sx, sy, px, py);
541    this->postConcat(m);
542}
543
544void SkMatrix::postSkew(SkScalar sx, SkScalar sy) {
545    SkMatrix    m;
546    m.setSkew(sx, sy);
547    this->postConcat(m);
548}
549
550///////////////////////////////////////////////////////////////////////////////
551
552bool SkMatrix::setRectToRect(const SkRect& src, const SkRect& dst, ScaleToFit align) {
553    if (src.isEmpty()) {
554        this->reset();
555        return false;
556    }
557
558    if (dst.isEmpty()) {
559        sk_bzero(fMat, 8 * sizeof(SkScalar));
560        fMat[kMPersp2] = 1;
561        this->setTypeMask(kScale_Mask | kRectStaysRect_Mask);
562    } else {
563        SkScalar    tx, sx = dst.width() / src.width();
564        SkScalar    ty, sy = dst.height() / src.height();
565        bool        xLarger = false;
566
567        if (align != kFill_ScaleToFit) {
568            if (sx > sy) {
569                xLarger = true;
570                sx = sy;
571            } else {
572                sy = sx;
573            }
574        }
575
576        tx = dst.fLeft - src.fLeft * sx;
577        ty = dst.fTop - src.fTop * sy;
578        if (align == kCenter_ScaleToFit || align == kEnd_ScaleToFit) {
579            SkScalar diff;
580
581            if (xLarger) {
582                diff = dst.width() - src.width() * sy;
583            } else {
584                diff = dst.height() - src.height() * sy;
585            }
586
587            if (align == kCenter_ScaleToFit) {
588                diff = SkScalarHalf(diff);
589            }
590
591            if (xLarger) {
592                tx += diff;
593            } else {
594                ty += diff;
595            }
596        }
597
598        this->setScaleTranslate(sx, sy, tx, ty);
599    }
600    return true;
601}
602
603///////////////////////////////////////////////////////////////////////////////
604
605static inline float muladdmul(float a, float b, float c, float d) {
606    return SkDoubleToFloat((double)a * b + (double)c * d);
607}
608
609static inline float rowcol3(const float row[], const float col[]) {
610    return row[0] * col[0] + row[1] * col[3] + row[2] * col[6];
611}
612
613static bool only_scale_and_translate(unsigned mask) {
614    return 0 == (mask & (SkMatrix::kAffine_Mask | SkMatrix::kPerspective_Mask));
615}
616
617void SkMatrix::setConcat(const SkMatrix& a, const SkMatrix& b) {
618    TypeMask aType = a.getType();
619    TypeMask bType = b.getType();
620
621    if (a.isTriviallyIdentity()) {
622        *this = b;
623    } else if (b.isTriviallyIdentity()) {
624        *this = a;
625    } else if (only_scale_and_translate(aType | bType)) {
626        this->setScaleTranslate(a.fMat[kMScaleX] * b.fMat[kMScaleX],
627                                a.fMat[kMScaleY] * b.fMat[kMScaleY],
628                                a.fMat[kMScaleX] * b.fMat[kMTransX] + a.fMat[kMTransX],
629                                a.fMat[kMScaleY] * b.fMat[kMTransY] + a.fMat[kMTransY]);
630    } else {
631        SkMatrix tmp;
632
633        if ((aType | bType) & kPerspective_Mask) {
634            tmp.fMat[kMScaleX] = rowcol3(&a.fMat[0], &b.fMat[0]);
635            tmp.fMat[kMSkewX]  = rowcol3(&a.fMat[0], &b.fMat[1]);
636            tmp.fMat[kMTransX] = rowcol3(&a.fMat[0], &b.fMat[2]);
637            tmp.fMat[kMSkewY]  = rowcol3(&a.fMat[3], &b.fMat[0]);
638            tmp.fMat[kMScaleY] = rowcol3(&a.fMat[3], &b.fMat[1]);
639            tmp.fMat[kMTransY] = rowcol3(&a.fMat[3], &b.fMat[2]);
640            tmp.fMat[kMPersp0] = rowcol3(&a.fMat[6], &b.fMat[0]);
641            tmp.fMat[kMPersp1] = rowcol3(&a.fMat[6], &b.fMat[1]);
642            tmp.fMat[kMPersp2] = rowcol3(&a.fMat[6], &b.fMat[2]);
643
644            normalize_perspective(tmp.fMat);
645            tmp.setTypeMask(kUnknown_Mask);
646        } else {
647            tmp.fMat[kMScaleX] = muladdmul(a.fMat[kMScaleX],
648                                           b.fMat[kMScaleX],
649                                           a.fMat[kMSkewX],
650                                           b.fMat[kMSkewY]);
651
652            tmp.fMat[kMSkewX]  = muladdmul(a.fMat[kMScaleX],
653                                           b.fMat[kMSkewX],
654                                           a.fMat[kMSkewX],
655                                           b.fMat[kMScaleY]);
656
657            tmp.fMat[kMTransX] = muladdmul(a.fMat[kMScaleX],
658                                           b.fMat[kMTransX],
659                                           a.fMat[kMSkewX],
660                                           b.fMat[kMTransY]) + a.fMat[kMTransX];
661
662            tmp.fMat[kMSkewY]  = muladdmul(a.fMat[kMSkewY],
663                                           b.fMat[kMScaleX],
664                                           a.fMat[kMScaleY],
665                                           b.fMat[kMSkewY]);
666
667            tmp.fMat[kMScaleY] = muladdmul(a.fMat[kMSkewY],
668                                           b.fMat[kMSkewX],
669                                           a.fMat[kMScaleY],
670                                           b.fMat[kMScaleY]);
671
672            tmp.fMat[kMTransY] = muladdmul(a.fMat[kMSkewY],
673                                           b.fMat[kMTransX],
674                                           a.fMat[kMScaleY],
675                                           b.fMat[kMTransY]) + a.fMat[kMTransY];
676
677            tmp.fMat[kMPersp0] = 0;
678            tmp.fMat[kMPersp1] = 0;
679            tmp.fMat[kMPersp2] = 1;
680            //SkDebugf("Concat mat non-persp type: %d\n", tmp.getType());
681            //SkASSERT(!(tmp.getType() & kPerspective_Mask));
682            tmp.setTypeMask(kUnknown_Mask | kOnlyPerspectiveValid_Mask);
683        }
684        *this = tmp;
685    }
686}
687
688void SkMatrix::preConcat(const SkMatrix& mat) {
689    // check for identity first, so we don't do a needless copy of ourselves
690    // to ourselves inside setConcat()
691    if(!mat.isIdentity()) {
692        this->setConcat(*this, mat);
693    }
694}
695
696void SkMatrix::postConcat(const SkMatrix& mat) {
697    // check for identity first, so we don't do a needless copy of ourselves
698    // to ourselves inside setConcat()
699    if (!mat.isIdentity()) {
700        this->setConcat(mat, *this);
701    }
702}
703
704///////////////////////////////////////////////////////////////////////////////
705
706/*  Matrix inversion is very expensive, but also the place where keeping
707    precision may be most important (here and matrix concat). Hence to avoid
708    bitmap blitting artifacts when walking the inverse, we use doubles for
709    the intermediate math, even though we know that is more expensive.
710 */
711
712static inline SkScalar scross_dscale(SkScalar a, SkScalar b,
713                                     SkScalar c, SkScalar d, double scale) {
714    return SkDoubleToScalar(scross(a, b, c, d) * scale);
715}
716
717static inline double dcross(double a, double b, double c, double d) {
718    return a * b - c * d;
719}
720
721static inline SkScalar dcross_dscale(double a, double b,
722                                     double c, double d, double scale) {
723    return SkDoubleToScalar(dcross(a, b, c, d) * scale);
724}
725
726static double sk_inv_determinant(const float mat[9], int isPerspective) {
727    double det;
728
729    if (isPerspective) {
730        det = mat[SkMatrix::kMScaleX] *
731              dcross(mat[SkMatrix::kMScaleY], mat[SkMatrix::kMPersp2],
732                     mat[SkMatrix::kMTransY], mat[SkMatrix::kMPersp1])
733              +
734              mat[SkMatrix::kMSkewX]  *
735              dcross(mat[SkMatrix::kMTransY], mat[SkMatrix::kMPersp0],
736                     mat[SkMatrix::kMSkewY],  mat[SkMatrix::kMPersp2])
737              +
738              mat[SkMatrix::kMTransX] *
739              dcross(mat[SkMatrix::kMSkewY],  mat[SkMatrix::kMPersp1],
740                     mat[SkMatrix::kMScaleY], mat[SkMatrix::kMPersp0]);
741    } else {
742        det = dcross(mat[SkMatrix::kMScaleX], mat[SkMatrix::kMScaleY],
743                     mat[SkMatrix::kMSkewX], mat[SkMatrix::kMSkewY]);
744    }
745
746    // Since the determinant is on the order of the cube of the matrix members,
747    // compare to the cube of the default nearly-zero constant (although an
748    // estimate of the condition number would be better if it wasn't so expensive).
749    if (SkScalarNearlyZero((float)det, SK_ScalarNearlyZero * SK_ScalarNearlyZero * SK_ScalarNearlyZero)) {
750        return 0;
751    }
752    return 1.0 / det;
753}
754
755void SkMatrix::SetAffineIdentity(SkScalar affine[6]) {
756    affine[kAScaleX] = 1;
757    affine[kASkewY] = 0;
758    affine[kASkewX] = 0;
759    affine[kAScaleY] = 1;
760    affine[kATransX] = 0;
761    affine[kATransY] = 0;
762}
763
764bool SkMatrix::asAffine(SkScalar affine[6]) const {
765    if (this->hasPerspective()) {
766        return false;
767    }
768    if (affine) {
769        affine[kAScaleX] = this->fMat[kMScaleX];
770        affine[kASkewY] = this->fMat[kMSkewY];
771        affine[kASkewX] = this->fMat[kMSkewX];
772        affine[kAScaleY] = this->fMat[kMScaleY];
773        affine[kATransX] = this->fMat[kMTransX];
774        affine[kATransY] = this->fMat[kMTransY];
775    }
776    return true;
777}
778
779bool SkMatrix::invertNonIdentity(SkMatrix* inv) const {
780    SkASSERT(!this->isIdentity());
781
782    TypeMask mask = this->getType();
783
784    if (0 == (mask & ~(kScale_Mask | kTranslate_Mask))) {
785        bool invertible = true;
786        if (inv) {
787            if (mask & kScale_Mask) {
788                SkScalar invX = fMat[kMScaleX];
789                SkScalar invY = fMat[kMScaleY];
790                if (0 == invX || 0 == invY) {
791                    return false;
792                }
793                invX = SkScalarInvert(invX);
794                invY = SkScalarInvert(invY);
795
796                // Must be careful when writing to inv, since it may be the
797                // same memory as this.
798
799                inv->fMat[kMSkewX] = inv->fMat[kMSkewY] =
800                inv->fMat[kMPersp0] = inv->fMat[kMPersp1] = 0;
801
802                inv->fMat[kMScaleX] = invX;
803                inv->fMat[kMScaleY] = invY;
804                inv->fMat[kMPersp2] = 1;
805                inv->fMat[kMTransX] = -fMat[kMTransX] * invX;
806                inv->fMat[kMTransY] = -fMat[kMTransY] * invY;
807
808                inv->setTypeMask(mask | kRectStaysRect_Mask);
809            } else {
810                // translate only
811                inv->setTranslate(-fMat[kMTransX], -fMat[kMTransY]);
812            }
813        } else {    // inv is NULL, just check if we're invertible
814            if (!fMat[kMScaleX] || !fMat[kMScaleY]) {
815                invertible = false;
816            }
817        }
818        return invertible;
819    }
820
821    int    isPersp = mask & kPerspective_Mask;
822    double scale = sk_inv_determinant(fMat, isPersp);
823
824    if (scale == 0) { // underflow
825        return false;
826    }
827
828    if (inv) {
829        SkMatrix tmp;
830        if (inv == this) {
831            inv = &tmp;
832        }
833
834        if (isPersp) {
835            inv->fMat[kMScaleX] = scross_dscale(fMat[kMScaleY], fMat[kMPersp2], fMat[kMTransY], fMat[kMPersp1], scale);
836            inv->fMat[kMSkewX]  = scross_dscale(fMat[kMTransX], fMat[kMPersp1], fMat[kMSkewX],  fMat[kMPersp2], scale);
837            inv->fMat[kMTransX] = scross_dscale(fMat[kMSkewX],  fMat[kMTransY], fMat[kMTransX], fMat[kMScaleY], scale);
838
839            inv->fMat[kMSkewY]  = scross_dscale(fMat[kMTransY], fMat[kMPersp0], fMat[kMSkewY],  fMat[kMPersp2], scale);
840            inv->fMat[kMScaleY] = scross_dscale(fMat[kMScaleX], fMat[kMPersp2], fMat[kMTransX], fMat[kMPersp0], scale);
841            inv->fMat[kMTransY] = scross_dscale(fMat[kMTransX], fMat[kMSkewY],  fMat[kMScaleX], fMat[kMTransY], scale);
842
843            inv->fMat[kMPersp0] = scross_dscale(fMat[kMSkewY],  fMat[kMPersp1], fMat[kMScaleY], fMat[kMPersp0], scale);
844            inv->fMat[kMPersp1] = scross_dscale(fMat[kMSkewX],  fMat[kMPersp0], fMat[kMScaleX], fMat[kMPersp1], scale);
845            inv->fMat[kMPersp2] = scross_dscale(fMat[kMScaleX], fMat[kMScaleY], fMat[kMSkewX],  fMat[kMSkewY],  scale);
846        } else {   // not perspective
847            inv->fMat[kMScaleX] = SkDoubleToScalar(fMat[kMScaleY] * scale);
848            inv->fMat[kMSkewX]  = SkDoubleToScalar(-fMat[kMSkewX] * scale);
849            inv->fMat[kMTransX] = dcross_dscale(fMat[kMSkewX], fMat[kMTransY], fMat[kMScaleY], fMat[kMTransX], scale);
850
851            inv->fMat[kMSkewY]  = SkDoubleToScalar(-fMat[kMSkewY] * scale);
852            inv->fMat[kMScaleY] = SkDoubleToScalar(fMat[kMScaleX] * scale);
853            inv->fMat[kMTransY] = dcross_dscale(fMat[kMSkewY], fMat[kMTransX], fMat[kMScaleX], fMat[kMTransY], scale);
854
855            inv->fMat[kMPersp0] = 0;
856            inv->fMat[kMPersp1] = 0;
857            inv->fMat[kMPersp2] = 1;
858        }
859
860        inv->setTypeMask(fTypeMask);
861
862        if (inv == &tmp) {
863            *(SkMatrix*)this = tmp;
864        }
865    }
866    return true;
867}
868
869///////////////////////////////////////////////////////////////////////////////
870
871void SkMatrix::Identity_pts(const SkMatrix& m, SkPoint dst[], const SkPoint src[], int count) {
872    SkASSERT(m.getType() == 0);
873
874    if (dst != src && count > 0) {
875        memcpy(dst, src, count * sizeof(SkPoint));
876    }
877}
878
879void SkMatrix::Trans_pts(const SkMatrix& m, SkPoint dst[], const SkPoint src[], int count) {
880    SkASSERT(m.getType() <= kTranslate_Mask);
881
882    if (count > 0) {
883        SkScalar tx = m.getTranslateX();
884        SkScalar ty = m.getTranslateY();
885        if (count & 1) {
886            dst->fX = src->fX + tx;
887            dst->fY = src->fY + ty;
888            src += 1;
889            dst += 1;
890        }
891        Sk4s trans4(tx, ty, tx, ty);
892        count >>= 1;
893        if (count & 1) {
894            (Sk4s::Load(&src->fX) + trans4).store(&dst->fX);
895            src += 2;
896            dst += 2;
897        }
898        count >>= 1;
899        for (int i = 0; i < count; ++i) {
900            (Sk4s::Load(&src[0].fX) + trans4).store(&dst[0].fX);
901            (Sk4s::Load(&src[2].fX) + trans4).store(&dst[2].fX);
902            src += 4;
903            dst += 4;
904        }
905    }
906}
907
908void SkMatrix::Scale_pts(const SkMatrix& m, SkPoint dst[], const SkPoint src[], int count) {
909    SkASSERT(m.getType() <= (kScale_Mask | kTranslate_Mask));
910
911    if (count > 0) {
912        SkScalar tx = m.getTranslateX();
913        SkScalar ty = m.getTranslateY();
914        SkScalar sx = m.getScaleX();
915        SkScalar sy = m.getScaleY();
916        if (count & 1) {
917            dst->fX = src->fX * sx + tx;
918            dst->fY = src->fY * sy + ty;
919            src += 1;
920            dst += 1;
921        }
922        Sk4s trans4(tx, ty, tx, ty);
923        Sk4s scale4(sx, sy, sx, sy);
924        count >>= 1;
925        if (count & 1) {
926            (Sk4s::Load(&src->fX) * scale4 + trans4).store(&dst->fX);
927            src += 2;
928            dst += 2;
929        }
930        count >>= 1;
931        for (int i = 0; i < count; ++i) {
932            (Sk4s::Load(&src[0].fX) * scale4 + trans4).store(&dst[0].fX);
933            (Sk4s::Load(&src[2].fX) * scale4 + trans4).store(&dst[2].fX);
934            src += 4;
935            dst += 4;
936        }
937    }
938}
939
940void SkMatrix::Persp_pts(const SkMatrix& m, SkPoint dst[],
941                         const SkPoint src[], int count) {
942    SkASSERT(m.hasPerspective());
943
944    if (count > 0) {
945        do {
946            SkScalar sy = src->fY;
947            SkScalar sx = src->fX;
948            src += 1;
949
950            SkScalar x = sdot(sx, m.fMat[kMScaleX], sy, m.fMat[kMSkewX])  + m.fMat[kMTransX];
951            SkScalar y = sdot(sx, m.fMat[kMSkewY],  sy, m.fMat[kMScaleY]) + m.fMat[kMTransY];
952#ifdef SK_LEGACY_MATRIX_MATH_ORDER
953            SkScalar z = sx * m.fMat[kMPersp0] + (sy * m.fMat[kMPersp1] + m.fMat[kMPersp2]);
954#else
955            SkScalar z = sdot(sx, m.fMat[kMPersp0], sy, m.fMat[kMPersp1]) + m.fMat[kMPersp2];
956#endif
957            if (z) {
958                z = SkScalarFastInvert(z);
959            }
960
961            dst->fY = y * z;
962            dst->fX = x * z;
963            dst += 1;
964        } while (--count);
965    }
966}
967
968void SkMatrix::Affine_vpts(const SkMatrix& m, SkPoint dst[], const SkPoint src[], int count) {
969    SkASSERT(m.getType() != kPerspective_Mask);
970
971    if (count > 0) {
972        SkScalar tx = m.getTranslateX();
973        SkScalar ty = m.getTranslateY();
974        SkScalar sx = m.getScaleX();
975        SkScalar sy = m.getScaleY();
976        SkScalar kx = m.getSkewX();
977        SkScalar ky = m.getSkewY();
978        if (count & 1) {
979            dst->set(src->fX * sx + src->fY * kx + tx,
980                     src->fX * ky + src->fY * sy + ty);
981            src += 1;
982            dst += 1;
983        }
984        Sk4s trans4(tx, ty, tx, ty);
985        Sk4s scale4(sx, sy, sx, sy);
986        Sk4s  skew4(kx, ky, kx, ky);    // applied to swizzle of src4
987        count >>= 1;
988        for (int i = 0; i < count; ++i) {
989            Sk4s src4 = Sk4s::Load(&src->fX);
990            Sk4s swz4(src[0].fY, src[0].fX, src[1].fY, src[1].fX);  // need ABCD -> BADC
991            (src4 * scale4 + swz4 * skew4 + trans4).store(&dst->fX);
992            src += 2;
993            dst += 2;
994        }
995    }
996}
997
998const SkMatrix::MapPtsProc SkMatrix::gMapPtsProcs[] = {
999    SkMatrix::Identity_pts, SkMatrix::Trans_pts,
1000    SkMatrix::Scale_pts,   SkMatrix::Scale_pts,
1001    SkMatrix::Affine_vpts,  SkMatrix::Affine_vpts,
1002    SkMatrix::Affine_vpts,  SkMatrix::Affine_vpts,
1003    // repeat the persp proc 8 times
1004    SkMatrix::Persp_pts,    SkMatrix::Persp_pts,
1005    SkMatrix::Persp_pts,    SkMatrix::Persp_pts,
1006    SkMatrix::Persp_pts,    SkMatrix::Persp_pts,
1007    SkMatrix::Persp_pts,    SkMatrix::Persp_pts
1008};
1009
1010///////////////////////////////////////////////////////////////////////////////
1011
1012void SkMatrix::mapHomogeneousPoints(SkScalar dst[], const SkScalar src[], int count) const {
1013    SkASSERT((dst && src && count > 0) || 0 == count);
1014    // no partial overlap
1015    SkASSERT(src == dst || &dst[3*count] <= &src[0] || &src[3*count] <= &dst[0]);
1016
1017    if (count > 0) {
1018        if (this->isIdentity()) {
1019            memcpy(dst, src, 3*count*sizeof(SkScalar));
1020            return;
1021        }
1022        do {
1023            SkScalar sx = src[0];
1024            SkScalar sy = src[1];
1025            SkScalar sw = src[2];
1026            src += 3;
1027
1028            SkScalar x = sdot(sx, fMat[kMScaleX], sy, fMat[kMSkewX],  sw, fMat[kMTransX]);
1029            SkScalar y = sdot(sx, fMat[kMSkewY],  sy, fMat[kMScaleY], sw, fMat[kMTransY]);
1030            SkScalar w = sdot(sx, fMat[kMPersp0], sy, fMat[kMPersp1], sw, fMat[kMPersp2]);
1031
1032            dst[0] = x;
1033            dst[1] = y;
1034            dst[2] = w;
1035            dst += 3;
1036        } while (--count);
1037    }
1038}
1039
1040///////////////////////////////////////////////////////////////////////////////
1041
1042void SkMatrix::mapVectors(SkPoint dst[], const SkPoint src[], int count) const {
1043    if (this->hasPerspective()) {
1044        SkPoint origin;
1045
1046        MapXYProc proc = this->getMapXYProc();
1047        proc(*this, 0, 0, &origin);
1048
1049        for (int i = count - 1; i >= 0; --i) {
1050            SkPoint tmp;
1051
1052            proc(*this, src[i].fX, src[i].fY, &tmp);
1053            dst[i].set(tmp.fX - origin.fX, tmp.fY - origin.fY);
1054        }
1055    } else {
1056        SkMatrix tmp = *this;
1057
1058        tmp.fMat[kMTransX] = tmp.fMat[kMTransY] = 0;
1059        tmp.clearTypeMask(kTranslate_Mask);
1060        tmp.mapPoints(dst, src, count);
1061    }
1062}
1063
1064bool SkMatrix::mapRect(SkRect* dst, const SkRect& src) const {
1065    SkASSERT(dst);
1066
1067    if (this->rectStaysRect()) {
1068        this->mapPoints((SkPoint*)dst, (const SkPoint*)&src, 2);
1069        dst->sort();
1070        return true;
1071    } else {
1072        SkPoint quad[4];
1073
1074        src.toQuad(quad);
1075        this->mapPoints(quad, quad, 4);
1076        dst->set(quad, 4);
1077        return false;
1078    }
1079}
1080
1081SkScalar SkMatrix::mapRadius(SkScalar radius) const {
1082    SkVector    vec[2];
1083
1084    vec[0].set(radius, 0);
1085    vec[1].set(0, radius);
1086    this->mapVectors(vec, 2);
1087
1088    SkScalar d0 = vec[0].length();
1089    SkScalar d1 = vec[1].length();
1090
1091    // return geometric mean
1092    return SkScalarSqrt(d0 * d1);
1093}
1094
1095///////////////////////////////////////////////////////////////////////////////
1096
1097void SkMatrix::Persp_xy(const SkMatrix& m, SkScalar sx, SkScalar sy,
1098                        SkPoint* pt) {
1099    SkASSERT(m.hasPerspective());
1100
1101    SkScalar x = sdot(sx, m.fMat[kMScaleX], sy, m.fMat[kMSkewX])  + m.fMat[kMTransX];
1102    SkScalar y = sdot(sx, m.fMat[kMSkewY],  sy, m.fMat[kMScaleY]) + m.fMat[kMTransY];
1103    SkScalar z = sdot(sx, m.fMat[kMPersp0], sy, m.fMat[kMPersp1]) + m.fMat[kMPersp2];
1104    if (z) {
1105        z = SkScalarFastInvert(z);
1106    }
1107    pt->fX = x * z;
1108    pt->fY = y * z;
1109}
1110
1111void SkMatrix::RotTrans_xy(const SkMatrix& m, SkScalar sx, SkScalar sy,
1112                           SkPoint* pt) {
1113    SkASSERT((m.getType() & (kAffine_Mask | kPerspective_Mask)) == kAffine_Mask);
1114
1115#ifdef SK_LEGACY_MATRIX_MATH_ORDER
1116    pt->fX = sx * m.fMat[kMScaleX] + (sy * m.fMat[kMSkewX]  +  m.fMat[kMTransX]);
1117    pt->fY = sx * m.fMat[kMSkewY]  + (sy * m.fMat[kMScaleY] + m.fMat[kMTransY]);
1118#else
1119    pt->fX = sdot(sx, m.fMat[kMScaleX], sy, m.fMat[kMSkewX])  + m.fMat[kMTransX];
1120    pt->fY = sdot(sx, m.fMat[kMSkewY],  sy, m.fMat[kMScaleY]) + m.fMat[kMTransY];
1121#endif
1122}
1123
1124void SkMatrix::Rot_xy(const SkMatrix& m, SkScalar sx, SkScalar sy,
1125                      SkPoint* pt) {
1126    SkASSERT((m.getType() & (kAffine_Mask | kPerspective_Mask))== kAffine_Mask);
1127    SkASSERT(0 == m.fMat[kMTransX]);
1128    SkASSERT(0 == m.fMat[kMTransY]);
1129
1130#ifdef SK_LEGACY_MATRIX_MATH_ORDER
1131    pt->fX = sx * m.fMat[kMScaleX] + (sy * m.fMat[kMSkewX]  + m.fMat[kMTransX]);
1132    pt->fY = sx * m.fMat[kMSkewY]  + (sy * m.fMat[kMScaleY] + m.fMat[kMTransY]);
1133#else
1134    pt->fX = sdot(sx, m.fMat[kMScaleX], sy, m.fMat[kMSkewX])  + m.fMat[kMTransX];
1135    pt->fY = sdot(sx, m.fMat[kMSkewY],  sy, m.fMat[kMScaleY]) + m.fMat[kMTransY];
1136#endif
1137}
1138
1139void SkMatrix::ScaleTrans_xy(const SkMatrix& m, SkScalar sx, SkScalar sy,
1140                             SkPoint* pt) {
1141    SkASSERT((m.getType() & (kScale_Mask | kAffine_Mask | kPerspective_Mask))
1142             == kScale_Mask);
1143
1144    pt->fX = sx * m.fMat[kMScaleX] + m.fMat[kMTransX];
1145    pt->fY = sy * m.fMat[kMScaleY] + m.fMat[kMTransY];
1146}
1147
1148void SkMatrix::Scale_xy(const SkMatrix& m, SkScalar sx, SkScalar sy,
1149                        SkPoint* pt) {
1150    SkASSERT((m.getType() & (kScale_Mask | kAffine_Mask | kPerspective_Mask))
1151             == kScale_Mask);
1152    SkASSERT(0 == m.fMat[kMTransX]);
1153    SkASSERT(0 == m.fMat[kMTransY]);
1154
1155    pt->fX = sx * m.fMat[kMScaleX];
1156    pt->fY = sy * m.fMat[kMScaleY];
1157}
1158
1159void SkMatrix::Trans_xy(const SkMatrix& m, SkScalar sx, SkScalar sy,
1160                        SkPoint* pt) {
1161    SkASSERT(m.getType() == kTranslate_Mask);
1162
1163    pt->fX = sx + m.fMat[kMTransX];
1164    pt->fY = sy + m.fMat[kMTransY];
1165}
1166
1167void SkMatrix::Identity_xy(const SkMatrix& m, SkScalar sx, SkScalar sy,
1168                           SkPoint* pt) {
1169    SkASSERT(0 == m.getType());
1170
1171    pt->fX = sx;
1172    pt->fY = sy;
1173}
1174
1175const SkMatrix::MapXYProc SkMatrix::gMapXYProcs[] = {
1176    SkMatrix::Identity_xy, SkMatrix::Trans_xy,
1177    SkMatrix::Scale_xy,    SkMatrix::ScaleTrans_xy,
1178    SkMatrix::Rot_xy,      SkMatrix::RotTrans_xy,
1179    SkMatrix::Rot_xy,      SkMatrix::RotTrans_xy,
1180    // repeat the persp proc 8 times
1181    SkMatrix::Persp_xy,    SkMatrix::Persp_xy,
1182    SkMatrix::Persp_xy,    SkMatrix::Persp_xy,
1183    SkMatrix::Persp_xy,    SkMatrix::Persp_xy,
1184    SkMatrix::Persp_xy,    SkMatrix::Persp_xy
1185};
1186
1187///////////////////////////////////////////////////////////////////////////////
1188
1189// if its nearly zero (just made up 26, perhaps it should be bigger or smaller)
1190#define PerspNearlyZero(x)  SkScalarNearlyZero(x, (1.0f / (1 << 26)))
1191
1192bool SkMatrix::fixedStepInX(SkScalar y, SkFixed* stepX, SkFixed* stepY) const {
1193    if (PerspNearlyZero(fMat[kMPersp0])) {
1194        if (stepX || stepY) {
1195            if (PerspNearlyZero(fMat[kMPersp1]) &&
1196                    PerspNearlyZero(fMat[kMPersp2] - 1)) {
1197                if (stepX) {
1198                    *stepX = SkScalarToFixed(fMat[kMScaleX]);
1199                }
1200                if (stepY) {
1201                    *stepY = SkScalarToFixed(fMat[kMSkewY]);
1202                }
1203            } else {
1204                SkScalar z = y * fMat[kMPersp1] + fMat[kMPersp2];
1205                if (stepX) {
1206                    *stepX = SkScalarToFixed(fMat[kMScaleX] / z);
1207                }
1208                if (stepY) {
1209                    *stepY = SkScalarToFixed(fMat[kMSkewY] / z);
1210                }
1211            }
1212        }
1213        return true;
1214    }
1215    return false;
1216}
1217
1218///////////////////////////////////////////////////////////////////////////////
1219
1220#include "SkPerspIter.h"
1221
1222SkPerspIter::SkPerspIter(const SkMatrix& m, SkScalar x0, SkScalar y0, int count)
1223        : fMatrix(m), fSX(x0), fSY(y0), fCount(count) {
1224    SkPoint pt;
1225
1226    SkMatrix::Persp_xy(m, x0, y0, &pt);
1227    fX = SkScalarToFixed(pt.fX);
1228    fY = SkScalarToFixed(pt.fY);
1229}
1230
1231int SkPerspIter::next() {
1232    int n = fCount;
1233
1234    if (0 == n) {
1235        return 0;
1236    }
1237    SkPoint pt;
1238    SkFixed x = fX;
1239    SkFixed y = fY;
1240    SkFixed dx, dy;
1241
1242    if (n >= kCount) {
1243        n = kCount;
1244        fSX += SkIntToScalar(kCount);
1245        SkMatrix::Persp_xy(fMatrix, fSX, fSY, &pt);
1246        fX = SkScalarToFixed(pt.fX);
1247        fY = SkScalarToFixed(pt.fY);
1248        dx = (fX - x) >> kShift;
1249        dy = (fY - y) >> kShift;
1250    } else {
1251        fSX += SkIntToScalar(n);
1252        SkMatrix::Persp_xy(fMatrix, fSX, fSY, &pt);
1253        fX = SkScalarToFixed(pt.fX);
1254        fY = SkScalarToFixed(pt.fY);
1255        dx = (fX - x) / n;
1256        dy = (fY - y) / n;
1257    }
1258
1259    SkFixed* p = fStorage;
1260    for (int i = 0; i < n; i++) {
1261        *p++ = x; x += dx;
1262        *p++ = y; y += dy;
1263    }
1264
1265    fCount -= n;
1266    return n;
1267}
1268
1269///////////////////////////////////////////////////////////////////////////////
1270
1271static inline bool checkForZero(float x) {
1272    return x*x == 0;
1273}
1274
1275static inline bool poly_to_point(SkPoint* pt, const SkPoint poly[], int count) {
1276    float   x = 1, y = 1;
1277    SkPoint pt1, pt2;
1278
1279    if (count > 1) {
1280        pt1.fX = poly[1].fX - poly[0].fX;
1281        pt1.fY = poly[1].fY - poly[0].fY;
1282        y = SkPoint::Length(pt1.fX, pt1.fY);
1283        if (checkForZero(y)) {
1284            return false;
1285        }
1286        switch (count) {
1287            case 2:
1288                break;
1289            case 3:
1290                pt2.fX = poly[0].fY - poly[2].fY;
1291                pt2.fY = poly[2].fX - poly[0].fX;
1292                goto CALC_X;
1293            default:
1294                pt2.fX = poly[0].fY - poly[3].fY;
1295                pt2.fY = poly[3].fX - poly[0].fX;
1296            CALC_X:
1297                x = sdot(pt1.fX, pt2.fX, pt1.fY, pt2.fY) / y;
1298                break;
1299        }
1300    }
1301    pt->set(x, y);
1302    return true;
1303}
1304
1305bool SkMatrix::Poly2Proc(const SkPoint srcPt[], SkMatrix* dst,
1306                         const SkPoint& scale) {
1307    float invScale = 1 / scale.fY;
1308
1309    dst->fMat[kMScaleX] = (srcPt[1].fY - srcPt[0].fY) * invScale;
1310    dst->fMat[kMSkewY] = (srcPt[0].fX - srcPt[1].fX) * invScale;
1311    dst->fMat[kMPersp0] = 0;
1312    dst->fMat[kMSkewX] = (srcPt[1].fX - srcPt[0].fX) * invScale;
1313    dst->fMat[kMScaleY] = (srcPt[1].fY - srcPt[0].fY) * invScale;
1314    dst->fMat[kMPersp1] = 0;
1315    dst->fMat[kMTransX] = srcPt[0].fX;
1316    dst->fMat[kMTransY] = srcPt[0].fY;
1317    dst->fMat[kMPersp2] = 1;
1318    dst->setTypeMask(kUnknown_Mask);
1319    return true;
1320}
1321
1322bool SkMatrix::Poly3Proc(const SkPoint srcPt[], SkMatrix* dst,
1323                         const SkPoint& scale) {
1324    float invScale = 1 / scale.fX;
1325    dst->fMat[kMScaleX] = (srcPt[2].fX - srcPt[0].fX) * invScale;
1326    dst->fMat[kMSkewY] = (srcPt[2].fY - srcPt[0].fY) * invScale;
1327    dst->fMat[kMPersp0] = 0;
1328
1329    invScale = 1 / scale.fY;
1330    dst->fMat[kMSkewX] = (srcPt[1].fX - srcPt[0].fX) * invScale;
1331    dst->fMat[kMScaleY] = (srcPt[1].fY - srcPt[0].fY) * invScale;
1332    dst->fMat[kMPersp1] = 0;
1333
1334    dst->fMat[kMTransX] = srcPt[0].fX;
1335    dst->fMat[kMTransY] = srcPt[0].fY;
1336    dst->fMat[kMPersp2] = 1;
1337    dst->setTypeMask(kUnknown_Mask);
1338    return true;
1339}
1340
1341bool SkMatrix::Poly4Proc(const SkPoint srcPt[], SkMatrix* dst,
1342                         const SkPoint& scale) {
1343    float   a1, a2;
1344    float   x0, y0, x1, y1, x2, y2;
1345
1346    x0 = srcPt[2].fX - srcPt[0].fX;
1347    y0 = srcPt[2].fY - srcPt[0].fY;
1348    x1 = srcPt[2].fX - srcPt[1].fX;
1349    y1 = srcPt[2].fY - srcPt[1].fY;
1350    x2 = srcPt[2].fX - srcPt[3].fX;
1351    y2 = srcPt[2].fY - srcPt[3].fY;
1352
1353    /* check if abs(x2) > abs(y2) */
1354    if ( x2 > 0 ? y2 > 0 ? x2 > y2 : x2 > -y2 : y2 > 0 ? -x2 > y2 : x2 < y2) {
1355        float denom = SkScalarMulDiv(x1, y2, x2) - y1;
1356        if (checkForZero(denom)) {
1357            return false;
1358        }
1359        a1 = (SkScalarMulDiv(x0 - x1, y2, x2) - y0 + y1) / denom;
1360    } else {
1361        float denom = x1 - SkScalarMulDiv(y1, x2, y2);
1362        if (checkForZero(denom)) {
1363            return false;
1364        }
1365        a1 = (x0 - x1 - SkScalarMulDiv(y0 - y1, x2, y2)) / denom;
1366    }
1367
1368    /* check if abs(x1) > abs(y1) */
1369    if ( x1 > 0 ? y1 > 0 ? x1 > y1 : x1 > -y1 : y1 > 0 ? -x1 > y1 : x1 < y1) {
1370        float denom = y2 - SkScalarMulDiv(x2, y1, x1);
1371        if (checkForZero(denom)) {
1372            return false;
1373        }
1374        a2 = (y0 - y2 - SkScalarMulDiv(x0 - x2, y1, x1)) / denom;
1375    } else {
1376        float denom = SkScalarMulDiv(y2, x1, y1) - x2;
1377        if (checkForZero(denom)) {
1378            return false;
1379        }
1380        a2 = (SkScalarMulDiv(y0 - y2, x1, y1) - x0 + x2) / denom;
1381    }
1382
1383    float invScale = SkScalarInvert(scale.fX);
1384    dst->fMat[kMScaleX] = (a2 * srcPt[3].fX + srcPt[3].fX - srcPt[0].fX) * invScale;
1385    dst->fMat[kMSkewY]  = (a2 * srcPt[3].fY + srcPt[3].fY - srcPt[0].fY) * invScale;
1386    dst->fMat[kMPersp0] = a2 * invScale;
1387
1388    invScale = SkScalarInvert(scale.fY);
1389    dst->fMat[kMSkewX]  = (a1 * srcPt[1].fX + srcPt[1].fX - srcPt[0].fX) * invScale;
1390    dst->fMat[kMScaleY] = (a1 * srcPt[1].fY + srcPt[1].fY - srcPt[0].fY) * invScale;
1391    dst->fMat[kMPersp1] = a1 * invScale;
1392
1393    dst->fMat[kMTransX] = srcPt[0].fX;
1394    dst->fMat[kMTransY] = srcPt[0].fY;
1395    dst->fMat[kMPersp2] = 1;
1396    dst->setTypeMask(kUnknown_Mask);
1397    return true;
1398}
1399
1400typedef bool (*PolyMapProc)(const SkPoint[], SkMatrix*, const SkPoint&);
1401
1402/*  Taken from Rob Johnson's original sample code in QuickDraw GX
1403*/
1404bool SkMatrix::setPolyToPoly(const SkPoint src[], const SkPoint dst[],
1405                             int count) {
1406    if ((unsigned)count > 4) {
1407        SkDebugf("--- SkMatrix::setPolyToPoly count out of range %d\n", count);
1408        return false;
1409    }
1410
1411    if (0 == count) {
1412        this->reset();
1413        return true;
1414    }
1415    if (1 == count) {
1416        this->setTranslate(dst[0].fX - src[0].fX, dst[0].fY - src[0].fY);
1417        return true;
1418    }
1419
1420    SkPoint scale;
1421    if (!poly_to_point(&scale, src, count) ||
1422            SkScalarNearlyZero(scale.fX) ||
1423            SkScalarNearlyZero(scale.fY)) {
1424        return false;
1425    }
1426
1427    static const PolyMapProc gPolyMapProcs[] = {
1428        SkMatrix::Poly2Proc, SkMatrix::Poly3Proc, SkMatrix::Poly4Proc
1429    };
1430    PolyMapProc proc = gPolyMapProcs[count - 2];
1431
1432    SkMatrix tempMap, result;
1433    tempMap.setTypeMask(kUnknown_Mask);
1434
1435    if (!proc(src, &tempMap, scale)) {
1436        return false;
1437    }
1438    if (!tempMap.invert(&result)) {
1439        return false;
1440    }
1441    if (!proc(dst, &tempMap, scale)) {
1442        return false;
1443    }
1444    this->setConcat(tempMap, result);
1445    return true;
1446}
1447
1448///////////////////////////////////////////////////////////////////////////////
1449
1450enum MinMaxOrBoth {
1451    kMin_MinMaxOrBoth,
1452    kMax_MinMaxOrBoth,
1453    kBoth_MinMaxOrBoth
1454};
1455
1456template <MinMaxOrBoth MIN_MAX_OR_BOTH> bool get_scale_factor(SkMatrix::TypeMask typeMask,
1457                                                              const SkScalar m[9],
1458                                                              SkScalar results[/*1 or 2*/]) {
1459    if (typeMask & SkMatrix::kPerspective_Mask) {
1460        return false;
1461    }
1462    if (SkMatrix::kIdentity_Mask == typeMask) {
1463        results[0] = SK_Scalar1;
1464        if (kBoth_MinMaxOrBoth == MIN_MAX_OR_BOTH) {
1465            results[1] = SK_Scalar1;
1466        }
1467        return true;
1468    }
1469    if (!(typeMask & SkMatrix::kAffine_Mask)) {
1470        if (kMin_MinMaxOrBoth == MIN_MAX_OR_BOTH) {
1471             results[0] = SkMinScalar(SkScalarAbs(m[SkMatrix::kMScaleX]),
1472                                      SkScalarAbs(m[SkMatrix::kMScaleY]));
1473        } else if (kMax_MinMaxOrBoth == MIN_MAX_OR_BOTH) {
1474             results[0] = SkMaxScalar(SkScalarAbs(m[SkMatrix::kMScaleX]),
1475                                      SkScalarAbs(m[SkMatrix::kMScaleY]));
1476        } else {
1477            results[0] = SkScalarAbs(m[SkMatrix::kMScaleX]);
1478            results[1] = SkScalarAbs(m[SkMatrix::kMScaleY]);
1479             if (results[0] > results[1]) {
1480                 SkTSwap(results[0], results[1]);
1481             }
1482        }
1483        return true;
1484    }
1485    // ignore the translation part of the matrix, just look at 2x2 portion.
1486    // compute singular values, take largest or smallest abs value.
1487    // [a b; b c] = A^T*A
1488    SkScalar a = sdot(m[SkMatrix::kMScaleX], m[SkMatrix::kMScaleX],
1489                      m[SkMatrix::kMSkewY],  m[SkMatrix::kMSkewY]);
1490    SkScalar b = sdot(m[SkMatrix::kMScaleX], m[SkMatrix::kMSkewX],
1491                      m[SkMatrix::kMScaleY], m[SkMatrix::kMSkewY]);
1492    SkScalar c = sdot(m[SkMatrix::kMSkewX],  m[SkMatrix::kMSkewX],
1493                      m[SkMatrix::kMScaleY], m[SkMatrix::kMScaleY]);
1494    // eigenvalues of A^T*A are the squared singular values of A.
1495    // characteristic equation is det((A^T*A) - l*I) = 0
1496    // l^2 - (a + c)l + (ac-b^2)
1497    // solve using quadratic equation (divisor is non-zero since l^2 has 1 coeff
1498    // and roots are guaranteed to be pos and real).
1499    SkScalar bSqd = b * b;
1500    // if upper left 2x2 is orthogonal save some math
1501    if (bSqd <= SK_ScalarNearlyZero*SK_ScalarNearlyZero) {
1502        if (kMin_MinMaxOrBoth == MIN_MAX_OR_BOTH) {
1503            results[0] = SkMinScalar(a, c);
1504        } else if (kMax_MinMaxOrBoth == MIN_MAX_OR_BOTH) {
1505            results[0] = SkMaxScalar(a, c);
1506        } else {
1507            results[0] = a;
1508            results[1] = c;
1509            if (results[0] > results[1]) {
1510                SkTSwap(results[0], results[1]);
1511            }
1512        }
1513    } else {
1514        SkScalar aminusc = a - c;
1515        SkScalar apluscdiv2 = SkScalarHalf(a + c);
1516        SkScalar x = SkScalarHalf(SkScalarSqrt(aminusc * aminusc + 4 * bSqd));
1517        if (kMin_MinMaxOrBoth == MIN_MAX_OR_BOTH) {
1518            results[0] = apluscdiv2 - x;
1519        } else if (kMax_MinMaxOrBoth == MIN_MAX_OR_BOTH) {
1520            results[0] = apluscdiv2 + x;
1521        } else {
1522            results[0] = apluscdiv2 - x;
1523            results[1] = apluscdiv2 + x;
1524        }
1525    }
1526    SkASSERT(results[0] >= 0);
1527    results[0] = SkScalarSqrt(results[0]);
1528    if (kBoth_MinMaxOrBoth == MIN_MAX_OR_BOTH) {
1529        SkASSERT(results[1] >= 0);
1530        results[1] = SkScalarSqrt(results[1]);
1531    }
1532    return true;
1533}
1534
1535SkScalar SkMatrix::getMinScale() const {
1536    SkScalar factor;
1537    if (get_scale_factor<kMin_MinMaxOrBoth>(this->getType(), fMat, &factor)) {
1538        return factor;
1539    } else {
1540        return -1;
1541    }
1542}
1543
1544SkScalar SkMatrix::getMaxScale() const {
1545    SkScalar factor;
1546    if (get_scale_factor<kMax_MinMaxOrBoth>(this->getType(), fMat, &factor)) {
1547        return factor;
1548    } else {
1549        return -1;
1550    }
1551}
1552
1553bool SkMatrix::getMinMaxScales(SkScalar scaleFactors[2]) const {
1554    return get_scale_factor<kBoth_MinMaxOrBoth>(this->getType(), fMat, scaleFactors);
1555}
1556
1557namespace {
1558
1559struct PODMatrix {
1560    SkScalar matrix[9];
1561    uint32_t typemask;
1562
1563    const SkMatrix& asSkMatrix() const { return *reinterpret_cast<const SkMatrix*>(this); }
1564};
1565SK_COMPILE_ASSERT(sizeof(PODMatrix) == sizeof(SkMatrix), PODMatrixSizeMismatch);
1566
1567}  // namespace
1568
1569const SkMatrix& SkMatrix::I() {
1570    SK_COMPILE_ASSERT(offsetof(SkMatrix, fMat)      == offsetof(PODMatrix, matrix),   BadfMat);
1571    SK_COMPILE_ASSERT(offsetof(SkMatrix, fTypeMask) == offsetof(PODMatrix, typemask), BadfTypeMask);
1572
1573    static const PODMatrix identity = { {SK_Scalar1, 0, 0,
1574                                         0, SK_Scalar1, 0,
1575                                         0, 0, SK_Scalar1 },
1576                                       kIdentity_Mask | kRectStaysRect_Mask};
1577    SkASSERT(identity.asSkMatrix().isIdentity());
1578    return identity.asSkMatrix();
1579}
1580
1581const SkMatrix& SkMatrix::InvalidMatrix() {
1582    SK_COMPILE_ASSERT(offsetof(SkMatrix, fMat)      == offsetof(PODMatrix, matrix),   BadfMat);
1583    SK_COMPILE_ASSERT(offsetof(SkMatrix, fTypeMask) == offsetof(PODMatrix, typemask), BadfTypeMask);
1584
1585    static const PODMatrix invalid =
1586        { {SK_ScalarMax, SK_ScalarMax, SK_ScalarMax,
1587           SK_ScalarMax, SK_ScalarMax, SK_ScalarMax,
1588           SK_ScalarMax, SK_ScalarMax, SK_ScalarMax },
1589         kTranslate_Mask | kScale_Mask | kAffine_Mask | kPerspective_Mask };
1590    return invalid.asSkMatrix();
1591}
1592
1593bool SkMatrix::decomposeScale(SkSize* scale, SkMatrix* remaining) const {
1594    if (this->hasPerspective()) {
1595        return false;
1596    }
1597
1598    const SkScalar sx = SkVector::Length(this->getScaleX(), this->getSkewY());
1599    const SkScalar sy = SkVector::Length(this->getSkewX(), this->getScaleY());
1600    if (!SkScalarIsFinite(sx) || !SkScalarIsFinite(sy) ||
1601        SkScalarNearlyZero(sx) || SkScalarNearlyZero(sy)) {
1602        return false;
1603    }
1604
1605    if (scale) {
1606        scale->set(sx, sy);
1607    }
1608    if (remaining) {
1609        *remaining = *this;
1610        remaining->postScale(SkScalarInvert(sx), SkScalarInvert(sy));
1611    }
1612    return true;
1613}
1614
1615///////////////////////////////////////////////////////////////////////////////
1616
1617size_t SkMatrix::writeToMemory(void* buffer) const {
1618    // TODO write less for simple matrices
1619    static const size_t sizeInMemory = 9 * sizeof(SkScalar);
1620    if (buffer) {
1621        memcpy(buffer, fMat, sizeInMemory);
1622    }
1623    return sizeInMemory;
1624}
1625
1626size_t SkMatrix::readFromMemory(const void* buffer, size_t length) {
1627    static const size_t sizeInMemory = 9 * sizeof(SkScalar);
1628    if (length < sizeInMemory) {
1629        return 0;
1630    }
1631    if (buffer) {
1632        memcpy(fMat, buffer, sizeInMemory);
1633        this->setTypeMask(kUnknown_Mask);
1634    }
1635    return sizeInMemory;
1636}
1637
1638void SkMatrix::dump() const {
1639    SkString str;
1640    this->toString(&str);
1641    SkDebugf("%s\n", str.c_str());
1642}
1643
1644void SkMatrix::toString(SkString* str) const {
1645    str->appendf("[%8.4f %8.4f %8.4f][%8.4f %8.4f %8.4f][%8.4f %8.4f %8.4f]",
1646             fMat[0], fMat[1], fMat[2], fMat[3], fMat[4], fMat[5],
1647             fMat[6], fMat[7], fMat[8]);
1648}
1649
1650///////////////////////////////////////////////////////////////////////////////
1651
1652#include "SkMatrixUtils.h"
1653
1654bool SkTreatAsSprite(const SkMatrix& mat, int width, int height,
1655                     unsigned subpixelBits) {
1656    // quick reject on affine or perspective
1657    if (mat.getType() & ~(SkMatrix::kScale_Mask | SkMatrix::kTranslate_Mask)) {
1658        return false;
1659    }
1660
1661    // quick success check
1662    if (!subpixelBits && !(mat.getType() & ~SkMatrix::kTranslate_Mask)) {
1663        return true;
1664    }
1665
1666    // mapRect supports negative scales, so we eliminate those first
1667    if (mat.getScaleX() < 0 || mat.getScaleY() < 0) {
1668        return false;
1669    }
1670
1671    SkRect dst;
1672    SkIRect isrc = { 0, 0, width, height };
1673
1674    {
1675        SkRect src;
1676        src.set(isrc);
1677        mat.mapRect(&dst, src);
1678    }
1679
1680    // just apply the translate to isrc
1681    isrc.offset(SkScalarRoundToInt(mat.getTranslateX()),
1682                SkScalarRoundToInt(mat.getTranslateY()));
1683
1684    if (subpixelBits) {
1685        isrc.fLeft <<= subpixelBits;
1686        isrc.fTop <<= subpixelBits;
1687        isrc.fRight <<= subpixelBits;
1688        isrc.fBottom <<= subpixelBits;
1689
1690        const float scale = 1 << subpixelBits;
1691        dst.fLeft *= scale;
1692        dst.fTop *= scale;
1693        dst.fRight *= scale;
1694        dst.fBottom *= scale;
1695    }
1696
1697    SkIRect idst;
1698    dst.round(&idst);
1699    return isrc == idst;
1700}
1701
1702// A square matrix M can be decomposed (via polar decomposition) into two matrices --
1703// an orthogonal matrix Q and a symmetric matrix S. In turn we can decompose S into U*W*U^T,
1704// where U is another orthogonal matrix and W is a scale matrix. These can be recombined
1705// to give M = (Q*U)*W*U^T, i.e., the product of two orthogonal matrices and a scale matrix.
1706//
1707// The one wrinkle is that traditionally Q may contain a reflection -- the
1708// calculation has been rejiggered to put that reflection into W.
1709bool SkDecomposeUpper2x2(const SkMatrix& matrix,
1710                         SkPoint* rotation1,
1711                         SkPoint* scale,
1712                         SkPoint* rotation2) {
1713
1714    SkScalar A = matrix[SkMatrix::kMScaleX];
1715    SkScalar B = matrix[SkMatrix::kMSkewX];
1716    SkScalar C = matrix[SkMatrix::kMSkewY];
1717    SkScalar D = matrix[SkMatrix::kMScaleY];
1718
1719    if (is_degenerate_2x2(A, B, C, D)) {
1720        return false;
1721    }
1722
1723    double w1, w2;
1724    SkScalar cos1, sin1;
1725    SkScalar cos2, sin2;
1726
1727    // do polar decomposition (M = Q*S)
1728    SkScalar cosQ, sinQ;
1729    double Sa, Sb, Sd;
1730    // if M is already symmetric (i.e., M = I*S)
1731    if (SkScalarNearlyEqual(B, C)) {
1732        cosQ = 1;
1733        sinQ = 0;
1734
1735        Sa = A;
1736        Sb = B;
1737        Sd = D;
1738    } else {
1739        cosQ = A + D;
1740        sinQ = C - B;
1741        SkScalar reciplen = SkScalarInvert(SkScalarSqrt(cosQ*cosQ + sinQ*sinQ));
1742        cosQ *= reciplen;
1743        sinQ *= reciplen;
1744
1745        // S = Q^-1*M
1746        // we don't calc Sc since it's symmetric
1747        Sa = A*cosQ + C*sinQ;
1748        Sb = B*cosQ + D*sinQ;
1749        Sd = -B*sinQ + D*cosQ;
1750    }
1751
1752    // Now we need to compute eigenvalues of S (our scale factors)
1753    // and eigenvectors (bases for our rotation)
1754    // From this, should be able to reconstruct S as U*W*U^T
1755    if (SkScalarNearlyZero(SkDoubleToScalar(Sb))) {
1756        // already diagonalized
1757        cos1 = 1;
1758        sin1 = 0;
1759        w1 = Sa;
1760        w2 = Sd;
1761        cos2 = cosQ;
1762        sin2 = sinQ;
1763    } else {
1764        double diff = Sa - Sd;
1765        double discriminant = sqrt(diff*diff + 4.0*Sb*Sb);
1766        double trace = Sa + Sd;
1767        if (diff > 0) {
1768            w1 = 0.5*(trace + discriminant);
1769            w2 = 0.5*(trace - discriminant);
1770        } else {
1771            w1 = 0.5*(trace - discriminant);
1772            w2 = 0.5*(trace + discriminant);
1773        }
1774
1775        cos1 = SkDoubleToScalar(Sb); sin1 = SkDoubleToScalar(w1 - Sa);
1776        SkScalar reciplen = SkScalarInvert(SkScalarSqrt(cos1*cos1 + sin1*sin1));
1777        cos1 *= reciplen;
1778        sin1 *= reciplen;
1779
1780        // rotation 2 is composition of Q and U
1781        cos2 = cos1*cosQ - sin1*sinQ;
1782        sin2 = sin1*cosQ + cos1*sinQ;
1783
1784        // rotation 1 is U^T
1785        sin1 = -sin1;
1786    }
1787
1788    if (scale) {
1789        scale->fX = SkDoubleToScalar(w1);
1790        scale->fY = SkDoubleToScalar(w2);
1791    }
1792    if (rotation1) {
1793        rotation1->fX = cos1;
1794        rotation1->fY = sin1;
1795    }
1796    if (rotation2) {
1797        rotation2->fX = cos2;
1798        rotation2->fY = sin2;
1799    }
1800
1801    return true;
1802}
1803