affine.cpp

Go to the documentation of this file.

/*
 * Authors:
 *   Lauris Kaplinski <lauris@kaplinski.com>
 *   Michael G. Sloan <mgsloan@gmail.com>
 *
 * This code is in public domain
 */
 
#include <2geom/affine.h>
#include <2geom/point.h>
#include <2geom/polynomial.h>
#include <2geom/utils.h>
 
namespace Geom {
 
//NOTE: Inkscape's version is broken, so when including this version, you'll have to search for code with this func
Affine from_basis(Point const &x_basis, Point const &y_basis, Point const &offset) {
    return Affine(x_basis[X], x_basis[Y],
                  y_basis[X], y_basis[Y],
                  offset [X], offset [Y]);
}
 
Point Affine::xAxis() const {
    return Point(_c[0], _c[1]);
}
 
Point Affine::yAxis() const {
    return Point(_c[2], _c[3]);
}
 
Point Affine::translation() const {
    return Point(_c[4], _c[5]);
}
 
void Affine::setXAxis(Point const &vec) {
    for(int i = 0; i < 2; i++)
        _c[i] = vec[i];
}
 
void Affine::setYAxis(Point const &vec) {
    for(int i = 0; i < 2; i++)
        _c[i + 2] = vec[i];
}
 
void Affine::setTranslation(Point const &loc) {
    for(int i = 0; i < 2; i++)
        _c[i + 4] = loc[i];
}
 
double Affine::expansionX() const {
    return sqrt(_c[0] * _c[0] + _c[1] * _c[1]);
}
 
double Affine::expansionY() const {
    return sqrt(_c[2] * _c[2] + _c[3] * _c[3]);
}
 
void Affine::setExpansionX(double val) {
    double exp_x = expansionX();
    if (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 Affine::setExpansionY(double val) {
    double exp_y = expansionY();
    if (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;
        }
    }
}
 
void Affine::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;
}
 
bool Affine::isIdentity(Coord 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);
}
 
bool Affine::isTranslation(Coord 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);
}
bool Affine::isNonzeroTranslation(Coord 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));
}
 
bool Affine::isScale(Coord eps) const {
    if (isSingular(eps)) return false;
    return 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);
}
 
bool Affine::isNonzeroScale(Coord eps) const {
    if (isSingular(eps)) return false;
    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);
}
 
