kig/objects/cubic_imp.cc

   1 // Copyright (C)  2003  Dominique Devriese <devriese@kde.org>
   2
   3 // This program is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU General Public License
   5 // as published by the Free Software Foundation; either version 2
   6 // of the License, or (at your option) any later version.
   7
   8 // This program is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 // GNU General Public License for more details.
  12
  13 // You should have received a copy of the GNU General Public License
  14 // along with this program; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  16 // 02111-1307, USA.
  17
  18 #include "cubic_imp.h"
  19
  20 #include "bogus_imp.h"
  21
  22 #include "../misc/kigpainter.h"
  23 #include "../misc/screeninfo.h"
  24 #include "../misc/kignumerics.h"
  25 #include "../misc/common.h"
  26 #include "../kig/kig_view.h"
  27
  28 #include <math.h>
  29 #include <klocale.h>
  30
  31 CubicImp::CubicImp( const CubicCartesianData& data )
  32   : CurveImp(), mdata( data )
  33 {
  34 }
  35
  36 CubicImp::~CubicImp()
  37 {
  38 }
  39
  40 ObjectImp* CubicImp::transform( const Transformation& t ) const
  41 {
  42   bool valid = true;
  43   CubicCartesianData d = calcCubicTransformation( data(), t, valid );
  44   if ( valid ) return new CubicImp( d );
  45   else return new InvalidImp;
  46 }
  47
  48 void CubicImp::draw( KigPainter& p ) const
  49 {
  50   p.drawCurve( this );
  51 }
  52
  53 bool CubicImp::contains( const Coordinate& o, int width, const KigWidget& w ) const
  54 {
  55   return internalContainsPoint( o, w.screenInfo().normalMiss( width ) );
  56 }
  57
  58 bool CubicImp::inRect( const Rect&, int, const KigWidget& ) const
  59 {
  60   // TODO ?
  61   return false;
  62 }
  63
  64 CubicImp* CubicImp::copy() const
  65 {
  66   return new CubicImp( mdata );
  67 }
  68
  69 double CubicImp::getParam( const Coordinate& p, const KigDocument& ) const
  70 {
  71   double x = p.x;
  72   double y = p.y;
  73   double t;
  74
  75   double a000 = mdata.coeffs[0];
  76   double a001 = mdata.coeffs[1];
  77   double a002 = mdata.coeffs[2];
  78   double a011 = mdata.coeffs[3];
  79   double a012 = mdata.coeffs[4];
  80   double a022 = mdata.coeffs[5];
  81   double a111 = mdata.coeffs[6];
  82   double a112 = mdata.coeffs[7];
  83   double a122 = mdata.coeffs[8];
  84   double a222 = mdata.coeffs[9];
  85
  86   /*
  87    * first project p onto the cubic. This is done by computing the
  88    * line through p in the direction of the gradient
  89    */
  90
  91   double f = a000 + a001*x + a002*y + a011*x*x + a012*x*y + a022*y*y +
  92              a111*x*x*x + a112*x*x*y + a122*x*y*y + a222*y*y*y;
  93   if ( f != 0 )
  94   {
  95     double fx = a001 + 2*a011*x + a012*y + 3*a111*x*x + 2*a112*x*y + a122*y*y;
  96     double fy = a002 + 2*a022*y + a012*x + 3*a222*y*y + 2*a122*x*y + a112*x*x;
  97     Coordinate v = Coordinate (fx, fy);
  98     if ( f < 0 ) v = -v;   // the line points away from the intersection
  99     double a, b, c, d;
 100     calcCubicLineRestriction ( mdata, p, v, a, b, c, d );
 101     if ( a < 0 )
 102     {
 103       a *= -1;
 104       b *= -1;
 105       c *= -1;
 106       d *= -1;
 107     }
 108
 109     // computing the coefficients of the Sturm sequence
 110     double p1a = 2*b*b - 6*a*c;
 111     double p1b = b*c - 9*a*d;
 112     double p0a = c*p1a*p1a + p1b*(3*a*p1b - 2*b*p1a);
 113     // compute the number of roots for negative lambda
 114     int variations = calcCubicVariations ( 0, a, b, c, d, p1a, p1b, p0a );
 115     bool valid;
 116     int numroots;
 117     double lambda = calcCubicRoot ( -1e10, 1e10, a, b, c, d, variations, valid,
 118                                     numroots );
 119     if ( valid )
 120     {
 121       Coordinate pnew = p + lambda*v;
 122       x = pnew.x;
 123       y = pnew.y;
 124     }
 125   }
 126
 127   if (x > 0) t = x/(1+x);
 128   else t = x/(1-x);
 129   t = 0.5*(t + 1);
 130   t /= 3;
 131
 132   Coordinate p1 = getPoint ( t );
 133   Coordinate p2 = getPoint ( t + 1.0/3.0 );
 134   Coordinate p3 = getPoint ( t + 2.0/3.0 );
 135
 136   double mint = t;
 137   double mindist = p1.valid() ? fabs ( y - p1.y ) : double_inf;
 138   if ( p2.valid() && fabs ( y - p2.y ) < mindist )
 139   {
 140     mint = t + 1.0/3.0;
 141     mindist = fabs ( y - p2.y );
 142   }
 143   if ( p3.valid() && fabs ( y - p3.y ) < mindist )
 144   {
 145     mint = t + 2.0/3.0;
 146   }
 147
 148   return mint;
 149 }
 150
 151 const Coordinate CubicImp::getPoint( double p, const KigDocument& ) const
 152 {
 153   return getPoint( p );
 154 }
 155
 156 const Coordinate CubicImp::getPoint( double p ) const
 157 {
 158   /*
 159    * this isn't really elegant...
 160    * the magnitude of p tells which one of the maximum 3 intersections
 161    * of the vertical line with the cubic to take.
 162    */
 163
 164   p *= 3;
 165   int root = (int) p;
 166   assert ( root >= 0 );
 167   assert ( root <= 3 );
 168   if ( root == 3 ) root = 2;
 169
 170   p -= root;
 171
 172   assert ( 0 <= p && p <= 1 );
 173   if ( p <= 0. ) p = 1e-6;
 174   if ( p >= 1. ) p = 1 - 1e-6;
 175   root++;
 176   p = 2*p - 1;
 177   double x;
 178   if (p > 0) x = p/(1 - p);
 179   else x = p/(1 + p);
 180
 181   // calc the third degree polynomial:
 182   // compute the third degree polinomial:
 183 //  double a000 = mdata.coeffs[0];
 184 //  double a001 = mdata.coeffs[1];
 185 //  double a002 = mdata.coeffs[2];
 186 //  double a011 = mdata.coeffs[3];
 187 //  double a012 = mdata.coeffs[4];
 188 //  double a022 = mdata.coeffs[5];
 189 //  double a111 = mdata.coeffs[6];
 190 //  double a112 = mdata.coeffs[7];
 191 //  double a122 = mdata.coeffs[8];
 192 //  double a222 = mdata.coeffs[9];
 193 //
 194 //  // first the y^3 coefficient, it coming only from a222:
 195 //  double a = a222;
 196 //  // next the y^2 coefficient (from a122 and a022):
 197 //  double b = a122*x + a022;
 198 //  // next the y coefficient (from a112, a012 and a002):
 199 //  double c = a112*x*x + a012*x + a002;
 200 //  // finally the constant coefficient (from a111, a011, a001 and a000):
 201 //  double d = a111*x*x*x + a011*x*x + a001*x + a000;
 202
 203 // commented out, since the bound is already computed when passing a huge
 204 // interval; the normalization is not needed also
 205
 206   // renormalize: positive a
 207 //  if ( a < 0 )
 208 //  {
 209 //    a *= -1;
 210 //    b *= -1;
 211 //    c *= -1;
 212 //    d *= -1;
 213 //  }
 214 //
 215 //  const double small = 1e-7;
 216 //  int degree = 3;
 217 //  if ( fabs(a) < small*fabs(b) ||
 218 //       fabs(a) < small*fabs(c) ||
 219 //       fabs(a) < small*fabs(d) )
 220 //  {
 221 //    degree = 2;
 222 //    if ( fabs(b) < small*fabs(c) ||
 223 //         fabs(b) < small*fabs(d) )
 224 //    {
 225 //      degree = 1;
 226 //    }
 227 //  }
 228
 229 // and a bound for all the real roots:
 230
 231 //  double bound;
 232 //  switch (degree)
 233 //  {
 234 //    case 3:
 235 //    bound = fabs(d/a);
 236 //    if ( fabs(c/a) + 1 > bound ) bound = fabs(c/a) + 1;
 237 //    if ( fabs(b/a) + 1 > bound ) bound = fabs(b/a) + 1;
 238 //    break;
 239 //
 240 //    case 2:
 241 //    bound = fabs(d/b);
 242 //    if ( fabs(c/b) + 1 > bound ) bound = fabs(c/b) + 1;
 243 //    break;
 244 //
 245 //    case 1:
 246 //    default:
 247 //    bound = fabs(d/c) + 1;
 248 //    break;
 249 //  }
 250
 251   int numroots;
 252   bool valid = true;
 253   double y = calcCubicYvalue ( x, -double_inf, double_inf, root, mdata, valid,
 254                                numroots );
 255   if ( ! valid ) return Coordinate::invalidCoord();
 256   else return Coordinate(x,y);
 257 //  if ( valid ) return Coordinate(x,y);
 258 //  root--; if ( root <= 0) root += 3;
 259 //  y = calcCubicYvalue ( x, -bound, bound, root, mdata, valid,
 260 //                        numroots );
 261 //  if ( valid ) return Coordinate(x,y);
 262 //  root--; if ( root <= 0) root += 3;
 263 //  y = calcCubicYvalue ( x, -bound, bound, root, mdata, valid,
 264 //                        numroots );
 265 //  assert ( valid );
 266 //  return Coordinate(x,y);
 267 }
 268
 269 const uint CubicImp::numberOfProperties() const
 270 {
 271   return Parent::numberOfProperties() + 1;
 272 }
 273
 274 const QCStringList CubicImp::propertiesInternalNames() const
 275 {
 276         QCStringList l = Parent::propertiesInternalNames();
 277   l << "cartesian-equation";
 278   assert( l.size() == CubicImp::numberOfProperties() );
 279   return l;
 280
 281 }
 282
 283 /*
 284  * cartesian equation property contributed by Carlo Sardi
 285  */
 286
 287 const QCStringList CubicImp::properties() const
 288 {
 289         QCStringList l = Parent::properties();
 290   l << I18N_NOOP( "Cartesian Equation" );
 291   assert( l.size() == CubicImp::numberOfProperties() );
 292   return l;
 293
 294 }
 295
 296 const ObjectImpType* CubicImp::impRequirementForProperty( uint which ) const
 297 {
 298         if ( which < Parent::numberOfProperties() )
 299     return Parent::impRequirementForProperty( which );
 300   else return CubicImp::stype();
 301 }
 302
 303 const char* CubicImp::iconForProperty( uint which ) const
 304 {
 305         int pnum = 0;
 306   if ( which < Parent::numberOfProperties() )
 307     return Parent::iconForProperty( which );
 308   if ( which == Parent::numberOfProperties() + pnum++ )
 309     return "text"; // cartesian equation string
 310   else
 311                 assert( false );
 312   return "";
 313 }
 314
 315 ObjectImp* CubicImp::property( uint which, const KigDocument& w ) const
 316 {
 317         int pnum = 0;
 318
 319   if ( which < Parent::numberOfProperties() )
 320     return Parent::property( which, w );
 321         if ( which == Parent::numberOfProperties() + pnum++ )
 322     return new StringImp( cartesianEquationString( w ) );
 323         else
 324                 assert( false );
 325         return new InvalidImp;
 326 }
 327
 328 const CubicCartesianData CubicImp::data() const
 329 {
 330   return mdata;
 331 }
 332
 333 void CubicImp::visit( ObjectImpVisitor* vtor ) const
 334 {
 335   vtor->visit( this );
 336 }
 337
 338 bool CubicImp::equals( const ObjectImp& rhs ) const
 339 {
 340   return rhs.inherits( CubicImp::stype() ) &&
 341     static_cast<const CubicImp&>( rhs ).data() == data();
 342 }
 343
 344 const ObjectImpType* CubicImp::type() const
 345 {
 346   return CubicImp::stype();
 347 }
 348
 349 const ObjectImpType* CubicImp::stype()
 350 {
 351   static const ObjectImpType t(
 352     Parent::stype(), "cubic",
 353     I18N_NOOP( "cubic curve" ),
 354     I18N_NOOP( "Select this cubic curve" ),
 355     I18N_NOOP( "Select cubic curve %1" ),
 356     I18N_NOOP( "Remove a Cubic Curve" ),
 357     I18N_NOOP( "Add a Cubic Curve" ),
 358     I18N_NOOP( "Move a Cubic Curve" ),
 359     I18N_NOOP( "Attach to this cubic curve" ),
 360     I18N_NOOP( "Show a Cubic Curve" ),
 361     I18N_NOOP( "Hide a Cubic Curve" )
 362     );
 363   return &t;
 364 }
 365
 366 bool CubicImp::containsPoint( const Coordinate& p, const KigDocument& ) const
 367 {
 368   return internalContainsPoint( p, test_threshold );
 369 }
 370
 371 bool CubicImp::internalContainsPoint( const Coordinate& p, double threshold ) const
 372 {
 373   double a000 = mdata.coeffs[0];
 374   double a001 = mdata.coeffs[1];
 375   double a002 = mdata.coeffs[2];
 376   double a011 = mdata.coeffs[3];
 377   double a012 = mdata.coeffs[4];
 378   double a022 = mdata.coeffs[5];
 379   double a111 = mdata.coeffs[6];
 380   double a112 = mdata.coeffs[7];
 381   double a122 = mdata.coeffs[8];
 382   double a222 = mdata.coeffs[9];
 383
 384   double x = p.x;
 385   double y = p.y;
 386
 387   double f = a000 + a001*x + a002*y + a011*x*x + a012*x*y + a022*y*y +
 388              a111*x*x*x + a112*x*x*y + a122*x*y*y + a222*y*y*y;
 389   double fx = a001 + 2*a011*x + a012*y + 3*a111*x*x + 2*a112*x*y + a122*y*y;
 390   double fy = a002 + a012*x + 2*a022*y + a112*x*x + 2*a122*x*y + 3*a222*y*y;
 391
 392   double dist = fabs(f)/(fabs(fx) + fabs(fy));
 393
 394   return dist <= threshold;
 395 }
 396
 397 bool CubicImp::isPropertyDefinedOnOrThroughThisImp( uint which ) const
 398 {
 399   return Parent::isPropertyDefinedOnOrThroughThisImp( which );
 400 }
 401
 402 Rect CubicImp::surroundingRect() const
 403 {
 404   // it's probably possible to calculate this if it exists, but for
 405   // now we don't.
 406   return Rect::invalidRect();
 407 }
 408
 409 QString CubicImp::cartesianEquationString( const KigDocument& ) const
 410 {
 411   QString ret = i18n( "%7 x³+ %10 y³ + %8 x²y + %9 xy² +%6 y² + %4 x² + %5 xy + %2 x + %3 y + %1 = 0" );
 412   ret = ret.arg( mdata.coeffs[0], 0, 'g', 3 );
 413   ret = ret.arg( mdata.coeffs[1], 0, 'g', 3 );
 414   ret = ret.arg( mdata.coeffs[2], 0, 'g', 3 );
 415   ret = ret.arg( mdata.coeffs[3], 0, 'g', 3 );
 416   ret = ret.arg( mdata.coeffs[4], 0, 'g', 3 );
 417   ret = ret.arg( mdata.coeffs[5], 0, 'g', 3 );
 418         ret = ret.arg( mdata.coeffs[6], 0, 'g', 3 );
 419   ret = ret.arg( mdata.coeffs[7], 0, 'g', 3 );
 420   ret = ret.arg( mdata.coeffs[8], 0, 'g', 3 );
 421   ret = ret.arg( mdata.coeffs[9], 0, 'g', 3 );
 422
 423   return ret;
 424 /*      double a000 = mdata.coeffs[0];
 425   double a001 = mdata.coeffs[1];
 426   double a002 = mdata.coeffs[2];
 427   double a011 = mdata.coeffs[3];
 428   double a012 = mdata.coeffs[4];
 429   double a022 = mdata.coeffs[5];
 430   double a111 = mdata.coeffs[6];
 431   double a112 = mdata.coeffs[7];
 432   double a122 = mdata.coeffs[8];
 433   double a222 = mdata.coeffs[9];*/
 434         //
 435 //  // first the y^3 coefficient, it coming only from a222:
 436 //  double a = a222;
 437 //  // next the y^2 coefficient (from a122 and a022):
 438 //  double b = a122*x + a022;
 439 //  // next the y coefficient (from a112, a012 and a002):
 440 //  double c = a112*x*x + a012*x + a002;
 441 //  // finally the constant coefficient (from a111, a011, a001 and a000):
 442 //  double d = a111*x*x*x + a011*x*x + a001*x + a000;
 443
 444 }