edlgeom.h

Go to the documentation of this file.

/*
 * edlgeom.h
 *
 * Copyright (C) 2007-2026 Hybrid Software Helix Ltd. All rights reserved.
 */

 
#ifndef EDLGEOM_H
#define EDLGEOM_H
 
#include <algorithm>
#include <edl/edltypes.h>
#include <edl/edlmath.h>
#include <edl/edlvector.h>
#include <edl/edlerrors.h>
 
_BEGIN_EDL_NAMESPACE

template <typename PointType>
class PointTmpl
{
  public:
    PointTmpl() : x(0), y(0)
    {
    }
 
    PointTmpl(PointType ax, PointType ay) : x(ax), y(ay)
    {
    }
 
    PointTmpl(const PointTmpl<PointType>& p) = default;
 
    PointTmpl(PointTmpl<PointType>&& p) = default;
 
    PointTmpl<PointType>& operator=(const PointTmpl<PointType>& other) = default;
 
    PointTmpl<PointType>& operator += (const PointTmpl<PointType>& pt)
    {
        x += pt.x;
        y += pt.y;
        return *this;
    }
 
    PointTmpl<PointType> &operator -= (const PointTmpl<PointType> &pt)
    {
        x -= pt.x;
        y -= pt.y;
        return *this;
    }
 
    bool equal(const PointTmpl& point) const
    {
        return ((x == point.x) && (y == point.y));
    }
 
    bool operator == (const PointTmpl& point) const
    {
        return ((x == point.x) && (y == point.y));
    }
 
    bool operator!= (const PointTmpl& point) const
    {
        return !(*this == point);
    }

    PointTmpl<PointType> getMidPoint(const PointTmpl<PointType>& pt) const
    {
        return PointTmpl<PointType>((PointType)(x + (pt.x - x) / 2.0), (PointType)(y + (pt.y - y) / 2.0));
    }

    double getDistance(const PointTmpl<PointType>& pt) const
    {
        return hypot(pt.x - x, pt.y - y);
    }
 
    PointType x;
    PointType y;
};
 
template <typename PointType>
const PointTmpl<PointType> operator +(const PointTmpl<PointType>& lhp, const PointTmpl<PointType>& rhp)
{
    return PointTmpl<PointType>(lhp.x + rhp.x, lhp.y + rhp.y);
}
template <typename PointType>
const PointTmpl<PointType> operator -(const PointTmpl<PointType>& lhp, const PointTmpl<PointType>& rhp)
{
    return PointTmpl<PointType>(lhp.x - rhp.x, lhp.y - rhp.y);
}
 
typedef PointTmpl<double>               FPoint;
typedef PointTmpl<int32>                IntPoint;
typedef PointTmpl<uint32>               UIntPoint;
typedef PointTmpl<int64>                Int64Point;
typedef CEDLVector<FPoint>              CFPointVect;
 
// Point type may not be an unsigned type
template <typename PointType>
class RectTmpl
{
  public:
    RectTmpl() :
        x(0), y(0), dX(-1), dY(-1)
    {}
 
    RectTmpl(PointType ax, PointType ay, PointType adX, PointType adY) :
        x(ax), y(ay), dX(adX), dY(adY)
    {}
 
    RectTmpl(const RectTmpl<PointType>& r) = default;
 
    RectTmpl(RectTmpl<PointType>&& p) = default;
 
    RectTmpl<PointType>& operator=(const RectTmpl<PointType>& other) = default;
 
    PointType x;
    PointType y;
    PointType dX;
    PointType dY;
 
    void setEmpty()
    {
        dX = (PointType)-1;
        dY = (PointType)-1;
    }
 
    bool isEmpty() const
    {
        return (dX < (PointType)0 || dY < (PointType)0);
    }
 
    bool equal(const RectTmpl& rect) const
    {
        return ((x == rect.x) && (y == rect.y) &&
            (dX == rect.dX) && (dY == rect.dY));
    }
 
    bool operator== (const RectTmpl& rect) const
    {
        return ((x == rect.x) && (y == rect.y) &&
            (dX == rect.dX) && (dY == rect.dY));
    }
 
    bool operator!= (const RectTmpl& rect) const
    {
        return !(*this == rect);
    }
 
