Add missing 2geom files

(bzr r10027)

3 files changed, 850 insertions, 0 deletions

diff --git a/src/2geom/affine.cpp b/src/2geom/affine.cpp
new file mode 100644
index 000000000..925f43820
--- /dev/null
+++ b/src/2geom/affine.cpp
@@ -0,0 +1,489 @@
+#define __Geom_MATRIX_C__
+
+/** \file
+ * Various matrix routines.  Currently includes some Geom::Rotate etc. routines too.
+ */
+
+/*
+ * Authors:
+ *   Lauris Kaplinski <lauris@kaplinski.com>
+ *   Michael G. Sloan <mgsloan@gmail.com>
+ *
+ * This code is in public domain
+ */
+
+#include <2geom/utils.h>
+#include <2geom/affine.h>
+#include <2geom/point.h>
+
+namespace Geom {
+
+/** Creates a Affine given an axis and origin point.
+ *  The axis is represented as two vectors, which represent skew, rotation, and scaling in two dimensions.
+ *  from_basis(Point(1, 0), Point(0, 1), Point(0, 0)) would return the identity matrix.
+
+ \param x_basis the vector for the x-axis.
+ \param y_basis the vector for the y-axis.
+ \param offset the translation applied by the matrix.
+ \return The new Affine.
+ */
+//NOTE: Inkscape's version is broken, so when including this version, you'll have to search for code with this func
+Affine from_basis(Point const x_basis, Point const y_basis, Point const offset) {
+    return Affine(x_basis[X], x_basis[Y],
+                  y_basis[X], y_basis[Y],
+                  offset [X], offset [Y]);
+}
+
+Point Affine::xAxis() const {
+    return Point(_c[0], _c[1]);
+}
+
+Point Affine::yAxis() const {
+    return Point(_c[2], _c[3]);
+}
+
+/** Gets the translation imparted by the Affine.
+ */
+Point Affine::translation() const {
+    return Point(_c[4], _c[5]);
+}
+
+void Affine::setXAxis(Point const &vec) {
+    for(int i = 0; i < 2; i++)
+        _c[i] = vec[i];
+}
+
+void Affine::setYAxis(Point const &vec) {
+    for(int i = 0; i < 2; i++)
+        _c[i + 2] = vec[i];
+}
+
+/** Sets the translation imparted by the Affine.
+ */
+void Affine::setTranslation(Point const &loc) {
+    for(int i = 0; i < 2; i++)
+        _c[i + 4] = loc[i];
+}
+
+/** Calculates the amount of x-scaling imparted by the Affine.  This is the scaling applied to
+ *  the original x-axis region.  It is \emph{not} the overall x-scaling of the transformation.
+ *  Equivalent to L2(m.xAxis())
+ */
+double Affine::expansionX() const {
+    return sqrt(_c[0] * _c[0] + _c[1] * _c[1]);
+}
+
+/** Calculates the amount of y-scaling imparted by the Affine.  This is the scaling applied before
+ *  the other transformations.  It is \emph{not} the overall y-scaling of the transformation. 
+ *  Equivalent to L2(m.yAxis())
+ */
+double Affine::expansionY() const {
+    return sqrt(_c[2] * _c[2] + _c[3] * _c[3]);
+}
+
+void Affine::setExpansionX(double val) {
+    double exp_x = expansionX();
+    if(!are_near(exp_x, 0.0)) {  //TODO: best way to deal with it is to skip op?
+        double coef = val / expansionX();
+        for(unsigned i=0;i<2;i++) _c[i] *= coef;
+    }
+}
+
+void Affine::setExpansionY(double val) {
+    double exp_y = expansionY();
+    if(!are_near(exp_y, 0.0)) {  //TODO: best way to deal with it is to skip op?
+        double coef = val / expansionY();
+        for(unsigned i=2; i<4; i++) _c[i] *= coef;
+    }
+}
+
+/** Sets this matrix to be the Identity Affine. */
+void Affine::setIdentity() {
+    _c[0] = 1.0; _c[1] = 0.0;
+    _c[2] = 0.0; _c[3] = 1.0;
+    _c[4] = 0.0; _c[5] = 0.0;
+}
+
+/** @brief Check whether this matrix is an identity matrix.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           1 & 0 & 0 \\
+           0 & 1 & 0 \\
+           0 & 0 & 1 \end{array}\right]\f$ */
+bool Affine::isIdentity(Coord eps) const {
+    return are_near(_c[0], 1.0, eps) && are_near(_c[1], 0.0, eps) &&
+           are_near(_c[2], 0.0, eps) && are_near(_c[3], 1.0, eps) &&
+           are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps);
+}
+
+/** @brief Check whether this matrix represents a pure translation.
+ * Will return true for the identity matrix, which represents a zero translation.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           1 & 0 & 0 \\
+           0 & 1 & 0 \\
+           a & b & 1 \end{array}\right]\f$ */
+bool Affine::isTranslation(Coord eps) const {
+    return are_near(_c[0], 1.0, eps) && are_near(_c[1], 0.0, eps) &&
+           are_near(_c[2], 0.0, eps) && are_near(_c[3], 1.0, eps);
+}
+/** @brief Check whether this matrix represents a pure nonzero translation.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           1 & 0 & 0 \\
+           0 & 1 & 0 \\
+           a & b & 1 \end{array}\right]\f$ and \f$a, b \neq 0\f$ */
+bool Affine::isNonzeroTranslation(Coord eps) const {
+    return are_near(_c[0], 1.0, eps) && are_near(_c[1], 0.0, eps) &&
+           are_near(_c[2], 0.0, eps) && are_near(_c[3], 1.0, eps) &&
+           (!are_near(_c[4], 0.0, eps) || !are_near(_c[5], 0.0, eps));
+}
+
+/** @brief Check whether this matrix represents pure scaling.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           a & 0 & 0 \\
+           0 & b & 0 \\
+           0 & 0 & 1 \end{array}\right]\f$. */
+bool Affine::isScale(Coord eps) const {
+    return are_near(_c[1], 0.0, eps) && are_near(_c[2], 0.0, eps) && 
+           are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps);
+}
+
+/** @brief Check whether this matrix represents pure, nonzero scaling.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           a & 0 & 0 \\
+           0 & b & 0 \\
+           0 & 0 & 1 \end{array}\right]\f$ and \f$a, b \neq 1\f$. */
+bool Affine::isNonzeroScale(Coord eps) const {
+    return (!are_near(_c[0], 1.0, eps) || !are_near(_c[3], 1.0, eps)) &&  //NOTE: these are the diags, and the next line opposite diags
+           are_near(_c[1], 0.0, eps) && are_near(_c[2], 0.0, eps) && 
+           are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps);
+}
+
+/** @brief Check whether this matrix represents pure uniform scaling.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           a & 0 & 0 \\
+           0 & a & 0 \\
+           0 & 0 & 1 \end{array}\right]\f$. */
+bool Affine::isUniformScale(Coord eps) const {
+    return are_near(_c[0], _c[3], eps) &&
+           are_near(_c[1], 0.0, eps) && are_near(_c[2], 0.0, eps) &&  
+           are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps);
+}
+
+/** @brief Check whether this matrix represents pure, nonzero uniform scaling.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           a & 0 & 0 \\
+           0 & a & 0 \\
+           0 & 0 & 1 \end{array}\right]\f$ and \f$a \neq 1\f$. */
+bool Affine::isNonzeroUniformScale(Coord eps) const {
+    return !are_near(_c[0], 1.0, eps) && are_near(_c[0], _c[3], eps) &&
+           are_near(_c[1], 0.0, eps) && are_near(_c[2], 0.0, eps) &&  
+           are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps);
+}
+
+/** @brief Check whether this matrix represents a pure rotation.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           a & b & 0 \\
+           -b & a & 0 \\
+           0 & 0 & 1 \end{array}\right]\f$ and \f$a^2 + b^2 = 1\f$. */
+bool Affine::isRotation(Coord eps) const {
+    return are_near(_c[0], _c[3], eps) && are_near(_c[1], -_c[2], eps) &&
+           are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps) &&
+           are_near(_c[0]*_c[0] + _c[1]*_c[1], 1.0, eps);
+}
+
+/** @brief Check whether this matrix represents a pure, nonzero rotation.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           a & b & 0 \\
+           -b & a & 0 \\
+           0 & 0 & 1 \end{array}\right]\f$, \f$a^2 + b^2 = 1\f$ and \f$a \neq 1\f$. */
+bool Affine::isNonzeroRotation(Coord eps) const {
+    return !are_near(_c[0], 1.0, eps) &&
+           are_near(_c[0], _c[3], eps) && are_near(_c[1], -_c[2], eps) &&
+           are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps) &&
+           are_near(_c[0]*_c[0] + _c[1]*_c[1], 1.0, eps);
+}
+
+/** @brief Check whether this matrix represents pure horizontal shearing.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           1 & 0 & 0 \\
+           k & 1 & 0 \\
+           0 & 0 & 1 \end{array}\right]\f$. */
+bool Affine::isHShear(Coord eps) const {
+    return are_near(_c[0], 1.0, eps) && are_near(_c[1], 0.0, eps) &&
+           are_near(_c[3], 1.0, eps) && are_near(_c[4], 0.0, eps) &&
+           are_near(_c[5], 0.0, eps);
+}
+/** @brief Check whether this matrix represents pure, nonzero horizontal shearing.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           1 & 0 & 0 \\
+           k & 1 & 0 \\
+           0 & 0 & 1 \end{array}\right]\f$ and \f$k \neq 0\f$. */
+bool Affine::isNonzeroHShear(Coord eps) const {
+    return are_near(_c[0], 1.0, eps) && are_near(_c[1], 0.0, eps) &&
+          !are_near(_c[2], 0.0, eps) && are_near(_c[3], 1.0, eps) &&
+           are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps);
+}
+
+/** @brief Check whether this matrix represents pure vertical shearing.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           1 & k & 0 \\
+           0 & 1 & 0 \\
+           0 & 0 & 1 \end{array}\right]\f$. */
+bool Affine::isVShear(Coord eps) const {
+    return are_near(_c[0], 1.0, eps) && are_near(_c[2], 0.0, eps) &&
+           are_near(_c[3], 1.0, eps) && are_near(_c[4], 0.0, eps) &&
+           are_near(_c[5], 0.0, eps);
+}
+
+/** @brief Check whether this matrix represents pure, nonzero vertical shearing.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           1 & k & 0 \\
+           0 & 1 & 0 \\
+           0 & 0 & 1 \end{array}\right]\f$ and \f$k \neq 0\f$. */
+bool Affine::isNonzeroVShear(Coord eps) const {
+    return are_near(_c[0], 1.0, eps) && !are_near(_c[1], 0.0, eps) &&
+           are_near(_c[2], 0.0, eps) && are_near(_c[3], 1.0, eps) &&
+           are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps);
+}
+
+/** @brief Check whether this matrix represents zooming.
+ * Zooming is any combination of translation and uniform scaling. It preserves angles, ratios
+ * of distances between arbitrary points and unit vectors of line segments.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           a & 0 & 0 \\
+           0 & a & 0 \\
+           b & c & 1 \end{array}\right]\f$. */
+bool Affine::isZoom(Coord eps) const {
+    return are_near(_c[0], _c[3], eps) && are_near(_c[1], 0, eps) && are_near(_c[2], 0, eps);
+}
+
+/** @brief Check whether the transformation preserves areas of polygons.
+ * This means that the transformation can be any combination of translation, rotation,
+ * shearing and squeezing (non-uniform scaling such that the absolute value of the product
+ * of Y-scale and X-scale is 1).
+ * @param eps Numerical tolerance
+ * @return True iff \f$|\det A| = 1\f$. */
+bool Affine::preservesArea(Coord eps) const
+{
+    return are_near(descrim2(), 1.0, eps);
+}
+
+/** @brief Check whether the transformation preserves angles between lines.
+ * This means that the transformation can be any combination of translation, uniform scaling
+ * and rotation.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           a & b & 0 \\
+           -b & a & 0 \\
+           c & d & 1 \end{array}\right]\f$. */
+bool Affine::preservesAngles(Coord eps) const
+{
+    return are_near(_c[0], _c[3], eps) && are_near(_c[1], -_c[2], eps);
+}
+
+/** @brief Check whether the transformation preserves distances between points.
+ * This means that the transformation can be any combination of translation and rotation.
+ * @param eps Numerical tolerance
+ * @return True iff the matrix is of the form
+ *         \f$\left[\begin{array}{ccc}
+           a & b & 0 \\
+           -b & a & 0 \\
+           c & d & 1 \end{array}\right]\f$ and \f$a^2 + b^2 = 1\f$. */
+bool Affine::preservesDistances(Coord eps) const
+{
+    return are_near(_c[0], _c[3], eps) && are_near(_c[1], -_c[2], eps) &&
+           are_near(_c[0] * _c[0] + _c[1] * _c[1], 1.0, eps);
+}
+
+/** @brief Check whether this transformation flips objects.
+ * A transformation flips objects if it has a negative scaling component. */
+bool Affine::flips() const {
+    // TODO shouldn't this be det() < 0?
+    return cross(xAxis(), yAxis()) > 0;
+}
+
+/** @brief Check whether this matrix is singular.
+ * Singular matrices have no inverse, which means that applying them to a set of points
+ * results in a loss of information.
+ * @param eps Numerical tolerance
+ * @return True iff the determinant is near zero. */
+bool Affine::isSingular(Coord eps) const {
+    return are_near(det(), 0.0, eps);
+}
+
+/** @brief Compute the inverse matrix.
+ * Inverse is a matrix (denoted \f$A^{-1}) such that \f$AA^{-1} = A^{-1}A = I\f$.
+ * Singular matrices have no inverse (for example a matrix that has two of its columns equal).
+ * For such matrices, the identity matrix will be returned instead.
+ * @param eps Numerical tolerance
+ * @return Inverse of the matrix, or the identity matrix if the inverse is undefined.
+ * @post (m * m.inverse()).isIdentity() == true */
+Affine Affine::inverse() const {
+    Affine d;
+    
+    double mx = std::max(fabs(_c[0]) + fabs(_c[1]), 
+                         fabs(_c[2]) + fabs(_c[3])); // a random matrix norm (either l1 or linfty
+    if(mx > 0) {
+        Geom::Coord const determ = det();
+        if (!rel_error_bound(determ, mx*mx)) {
+            Geom::Coord const ideterm = 1.0 / (determ);
+            
+            d._c[0] =  _c[3] * ideterm;
+            d._c[1] = -_c[1] * ideterm;
+            d._c[2] = -_c[2] * ideterm;
+            d._c[3] =  _c[0] * ideterm;
+            d._c[4] = (-_c[4] * d._c[0] - _c[5] * d._c[2]);
+            d._c[5] = (-_c[4] * d._c[1] - _c[5] * d._c[3]);
+        } else {
+            d.setIdentity();
+        }
+    } else {
+        d.setIdentity();
+    }
+
+    return d;
+}
+
+/** @brief Calculate the determinant.
+ * @return \f$\det A\f$. */
+Coord Affine::det() const {
+    // TODO this can overflow
+    return _c[0] * _c[3] - _c[1] * _c[2];
+}
+
+/** @brief Calculate the square of the descriminant.
+ * This is simply the absolute value of the determinant.
+ * @return \f$|\det A|\f$. */
+Coord Affine::descrim2() const {
+    return fabs(det());
+}
+
+/** @brief Calculate the descriminant.
+ * If the matrix doesn't contain a non-uniform scaling or shearing component, this value says
+ * how will the length any line segment change after applying this transformation
+ * to arbitrary objects on a plane (the new length will be
+ * @code line_seg.length() * m.descrim()) @endcode.
+ * @return \f$\sqrt{|\det A|}\f$. */
+Coord Affine::descrim() const {
+    return sqrt(descrim2());
+}
+
+/** @brief Combine this transformation with another one.
+ * After this operation, the matrix will correspond to the transformation
+ * obtained by first applying the original version of this matrix, and then
+ * applying @a m. */
+Affine &Affine::operator*=(Affine const &o) {
+    Coord nc[6];
+    for(int a = 0; a < 5; a += 2) {
+        for(int b = 0; b < 2; b++) {
+            nc[a + b] = _c[a] * o._c[b] + _c[a + 1] * o._c[b + 2];
+        }
+    }
+    for(int a = 0; a < 6; ++a) {
+        _c[a] = nc[a];
+    }
+    _c[4] += o._c[4];
+    _c[5] += o._c[5];
+    return *this;
+}
+
+//TODO: What's this!?!
+Affine elliptic_quadratic_form(Affine const &m) {
+    double od = m[0] * m[1]  +  m[2] * m[3];
+    Affine ret (m[0]*m[0] + m[1]*m[1], od,
+                od, m[2]*m[2] + m[3]*m[3],
+                0, 0);
+    return ret; // allow NRVO
+}
+
+Eigen::Eigen(Affine const &m) {
+    double const B = -m[0] - m[3];
+    double const C = m[0]*m[3] - m[1]*m[2];
+    double const center = -B/2.0;
+    double const delta = sqrt(B*B-4*C)/2.0;
+    values[0] = center + delta; values[1] = center - delta;
+    for (int i = 0; i < 2; i++) {
+        vectors[i] = unit_vector(rot90(Point(m[0]-values[i], m[1])));
+    }
+}
+
+static void quadratic_roots(double q0, double q1, double q2, int &n, double&r0, double&r1) {
+    std::vector<double> r;
+    if(q2 == 0) {
+        if(q1 == 0) { // zero or infinite roots
+            n = 0;
+        } else {
+            n = 1;
+            r0 = -q0/q1;
+        }
+    } else {
+        double desc = q1*q1 - 4*q2*q0;
+        if (desc < 0)
+            n = 0;
+        else if (desc == 0) {
+            n = 1;
+            r0 = -q1/(2*q2);
+        } else {
+            n = 2;
+            desc = std::sqrt(desc);
+            double t = -0.5*(q1+sgn(q1)*desc);
+            r0 = t/q2;
+            r1 = q0/t;
+        }
+    }
+}
+
+Eigen::Eigen(double m[2][2]) {
+    double const B = -m[0][0] - m[1][1];
+    double const C = m[0][0]*m[1][1] - m[1][0]*m[0][1];
+    //double const desc = B*B-4*C;
+    //double t = -0.5*(B+sgn(B)*desc);
+    int n;
+    values[0] = values[1] = 0;
+    quadratic_roots(C, B, 1, n, values[0], values[1]);
+    for (int i = 0; i < n; i++)
+        vectors[i] = unit_vector(rot90(Point(m[0][0]-values[i], m[0][1])));
+    for (int i = n; i < 2; i++) 
+        vectors[i] = Point(0,0);
+}
+
+}  //namespace Geom
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/affine.h b/src/2geom/affine.h
new file mode 100644
index 000000000..9a6352f8b
--- /dev/null
+++ b/src/2geom/affine.h
@@ -0,0 +1,253 @@
+/** \file
+ *  \brief 3x3 affine transformation matrix.
+ *//*
+ * Main authors:
+ *   Lauris Kaplinski <lauris@kaplinski.com> (Original NRAffine definition and related macros)
+ *   Nathan Hurst <njh@mail.csse.monash.edu.au> (Geom::Affine class version of the above)
+ *   Michael G. Sloan <mgsloan@gmail.com> (reorganization and additions)
+ *   Krzysztof Kosiński <tweenk.pl@gmail.com> (removal of boilerplate, docs)
+ *
+ * This code is in public domain.
+ */
+
+#ifndef SEEN_LIB2GEOM_MATRIX_H
+#define SEEN_LIB2GEOM_MATRIX_H
+
+#include <boost/operators.hpp>
+#include <2geom/forward.h>
+#include <2geom/point.h>
+
+namespace Geom {
+
+/**
+ * @brief 3x3 matrix representing an affine transformation.
+ *
+ * Affine transformations on elements of a vector space are transformations which can be
+ * expressed in terms of matrix multiplication followed by addition
+ * (\f$x \mapsto A x + b\f$). They can be thought of as generalizations of linear functions
+ * (\f$y = a x + b\f$) to vector spaces. Affine transformations of points on a 2D plane preserve
+ * the following properties:
+ *
+ * - Colinearity: if three points lie on the same line, they will still be colinear after
+ *   an affine transformation.
+ * - Ratios of distances between points on the same line are preserved
+ * - Parallel lines remain parallel.
+ *
+ * All affine transformations on 2D points can be written as combinations of scaling, rotation,
+ * shearing and translation. They can be represented as a combination of a vector and a 2x2 matrix,
+ * but this form is inconvenient to work with. A better solution is to represent all affine
+ * transformations on the 2D plane as 3x3 matrices, where the last column has fixed values.
+ * \f[ A = \left[ \begin{array}{ccc}
+    c_0 & c_1 & 0 \\
+    c_2 & c_3 & 0 \\
+    c_4 & c_5 & 1 \end{array} \right]\f]
+ *
+ * We then interpret points as row vectors of the form \f$[p_X, p_Y, 1]\f$. Applying a
+ * transformation to a point can be written as right-multiplication by a 3x3 matrix
+ * (\f$p' = pA\f$). This subset of matrices is closed under multiplication - combination
+ * of any two transforms can be expressed as the multiplication of their matrices.
+ * In this representation, the \f$c_4\f$ and \f$c_5\f$ coefficients represent
+ * the translation component of the transformation.
+ *
+ * Matrices can be multiplied by other more specific transformations. When multiplying,
+ * the transformations are applied from left to right, so for example <code>m = a * b</code>
+ * means: @a m first transforms by a, then by b.
+ *
+ * @ingroup Transforms
+ */
+class Affine
+    : boost::equality_comparable< Affine // generates operator!= from operator==
+    , boost::multipliable< Affine
+    , boost::multipliable< Affine, Translate
+    , boost::multipliable< Affine, Scale
+    , boost::multipliable< Affine, Rotate
+    , boost::multipliable< Affine, HShear
+    , boost::multipliable< Affine, VShear
+      > > > > > > >
+    // boost::multipliable< A, B > generates operator*(A const &, B const &)
+    // and operator*(B const &, A const &) from A::operator*=(B const &)
+{
+    Coord _c[6];
+public:
+    Affine() {
+        _c[0] = _c[3] = 1;
+        _c[1] = _c[2] = _c[4] = _c[5] = 0;
+    }
+
+    Affine(Affine const &m) {
+        for(int i = 0; i < 6; i++) {
+            _c[i] = m[i];
+        }
+    }
+
+    /** @brief Create a matrix from its coefficient values.
+     * It's rather inconvenient to directly create matrices in this way. Use transform classes
+     * if your transformation has a geometric interpretation.
+     * @see Translate
+     * @see Scale
+     * @see Rotate
+     * @see HShear
+     * @see VShear */
+    Affine(Coord c0, Coord c1, Coord c2, Coord c3, Coord c4, Coord c5) {
+        _c[0] = c0; _c[1] = c1;
+        _c[2] = c2; _c[3] = c3;
+        _c[4] = c4; _c[5] = c5;
+    }
+
+    Affine &operator=(Affine const &m) {
+        for(int i = 0; i < 6; i++)
+            _c[i] = m._c[i];
+        return *this;
+    }
+
+    /** @brief Access a coefficient by its index. */
+    inline Coord operator[](unsigned i) const { return _c[i]; }
+    inline Coord &operator[](unsigned i) { return _c[i]; }
+
+    /// @name Combine with other transformations
+    /// @{
+    Affine &operator*=(Affine const &m);
+    // implemented in transforms.cpp
+    Affine &operator*=(Translate const &t);
+    Affine &operator*=(Scale const &s);
+    Affine &operator*=(Rotate const &r);
+    Affine &operator*=(HShear const &h);
+    Affine &operator*=(VShear const &v);
+    /// @}
+
+    bool operator==(Affine const &o) const {
+        for(unsigned i = 0; i < 6; ++i) {
+            if ( _c[i] != o._c[i] ) return false;
+        }
+        return true;
+    }
+
+    /// @name Get the parameters of the matrix's transform
+    /// @{
+    Point xAxis() const;
+    Point yAxis() const;
+    Point translation() const;
+    Coord expansionX() const;
+    Coord expansionY() const;
+    Point expansion() const { return Point(expansionX(), expansionY()); }
+    /// @}
+    
+    /// @name Modify the matrix
+    /// @{
+    void setXAxis(Point const &vec);
+    void setYAxis(Point const &vec);
+
+    void setTranslation(Point const &loc);
+
+    void setExpansionX(Coord val);
+    void setExpansionY(Coord val);
+    void setIdentity();
+    /// @}
+
+    /// @name Inspect the matrix's transform
+    /// @{
+    bool isIdentity(Coord eps = EPSILON) const;
+
+    bool isTranslation(Coord eps = EPSILON) const;
+    bool isScale(Coord eps = EPSILON) const;
+    bool isUniformScale(Coord eps = EPSILON) const;
+    bool isRotation(Coord eps = EPSILON) const;
+    bool isHShear(Coord eps = EPSILON) const;
+    bool isVShear(Coord eps = EPSILON) const;
+
+    bool isNonzeroTranslation(Coord eps = EPSILON) const;
+    bool isNonzeroScale(Coord eps = EPSILON) const;
+    bool isNonzeroUniformScale(Coord eps = EPSILON) const;
+    bool isNonzeroRotation(Coord eps = EPSILON) const;
+    bool isNonzeroHShear(Coord eps = EPSILON) const;
+    bool isNonzeroVShear(Coord eps = EPSILON) const;
+
+    bool isZoom(Coord eps = EPSILON) const;
+    bool preservesArea(Coord eps = EPSILON) const;
+    bool preservesAngles(Coord eps = EPSILON) const;
+    bool preservesDistances(Coord eps = EPSILON) const;
+    bool flips() const;
+
+    bool isSingular(Coord eps = EPSILON) const;
+    /// @}
+
+    /// @name Compute other matrices
+    /// @{
+    Affine withoutTranslation() const {
+        Affine ret(*this);
+        ret.setTranslation(Point(0,0));
+        return ret;
+    }
+    Affine inverse() const;
+    /// @}
+
+    /// @name Compute scalar values
+    /// @{
+    Coord det() const;
+    Coord descrim2() const;
+    Coord descrim() const;
+    /// @}
+    inline static Affine identity();
+};
+
+/** @brief Print out the Affine (for debugging).
+ * @relates Affine */
+inline std::ostream &operator<< (std::ostream &out_file, const Geom::Affine &m) {
+    out_file << "A: " << m[0] << "  C: " << m[2] << "  E: " << m[4] << "\n";
+    out_file << "B: " << m[1] << "  D: " << m[3] << "  F: " << m[5] << "\n";
+    return out_file;
+}
+
+/** Given a matrix m such that unit_circle = m*x, this returns the
+ * quadratic form x*A*x = 1.
+ * @relates Affine */
+Affine elliptic_quadratic_form(Affine const &m);
+
+/** Given a matrix (ignoring the translation) this returns the eigen
+ * values and vectors. */
+class Eigen{
+public:
+    Point vectors[2];
+    double values[2];
+    Eigen(Affine const &m);
+    Eigen(double M[2][2]);
+};
+
+// Affine factories
+Affine from_basis(const Point x_basis, const Point y_basis, const Point offset=Point(0,0));
+
+/** @brief Create an identity matrix.
+ * This is a convenience function identical to Affine::identity(). */
+inline Affine identity() {
+    Affine ret(Affine::identity());
+    return ret; // allow NRVO
+}
+
+/** @brief Create an identity matrix.
+ * @return The matrix
+ *         \f$\left[\begin{array}{ccc}
+           1 & 0 & 0 \\
+           0 & 1 & 0 \\
+           0 & 0 & 1 \end{array}\right]\f$.
+ * @relates Affine */
+inline Affine Affine::identity() {
+    Affine ret(1.0, 0.0,
+               0.0, 1.0,
+               0.0, 0.0);
+    return ret; // allow NRVO
+}
+
+} /* namespace Geom */
+
+#endif /* !__Geom_MATRIX_H__ */
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/math-utils.h b/src/2geom/math-utils.h
new file mode 100644
index 000000000..2c348f54b
--- /dev/null
+++ b/src/2geom/math-utils.h
@@ -0,0 +1,108 @@
+#ifndef LIB2GEOM_MATH_UTILS_HEADER
+#define LIB2GEOM_MATH_UTILS_HEADER
+
+/**
+ * \file
+ * \brief Low level math functions and compatibility wrappers
+ *//*
+ * Authors:
+ *   Johan Engelen <goejendaagh@zonnet.nl>
+ *   Michael G. Sloan <mgsloan@gmail.com>
+ *   Krzysztof Kosiński <tweenk.pl@gmail.com>
+ * Copyright 2006-2009 Authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ *
+ */
+
+#include "config.h"
+#include <math.h> // sincos is usually only available in math.h
+#include <cmath>
+#include <utility> // for std::pair
+
+namespace Geom {
+
+/** @brief Sign function - indicates the sign of a numeric type.
+ * Mathsy people will know this is basically the derivative of abs, except for the fact
+ * that it is defined on 0.
+ * @return -1 when x is negative, 1 when positive, and 0 if equal to 0. */
+template <class T> inline int sgn(const T& x) {
+    return (x < 0 ? -1 : (x > 0 ? 1 : 0) );
+}
+
+template <class T> inline T sqr(const T& x) {return x * x;}
+template <class T> inline T cube(const T& x) {return x * x * x;}
+
+/** Between function - returns true if a number x is within a range: (min < x) && (max > x).
+ * The values delimiting the range and the number must have the same type.
+ */
+template <class T> inline const T& between (const T& min, const T& max, const T& x)
+    { return (min < x) && (max > x); }
+
+/** @brief Returns @a x rounded to the nearest multiple of \f$10^{p}\f$.
+
+    Implemented in terms of round, i.e. we make no guarantees as to what happens if x is
+    half way between two rounded numbers.
+    
+    Note: places is the number of decimal places without using scientific (e) notation, not the
+    number of significant figures.  This function may not be suitable for values of x whose
+    magnitude is so far from 1 that one would want to use scientific (e) notation.
+
+    places may be negative: e.g. places = -2 means rounding to a multiple of .01
+**/
+inline double decimal_round(double x, int p) {
+    //TODO: possibly implement with modulus instead?
+    double const multiplier = ::pow(10.0, p);
+    return ::round( x * multiplier ) / multiplier;
+}
+
+/** @brief Simultaneously compute a sine and a cosine of the same angle.
+ * This function can be up to 2 times faster than separate computation, depending
+ * on the platform. It uses the standard library function sincos() if available.
+ * @param angle Angle
+ * @param sin_ Variable that will store the sine
+ * @param cos_ Variable that will store the cosine */
+inline void sincos(double angle, double &sin_, double &cos_) {
+#ifdef HAVE_SINCOS
+    ::sincos(angle, &sin_, &cos_);
+#else
+    sin_ = ::sin(angle);
+    cos_ = ::cos(angle);
+#endif
+}
+
+}
+
+#endif
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :