Repository navigation
Expand file tree
/
Copy pathnurbs.scad
More file actions
5651 lines (5280 loc) · 310 KB
/
Copy pathnurbs.scad
File metadata and controls
5651 lines (5280 loc) · 310 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/////////////////////////////////////////////////////////////////////
// LibFile: nurbs.scad
// B-Splines and Non-uniform Rational B-Splines (NURBS) are a way to represent smooth curves and smoothly curving
// surfaces with a set of control points. The curve or surface is defined by
// the control points and a set of "knot" points. The NURBS can be "clamped" in which case the curve passes through
// the first and last point, or they can be "closed" in which case the first and last point are coincident. Also possible
// are "open" curves which do not necessarily pass through any of their control points. Unlike Bezier curves, a NURBS
// can have an unlimited number of control points and changes to the control points only affect the curve locally.
//
// Includes:
// include <BOSL2/std.scad>
// include <BOSL2/nurbs.scad>
// FileGroup: Advanced Modeling
// FileSummary: NURBS and B-spline curves and surfaces.
//////////////////////////////////////////////////////////////////////
_BOSL2_NURBS = is_undef(_BOSL2_STD) && (is_undef(BOSL2_NO_STD_WARNING) || !BOSL2_NO_STD_WARNING) ?
echo("Warning: nurbs.scad included without std.scad; dependencies may be missing\nSet BOSL2_NO_STD_WARNING = true to mute this warning.") true : true;
// Section: NURBS Curves
// Private helper: compute the internal control list, multiplicities,
// expanded knot vector, and public-u domain mapping used by nurbs_curve().
// Returns [control2, mult2, knot2, bound].
function _nurbs_param_info(control, degree, mult, type="clamped", knots) =
let(
uniform = is_undef(knots),
dum = assert(in_list(type, ["closed","open","clamped"]), str("Unknown nurbs spline type: ", type))
assert(type=="closed" || len(control)>=degree+1, str(type," nurbs requires at least degree+1 control points"))
assert(is_undef(mult) || is_vector(mult), "mult must be a vector")
0,
badmult = is_undef(mult) ? []
: [for(i=idx(mult)) if (!(
is_int(mult[i])
&& mult[i]>0
&& (mult[i]<=degree
|| (type!="closed"
&& mult[i]==degree+1
&& (i==0 || i==len(mult)-1)
)
)
)) i],
dummy0 = assert(badmult==[], str("mult vector should contain positive integers no larger than the degree, except at ends of open splines, ",
"where degree+1 is allowed. The mult vector has bad values at indices: ",badmult))
assert(is_undef(knots) || is_undef(mult) || len(mult)==len(knots), "If both mult and knots are given they must be vectors of the same length")
assert(is_undef(mult) || type!="clamped" || sum(mult)==len(control)-degree+1,
str("For clamped spline knot count (sum of multiplicity vector) must be ",len(control)-degree+1," but is instead ",
is_undef(mult)?0:sum(mult)))
assert(is_undef(mult) || type!="closed" || sum(mult)==len(control)+1,
str("For closed spline knot count (sum of multiplicity vector) must be ",len(control)+1," but is instead ",
is_undef(mult)?0:sum(mult)))
assert(is_undef(mult) || type!="open" || sum(mult)==len(control)+degree+1,
str("For open spline knot count (sum of multiplicity vector) must be ",len(control)+degree+1," but is instead ",
is_undef(mult)?0:sum(mult)))
assert(uniform || is_increasing(knots), "Knot vector must be increasing")
0,
control2 = type=="open" ? control
: type=="clamped" ? control
: concat(control, select(control,count(degree))),
mult2 = !uniform ? mult
: type=="clamped" ? assert(is_undef(mult) || mult[0]==1 && last(mult)==1,"For clamped b-splines, first and last multiplicity must be 1")
[degree+1,each slice(default(mult, repeat(1,len(control2)-degree+1)),1,-2),degree+1]
: is_undef(mult) ? repeat(1,len(control2)+degree+1)
: type=="open" ? mult
: let(
lastmult = last(mult)+mult[0]-1,
dummy=assert(lastmult<=degree, "For closed spline, first and last knot multiplicity cannot total more than the degree+1"),
adjlast = [each select(mult,0,-2), lastmult]
)
_extend_knot_mult(adjlast,1,len(control2)+degree+1),
// IMPORTANT: this must use mult2, not the raw input mult.
knot2 = uniform && is_undef(mult2) ? lerpn(0,1,len(control2)+degree+1)
: uniform ? [for(i=idx(mult2)) each repeat(i/(len(mult2)-1),mult2[i])]
: let(
xknots = is_undef(mult2)? knots
: assert(len(mult2) == len(knots), "Knot vector and mult vector must be the same length")
[for(i=idx(mult2)) each repeat(knots[i], mult2[i])]
)
type=="open" ? assert(len(xknots)==len(control2)+degree+1, str("For open spline, knot vector with multiplicity must have length ",
len(control2)+degree+1," but has length ", len(xknots)))
xknots
: type=="clamped" ? assert(len(xknots) == len(control2)+1-degree, str("For clamped spline of degree ",degree,", knot vector with multiplicity must have length ",
len(control2)+1-degree," but has length ", len(xknots)))
assert(xknots[0]!=xknots[1] && last(xknots)!=select(xknots,-2),
"For clamped spline, first and last knots cannot repeat (must have multiplicity one)")
concat(repeat(xknots[0],degree), xknots, repeat(last(xknots),degree))
: assert(len(xknots) == len(control2)+1-degree, str("For closed spline, knot vector (including multiplicity) must have length ",
len(control2)+1-degree," but has length ", len(xknots),control2))
let(gmult=_calc_mult(xknots))
assert(gmult[0]+last(gmult)<=degree+1, "For closed spline, first and last knot multiplicity together cannot total more than the degree+1")
_extend_knot_vector(xknots,0,len(control2)+degree+1),
bound = type=="clamped" ? undef
: [knot2[degree], knot2[len(control2)]]
)
[control2, mult2, knot2, bound];
// Private helper: return the public-u breakpoints from _nurbs_param_info().
function _nurbs_breaks(degree, pinfo) =
let(
control = pinfo[0],
knot = pinfo[2],
bound = pinfo[3],
tvals = [
for (i=[degree:1:len(control)-1])
if (!approx(knot[i],knot[i+1])) knot[i],
knot[len(control)]
]
)
is_undef(bound) ? tvals
: [for (t=tvals) (t-bound[0])/(bound[1]-bound[0])];
// Function: nurbs_curve_breaks()
// Synopsis: Returns the parameter-space breakpoints of a NURBS curve.
// Topics: NURBS Curves
// See Also: nurbs_curve(), nurbs_length()
// Usage:
// breaks = nurbs_curve_breaks(control, degree, [mult=], [type=], [knots=]);
// breaks = nurbs_curve_breaks(nurbs_param_list);
// Description:
// Returns the breakpoints of the curve in the public `u` parameter space `[0,1]`.
// These are the parameter values where knot spans meet. For `type="closed"` both
// `0` and `1` are returned; they represent the same geometric point at the seam.
// Arguments:
// control = list of control points in any dimension or a NURBS parameter list
// degree = degree of NURBS
// ---
// mult = list of multiplicities of the knots. Default: all 1
// type = One of "clamped", "closed" or "open" to define end point handling of the spline. Default: "clamped"
// knots = List of knot values. Default: uniform
// Example(2D,NoAxes): Show the breakpoints of a curve in parameter space
// pts = [[13,43],[30,52],[49,22],[24,3]];
// echo(nurbs_curve_breaks(pts, 2));
function nurbs_curve_breaks(control, degree, mult, type="clamped", knots) =
is_list(control) && in_list(control[0], ["closed","open","clamped"]) ?
assert(len(control)>=6, "Invalid NURBS parameter list")
assert(num_defined([degree,mult,knots])==0,
"Cannot give degree, mult or knots when you provide a NURBS parameter list")
nurbs_curve_breaks(control[2], control[1], mult=control[4], type=control[0], knots=control[3])
: _nurbs_breaks(degree, _nurbs_param_info(control, degree, mult, type, knots));
// Function: nurbs_curve()
// Synopsis: Computes one or more points (or derivatives) on a NURBS curve.
// SynTags: Path
// Topics: NURBS Curves
// See Also: debug_nurbs()
// Usage:
// pts = nurbs_curve(control, degree, splinesteps, [mult=], [weights=], [type=], [knots=], [close_loop=]);
// pts = nurbs_curve(control, degree, u=, [mult=], [weights=], [type=], [knots=], [close_loop=]);
// dpts = nurbs_curve(control, degree, splinesteps, deriv=d, [two_sided=], [close_loop=], ...);
// dpts = nurbs_curve(control, degree, u=, deriv=d, [two_sided=], [close_loop=], ...);
// list = nurbs_curve(control, degree, splinesteps, deriv=[d1,d2,...], [two_sided=], [close_loop=], ...);
// list = nurbs_curve(control, degree, u=, deriv=[d1,d2,...], [two_sided=], [close_loop=], ...);
// Description:
// Compute the points specified by a NURBS curve. You specify the NURBS by supplying the control points, knots and weights.
// Only the control points are required. The knots and weights default to uniform, in which case you get a uniform B-spline.
// The length of `weights`, if given, must match the length of `control`.
// You can specify endpoint behavior using the `type` parameter. The default, "clamped", gives a curve which starts and
// ends at the first and last control points and moves in the tangent direction to the first and last control point segments.
// A "closed" curve is a one that starts where it ends. An "open" spline is a generic curve that starts somewhere
// in the middle of the control points. The "open" curve is less common; you only need this if you are managing the
// knots and control points yourself to create your own clamped or closed curve, so avoid this type unless you know what you're doing.
// Each of these types of curve require a different number of knots as described below.
// .
// The control points are the most important control over the shape
// of the curve. You must have at least degree+1 control points for clamped and open NURBS. Don't confuse the degree of a
// NURBS with its *order*: the order of a NURBS, often called $p$, is degree+1. Unlike a bezier, there is no maximum
// number of control points. A single NURBS is more like a bezier **path** than like a single bezier spline.
// .
// A NURBS or B-spline is a curve made from a moving average of several Bezier curves. The knots specify when one Bezier fades
// away to be replaced by the next one. The knot list is a non-decreasing list of values that you specify using two parameters,
// `knots` and `mult`. In practice changing the knot values doesn't have a strong effect on the curve, so it usually suffices
// to use a uniform knot vector, which is the default. The major exception to this is repeated knot values.
// At generic points in the NURBS, the curve is infinitely differentiable, but at a point that
// corresponds to a knot, a NURBS with degree $d$ will have a $(d-1)\mathrm{th}$ derivative that is continuous.
// However, if a value repeats in the knot vector that creates a knot with a multiplicity larger than 1, and each
// repetition decreases the smoothness of the curve at the corresponding NURBS point by 1. This means that
// if the multiplicity equals the degree then the curve is not differentiable: it has a corner at the knot point. Using the `mult` parameter
// without giving `knots` allows you to give a vector of multiplicities, which produces a knot vector that is uniform
// except it has some repeated knots. A value of 1 in the `mult` vector means the knot is not repeated; a value of 2 means it is
// repeated twice. The multiplicity can be as large as the degree but no larger. (A special exception is at the ends for open
// NURBS, where multiplicity degree+1 is permitted.) When you specify the multiplicity vector the total number of knots is the sum of that vector. You can also list
// the knots explicitly yourself. The knots exist in the parameter space of the NURBS, but the knot values you give can cover any range;
// they will be scaled to correspond properly to the NURBS parameter space: regardless of the knot values you give, the domain of evaluation
// for u is always the interval [0,1], and it will be scaled to give the entire valid portion of the curve you have chosen.
// .
// For an open spline the number of knots must be `len(control)+degree+1`. For a clamped spline the number of knots is `len(control)-degree+1`,
// and for a closed spline you need `len(control)+1` knots. If you are using the default uniform knots then the way to
// ensure that you have the right number is to check that mult is not set or `sum(mult)` equals the correct value.
// .
// You can use this function to evaluate the NURBS at `u`, which can be a single point or a list of points. You can also
// use it to evaluate the NURBS over its entire domain by giving a splinesteps value. This specifies the number of segments
// to use between each knot and guarantees a point exactly at each knot. This may be important if you set the knot multiplicity
// to the degree somewhere in your curve, which creates a corner at the knot, because it guarantees a sharp corner regardless
// of the number of points. If you don't give `u` or `splinesteps` then `splinesteps=16` is used as the default evaluation.
// .
// Instead of providing separate parameters you can give a first parameter of the form of a NURBS parameter list: `[type, degree, control, knots, mult, weights]`.
// .
// **Derivatives:** The `deriv` parameter requests curve derivatives in addition to, or instead of, curve points.
// - `deriv=0` (default) — returns the curve points, same as without the parameter. The output is a flat list of points, backward-compatible with code that does not use `deriv`.
// - `deriv=d` (positive integer) — returns a flat list of d-th derivative entries at each evaluation point, one entry per point.
// - `deriv=[d1,d2,...]` (list of integers) — returns a list of lists. Each element of the outer list corresponds to one entry of the `deriv` list, and contains the derivative entries of that order at every evaluation point. The output order matches the order of the `deriv` list, so `deriv=[0,1,2]` gives `[curve_pts, first_derivs, second_derivs]`.
// .
// The derivative order must be non-negative and cannot exceed the curve degree. When `u` is a single scalar and `deriv` is a list, the return value is a list of single entries (one per requested order) rather than a list of lists.
// .
// A NURBS curve need not be differentiable everywhere: at a knot of multiplicity m, derivatives
// of order larger than degree−m generally do not exist because the left and right limits differ.
// At a corner (m equal to the degree) already the first derivative jumps. Because `splinesteps`
// sampling places an evaluation point exactly at every knot, such points are hit routinely; in
// particular `deriv=degree` does not exist at any interior knot unless the curve happens to be
// smoother there than the multiplicity guarantees. The `two_sided` parameter selects what is
// returned at these points:
// - `two_sided=false` (default) — every derivative entry is a single vector. Where the requested
// derivative does not exist, the entry is `NAN`. You can detect these entries with `is_nan()`.
// Note that NAN propagates silently through further arithmetic, so check for it before using
// derivative values from a curve that may have corners.
// - `two_sided=true` — every derivative entry is a list: a singleton `[d]` containing the
// derivative vector at differentiable points, and a pair `[dleft, dright]` giving the one-sided
// left and right derivatives at non-differentiable points. Test `len(entry)==2` to detect the
// non-differentiable points.
// .
// Whether a derivative exists is decided by comparing the two one-sided values, not from the knot
// multiplicity alone, so a curve that is accidentally smoother than its knot multiplicity
// guarantees is treated as differentiable. At the domain endpoints of clamped and open curves
// only one side of the curve exists, and the one-sided derivative is returned as the derivative
// there (never NAN). Closed curves are periodic, so both sides exist everywhere, including at
// the seam. A clamped or open curve whose first and last control points coincide is also treated
// as periodic at its domain boundary — this is the case for a clamped curve built by gluing
// segments end to end (as in the circle example below), and for the curve
// {{nurbs_interp()}}(`closed=true`, `corners=`) returns, whose corners force a "clamped" type
// NURBS that nonetheless starts and ends at the same corner point. The `two_sided` parameter
// affects only derivative orders 1 and higher: order-0 entries (curve points) are always plain
// points.
// .
// For unweighted B-splines the derivatives are computed exactly using the difference-control-point method (Piegl & Tiller, "The NURBS Book", Algorithm A3.3). For rational NURBS (when `weights` is given) the geometric derivatives are obtained from the homogeneous B-spline derivatives via the quotient-rule formula (Piegl & Tiller, Eq. 4.8 / Algorithm A4.2).
// .
// **Closed curves and point repetition:** For `type="closed"` with `splinesteps`, the default behavior (`close_loop=false`) omits the final curve point because the curve is periodic: u=0 and u=1 map to the same point. This is consistent with BOSL2's general convention of not repeating closing points. Set `close_loop=true` to include the final point explicitly, making the output path suitable for direct use with `polygon()`. This parameter only affects the output when using `splinesteps`; it has no effect when specifying explicit `u` values. For derivatives with `two_sided=true`, the seam point uses wraparound logic to obtain two-sided derivatives, so no information is lost even though the final point is not repeated in the default output.
// Arguments:
// control = list of control points in any dimension or a NURBS parameter list
// degree = degree of NURBS
// splinesteps = evaluate whole spline with this number of segments between each pair of knots. Default: 16 if `u` is not given
// ---
// u = list of values or range in the interval [0,1] where the NURBS should be evaluated
// mult = list of multiplicities of the knots. Default: all 1
// weights = vector whose length is the same as control giving weights at each control point. Default: all 1
// type = One of "clamped", "closed" or "open" to define end point handling of the spline. Default: "clamped"
// knots = List of knot values. Default: uniform
// deriv = Integer or list of integers selecting which derivative orders to return. 0 = curve points. Default: 0
// two_sided = If true, return each derivative entry as a singleton `[d]` at differentiable points and as a pair `[dleft,dright]` of one-sided derivatives at non-differentiable points. If false, return plain derivative vectors, with NAN where the derivative does not exist. Default: false
// close_loop = For `type="closed"`, if true include the final curve point (which coincides with the first point, making the output path explicitly closed). If false, omit the final point. Default: false
// Example(2D,NoAxes): Compute some points and draw a curve and also some specific points:
// control = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// curve = nurbs_curve(control,2,splinesteps=16);
// pts = nurbs_curve(control,2,u=[0.4,0.8]);
// stroke(curve);
// color("red")move_copies(pts) circle(r=1.5,$fn=16);
// Example(2D,NoAxes): Compute NURBS points and make a polygon
// control = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// curve = nurbs_curve(control,2,splinesteps=16,type="closed");
// polygon(curve);
// Example(2D,NoAxes): Simple quadratic uniform clamped b-spline with some points computed using splinesteps.
// pts = [[13,43],[30,52],[49,22],[24,3]];
// debug_nurbs(pts,2);
// npts = nurbs_curve(pts, 2, splinesteps=3);
// color("red")move_copies(npts) circle(r=1);
// Example(2D,NoAxes): Simple quadratic uniform clamped b-spline with some points computed using the u parameter. Note that a uniform u parameter doesn't necessarily sample the curve uniformly.
// pts = [[13,43],[30,52],[49,22],[24,3]];
// debug_nurbs(pts,2);
// npts = nurbs_curve(pts, 2, u=[0:.2:1]);
// color("red")move_copies(npts) circle(r=1);
// Example(2D,NoAxes): Same control points, but cubic
// pts = [[13,43],[30,52],[49,22],[24,3]];
// debug_nurbs(pts,3);
// Example(2D,NoAxes): Same control points, quadratic and closed
// pts = [[13,43],[30,52],[49,22],[24,3]];
// debug_nurbs(pts,2,type="closed");
// Example(2D,NoAxes): Same control points, cubic and closed
// pts = [[13,43],[30,52],[49,22],[24,3]];
// debug_nurbs(pts,3,type="closed");
// Example(2D,NoAxes): Ten control points, quadratic, clamped
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// debug_nurbs(pts,2);
// Example(2D,NoAxes): Same thing, degree 4
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// debug_nurbs(pts,4);
// Example(2D,NoAxes): Same control points, degree 2, open. Note it doesn't reach the ends
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// debug_nurbs(pts,2, type="open");
// Example(2D,NoAxes): Same control points, degree 4, open. Note it starts farther from the ends
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// debug_nurbs(pts,4,type="open");
// Example(2D,NoAxes): Same control points, degree 2, closed
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// debug_nurbs(pts,2,type="closed");
// Example(2D,NoAxes): Same control points, degree 4, closed
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// debug_nurbs(pts,4,type="closed");
// Example(2D,Med,NoAxes): Adding weights
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// weights = [1,1,1,3,1,1,3,1,1,1];
// debug_nurbs(pts,4,type="clamped",weights=weights);
// Example(2D,NoAxes): Using knot multiplicity with quadratic clamped case. Knot count is len(control)-degree+1 = 9. The multiplicity 2 knot creates a corner for a quadratic.
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// mult = [1,1,1,2,1,1,1,1];
// debug_nurbs(pts,2,mult=mult,show_knots=true);
// Example(2D,NoAxes): Using knot multiplicity with quadratic clamped case. Two knots of multiplicity 2 gives two corners
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// mult = [1,1,1,2,2,1,1];
// debug_nurbs(pts,2,mult=mult,show_knots=true);
// Example(2D,NoAxes): Using knot multiplicity with cubic clamped case. Knot count is now 8. We need multiplicity equal to degree (3) to create a corner.
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// mult = [1,3,1,1,1,1];
// debug_nurbs(pts,3,mult=mult,show_knots=true);
// Example(2D,NoAxes): Using knot multiplicity with cubic closed case. Knot count is now len(control)+1=11. Here are three corners.
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// mult = [1,3,1,3,3];
// debug_nurbs(pts,3,mult=mult,type="closed",show_knots=true);
// Example(2D,NoAxes): Explicitly specified knots only change the quadratic clamped curve slightly. Knot count is len(control)-degree+1 = 9.
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// knots = [0,1,3,5,9,13,14,19,21];
// debug_nurbs(pts,2,knots=knots);
// Example(2D,NoAxes): Combining explicit knots with mult for the quadratic curve to add a corner
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// knots = [0,1,3,9,13,14,19,21];
// mult = [1,1,1,2,1,1,1,1];
// debug_nurbs(pts,2,knots=knots,mult=mult);
// Example(2D,NoAxes): Directly repeating a knot in the knot list to create a corner for a cubic spline
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// knots = [0,1,3,13,13,13,19,21];
// debug_nurbs(pts,3,knots=knots);
// Example(2D,NoAxes): Open cubic spline with explicit knots
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// knots = [0,1,3,13,13,13,19,21,27,28,29,40,42,44];
// debug_nurbs(pts,3,knots=knots,type="open");
// Example(2D,NoAxes): Closed quintic spline with explicit knots
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// knots = [0,1,3,13,13,13,19,21,27,28,33];
// debug_nurbs(pts,5,knots=knots,type="closed");
// Example(2D,Med,NoAxes): Closed quintic spline with explicit knots and weights
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// weights = [1,2,3,4,5,6,7,6,5,4];
// knots = [0,1,3,13,13,13,19,21,27,28,33];
// debug_nurbs(pts,5,knots=knots,weights=weights,type="closed");
// Example(2D,NoAxes): Circular arcs are possible with NURBS. This example gives a semi-circle
// control = [[1,0],[1,2],[-1,2],[-1,0]];
// w = [1,1/3,1/3,1];
// debug_nurbs(control, 3, weights=w, width=0.1, size=.2);
// Example(2D,NoAxes): Gluing two semi-circles together gives a whole circle. Note that this is a clamped not a closed NURBS. The interface uses a knot of multiplicity 3 where the clamped ends of the semi-circles meet.
// control = [[1,0],[1,2],[-1,2],[-1,0],[-1,-2],[1,-2],[1,0]];
// w = [1,1/3,1/3,1,1/3,1/3,1];
// debug_nurbs(control, 3, splinesteps=16,weights=w,mult=[1,3,1],width=.1,size=.2);
// Example(2D,NoAxes): Circle constructed with type="closed"
// control = [[1,0],[1,2],[-1,2],[-1,0],[-1,-2],[1,-2]];
// w = [1,1/3,1/3,1,1/3,1/3];
// debug_nurbs(control, 3, splinesteps=16,weights=w,mult=[1,3,3],width=.1,size=.2,type="closed",show_knots=true);
// Example(2D,NoAxes): One-sided derivatives at a corner. This quadratic has a corner at the multiplicity-2 knot (u=3/7). The default `two_sided=false` returns NAN for the first derivative there; `two_sided=true` returns both one-sided tangents, drawn here as arrows into (red) and out of (blue) the corner.
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// mult = [1,1,1,2,1,1,1,1];
// stroke(nurbs_curve(pts,2,mult=mult,splinesteps=16));
// corner = nurbs_curve(pts,2,mult=mult,u=3/7);
// tang = nurbs_curve(pts,2,mult=mult,u=3/7,deriv=1,two_sided=true);
// color("red") stroke([corner-15*unit(tang[0]), corner], endcap2="arrow2", width=1);
// color("blue") stroke([corner, corner+15*unit(tang[1])], endcap2="arrow2", width=1);
function nurbs_curve(control,degree,splinesteps,u, mult,weights,type="clamped",knots,deriv=0,two_sided=false,close_loop=false) =
let(
splinesteps = !any_defined([splinesteps,u]) ? 16 : splinesteps
)
is_list(control) && in_list(control[0], ["closed","open","clamped"]) ?
assert(len(control)>=6, "Invalid NURBS parameter list")
assert(num_defined([degree,mult,weights,knots])==0,
"Cannot give degree, mult, weights or knots when you provide a NURBS parameter list")
nurbs_curve(control[2], control[1], splinesteps, u, weights=control[5],mult=control[4], type=control[0], knots=control[3], deriv=deriv, two_sided=two_sided, close_loop=close_loop)
: assert(num_defined([splinesteps,u])==1, "Must define exactly one of u and splinesteps")
is_finite(u) ?
let(r=nurbs_curve(control,degree,u=[u],mult=mult,weights=weights,knots=knots,type=type,deriv=deriv,two_sided=two_sided,close_loop=close_loop))
!is_list(deriv) ? r[0] : [for(d=r) d[0]]
: assert(is_undef(splinesteps) || (is_int(splinesteps) && splinesteps>0), "splinesteps must be a positive integer")
let(u=is_range(u) ? list(u) : u)
assert(is_undef(u) || (is_vector(u) && min(u)>=0 && max(u)<=1), "u must be a list of points on the interval [0,1] or a range contained in that interval")
assert(is_bool(two_sided), "two_sided must be a boolean")
is_def(weights) ? assert(is_vector(weights, len(control)), "Weights should be a vector whose length is the number of control points")
let(
dim = len(control[0]),
hctrl = [for(i=idx(control)) [each control[i]*weights[i],weights[i]]]
)
deriv == 0 ?
let(curve = nurbs_curve(hctrl,degree,u=u,splinesteps=splinesteps, mult=mult, knots=knots, type=type, close_loop=close_loop))
[for(pt=curve) select(pt,0,-2)/last(pt)]
:
assert(is_list(deriv) ? min(deriv)>=0 : deriv>=0, "Derivative orders must be non-negative")
assert(is_list(deriv) ? max(deriv)<=degree : deriv<=degree,
str("Derivative order exceeds curve degree ",degree))
let(
max_d = is_list(deriv) ? max(deriv) : deriv,
// Request both one-sided homogeneous derivatives so the geometric
// left/right derivatives can be formed at points where the curve
// may not be differentiable. Entries of hderivs[k] for k>=1 are
// singletons [v] or pairs [left,right]; hderivs[0] is plain points.
hderivs = nurbs_curve(hctrl, degree, u=u, splinesteps=splinesteps,
mult=mult, knots=knots, type=type,
deriv=[for(k=[0:1:max_d]) k], two_sided=true, close_loop=close_loop),
n_u = len(hderivs[0]),
geo_lr = [for(i=[0:1:n_u-1])
let(
anypair = [for(k=[1:1:max_d]) if (len(hderivs[k][i])==2) 1] != [],
geo_l = _rational_derivs_at_u(
[hderivs[0][i], for(k=[1:1:max_d]) hderivs[k][i][0]],
dim, max_d),
geo_r = !anypair ? geo_l
: _rational_derivs_at_u(
[hderivs[0][i], for(k=[1:1:max_d]) last(hderivs[k][i])],
dim, max_d)
)
[geo_l, geo_r]
]
)
is_list(deriv) ?
[for(d=deriv) [for(i=[0:1:n_u-1]) _sided_out(geo_lr[i][0][d], geo_lr[i][1][d], d, two_sided)]]
:
[for(i=[0:1:n_u-1]) _sided_out(geo_lr[i][0][deriv], geo_lr[i][1][deriv], deriv, two_sided)]
:
let(
uniform = is_undef(knots),
pinfo = _nurbs_param_info(control, degree, mult, type, knots),
control = pinfo[0],
mult = pinfo[1],
knot = pinfo[2],
bound = pinfo[3],
dummy_deriv = is_list(deriv) ?
assert(min(deriv)>=0, "Derivative orders must be non-negative")
assert(max(deriv)<=degree, str("Maximum derivative order ",max(deriv)," exceeds curve degree ",degree))
0
: deriv != 0 ?
assert(deriv>=0, "Derivative order must be non-negative")
assert(deriv<=degree, str("Derivative order ",deriv," exceeds curve degree ",degree))
0 : 0,
adjusted_u_orig = !is_undef(splinesteps) ?
[for(i=[degree:1:len(control)-1])
each
if (!approx(knot[i],knot[i+1]))
lerpn(knot[i],knot[i+1],splinesteps, endpoint=false),
if (type!="closed" || close_loop) knot[len(control)]
]
: is_undef(bound) ? u
: add_scalar((bound[1]-bound[0])*u,bound[0]),
reorder = is_undef(splinesteps) && !is_increasing(adjusted_u_orig) ?
let(ind = sortidx(adjusted_u_orig))
[ind,sortidx(ind)]
: false,
// The u list needs to be sorted for the algorithm to identify the knot spans, so sort it if necessary
adjusted_u = reorder ? select(adjusted_u_orig,reorder[0]) : adjusted_u_orig,
nurbs_pts =
deriv != 0 ?
let(
msum_d = uniform ? cumsum(mult) : undef,
kmult_d = !uniform ? _calc_mult(knot) : undef,
knotidx_list = uniform ?
[for(uval=adjusted_u)
let(
mind = floor(uval*(len(mult)-1)),
knotidxR = msum_d[mind]-1
)
knotidxR<len(control) ? knotidxR : knotidxR - mult[mind]
]
: [for(
kind=kmult_d[0]-1, uind=0, kmultind=1, output=undef, done=false
;
!done
;
output = (uind<len(adjusted_u) && approx(adjusted_u[uind],knot[kind]) && kind>kmult_d[0]-1
&& ((kmultind>=len(kmult_d)-1 || kind+kmult_d[kmultind]>=len(control))))
?kind-kmult_d[kmultind-1]
: (uind<len(adjusted_u) && adjusted_u[uind]>=knot[kind] && adjusted_u[uind]>=knot[kind]
&& adjusted_u[uind]<knot[kind+kmult_d[kmultind]]) ? kind
: undef,
done = uind==len(adjusted_u),
uind = is_def(output) ? uind+1 : uind,
inc_k = uind<len(adjusted_u) && adjusted_u[uind]>=knot[kind+kmult_d[kmultind]],
kind = inc_k ? kind+kmult_d[kmultind] : kind,
kmultind = inc_k ? kmultind+1 : kmultind
)
if (is_def(output)) output]
)
// Chain rule: t = bound[0] + (bound[1]-bound[0])*u, so
// d^k/du^k = (bound[1]-bound[0])^k * d^k/dt^k.
// For clamped, bound=undef and the mapping is identity (no scaling needed).
let(
// Knot-run info for every evaluation point: [first_index, multiplicity]
// when the point lies (approx) on a knot, undef otherwise. Used to
// detect points where the requested derivative order may not exist.
runinfo = [for(uval=adjusted_u) _knot_run_info(knot, uval)]
)
is_list(deriv) ?
[for(d=deriv)
let(scale = is_undef(bound) || d==0 ? 1 : pow(bound[1]-bound[0], d))
[for(i=idx(adjusted_u))
let(base = scale * _nurbs_eval_deriv(knot,
select(control, knotidx_list[i]-degree, knotidx_list[i]),
adjusted_u[i], knotidx_list[i], degree, d))
d==0 ? base
: _deriv_sided(base, d, runinfo[i], adjusted_u[i], knot, control,
degree, type, scale, two_sided)
]
]
: // integer deriv > 0: flat list of derivative entries
let(scale = is_undef(bound) ? 1 : pow(bound[1]-bound[0], deriv))
[for(i=idx(adjusted_u))
let(base = scale * _nurbs_eval_deriv(knot,
select(control, knotidx_list[i]-degree, knotidx_list[i]),
adjusted_u[i], knotidx_list[i], degree, deriv))
_deriv_sided(base, deriv, runinfo[i], adjusted_u[i], knot, control,
degree, type, scale, two_sided)
]
: uniform?
let(
msum = cumsum(mult)
)
[for(uval=adjusted_u)
let(
mind = floor(uval*(len(mult)-1)),
knotidxR=msum[mind]-1,
knotidx = knotidxR<len(control) ? knotidxR : knotidxR - mult[mind]
)
_nurbs_pt(knot,select(control,knotidx-degree,knotidx),uval,1,degree,knotidx)
]
: let(
kmult = _calc_mult(knot),
knotidx =
[for(
kind = kmult[0]-1,
uind=0,
kmultind=1,
output=undef,
done=false
;
!done
;
output = (uind<len(adjusted_u) && approx(adjusted_u[uind],knot[kind]) && kind>kmult[0]-1
&& ((kmultind>=len(kmult)-1 || kind+kmult[kmultind]>=len(control))))
?kind-kmult[kmultind-1]
: (uind<len(adjusted_u) && adjusted_u[uind]>=knot[kind] && adjusted_u[uind]>=knot[kind]
&& adjusted_u[uind]<knot[kind+kmult[kmultind]]) ? kind
: undef,
done = uind==len(adjusted_u),
uind = is_def(output) ? uind+1 : uind,
inc_k = uind<len(adjusted_u) && adjusted_u[uind]>=knot[kind+kmult[kmultind]],
kind = inc_k ? kind+kmult[kmultind] : kind,
kmultind = inc_k ? kmultind+1 : kmultind
)
if (is_def(output)) output]
)
[for(i=idx(adjusted_u))
_nurbs_pt(knot,slice(control, knotidx[i]-degree,knotidx[i]), adjusted_u[i], 1, degree, knotidx[i])
]
)
!reorder ? nurbs_pts
: !is_list(deriv) ? select(nurbs_pts,reorder[1])
: [for(dr=nurbs_pts) select(dr,reorder[1])];
function _nurbs_pt(knot, control, u, r, p, k) =
r>p ? control[0]
:
let(
ctrl_new = [for(i=[k-p+r:1:k])
let(
alpha = (u-knot[i]) / (knot[i+p-r+1]-knot[i])
)
(1-alpha) * control[i-1-(k-p)-r+1] + alpha*control[i-(k-p)-r+1]
]
)
_nurbs_pt(knot,ctrl_new,u,r+1,p,k);
// Piegl & Tiller Algorithm A3.3: build difference control points for the d-th derivative.
// Each call reduces the local control polygon by one point and increments d_so_far.
// A zero knot span (repeated knot) produces a zero vector, giving zero derivative there.
function _bspline_deriv_cpts(local_ctrl, knot, k, p, d_so_far, target_d) =
d_so_far >= target_d ? local_ctrl
: let(
step = d_so_far + 1,
coef = p - d_so_far,
new_cpts = [for(i=[0:1:len(local_ctrl)-2])
let(denom = knot[k+i+1] - knot[k-p+i+step])
denom == 0 ? 0*local_ctrl[0]
: coef * (local_ctrl[i+1] - local_ctrl[i]) / denom
]
)
_bspline_deriv_cpts(new_cpts, knot, k, p, step, target_d);
// Evaluate the d-th derivative of a B-spline at u in knot span k.
// Reduces to a degree-(p-d) B-spline evaluation on the difference control points.
function _nurbs_eval_deriv(knot, local_ctrl, u, k, p, d) =
d == 0 ? _nurbs_pt(knot, local_ctrl, u, 1, p, k)
: _nurbs_pt(knot, _bspline_deriv_cpts(local_ctrl, knot, k, p, 0, d), u, 1, p-d, k);
// Locate t in the knot vector: returns [first_index, multiplicity] of the run of
// knots (approx) equal to t, or undef when t does not lie on a knot.
function _knot_run_info(knot, t) =
let(idxs = [for(i=idx(knot)) if (approx(knot[i], t)) i])
idxs == [] ? undef : [idxs[0], len(idxs)];
// Package the d-th derivative at one evaluation point from its one-sided values.
// Existence is decided by comparing the sides with a magnitude-relative tolerance
// (derivative magnitudes scale with the curve size, so the absolute EPSILON of
// approx() alone would misclassify large curves).
function _sided_out(left, right, d, two_sided) =
d == 0 ? left
: approx(left, right, EPSILON * (1 + norm(left) + norm(right)))
? (two_sided ? [left] : left)
: two_sided ? [left, right] : NAN;
// Resolve the d-th derivative at an evaluation point that may fall on a knot
// where that derivative does not exist.
//
// base = one-sided value the span logic produced (right-sided, except
// left-sided at the domain end), already chain-rule scaled.
// run = [first_knot_index, multiplicity] if the point lies on a knot, else undef.
//
// At a knot of multiplicity m the curve is guaranteed C^(degree-m), so for
// d <= degree-m the derivative exists and base is returned. Otherwise both
// one-sided derivatives are computed and compared: if they agree the curve is
// smoother there than the multiplicity bound guarantees and the common value is
// returned; if they disagree the result is NAN (two_sided=false) or the pair
// [left, right] (two_sided=true). At true clamped/open domain endpoints only
// one side exists and it is returned as the derivative; closed curves are
// periodic, so the side missing at one end of the seam is evaluated at the
// other end.
//
// A clamped/open curve whose first and last control points coincide is also
// periodic at its domain boundary even though its type is not "closed" --
// this is BOSL2's own documented pattern for gluing curve segments into a
// loop via a repeated endpoint (see the nurbs_curve() "gluing two
// semi-circles together" example), and it is exactly what
// nurbs_interp(closed=true, corners=...) produces: corners force a
// "clamped" assembly, but the assembled curve starts and ends at the same
// corner point, so u=0 and u=1 are the same point in space with a real
// corner there. The wraparound math below only depends on knot[degree]
// being the domain-start value and knot[len(control)] the domain-end value
// (both true for clamped and open knot vectors by construction, independent
// of any array extension), so the same formula that handles the "closed"
// type's pre-extended arrays works unchanged on the plain unextended arrays
// of such a loop.
function _deriv_sided(base, d, run, t, knot, control, degree, type, scale, two_sided) =
is_undef(run) || d <= degree - run[1] ? (two_sided ? [base] : base)
: let(
n_ctrl = len(control),
k_l = run[0] - 1,
k_r = run[0] + run[1] - 1,
valid_l = k_l >= degree,
valid_r = k_r <= n_ctrl - 1,
loop_ends = type != "closed"
&& approx(control[0], last(control),
EPSILON * (1 + norm(control[0]) + norm(last(control))))
)
type != "closed" && !loop_ends && (!valid_l || !valid_r) ? (two_sided ? [base] : base)
: let(
lrun = valid_l ? run : _knot_run_info(knot, knot[n_ctrl]),
lt = valid_l ? t : knot[n_ctrl],
lk = valid_l ? k_l : lrun[0] - 1,
rrun = valid_r ? run : _knot_run_info(knot, knot[degree]),
rt = valid_r ? t : knot[degree],
rk = valid_r ? k_r : rrun[0] + rrun[1] - 1,
left = scale * _nurbs_eval_deriv(knot, select(control, lk-degree, lk), lt, lk, degree, d),
right = scale * _nurbs_eval_deriv(knot, select(control, rk-degree, rk), rt, rk, degree, d)
)
_sided_out(left, right, d, two_sided);
// Binomial coefficient C(n,k), computed recursively via Pascal's triangle.
function _binom(n, k) =
k == 0 || k == n ? 1 : _binom(n-1, k-1) + _binom(n-1, k);
// Piegl & Tiller Eq. 4.8 / Algorithm A4.2: convert homogeneous B-spline derivatives to
// geometric NURBS derivatives at a single u. Aw[k] is the k-th homogeneous derivative
// (a dim+1 vector); returns [C^0, C^1, ..., C^max_d] in geometric (un-homogenised) space.
function _rational_derivs_at_u(Aw, dim, max_d, k=0, Cders=[]) =
k > max_d ? Cders
: let(
w0 = last(Aw[0]),
Ak = select(Aw[k], 0, dim-1),
corr = k == 0 ? repeat(0, dim)
: sum([for(j=[1:1:k]) _binom(k,j) * last(Aw[j]) * Cders[k-j]]),
Ck = (Ak - corr) / w0
)
_rational_derivs_at_u(Aw, dim, max_d, k+1, [each Cders, Ck]);
function _extend_knot_mult(mult, next, len) =
let(total = sum(mult))
total == len ? mult
: total>len ? [ each select(mult,0,-2), last(mult)-(total-len) ]
: _extend_knot_mult([each mult,mult[next]], next+1, len);
function _extend_knot_vector(knots,next,len) =
len(knots)==len ? knots
: _extend_knot_vector([each knots, last(knots)+knots[next+1]-knots[next]], next+1, len);
function _calc_mult(knots) =
let(
ind=[ 0,
for(i=[1:len(knots)-1])
if (!approx(knots[i],knots[i-1])) i,
len(knots)
]
)
deltas(ind);
// Module: debug_nurbs()
// Synopsis: Shows a NURBS curve and its control points, knots and weights
// SynTags: Geom
// Topics: NURBS, Debugging
// See Also: nurbs_curve()
// Usage:
// debug_nurbs(control, degree, [width], [splinesteps=], [type=], [mult=], [knots=], [size=], [show_weights=], [show_knots=], [show_index=]);
// Description:
// Displays a 2D or 3D NURBS and the associated control points to help debug NURBS curves. You can display the
// control point indices and weights, and can also display the knot points.
// Instead of providing separate parameters you can give a first parameter of the form of a NURBS parameter list: `[type, degree, control, knots, mult, weights]`.
// Arguments:
// control = list of control points in any dimension or a NURBS parameter list
// degree = degree of NURBS
// splinesteps = number of segments between each pair of knots. Default: 16
// width = width of the line. Default: 1
// size = size of text annotations. Default: 3 times the width
// mult = multiplicity vector for NURBS
// weights = weight vector for NURBS
// type = NURBS type, one of "clamped", "open" or "closed". Default: "clamped"
// show_index = if true then display index of each control point vertex. Default: true
// show_weights = if true then display any non-unity weights. Default: true if weights vector is supplied, false otherwise
// show_knots = If true then show the knots on the spline curve. Default: false
// show_control = If true then show the control points and its polygon. Default: true
//
// Example(2D,Med,NoAxes): If you want to see the knots set `show_knots=true`:
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// debug_nurbs(pts,4,type="clamped",show_knots=true);
//
// Example(2D,Med,NoAxes): Non-unity weights are displayed if you give a weight vector
// pts = [[5,0],[0,20],[33,43],[37,88],[60,62],[44,22],[77,44],[79,22],[44,3],[22,7]];
// weights = [1,1,1,7,1,1,7,1,1,1];
// debug_nurbs(pts,4,type="closed",weights=weights);
//
module debug_nurbs(control,degree,splinesteps=16,width=1, size, mult,weights,type="clamped",knots, show_weights, show_knots=false, show_index=true, show_control=true)
{
if (is_list(control) && in_list(control[0], ["closed","open","clamped"])) {
assert(len(control)>=6, "Invalid NURBS parameter list")
assert(num_defined([degree,mult,weights,knots])==0,
"Cannot give degree, mult, weights or knots when you provide a NURBS parameter list")
debug_nurbs(control[2], control[1], splinesteps, width, size, weights=control[5],mult=control[4], type=control[0], knots=control[3],
show_weights=show_weights, show_knots=show_knots, show_index=show_index, show_control=show_control);
}
else {
$fn=8;
size = default(size, 3*width);
show_weights = default(show_weights, is_def(weights));
N=len(control);
twodim = len(control[0])==2;
curve = nurbs_curve(control=control,degree=degree,splinesteps=splinesteps, mult=mult,weights=weights, type=type, knots=knots);
stroke(curve, width=width, closed=type=="closed");//, color="green");
if (show_control)
stroke(control, width=width/2, color="white", closed=type=="closed");
if (show_knots){
knotpts = nurbs_curve(control=control, degree=degree, splinesteps=1, mult=mult, weights=weights, type=type, knots=knots);
color([0,.8,0])
move_copies(knotpts)
if (twodim) union(){rect([3*width,width]); rect([width,3*width]);} //circle(r=width);
else for(i=[0:2]) cuboid(list_rotate([3*width,width,width],i));
}
color("blue")
if (show_index && show_control)
move_copies(control){
let(label = str($idx),
anch = show_weights && is_def(weights[$idx]) && weights[$idx]!=1 ? FWD : CENTER)
if (twodim) text(text=label, size=size, anchor=anch);
else rot($vpr) text3d(text=label, size=size, anchor=anch);
}
color("blue")
if ( show_weights)
move_copies(control){
if(is_def(weights[$idx]) && weights[$idx]!=1)
let(label = str("w=",weights[$idx]),
anch = show_index ? BACK : CENTER
)
if (twodim) text(text=label, size=size, anchor=anch);
else rot($vpr) text3d(text=label, size=size, anchor=anch);
}
}
}
// Function: nurbs_interp()
// Synopsis: Finds a NURBS curve passing through a point list with optional derivative constraints.
// Topics: NURBS Curves, Interpolation
// See Also: nurbs_curve(), debug_nurbs(), debug_nurbs_interp()
//
// Usage:
// nurbs_param = nurbs_interp(points, degree, [method=], [closed=], [start_deriv=], [end_deriv=], [curvature=], [start_curvature=], [end_curvature=], [corners=], [deriv=], [extra_pts=], [smooth=]);
//
// Description:
// Given a list of data points and a NURBS degree, computes a curve of the specified degree
// that passes exactly through every data point. The computed curve always has
// uniform weights, but irregularly spaced knots, so it is actually a non-uniform B-spline.
// Data points may 2D or any higher dimension. Returns a NURBS parameter list of the form
// `[type, degree, control_points, knots, undef, undef, u]` that can be
// passed directly to {{nurbs_curve()}} and other NURBS functions. The extra return value `u`,
// described in detail below, enables you to locate your input points in the computed spline
// .
// When `closed=false` (the default) the output is a "clamped" NURBS.
// When `closed=true`, the interpolation treats the data points as a loop and produces a
// curve that is smooth at the closing point. The output will be a "closed" NURBS (unless you
// specify corners as described below).
// If you instead duplicate the closing point and set `closed=false` then the
// result will have a corner at the closing point.
// .
// Inserting a corner converts the output from a closed NURBS to to a clamped NURBS. Adding
// more than one corner converts the output to a piecewise sequence of clamped NURBS.
// .
// **Parameterization** (`method=`)
// .
// In order to solve the interpolation problem, the algorithm first chooses
// the NURBS parameter value `u[k]` that will correspond to each `points[k]`.
// This parametrization step significantly affects the shape of the output curve, particularly when the
// data points are not evenly spaced. The following methods are supported:
// .
// - `"length"` — Base parameters values on the chord length, which is distance between the consecutive data points.
// Best when data points are fairly evenly spaced.
// - `"centripetal"` (default) — Base parameters values on the square root of the chord length. (Lee 1989).
// - `"dynamic"` — like centripetal, but the exponent 0.5 is replaced
// by a per-chord value chosen based on local spacing variation. Long chords
// get a smaller exponent and short chords a larger one, compressing the
// influence of outliers. Chord lengths are normalized, which makes the method scale
// invariant and prevents misbehavior at extreme scales. Scaling is not given in the original reference. (Balta et al. 2020).
// - `"foley"` — centripetal base, augmented by corrections at each point that
// are proportional to the local turn angle. Sharp bends pull parameter values
// closer together, which tends to reduce overshoot at corners (Foley & Neilson 1987).
// - `"fang"` — centripetal base, augmented by a correction based on the radius
// of the osculating circle at each point. Said to handles mixed straight-and-curved
// segments particularly well. This method is NOT scale invariant, so results will
// change if you scale your input data. (Fang & Hung 2013).
// .
// The other required input to the interpolation is the location of the knots.
// We place knots using a moving average of `degree` consecutive parameter values, which links
// the knots to the local parameter spacing. A consequence of this process for selection
// of the parameters and knot locations is that even if your input data has symmetry it is
// likely that the symmetry will be broken in the output. For closed curves, another
// consequence is that the resulting curve will depend on which point is chosen as the
// starting point for the interpolation. The algorithm chooses a starting point
// that is expected to provide the best behaved interpolation curve. Examining the
// knot positions with {{debug_nurbs_interp()}} may help you understand unexpected behavior
// you observe in the output. If your curve does not
// behave as desired you may be able to adjust it by imposing additional constraints or
// by giving it more freedom using `extra_pts`.
// .
// **Derivative constraints** (`deriv=`, `start_deriv=`, `end_deriv=`)
// .
// `deriv[k]` specifies the tangent direction and speed the curve must have
// as it passes through `points[k]`. The length of `deriv[k]` gives the speed
// as a multiple of `path_length(points)` which means a unit vector gives a natural
// speed that is a good starting point.
// The speed has a big effect on the shape of the curve, so if the local shape is
// not as you desire you should try increasing it, which will make the curve around
// the point flatter or decreasing it, which will make the curve more pointy.
// Set `deriv[k] = undef` to leave point `k` unconstrained.
// If you only want to set the derivative at the ends of a "clamped" curve you can use
// `start_deriv=` and `end_deriv=`, which set
// `deriv[0]` and `last(deriv)` without the need to provide a list of undefs for all the interior points.
// .
// **Curvature constraints** (`curvature=`, `start_curvature=`, `end_curvature=`)
// .
// The curvature at a point measures how tightly a curve bends.
// When a point has curvature $\kappa$ then a circle with radius $1/\kappa$
// locally matches the curve at that point so both its first and second derivatives agree.
// This matched circle is called the osculating circle. When you set `curvature[k]` this
// constrains the curvature at `points[k]`. Every curvature-constrained point **must** also have a derivative constraint
// at the same index. Curvature constraints require a degree of at least 2.
// .
// In general curvature constraints require the curvature **vector**, which
// points in the direction of the osculating circle and has length equal to the curvature.
// The curvature vector must be orthogonal to the tangent vector at the point;
// when you specify a curvature vector any component parallel to the tangent is removed.
// The magnitude of the curvature is taken as the magnitude of your original input vector,
// even if subtracting the tangent component changes its length.
// For 2D curves you can also provide curvature as a scalar, with the sign indicating direction.
// (positive = left/CCW, negative = right/CW).
// .
// You can specify the curvature at the ends of "clamped" curves using
// `start_curvature=` and `end_curvature=`, which specify `curvature[0]`
// and `last(curvature)` without the need to create undefs for all the interior points.
// .
// **Corners** (`corners=`)
// .
// `corners=` is a list of interior point indices where the curve has
// a corner, a discontinuity in the derivative. You can also specify a corner
// at point `k` by setting `deriv[k]=NAN`. When you request corners, the
// algorithm chops up the input data into separate clamped splines that run from corner
// to corner. When `closed=true` this results in a "clamped" output spline, and the curve
// will start at one of your corner points.
// If you place corners close together, the effective degree of the short segment
// in between the corners may be reduced. These curve sections are assembled into a single
// NURBS so this process is transparent to the user. A limitation is that you cannot control
// the derivatives of the two segments that meet at a corner. If you need to do this you
// must construct your own sequence of clamped interpolations.
// .
// **Extra control points** (`extra_pts=`, `smooth=`)
// .
// By default, the solver uses exactly as many control points as are needed to
// satisfy the interpolation and constraint conditions, which gives a unique
// solution. This unique solution may be badly behaved, with undesirable oscillations.
// You can improve the behavior by requesting extra points.
// Specifying `extra_pts=N` inserts `N` additional control points and knots, making the
// system underdetermined: infinitely many curves pass through the data points and satisfy
// the constraints. The solver picks the one that satisfies
// a smoothness criterion specified by `smooth=`:
// .
// - `smooth=1` — minimises the sum of squared differences between consecutive
// control points. This tends to keep the control polygon short and reduces
// large-scale variation in the curve.
// - `smooth=2` — minimises the sum of squared second differences of the control
// points. This penalises bending in the control polygon, generally producing
// a fairer, less wiggly curve than `smooth=1`.
// - `smooth=3` (default) — minimises the integrated squared second derivative
// $\int \|\mathbf{C}''(t)\|^2 \, dt$, often called the *bending energy* of
// the curve. Unlike `smooth=2`, which only looks at the control polygon,
// this criterion acts directly on the curve shape and is the most
// mathematically principled choice for smooth interpolation. Requires
// `degree >= 2`.
// .
// The number of extra control points cannot exceed the number of knot spans.
// If you request too many, the number is capped and a warning is displayed.
// With `corners=`, the curve is split into independent clamped segments and
// the extra points are distributed across eligible segments proportionally
// to their control-point count, rounding up, so the total may
// exceed the requested number but will never be less. A segment is eligible when
// its effective degree is 3 or higher, or when it is degree 2 with `smooth=1`.
// .
// **Locating points in the spline** — In order to locate your original data
// points in the spline you need the `u` parameter value that you
// can pass to {{nurbs_curve()}}. The last return value `u` is a list
// where `u[k]` is the NURBS parameter at which the curve passes through
// `points[k]`.
// .
// **Smoothness** — The smoothness of B-splines is determined by the
// degree. If you request a degree $p$ spline then it will be $C^{p-1}$ at
// knot points and $C^\infty$ everywhere else. If you request corners then
// these are points where the curve is not differentiable; corners may
// also divide the curve into small segments that lack sufficient points
// to support an interpolation at your requested degree: a degree $p$ interpolation
// requires $p+1$ points. In this case, the interpolation is performed at a lower
// degree and elevated, which means it will be less smooth at knots.
//
// Arguments:
// points = List of data points to interpolate (2D or any higher dimension).
// degree = Degree of the NURBS. Degree 3 (cubic) is the most common choice.
// ---
// method = Parameterization method: `"length"`, `"centripetal"`, `"dynamic"`, `"foley"`, or `"fang"`. Default: `"centripetal"`
// closed = If true treat point list as a loop . Default: `false`
// start_deriv = If `closed=false`, gives the tangent vector at the first point
// end_deriv = If `closed=false`, gives tangent vector at the last point.
// deriv = List of tangent vector constraints for every point, NAN at corners or undef at unconstrained points. Cannot be combined with `start_deriv=`/`end_deriv=`.
// start_curvature = If `closed=false` gives curvature at first point. (Requires matching derivative.)
// end_curvature = If `closed=false` gives curvature at last point. (Requires matching derivative.)
// curvature = List of curvature constraints for every point, or undef at unconstrained points. Each curvature constraint must be paired with a derivative constraint at the same point. Cannot be combined with `start_curvature=`/`end_curvature=`.
// corners = List of interior point indices where corners are permitted. Equivalent to setting entries of `deriv` to NAN.
// extra_pts = Number of extra control points to add to provide additional freedom to control undesirable oscillations. Default: 0
// smooth = Smoothness criterion used with extra control points. Set to 1 (minimize control-polygon length), 2 (minimize control-polygon bending) or 3 (minimize curve bending energy). Default: 3
//
// Example(2D): A NURBS curve where closed = false. (default)
// data = [[0,0], [10,30], [25,15], [40,35], [60,10], [80,25]];
// debug_nurbs_interp(data, 3);
//
// Example(2D,VPT=[40,12,0],VPD = 200): A NURBS curve where closed = true - Do NOT repeat the first point at the end.
// data = [[0,0], [30,50], [60,40], [80,10], [50,-20], [20,-10]];
// debug_nurbs_interp(data, 3, closed = true);
//
// Example(2D, VPT=[40,12,0],VPD = 200): Closed polygon - All data points lie exactly on the polygon boundary.
// data = [[0,0], [30,50], [60,40], [80,10], [50,-20], [20,-10]];
// path = nurbs_curve(nurbs_interp(data, 3, closed=true), splinesteps=16);
// polygon(path);
// color("red") move_copies(data) circle(r=1, $fn=16);
//
// Example(3D,VPT=[0,0,10],VPR=[60,0,330],VPD=120): 3D closed curve