kig/misc/kigtransform.cpp

   1 /**
   2  This file is part of Kig, a KDE program for Interactive Geometry...
   3  Copyright (C) 2002  Maurizio Paolini <paolini@dmf.unicatt.it>
   4  Copyright (C) 2003  Dominique Devriese <devriese@kde.org>
   5
   6  This program is free software; you can redistribute it and/or modify
   7  it under the terms of the GNU General Public License as published by
   8  the Free Software Foundation; either version 2 of the License, or
   9  (at your option) any later version.
  10
  11  This program is distributed in the hope that it will be useful,
  12  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  GNU General Public License for more details.
  15
  16  You should have received a copy of the GNU General Public License
  17  along with this program; if not, write to the Free Software
  18  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
  19  USA
  20 **/
  21
  22 #include "kigtransform.h"
  23
  24 #include "kignumerics.h"
  25 #include "common.h"
  26
  27 #include <cmath>
  28
  29 #include <klocale.h>
  30 #include <kdebug.h>
  31
  32 // Transformation getProjectiveTransformation ( int argsnum,
  33 //     Object *transforms[], bool& valid )
  34 // {
  35 //   valid = true;
  36
  37 //   assert ( argsnum > 0 );
  38 //   int argn = 0;
  39 //   Object* transform = transforms[argn++];
  40 //   if (transform->toVector())
  41 //   {
  42 //     // translation
  43 //     assert (argn == argsnum);
  44 //     Vector* v = transform->toVector();
  45 //     Coordinate dir = v->getDir();
  46 //     return Transformation::translation( dir );
  47 //   }
  48
  49 //   if (transform->toPoint())
  50 //   {
  51 //     // point reflection ( or is point symmetry the correct term ? )
  52 //     assert (argn == argsnum);
  53 //     Point* p = transform->toPoint();
  54 //     return Transformation::pointReflection( p->getCoord() );
  55 //   }
  56
  57 //   if (transform->toLine())
  58 //   {
  59 //     // line reflection ( or is it line symmetry ? )
  60 //     Line* line = transform->toLine();
  61 //     assert (argn == argsnum);
  62 //     return Transformation::lineReflection( line->lineData() );
  63 //   }
  64
  65 //   if (transform->toRay())
  66 //   {
  67 //     // domi: sorry, but what kind of transformation does this do ?
  68 //     //       i'm guessing it's some sort of rotation, but i'm not
  69 //     //       really sure..
  70 //      Ray* line = transform->toRay();
  71 //      Coordinate d = line->direction().normalize();
  72 //      Coordinate t = line->p1();
  73 //      double alpha = 0.1*M_PI/2;  // a small angle for the DrawPrelim
  74 //      if (argn < argsnum)
  75 //      {
  76 //        Angle* angle = transforms[argn++]->toAngle();
  77 //        alpha = angle->size();
  78 //      }
  79 //      assert (argn == argsnum);
  80 //      return Transformation::projectiveRotation( alpha, d, t );
  81 //   }
  82
  83 //   if (transform->toAngle())
  84 //   {
  85 //     // rotation..
  86 //     Coordinate center = Coordinate( 0., 0. );
  87 //     if (argn < argsnum)
  88 //     {
  89 //       Object* arg = transforms[argn++];
  90 //       assert (arg->toPoint());
  91 //       center = arg->toPoint()->getCoord();
  92 //     }
  93 //     Angle* angle = transform->toAngle();
  94 //     double alpha = angle->size();
  95
  96 //     assert (argn == argsnum);
  97
  98 //     return Transformation::rotation( alpha, center );
  99 //   }
 100
 101 //   if (transform->toSegment())  // this is a scaling
 102 //   {
 103 //     Segment* segment = transform->toSegment();
 104 //     Coordinate p = segment->p2() - segment->p1();
 105 //     double s = p.length();
 106 //     if (argn < argsnum)
 107 //     {
 108 //       Object* arg = transforms[argn++];
 109 //       if (arg->toSegment()) // s is the length of the first segment
 110 //                             // divided by the length of the second..
 111 //       {
 112 //         Segment* segment = arg->toSegment();
 113 //         Coordinate p = segment->p2() - segment->p1();
 114 //         s /= p.length();
 115 //         if (argn < argsnum) arg = transforms[argn++];
 116 //       }
 117 //       if (arg->toPoint())   // scaling w.r. to a point
 118 //       {
 119 //         Point* p = arg->toPoint();
 120 //         assert (argn == argsnum);
 121 //         return Transformation::scaling( s, p->getCoord() );
 122 //       }
 123 //       if (arg->toLine())  // scaling w.r. to a line
 124 //       {
 125 //         Line* line = arg->toLine();
 126 //         assert( argn == argsnum );
 127 //         return Transformation::scaling( s, line->lineData() );
 128 //       }
 129 //     }
 130
 131 //     return Transformation::scaling( s, Coordinate( 0., 0. ) );
 132 //   }
 133
 134 //   valid = false;
 135 //   return Transformation::identity();
 136 // }
 137
 138 // tWantArgsResult WantTransformation ( Objects::const_iterator& i,
 139 //       const Objects& os )
 140 // {
 141 //   Object* o = *i++;
 142 //   if (o->toVector()) return tComplete;
 143 //   if (o->toPoint()) return tComplete;
 144 //   if (o->toLine()) return tComplete;
 145 //   if (o->toAngle())
 146 //   {
 147 //     if ( i == os.end() ) return tNotComplete;
 148 //     o = *i++;
 149 //     if (o->toPoint()) return tComplete;
 150 //     if (o->toLine()) return tComplete;
 151 //     return tNotGood;
 152 //   }
 153 //   if (o->toRay())
 154 //   {
 155 //     if ( i == os.end() ) return tNotComplete;
 156 //     o = *i++;
 157 //     if (o->toAngle()) return tComplete;
 158 //     return tNotGood;
 159 //   }
 160 //   if (o->toSegment())
 161 //   {
 162 //     if ( i == os.end() ) return tNotComplete;
 163 //     o = *i++;
 164 //     if ( o->toSegment() )
 165 //     {
 166 //       if ( i == os.end() ) return tNotComplete;
 167 //       o = *i++;
 168 //     }
 169 //     if (o->toPoint()) return tComplete;
 170 //     if (o->toLine()) return tComplete;
 171 //     return tNotGood;
 172 //   }
 173 //   return tNotGood;
 174 // }
 175
 176 // QString getTransformMessage ( const Objects& os, const Object *o )
 177 // {
 178 //   int size = os.size();
 179 //   switch (size)
 180 //   {
 181 //     case 1:
 182 //     if (o->toVector()) return i18n("translate by this vector");
 183 //     if (o->toPoint()) return i18n("central symmetry by this point. You"
 184 //      " can obtain different transformations by clicking on lines (mirror),"
 185 //      " vectors (translation), angles (rotation), segments (scaling) and rays"
 186 //      " (projective transformation)");
 187 //     if (o->toLine()) return i18n("reflect in this line");
 188 //     if (o->toAngle()) return i18n("rotate by this angle");
 189 //     if (o->toSegment()) return i18n("scale using the length of this vector");
 190 //     if (o->toRay()) return i18n("a projective transformation in the direction"
 191 //      " indicated by this ray, it is a rotation in the projective plane"
 192 //      " about a point at infinity");
 193 //     return i18n("Use this transformation");
 194
 195 //     case 2:   // we ask for the first parameter of the transformation
 196 //     case 3:
 197 //     if (os[1]->toAngle())
 198 //     {
 199 //       if (o->toPoint()) return i18n("about this point");
 200 //       assert (false);
 201 //     }
 202 //     if (os[1]->toSegment())
 203 //     {
 204 //       if (o->toSegment())
 205 //         return i18n("relative to the length of this other vector");
 206 //       if (o->toPoint())
 207 //         return i18n("about this point");
 208 //       if (o->toLine())
 209 //         return i18n("about this line");
 210 //     }
 211 //     if (os[1]->toRay())
 212 //     {
 213 //       if (o->toAngle()) return i18n("rotate by this angle in the projective"
 214 //        " plane");
 215 //     }
 216 //     return i18n("Using this object");
 217
 218 //     default: assert(false);
 219 //   }
 220
 221 //   return i18n("Use this transformation");
 222 // }
 223
 224
 225 /* domi: not necessary anymore, homotheticness is kept as a bool in
 226  *       the Transformation class..
 227  * keeping it here, in case a need for it arises some time in the
 228  * future...
 229  * decide if the given transformation is homotetic
 230  */
 231 // bool isHomoteticTransformation ( double transformation[3][3] )
 232 // {
 233 //   if (transformation[0][1] != 0 || transformation[0][2] != 0) return (false);
 234 //   // test the orthogonality of the matrix 2x2 of second and third rows
 235 //   // and columns
 236 //   if (fabs(fabs(transformation[1][1]) -
 237 //             fabs(transformation[2][2])) > 1e-8) return (false);
 238 //   if (fabs(fabs(transformation[1][2]) -
 239 //             fabs(transformation[2][1])) > 1e-8) return (false);
 240
 241 //   return transformation[1][2] * transformation[2][1] *
 242 //          transformation[1][1] * transformation[2][2] <= 0.;
 243 // }
 244
 245 const Transformation Transformation::identity()
 246 {
 247   Transformation ret;
 248   for ( int i = 0; i < 3; ++i )
 249     for ( int j = 0; j < 3; ++j )
 250       ret.mdata[i][j] = ( i == j ? 1 : 0 );
 251   ret.mIsHomothety = ret.mIsAffine = true;
 252   return ret;
 253 }
 254
 255 const Transformation Transformation::scalingOverPoint( double factor, const Coordinate& center )
 256 {
 257   Transformation ret;
 258   for ( int i = 0; i < 3; ++i )
 259     for ( int j = 0; j < 3; ++j )
 260       ret.mdata[i][j] = ( i == j ? factor : 0 );
 261   ret.mdata[0][0] = 1;
 262   ret.mdata[1][0] = center.x - factor * center.x;
 263   ret.mdata[2][0] = center.y - factor * center.y;
 264   ret.mIsHomothety = ret.mIsAffine = true;
 265   return ret;
 266 }
 267
 268 const Transformation Transformation::translation( const Coordinate& c )
 269 {
 270   Transformation ret = identity();
 271   ret.mdata[1][0] = c.x;
 272   ret.mdata[2][0] = c.y;
 273
 274   // this is already set in the identity() constructor, but just for
 275   // clarity..
 276   ret.mIsHomothety = ret.mIsAffine = true;
 277   return ret;
 278 }
 279
 280 const Transformation Transformation::pointReflection( const Coordinate& c )
 281 {
 282   Transformation ret = scalingOverPoint( -1, c );
 283   ret.mIsHomothety = ret.mIsAffine = true;
 284   return ret;
 285 }
 286
 287 const Transformation operator*( const Transformation& a, const Transformation& b )
 288 {
 289   // just multiply the two matrices..
 290   Transformation ret;
 291
 292   for ( int i = 0; i < 3; ++i )
 293     for ( int j = 0; j < 3; ++j )
 294     {
 295       ret.mdata[i][j] = 0;
 296       for ( int k = 0; k < 3; ++k )
 297         ret.mdata[i][j] += a.mdata[i][k] * b.mdata[k][j];
 298     };
 299
 300   // combination of two homotheties is a homothety..
 301
 302   ret.mIsHomothety = a.mIsHomothety && b.mIsHomothety;
 303
 304   // combination of two affinities is affine..
 305
 306   ret.mIsAffine = a.mIsAffine && b.mIsAffine;
 307
 308   return ret;
 309 }
 310
 311 const Transformation Transformation::lineReflection( const LineData& l )
 312 {
 313   Transformation ret = scalingOverLine( -1, l );
 314   // a reflection is a homothety...
 315   ret.mIsHomothety = ret.mIsAffine = true;
 316   return ret;
 317 }
 318
 319 const Transformation Transformation::scalingOverLine( double factor, const LineData& l )
 320 {
 321   Transformation ret = identity();
 322
 323   Coordinate a = l.a;
 324   Coordinate d = l.dir();
 325   double dirnormsq = d.squareLength();
 326   ret.mdata[1][1] = (d.x*d.x + factor*d.y*d.y)/dirnormsq;
 327   ret.mdata[2][2] = (d.y*d.y + factor*d.x*d.x)/dirnormsq;
 328   ret.mdata[1][2] = ret.mdata[2][1] = (d.x*d.y - factor*d.x*d.y)/dirnormsq;
 329
 330   ret.mdata[1][0] = a.x - ret.mdata[1][1]*a.x - ret.mdata[1][2]*a.y;
 331   ret.mdata[2][0] = a.y - ret.mdata[2][1]*a.x - ret.mdata[2][2]*a.y;
 332
 333   // domi: is 1e-8 a good value ?
 334   ret.mIsHomothety = ( fabs( factor - 1 ) < 1e-8 || fabs ( factor + 1 ) < 1e-8 );
 335   ret.mIsAffine = true;
 336   return ret;
 337 }
 338
 339 const Transformation Transformation::harmonicHomology(
 340   const Coordinate& center, const LineData& axis )
 341 {
 342   // this is a well known projective transformation.  We find it by first
 343   // computing the homogeneous equation of the axis ax + by + cz = 0
 344   // then a straightforward computation shows that the 3x3 matrix describing
 345   // the transformation is of the form:
 346   //
 347   // (r . C) Id - 2 (C tensor r)
 348   //
 349   // where r = [c, a, b], C = [1, Cx, Cy], Cx and Cy are the coordinates of
 350   // the center, '.' denotes the scalar product, Id is the identity matrix,
 351   // 'tensor' is the tensor product producing a 3x3 matrix.
 352   //
 353   // note: here we decide to use coordinate '0' in place of the third coordinate
 354   // in homogeneous notation; e.g. C = [1, cx, cy]
 355
 356   Coordinate pointa = axis.a;
 357   Coordinate pointb = axis.b;
 358
 359   double a = pointa.y - pointb.y;
 360   double b = pointb.x - pointa.x;
 361   double c = pointa.x*pointb.y - pointa.y*pointb.x;
 362
 363   double cx = center.x;
 364   double cy = center.y;
 365
 366   double scalprod = a*cx + b*cy + c;
 367   scalprod *= 0.5;
 368   Transformation ret;
 369
 370   ret.mdata[0][0]  = c - scalprod;
 371   ret.mdata[0][1]  = a;
 372   ret.mdata[0][2]  = b;
 373
 374   ret.mdata[1][0]  = c*cx;
 375   ret.mdata[1][1]  = a*cx - scalprod;
 376   ret.mdata[1][2]  = b*cx;
 377
 378   ret.mdata[2][0]  = c*cy;
 379   ret.mdata[2][1]  = a*cy;
 380   ret.mdata[2][2]  = b*cy - scalprod;
 381
 382   ret.mIsHomothety = ret.mIsAffine = false;
 383   return ret;
 384 }
 385
 386 const Transformation Transformation::affinityGI3P(
 387   const std::vector<Coordinate>& FromPoints,
 388   const std::vector<Coordinate>& ToPoints,
 389   bool& valid )
 390 {
 391   // construct the (generically) unique affinity that transforms 3 given
 392   // point into 3 other given points; i.e. it depends on the coordinates of
 393   // a total of 6 points.  This actually amounts in solving a 6x6 linear
 394   // system to find the entries of a 2x2 linear transformation matrix T
 395   // and of a translation vector t.
 396   // If Pi denotes one of the starting points and Qi the corresponding
 397   // final position we actually have to solve: Qi = t + T Pi, for i=1,2,3
 398   // (each one is a vector equation, so that it really gives 2 equations).
 399   // In our context T and t are used to build a 3x3 projective transformation
 400   // as follows:
 401   //
 402   //    [  1  0   0  ]
 403   //    [ t1 T11 T12 ]
 404   //    [ t2 T21 T22 ]
 405   //
 406   // In order to take advantage of the two functions "GaussianElimination"
 407   // and "BackwardSubstitution", which are specifically aimed at solving
 408   // homogeneous underdetermined linear systems, we just add a further
 409   // unknown m and solve for t + T Pi - m Qi = 0.  Since our functions
 410   // returns a nonzero solution we shall have a nonzero 'm' in the end and
 411   // can build the 3x3 matrix as follows:
 412   //
 413   //    [  m  0   0  ]
 414   //    [ t1 T11 T12 ]
 415   //    [ t2 T21 T22 ]
 416   //
 417   // we order the unknowns as follows: m, t1, t2, T11, T12, T21, T22
 418
 419   double row0[7], row1[7], row2[7], row3[7], row4[7], row5[7];
 420
 421   double *matrix[6] = {row0, row1, row2, row3, row4, row5};
 422
 423   double solution[7];
 424   int scambio[7];
 425
 426   assert (FromPoints.size() == 3);
 427   assert (ToPoints.size() == 3);
 428
 429   // fill in the matrix elements
 430   for ( int i = 0; i < 6; i++ )
 431   {
 432     for ( int j = 0; j < 7; j++ )
 433     {
 434       matrix[i][j] = 0.0;
 435     }
 436   }
 437
 438   for ( int i = 0; i < 3; i++ )
 439   {
 440     Coordinate p = FromPoints[i];
 441     Coordinate q = ToPoints[i];
 442     matrix[i][0] = -q.x;
 443     matrix[i][1] = 1.0;
 444     matrix[i][3] = p.x;
 445     matrix[i][4] = p.y;
 446     matrix[i+3][0] = -q.y;
 447     matrix[i+3][2] = 1.0;
 448     matrix[i+3][5] = p.x;
 449     matrix[i+3][6] = p.y;
 450   }
 451
 452   Transformation ret;
 453   valid = true;
 454   if ( ! GaussianElimination( matrix, 6, 7, scambio ) )
 455     { valid = false; return ret; }
 456
 457   // fine della fase di eliminazione
 458   BackwardSubstitution( matrix, 6, 7, scambio, solution );
 459
 460   // now we can build the 3x3 transformation matrix; remember that
 461   // unknown 0 is the multiplicator 'm'
 462
 463   ret.mdata[0][0] = solution[0];
 464   ret.mdata[0][1] = ret.mdata[0][2] = 0.0;
 465   ret.mdata[1][0] = solution[1];
 466   ret.mdata[2][0] = solution[2];
 467   ret.mdata[1][1] = solution[3];
 468   ret.mdata[1][2] = solution[4];
 469   ret.mdata[2][1] = solution[5];
 470   ret.mdata[2][2] = solution[6];
 471
 472   ret.mIsHomothety = false;
 473   ret.mIsAffine = true;
 474   return ret;
 475 }
 476
 477 const Transformation Transformation::projectivityGI4P(
 478   const std::vector<Coordinate>& FromPoints,
 479   const std::vector<Coordinate>& ToPoints,
 480   bool& valid )
 481 {
 482   // construct the (generically) unique projectivity that transforms 4 given
 483   // point into 4 other given points; i.e. it depends on the coordinates of
 484   // a total of 8 points.  This actually amounts in solving an underdetermined
 485   // homogeneous linear system.
 486
 487   double
 488     row0[13], row1[13], row2[13], row3[13], row4[13], row5[13], row6[13], row7[13],
 489     row8[13], row9[13], row10[13], row11[13];
 490
 491   double *matrix[12] = {row0, row1, row2, row3, row4, row5, row6, row7,
 492                         row8, row9, row10, row11};
 493
 494   double solution[13];
 495   int scambio[13];
 496
 497   assert (FromPoints.size() == 4);
 498   assert (ToPoints.size() == 4);
 499
 500   // fill in the matrix elements
 501   for ( int i = 0; i < 12; i++ )
 502   {
 503     for ( int j = 0; j < 13; j++ )
 504     {
 505       matrix[i][j] = 0.0;
 506     }
 507   }
 508
 509   for ( int i = 0; i < 4; i++ )
 510   {
 511     Coordinate p = FromPoints[i];
 512     Coordinate q = ToPoints[i];
 513     matrix[i][0] = matrix[4+i][3] = matrix[8+i][6] = 1.0;
 514     matrix[i][1] = matrix[4+i][4] = matrix[8+i][7] = p.x;
 515     matrix[i][2] = matrix[4+i][5] = matrix[8+i][8] = p.y;
 516     matrix[i][9+i] = -1.0;
 517     matrix[4+i][9+i] = -q.x;
 518     matrix[8+i][9+i] = -q.y;
 519   }
 520
 521   Transformation ret;
 522   valid = true;
 523   if ( ! GaussianElimination( matrix, 12, 13, scambio ) )
 524     { valid = false; return ret; }
 525
 526   // fine della fase di eliminazione
 527   BackwardSubstitution( matrix, 12, 13, scambio, solution );
 528
 529   // now we can build the 3x3 transformation matrix; remember that
 530   // unknowns from 9 to 13 are just multiplicators that we don't need here
 531
 532   int k = 0;
 533   for ( int i = 0; i < 3; i++ )
 534   {
 535     for ( int j = 0; j < 3; j++ )
 536     {
 537       ret.mdata[i][j]  = solution[k++];
 538     }
 539   }
 540
 541   ret.mIsHomothety = ret.mIsAffine = false;
 542   return ret;
 543 }
 544
 545 const Transformation Transformation::castShadow(
 546   const Coordinate& lightsrc, const LineData& l )
 547 {
 548   // first deal with the line l, I need to find an appropriate reflection
 549   // that transforms l onto the x-axis
 550
 551   Coordinate d = l.dir();
 552   Coordinate a = l.a;
 553   double k = d.length();
 554   if ( d.x < 0 ) k *= -1;         // for numerical stability
 555   Coordinate w = d + Coordinate( k, 0 );
 556   // w /= w.length();
 557   // w defines a Householder transformation, but we don't need to normalize
 558   // it here.
 559   // warning: this w is the orthogonal of the w of the textbooks!
 560   // this is fine for us since in this way it indicates the line direction
 561   Coordinate ra = Coordinate ( a.x + w.y*a.y/(2*w.x), a.y/2 );
 562   Transformation sym = lineReflection ( LineData( ra, ra + w ) );
 563
 564   // in the new coordinates the line is the x-axis
 565   // I must transform the point
 566
 567   Coordinate modlightsrc = sym.apply ( lightsrc );
 568   Transformation ret = identity();
 569   // parameter t indicates the distance of the light source from
 570   // the plane of the drawing. A negative value means that the light
 571   // source is behind the plane.
 572   double t = -1.0;
 573   //  double t = -modlightsrc.y;    <-- this gives the old transformation!
 574   double e = modlightsrc.y - t;
 575   ret.mdata[0][0] =  e;
 576   ret.mdata[0][2] = -1;
 577   ret.mdata[1][1] =  e;
 578   ret.mdata[1][2] = -modlightsrc.x;
 579   ret.mdata[2][2] =  -t;
 580
 581   ret.mIsHomothety = ret.mIsAffine = false;
 582   return sym*ret*sym;
 583 //  return translation( t )*ret*translation( -t );
 584 }
 585
 586 const Transformation Transformation::projectiveRotation(
 587   double alpha, const Coordinate& d, const Coordinate& t )
 588 {
 589   Transformation ret;
 590   double cosalpha = cos( alpha );
 591   double sinalpha = sin( alpha );
 592   ret.mdata[0][0] =  cosalpha;
 593   ret.mdata[1][1] =  cosalpha*d.x*d.x + d.y*d.y;
 594   ret.mdata[0][1] = -sinalpha*d.x;
 595   ret.mdata[1][0] =  sinalpha*d.x;
 596   ret.mdata[0][2] = -sinalpha*d.y;
 597   ret.mdata[2][0] =  sinalpha*d.y;
 598   ret.mdata[1][2] =  cosalpha*d.x*d.y - d.x*d.y;
 599   ret.mdata[2][1] =  cosalpha*d.x*d.y - d.x*d.y;
 600   ret.mdata[2][2] =  cosalpha*d.y*d.y + d.x*d.x;
 601
 602   ret.mIsHomothety = ret.mIsAffine = false;
 603   return translation( t )*ret*translation( -t );
 604 }
 605
 606 const Coordinate Transformation::apply( const double x0,
 607                                         const double x1,
 608                                         const double x2) const
 609 {
 610   double phom[3] = {x0, x1, x2};
 611   double rhom[3] = {0., 0., 0.};
 612
 613
 614   for (int i = 0; i < 3; i++)
 615   {
 616     for (int j = 0; j < 3; j++)
 617     {
 618       rhom[i] += mdata[i][j]*phom[j];
 619     }
 620   }
 621
 622   if (rhom[0] == 0.)
 623     return Coordinate::invalidCoord();
 624
 625   return Coordinate (rhom[1]/rhom[0], rhom[2]/rhom[0]);
 626 }
 627
 628 const Coordinate Transformation::apply( const Coordinate& p ) const
 629 {
 630   return apply( 1., p.x, p.y );
 631 //  double phom[3] = {1., p.x, p.y};
 632 //  double rhom[3] = {0., 0., 0.};
 633 //
 634 //  for (int i = 0; i < 3; i++)
 635 //  {
 636 //    for (int j = 0; j < 3; j++)
 637 //    {
 638 //      rhom[i] += mdata[i][j]*phom[j];
 639 //    }
 640 //  }
 641 //
 642 //  if (rhom[0] == 0.)
 643 //    return Coordinate::invalidCoord();
 644 //
 645 //  return Coordinate (rhom[1]/rhom[0], rhom[2]/rhom[0]);
 646 }
 647
 648 const Coordinate Transformation::apply0( const Coordinate& p ) const
 649 {
 650   return apply( 0., p.x, p.y );
 651 }
 652
 653 const Transformation Transformation::rotation( double alpha, const Coordinate& center )
 654 {
 655   Transformation ret = identity();
 656
 657   double x = center.x;
 658   double y = center.y;
 659
 660   double cosalpha = cos( alpha );
 661   double sinalpha = sin( alpha );
 662
 663   ret.mdata[1][1] = ret.mdata[2][2] = cosalpha;
 664   ret.mdata[1][2] = -sinalpha;
 665   ret.mdata[2][1] = sinalpha;
 666   ret.mdata[1][0] = x - ret.mdata[1][1]*x - ret.mdata[1][2]*y;
 667   ret.mdata[2][0] = y - ret.mdata[2][1]*x - ret.mdata[2][2]*y;
 668
 669   // this is already set in the identity() constructor, but just for
 670   // clarity..
 671   ret.mIsHomothety = ret.mIsAffine = true;
 672
 673   return ret;
 674 }
 675
 676 bool Transformation::isHomothetic() const
 677 {
 678   return mIsHomothety;
 679 }
 680
 681 bool Transformation::isAffine() const
 682 {
 683   return mIsAffine;
 684 }
 685
 686 /*
 687  *mp:
 688  * this function has the property that it changes sign if computed
 689  * on two points that lie on either sides with respect to the critical
 690  * line (this is the line that goes to the line at infinity).
 691  * For affine transformations the result has always the same sign.
 692  * NOTE: the result is *not* invariant under rescaling of all elements
 693  * of the transformation matrix.
 694  * The typical use is to determine whether a segment is transformed
 695  * into a segment or a couple of half-lines.
 696  */
 697
 698 double Transformation::getProjectiveIndicator( const Coordinate& c ) const
 699 {
 700   return mdata[0][0] + mdata[0][1]*c.x + mdata[0][2]*c.y;
 701 }
 702
 703 // assuming that this is an affine transformation, return its
 704 // determinant.  What is really important here is just the sign
 705 // of the determinant.
 706 double Transformation::getAffineDeterminant() const
 707 {
 708   return mdata[1][1]*mdata[2][2] - mdata[1][2]*mdata[2][1];
 709 }
 710
 711 // this assumes that the 2x2 affine part of the matrix is of the
 712 // form [ cos a, sin a; -sin a, cos a] or a multiple
 713 double Transformation::getRotationAngle() const
 714 {
 715   return atan2( mdata[1][2], mdata[1][1] );
 716 }
 717
 718 const Coordinate Transformation::apply2by2only( const Coordinate& p ) const
 719 {
 720   double x = p.x;
 721   double y = p.y;
 722   double nx = mdata[1][1]*x + mdata[1][2]*y;
 723   double ny = mdata[2][1]*x + mdata[2][2]*y;
 724   return Coordinate( nx, ny );
 725 }
 726
 727 double Transformation::data( int r, int c ) const
 728 {
 729   return mdata[r][c];
 730 }
 731
 732 const Transformation Transformation::inverse( bool& valid ) const
 733 {
 734   Transformation ret;
 735
 736   valid = Invert3by3matrix( mdata, ret.mdata );
 737
 738   // the inverse of a homothety is a homothety, same for affinities..
 739   ret.mIsHomothety = mIsHomothety;
 740   ret.mIsAffine = mIsAffine;
 741
 742   return ret;
 743 }
 744
 745 Transformation::Transformation()
 746 {
 747   // this is the constructor used by the static Transformation
 748   // creation functions, so mIsHomothety is in general false
 749   mIsHomothety = mIsAffine = false;
 750   for ( int i = 0; i < 3; ++i )
 751     for ( int j = 0; j < 3; ++j )
 752       mdata[i][j] = ( i == j ) ? 1 : 0;
 753 }
 754
 755 Transformation::~Transformation()
 756 {
 757 }
 758
 759 double Transformation::apply( double length ) const
 760 {
 761   assert( isHomothetic() );
 762   double det = mdata[1][1]*mdata[2][2] -
 763                mdata[1][2]*mdata[2][1];
 764   return sqrt( fabs( det ) ) * length;
 765 }
 766
 767 Transformation::Transformation( double data[3][3], bool ishomothety )
 768   : mIsHomothety( ishomothety )
 769 {
 770   for ( int i = 0; i < 3; ++i )
 771     for ( int j = 0; j < 3; ++j )
 772       mdata[i][j] = data[i][j];
 773
 774   //mp: a test for affinity is used to initialize mIsAffine...
 775   mIsAffine = false;
 776   if ( fabs(mdata[0][1]) + fabs(mdata[0][2]) < 1e-8 * fabs(mdata[0][0]) )
 777     mIsAffine = true;
 778 }
 779
 780 bool operator==( const Transformation& lhs, const Transformation& rhs )
 781 {
 782   for ( int i = 0; i < 3; ++i )
 783     for ( int j = 0; j < 3; ++j )
 784       if ( lhs.data( i, j ) != rhs.data( i, j ) )
 785         return false;
 786   return true;
 787 }
 788
 789 const Transformation Transformation::similitude(
 790   const Coordinate& center, double theta, double factor )
 791 {
 792   //kdDebug() << k_funcinfo << "theta: " << theta << " factor: " << factor << endl;
 793   Transformation ret;
 794   ret.mIsHomothety = true;
 795   double costheta = cos( theta );
 796   double sintheta = sin( theta );
 797   ret.mdata[0][0] = 1;
 798   ret.mdata[0][1] = 0;
 799   ret.mdata[0][2] = 0;
 800   ret.mdata[1][0] = ( 1 - factor*costheta )*center.x + factor*sintheta*center.y;
 801   ret.mdata[1][1] = factor*costheta;
 802   ret.mdata[1][2] = -factor*sintheta;
 803   ret.mdata[2][0] = -factor*sintheta*center.x + ( 1 - factor*costheta )*center.y;
 804   ret.mdata[2][1] = factor*sintheta;
 805   ret.mdata[2][2] = factor*costheta;
 806   // fails for factor == infinity
 807   //assert( ( ret.apply( center ) - center ).length() < 1e-5 );
 808   ret.mIsHomothety = ret.mIsAffine = true;
 809   return ret;
 810 }