    bool similar(const RectTmpl& rect, float epsilon) const
    {
        return (fabs(x - rect.x) < epsilon &&
            fabs(y - rect.y) < epsilon &&
            fabs(dX - rect.dX) < epsilon &&
            fabs(dY - rect.dY) < epsilon);
    }
 
    RectTmpl& offset(PointType offX, PointType offY)
    {
        if (!isEmpty())
        {
            x += offX;
            y += offY;
        }
        return *this;
    }

    RectTmpl& inset(PointType insetX, PointType insetY)
    {
        if (!isEmpty())
        {
            if (dX < insetX * 2)
            {
                x += dX / 2;
                dX = 0;
            }
            else
            {
                dX -= insetX * 2;
                x += insetX;
            }
            if (dY < insetY * 2)
            {
                y += dY / 2;
                dY = 0;
            }
            else
            {
                dY -= insetY * 2;
                y += insetY;
            }
        }
        return *this;
    }

    bool expandToPoint(const PointTmpl<PointType>& point)
    {
        bool result = false;
        if (isEmpty())
        {
            x = point.x; y = point.y;
            dX = 0; dY = 0;
            result = true;
        }
        else
        {
            if (point.x < x)
            {
                dX += (x - point.x);
                x = point.x;
                result = true;
            }
            if (point.y < y)
            {
                dY += (y - point.y);
                y = point.y;
                result = true;
            }
            if (point.x > x + dX)
            {
                dX = point.x - x;
                result = true;
            }
            if (point.y > y + dY)
            {
                dY = point.y - y;
                result = true;
            }
        }
        return result;
    }

    void intersectRect(const RectTmpl<PointType>& rect)
    {
        if (!isEmpty() && !rect.isEmpty())
        {
            PointType llX = std::max<PointType>(x, rect.x);
            PointType urX = std::min<PointType>(x + dX, rect.x + rect.dX);
            PointType llY = std::max<PointType>(y, rect.y);
            PointType urY = std::min<PointType>(y + dY, rect.y + rect.dY);
            x = llX; y = llY; dX = urX - x; dY = urY - y;
        }
        else
        {
            setEmpty(); // Empty result
        }
    }

    void unionRect(const RectTmpl<PointType>& rect)
    {
        if (!isEmpty())
        {
            if (containsRect(rect))
            {
                // We completely contain the input rect; no need to grow
            }
            else if (!rect.isEmpty())
            {
                PointType llX = std::min<PointType>(x, rect.x);
                PointType urX = std::max<PointType>(x + dX, rect.x + rect.dX);
                PointType llY = std::min<PointType>(y, rect.y);
                PointType urY = std::max<PointType>(y + dY, rect.y + rect.dY);
                x = llX; y = llY; dX = urX - x; dY = urY - y;
            }
        }
        else
        {
            x = rect.x; y = rect.y; dX = rect.dX; dY = rect.dY;
        }
    }

    bool intersectsWithRect(const RectTmpl<PointType>& rect) const
    {
        RectTmpl<PointType> copy = rect;
        copy.intersectRect(*this);
        return !copy.isEmpty();
    }

    bool containsRect(const RectTmpl<PointType>& rect) const
    {
        if (rect.isEmpty() || isEmpty())
        {
            return false;
        }
        if (x <= rect.x &&                    // We start at or to the left of the input rect
            y <= rect.y &&                    // We start at or above the top of the input rect
            (x + dX) >= (rect.x + rect.dX) && // We extend to or beyond the rect to the right
            (y + dY) >= (rect.y + rect.dY))   // We extend to or beyond the rect to the bottom
        {
            return true;
        }
        else
        {
            return false;
        }
    }
 
    bool containsPoint(PointType px, PointType py) const
    {
        return x <= px && px <= getRight() && y <= py && py <= getBottom();
    }
 
    PointType getRight() const
    {
        return x + dX;
    }
 
    PointType getBottom() const
    {
        return y + dY;
    }
};
 
typedef RectTmpl<double>                FRect;
typedef RectTmpl<int32>                 IntRect;
typedef RectTmpl<int64>                 Int64Rect;

template <typename PointType>
class BoxTmpl
{
  public:
    BoxTmpl() : left(0), bottom(0), right(0), top(0)
    {}
 
    BoxTmpl(PointType _left, PointType _bottom, PointType _right, PointType _top) :
        left(_left), bottom(_bottom), right(_right), top(_top)
    {}
 
