path: root/gfx/2d/Matrix.h

/* -*- Mode: C++; tab-width: 20; indent-tabs-mode: nil; c-basic-offset: 2 -*-
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */

#ifndef MOZILLA_GFX_MATRIX_H_
#define MOZILLA_GFX_MATRIX_H_

#include "Types.h"
#include "Triangle.h"
#include "Rect.h"
#include "Point.h"
#include "Quaternion.h"
#include <iosfwd>
#include <math.h>
#include "mozilla/Attributes.h"
#include "mozilla/DebugOnly.h"
#include "mozilla/FloatingPoint.h"

namespace mozilla {
namespace gfx {

static bool FuzzyEqual(Float aV1, Float aV2) {
  // XXX - Check if fabs does the smart thing and just negates the sign bit.
  return fabs(aV2 - aV1) < 1e-6;
}

class Matrix
{
public:
  Matrix()
    : _11(1.0f), _12(0)
    , _21(0), _22(1.0f)
    , _31(0), _32(0)
  {}
  Matrix(Float a11, Float a12, Float a21, Float a22, Float a31, Float a32)
    : _11(a11), _12(a12)
    , _21(a21), _22(a22)
    , _31(a31), _32(a32)
  {}
  union {
    struct {
      Float _11, _12;
      Float _21, _22;
      Float _31, _32;
    };
    Float components[6];
  };

  MOZ_ALWAYS_INLINE Matrix Copy() const
  {
    return Matrix(*this);
  }

  friend std::ostream& operator<<(std::ostream& aStream, const Matrix& aMatrix);

  Point TransformPoint(const Point &aPoint) const
  {
    Point retPoint;

    retPoint.x = aPoint.x * _11 + aPoint.y * _21 + _31;
    retPoint.y = aPoint.x * _12 + aPoint.y * _22 + _32;

    return retPoint;
  }

  Size TransformSize(const Size &aSize) const
  {
    Size retSize;

    retSize.width = aSize.width * _11 + aSize.height * _21;
    retSize.height = aSize.width * _12 + aSize.height * _22;

    return retSize;
  }

  GFX2D_API Rect TransformBounds(const Rect& rect) const;

  static Matrix Translation(Float aX, Float aY)
  {
    return Matrix(1.0f, 0.0f, 0.0f, 1.0f, aX, aY);
  }

  static Matrix Translation(Point aPoint)
  {
    return Translation(aPoint.x, aPoint.y);
  }

  /**
   * Apply a translation to this matrix.
   *
   * The "Pre" in this method's name means that the translation is applied
   * -before- this matrix's existing transformation. That is, any vector that
   * is multiplied by the resulting matrix will first be translated, then be
   * transformed by the original transform.
   *
   * Calling this method will result in this matrix having the same value as
   * the result of:
   *
   *   Matrix::Translation(x, y) * this
   *
   * (Note that in performance critical code multiplying by the result of a
   * Translation()/Scaling() call is not recommended since that results in a
   * full matrix multiply involving 12 floating-point multiplications. Calling
   * this method would be preferred since it only involves four floating-point
   * multiplications.)
   */
  Matrix &PreTranslate(Float aX, Float aY)
  {
    _31 += _11 * aX + _21 * aY;
    _32 += _12 * aX + _22 * aY;

    return *this;
  }

  Matrix &PreTranslate(const Point &aPoint)
  {
    return PreTranslate(aPoint.x, aPoint.y);
  }

  /**
   * Similar to PreTranslate, but the translation is applied -after- this
   * matrix's existing transformation instead of before it.
   *
   * This method is generally less used than PreTranslate since typically code
   * want to adjust an existing user space to device space matrix to create a
   * transform to device space from a -new- user space (translated from the
   * previous user space). In that case consumers will need to use the Pre*
   * variants of the matrix methods rather than using the Post* methods, since
   * the Post* methods add a transform to the device space end of the
   * transformation.
   */
  Matrix &PostTranslate(Float aX, Float aY)
  {
    _31 += aX;
    _32 += aY;
    return *this;
  }

  Matrix &PostTranslate(const Point &aPoint)
  {
    return PostTranslate(aPoint.x, aPoint.y);
  }

  static Matrix Scaling(Float aScaleX, Float aScaleY)
  {
    return Matrix(aScaleX, 0.0f, 0.0f, aScaleY, 0.0f, 0.0f);
  }
  
  /**
   * Similar to PreTranslate, but applies a scale instead of a translation.
   */
  Matrix &PreScale(Float aX, Float aY)
  {
    _11 *= aX;
    _12 *= aX;
    _21 *= aY;
    _22 *= aY;

    return *this;
  }

  /**
   * Similar to PostTranslate, but applies a scale instead of a translation.
   */
  Matrix &PostScale(Float aScaleX, Float aScaleY)
  {
    _11 *= aScaleX;
    _12 *= aScaleY;
    _21 *= aScaleX;
    _22 *= aScaleY;
    _31 *= aScaleX;
    _32 *= aScaleY;

    return *this;
  }

  GFX2D_API static Matrix Rotation(Float aAngle);

  /**
   * Similar to PreTranslate, but applies a rotation instead of a translation.
   */
  Matrix &PreRotate(Float aAngle)
  {
    return *this = Matrix::Rotation(aAngle) * *this;
  }

  bool Invert()
  {
    // Compute co-factors.
    Float A = _22;
    Float B = -_21;
    Float C = _21 * _32 - _22 * _31;
    Float D = -_12;
    Float E = _11;
    Float F = _31 * _12 - _11 * _32;

    Float det = Determinant();

    if (!det) {
      return false;
    }

    Float inv_det = 1 / det;

    _11 = inv_det * A;
    _12 = inv_det * D;
    _21 = inv_det * B;
    _22 = inv_det * E;
    _31 = inv_det * C;
    _32 = inv_det * F;

    return true;
  }

  Matrix Inverse() const
  {
    Matrix clone = *this;
    DebugOnly<bool> inverted = clone.Invert();
    MOZ_ASSERT(inverted, "Attempted to get the inverse of a non-invertible matrix");
    return clone;
  }

  Float Determinant() const
  {
    return _11 * _22 - _12 * _21;
  }

  Matrix operator*(const Matrix &aMatrix) const
  {
    Matrix resultMatrix;

    resultMatrix._11 = this->_11 * aMatrix._11 + this->_12 * aMatrix._21;
    resultMatrix._12 = this->_11 * aMatrix._12 + this->_12 * aMatrix._22;
    resultMatrix._21 = this->_21 * aMatrix._11 + this->_22 * aMatrix._21;
    resultMatrix._22 = this->_21 * aMatrix._12 + this->_22 * aMatrix._22;
    resultMatrix._31 = this->_31 * aMatrix._11 + this->_32 * aMatrix._21 + aMatrix._31;
    resultMatrix._32 = this->_31 * aMatrix._12 + this->_32 * aMatrix._22 + aMatrix._32;

    return resultMatrix;
  }

  Matrix& operator*=(const Matrix &aMatrix)
  {
    *this = *this * aMatrix;
    return *this;
  }

  /**
   * Multiplies in the opposite order to operator=*.
   */
  Matrix &PreMultiply(const Matrix &aMatrix)
  {
    *this = aMatrix * *this;
    return *this;
  }

  /* Returns true if the other matrix is fuzzy-equal to this matrix.
   * Note that this isn't a cheap comparison!
   */
  bool operator==(const Matrix& other) const
  {
    return FuzzyEqual(_11, other._11) && FuzzyEqual(_12, other._12) &&
           FuzzyEqual(_21, other._21) && FuzzyEqual(_22, other._22) &&
           FuzzyEqual(_31, other._31) && FuzzyEqual(_32, other._32);
  }

  bool operator!=(const Matrix& other) const
  {
    return !(*this == other);
  }

  bool ExactlyEquals(const Matrix& o) const
  {
    return _11 == o._11 && _12 == o._12 &&
           _21 == o._21 && _22 == o._22 &&
           _31 == o._31 && _32 == o._32;
  }

  /* Verifies that the matrix contains no Infs or NaNs. */
  bool IsFinite() const
  {
    return mozilla::IsFinite(_11) && mozilla::IsFinite(_12) &&
           mozilla::IsFinite(_21) && mozilla::IsFinite(_22) &&
           mozilla::IsFinite(_31) && mozilla::IsFinite(_32);
  }

  /* Returns true if the matrix is a rectilinear transformation (i.e.
   * grid-aligned rectangles are transformed to grid-aligned rectangles)
   */
  bool IsRectilinear() const {
    if (FuzzyEqual(_12, 0) && FuzzyEqual(_21, 0)) {
      return true;
    } else if (FuzzyEqual(_22, 0) && FuzzyEqual(_11, 0)) {
      return true;
    }

    return false;
  }

