← Back to branch summary

~librecad-dev/librecad/librecad

« back to all changes in this revision

Viewing changes to librecad/src/lib/engine/rs_spline.cpp

  • Committer: Scott Howard
  • Date: 2014-02-21 19:07:55 UTC
  • Revision ID: showard@debian.org-20140221190755-csjax9wb146hgdq4
first commit

Show diffs side-by-side

added added

removed removed

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