    BoxTmpl(const BoxTmpl<PointType>& b) :
        left(b.left), bottom(b.bottom), right(b.right), top(b.top)
    {}
 
    BoxTmpl(const RectTmpl<PointType>& rect)
    {
        if (rect.isEmpty())
        {
            left = rect.x;
            bottom = rect.y;
            right = rect.x;
            top = rect.y;
        }
        else
        {
            left = rect.x;
            bottom = rect.y;
            right = rect.x + rect.dX;
            top = rect.y + rect.dY;
        }
    }
 
    BoxTmpl(BoxTmpl<PointType>&& p) = default;
 
    BoxTmpl<PointType>& operator=(const BoxTmpl<PointType>& other) = default;
 
    RectTmpl<PointType> asRect() const
    {
        RectTmpl<PointType> result;
        result.x = left < right ? left : right;
        result.y = top < bottom ? top : bottom;
        result.dX = left < right ? (right - left) : (left - right);
        result.dY = top < bottom ? (bottom - top) : (top - bottom);
        return result;
    }
 
    void normalize()
    {
        if (left > right)
        {
            PointType tmp = left;
            left = right;
            right = tmp;
        }
        if (bottom > top)
        {
            PointType tmp = bottom;
            bottom = top;
            top = tmp;
        }
    }
 
    void offset(PointType x, PointType y)
    {
        left += x;
        right += x;
        top += y;
        bottom += y;
    }
 
    void scale(PointType s)
    {
        left *= s;
        right *= s;
        top *= s;
        bottom *= s;
    }
 
    bool equal(const BoxTmpl<PointType>& other) const
    {
        return (left == other.left)
            && (bottom == other.bottom)
            && (top == other.top)
            && (right == other.right);
    }
 
    PointType left;
    PointType bottom;
    PointType top;
    PointType right;
};
 
typedef BoxTmpl<double>                FBox;
typedef BoxTmpl<int32>                 IntBox;
typedef BoxTmpl<int64>                 Int64Box;

template <typename TItem>
class CTransformMatrix
{
  public:
    typedef enum {
        eDoesTranslate = 0x1, 
        eDoesScale     = 0x2, 
        eDoesRotate    = 0x4, 
        eIsComplex     = 0x8  
    } eOperationTypes;

    CTransformMatrix()
    {
        set();
    }

    CTransformMatrix(TItem _xx, TItem _xy, TItem _yx, TItem _yy, TItem _dx, TItem _dy)
    {
        set(_xx, _xy, _yx, _yy, _dx, _dy);
    }

    CTransformMatrix(const CTransformMatrix<TItem>& m)
    {
        *this = m;
    }

    CTransformMatrix& operator=(const CTransformMatrix<TItem>& m)
    {
        if (this != &m)
        {
            m_xx = m.m_xx;
            m_xy = m.m_xy;
            m_yx = m.m_yx;
            m_yy = m.m_yy;
            m_dx = m.m_dx;
            m_dy = m.m_dy;
        }
        return *this;
    }

    CTransformMatrix(const RectTmpl<TItem>& sourceRect, const RectTmpl<TItem>& destRect, bool allowZeroWidthRect = false)
    {
        set();
        if (sourceRect.dX > 0 && destRect.dX > 0 && sourceRect.dY > 0 && destRect.dY > 0)
        {
            m_xx = destRect.dX / sourceRect.dX;
            m_yy = destRect.dY / sourceRect.dY;
            m_dx = destRect.x - sourceRect.x * m_xx;
            m_dy = destRect.y - sourceRect.y * m_yy;
        }
        else if (allowZeroWidthRect &&
                 sourceRect.dX > 0 && destRect.dX > 0 && sourceRect.dY == 0 && destRect.dY == 0)
        {
            m_xx = destRect.dX / sourceRect.dX;
            m_yy = 1.0;
            m_dx = destRect.x - sourceRect.x * m_xx;
            m_dy = destRect.y - sourceRect.y * m_yy;
        }
        else if (allowZeroWidthRect &&
                 sourceRect.dX == 0 && destRect.dX == 0 && sourceRect.dY > 0 && destRect.dY > 0)
        {
            m_xx = 1.0;
            m_yy = destRect.dY / sourceRect.dY;
            m_dx = destRect.x - sourceRect.x * m_xx;
            m_dy = destRect.y - sourceRect.y * m_yy;
        }
    }