bool Affine::isUniformScale(Coord eps) const {
    if (isSingular(eps)) return false;
    return are_near(fabs(_c[0]), fabs(_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);
}
 
bool Affine::isNonzeroUniformScale(Coord eps) const {
    if (isSingular(eps)) return false;
    // we need to test both c0 and c3 to handle the case of flips,
    // which should be treated as nonzero uniform scales
    return !(are_near(_c[0], 1.0, eps) && are_near(_c[3], 1.0, eps)) &&
           are_near(fabs(_c[0]), fabs(_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);
}
 
bool Affine::isRotation(Coord 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 Affine::isNonzeroRotation(Coord eps) const {
    return !are_near(_c[0], 1.0, eps) &&
           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 Affine::isNonzeroNonpureRotation(Coord eps) const {
    return !are_near(_c[0], 1.0, eps) &&
           are_near(_c[0], _c[3], eps) && are_near(_c[1], -_c[2], eps) &&
           are_near(_c[0]*_c[0] + _c[1]*_c[1], 1.0, eps);
}
 
Point Affine::rotationCenter() const {
    Coord x = (_c[2]*_c[5]+_c[4]-_c[4]*_c[3]) / (1-_c[3]-_c[0]+_c[0]*_c[3]-_c[2]*_c[1]);
    Coord y = (_c[1]*x + _c[5]) / (1 - _c[3]);
    return Point(x,y);
};
 
bool Affine::isHShear(Coord eps) const {
    return are_near(_c[0], 1.0, eps) && are_near(_c[1], 0.0, eps) &&
           are_near(_c[3], 1.0, eps) && are_near(_c[4], 0.0, eps) &&
           are_near(_c[5], 0.0, eps);
}
bool Affine::isNonzeroHShear(Coord 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);
}
 
bool Affine::isVShear(Coord eps) const {
    return are_near(_c[0], 1.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);
}
 
bool Affine::isNonzeroVShear(Coord 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);
}
 
bool Affine::isZoom(Coord eps) const {
    if (isSingular(eps)) return false;
    return are_near(_c[0], _c[3], eps) && are_near(_c[1], 0, eps) && are_near(_c[2], 0, eps);
}
 
bool Affine::preservesArea(Coord eps) const
{
    return are_near(descrim2(), 1.0, eps);
}
 
bool Affine::preservesAngles(Coord eps) const
{
    if (isSingular(eps)) return false;
    return (are_near(_c[0], _c[3], eps) && are_near(_c[1], -_c[2], eps)) ||
           (are_near(_c[0], -_c[3], eps) && are_near(_c[1], _c[2], eps));
}
 
bool Affine::preservesDistances(Coord eps) const
{
    return ((are_near(_c[0], _c[3], eps) && are_near(_c[1], -_c[2], eps)) ||
            (are_near(_c[0], -_c[3], eps) && are_near(_c[1], _c[2], eps))) &&
           are_near(_c[0] * _c[0] + _c[1] * _c[1], 1.0, eps);
}
 
bool Affine::flips() const {
    return det() < 0;
}
 
bool Affine::isSingular(Coord eps) const {
    return are_near(det(), 0.0, eps);
}
 
Affine Affine::inverse() const {
    Affine d;
    
    double mx = std::max(fabs(_c[0]) + fabs(_c[1]), 
                         fabs(_c[2]) + fabs(_c[3])); // a random matrix norm (either l1 or linfty
    if(mx > 0) {
        Geom::Coord const determ = det();
        if (!rel_error_bound(std::sqrt(fabs(determ)), mx)) {
            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();
        }
    } else {
        d.setIdentity();
    }
 
    return d;
}
 
Coord Affine::det() const {
    // TODO this can overflow
    return _c[0] * _c[3] - _c[1] * _c[2];
}
 
Coord Affine::descrim2() const {
    return fabs(det());
}
 
Coord Affine::descrim() const {
    return sqrt(descrim2());
}
 
Affine &Affine::operator*=(Affine const &o) {
    Coord nc[6];
    for(int a = 0; a < 5; a += 2) {
        for(int b = 0; b < 2; b++) {
            nc[a + b] = _c[a] * o._c[b] + _c[a + 1] * o._c[b + 2];
        }
    }
    for(int a = 0; a < 6; ++a) {
        _c[a] = nc[a];
    }
    _c[4] += o._c[4];
    _c[5] += o._c[5];
    return *this;
}
 
//TODO: What's this!?!
Affine elliptic_quadratic_form(Affine const &m) {
    double od = m[0] * m[1]  +  m[2] * m[3];
    Affine ret (m[0]*m[0] + m[1]*m[1], od,
                od, m[2]*m[2] + m[3]*m[3],
                0, 0);
    return ret; // allow NRVO
}
 
Eigen::Eigen(Affine const &m) {
    double const B = -m[0] - m[3];
    double const C = m[0]*m[3] - m[1]*m[2];
 
    std::vector<double> v = solve_quadratic(1, B, C);
 
    for (unsigned i = 0; i < v.size(); ++i) {
        values[i] = v[i];
        vectors[i] = unit_vector(rot90(Point(m[0] - values[i], m[1])));
    }
    for (unsigned i = v.size(); i < 2; ++i) {
        values[i] = 0;
        vectors[i] = Point(0,0);
    }
}
 
Eigen::Eigen(double m[2][2]) {
    double const B = -m[0][0] - m[1][1];
    double const C = m[0][0]*m[1][1] - m[1][0]*m[0][1];
 
    std::vector<double> v = solve_quadratic(1, B, C);
 
    for (unsigned i = 0; i < v.size(); ++i) {
        values[i] = v[i];
        vectors[i] = unit_vector(rot90(Point(m[0][0] - values[i], m[0][1])));
    }
    for (unsigned i = v.size(); i < 2; ++i) {
        values[i] = 0;
        vectors[i] = Point(0,0);
    }
}
 
bool are_near(Affine const &a, Affine const &b, Coord eps)
{
    return are_near(a[0], b[0], eps) && are_near(a[1], b[1], eps) &&
           are_near(a[2], b[2], eps) && are_near(a[3], b[3], eps) &&
           are_near(a[4], b[4], eps) && are_near(a[5], b[5], eps);
}
 
}  //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:fileencoding=utf-8:textwidth=99 :