1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
|
#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 "utils.h"
#include "matrix.h"
#include "point.h"
namespace Geom {
/** Creates a Matrix 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 Matrix.
*/
//NOTE: Inkscape's version is broken, so when including this version, you'll have to search for code with this func
//TODO: move to Matrix::from_basis
Matrix from_basis(Point const x_basis, Point const y_basis, Point const offset) {
return Matrix(x_basis[X], x_basis[Y],
y_basis[X], y_basis[Y],
offset [X], offset [Y]);
}
Point Matrix::xAxis() const {
return Point(_c[0], _c[1]);
}
Point Matrix::yAxis() const {
return Point(_c[2], _c[3]);
}
/** Gets the translation imparted by the Matrix.
*/
Point Matrix::translation() const {
return Point(_c[4], _c[5]);
}
void Matrix::setXAxis(Point const &vec) {
for(int i = 0; i < 2; i++)
_c[i] = vec[i];
}
void Matrix::setYAxis(Point const &vec) {
for(int i = 0; i < 2; i++)
_c[i + 2] = vec[i];
}
/** Sets the translation imparted by the Matrix.
*/
void Matrix::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 Matrix. 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 Matrix::expansionX() const {
return sqrt(_c[0] * _c[0] + _c[1] * _c[1]);
}
/** Calculates the amount of y-scaling imparted by the Matrix. 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 Matrix::expansionY() const {
return sqrt(_c[2] * _c[2] + _c[3] * _c[3]);
}
void Matrix::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 Matrix::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 Matrix. */
void Matrix::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;
}
//TODO: use eps
bool Matrix::isIdentity(Coord const eps) const {
return are_near(_c[0], 1.0) && are_near(_c[1], 0.0) &&
are_near(_c[2], 0.0) && are_near(_c[3], 1.0) &&
are_near(_c[4], 0.0) && are_near(_c[5], 0.0);
}
/** Answers the question "Does this matrix perform a translation, and \em{only} a translation?"
\param eps an epsilon value defaulting to EPSILON
\return A bool representing yes/no.
*/
bool Matrix::isTranslation(Coord const eps) const {
return are_near(_c[0], 1.0) && are_near(_c[1], 0.0) &&
are_near(_c[2], 0.0) && are_near(_c[3], 1.0) &&
(!are_near(_c[4], 0.0) || !are_near(_c[5], 0.0));
}
/** Answers the question "Does this matrix perform a scale, and \em{only} a Scale?"
\param eps an epsilon value defaulting to EPSILON
\return A bool representing yes/no.
*/
bool Matrix::isScale(Coord const eps) const {
return !are_near(_c[0], 1.0) || !are_near(_c[3], 1.0) && //NOTE: these are the diags, and the next line opposite diags
are_near(_c[1], 0.0) && are_near(_c[2], 0.0) &&
are_near(_c[4], 0.0) && are_near(_c[5], 0.0);
}
/** Answers the question "Does this matrix perform a uniform scale, and \em{only} a uniform scale?"
\param eps an epsilon value defaulting to EPSILON
\return A bool representing yes/no.
*/
bool Matrix::isUniformScale(Coord const eps) const {
return !are_near(_c[0], 1.0) && are_near(_c[0], _c[3]) &&
are_near(_c[1], 0.0) && are_near(_c[2], 0.0) &&
are_near(_c[4], 0.0) && are_near(_c[5], 0.0);
}
/** Answers the question "Does this matrix perform a rotation, and \em{only} a rotation?"
\param eps an epsilon value defaulting to EPSILON
\return A bool representing yes/no.
*/
bool Matrix::isRotation(Coord const eps) const {
return !are_near(_c[0], _c[3]) && are_near(_c[1], -_c[2]) &&
are_near(_c[4], 0.0) && are_near(_c[5], 0.0) &&
are_near(_c[0]*_c[0] + _c[1]*_c[1], 1.0);
}
bool Matrix::onlyScaleAndTranslation(Coord const eps) const {
return are_near(_c[0], _c[3]) && are_near(_c[1], 0) && are_near(_c[2], 0);
}
bool Matrix::flips() const {
return cross(xAxis(), yAxis()) > 0;
}
/** Returns the Scale/Rotate/skew part of the matrix without the translation part. */
Matrix Matrix::without_translation() const {
return Matrix(_c[0], _c[1], _c[2], _c[3], 0, 0);
}
/** Attempts to calculate the inverse of a matrix.
* This is a Matrix such that m * m.inverse() is very near (hopefully < epsilon difference) the identity Matrix.
* \textbf{The Identity Matrix is returned if the matrix has no inverse.}
\return The inverse of the Matrix if defined, otherwise the Identity Matrix.
*/
Matrix Matrix::inverse() const {
Matrix d;
Geom::Coord const determ = det();
if (!are_near(determ, 0.0)) {
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();
}
return d;
}
/** Calculates the determinant of a Matrix. */
Geom::Coord Matrix::det() const {
return _c[0] * _c[3] - _c[1] * _c[2];
}
/** Calculates the scalar of the descriminant of the Matrix.
* This is simply the absolute value of the determinant.
*/
Geom::Coord Matrix::descrim2() const {
return fabs(det());
}
/** Calculates the descriminant of the Matrix. */
Geom::Coord Matrix::descrim() const {
return sqrt(descrim2());
}
Matrix operator*(Matrix const &m0, Matrix const &m1) {
Matrix ret;
for(int a = 0; a < 5; a += 2) {
for(int b = 0; b < 2; b++) {
ret[a + b] = m0[a] * m1[b] + m0[a + 1] * m1[b + 2];
}
}
ret[4] += m1[4];
ret[5] += m1[5];
return ret;
}
//TODO: What's this!?!
Matrix elliptic_quadratic_form(Matrix const &m) {
double const od = m[0] * m[1] + m[2] * m[3];
return Matrix(m[0]*m[0] + m[1]*m[1], od,
od, m[2]*m[2] + m[3]*m[3],
0, 0);
}
Eigen::Eigen(Matrix 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])));
}
}
} //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:encoding=utf-8:textwidth=99 :
|