    void set(TItem _xx = 1, TItem _xy = 0, TItem _yx = 0, TItem _yy = 1, TItem _dx = 0, TItem _dy = 0)
    {
        m_xx = _xx; m_xy = _xy; m_yx = _yx; m_yy = _yy; m_dx = _dx; m_dy = _dy;
    }

    TItem xx() const { return m_xx; }

    TItem xy() const { return m_xy; }

    TItem yx() const { return m_yx; }

    TItem yy() const { return m_yy; }

    TItem dx() const { return m_dx; }

    TItem dy() const { return m_dy; }

    void setXX(TItem x) { m_xx = x; }

    void setXY(TItem x) { m_xy = x; }

    void setYX(TItem x) { m_yx = x; }

    void setYY(TItem x) { m_yy = x; }

    void setDX(TItem x) { m_dx = x; }

    void setDY(TItem x) { m_dy = x; }

    bool equal(const CTransformMatrix<TItem>& matrix, bool ignoreDXDY = false) const
    {
        return  (m_xx == matrix.xx()) &&
            (m_xy == matrix.xy()) &&
            (m_yx == matrix.yx()) &&
            (m_yy == matrix.yy()) &&
            (ignoreDXDY || (
                (m_dx == matrix.dx()) &&
                (m_dy == matrix.dy())));
    }

    bool identity(bool ignoreDXDY = false) const
    {
        return  (m_xx == 1.0) &&
            (m_xy == 0.0) &&
            (m_yx == 0.0) &&
            (m_yy == 1.0) &&
            (ignoreDXDY || (
                (m_dx == 0.0) &&
                (m_dy == 0.0)));
    }

    bool ortho(bool &rotated) const
    {
        rotated = false; // Until proven otherwise
        if (m_xy == 0 && m_yx == 0)
        {
            return true;
        }
        else if (m_xx == 0 && m_yy == 0)
        {
            rotated = true;
            return true;
        }
        else
        {
            return false;
        }
    }

    CTransformMatrix<TItem>& preMul(const CTransformMatrix<TItem>& matrix)
    {
        TItem a = matrix.xx() * m_xx + matrix.xy() * m_yx;
        TItem b = matrix.xx() * m_xy + matrix.xy() * m_yy;
        TItem c = matrix.yx() * m_xx + matrix.yy() * m_yx;
        TItem d = matrix.yx() * m_xy + matrix.yy() * m_yy;
        TItem e = matrix.dx() * m_xx + matrix.dy() * m_yx + m_dx;
        TItem f = matrix.dx() * m_xy + matrix.dy() * m_yy + m_dy;
        m_xx = a; m_xy = b; m_yx = c; m_yy = d; m_dx = e; m_dy = f;
        return *this;
    }

    CTransformMatrix<TItem>& postMul(const CTransformMatrix<TItem>& matrix) {
        TItem a = m_xx * matrix.xx() + m_xy * matrix.yx();
        TItem b = m_xx * matrix.xy() + m_xy * matrix.yy();
        TItem c = m_yx * matrix.xx() + m_yy * matrix.yx();
        TItem d = m_yx * matrix.xy() + m_yy * matrix.yy();
        TItem e = m_dx * matrix.xx() + m_dy * matrix.yx() + matrix.dx();
        TItem f = m_dx * matrix.xy() + m_dy * matrix.yy() + matrix.dy();
        m_xx = a; m_xy = b; m_yx = c; m_yy = d; m_dx = e; m_dy = f;
        return *this;
    }

    bool degenerate () const
    {
        // Find the determinant
        TItem det = m_xx * m_yy - m_yx * m_xy;
 
        // A zero determinant indicates a degenerate matrix.
        return (det == 0.0);
    }

    bool invert()
    {
        TItem det;
 
        // Find the determinant
        det = m_xx * m_yy - m_yx * m_xy;
        if (det == 0.0)
        {
            // This matrix cannot be inverted
            return false;
        }
 
        TItem a = m_yy / det;
        TItem b = -m_xy / det;
        TItem c = -m_yx / det;
        TItem d = m_xx / det;
        TItem e = -(m_dx * m_yy - m_dy * m_yx) / det;
        TItem f = (m_dx * m_xy - m_dy * m_xx) / det;
        m_xx = a; m_xy = b; m_yx = c; m_yy = d; m_dx = e; m_dy = f;
        return true;
    }


