src/libnr/nr-matrix.cpp

   1 #define __NR_MATRIX_C__
   2
   3 /** \file
   4  * Various matrix routines.  Currently includes some NR::rotate etc. routines too.
   5  */
   6
   7 /*
   8  * Authors:
   9  *   Lauris Kaplinski <lauris@kaplinski.com>
  10  *
  11  * This code is in public domain
  12  */
  13
  14 #include <cstdlib>
  15 #include "nr-matrix.h"
  16 #include "nr-matrix-fns.h"
  17
  18
  19
  20 /**
  21  *  Implement NR functions and methods
  22  */
  23 namespace NR {
  24
  25
  26
  27
  28
  29 /**
  30  *  Multiply two matrices together
  31  */
  32 Matrix operator*(Matrix const &m0, Matrix const &m1)
  33 {
  34     NR::Coord const d0 = m0[0] * m1[0]  +  m0[1] * m1[2];
  35     NR::Coord const d1 = m0[0] * m1[1]  +  m0[1] * m1[3];
  36     NR::Coord const d2 = m0[2] * m1[0]  +  m0[3] * m1[2];
  37     NR::Coord const d3 = m0[2] * m1[1]  +  m0[3] * m1[3];
  38     NR::Coord const d4 = m0[4] * m1[0]  +  m0[5] * m1[2]  +  m1[4];
  39     NR::Coord const d5 = m0[4] * m1[1]  +  m0[5] * m1[3]  +  m1[5];
  40
  41     Matrix ret( d0, d1, d2, d3, d4, d5 );
  42
  43     return ret;
  44 }
  45
  46
  47
  48
  49
  50 /**
  51  *  Return the inverse of this matrix.  If an inverse is not defined,
  52  *  then return the identity matrix.
  53  */
  54 Matrix Matrix::inverse() const
  55 {
  56     Matrix d(0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
  57
  58     NR::Coord const det = _c[0] * _c[3] - _c[1] * _c[2];
  59     if (!NR_DF_TEST_CLOSE(det, 0.0, NR_EPSILON)) {
  60
  61         NR::Coord const idet = 1.0 / det;
  62         NR::Coord *dest = d._c;
  63
  64         /*0*/ *dest++ =  _c[3] * idet;
  65         /*1*/ *dest++ = -_c[1] * idet;
  66         /*2*/ *dest++ = -_c[2] * idet;
  67         /*3*/ *dest++ =  _c[0] * idet;
  68         /*4*/ *dest++ = -_c[4] * d._c[0] - _c[5] * d._c[2];
  69         /*5*/ *dest   = -_c[4] * d._c[1] - _c[5] * d._c[3];
  70
  71     } else {
  72         d.set_identity();
  73     }
  74
  75     return d;
  76 }
  77
  78
  79
  80
  81
  82 /**
  83  *  Set this matrix to Identity
  84  */
  85 void Matrix::set_identity()
  86 {
  87     NR::Coord *dest = _c;
  88
  89     *dest++ = 1.0; //0
  90     *dest++ = 0.0; //1
  91     *dest++ = 0.0; //2
  92     *dest++ = 1.0; //3
  93     // translation
  94     *dest++ = 0.0; //4
  95     *dest   = 0.0; //5
  96 }
  97
  98
  99
 100
 101
 102 /**
 103  *  return an Identity matrix
 104  */
 105 Matrix identity()
 106 {
 107     Matrix ret(1.0, 0.0,
 108                0.0, 1.0,
 109                0.0, 0.0);
 110     return ret;
 111 }
 112
 113
 114
 115
 116
 117 /**
 118  *
 119  */
 120 Matrix from_basis(Point const x_basis, Point const y_basis, Point const offset)
 121 {
 122     Matrix const ret(x_basis[X], y_basis[X],
 123                      x_basis[Y], y_basis[Y],
 124                      offset[X], offset[Y]);
 125     return ret;
 126 }
 127
 128
 129
 130
 131 /**
 132  * Returns a rotation matrix corresponding by the specified angle (in radians) about the origin.
 133  *
 134  * \see NR::rotate_degrees
 135  *
 136  * Angle direction in Inkscape code: If you use the traditional mathematics convention that y
 137  * increases upwards, then positive angles are anticlockwise as per the mathematics convention.  If
 138  * you take the common non-mathematical convention that y increases downwards, then positive angles
 139  * are clockwise, as is common outside of mathematics.
 140  */
 141 rotate::rotate(NR::Coord const theta) :
 142     vec(cos(theta),
 143         sin(theta))
 144 {
 145 }
 146
 147
 148
 149
 150
 151 /**
 152  *  Return the determinant of the Matrix
 153  */
 154 NR::Coord Matrix::det() const
 155 {
 156     return _c[0] * _c[3] - _c[1] * _c[2];
 157 }
 158
 159
 160
 161
 162
 163 /**
 164  * Return the scalar of the descriminant of the Matrix
 165  */
 166 NR::Coord Matrix::descrim2() const
 167 {
 168     return fabs(det());
 169 }
 170
 171
 172
 173
 174
 175 /**
 176  *  Return the descriminant of the Matrix
 177  */
 178 NR::Coord Matrix::descrim() const
 179 {
 180     return sqrt(descrim2());
 181 }
 182
 183
 184
 185
 186
 187 /**
 188  *
 189  */
 190 bool Matrix::is_translation(Coord const eps) const {
 191     return ( fabs(_c[0] - 1.0) < eps &&
 192              fabs(_c[3] - 1.0) < eps &&
 193              fabs(_c[1])       < eps &&
 194              fabs(_c[2])       < eps   );
 195 }
 196
 197
 198 /**
 199  *
 200  */
 201 bool Matrix::is_scale(Coord const eps) const {
 202     return ( (fabs(_c[0] - 1.0) > eps || fabs(_c[3] - 1.0) > eps) &&
 203              fabs(_c[1])       < eps &&
 204              fabs(_c[2])       < eps   );
 205 }
 206
 207
 208 /**
 209  *
 210  */
 211 bool Matrix::is_rotation(Coord const eps) const {
 212     return ( fabs(_c[1]) > eps &&
 213              fabs(_c[2]) > eps &&
 214              fabs(_c[1] + _c[2]) < 2 * eps);
 215 }
 216
 217
 218
 219
 220
 221 /**
 222  *  test whether the matrix is the identity matrix (true).  (2geom's Matrix::isIdentity() does the same)
 223  */
 224 bool Matrix::test_identity() const {
 225     return matrix_equalp(*this, NR_MATRIX_IDENTITY, NR_EPSILON);
 226 }
 227
 228
 229
 230
 231
 232 /**
 233  * calculates the descriminant of the matrix. (Geom::Coord Matrix::descrim() does the same)
 234  */
 235 double expansion(Matrix const &m) {
 236     return sqrt(fabs(m.det()));
 237 }
 238
 239
 240
 241
 242
 243 /**
 244  *
 245  */
 246 bool transform_equalp(Matrix const &m0, Matrix const &m1, NR::Coord const epsilon) {
 247     return
 248         NR_DF_TEST_CLOSE(m0[0], m1[0], epsilon) &&
 249         NR_DF_TEST_CLOSE(m0[1], m1[1], epsilon) &&
 250         NR_DF_TEST_CLOSE(m0[2], m1[2], epsilon) &&
 251         NR_DF_TEST_CLOSE(m0[3], m1[3], epsilon);
 252 }
 253
 254
 255
 256
 257
 258 /**
 259  *
 260  */
 261 bool translate_equalp(Matrix const &m0, Matrix const &m1, NR::Coord const epsilon) {
 262     return NR_DF_TEST_CLOSE(m0[4], m1[4], epsilon) && NR_DF_TEST_CLOSE(m0[5], m1[5], epsilon);
 263 }
 264
 265
 266
 267
 268
 269 /**
 270  *
 271  */
 272 bool matrix_equalp(Matrix const &m0, Matrix const &m1, NR::Coord const epsilon) {
 273     return transform_equalp(m0, m1, epsilon) && translate_equalp(m0, m1, epsilon);
 274 }
 275
 276
 277
 278 }  //namespace NR
 279
 280
 281
 282 /*
 283   Local Variables:
 284   mode:c++
 285   c-file-style:"stroustrup"
 286   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 287   indent-tabs-mode:nil
 288   fill-column:99
 289   End:
 290 */
 291 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :