1 #define __Geom_MATRIX_C__
3 /** \file
4 * Various matrix routines. Currently includes some Geom::Rotate etc. routines too.
5 */
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 */
15 #include "utils.h"
16 #include "matrix.h"
17 #include "point.h"
19 namespace Geom {
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.
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 }
38 Point Matrix::xAxis() const {
39 return Point(_c[0], _c[1]);
40 }
42 Point Matrix::yAxis() const {
43 return Point(_c[2], _c[3]);
44 }
46 /** Gets the translation imparted by the Matrix.
47 */
48 Point Matrix::translation() const {
49 return Point(_c[4], _c[5]);
50 }
52 void Matrix::setXAxis(Point const &vec) {
53 for(int i = 0; i < 2; i++)
54 _c[i] = vec[i];
55 }
57 void Matrix::setYAxis(Point const &vec) {
58 for(int i = 0; i < 2; i++)
59 _c[i + 2] = vec[i];
60 }
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 }
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 }
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 }
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 }
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 }
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 }
108 //TODO: use eps
110 bool Matrix::isIdentity(Coord const eps) const {
111 return are_near(_c[0], 1.0) && are_near(_c[1], 0.0) &&
112 are_near(_c[2], 0.0) && are_near(_c[3], 1.0) &&
113 are_near(_c[4], 0.0) && are_near(_c[5], 0.0);
114 }
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) && are_near(_c[1], 0.0) &&
122 are_near(_c[2], 0.0) && are_near(_c[3], 1.0) &&
123 (!are_near(_c[4], 0.0) || !are_near(_c[5], 0.0));
124 }
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) || !are_near(_c[3], 1.0) && //NOTE: these are the diags, and the next line opposite diags
132 are_near(_c[1], 0.0) && are_near(_c[2], 0.0) &&
133 are_near(_c[4], 0.0) && are_near(_c[5], 0.0);
134 }
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) && are_near(_c[0], _c[3]) &&
142 are_near(_c[1], 0.0) && are_near(_c[2], 0.0) &&
143 are_near(_c[4], 0.0) && are_near(_c[5], 0.0);
144 }
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]) && are_near(_c[1], -_c[2]) &&
152 are_near(_c[4], 0.0) && are_near(_c[5], 0.0) &&
153 are_near(_c[0]*_c[0] + _c[1]*_c[1], 1.0);
154 }
156 bool Matrix::onlyScaleAndTranslation(Coord const eps) const {
157 return are_near(_c[0], _c[3]) && are_near(_c[1], 0) && are_near(_c[2], 0);
158 }
160 bool Matrix::flips() const {
161 return cross(xAxis(), yAxis()) > 0;
162 }
164 /** Returns the Scale/Rotate/skew part of the matrix without the translation part. */
165 Matrix Matrix::without_translation() const {
166 return Matrix(_c[0], _c[1], _c[2], _c[3], 0, 0);
167 }
169 /** Attempts to calculate the inverse of a matrix.
170 * This is a Matrix such that m * m.inverse() is very near (hopefully < epsilon difference) the identity Matrix.
171 * \textbf{The Identity Matrix is returned if the matrix has no inverse.}
172 \return The inverse of the Matrix if defined, otherwise the Identity Matrix.
173 */
174 Matrix Matrix::inverse() const {
175 Matrix d;
177 Geom::Coord const determ = det();
178 if (!are_near(determ, 0.0)) {
179 Geom::Coord const ideterm = 1.0 / determ;
181 d._c[0] = _c[3] * ideterm;
182 d._c[1] = -_c[1] * ideterm;
183 d._c[2] = -_c[2] * ideterm;
184 d._c[3] = _c[0] * ideterm;
185 d._c[4] = -_c[4] * d._c[0] - _c[5] * d._c[2];
186 d._c[5] = -_c[4] * d._c[1] - _c[5] * d._c[3];
187 } else {
188 d.setIdentity();
189 }
191 return d;
192 }
194 /** Calculates the determinant of a Matrix. */
195 Geom::Coord Matrix::det() const {
196 return _c[0] * _c[3] - _c[1] * _c[2];
197 }
199 /** Calculates the scalar of the descriminant of the Matrix.
200 * This is simply the absolute value of the determinant.
201 */
202 Geom::Coord Matrix::descrim2() const {
203 return fabs(det());
204 }
206 /** Calculates the descriminant of the Matrix. */
207 Geom::Coord Matrix::descrim() const {
208 return sqrt(descrim2());
209 }
211 Matrix operator*(Matrix const &m0, Matrix const &m1) {
212 Matrix ret;
213 for(int a = 0; a < 5; a += 2) {
214 for(int b = 0; b < 2; b++) {
215 ret[a + b] = m0[a] * m1[b] + m0[a + 1] * m1[b + 2];
216 }
217 }
218 ret[4] += m1[4];
219 ret[5] += m1[5];
220 return ret;
221 }
223 //TODO: What's this!?!
224 Matrix elliptic_quadratic_form(Matrix const &m) {
225 double const od = m[0] * m[1] + m[2] * m[3];
226 return Matrix(m[0]*m[0] + m[1]*m[1], od,
227 od, m[2]*m[2] + m[3]*m[3],
228 0, 0);
229 }
231 Eigen::Eigen(Matrix const &m) {
232 double const B = -m[0] - m[3];
233 double const C = m[0]*m[3] - m[1]*m[2];
234 double const center = -B/2.0;
235 double const delta = sqrt(B*B-4*C)/2.0;
236 values[0] = center + delta; values[1] = center - delta;
237 for (int i = 0; i < 2; i++) {
238 vectors[i] = unit_vector(rot90(Point(m[0]-values[i], m[1])));
239 }
240 }
242 } //namespace Geom
244 /*
245 Local Variables:
246 mode:c++
247 c-file-style:"stroustrup"
248 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
249 indent-tabs-mode:nil
250 fill-column:99
251 End:
252 */
253 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :