src/2geom/matrix.cpp

   1 #define __Geom_MATRIX_C__
   2
   3 /** \file
   4  * Various matrix routines.  Currently includes some Geom::Rotate etc. routines too.
   5  */
   6
   7 /*
   8  * Authors:
   9  *   Lauris Kaplinski <lauris@kaplinski.com>
  10  *   Michael G. Sloan <mgsloan@gmail.com>
  11  *
  12  * This code is in public domain
  13  */
  14
  15 #include <2geom/utils.h>
  16 #include <2geom/matrix.h>
  17 #include <2geom/point.h>
  18
  19 namespace Geom {
  20
  21 /** Creates a Matrix given an axis and origin point.
  22  *  The axis is represented as two vectors, which represent skew, rotation, and scaling in two dimensions.
  23  *  from_basis(Point(1, 0), Point(0, 1), Point(0, 0)) would return the identity matrix.
  24
  25  \param x_basis the vector for the x-axis.
  26  \param y_basis the vector for the y-axis.
  27  \param offset the translation applied by the matrix.
  28  \return The new Matrix.
  29  */
  30 //NOTE: Inkscape's version is broken, so when including this version, you'll have to search for code with this func
  31 //TODO: move to Matrix::from_basis
  32 Matrix from_basis(Point const x_basis, Point const y_basis, Point const offset) {
  33     return Matrix(x_basis[X], x_basis[Y],
  34                   y_basis[X], y_basis[Y],
  35                   offset [X], offset [Y]);
  36 }
  37
  38 Point Matrix::xAxis() const {
  39     return Point(_c[0], _c[1]);
  40 }
  41
  42 Point Matrix::yAxis() const {
  43     return Point(_c[2], _c[3]);
  44 }
  45
  46 /** Gets the translation imparted by the Matrix.
  47  */
  48 Point Matrix::translation() const {
  49     return Point(_c[4], _c[5]);
  50 }
  51
  52 void Matrix::setXAxis(Point const &vec) {
  53     for(int i = 0; i < 2; i++)
  54         _c[i] = vec[i];
  55 }
  56
  57 void Matrix::setYAxis(Point const &vec) {
  58     for(int i = 0; i < 2; i++)
  59         _c[i + 2] = vec[i];
  60 }
  61
  62 /** Sets the translation imparted by the Matrix.
  63  */
  64 void Matrix::setTranslation(Point const &loc) {
  65     for(int i = 0; i < 2; i++)
  66         _c[i + 4] = loc[i];
  67 }
  68
  69 /** Calculates the amount of x-scaling imparted by the Matrix.  This is the scaling applied to
  70  *  the original x-axis region.  It is \emph{not} the overall x-scaling of the transformation.
  71  *  Equivalent to L2(m.xAxis())
  72  */
  73 double Matrix::expansionX() const {
  74     return sqrt(_c[0] * _c[0] + _c[1] * _c[1]);
  75 }
  76
  77 /** Calculates the amount of y-scaling imparted by the Matrix.  This is the scaling applied before
  78  *  the other transformations.  It is \emph{not} the overall y-scaling of the transformation.
  79  *  Equivalent to L2(m.yAxis())
  80  */
  81 double Matrix::expansionY() const {
  82     return sqrt(_c[2] * _c[2] + _c[3] * _c[3]);
  83 }
  84
  85 void Matrix::setExpansionX(double val) {
  86     double exp_x = expansionX();
  87     if(!are_near(exp_x, 0.0)) {  //TODO: best way to deal with it is to skip op?
  88         double coef = val / expansionX();
  89         for(unsigned i=0;i<2;i++) _c[i] *= coef;
  90     }
  91 }
  92
  93 void Matrix::setExpansionY(double val) {
  94     double exp_y = expansionY();
  95     if(!are_near(exp_y, 0.0)) {  //TODO: best way to deal with it is to skip op?
  96         double coef = val / expansionY();
  97         for(unsigned i=2; i<4; i++) _c[i] *= coef;
  98     }
  99 }
 100
 101 /** Sets this matrix to be the Identity Matrix. */
 102 void Matrix::setIdentity() {
 103     _c[0] = 1.0; _c[1] = 0.0;
 104     _c[2] = 0.0; _c[3] = 1.0;
 105     _c[4] = 0.0; _c[5] = 0.0;
 106 }
 107
 108 //TODO: use eps
 109
 110 bool Matrix::isIdentity(Coord const eps) const {
 111     return are_near(_c[0], 1.0, eps) && are_near(_c[1], 0.0, eps) &&
 112            are_near(_c[2], 0.0, eps) && are_near(_c[3], 1.0, eps) &&
 113            are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps);
 114 }
 115
 116 /** Answers the question "Does this matrix perform a translation, and \em{only} a translation?"
 117  \param eps an epsilon value defaulting to EPSILON
 118  \return A bool representing yes/no.
 119  */
 120 bool Matrix::isTranslation(Coord const eps) const {
 121     return are_near(_c[0], 1.0, eps) && are_near(_c[1], 0.0, eps) &&
 122            are_near(_c[2], 0.0, eps) && are_near(_c[3], 1.0, eps) &&
 123            (!are_near(_c[4], 0.0, eps) || !are_near(_c[5], 0.0, eps));
 124 }
 125
 126 /** Answers the question "Does this matrix perform a scale, and \em{only} a Scale?"
 127  \param eps an epsilon value defaulting to EPSILON
 128  \return A bool representing yes/no.
 129  */
 130 bool Matrix::isScale(Coord const eps) const {
 131     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
 132            are_near(_c[1], 0.0, eps) && are_near(_c[2], 0.0, eps) &&
 133            are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps);
 134 }
 135
 136 /** Answers the question "Does this matrix perform a uniform scale, and \em{only} a uniform scale?"
 137  \param eps an epsilon value defaulting to EPSILON
 138  \return A bool representing yes/no.
 139  */
 140 bool Matrix::isUniformScale(Coord const eps) const {
 141     return !are_near(_c[0], 1.0, eps) && are_near(_c[0], _c[3], eps) &&
 142            are_near(_c[1], 0.0, eps) && are_near(_c[2], 0.0, eps) &&
 143            are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps);
 144 }
 145
 146 /** Answers the question "Does this matrix perform a rotation, and \em{only} a rotation?"
 147  \param eps an epsilon value defaulting to EPSILON
 148  \return A bool representing yes/no.
 149  */
 150 bool Matrix::isRotation(Coord const eps) const {
 151     return are_near(_c[0], _c[3], eps) && are_near(_c[1], -_c[2], eps) &&
 152            are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps) &&
 153            are_near(_c[0]*_c[0] + _c[1]*_c[1], 1.0, eps);
 154 }
 155
 156 bool Matrix::onlyScaleAndTranslation(Coord const eps) const {
 157     return are_near(_c[0], _c[3], eps) && are_near(_c[1], 0, eps) && are_near(_c[2], 0, eps);
 158 }
 159
 160 bool Matrix::isSingular(Coord const eps) const {
 161     return are_near(det(), 0.0, eps);
 162 }
 163
 164 bool Matrix::flips() const {
 165     return cross(xAxis(), yAxis()) > 0;
 166 }
 167
 168 /** Returns the Scale/Rotate/skew part of the matrix without the translation part. */
 169 Matrix Matrix::without_translation() const {
 170     return Matrix(_c[0], _c[1], _c[2], _c[3], 0, 0);
 171 }
 172
 173 /** Attempts to calculate the inverse of a matrix.
 174  *  This is a Matrix such that m * m.inverse() is very near (hopefully < epsilon difference) the identity Matrix.
 175  *  \textbf{The Identity Matrix is returned if the matrix has no inverse.}
 176  \return The inverse of the Matrix if defined, otherwise the Identity Matrix.
 177  */
 178 Matrix Matrix::inverse() const {
 179     Matrix d;
 180
 181     Geom::Coord const determ = det();
 182     // the numerical precision of the determinant must be significant
 183     if (fabs(determ) > 1e-18) {
 184         Geom::Coord const ideterm = 1.0 / determ;
 185
 186         d._c[0] =  _c[3] * ideterm;
 187         d._c[1] = -_c[1] * ideterm;
 188         d._c[2] = -_c[2] * ideterm;
 189         d._c[3] =  _c[0] * ideterm;
 190         d._c[4] = -_c[4] * d._c[0] - _c[5] * d._c[2];
 191         d._c[5] = -_c[4] * d._c[1] - _c[5] * d._c[3];
 192     } else {
 193         d.setIdentity();
 194     }
 195
 196     return d;
 197 }
 198
 199 /** Calculates the determinant of a Matrix. */
 200 Geom::Coord Matrix::det() const {
 201     return _c[0] * _c[3] - _c[1] * _c[2];
 202 }
 203
 204 /** Calculates the scalar of the descriminant of the Matrix.
 205  *  This is simply the absolute value of the determinant.
 206  */
 207 Geom::Coord Matrix::descrim2() const {
 208     return fabs(det());
 209 }
 210
 211 /** Calculates the descriminant of the Matrix. */
 212 Geom::Coord Matrix::descrim() const {
 213     return sqrt(descrim2());
 214 }
 215
 216 Matrix operator*(Matrix const &m0, Matrix const &m1) {
 217     Matrix ret;
 218     for(int a = 0; a < 5; a += 2) {
 219         for(int b = 0; b < 2; b++) {
 220             ret[a + b] = m0[a] * m1[b] + m0[a + 1] * m1[b + 2];
 221         }
 222     }
 223     ret[4] += m1[4];
 224     ret[5] += m1[5];
 225     return ret;
 226 }
 227
 228 //TODO: What's this!?!
 229 Matrix elliptic_quadratic_form(Matrix const &m) {
 230     double const od = m[0] * m[1]  +  m[2] * m[3];
 231     return Matrix(m[0]*m[0] + m[1]*m[1], od,
 232                   od, m[2]*m[2] + m[3]*m[3],
 233                   0, 0);
 234 }
 235
 236 Eigen::Eigen(Matrix const &m) {
 237     double const B = -m[0] - m[3];
 238     double const C = m[0]*m[3] - m[1]*m[2];
 239     double const center = -B/2.0;
 240     double const delta = sqrt(B*B-4*C)/2.0;
 241     values[0] = center + delta; values[1] = center - delta;
 242     for (int i = 0; i < 2; i++) {
 243         vectors[i] = unit_vector(rot90(Point(m[0]-values[i], m[1])));
 244     }
 245 }
 246
 247 }  //namespace Geom
 248
 249 /*
 250   Local Variables:
 251   mode:c++
 252   c-file-style:"stroustrup"
 253   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 254   indent-tabs-mode:nil
 255   fill-column:99
 256   End:
 257 */
 258 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :