Shamusworld >> Repos - architektonas/blob - src/base/spline.cpp

   1 // spline.cpp
   2 //
   3 // Part of the Architektonas Project
   4 // Originally part of QCad Community Edition by Andrew Mustun
   5 // Extensively rewritten and refactored by James L. Hammons
   6 // Portions copyright (C) 2001-2003 RibbonSoft
   7 // Copyright (C) 2010 Underground Software
   8 // See the README and GPLv2 files for licensing and warranty information
   9 //
  10 // JLH = James L. Hammons <jlhamm@acm.org>
  11 //
  12 // Who  When        What
  13 // ---  ----------  -----------------------------------------------------------
  14 // JLH  06/02/2010  Added this text. :-)
  15 //
  16
  17 #include "spline.h"
  18
  19 #include "debug.h"
  20 #include "graphicview.h"
  21 #include "drawing.h"
  22 #include "paintinterface.h"
  23
  24 /**
  25  * Constructor.
  26  */
  27 Spline::Spline(EntityContainer * parent, const SplineData & d):
  28         EntityContainer(parent), data(d)
  29 {
  30         calculateBorders();
  31 }
  32
  33 /**
  34  * Destructor.
  35  */
  36 Spline::~Spline()
  37 {
  38 }
  39
  40 Entity * Spline::clone()
  41 {
  42         Spline * l = new Spline(*this);
  43 #warning "!!! Need to deal with setAutoDelete() Qt3->Qt4 !!!"
  44 //      l->entities.setAutoDelete(entities.autoDelete());
  45         l->initId();
  46         l->detach();
  47         return l;
  48 }
  49
  50 /**     @return RS2::EntitySpline */
  51 /*virtual*/ RS2::EntityType Spline::rtti() const
  52 {
  53         return RS2::EntitySpline;
  54 }
  55
  56 /** @return false */
  57 /*virtual*/ bool Spline::isEdge() const
  58 {
  59         return false;
  60 }
  61
  62 /** @return Copy of data that defines the spline. */
  63 SplineData Spline::getData() const
  64 {
  65         return data;
  66 }
  67
  68 /** Sets the splines degree (1-3). */
  69 void Spline::setDegree(int deg)
  70 {
  71         if (deg >= 1 && deg <= 3)
  72                 data.degree = deg;
  73 }
  74
  75 /** @return Degree of this spline curve (1-3).*/
  76 int Spline::getDegree()
  77 {
  78         return data.degree;
  79 }
  80
  81 /** @return 0. */
  82 int Spline::getNumberOfKnots()
  83 {
  84         return 0;
  85 }
  86
  87 /** @return Number of control points. */
  88 int Spline::getNumberOfControlPoints()
  89 {
  90         return data.controlPoints.count();
  91 }
  92
  93 /**
  94  * @retval true if the spline is closed.
  95  * @retval false otherwise.
  96  */
  97 bool Spline::isClosed()
  98 {
  99         return data.closed;
 100 }
 101
 102 /**
 103  * Sets the closed falg of this spline.
 104  */
 105 void Spline::setClosed(bool c)
 106 {
 107         data.closed = c;
 108         update();
 109 }
 110
 111 void Spline::calculateBorders()
 112 {
 113     /*minV = Vector::minimum(data.startpoint, data.endpoint);
 114     maxV = Vector::maximum(data.startpoint, data.endpoint);
 115
 116     Q3ValueList<Vector>::iterator it;
 117     for (it = data.controlPoints.begin();
 118     it!=data.controlPoints.end(); ++it) {
 119
 120     minV = Vector::minimum(*it, minV);
 121     maxV = Vector::maximum(*it, maxV);
 122 }
 123     */
 124 }
 125
 126 VectorSolutions Spline::getRefPoints()
 127 {
 128         VectorSolutions ret(data.controlPoints.count());
 129
 130         int i = 0;
 131         QList<Vector>::iterator it;
 132
 133         for(it=data.controlPoints.begin(); it!=data.controlPoints.end(); ++it, ++i)
 134         {
 135                 ret.set(i, (*it));
 136         }
 137
 138         return ret;
 139 }
 140
 141 Vector Spline::getNearestRef(const Vector & coord, double * dist)
 142 {
 143         //return getRefPoints().getClosest(coord, dist);
 144         return Entity::getNearestRef(coord, dist);
 145 }
 146
 147 Vector Spline::getNearestSelectedRef(const Vector & coord, double * dist)
 148 {
 149         //return getRefPoints().getClosest(coord, dist);
 150         return Entity::getNearestSelectedRef(coord, dist);
 151 }
 152
 153 /**
 154  * Updates the internal polygon of this spline. Called when the
 155  * spline or it's data, position, .. changes.
 156  */
 157 void Spline::update()
 158 {
 159         DEBUG->print("Spline::update");
 160
 161         clear();
 162
 163         if (isUndone())
 164                 return;
 165
 166         if (data.degree < 1 || data.degree > 3)
 167         {
 168                 DEBUG->print("Spline::update: invalid degree: %d", data.degree);
 169                 return;
 170         }
 171
 172         if (data.controlPoints.count() < (uint)data.degree + 1)
 173         {
 174                 DEBUG->print("Spline::update: not enough control points");
 175                 return;
 176         }
 177
 178         resetBorders();
 179
 180 //      Q3ValueList<Vector> tControlPoints = data.controlPoints;
 181         QList<Vector> tControlPoints = data.controlPoints;
 182
 183         if (data.closed)
 184         {
 185                 for(int i=0; i<data.degree; ++i)
 186                         tControlPoints.append(data.controlPoints[i]);
 187         }
 188
 189         int i;
 190         int npts = tControlPoints.count();
 191         // order:
 192         int k = data.degree + 1;
 193         // resolution:
 194         int p1 = getGraphicVariableInt("$SPLINESEGS", 8) * npts;
 195
 196         double * b = new double[npts * 3 + 1];
 197         double * h = new double[npts + 1];
 198         double * p = new double[p1 * 3 + 1];
 199
 200 //      Q3ValueList<Vector>::iterator it;
 201         QList<Vector>::iterator it;
 202         i = 1;
 203
 204         for(it=tControlPoints.begin(); it!=tControlPoints.end(); ++it)
 205         {
 206                 b[i + 0] = (*it).x;
 207                 b[i + 1] = (*it).y;
 208                 b[i + 2] = 0.0;
 209
 210                 DEBUG->print("Spline::update: b[%d]: %f/%f", i, b[i], b[i + 1]);
 211                 i += 3;
 212         }
 213
 214         // set all homogeneous weighting factors to 1.0
 215         for(i=1; i<=npts; i++)
 216                 h[i] = 1.0;
 217
 218         for(i=1; i<=3*p1; i++)
 219                 p[i] = 0.0;
 220
 221         if (data.closed)
 222                 rbsplinu(npts, k, p1, b, h, p);
 223         else
 224                 rbspline(npts, k, p1, b, h, p);
 225
 226         Vector prev(false);
 227
 228         for(i=1; i<=3*p1; i=i+3)
 229         {
 230                 if (prev.valid)
 231                 {
 232                         Line * line = new Line(this, LineData(prev, Vector(p[i], p[i + 1])));
 233                         line->setLayer(NULL);
 234                         line->setPen(Pen(RS2::FlagInvalid));
 235                         addEntity(line);
 236                 }
 237
 238                 prev = Vector(p[i], p[i+1]);
 239                 minV = Vector::minimum(prev, minV);
 240                 maxV = Vector::maximum(prev, maxV);
 241         }
 242
 243         delete[] b;
 244         delete[] h;
 245         delete[] p;
 246 }
 247
 248 Vector Spline::getNearestEndpoint(const Vector & coord, double * dist)
 249 {
 250     double minDist = RS_MAXDOUBLE;
 251     double d;
 252     Vector ret(false);
 253
 254     for (uint i=0; i<data.controlPoints.count(); i++) {
 255         d = data.controlPoints[i].distanceTo(coord);
 256
 257         if (d<minDist) {
 258             minDist = d;
 259             ret = data.controlPoints[i];
 260         }
 261     }
 262     if (dist!=NULL) {
 263         *dist = minDist;
 264     }
 265         return ret;
 266 }
 267
 268 /*
 269 // The default implementation of EntityContainer is inaccurate but
 270 //   has to do for now..
 271 Vector Spline::getNearestPointOnEntity(const Vector& coord,
 272         bool onEntity, double* dist, Entity** entity) {
 273 }
 274 */
 275
 276 Vector Spline::getNearestCenter(const Vector & /*coord*/, double * dist)
 277 {
 278         if (dist != NULL)
 279                 *dist = RS_MAXDOUBLE;
 280
 281         return Vector(false);
 282 }
 283
 284 Vector Spline::getNearestMiddle(const Vector & /*coord*/, double * dist)
 285 {
 286         if (dist!=NULL)
 287                 *dist = RS_MAXDOUBLE;
 288
 289         return Vector(false);
 290 }
 291
 292 Vector Spline::getNearestDist(double /*distance*/, const Vector & /*coord*/, double * dist)
 293 {
 294         if (dist != NULL)
 295                 *dist = RS_MAXDOUBLE;
 296
 297         return Vector(false);
 298 }
 299
 300 void Spline::move(Vector offset)
 301 {
 302 //      Q3ValueList<Vector>::iterator it;
 303         QList<Vector>::iterator it;
 304
 305         for(it=data.controlPoints.begin(); it!=data.controlPoints.end(); ++it)
 306                 (*it).move(offset);
 307
 308         update();
 309 }
 310
 311 void Spline::rotate(Vector center, double angle)
 312 {
 313 //      Q3ValueList<Vector>::iterator it;
 314         QList<Vector>::iterator it;
 315
 316         for(it=data.controlPoints.begin(); it!=data.controlPoints.end(); ++it)
 317                 (*it).rotate(center, angle);
 318
 319         update();
 320 }
 321
 322 void Spline::scale(Vector center, Vector factor)
 323 {
 324 //      Q3ValueList<Vector>::iterator it;
 325         QList<Vector>::iterator it;
 326
 327         for(it=data.controlPoints.begin(); it!=data.controlPoints.end(); ++it)
 328                 (*it).scale(center, factor);
 329
 330         update();
 331 }
 332
 333 void Spline::mirror(Vector axisPoint1, Vector axisPoint2)
 334 {
 335 //      Q3ValueList<Vector>::iterator it;
 336         QList<Vector>::iterator it;
 337
 338         for(it=data.controlPoints.begin(); it!=data.controlPoints.end(); ++it)
 339                 (*it).mirror(axisPoint1, axisPoint2);
 340
 341         update();
 342 }
 343
 344 void Spline::moveRef(const Vector & ref, const Vector & offset)
 345 {
 346 //      Q3ValueList<Vector>::iterator it;
 347         QList<Vector>::iterator it;
 348
 349         for(it=data.controlPoints.begin(); it!=data.controlPoints.end(); ++it)
 350         {
 351                 if (ref.distanceTo(*it) < 1.0e-4)
 352                         (*it).move(offset);
 353         }
 354
 355         update();
 356 }
 357
 358 void Spline::draw(PaintInterface * painter, GraphicView * view, double /*patternOffset*/)
 359 {
 360         if (painter == NULL || view == NULL)
 361                 return;
 362
 363         Entity * e = firstEntity(RS2::ResolveNone);
 364         double offset = 0.0;
 365
 366         if (e != NULL)
 367         {
 368                 view->drawEntity(e);
 369                 offset += e->getLength();
 370                 //DEBUG->print("offset: %f\nlength was: %f", offset, e->getLength());
 371         }
 372
 373         for (Entity * e=nextEntity(RS2::ResolveNone); e!=NULL; e = nextEntity(RS2::ResolveNone))
 374         {
 375                 view->drawEntityPlain(e, -offset);
 376                 offset += e->getLength();
 377                 //DEBUG->print("offset: %f\nlength was: %f", offset, e->getLength());
 378         }
 379 }
 380
 381 /**
 382  * Todo: draw the spline, user patterns.
 383  */
 384 /*
 385 void Spline::draw(RS_Painter* painter, GraphicView* view) {
 386    if (painter==NULL || view==NULL) {
 387        return;
 388    }
 389
 390    / *
 391       if (data.controlPoints.count()>0) {
 392           Vector prev(false);
 393           Q3ValueList<Vector>::iterator it;
 394           for (it = data.controlPoints.begin(); it!=data.controlPoints.end(); ++it) {
 395               if (prev.valid) {
 396                   painter->drawLine(view->toGui(prev),
 397                                     view->toGui(*it));
 398               }
 399               prev = (*it);
 400           }
 401       }
 402    * /
 403
 404    int i;
 405    int npts = data.controlPoints.count();
 406    // order:
 407    int k = 4;
 408    // resolution:
 409    int p1 = 100;
 410
 411    double* b = new double[npts*3+1];
 412    double* h = new double[npts+1];
 413    double* p = new double[p1*3+1];
 414
 415    Q3ValueList<Vector>::iterator it;
 416    i = 1;
 417    for (it = data.controlPoints.begin(); it!=data.controlPoints.end(); ++it) {
 418        b[i] = (*it).x;
 419        b[i+1] = (*it).y;
 420        b[i+2] = 0.0;
 421
 422         DEBUG->print("Spline::draw: b[%d]: %f/%f", i, b[i], b[i+1]);
 423         i+=3;
 424    }
 425
 426    // set all homogeneous weighting factors to 1.0
 427    for (i=1; i <= npts; i++) {
 428        h[i] = 1.0;
 429    }
 430
 431    //
 432    for (i = 1; i <= 3*p1; i++) {
 433        p[i] = 0.0;
 434    }
 435
 436    rbspline(npts,k,p1,b,h,p);
 437
 438    Vector prev(false);
 439    for (i = 1; i <= 3*p1; i=i+3) {
 440        if (prev.valid) {
 441            painter->drawLine(view->toGui(prev),
 442                              view->toGui(Vector(p[i], p[i+1])));
 443        }
 444        prev = Vector(p[i], p[i+1]);
 445    }
 446 }
 447 */
 448
 449 /**
 450  * @return The reference points of the spline.
 451  */
 452 //Q3ValueList<Vector> Spline::getControlPoints()
 453 QList<Vector> Spline::getControlPoints()
 454 {
 455         return data.controlPoints;
 456 }
 457
 458 /**
 459  * Appends the given point to the control points.
 460  */
 461 void Spline::addControlPoint(const Vector & v)
 462 {
 463         data.controlPoints.append(v);
 464 }
 465
 466 /**
 467  * Removes the control point that was last added.
 468  */
 469 void Spline::removeLastControlPoint()
 470 {
 471         data.controlPoints.pop_back();
 472 }
 473
 474 /**
 475  * Generates B-Spline open knot vector with multiplicity
 476  * equal to the order at the ends.
 477  */
 478 void Spline::knot(int num, int order, int knotVector[])
 479 {
 480         knotVector[1] = 0;
 481
 482         for(int i=2; i<=num+order; i++)
 483         {
 484                 if ((i > order) && (i < num + 2))
 485                 {
 486                         knotVector[i] = knotVector[i - 1] + 1;
 487                 }
 488                 else
 489                 {
 490                         knotVector[i] = knotVector[i - 1];
 491                 }
 492         }
 493 }
 494
 495 /**
 496  * Generates rational B-spline basis functions for an open knot vector.
 497  */
 498 void Spline::rbasis(int c, double t, int npts, int x[], double h[], double r[])
 499 {
 500         int nplusc;
 501         int i, k;
 502         double d, e;
 503         double sum;
 504         //double temp[36];
 505
 506         nplusc = npts + c;
 507
 508         double * temp = new double[nplusc + 1];
 509
 510         // calculate the first order nonrational basis functions n[i]
 511         for(i=1; i<=nplusc-1; i++)
 512         {
 513                 if ((t >= x[i]) && (t < x[i + 1]))
 514                         temp[i] = 1;
 515                 else
 516                         temp[i] = 0;
 517         }
 518
 519         /* calculate the higher order nonrational basis functions */
 520
 521         for(k=2; k<=c; k++)
 522         {
 523                 for(i=1; i<=nplusc-k; i++)
 524                 {
 525                         // if the lower order basis function is zero skip the calculation
 526                         if (temp[i] != 0)
 527                                 d = ((t - x[i]) * temp[i]) / (x[i + k - 1] - x[i]);
 528                         else
 529                                 d = 0;
 530
 531                         // if the lower order basis function is zero skip the calculation
 532                         if (temp[i + 1] != 0)
 533                                 e = ((x[i + k] - t) * temp[i + 1]) / (x[i + k] - x[i + 1]);
 534                         else
 535                                 e = 0;
 536
 537                         temp[i] = d + e;
 538                 }
 539         }
 540
 541         // pick up last point
 542         if (t == (double)x[nplusc])
 543                 temp[npts] = 1;
 544
 545         // calculate sum for denominator of rational basis functions
 546         sum = 0.;
 547
 548         for(i=1; i<=npts; i++)
 549                 sum = sum + temp[i] * h[i];
 550
 551         // form rational basis functions and put in r vector
 552         for(i=1; i<=npts; i++)
 553         {
 554                 if (sum != 0)
 555                         r[i] = (temp[i] * h[i]) / (sum);
 556                 else
 557                         r[i] = 0;
 558         }
 559
 560         delete[] temp;
 561 }
 562
 563 /**
 564  * Generates a rational B-spline curve using a uniform open knot vector.
 565  */
 566 void Spline::rbspline(int npts, int k, int p1, double b[], double h[], double p[])
 567 {
 568         int i, j, icount, jcount;
 569         int i1;
 570         //int x[30]; /* allows for 20 data points with basis function of order 5 */
 571         int nplusc;
 572
 573         double step;
 574         double t;
 575         //double nbasis[20];
 576         double temp;
 577
 578         nplusc = npts + k;
 579
 580         int * x = new int[nplusc + 1];
 581         double * nbasis = new double[npts + 1];
 582
 583         // zero and redimension the knot vector and the basis array
 584
 585         for(i = 0; i <= npts; i++)
 586                 nbasis[i] = 0.0;
 587
 588         for(i = 0; i <= nplusc; i++)
 589                 x[i] = 0;
 590
 591         // generate the uniform open knot vector
 592         knot(npts, k, x);
 593
 594         icount = 0;
 595
 596         // calculate the points on the rational B-spline curve
 597         t = 0;
 598         step = ((double)x[nplusc]) / ((double)(p1 - 1));
 599
 600         for(i1=1; i1<= p1; i1++)
 601         {
 602                 if ((double)x[nplusc] - t < 5e-6)
 603                         t = (double)x[nplusc];
 604
 605                 // generate the basis function for this value of t
 606                 rbasis(k, t, npts, x, h, nbasis);
 607
 608                 // generate a point on the curve
 609                 for(j=1; j<=3; j++)
 610                 {
 611                         jcount = j;
 612                         p[icount + j] = 0.0;
 613
 614                         // Do local matrix multiplication
 615                         for(i=1; i<=npts; i++)
 616                         {
 617                                 temp = nbasis[i] * b[jcount];
 618                                 p[icount + j] = p[icount + j] + temp;
 619                                 jcount = jcount + 3;
 620                         }
 621                 }
 622
 623                 icount = icount + 3;
 624                 t = t + step;
 625         }
 626
 627         delete[] x;
 628         delete[] nbasis;
 629 }
 630
 631 void Spline::knotu(int num, int order, int knotVector[])
 632 {
 633         int nplusc = num + order;
 634         int nplus2 = num + 2;
 635         knotVector[1] = 0;
 636
 637         for(int i=2; i<=nplusc; i++)
 638                 knotVector[i] = i - 1;
 639 }
 640
 641 void Spline::rbsplinu(int npts, int k, int p1, double b[], double h[], double p[])
 642 {
 643         int i, j, icount, jcount;
 644         int i1;
 645         //int x[30];            /* allows for 20 data points with basis function of order 5 */
 646         int nplusc;
 647
 648         double step;
 649         double t;
 650         //double nbasis[20];
 651         double temp;
 652
 653         nplusc = npts + k;
 654
 655         int * x = new int[nplusc + 1];
 656         double * nbasis = new double[npts + 1];
 657
 658         /*  zero and redimension the knot vector and the basis array */
 659
 660         for(i=0; i<=npts; i++)
 661                 nbasis[i] = 0.0;
 662
 663         for(i=0; i<=nplusc; i++)
 664                 x[i] = 0;
 665
 666         /* generate the uniform periodic knot vector */
 667
 668         knotu(npts, k, x);
 669
 670         /*
 671                 printf("The knot vector is ");
 672                 for (i = 1; i <= nplusc; i++){
 673                         printf(" %d ", x[i]);
 674                 }
 675                 printf("\n");
 676
 677                 printf("The usable parameter range is ");
 678                 for (i = k; i <= npts+1; i++){
 679                         printf(" %d ", x[i]);
 680                 }
 681                 printf("\n");
 682         */
 683
 684         icount = 0;
 685
 686         /*    calculate the points on the rational B-spline curve */
 687
 688         t = k - 1;
 689         step = ((double)((npts) - (k - 1))) / ((double)(p1 - 1));
 690
 691         for(i1=1; i1<=p1; i1++)
 692         {
 693                 if ((double)x[nplusc] - t < 5e-6)
 694                         t = (double)x[nplusc];
 695
 696                 rbasis(k, t, npts, x, h, nbasis);      /* generate the basis function for this value of t */
 697                 /*
 698                                 printf("t = %f \n",t);
 699                                 printf("nbasis = ");
 700                                 for (i = 1; i <= npts; i++){
 701                                         printf("%f  ",nbasis[i]);
 702                                 }
 703                                 printf("\n");
 704                 */
 705                 for(j=1; j<=3; j++)      /* generate a point on the curve */
 706                 {
 707                         jcount = j;
 708                         p[icount + j] = 0.0;
 709
 710                         for(i=1; i<=npts; i++)  /* Do local matrix multiplication */
 711                         {
 712                                 temp = nbasis[i] * b[jcount];
 713                                 p[icount + j] = p[icount + j] + temp;
 714                                 /*
 715                                                                 printf("jcount,nbasis,b,nbasis*b,p = %d %f %f %f %f\n",jcount,nbasis[i],b[jcount],temp,p[icount+j]);
 716                                 */
 717                                 jcount = jcount + 3;
 718                         }
 719                 }
 720                 /*
 721                                 printf("icount, p %d %f %f %f \n",icount,p[icount+1],p[icount+2],p[icount+3]);
 722                 */
 723                 icount = icount + 3;
 724                 t = t + step;
 725         }
 726
 727         delete[] x;
 728         delete[] nbasis;
 729 }
 730
 731 /**
 732  * Dumps the spline's data to stdout.
 733  */
 734 std::ostream & operator<<(std::ostream & os, const Spline & l)
 735 {
 736         os << " Spline: " << l.getData() << "\n";
 737         return os;
 738 }