  /**
   * Returns true if the matrix is anything other than a straight
   * translation by integers.
  */
  bool HasNonIntegerTranslation() const {
    return HasNonTranslation() ||
      !FuzzyEqual(_31, floor(_31 + Float(0.5))) ||
      !FuzzyEqual(_32, floor(_32 + Float(0.5)));
  }

  /**
   * Returns true if the matrix only has an integer translation.
   */
  bool HasOnlyIntegerTranslation() const {
    return !HasNonIntegerTranslation();
  }

  /**
   * Returns true if the matrix has any transform other
   * than a straight translation.
   */
  bool HasNonTranslation() const {
    return !FuzzyEqual(_11, 1.0) || !FuzzyEqual(_22, 1.0) ||
           !FuzzyEqual(_12, 0.0) || !FuzzyEqual(_21, 0.0);
  }

  /**
   * Returns true if the matrix has any transform other
   * than a translation or a -1 y scale (y axis flip)
   */
  bool HasNonTranslationOrFlip() const {
      return !FuzzyEqual(_11, 1.0) ||
             (!FuzzyEqual(_22, 1.0) && !FuzzyEqual(_22, -1.0)) ||
             !FuzzyEqual(_21, 0.0) || !FuzzyEqual(_12, 0.0);
  }

  /* Returns true if the matrix is an identity matrix.
   */
  bool IsIdentity() const
  {
    return _11 == 1.0f && _12 == 0.0f &&
           _21 == 0.0f && _22 == 1.0f &&
           _31 == 0.0f && _32 == 0.0f;
  }

  /* Returns true if the matrix is singular.
   */
  bool IsSingular() const
  {
    Float det = Determinant();
    return !mozilla::IsFinite(det) || det == 0;
  }

  GFX2D_API Matrix &NudgeToIntegers();

  bool IsTranslation() const
  {
    return FuzzyEqual(_11, 1.0f) && FuzzyEqual(_12, 0.0f) &&
           FuzzyEqual(_21, 0.0f) && FuzzyEqual(_22, 1.0f);
  }

  static bool FuzzyIsInteger(Float aValue)
  {
    return FuzzyEqual(aValue, floorf(aValue + 0.5f));
  }

  bool IsIntegerTranslation() const
  {
    return IsTranslation() && FuzzyIsInteger(_31) && FuzzyIsInteger(_32);
  }

  bool IsAllIntegers() const
  {
    return FuzzyIsInteger(_11) && FuzzyIsInteger(_12) &&
           FuzzyIsInteger(_21) && FuzzyIsInteger(_22) &&
           FuzzyIsInteger(_31) && FuzzyIsInteger(_32);
  }

  Point GetTranslation() const {
    return Point(_31, _32);
  }

  /**
   * Returns true if matrix is multiple of 90 degrees rotation with flipping,
   * scaling and translation.
   */
  bool PreservesAxisAlignedRectangles() const {
      return ((FuzzyEqual(_11, 0.0) && FuzzyEqual(_22, 0.0))
          || (FuzzyEqual(_12, 0.0) && FuzzyEqual(_21, 0.0)));
  }

  /**
   * Returns true if the matrix has any transform other
   * than a translation or scale; this is, if there is
   * rotation.
   */
  bool HasNonAxisAlignedTransform() const {
      return !FuzzyEqual(_21, 0.0) || !FuzzyEqual(_12, 0.0);
  }

  /**
   * Returns true if the matrix has negative scaling (i.e. flip).
   */
  bool HasNegativeScaling() const {
      return (_11 < 0.0) || (_22 < 0.0);
  }
};

// Helper functions used by Matrix4x4Typed defined in Matrix.cpp
double
SafeTangent(double aTheta);
double
FlushToZero(double aVal);

template<class Units, class F>
Point4DTyped<Units, F>
ComputePerspectivePlaneIntercept(const Point4DTyped<Units, F>& aFirst,
                                 const Point4DTyped<Units, F>& aSecond)
{
  // This function will always return a point with a w value of 0.
  // The X, Y, and Z components will point towards an infinite vanishing
  // point.

  // We want to interpolate aFirst and aSecond to find the point intersecting
  // with the w=0 plane.

  // Since we know what we want the w component to be, we can rearrange the
  // interpolation equation and solve for t.
  float t = -aFirst.w / (aSecond.w - aFirst.w);

  // Use t to find the remainder of the components
  return aFirst + (aSecond - aFirst) * t;
}


template <typename SourceUnits, typename TargetUnits>
class Matrix4x4Typed
{
public:
  typedef PointTyped<SourceUnits> SourcePoint;
  typedef PointTyped<TargetUnits> TargetPoint;
  typedef Point3DTyped<SourceUnits> SourcePoint3D;
  typedef Point3DTyped<TargetUnits> TargetPoint3D;
  typedef Point4DTyped<SourceUnits> SourcePoint4D;
  typedef Point4DTyped<TargetUnits> TargetPoint4D;
  typedef RectTyped<SourceUnits> SourceRect;
  typedef RectTyped<TargetUnits> TargetRect;

  Matrix4x4Typed()
    : _11(1.0f), _12(0.0f), _13(0.0f), _14(0.0f)
    , _21(0.0f), _22(1.0f), _23(0.0f), _24(0.0f)
    , _31(0.0f), _32(0.0f), _33(1.0f), _34(0.0f)
    , _41(0.0f), _42(0.0f), _43(0.0f), _44(1.0f)
  {}

  Matrix4x4Typed(Float a11, Float a12, Float a13, Float a14,
                 Float a21, Float a22, Float a23, Float a24,
                 Float a31, Float a32, Float a33, Float a34,
                 Float a41, Float a42, Float a43, Float a44)
    : _11(a11), _12(a12), _13(a13), _14(a14)
    , _21(a21), _22(a22), _23(a23), _24(a24)
    , _31(a31), _32(a32), _33(a33), _34(a34)
    , _41(a41), _42(a42), _43(a43), _44(a44)
  {}

  explicit Matrix4x4Typed(const Float aArray[16])
  {
    memcpy(components, aArray, sizeof(components));
  }

  Matrix4x4Typed(const Matrix4x4Typed& aOther)
  {
    memcpy(this, &aOther, sizeof(*this));
  }

  union {
    struct {
      Float _11, _12, _13, _14;
      Float _21, _22, _23, _24;
      Float _31, _32, _33, _34;
      Float _41, _42, _43, _44;
    };
    Float components[16];
  };

  friend std::ostream& operator<<(std::ostream& aStream, const Matrix4x4Typed& aMatrix)
  {
    const Float *f = &aMatrix._11;
    aStream << "[ " << f[0] << " "  << f[1] << " " << f[2] << " " << f[3] << " ;" << std::endl; f += 4;
    aStream << "  " << f[0] << " "  << f[1] << " " << f[2] << " " << f[3] << " ;" << std::endl; f += 4;
    aStream << "  " << f[0] << " "  << f[1] << " " << f[2] << " " << f[3] << " ;" << std::endl; f += 4;
    aStream << "  " << f[0] << " "  << f[1] << " " << f[2] << " " << f[3] << " ]" << std::endl;
    return aStream;
  }

  Point4D& operator[](int aIndex)
  {
      MOZ_ASSERT(aIndex >= 0 && aIndex <= 3, "Invalid matrix array index");
      return *reinterpret_cast<Point4D*>((&_11)+4*aIndex);
  }
  const Point4D& operator[](int aIndex) const
  {
      MOZ_ASSERT(aIndex >= 0 && aIndex <= 3, "Invalid matrix array index");
      return *reinterpret_cast<const Point4D*>((&_11)+4*aIndex);
  }

  /**
   * Returns true if the matrix is isomorphic to a 2D affine transformation.
   */
  bool Is2D() const
  {
    if (_13 != 0.0f || _14 != 0.0f ||
        _23 != 0.0f || _24 != 0.0f ||
        _31 != 0.0f || _32 != 0.0f || _33 != 1.0f || _34 != 0.0f ||
        _43 != 0.0f || _44 != 1.0f) {
      return false;
    }
    return true;
  }

  bool Is2D(Matrix* aMatrix) const {
    if (!Is2D()) {
      return false;
    }
    if (aMatrix) {
      aMatrix->_11 = _11;
      aMatrix->_12 = _12;
      aMatrix->_21 = _21;
      aMatrix->_22 = _22;
      aMatrix->_31 = _41;
      aMatrix->_32 = _42;
    }
    return true;
  }

