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

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