    void transform(PointTmpl<TItem>& result, const PointTmpl<TItem>& point, bool ignoreDXDY = false) const
    {
        TItem x = point.x * m_xx + point.y * m_yx;
        TItem y = point.x * m_xy + point.y * m_yy;
        if (!ignoreDXDY)
        {
            x += m_dx;
            y += m_dy;
        }
        result.x = x;
        result.y = y;
    }

    PointTmpl<TItem> transform(const PointTmpl<TItem>& point, bool ignoreDXDY = false) const
    {
        PointTmpl<TItem> result;
 
        TItem x = point.x * m_xx + point.y * m_yx;
        TItem y = point.x * m_xy + point.y * m_yy;
        if (!ignoreDXDY)
        {
            x += m_dx;
            y += m_dy;
        }
        result.x = x;
        result.y = y;
 
        return result;
    }

    bool iTransform(PointTmpl<TItem>& result, const PointTmpl<TItem>& point, bool ignoreDXDY = false) const
    {
        TItem det = m_xx * m_yy - m_yx * m_xy;
        if (det == 0.0)
        {
            // This matrix cannot be inverted
            return false;
        }
        TItem iDet = 1 / det;
 
        TItem x = point.x;
        TItem y = point.y;
        if (!ignoreDXDY)
        {
            x -= m_dx;
            y -= m_dy;
        }
 
        result.x = ((x * m_yy - y * m_yx) * iDet);
        result.y = ((-x * m_xy + y * m_xx) * iDet);
 
        return true;
    }

    std::pair<bool, PointTmpl<TItem> > iTransform(const PointTmpl<TItem>& point, bool ignoreDXDY = false) const
    {
        PointTmpl<TItem> result;
 
        TItem det = m_xx * m_yy - m_yx * m_xy;
        if (det == 0.0)
        {
            // This matrix cannot be inverted
            return std::pair<bool, PointTmpl<TItem> >(false, result);
        }
        TItem iDet = 1 / det;
 
        TItem x = point.x;
        TItem y = point.y;
        if (!ignoreDXDY)
        {
            x -= m_dx;
            y -= m_dy;
        }
 
        result.x = ((x * m_yy - y * m_yx) * iDet);
        result.y = ((-x * m_xy + y * m_xx) * iDet);
 
        return std::pair<bool, PointTmpl<TItem> >(true, result);
    }

    void rotate(double radians)
    {
        constexpr double rad_90 = PI / 2.0;
        constexpr double rad_180 = PI;
        constexpr double rad_270 = 3.0 * PI / 2.0;
        constexpr double rad_360 = 2.0 * PI;
        constexpr double eps = 1.0e-6;
 
        // Clamp to between 0 and 2*PI
        while (radians >= rad_360)
            radians -= rad_360;
        while (radians < 0.0)
            radians += rad_360;
 
        CTransformMatrix rotate;
 
        if (fabs(radians) < 1.0e-6);                                    // 0 degrees
        else if (fabs(radians - rad_90) < eps)                          // 90 degrees
            rotate = CTransformMatrix(0.0, 1.0, -1.0, 0.0, 0.0, 0.0);
        else if (fabs(radians - rad_180) < eps)                         // 180 degrees
            rotate = CTransformMatrix(-1.0, 0.0, 0.0, -1.0, 0.0, 0.0);
        else if (fabs(radians - rad_270) < eps)                         // 270 degrees
            rotate = CTransformMatrix(0.0, -1.0, 1.0, 0.0, 0.0, 0.0);
        else
            rotate = CTransformMatrix(cos(radians), sin(radians), -sin(radians), cos(radians), 0.0, 0.0);
 
        preMul(rotate);
    }

    void scale(TItem xscale, TItem yscale)
    {
        CTransformMatrix mat(xscale, 0.0, 0.0, yscale, 0.0, 0.0);
        preMul(mat);
    }

    void translate(TItem dx, TItem dy)
    {
        CTransformMatrix mat(1.0, 0.0, 0.0, 1.0, dx, dy);
        preMul(mat);
    }

    TItem determinant() const
    {
        return m_xx * m_yy - m_yx * m_xy;
    }