  Matrix As2D() const
  {
    MOZ_ASSERT(Is2D(), "Matrix is not a 2D affine transform");

    return Matrix(_11, _12, _21, _22, _41, _42);
  }

  bool CanDraw2D(Matrix* aMatrix = nullptr) const {
    if (_14 != 0.0f ||
        _24 != 0.0f ||
        _44 != 1.0f) {
      return false;
    }
    if (aMatrix) {
      aMatrix->_11 = _11;
      aMatrix->_12 = _12;
      aMatrix->_21 = _21;
      aMatrix->_22 = _22;
      aMatrix->_31 = _41;
      aMatrix->_32 = _42;
    }
    return true;
  }

  Matrix4x4Typed& ProjectTo2D() {
    _31 = 0.0f;
    _32 = 0.0f;
    _13 = 0.0f;
    _23 = 0.0f;
    _33 = 1.0f;
    _43 = 0.0f;
    _34 = 0.0f;
    // Some matrices, such as those derived from perspective transforms,
    // can modify _44 from 1, while leaving the rest of the fourth column
    // (_14, _24) at 0. In this case, after resetting the third row and
    // third column above, the value of _44 functions only to scale the
    // coordinate transform divide by W. The matrix can be converted to
    // a true 2D matrix by normalizing out the scaling effect of _44 on
    // the remaining components ahead of time.
    if (_14 == 0.0f && _24 == 0.0f &&
        _44 != 1.0f && _44 != 0.0f) {
      Float scale = 1.0f / _44;
      _11 *= scale;
      _12 *= scale;
      _21 *= scale;
      _22 *= scale;
      _41 *= scale;
      _42 *= scale;
      _44 = 1.0f;
    }
    return *this;
  }

  template<class F>
  Point4DTyped<TargetUnits, F>
  ProjectPoint(const PointTyped<SourceUnits, F>& aPoint) const {
    // Find a value for z that will transform to 0.

    // The transformed value of z is computed as:
    // z' = aPoint.x * _13 + aPoint.y * _23 + z * _33 + _43;

    // Solving for z when z' = 0 gives us:
    F z = -(aPoint.x * _13 + aPoint.y * _23 + _43) / _33;

    // Compute the transformed point
    return this->TransformPoint(Point4DTyped<SourceUnits, F>(aPoint.x, aPoint.y, z, 1));
  }

  template<class F>
  RectTyped<TargetUnits, F>
  ProjectRectBounds(const RectTyped<SourceUnits, F>& aRect, const RectTyped<TargetUnits, F>& aClip) const
  {
    // This function must never return std::numeric_limits<Float>::max() or any
    // other arbitrary large value in place of inifinity.  This often occurs when
    // aRect is an inversed projection matrix or when aRect is transformed to be
    // partly behind and in front of the camera (w=0 plane in homogenous
    // coordinates) - See Bug 1035611

    // Some call-sites will call RoundGfxRectToAppRect which clips both the
    // extents and dimensions of the rect to be bounded by nscoord_MAX.
    // If we return a Rect that, when converted to nscoords, has a width or height
    // greater than nscoord_MAX, RoundGfxRectToAppRect will clip the overflow
    // off both the min and max end of the rect after clipping the extents of the
    // rect, resulting in a translation of the rect towards the infinite end.

    // The bounds returned by ProjectRectBounds are expected to be clipped only on
    // the edges beyond the bounds of the coordinate system; otherwise, the
    // clipped bounding box would be smaller than the correct one and result
    // bugs such as incorrect culling (eg. Bug 1073056)

    // To address this without requiring all code to work in homogenous
    // coordinates or interpret infinite values correctly, a specialized
    // clipping function is integrated into ProjectRectBounds.

    // Callers should pass an aClip value that represents the extents to clip
    // the result to, in the same coordinate system as aRect.
    Point4DTyped<TargetUnits, F> points[4];

    points[0] = ProjectPoint(aRect.TopLeft());
    points[1] = ProjectPoint(aRect.TopRight());
    points[2] = ProjectPoint(aRect.BottomRight());
    points[3] = ProjectPoint(aRect.BottomLeft());

    F min_x = std::numeric_limits<F>::max();
    F min_y = std::numeric_limits<F>::max();
    F max_x = -std::numeric_limits<F>::max();
    F max_y = -std::numeric_limits<F>::max();

    for (int i=0; i<4; i++) {
      // Only use points that exist above the w=0 plane
      if (points[i].HasPositiveWCoord()) {
        PointTyped<TargetUnits, F> point2d = aClip.ClampPoint(points[i].As2DPoint());
        min_x = std::min<F>(point2d.x, min_x);
        max_x = std::max<F>(point2d.x, max_x);
        min_y = std::min<F>(point2d.y, min_y);
        max_y = std::max<F>(point2d.y, max_y);
      }

      int next = (i == 3) ? 0 : i + 1;
      if (points[i].HasPositiveWCoord() != points[next].HasPositiveWCoord()) {
        // If the line between two points crosses the w=0 plane, then interpolate
        // to find the point of intersection with the w=0 plane and use that
        // instead.
        Point4DTyped<TargetUnits, F> intercept =
          ComputePerspectivePlaneIntercept(points[i], points[next]);
        // Since intercept.w will always be 0 here, we interpret x,y,z as a
        // direction towards an infinite vanishing point.
        if (intercept.x < 0.0f) {
          min_x = aClip.x;
        } else if (intercept.x > 0.0f) {
          max_x = aClip.XMost();
        }
        if (intercept.y < 0.0f) {
          min_y = aClip.y;
        } else if (intercept.y > 0.0f) {
          max_y = aClip.YMost();
        }
      }
    }

    if (max_x < min_x || max_y < min_y) {
      return RectTyped<TargetUnits, F>(0, 0, 0, 0);
    }

    return RectTyped<TargetUnits, F>(min_x, min_y, max_x - min_x, max_y - min_y);
  }

  /**
   * TransformAndClipBounds transforms aRect as a bounding box, while clipping
   * the transformed bounds to the extents of aClip.
   */
  template<class F>
  RectTyped<TargetUnits, F> TransformAndClipBounds(const RectTyped<SourceUnits, F>& aRect,
                                                   const RectTyped<TargetUnits, F>& aClip) const
  {
    PointTyped<UnknownUnits, F> verts[kTransformAndClipRectMaxVerts];
    size_t vertCount = TransformAndClipRect(aRect, aClip, verts);

    F min_x = std::numeric_limits<F>::max();
    F min_y = std::numeric_limits<F>::max();
    F max_x = -std::numeric_limits<F>::max();
    F max_y = -std::numeric_limits<F>::max();
    for (size_t i=0; i < vertCount; i++) {
      min_x = std::min(min_x, verts[i].x);
      max_x = std::max(max_x, verts[i].x);
      min_y = std::min(min_y, verts[i].y);
      max_y = std::max(max_y, verts[i].y);
    }

    if (max_x < min_x || max_y < min_y) {
      return RectTyped<TargetUnits, F>(0, 0, 0, 0);
    }

    return RectTyped<TargetUnits, F>(min_x, min_y, max_x - min_x, max_y - min_y);
  }

  template<class F>
  RectTyped<TargetUnits, F> TransformAndClipBounds(const TriangleTyped<SourceUnits, F>& aTriangle,
                                                   const RectTyped<TargetUnits, F>& aClip) const
  {
    return TransformAndClipBounds(aTriangle.BoundingBox(), aClip);
  }

