matrix.cpp revision a7db6de6cc11de8b084ef3b9061a4b40a9c08a2b
#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/matrix.h>
#include <2geom/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, 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);
}
/** 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, 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));
}
/** 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, 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);
}
/** 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, 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);
}
/** 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], 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);
}
bool Matrix::onlyScaleAndTranslation(Coord const eps) const {
return are_near(_c[0], _c[3], eps) && are_near(_c[1], 0, eps) && are_near(_c[2], 0, eps);
}
bool Matrix::isSingular(Coord const eps) const {
return are_near(det(), 0.0, eps);
}
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 :