    void transformRect(RectTmpl<TItem>& rect, bool ignoreDXDY = false) const
    {
        TItem llX, llY, urX, urY;
        PointTmpl<TItem> point, transformedPoint;
 
        if (rect.isEmpty())
            return;
#define EDLTMIN(a, b) (a) < (b) ? (a) : (b)
#define EDLTMAX(a, b) (a) < (b) ? (b) : (a)
        // Transform each of the points that make up the corners,
        // and find the extremes.
        point.x = rect.x;
        point.y = rect.y;
        transform(transformedPoint, point, ignoreDXDY);
        llX = urX = transformedPoint.x;
        llY = urY = transformedPoint.y;
        point.x = rect.x + rect.dX;
        transform(transformedPoint, point, ignoreDXDY);
        llX = EDLTMIN(llX, transformedPoint.x);
        llY = EDLTMIN(llY, transformedPoint.y);
        urX = EDLTMAX(urX, transformedPoint.x);
        urY = EDLTMAX(urY, transformedPoint.y);
        point.y = rect.y + rect.dY;
        transform(transformedPoint, point, ignoreDXDY);
        llX = EDLTMIN(llX, transformedPoint.x);
        llY = EDLTMIN(llY, transformedPoint.y);
        urX = EDLTMAX(urX, transformedPoint.x);
        urY = EDLTMAX(urY, transformedPoint.y);
        point.x = rect.x;
        transform(transformedPoint, point, ignoreDXDY);
        llX = EDLTMIN(llX, transformedPoint.x);
        llY = EDLTMIN(llY, transformedPoint.y);
        urX = EDLTMAX(urX, transformedPoint.x);
        urY = EDLTMAX(urY, transformedPoint.y);
#undef EDLTMIN
#undef EDLTMAX
        rect.x = llX; rect.y = llY; rect.dX = urX - llX; rect.dY = urY - llY;
    }

    uint32 classify() const
    {
        uint32 operationFlags = 0;
 
        if ((m_dx != 0.0) || (m_dy != 0))
            operationFlags |= eDoesTranslate;
 
        if ((m_xy == 0) && (m_yx == 0))
            operationFlags |= eDoesScale;
 
        if ((m_xx == m_yy) && (m_xy == -m_yx)) {
            operationFlags |= eDoesRotate;
 
            double det = m_xx * m_yy - m_yx * m_xy;
            double eps = 1.0e-06;
 
            if (fabs(det - 1.0) > eps)
                operationFlags |= eDoesScale;
 
        }
        else {
 
            operationFlags |= eIsComplex;
 
        }
 
        return operationFlags;
    }
 

    void decompose(PointTmpl<TItem>& translate, FPoint& scale, FPoint& shear, double& rotationAngle, double eps = 1.0e-06) const
    {
        // unit square decomposition of matrix ...
        // A={a,c}, B={b,d}, Q=angle between vectors (ie PI/2-shear angle)!
        // |A|=sqrt(aa+cc)
        // |B|=sqrt(bb+dd)
        // AxB=|A||B|sinQ=ad-bc=(new area of unit square)
        // A.B=|A||B|cosQ=ab+cd
 
        // ScaleX   = |A|
        // ScaleY   = |B|
        // Scale    = ScaleY*cos(PI/2-Q) = ScaleY*sinQ = AxB/|A|
        // AbsScale = abs(Scale) // effectively remove any mirror
        // Aspect   = ScaleX/abs(Scale)
        // Shear    = PI/2-Q = PI/2-acos((A.B)/(|A||B|))
        // Rotate   = atan2(c,a)
 
        double a = m_xx;
        double b = m_xy;
        double c = m_yx;
        double d = m_yy;
 
        // Determine the cross product (determinant), modulus (length) of A and scale
        double AdotB = a * b + c * d;
        double ModA = ::sqrt(a * a + c * c);
        double ModB = ::sqrt(b * b + d * d);
        double pi = PI;
        double Shear = (pi / 2 - ::acos(AdotB / (ModB * ModA)));
 
        // Set the translate
        translate.x = m_dx;
        translate.y = m_dy;
 
        // Set the rotation angle
        if ((b == 0) || (c == 0))
            rotationAngle = 0;
        else {
            rotationAngle = ::atan2(c, a);
            if (fabs(rotationAngle) < eps)
                rotationAngle = 0.0;
        }
 
        //Set the scale
        scale.x = ((a != 0) && (a > 0)) || ((a == 0) && (c * b <= 0)) ? ModA : -ModA;
        scale.y = ((d != 0) && (d > 0)) || ((d == 0) && (c * b <= 0)) ? ModB : -ModB;
 
        // Set shear
        if (fabs(Shear) < eps) {
 
            // Case: no shear
            shear.x = 0;
            shear.y = 0;
 
        }
        else {
 
            // Case: shear present.
            shear.x = -tan(Shear);
            shear.y = 0;
 
            if (a == d) {
 
                // Shear will handle differential x and y scaling
                scale.x = a;
                scale.y = d;
 
            }
        }
 
    }
 