  /**
   * TransformAndClipRect projects a rectangle and clips against view frustum
   * clipping planes in homogenous space so that its projected vertices are
   * constrained within the 2d rectangle passed in aClip.
   * The resulting vertices are populated in aVerts.  aVerts must be
   * pre-allocated to hold at least kTransformAndClipRectMaxVerts Points.
   * The vertex count is returned by TransformAndClipRect.  It is possible to
   * emit fewer that 3 vertices, indicating that aRect will not be visible
   * within aClip.
   */
  template<class F>
  size_t TransformAndClipRect(const RectTyped<SourceUnits, F>& aRect,
                              const RectTyped<TargetUnits, F>& aClip,
                              PointTyped<TargetUnits, F>* aVerts) const
  {
    // Initialize a double-buffered array of points in homogenous space with
    // the input rectangle, aRect.
    Point4DTyped<UnknownUnits, F> points[2][kTransformAndClipRectMaxVerts];
    Point4DTyped<UnknownUnits, F>* dstPoint = points[0];

    *dstPoint++ = TransformPoint(Point4DTyped<UnknownUnits, F>(aRect.x, aRect.y, 0, 1));
    *dstPoint++ = TransformPoint(Point4DTyped<UnknownUnits, F>(aRect.XMost(), aRect.y, 0, 1));
    *dstPoint++ = TransformPoint(Point4DTyped<UnknownUnits, F>(aRect.XMost(), aRect.YMost(), 0, 1));
    *dstPoint++ = TransformPoint(Point4DTyped<UnknownUnits, F>(aRect.x, aRect.YMost(), 0, 1));

    // View frustum clipping planes are described as normals originating from
    // the 0,0,0,0 origin.
    Point4DTyped<UnknownUnits, F> planeNormals[4];
    planeNormals[0] = Point4DTyped<UnknownUnits, F>(1.0, 0.0, 0.0, -aClip.x);
    planeNormals[1] = Point4DTyped<UnknownUnits, F>(-1.0, 0.0, 0.0, aClip.XMost());
    planeNormals[2] = Point4DTyped<UnknownUnits, F>(0.0, 1.0, 0.0, -aClip.y);
    planeNormals[3] = Point4DTyped<UnknownUnits, F>(0.0, -1.0, 0.0, aClip.YMost());

    // Iterate through each clipping plane and clip the polygon.
    // In each pass, we double buffer, alternating between points[0] and
    // points[1].
    for (int plane=0; plane < 4; plane++) {
      planeNormals[plane].Normalize();
      Point4DTyped<UnknownUnits, F>* srcPoint = points[plane & 1];
      Point4DTyped<UnknownUnits, F>* srcPointEnd = dstPoint;

      dstPoint = points[~plane & 1];
      Point4DTyped<UnknownUnits, F>* dstPointStart = dstPoint;

      Point4DTyped<UnknownUnits, F>* prevPoint = srcPointEnd - 1;
      F prevDot = planeNormals[plane].DotProduct(*prevPoint);
      while (srcPoint < srcPointEnd && ((dstPoint - dstPointStart) < kTransformAndClipRectMaxVerts)) {
        F nextDot = planeNormals[plane].DotProduct(*srcPoint);

        if ((nextDot >= 0.0) != (prevDot >= 0.0)) {
          // An intersection with the clipping plane has been detected.
          // Interpolate to find the intersecting point and emit it.
          F t = -prevDot / (nextDot - prevDot);
          *dstPoint++ = *srcPoint * t + *prevPoint * (1.0 - t);
        }

        if (nextDot >= 0.0) {
          // Emit any source points that are on the positive side of the
          // clipping plane.
          *dstPoint++ = *srcPoint;
        }

        prevPoint = srcPoint++;
        prevDot = nextDot;
      }

      if (dstPoint == dstPointStart) {
        break;
      }
    }

    size_t dstPointCount = 0;
    size_t srcPointCount = dstPoint - points[0];
    for (Point4DTyped<UnknownUnits, F>* srcPoint = points[0]; srcPoint < points[0] + srcPointCount; srcPoint++) {

      PointTyped<TargetUnits, F> p;
      if (srcPoint->w == 0.0) {
        // If a point lies on the intersection of the clipping planes at
        // (0,0,0,0), we must avoid a division by zero w component.
        p = PointTyped<TargetUnits, F>(0.0, 0.0);
      } else {
        p = srcPoint->As2DPoint();
      }
      // Emit only unique points
      if (dstPointCount == 0 || p != aVerts[dstPointCount - 1]) {
        aVerts[dstPointCount++] = p;
      }
    }

    return dstPointCount;
  }

  static const int kTransformAndClipRectMaxVerts = 32;

  static Matrix4x4Typed From2D(const Matrix &aMatrix) {
    Matrix4x4Typed matrix;
    matrix._11 = aMatrix._11;
    matrix._12 = aMatrix._12;
    matrix._21 = aMatrix._21;
    matrix._22 = aMatrix._22;
    matrix._41 = aMatrix._31;
    matrix._42 = aMatrix._32;
    return matrix;
  }

  bool Is2DIntegerTranslation() const
  {
    return Is2D() && As2D().IsIntegerTranslation();
  }

  TargetPoint4D TransposeTransform4D(const SourcePoint4D& aPoint) const
  {
      Float x = aPoint.x * _11 + aPoint.y * _12 + aPoint.z * _13 + aPoint.w * _14;
      Float y = aPoint.x * _21 + aPoint.y * _22 + aPoint.z * _23 + aPoint.w * _24;
      Float z = aPoint.x * _31 + aPoint.y * _32 + aPoint.z * _33 + aPoint.w * _34;
      Float w = aPoint.x * _41 + aPoint.y * _42 + aPoint.z * _43 + aPoint.w * _44;

      return TargetPoint4D(x, y, z, w);
  }

  template<class F>
  Point4DTyped<TargetUnits, F> TransformPoint(const Point4DTyped<SourceUnits, F>& aPoint) const
  {
    Point4DTyped<TargetUnits, F> retPoint;

    retPoint.x = aPoint.x * _11 + aPoint.y * _21 + aPoint.z * _31 + aPoint.w * _41;
    retPoint.y = aPoint.x * _12 + aPoint.y * _22 + aPoint.z * _32 + aPoint.w * _42;
    retPoint.z = aPoint.x * _13 + aPoint.y * _23 + aPoint.z * _33 + aPoint.w * _43;
    retPoint.w = aPoint.x * _14 + aPoint.y * _24 + aPoint.z * _34 + aPoint.w * _44;

    return retPoint;
  }

  template<class F>
  Point3DTyped<TargetUnits, F> TransformPoint(const Point3DTyped<SourceUnits, F>& aPoint) const
  {
    Point3DTyped<TargetUnits, F> result;
    result.x = aPoint.x * _11 + aPoint.y * _21 + aPoint.z * _31 + _41;
    result.y = aPoint.x * _12 + aPoint.y * _22 + aPoint.z * _32 + _42;
    result.z = aPoint.x * _13 + aPoint.y * _23 + aPoint.z * _33 + _43;

    result /= (aPoint.x * _14 + aPoint.y * _24 + aPoint.z * _34 + _44);

    return result;
  }

  template<class F>
  PointTyped<TargetUnits, F> TransformPoint(const PointTyped<SourceUnits, F> &aPoint) const
  {
    Point4DTyped<SourceUnits, F> temp(aPoint.x, aPoint.y, 0, 1);
    return TransformPoint(temp).As2DPoint();
  }

  template<class F>
  GFX2D_API RectTyped<TargetUnits, F> TransformBounds(const RectTyped<SourceUnits, F>& aRect) const
  {
    Point4DTyped<TargetUnits, F> verts[4];
    verts[0] = TransformPoint(Point4DTyped<SourceUnits, F>(aRect.x, aRect.y, 0.0, 1.0));
    verts[1] = TransformPoint(Point4DTyped<SourceUnits, F>(aRect.XMost(), aRect.y, 0.0, 1.0));
    verts[2] = TransformPoint(Point4DTyped<SourceUnits, F>(aRect.XMost(), aRect.YMost(), 0.0, 1.0));
    verts[3] = TransformPoint(Point4DTyped<SourceUnits, F>(aRect.x, aRect.YMost(), 0.0, 1.0));

    PointTyped<TargetUnits, F> quad[4];
    F min_x, max_x;
    F min_y, max_y;

    quad[0] = TransformPoint(aRect.TopLeft());
    quad[1] = TransformPoint(aRect.TopRight());
    quad[2] = TransformPoint(aRect.BottomLeft());
    quad[3] = TransformPoint(aRect.BottomRight());

    min_x = max_x = quad[0].x;
    min_y = max_y = quad[0].y;

    for (int i = 1; i < 4; i++) {
      if (quad[i].x < min_x) {
        min_x = quad[i].x;
      }
      if (quad[i].x > max_x) {
        max_x = quad[i].x;
      }

      if (quad[i].y < min_y) {
        min_y = quad[i].y;
      }
      if (quad[i].y > max_y) {
        max_y = quad[i].y;
      }
    }

    return RectTyped<TargetUnits, F>(min_x, min_y, max_x - min_x, max_y - min_y);
  }

  static Matrix4x4Typed Translation(Float aX, Float aY, Float aZ)
  {
    return Matrix4x4Typed(1.0f, 0.0f, 0.0f, 0.0f,
                          0.0f, 1.0f, 0.0f, 0.0f,
                          0.0f, 0.0f, 1.0f, 0.0f,
                          aX,   aY,   aZ, 1.0f);
  }

  static Matrix4x4Typed Translation(const TargetPoint3D& aP)
  {
    return Translation(aP.x, aP.y, aP.z);
  }

  static Matrix4x4Typed Translation(const TargetPoint& aP)
  {
    return Translation(aP.x, aP.y, 0);
  }

