diff options
Diffstat (limited to 'kig/misc/kigtransform.cpp')
| -rw-r--r-- | kig/misc/kigtransform.cpp | 810 |
1 files changed, 810 insertions, 0 deletions
diff --git a/kig/misc/kigtransform.cpp b/kig/misc/kigtransform.cpp new file mode 100644 index 00000000..a297ed6e --- /dev/null +++ b/kig/misc/kigtransform.cpp @@ -0,0 +1,810 @@ +/** + This file is part of Kig, a KDE program for Interactive Geometry... + Copyright (C) 2002 Maurizio Paolini <paolini@dmf.unicatt.it> + Copyright (C) 2003 Dominique Devriese <devriese@kde.org> + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 + USA +**/ + +#include "kigtransform.h" + +#include "kignumerics.h" +#include "common.h" + +#include <cmath> + +#include <klocale.h> +#include <kdebug.h> + +// Transformation getProjectiveTransformation ( int argsnum, +// Object *transforms[], bool& valid ) +// { +// valid = true; + +// assert ( argsnum > 0 ); +// int argn = 0; +// Object* transform = transforms[argn++]; +// if (transform->toVector()) +// { +// // translation +// assert (argn == argsnum); +// Vector* v = transform->toVector(); +// Coordinate dir = v->getDir(); +// return Transformation::translation( dir ); +// } + +// if (transform->toPoint()) +// { +// // point reflection ( or is point symmetry the correct term ? ) +// assert (argn == argsnum); +// Point* p = transform->toPoint(); +// return Transformation::pointReflection( p->getCoord() ); +// } + +// if (transform->toLine()) +// { +// // line reflection ( or is it line symmetry ? ) +// Line* line = transform->toLine(); +// assert (argn == argsnum); +// return Transformation::lineReflection( line->lineData() ); +// } + +// if (transform->toRay()) +// { +// // domi: sorry, but what kind of transformation does this do ? +// // i'm guessing it's some sort of rotation, but i'm not +// // really sure.. +// Ray* line = transform->toRay(); +// Coordinate d = line->direction().normalize(); +// Coordinate t = line->p1(); +// double alpha = 0.1*M_PI/2; // a small angle for the DrawPrelim +// if (argn < argsnum) +// { +// Angle* angle = transforms[argn++]->toAngle(); +// alpha = angle->size(); +// } +// assert (argn == argsnum); +// return Transformation::projectiveRotation( alpha, d, t ); +// } + +// if (transform->toAngle()) +// { +// // rotation.. +// Coordinate center = Coordinate( 0., 0. ); +// if (argn < argsnum) +// { +// Object* arg = transforms[argn++]; +// assert (arg->toPoint()); +// center = arg->toPoint()->getCoord(); +// } +// Angle* angle = transform->toAngle(); +// double alpha = angle->size(); + +// assert (argn == argsnum); + +// return Transformation::rotation( alpha, center ); +// } + +// if (transform->toSegment()) // this is a scaling +// { +// Segment* segment = transform->toSegment(); +// Coordinate p = segment->p2() - segment->p1(); +// double s = p.length(); +// if (argn < argsnum) +// { +// Object* arg = transforms[argn++]; +// if (arg->toSegment()) // s is the length of the first segment +// // divided by the length of the second.. +// { +// Segment* segment = arg->toSegment(); +// Coordinate p = segment->p2() - segment->p1(); +// s /= p.length(); +// if (argn < argsnum) arg = transforms[argn++]; +// } +// if (arg->toPoint()) // scaling w.r. to a point +// { +// Point* p = arg->toPoint(); +// assert (argn == argsnum); +// return Transformation::scaling( s, p->getCoord() ); +// } +// if (arg->toLine()) // scaling w.r. to a line +// { +// Line* line = arg->toLine(); +// assert( argn == argsnum ); +// return Transformation::scaling( s, line->lineData() ); +// } +// } + +// return Transformation::scaling( s, Coordinate( 0., 0. ) ); +// } + +// valid = false; +// return Transformation::identity(); +// } + +// tWantArgsResult WantTransformation ( Objects::const_iterator& i, +// const Objects& os ) +// { +// Object* o = *i++; +// if (o->toVector()) return tComplete; +// if (o->toPoint()) return tComplete; +// if (o->toLine()) return tComplete; +// if (o->toAngle()) +// { +// if ( i == os.end() ) return tNotComplete; +// o = *i++; +// if (o->toPoint()) return tComplete; +// if (o->toLine()) return tComplete; +// return tNotGood; +// } +// if (o->toRay()) +// { +// if ( i == os.end() ) return tNotComplete; +// o = *i++; +// if (o->toAngle()) return tComplete; +// return tNotGood; +// } +// if (o->toSegment()) +// { +// if ( i == os.end() ) return tNotComplete; +// o = *i++; +// if ( o->toSegment() ) +// { +// if ( i == os.end() ) return tNotComplete; +// o = *i++; +// } +// if (o->toPoint()) return tComplete; +// if (o->toLine()) return tComplete; +// return tNotGood; +// } +// return tNotGood; +// } + +// QString getTransformMessage ( const Objects& os, const Object *o ) +// { +// int size = os.size(); +// switch (size) +// { +// case 1: +// if (o->toVector()) return i18n("translate by this vector"); +// if (o->toPoint()) return i18n("central symmetry by this point. You" +// " can obtain different transformations by clicking on lines (mirror)," +// " vectors (translation), angles (rotation), segments (scaling) and rays" +// " (projective transformation)"); +// if (o->toLine()) return i18n("reflect in this line"); +// if (o->toAngle()) return i18n("rotate by this angle"); +// if (o->toSegment()) return i18n("scale using the length of this vector"); +// if (o->toRay()) return i18n("a projective transformation in the direction" +// " indicated by this ray, it is a rotation in the projective plane" +// " about a point at infinity"); +// return i18n("Use this transformation"); + +// case 2: // we ask for the first parameter of the transformation +// case 3: +// if (os[1]->toAngle()) +// { +// if (o->toPoint()) return i18n("about this point"); +// assert (false); +// } +// if (os[1]->toSegment()) +// { +// if (o->toSegment()) +// return i18n("relative to the length of this other vector"); +// if (o->toPoint()) +// return i18n("about this point"); +// if (o->toLine()) +// return i18n("about this line"); +// } +// if (os[1]->toRay()) +// { +// if (o->toAngle()) return i18n("rotate by this angle in the projective" +// " plane"); +// } +// return i18n("Using this object"); + +// default: assert(false); +// } + +// return i18n("Use this transformation"); +// } + + +/* domi: not necessary anymore, homotheticness is kept as a bool in + * the Transformation class.. + * keeping it here, in case a need for it arises some time in the + * future... + * decide if the given transformation is homotetic + */ +// bool isHomoteticTransformation ( double transformation[3][3] ) +// { +// if (transformation[0][1] != 0 || transformation[0][2] != 0) return (false); +// // test the orthogonality of the matrix 2x2 of second and third rows +// // and columns +// if (fabs(fabs(transformation[1][1]) - +// fabs(transformation[2][2])) > 1e-8) return (false); +// if (fabs(fabs(transformation[1][2]) - +// fabs(transformation[2][1])) > 1e-8) return (false); + +// return transformation[1][2] * transformation[2][1] * +// transformation[1][1] * transformation[2][2] <= 0.; +// } + +const Transformation Transformation::identity() +{ + Transformation ret; + for ( int i = 0; i < 3; ++i ) + for ( int j = 0; j < 3; ++j ) + ret.mdata[i][j] = ( i == j ? 1 : 0 ); + ret.mIsHomothety = ret.mIsAffine = true; + return ret; +} + +const Transformation Transformation::scalingOverPoint( double factor, const Coordinate& center ) +{ + Transformation ret; + for ( int i = 0; i < 3; ++i ) + for ( int j = 0; j < 3; ++j ) + ret.mdata[i][j] = ( i == j ? factor : 0 ); + ret.mdata[0][0] = 1; + ret.mdata[1][0] = center.x - factor * center.x; + ret.mdata[2][0] = center.y - factor * center.y; + ret.mIsHomothety = ret.mIsAffine = true; + return ret; +} + +const Transformation Transformation::translation( const Coordinate& c ) +{ + Transformation ret = identity(); + ret.mdata[1][0] = c.x; + ret.mdata[2][0] = c.y; + + // this is already set in the identity() constructor, but just for + // clarity.. + ret.mIsHomothety = ret.mIsAffine = true; + return ret; +} + +const Transformation Transformation::pointReflection( const Coordinate& c ) +{ + Transformation ret = scalingOverPoint( -1, c ); + ret.mIsHomothety = ret.mIsAffine = true; + return ret; +} + +const Transformation operator*( const Transformation& a, const Transformation& b ) +{ + // just multiply the two matrices.. + Transformation ret; + + for ( int i = 0; i < 3; ++i ) + for ( int j = 0; j < 3; ++j ) + { + ret.mdata[i][j] = 0; + for ( int k = 0; k < 3; ++k ) + ret.mdata[i][j] += a.mdata[i][k] * b.mdata[k][j]; + }; + + // combination of two homotheties is a homothety.. + + ret.mIsHomothety = a.mIsHomothety && b.mIsHomothety; + + // combination of two affinities is affine.. + + ret.mIsAffine = a.mIsAffine && b.mIsAffine; + + return ret; +} + +const Transformation Transformation::lineReflection( const LineData& l ) +{ + Transformation ret = scalingOverLine( -1, l ); + // a reflection is a homothety... + ret.mIsHomothety = ret.mIsAffine = true; + return ret; +} + +const Transformation Transformation::scalingOverLine( double factor, const LineData& l ) +{ + Transformation ret = identity(); + + Coordinate a = l.a; + Coordinate d = l.dir(); + double dirnormsq = d.squareLength(); + ret.mdata[1][1] = (d.x*d.x + factor*d.y*d.y)/dirnormsq; + ret.mdata[2][2] = (d.y*d.y + factor*d.x*d.x)/dirnormsq; + ret.mdata[1][2] = ret.mdata[2][1] = (d.x*d.y - factor*d.x*d.y)/dirnormsq; + + ret.mdata[1][0] = a.x - ret.mdata[1][1]*a.x - ret.mdata[1][2]*a.y; + ret.mdata[2][0] = a.y - ret.mdata[2][1]*a.x - ret.mdata[2][2]*a.y; + + // domi: is 1e-8 a good value ? + ret.mIsHomothety = ( fabs( factor - 1 ) < 1e-8 || fabs ( factor + 1 ) < 1e-8 ); + ret.mIsAffine = true; + return ret; +} + +const Transformation Transformation::harmonicHomology( + const Coordinate& center, const LineData& axis ) +{ + // this is a well known projective transformation. We find it by first + // computing the homogeneous equation of the axis ax + by + cz = 0 + // then a straightforward computation shows that the 3x3 matrix describing + // the transformation is of the form: + // + // (r . C) Id - 2 (C tensor r) + // + // where r = [c, a, b], C = [1, Cx, Cy], Cx and Cy are the coordinates of + // the center, '.' denotes the scalar product, Id is the identity matrix, + // 'tensor' is the tensor product producing a 3x3 matrix. + // + // note: here we decide to use coordinate '0' in place of the third coordinate + // in homogeneous notation; e.g. C = [1, cx, cy] + + Coordinate pointa = axis.a; + Coordinate pointb = axis.b; + + double a = pointa.y - pointb.y; + double b = pointb.x - pointa.x; + double c = pointa.x*pointb.y - pointa.y*pointb.x; + + double cx = center.x; + double cy = center.y; + + double scalprod = a*cx + b*cy + c; + scalprod *= 0.5; + Transformation ret; + + ret.mdata[0][0] = c - scalprod; + ret.mdata[0][1] = a; + ret.mdata[0][2] = b; + + ret.mdata[1][0] = c*cx; + ret.mdata[1][1] = a*cx - scalprod; + ret.mdata[1][2] = b*cx; + + ret.mdata[2][0] = c*cy; + ret.mdata[2][1] = a*cy; + ret.mdata[2][2] = b*cy - scalprod; + + ret.mIsHomothety = ret.mIsAffine = false; + return ret; +} + +const Transformation Transformation::affinityGI3P( + const std::vector<Coordinate>& FromPoints, + const std::vector<Coordinate>& ToPoints, + bool& valid ) +{ + // construct the (generically) unique affinity that transforms 3 given + // point into 3 other given points; i.e. it depends on the coordinates of + // a total of 6 points. This actually amounts in solving a 6x6 linear + // system to find the entries of a 2x2 linear transformation matrix T + // and of a translation vector t. + // If Pi denotes one of the starting points and Qi the corresponding + // final position we actually have to solve: Qi = t + T Pi, for i=1,2,3 + // (each one is a vector equation, so that it really gives 2 equations). + // In our context T and t are used to build a 3x3 projective transformation + // as follows: + // + // [ 1 0 0 ] + // [ t1 T11 T12 ] + // [ t2 T21 T22 ] + // + // In order to take advantage of the two functions "GaussianElimination" + // and "BackwardSubstitution", which are specifically aimed at solving + // homogeneous underdetermined linear systems, we just add a further + // unknown m and solve for t + T Pi - m Qi = 0. Since our functions + // returns a nonzero solution we shall have a nonzero 'm' in the end and + // can build the 3x3 matrix as follows: + // + // [ m 0 0 ] + // [ t1 T11 T12 ] + // [ t2 T21 T22 ] + // + // we order the unknowns as follows: m, t1, t2, T11, T12, T21, T22 + + double row0[7], row1[7], row2[7], row3[7], row4[7], row5[7]; + + double *matrix[6] = {row0, row1, row2, row3, row4, row5}; + + double solution[7]; + int scambio[7]; + + assert (FromPoints.size() == 3); + assert (ToPoints.size() == 3); + + // fill in the matrix elements + for ( int i = 0; i < 6; i++ ) + { + for ( int j = 0; j < 7; j++ ) + { + matrix[i][j] = 0.0; + } + } + + for ( int i = 0; i < 3; i++ ) + { + Coordinate p = FromPoints[i]; + Coordinate q = ToPoints[i]; + matrix[i][0] = -q.x; + matrix[i][1] = 1.0; + matrix[i][3] = p.x; + matrix[i][4] = p.y; + matrix[i+3][0] = -q.y; + matrix[i+3][2] = 1.0; + matrix[i+3][5] = p.x; + matrix[i+3][6] = p.y; + } + + Transformation ret; + valid = true; + if ( ! GaussianElimination( matrix, 6, 7, scambio ) ) + { valid = false; return ret; } + + // fine della fase di eliminazione + BackwardSubstitution( matrix, 6, 7, scambio, solution ); + + // now we can build the 3x3 transformation matrix; remember that + // unknown 0 is the multiplicator 'm' + + ret.mdata[0][0] = solution[0]; + ret.mdata[0][1] = ret.mdata[0][2] = 0.0; + ret.mdata[1][0] = solution[1]; + ret.mdata[2][0] = solution[2]; + ret.mdata[1][1] = solution[3]; + ret.mdata[1][2] = solution[4]; + ret.mdata[2][1] = solution[5]; + ret.mdata[2][2] = solution[6]; + + ret.mIsHomothety = false; + ret.mIsAffine = true; + return ret; +} + +const Transformation Transformation::projectivityGI4P( + const std::vector<Coordinate>& FromPoints, + const std::vector<Coordinate>& ToPoints, + bool& valid ) +{ + // construct the (generically) unique projectivity that transforms 4 given + // point into 4 other given points; i.e. it depends on the coordinates of + // a total of 8 points. This actually amounts in solving an underdetermined + // homogeneous linear system. + + double + row0[13], row1[13], row2[13], row3[13], row4[13], row5[13], row6[13], row7[13], + row8[13], row9[13], row10[13], row11[13]; + + double *matrix[12] = {row0, row1, row2, row3, row4, row5, row6, row7, + row8, row9, row10, row11}; + + double solution[13]; + int scambio[13]; + + assert (FromPoints.size() == 4); + assert (ToPoints.size() == 4); + + // fill in the matrix elements + for ( int i = 0; i < 12; i++ ) + { + for ( int j = 0; j < 13; j++ ) + { + matrix[i][j] = 0.0; + } + } + + for ( int i = 0; i < 4; i++ ) + { + Coordinate p = FromPoints[i]; + Coordinate q = ToPoints[i]; + matrix[i][0] = matrix[4+i][3] = matrix[8+i][6] = 1.0; + matrix[i][1] = matrix[4+i][4] = matrix[8+i][7] = p.x; + matrix[i][2] = matrix[4+i][5] = matrix[8+i][8] = p.y; + matrix[i][9+i] = -1.0; + matrix[4+i][9+i] = -q.x; + matrix[8+i][9+i] = -q.y; + } + + Transformation ret; + valid = true; + if ( ! GaussianElimination( matrix, 12, 13, scambio ) ) + { valid = false; return ret; } + + // fine della fase di eliminazione + BackwardSubstitution( matrix, 12, 13, scambio, solution ); + + // now we can build the 3x3 transformation matrix; remember that + // unknowns from 9 to 13 are just multiplicators that we don't need here + + int k = 0; + for ( int i = 0; i < 3; i++ ) + { + for ( int j = 0; j < 3; j++ ) + { + ret.mdata[i][j] = solution[k++]; + } + } + + ret.mIsHomothety = ret.mIsAffine = false; + return ret; +} + +const Transformation Transformation::castShadow( + const Coordinate& lightsrc, const LineData& l ) +{ + // first deal with the line l, I need to find an appropriate reflection + // that transforms l onto the x-axis + + Coordinate d = l.dir(); + Coordinate a = l.a; + double k = d.length(); + if ( d.x < 0 ) k *= -1; // for numerical stability + Coordinate w = d + Coordinate( k, 0 ); + // w /= w.length(); + // w defines a Householder transformation, but we don't need to normalize + // it here. + // warning: this w is the orthogonal of the w of the textbooks! + // this is fine for us since in this way it indicates the line direction + Coordinate ra = Coordinate ( a.x + w.y*a.y/(2*w.x), a.y/2 ); + Transformation sym = lineReflection ( LineData( ra, ra + w ) ); + + // in the new coordinates the line is the x-axis + // I must transform the point + + Coordinate modlightsrc = sym.apply ( lightsrc ); + Transformation ret = identity(); + // parameter t indicates the distance of the light source from + // the plane of the drawing. A negative value means that the light + // source is behind the plane. + double t = -1.0; + // double t = -modlightsrc.y; <-- this gives the old transformation! + double e = modlightsrc.y - t; + ret.mdata[0][0] = e; + ret.mdata[0][2] = -1; + ret.mdata[1][1] = e; + ret.mdata[1][2] = -modlightsrc.x; + ret.mdata[2][2] = -t; + + ret.mIsHomothety = ret.mIsAffine = false; + return sym*ret*sym; +// return translation( t )*ret*translation( -t ); +} + +const Transformation Transformation::projectiveRotation( + double alpha, const Coordinate& d, const Coordinate& t ) +{ + Transformation ret; + double cosalpha = cos( alpha ); + double sinalpha = sin( alpha ); + ret.mdata[0][0] = cosalpha; + ret.mdata[1][1] = cosalpha*d.x*d.x + d.y*d.y; + ret.mdata[0][1] = -sinalpha*d.x; + ret.mdata[1][0] = sinalpha*d.x; + ret.mdata[0][2] = -sinalpha*d.y; + ret.mdata[2][0] = sinalpha*d.y; + ret.mdata[1][2] = cosalpha*d.x*d.y - d.x*d.y; + ret.mdata[2][1] = cosalpha*d.x*d.y - d.x*d.y; + ret.mdata[2][2] = cosalpha*d.y*d.y + d.x*d.x; + + ret.mIsHomothety = ret.mIsAffine = false; + return translation( t )*ret*translation( -t ); +} + +const Coordinate Transformation::apply( const double x0, + const double x1, + const double x2) const +{ + double phom[3] = {x0, x1, x2}; + double rhom[3] = {0., 0., 0.}; + + + for (int i = 0; i < 3; i++) + { + for (int j = 0; j < 3; j++) + { + rhom[i] += mdata[i][j]*phom[j]; + } + } + + if (rhom[0] == 0.) + return Coordinate::invalidCoord(); + + return Coordinate (rhom[1]/rhom[0], rhom[2]/rhom[0]); +} + +const Coordinate Transformation::apply( const Coordinate& p ) const +{ + return apply( 1., p.x, p.y ); +// double phom[3] = {1., p.x, p.y}; +// double rhom[3] = {0., 0., 0.}; +// +// for (int i = 0; i < 3; i++) +// { +// for (int j = 0; j < 3; j++) +// { +// rhom[i] += mdata[i][j]*phom[j]; +// } +// } +// +// if (rhom[0] == 0.) +// return Coordinate::invalidCoord(); +// +// return Coordinate (rhom[1]/rhom[0], rhom[2]/rhom[0]); +} + +const Coordinate Transformation::apply0( const Coordinate& p ) const +{ + return apply( 0., p.x, p.y ); +} + +const Transformation Transformation::rotation( double alpha, const Coordinate& center ) +{ + Transformation ret = identity(); + + double x = center.x; + double y = center.y; + + double cosalpha = cos( alpha ); + double sinalpha = sin( alpha ); + + ret.mdata[1][1] = ret.mdata[2][2] = cosalpha; + ret.mdata[1][2] = -sinalpha; + ret.mdata[2][1] = sinalpha; + ret.mdata[1][0] = x - ret.mdata[1][1]*x - ret.mdata[1][2]*y; + ret.mdata[2][0] = y - ret.mdata[2][1]*x - ret.mdata[2][2]*y; + + // this is already set in the identity() constructor, but just for + // clarity.. + ret.mIsHomothety = ret.mIsAffine = true; + + return ret; +} + +bool Transformation::isHomothetic() const +{ + return mIsHomothety; +} + +bool Transformation::isAffine() const +{ + return mIsAffine; +} + +/* + *mp: + * this function has the property that it changes sign if computed + * on two points that lie on either sides with respect to the critical + * line (this is the line that goes to the line at infinity). + * For affine transformations the result has always the same sign. + * NOTE: the result is *not* invariant under rescaling of all elements + * of the transformation matrix. + * The typical use is to determine whether a segment is transformed + * into a segment or a couple of half-lines. + */ + +double Transformation::getProjectiveIndicator( const Coordinate& c ) const +{ + return mdata[0][0] + mdata[0][1]*c.x + mdata[0][2]*c.y; +} + +// assuming that this is an affine transformation, return its +// determinant. What is really important here is just the sign +// of the determinant. +double Transformation::getAffineDeterminant() const +{ + return mdata[1][1]*mdata[2][2] - mdata[1][2]*mdata[2][1]; +} + +// this assumes that the 2x2 affine part of the matrix is of the +// form [ cos a, sin a; -sin a, cos a] or a multiple +double Transformation::getRotationAngle() const +{ + return atan2( mdata[1][2], mdata[1][1] ); +} + +const Coordinate Transformation::apply2by2only( const Coordinate& p ) const +{ + double x = p.x; + double y = p.y; + double nx = mdata[1][1]*x + mdata[1][2]*y; + double ny = mdata[2][1]*x + mdata[2][2]*y; + return Coordinate( nx, ny ); +} + +double Transformation::data( int r, int c ) const +{ + return mdata[r][c]; +} + +const Transformation Transformation::inverse( bool& valid ) const +{ + Transformation ret; + + valid = Invert3by3matrix( mdata, ret.mdata ); + + // the inverse of a homothety is a homothety, same for affinities.. + ret.mIsHomothety = mIsHomothety; + ret.mIsAffine = mIsAffine; + + return ret; +} + +Transformation::Transformation() +{ + // this is the constructor used by the static Transformation + // creation functions, so mIsHomothety is in general false + mIsHomothety = mIsAffine = false; + for ( int i = 0; i < 3; ++i ) + for ( int j = 0; j < 3; ++j ) + mdata[i][j] = ( i == j ) ? 1 : 0; +} + +Transformation::~Transformation() +{ +} + +double Transformation::apply( double length ) const +{ + assert( isHomothetic() ); + double det = mdata[1][1]*mdata[2][2] - + mdata[1][2]*mdata[2][1]; + return sqrt( fabs( det ) ) * length; +} + +Transformation::Transformation( double data[3][3], bool ishomothety ) + : mIsHomothety( ishomothety ) +{ + for ( int i = 0; i < 3; ++i ) + for ( int j = 0; j < 3; ++j ) + mdata[i][j] = data[i][j]; + + //mp: a test for affinity is used to initialize mIsAffine... + mIsAffine = false; + if ( fabs(mdata[0][1]) + fabs(mdata[0][2]) < 1e-8 * fabs(mdata[0][0]) ) + mIsAffine = true; +} + +bool operator==( const Transformation& lhs, const Transformation& rhs ) +{ + for ( int i = 0; i < 3; ++i ) + for ( int j = 0; j < 3; ++j ) + if ( lhs.data( i, j ) != rhs.data( i, j ) ) + return false; + return true; +} + +const Transformation Transformation::similitude( + const Coordinate& center, double theta, double factor ) +{ + //kdDebug() << k_funcinfo << "theta: " << theta << " factor: " << factor << endl; + Transformation ret; + ret.mIsHomothety = true; + double costheta = cos( theta ); + double sintheta = sin( theta ); + ret.mdata[0][0] = 1; + ret.mdata[0][1] = 0; + ret.mdata[0][2] = 0; + ret.mdata[1][0] = ( 1 - factor*costheta )*center.x + factor*sintheta*center.y; + ret.mdata[1][1] = factor*costheta; + ret.mdata[1][2] = -factor*sintheta; + ret.mdata[2][0] = -factor*sintheta*center.x + ( 1 - factor*costheta )*center.y; + ret.mdata[2][1] = factor*sintheta; + ret.mdata[2][2] = factor*costheta; + // fails for factor == infinity + //assert( ( ret.apply( center ) - center ).length() < 1e-5 ); + ret.mIsHomothety = ret.mIsAffine = true; + return ret; +} |