    typedef struct
    {
        PointTmpl<TItem>    translate;
        FPoint              scale;
        FPoint              shear;
        double              rotationAngle;
    } DecomposeInfo;

    DecomposeInfo decompose(double eps = 1.0e-06) const
    {
        DecomposeInfo info;
 
        decompose(info.translate, info.scale, info.shear, info.rotationAngle, eps);
 
        return info;
    }

    double getScale() const
    {
        double a = m_xx;
        double b = m_xy;
        double c = m_yx;
        double d = m_yy;
 
        double AxB = a * d - b * c;
        double ModA = ::sqrt(a * a + c * c);
 
        return ::fabs(AxB / ModA);
    }

    bool isSimpleScaled(TItem& scale) const
    {
        bool isSimple = ((m_xx == m_yy) && (m_xy == 0) && (m_yx == 0));
 
        if (isSimple)
            scale = m_xx;
 
        return isSimple;
    }

    std::pair<bool, TItem> isSimpleScaled() const
    {
        TItem   scale;
        bool    val = isSimpleScaled(scale);
 
        return std::pair<bool, TItem>(val, scale);
    }

    void decompose(PointTmpl<TItem>& translate, double& rotationAngle, double& shearAngle, FPoint& scale) const
    {
        CTransformMatrix<TItem> r = *this;
        translate = PointTmpl<TItem>(r.dx(), r.dy());
        r.setDX(0);
        r.setDY(0);
        {
            PointTmpl<TItem> point(1, 0);
            r.transform(point, point);
            rotationAngle = atan2(point.y, point.x);
            {
                CTransformMatrix<TItem> rotate;
                rotate.rotate(-rotationAngle);
                r.postMul(rotate);
            }
            r.setXY(0);
        }
        {
            PointTmpl<TItem> point(0, 1);
            r.transform(point, point);
            shearAngle = (PI / 2.0) - atan2(point.y, point.x);
            if (shearAngle < -(PI / 2.0))
            {
                shearAngle += PI;
            }
            else if (shearAngle > (PI / 2.0))
            {
                shearAngle -= PI;
            }
            r.setYX(0);
        }
        scale = FPoint(r.xx(), r.yy());
    }

    void compose(const PointTmpl<TItem>& translateAmount, double rotationAngle, double shearAngle, const FPoint& scaleAmount)
    {
        // Reset
        set();
 
        // First translate
        translate(translateAmount.x, translateAmount.y);
 
        // Then rotate
        rotate(rotationAngle);
 
        // Then the shear. Must be less than +/- PI/2.
        if (shearAngle < -(PI / 2.0) || shearAngle >(PI / 2.0))
        {
            throwEDLError(EDL_ERR_BAD_ARGUMENTS);
        }
        double shearAmount = tan(shearAngle);
        CTransformMatrix<TItem> shearMat(1, 0, (TItem)shearAmount, 1, 0, 0);
        preMul(shearMat);
 
        // Finally scale
        scale((TItem)scaleAmount.x, (TItem)scaleAmount.y);
    }

    template <typename AType>
    void asArray(AType* array) const
    {
        if (!array)
        {
            throwEDLError(EDL_ERR_BAD_ARGUMENTS);
        }
        array[0] = (AType)m_xx;
        array[1] = (AType)m_xy;
        array[2] = (AType)m_yx;
        array[3] = (AType)m_yy;
        array[4] = (AType)m_dx;
        array[5] = (AType)m_dy;
    }
 
  private:
    TItem m_xx, m_xy, m_yx, m_yy, m_dx, m_dy;
};
 
typedef CTransformMatrix<double> FMatrix;
 
_END_EDL_NAMESPACE
 
#endif /* EDLGEOM_H */