  /**
   * Apply a translation to this matrix.
   *
   * The "Pre" in this method's name means that the translation is applied
   * -before- this matrix's existing transformation. That is, any vector that
   * is multiplied by the resulting matrix will first be translated, then be
   * transformed by the original transform.
   *
   * Calling this method will result in this matrix having the same value as
   * the result of:
   *
   *   Matrix4x4::Translation(x, y) * this
   *
   * (Note that in performance critical code multiplying by the result of a
   * Translation()/Scaling() call is not recommended since that results in a
   * full matrix multiply involving 64 floating-point multiplications. Calling
   * this method would be preferred since it only involves 12 floating-point
   * multiplications.)
   */
  Matrix4x4Typed &PreTranslate(Float aX, Float aY, Float aZ)
  {
    _41 += aX * _11 + aY * _21 + aZ * _31;
    _42 += aX * _12 + aY * _22 + aZ * _32;
    _43 += aX * _13 + aY * _23 + aZ * _33;
    _44 += aX * _14 + aY * _24 + aZ * _34;

    return *this;
  }

  Matrix4x4Typed &PreTranslate(const Point3D& aPoint) {
    return PreTranslate(aPoint.x, aPoint.y, aPoint.z);
  }

  /**
   * Similar to PreTranslate, but the translation is applied -after- this
   * matrix's existing transformation instead of before it.
   *
   * This method is generally less used than PreTranslate since typically code
   * wants to adjust an existing user space to device space matrix to create a
   * transform to device space from a -new- user space (translated from the
   * previous user space). In that case consumers will need to use the Pre*
   * variants of the matrix methods rather than using the Post* methods, since
   * the Post* methods add a transform to the device space end of the
   * transformation.
   */
  Matrix4x4Typed &PostTranslate(Float aX, Float aY, Float aZ)
  {
    _11 += _14 * aX;
    _21 += _24 * aX;
    _31 += _34 * aX;
    _41 += _44 * aX;
    _12 += _14 * aY;
    _22 += _24 * aY;
    _32 += _34 * aY;
    _42 += _44 * aY;
    _13 += _14 * aZ;
    _23 += _24 * aZ;
    _33 += _34 * aZ;
    _43 += _44 * aZ;

    return *this;
  }

  Matrix4x4Typed &PostTranslate(const TargetPoint3D& aPoint) {
    return PostTranslate(aPoint.x, aPoint.y, aPoint.z);
  }

  Matrix4x4Typed &PostTranslate(const TargetPoint& aPoint) {
    return PostTranslate(aPoint.x, aPoint.y, 0);
  }

  static Matrix4x4Typed Scaling(Float aScaleX, Float aScaleY, float aScaleZ)
  {
    return Matrix4x4Typed(aScaleX, 0.0f, 0.0f, 0.0f,
                          0.0f, aScaleY, 0.0f, 0.0f,
                          0.0f, 0.0f, aScaleZ, 0.0f,
                          0.0f, 0.0f, 0.0f, 1.0f);
  }

  /**
   * Similar to PreTranslate, but applies a scale instead of a translation.
   */
  Matrix4x4Typed &PreScale(Float aX, Float aY, Float aZ)
  {
    _11 *= aX;
    _12 *= aX;
    _13 *= aX;
    _14 *= aX;
    _21 *= aY;
    _22 *= aY;
    _23 *= aY;
    _24 *= aY;
    _31 *= aZ;
    _32 *= aZ;
    _33 *= aZ;
    _34 *= aZ;

    return *this;
  }

  /**
   * Similar to PostTranslate, but applies a scale instead of a translation.
   */
  Matrix4x4Typed &PostScale(Float aScaleX, Float aScaleY, Float aScaleZ)
  {
    _11 *= aScaleX;
    _21 *= aScaleX;
    _31 *= aScaleX;
    _41 *= aScaleX;
    _12 *= aScaleY;
    _22 *= aScaleY;
    _32 *= aScaleY;
    _42 *= aScaleY;
    _13 *= aScaleZ;
    _23 *= aScaleZ;
    _33 *= aScaleZ;
    _43 *= aScaleZ;

    return *this;
  }

  void SkewXY(Float aSkew)
  {
    (*this)[1] += (*this)[0] * aSkew;
  }

  void SkewXZ(Float aSkew)
  {
      (*this)[2] += (*this)[0] * aSkew;
  }

  void SkewYZ(Float aSkew)
  {
      (*this)[2] += (*this)[1] * aSkew;
  }

  Matrix4x4Typed &ChangeBasis(const Point3D& aOrigin)
  {
    return ChangeBasis(aOrigin.x, aOrigin.y, aOrigin.z);
  }

  Matrix4x4Typed &ChangeBasis(Float aX, Float aY, Float aZ)
  {
    // Translate to the origin before applying this matrix
    PreTranslate(-aX, -aY, -aZ);

    // Translate back into position after applying this matrix
    PostTranslate(aX, aY, aZ);

    return *this;
  }

  Matrix4x4Typed& Transpose() {
    std::swap(_12, _21);
    std::swap(_13, _31);
    std::swap(_14, _41);

    std::swap(_23, _32);
    std::swap(_24, _42);

    std::swap(_34, _43);

    return *this;
  }

  bool operator==(const Matrix4x4Typed& o) const
  {
    // XXX would be nice to memcmp here, but that breaks IEEE 754 semantics
    return _11 == o._11 && _12 == o._12 && _13 == o._13 && _14 == o._14 &&
           _21 == o._21 && _22 == o._22 && _23 == o._23 && _24 == o._24 &&
           _31 == o._31 && _32 == o._32 && _33 == o._33 && _34 == o._34 &&
           _41 == o._41 && _42 == o._42 && _43 == o._43 && _44 == o._44;
  }

  bool operator!=(const Matrix4x4Typed& o) const
  {
    return !((*this) == o);
  }

  template <typename NewTargetUnits>
  Matrix4x4Typed<SourceUnits, NewTargetUnits> operator*(const Matrix4x4Typed<TargetUnits, NewTargetUnits> &aMatrix) const
  {
    Matrix4x4Typed<SourceUnits, NewTargetUnits> matrix;

    matrix._11 = _11 * aMatrix._11 + _12 * aMatrix._21 + _13 * aMatrix._31 + _14 * aMatrix._41;
    matrix._21 = _21 * aMatrix._11 + _22 * aMatrix._21 + _23 * aMatrix._31 + _24 * aMatrix._41;
    matrix._31 = _31 * aMatrix._11 + _32 * aMatrix._21 + _33 * aMatrix._31 + _34 * aMatrix._41;
    matrix._41 = _41 * aMatrix._11 + _42 * aMatrix._21 + _43 * aMatrix._31 + _44 * aMatrix._41;
    matrix._12 = _11 * aMatrix._12 + _12 * aMatrix._22 + _13 * aMatrix._32 + _14 * aMatrix._42;
    matrix._22 = _21 * aMatrix._12 + _22 * aMatrix._22 + _23 * aMatrix._32 + _24 * aMatrix._42;
    matrix._32 = _31 * aMatrix._12 + _32 * aMatrix._22 + _33 * aMatrix._32 + _34 * aMatrix._42;
    matrix._42 = _41 * aMatrix._12 + _42 * aMatrix._22 + _43 * aMatrix._32 + _44 * aMatrix._42;
    matrix._13 = _11 * aMatrix._13 + _12 * aMatrix._23 + _13 * aMatrix._33 + _14 * aMatrix._43;
    matrix._23 = _21 * aMatrix._13 + _22 * aMatrix._23 + _23 * aMatrix._33 + _24 * aMatrix._43;
    matrix._33 = _31 * aMatrix._13 + _32 * aMatrix._23 + _33 * aMatrix._33 + _34 * aMatrix._43;
    matrix._43 = _41 * aMatrix._13 + _42 * aMatrix._23 + _43 * aMatrix._33 + _44 * aMatrix._43;
    matrix._14 = _11 * aMatrix._14 + _12 * aMatrix._24 + _13 * aMatrix._34 + _14 * aMatrix._44;
    matrix._24 = _21 * aMatrix._14 + _22 * aMatrix._24 + _23 * aMatrix._34 + _24 * aMatrix._44;
    matrix._34 = _31 * aMatrix._14 + _32 * aMatrix._24 + _33 * aMatrix._34 + _34 * aMatrix._44;
    matrix._44 = _41 * aMatrix._14 + _42 * aMatrix._24 + _43 * aMatrix._34 + _44 * aMatrix._44;

    return matrix;
  }

  Matrix4x4Typed& operator*=(const Matrix4x4Typed<TargetUnits, TargetUnits> &aMatrix)
  {
    *this = *this * aMatrix;
    return *this;
  }

  /* Returns true if the matrix is an identity matrix.
   */
  bool IsIdentity() const
  {
    return _11 == 1.0f && _12 == 0.0f && _13 == 0.0f && _14 == 0.0f &&
           _21 == 0.0f && _22 == 1.0f && _23 == 0.0f && _24 == 0.0f &&
           _31 == 0.0f && _32 == 0.0f && _33 == 1.0f && _34 == 0.0f &&
           _41 == 0.0f && _42 == 0.0f && _43 == 0.0f && _44 == 1.0f;
  }

  bool IsSingular() const
  {
    return Determinant() == 0.0;
  }

  Float Determinant() const
  {
    return _14 * _23 * _32 * _41
         - _13 * _24 * _32 * _41
         - _14 * _22 * _33 * _41
         + _12 * _24 * _33 * _41
         + _13 * _22 * _34 * _41
         - _12 * _23 * _34 * _41
         - _14 * _23 * _31 * _42
         + _13 * _24 * _31 * _42
         + _14 * _21 * _33 * _42
         - _11 * _24 * _33 * _42
         - _13 * _21 * _34 * _42
         + _11 * _23 * _34 * _42
         + _14 * _22 * _31 * _43
         - _12 * _24 * _31 * _43
         - _14 * _21 * _32 * _43
         + _11 * _24 * _32 * _43
         + _12 * _21 * _34 * _43
         - _11 * _22 * _34 * _43
         - _13 * _22 * _31 * _44
         + _12 * _23 * _31 * _44
         + _13 * _21 * _32 * _44
         - _11 * _23 * _32 * _44
         - _12 * _21 * _33 * _44
         + _11 * _22 * _33 * _44;
  }

  // Invert() is not unit-correct. Prefer Inverse() where possible.
  bool Invert()
  {
    Float det = Determinant();
    if (!det) {
      return false;
    }

    Matrix4x4Typed<SourceUnits, TargetUnits> result;
    result._11 = _23 * _34 * _42 - _24 * _33 * _42 + _24 * _32 * _43 - _22 * _34 * _43 - _23 * _32 * _44 + _22 * _33 * _44;
    result._12 = _14 * _33 * _42 - _13 * _34 * _42 - _14 * _32 * _43 + _12 * _34 * _43 + _13 * _32 * _44 - _12 * _33 * _44;
    result._13 = _13 * _24 * _42 - _14 * _23 * _42 + _14 * _22 * _43 - _12 * _24 * _43 - _13 * _22 * _44 + _12 * _23 * _44;
    result._14 = _14 * _23 * _32 - _13 * _24 * _32 - _14 * _22 * _33 + _12 * _24 * _33 + _13 * _22 * _34 - _12 * _23 * _34;
    result._21 = _24 * _33 * _41 - _23 * _34 * _41 - _24 * _31 * _43 + _21 * _34 * _43 + _23 * _31 * _44 - _21 * _33 * _44;
    result._22 = _13 * _34 * _41 - _14 * _33 * _41 + _14 * _31 * _43 - _11 * _34 * _43 - _13 * _31 * _44 + _11 * _33 * _44;
    result._23 = _14 * _23 * _41 - _13 * _24 * _41 - _14 * _21 * _43 + _11 * _24 * _43 + _13 * _21 * _44 - _11 * _23 * _44;
    result._24 = _13 * _24 * _31 - _14 * _23 * _31 + _14 * _21 * _33 - _11 * _24 * _33 - _13 * _21 * _34 + _11 * _23 * _34;
    result._31 = _22 * _34 * _41 - _24 * _32 * _41 + _24 * _31 * _42 - _21 * _34 * _42 - _22 * _31 * _44 + _21 * _32 * _44;
    result._32 = _14 * _32 * _41 - _12 * _34 * _41 - _14 * _31 * _42 + _11 * _34 * _42 + _12 * _31 * _44 - _11 * _32 * _44;
    result._33 = _12 * _24 * _41 - _14 * _22 * _41 + _14 * _21 * _42 - _11 * _24 * _42 - _12 * _21 * _44 + _11 * _22 * _44;
    result._34 = _14 * _22 * _31 - _12 * _24 * _31 - _14 * _21 * _32 + _11 * _24 * _32 + _12 * _21 * _34 - _11 * _22 * _34;
    result._41 = _23 * _32 * _41 - _22 * _33 * _41 - _23 * _31 * _42 + _21 * _33 * _42 + _22 * _31 * _43 - _21 * _32 * _43;
    result._42 = _12 * _33 * _41 - _13 * _32 * _41 + _13 * _31 * _42 - _11 * _33 * _42 - _12 * _31 * _43 + _11 * _32 * _43;
    result._43 = _13 * _22 * _41 - _12 * _23 * _41 - _13 * _21 * _42 + _11 * _23 * _42 + _12 * _21 * _43 - _11 * _22 * _43;
    result._44 = _12 * _23 * _31 - _13 * _22 * _31 + _13 * _21 * _32 - _11 * _23 * _32 - _12 * _21 * _33 + _11 * _22 * _33;

    result._11 /= det;
    result._12 /= det;
    result._13 /= det;
    result._14 /= det;
    result._21 /= det;
    result._22 /= det;
    result._23 /= det;
    result._24 /= det;
    result._31 /= det;
    result._32 /= det;
    result._33 /= det;
    result._34 /= det;
    result._41 /= det;
    result._42 /= det;
    result._43 /= det;
    result._44 /= det;
    *this = result;

    return true;
  }

  Matrix4x4Typed<TargetUnits, SourceUnits> Inverse() const
  {
    typedef Matrix4x4Typed<TargetUnits, SourceUnits> InvertedMatrix;
    InvertedMatrix clone = InvertedMatrix::FromUnknownMatrix(ToUnknownMatrix());
    DebugOnly<bool> inverted = clone.Invert();
    MOZ_ASSERT(inverted, "Attempted to get the inverse of a non-invertible matrix");
    return clone;
  }

  void Normalize()
  {
      for (int i = 0; i < 4; i++) {
          for (int j = 0; j < 4; j++) {
              (*this)[i][j] /= (*this)[3][3];
         }
      }
  }

  bool FuzzyEqual(const Matrix4x4Typed& o) const
  {
    return gfx::FuzzyEqual(_11, o._11) && gfx::FuzzyEqual(_12, o._12) &&
           gfx::FuzzyEqual(_13, o._13) && gfx::FuzzyEqual(_14, o._14) &&
           gfx::FuzzyEqual(_21, o._21) && gfx::FuzzyEqual(_22, o._22) &&
           gfx::FuzzyEqual(_23, o._23) && gfx::FuzzyEqual(_24, o._24) &&
           gfx::FuzzyEqual(_31, o._31) && gfx::FuzzyEqual(_32, o._32) &&
           gfx::FuzzyEqual(_33, o._33) && gfx::FuzzyEqual(_34, o._34) &&
           gfx::FuzzyEqual(_41, o._41) && gfx::FuzzyEqual(_42, o._42) &&
           gfx::FuzzyEqual(_43, o._43) && gfx::FuzzyEqual(_44, o._44);
  }

  bool FuzzyEqualsMultiplicative(const Matrix4x4Typed& o) const
  {
    return ::mozilla::FuzzyEqualsMultiplicative(_11, o._11) &&
           ::mozilla::FuzzyEqualsMultiplicative(_12, o._12) &&
           ::mozilla::FuzzyEqualsMultiplicative(_13, o._13) &&
           ::mozilla::FuzzyEqualsMultiplicative(_14, o._14) &&
           ::mozilla::FuzzyEqualsMultiplicative(_21, o._21) &&
           ::mozilla::FuzzyEqualsMultiplicative(_22, o._22) &&
           ::mozilla::FuzzyEqualsMultiplicative(_23, o._23) &&
           ::mozilla::FuzzyEqualsMultiplicative(_24, o._24) &&
           ::mozilla::FuzzyEqualsMultiplicative(_31, o._31) &&
           ::mozilla::FuzzyEqualsMultiplicative(_32, o._32) &&
           ::mozilla::FuzzyEqualsMultiplicative(_33, o._33) &&
           ::mozilla::FuzzyEqualsMultiplicative(_34, o._34) &&
           ::mozilla::FuzzyEqualsMultiplicative(_41, o._41) &&
           ::mozilla::FuzzyEqualsMultiplicative(_42, o._42) &&
           ::mozilla::FuzzyEqualsMultiplicative(_43, o._43) &&
           ::mozilla::FuzzyEqualsMultiplicative(_44, o._44);
  }

  bool IsBackfaceVisible() const
  {
    // Inverse()._33 < 0;
    Float det = Determinant();
    Float __33 = _12*_24*_41 - _14*_22*_41 +
                _14*_21*_42 - _11*_24*_42 -
                _12*_21*_44 + _11*_22*_44;
    return (__33 * det) < 0;
  }

  Matrix4x4Typed &NudgeToIntegersFixedEpsilon()
  {
    NudgeToInteger(&_11);
    NudgeToInteger(&_12);
    NudgeToInteger(&_13);
    NudgeToInteger(&_14);
    NudgeToInteger(&_21);
    NudgeToInteger(&_22);
    NudgeToInteger(&_23);
    NudgeToInteger(&_24);
    NudgeToInteger(&_31);
    NudgeToInteger(&_32);
    NudgeToInteger(&_33);
    NudgeToInteger(&_34);
    static const float error = 1e-5f;
    NudgeToInteger(&_41, error);
    NudgeToInteger(&_42, error);
    NudgeToInteger(&_43, error);
    NudgeToInteger(&_44, error);
    return *this;
  }

  Point4D TransposedVector(int aIndex) const
  {
      MOZ_ASSERT(aIndex >= 0 && aIndex <= 3, "Invalid matrix array index");
      return Point4D(*((&_11)+aIndex), *((&_21)+aIndex), *((&_31)+aIndex), *((&_41)+aIndex));
  }

  void SetTransposedVector(int aIndex, Point4D &aVector)
  {
      MOZ_ASSERT(aIndex >= 0 && aIndex <= 3, "Invalid matrix array index");
      *((&_11)+aIndex) = aVector.x;
      *((&_21)+aIndex) = aVector.y;
      *((&_31)+aIndex) = aVector.z;
      *((&_41)+aIndex) = aVector.w;
  }

  // Sets this matrix to a rotation matrix given by aQuat.
  // This quaternion *MUST* be normalized!
  // Implemented in Quaternion.cpp
  void SetRotationFromQuaternion(const Quaternion& q)
  {
    const Float x2 = q.x + q.x, y2 = q.y + q.y, z2 = q.z + q.z;
    const Float xx = q.x * x2, xy = q.x * y2, xz = q.x * z2;
    const Float yy = q.y * y2, yz = q.y * z2, zz = q.z * z2;
    const Float wx = q.w * x2, wy = q.w * y2, wz = q.w * z2;

    _11 = 1.0f - (yy + zz);
    _21 = xy + wz;
    _31 = xz - wy;
    _41 = 0.0f;

    _12 = xy - wz;
    _22 = 1.0f - (xx + zz);
    _32 = yz + wx;
    _42 = 0.0f;

    _13 = xz + wy;
    _23 = yz - wx;
    _33 = 1.0f - (xx + yy);
    _43 = 0.0f;

    _14 = _42 = _43 = 0.0f;
    _44 = 1.0f;
  }

  // Set all the members of the matrix to NaN
  void SetNAN()
  {
    _11 = UnspecifiedNaN<Float>();
    _21 = UnspecifiedNaN<Float>();
    _31 = UnspecifiedNaN<Float>();
    _41 = UnspecifiedNaN<Float>();
    _12 = UnspecifiedNaN<Float>();
    _22 = UnspecifiedNaN<Float>();
    _32 = UnspecifiedNaN<Float>();
    _42 = UnspecifiedNaN<Float>();
    _13 = UnspecifiedNaN<Float>();
    _23 = UnspecifiedNaN<Float>();
    _33 = UnspecifiedNaN<Float>();
    _43 = UnspecifiedNaN<Float>();
    _14 = UnspecifiedNaN<Float>();
    _24 = UnspecifiedNaN<Float>();
    _34 = UnspecifiedNaN<Float>();
    _44 = UnspecifiedNaN<Float>();
  }

  void SkewXY(double aXSkew, double aYSkew)
  {
    // XXX Is double precision really necessary here
    float tanX = SafeTangent(aXSkew);
    float tanY = SafeTangent(aYSkew);
    float temp;

    temp = _11;
    _11 += tanY * _21;
    _21 += tanX * temp;

    temp = _12;
    _12 += tanY * _22;
    _22 += tanX * temp;

    temp = _13;
    _13 += tanY * _23;
    _23 += tanX * temp;

    temp = _14;
    _14 += tanY * _24;
    _24 += tanX * temp;
  }

  void RotateX(double aTheta)
  {
    // XXX Is double precision really necessary here
    double cosTheta = FlushToZero(cos(aTheta));
    double sinTheta = FlushToZero(sin(aTheta));

    float temp;

    temp = _21;
    _21 = cosTheta * _21 + sinTheta * _31;
    _31 = -sinTheta * temp + cosTheta * _31;

    temp = _22;
    _22 = cosTheta * _22 + sinTheta * _32;
    _32 = -sinTheta * temp + cosTheta * _32;

    temp = _23;
    _23 = cosTheta * _23 + sinTheta * _33;
    _33 = -sinTheta * temp + cosTheta * _33;

    temp = _24;
    _24 = cosTheta * _24 + sinTheta * _34;
    _34 = -sinTheta * temp + cosTheta * _34;
  }

  void RotateY(double aTheta)
  {
    // XXX Is double precision really necessary here
    double cosTheta = FlushToZero(cos(aTheta));
    double sinTheta = FlushToZero(sin(aTheta));

    float temp;

    temp = _11;
    _11 = cosTheta * _11 + -sinTheta * _31;
    _31 = sinTheta * temp + cosTheta * _31;

    temp = _12;
    _12 = cosTheta * _12 + -sinTheta * _32;
    _32 = sinTheta * temp + cosTheta * _32;

    temp = _13;
    _13 = cosTheta * _13 + -sinTheta * _33;
    _33 = sinTheta * temp + cosTheta * _33;

    temp = _14;
    _14 = cosTheta * _14 + -sinTheta * _34;
    _34 = sinTheta * temp + cosTheta * _34;
  }

  void RotateZ(double aTheta)
  {
    // XXX Is double precision really necessary here
    double cosTheta = FlushToZero(cos(aTheta));
    double sinTheta = FlushToZero(sin(aTheta));

    float temp;

    temp = _11;
    _11 = cosTheta * _11 + sinTheta * _21;
    _21 = -sinTheta * temp + cosTheta * _21;

    temp = _12;
    _12 = cosTheta * _12 + sinTheta * _22;
    _22 = -sinTheta * temp + cosTheta * _22;

    temp = _13;
    _13 = cosTheta * _13 + sinTheta * _23;
    _23 = -sinTheta * temp + cosTheta * _23;

    temp = _14;
    _14 = cosTheta * _14 + sinTheta * _24;
    _24 = -sinTheta * temp + cosTheta * _24;
  }

  // Sets this matrix to a rotation matrix about a
  // vector [x,y,z] by angle theta. The vector is normalized
  // to a unit vector.
  // https://www.w3.org/TR/css3-3d-transforms/#Rotate3dDefined
  void SetRotateAxisAngle(double aX, double aY, double aZ, double aTheta)
  {
    Point3D vector(aX, aY, aZ);
    if (!vector.Length()) {
      return;
    }
    vector.Normalize();

    double x = vector.x;
    double y = vector.y;
    double z = vector.z;

    double cosTheta = FlushToZero(cos(aTheta));
    double sinTheta = FlushToZero(sin(aTheta));

    // sin(aTheta / 2) * cos(aTheta / 2)
    double sc = sinTheta / 2;
    // pow(sin(aTheta / 2), 2)
    double sq = (1 - cosTheta) / 2;

    _11 = 1 - 2 * (y * y + z * z) * sq;
    _12 = 2 * (x * y * sq + z * sc);
    _13 = 2 * (x * z * sq - y * sc);
    _14 = 0.0f;
    _21 = 2 * (x * y * sq - z * sc);
    _22 = 1 - 2 * (x * x + z * z) * sq;
    _23 = 2 * (y * z * sq + x * sc);
    _24 = 0.0f;
    _31 = 2 * (x * z * sq + y * sc);
    _32 = 2 * (y * z * sq - x * sc);
    _33 = 1 - 2 * (x * x + y * y) * sq;
    _34 = 0.0f;
    _41 = 0.0f;
    _42 = 0.0f;
    _43 = 0.0f;
    _44 = 1.0f;
  }

  void Perspective(float aDepth)
  {
    MOZ_ASSERT(aDepth > 0.0f, "Perspective must be positive!");
    _31 += -1.0/aDepth * _41;
    _32 += -1.0/aDepth * _42;
    _33 += -1.0/aDepth * _43;
    _34 += -1.0/aDepth * _44;
  }

  Point3D GetNormalVector() const
  {
    // Define a plane in transformed space as the transformations
    // of 3 points on the z=0 screen plane.
    Point3D a = TransformPoint(Point3D(0, 0, 0));
    Point3D b = TransformPoint(Point3D(0, 1, 0));
    Point3D c = TransformPoint(Point3D(1, 0, 0));

    // Convert to two vectors on the surface of the plane.
    Point3D ab = b - a;
    Point3D ac = c - a;

    return ac.CrossProduct(ab);
  }

  /**
   * Returns true if the matrix has any transform other
   * than a straight translation.
   */
  bool HasNonTranslation() const {
    return !gfx::FuzzyEqual(_11, 1.0) || !gfx::FuzzyEqual(_22, 1.0) ||
           !gfx::FuzzyEqual(_12, 0.0) || !gfx::FuzzyEqual(_21, 0.0) ||
           !gfx::FuzzyEqual(_13, 0.0) || !gfx::FuzzyEqual(_23, 0.0) ||
           !gfx::FuzzyEqual(_31, 0.0) || !gfx::FuzzyEqual(_32, 0.0) ||
           !gfx::FuzzyEqual(_33, 1.0);
  }

  /**
   * Returns true if the matrix is anything other than a straight
   * translation by integers.
  */
  bool HasNonIntegerTranslation() const {
    return HasNonTranslation() ||
      !gfx::FuzzyEqual(_41, floor(_41 + 0.5)) ||
      !gfx::FuzzyEqual(_42, floor(_42 + 0.5)) ||
      !gfx::FuzzyEqual(_43, floor(_43 + 0.5));
  }

  /**
   * Return true if the matrix is with perspective (w).
   */
  bool HasPerspectiveComponent() const {
    return _14 != 0 || _24 != 0 || _34 != 0 || _44 != 1;
  }

  /**
   * Convert between typed and untyped matrices.
   */
  Matrix4x4 ToUnknownMatrix() const {
    return Matrix4x4{_11, _12, _13, _14,
                     _21, _22, _23, _24,
                     _31, _32, _33, _34,
                     _41, _42, _43, _44};
  }
  static Matrix4x4Typed FromUnknownMatrix(const Matrix4x4& aUnknown) {
    return Matrix4x4Typed{aUnknown._11, aUnknown._12, aUnknown._13, aUnknown._14,
                          aUnknown._21, aUnknown._22, aUnknown._23, aUnknown._24,
                          aUnknown._31, aUnknown._32, aUnknown._33, aUnknown._34,
                          aUnknown._41, aUnknown._42, aUnknown._43, aUnknown._44};
  }
};

typedef Matrix4x4Typed<UnknownUnits, UnknownUnits> Matrix4x4;

class Matrix5x4
{
public:
  Matrix5x4()
    : _11(1.0f), _12(0), _13(0), _14(0)
    , _21(0), _22(1.0f), _23(0), _24(0)
    , _31(0), _32(0), _33(1.0f), _34(0)
    , _41(0), _42(0), _43(0), _44(1.0f)
    , _51(0), _52(0), _53(0), _54(0)
  {}
  Matrix5x4(Float a11, Float a12, Float a13, Float a14,
         Float a21, Float a22, Float a23, Float a24,
         Float a31, Float a32, Float a33, Float a34,
         Float a41, Float a42, Float a43, Float a44,
         Float a51, Float a52, Float a53, Float a54)
    : _11(a11), _12(a12), _13(a13), _14(a14)
    , _21(a21), _22(a22), _23(a23), _24(a24)
    , _31(a31), _32(a32), _33(a33), _34(a34)
    , _41(a41), _42(a42), _43(a43), _44(a44)
    , _51(a51), _52(a52), _53(a53), _54(a54)
  {}

  bool operator==(const Matrix5x4 &o) const
  {
    return _11 == o._11 && _12 == o._12 && _13 == o._13 && _14 == o._14 &&
           _21 == o._21 && _22 == o._22 && _23 == o._23 && _24 == o._24 &&
           _31 == o._31 && _32 == o._32 && _33 == o._33 && _34 == o._34 &&
           _41 == o._41 && _42 == o._42 && _43 == o._43 && _44 == o._44 &&
           _51 == o._51 && _52 == o._52 && _53 == o._53 && _54 == o._54;
  }

  bool operator!=(const Matrix5x4 &aMatrix) const
  {
    return !(*this == aMatrix);
  }

  Matrix5x4 operator*(const Matrix5x4 &aMatrix) const
  {
    Matrix5x4 resultMatrix;

    resultMatrix._11 = this->_11 * aMatrix._11 + this->_12 * aMatrix._21 + this->_13 * aMatrix._31 + this->_14 * aMatrix._41;
    resultMatrix._12 = this->_11 * aMatrix._12 + this->_12 * aMatrix._22 + this->_13 * aMatrix._32 + this->_14 * aMatrix._42;
    resultMatrix._13 = this->_11 * aMatrix._13 + this->_12 * aMatrix._23 + this->_13 * aMatrix._33 + this->_14 * aMatrix._43;
    resultMatrix._14 = this->_11 * aMatrix._14 + this->_12 * aMatrix._24 + this->_13 * aMatrix._34 + this->_14 * aMatrix._44;
    resultMatrix._21 = this->_21 * aMatrix._11 + this->_22 * aMatrix._21 + this->_23 * aMatrix._31 + this->_24 * aMatrix._41;
    resultMatrix._22 = this->_21 * aMatrix._12 + this->_22 * aMatrix._22 + this->_23 * aMatrix._32 + this->_24 * aMatrix._42;
    resultMatrix._23 = this->_21 * aMatrix._13 + this->_22 * aMatrix._23 + this->_23 * aMatrix._33 + this->_24 * aMatrix._43;
    resultMatrix._24 = this->_21 * aMatrix._14 + this->_22 * aMatrix._24 + this->_23 * aMatrix._34 + this->_24 * aMatrix._44;
    resultMatrix._31 = this->_31 * aMatrix._11 + this->_32 * aMatrix._21 + this->_33 * aMatrix._31 + this->_34 * aMatrix._41;
    resultMatrix._32 = this->_31 * aMatrix._12 + this->_32 * aMatrix._22 + this->_33 * aMatrix._32 + this->_34 * aMatrix._42;
    resultMatrix._33 = this->_31 * aMatrix._13 + this->_32 * aMatrix._23 + this->_33 * aMatrix._33 + this->_34 * aMatrix._43;
    resultMatrix._34 = this->_31 * aMatrix._14 + this->_32 * aMatrix._24 + this->_33 * aMatrix._34 + this->_34 * aMatrix._44;
    resultMatrix._41 = this->_41 * aMatrix._11 + this->_42 * aMatrix._21 + this->_43 * aMatrix._31 + this->_44 * aMatrix._41;
    resultMatrix._42 = this->_41 * aMatrix._12 + this->_42 * aMatrix._22 + this->_43 * aMatrix._32 + this->_44 * aMatrix._42;
    resultMatrix._43 = this->_41 * aMatrix._13 + this->_42 * aMatrix._23 + this->_43 * aMatrix._33 + this->_44 * aMatrix._43;
    resultMatrix._44 = this->_41 * aMatrix._14 + this->_42 * aMatrix._24 + this->_43 * aMatrix._34 + this->_44 * aMatrix._44;
    resultMatrix._51 = this->_51 * aMatrix._11 + this->_52 * aMatrix._21 + this->_53 * aMatrix._31 + this->_54 * aMatrix._41 + aMatrix._51;
    resultMatrix._52 = this->_51 * aMatrix._12 + this->_52 * aMatrix._22 + this->_53 * aMatrix._32 + this->_54 * aMatrix._42 + aMatrix._52;
    resultMatrix._53 = this->_51 * aMatrix._13 + this->_52 * aMatrix._23 + this->_53 * aMatrix._33 + this->_54 * aMatrix._43 + aMatrix._53;
    resultMatrix._54 = this->_51 * aMatrix._14 + this->_52 * aMatrix._24 + this->_53 * aMatrix._34 + this->_54 * aMatrix._44 + aMatrix._54;

    return resultMatrix;
  }

  Matrix5x4& operator*=(const Matrix5x4 &aMatrix)
  {
    *this = *this * aMatrix;
    return *this;
  }

  union {
    struct {
      Float _11, _12, _13, _14;
      Float _21, _22, _23, _24;
      Float _31, _32, _33, _34;
      Float _41, _42, _43, _44;
      Float _51, _52, _53, _54;
    };
    Float components[20];
  };
};

} // namespace gfx
} // namespace mozilla

#endif /* MOZILLA_GFX_MATRIX_H_ */