-
Notifications
You must be signed in to change notification settings - Fork 10
/
ArachNURBS.py
5517 lines (4792 loc) · 206 KB
/
ArachNURBS.py
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
# ArachNURBS
# (c) Edward Mills 2016-2018
# edwardvmills@gmail.com
#
# ArachNURBS is a library of functions and classes to manipulate NURBS
# curves, surfaces, and the associated control polygons and grids.
# ArachNURBS is built on top FreeCAD's standard NURBS functions.
#
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
from __future__ import division # allows floating point division from integers
import Part
import FreeCAD
from FreeCAD import Base
from FreeCAD import Gui
import math
import numpy as np
# test message to verify load and reloads
print ("importing ArachNURBS")
#
## Basic NURBS rules and terms
## Order >= 2
## Order: 2 = line, 3 = quadratic, 4 = cubic ...
## Degree = Order - 1
## Order = Degree + 1
## nPoles >= Order
## nPoles >= Degree + 1
## nKnots = nPoles + Order
## nKnots = nPoles + degree + 1
## knot vector strictly ascending
## Pinned knot vector: k=Order, first k knots are equal, last k knots are equal
####
#### SECTION 1: DIRECT FUNCTIONS - NO PARAMETRIC LINKING BETWEEN OBJECTS
#### SECTION 2: PYTHON FEATURE CLASSES - PARAMETRIC LINKING BETWEEN OBJECTS (start around line 364)
####
### SECTION 1: DIRECT FUNCTIONS - NO PARAMETRIC LINKING BETWEEN OBJECTS
## Bottom up view:
## poles = 3D points with weights, as [[x,y,z],w], or [x,y,z] (these are leftovers waiting to receive weights).
## These are the basic input data for all that follows.
## They are obtained from the FreeCAD functions .getPoles() and .getWeights()
## NOTE: Poles in FreeCAD, such as returned by .getPoles(), refer only to xyz coordinates of a control point,
## THROUGHOUT the following functions, pole means [[x,y,z],w]
## lists are probably not efficient, but until FreeCAD has fully integrated homogeneous coordinates
## for all NURBS functions, this is easier for me :)
## Right now, the computation of these scripts is ridiculously fast compared
## to the time taken to generate the surfaces using the FreeCAD Part.BSplineSurface() function
## direct functions actually used in the Classes / available through the Silk FreeCAD workbench:
def equalVectors(vector0,vector1,tol): # 3D point equality test
if (vector1-vector0).Length <= tol:
return 1
elif (vector1-vector0).Length > tol:
return 0
def ClosestPointOnLine(a, b, p):
ap = p-a
ab = b-a
result = a + ap.dot(ab)/ab.dot(ab) * ab
return result
def int_2l(la,lb):
pa1=la.StartPoint
pa2=la.EndPoint
pb1=lb.StartPoint
pb2=lb.EndPoint
va=pa2-pa1
vb=pb2-pb1
lab=Part.Line(pa1,pb1)
ab=lab.length()
van=va.normalize()
vbn=vb.normalize()
pa3=pa1+van.multiply(10*ab)
pb3=pb1+vbn.multiply(10*ab)
lax=Part.Line(pa1,pa3)
lbx=Part.Line(pb1,pb3)
pln=Part.Plane(pa1,pb1,pa2)
int_0_1= lax.intersect2d(lbx,pln) # works down to 5.73 degrees between the lines
if int_0_1==[]:
return 'intersection failed'
int_abs_coord=pln.value(int_0_1[0][0],int_0_1[0][1])
return int_abs_coord
def lineOrPoint(p0,p1):
if equalVectors(p0, p1, .000001):
return [0]
else:
return [1, Part.LineSegment(p0, p1)]
def drawGrid(poles, columns):
nbPoles = len(poles)
# print ('nbPoles = ', nbPoles)
# print ('columns = ', columns)
rows = int(len(poles) / columns)
# print ('rows = ', rows)
#legs_total = 2 * rows * columns - rows - columns
#print ('legs_total = ', legs_total)
#legs = [0] * legs_total
legs = []
leg_index = 0
# loop all rows
for i in range(0,rows):
# loop over a row
for j in range(0, columns-1):
# print ('in row ', i)
start = i*columns + j
# print ('start = ', start)
end = i*columns + j + 1
# print ('end = ', end)
leg = lineOrPoint(poles[start], poles[end])
if (leg[0] == 1):
legs.append(leg[1])
# loop all columns
for i in range(0,columns):
# loop over a column
for j in range(0, rows-1):
# print ('in column ', i)
start = j*columns + i
# print ('start = ', start)
end = j*columns + columns + i
# print ('end = ', end)
leg = lineOrPoint(poles[start], poles[end])
if (leg[0] == 1):
legs.append(leg[1])
return legs
def orient_a_to_b(polesa,polesb): # polesa and polesb are lists of poles that share one endpoint.
# if needed, this function reorders a so that a.end = b.start or b.end. b is never modified
if equalVectors(polesa[-1],polesb[0],0.000001): # last point of first curve is first point of second curve
# curve 1 is oriented properly
return polesa
elif equalVectors(polesa[-1],polesb[-1],0.000001): # last point of first curve is last point of second curve
# curve 1 is oriented properly
return polesa
elif equalVectors(polesa[0],polesb[0],0.000001): # first point of first curve is first point of second curve
# curve 1 is reversed
return polesa[::-1]
elif equalVectors(polesa[0],polesb[-1],0.000001): # first point of first curve is last point of second curve
# curve 1 is reversed
return polesa[::-1]
else:
print ('curves do not share endpoints')
return 0
def Cubic_Bezier_ddu(pole0, pole1): # cubic derivative at curve start (pole1) based on first
# two poles (no curve required). Weights not included yet
P0=Base.Vector(pole0)
P1=Base.Vector(pole1)
Cubic_Bezier_ddu = (P1 - P0)*3
return Cubic_Bezier_ddu
def Cubic_6P_ddu(pole0, pole1): # cubic derivative at curve start (pole1) based on first
# two poles (no curve required). Weights not included yet.
P0=Base.Vector(pole0)
P1=Base.Vector(pole1)
Cubic_6P_ddu = (P1 - P0)*9
return Cubic_6P_ddu
def Cubic_Bezier_d2du2(pole0, pole1, pole2): # cubic second derivative at curve start (pole1) based on first
# three poles (no curve required). Weights not included yet.
P0=Base.Vector(pole0)
P1=Base.Vector(pole1)
P2=Base.Vector(pole2)
Cubic_Bezier_d2du2 = (P0- P1*2 + P2)*6
return Cubic_Bezier_d2du2
def Cubic_6P_d2du2(pole0, pole1, pole2): # cubic second derivative at curve start (pole1) based on first
# three poles (no curve required). Weights not included yet.
P0=Base.Vector(pole0)
P1=Base.Vector(pole1)
P2=Base.Vector(pole2)
Cubic_6P_d2du2 = (P0*2- P1*3 + P2)*27
return Cubic_6P_d2du2
def Cubic_Bezier_curvature(pole0, pole1, pole2): # curvature at curve start (pole1) based on the first three
# poles (no curve required). Weights not included yet.
ddu = Cubic_Bezier_ddu(pole0, pole1)
d2du2 = Cubic_Bezier_d2du2(pole0, pole1, pole2)
Cubic_Bezier_curvature = ddu.cross(d2du2).Length/ddu.Length.__pow__(3)
return Cubic_Bezier_curvature
def Cubic_6P_curvature(pole0, pole1, pole2): # curvature at curve start (pole1) based on the first three
# poles (no curve required). Weights not included yet
ddu = Cubic_6P_ddu(pole0, pole1)
d2du2 = Cubic_6P_d2du2(pole0, pole1, pole2)
Cubic_6P_curvature = ddu.cross(d2du2).Length/ddu.Length.__pow__(3)
return Cubic_6P_curvature
def Bezier_Cubic_curve(poles): # pinned cubic rational B spline, 4 control points
# Part.BSplineCurve(), cubic bezier form
#draws a degree 3 rational bspline from first to last point,
# second and third act as tangents
# poles is a list: [[[x,y,z],w],[[x,y,z],w],[[x,y,z],w],[[x,y,z],w]]
## nKnot = 4 + 3 +1 = 8
## Order = 3 + 1 = 4
degree=3
nPoles=4
knot=[0,0,0,0,1,1,1,1]
bs=Part.BSplineCurve()
bs.increaseDegree(degree)
id=1
for i in range(0,len(knot)): #-1):
bs.insertKnot(knot[i],id,0.0000001)
i=0
for ii in range(0,nPoles):
bs.setPole(ii+1,poles[i][0],poles[i][1])
i=i+1;
return bs
def Bezier_Bicubic_surf(grid_44): # given a 4 x 4 control grid, build the bicubic bezier
# surface from a Part.BSplineSurface() in Bicubic Bezier form
# len(knot_u) := nNodes_u + degree_u + 1
# len(knot_v) := nNodes_v + degree_v + 1
degree_u=3
degree_v=3
nNodes_u=4
nNodes_v=4
knot_u=[0,0,0,0,1,1,1,1]
knot_v=[0,0,0,0,1,1,1,1]
Bezier_Bicubic_surf=Part.BSplineSurface()
Bezier_Bicubic_surf.increaseDegree(degree_u,degree_v)
id=1
for i in range(0,len(knot_u)): #-1):
Bezier_Bicubic_surf.insertUKnot(knot_u[i],id,0.0000001)
id=1
for i in range(0,len(knot_v)): #-1):
Bezier_Bicubic_surf.insertVKnot(knot_v[i],id,0.0000001)
i=0
for jj in range(0,nNodes_v):
for ii in range(0,nNodes_u):
Bezier_Bicubic_surf.setPole(ii+1,jj+1,grid_44[i][0], grid_44[i][1]);
i=i+1;
return Bezier_Bicubic_surf
def NURBS_Cubic_6P_curve(poles): # pinned cubic rational Bspline, 6 control points
# Part.BSplineCurve(), just enough to have independent endpoint curvature
# draws a degree 3 rational bspline from first to last point,
# second and third act as tangents
# poles is a list: [[[x,y,z],w],[[x,y,z],w],[[x,y,z],w],[[x,y,z],w],[[x,y,z],w],[[x,y,z],w]]
## nKnot = 6 + 3 +1 = 10
## Order = 3 + 1 = 4
degree=3
nPoles=6
knot=[0,0,0,0,1.0/3.0,2.0/3.0,1,1,1,1]
bs=Part.BSplineCurve()
bs.increaseDegree(degree)
id=1
for i in range(0,len(knot)): #-1):
bs.insertKnot(knot[i],id,0.0000001)
i=0
for ii in range(0,nPoles):
bs.setPole(ii+1,poles[i][0],poles[i][1])
i=i+1;
return bs
def NURBS_Cubic_66_surf(grid_66): # given a 6 x 6 control grid, build the cubic
# NURBS surface from a Part.BSplineSurface().
# len(knot_u) := nNodes_u + degree_u + 1
# len(knot_v) := nNodes_v + degree_v + 1
degree_u=3
degree_v=3
nNodes_u=6
nNodes_v=6
knot_u=[0,0,0,0,1.0/3.0,2.0/3.0,1,1,1,1]
knot_v=[0,0,0,0,1.0/3.0,2.0/3.0,1,1,1,1]
NURBS_Cubic_66_surf=Part.BSplineSurface()
NURBS_Cubic_66_surf.increaseDegree(degree_u,degree_v)
id=1
for i in range(0,len(knot_u)): #-1):
NURBS_Cubic_66_surf.insertUKnot(knot_u[i],id,0.0000001)
id=1
for i in range(0,len(knot_v)): #-1):
NURBS_Cubic_66_surf.insertVKnot(knot_v[i],id,0.0000001)
i=0
for jj in range(0,nNodes_v):
for ii in range(0,nNodes_u):
NURBS_Cubic_66_surf.setPole(ii+1,jj+1,grid_66[i][0],grid_66[i][1]);
i=i+1;
return NURBS_Cubic_66_surf
def NURBS_Cubic_64_surf(grid_64): # given a 6 x 4 control grid, build the cubic
# NURBS surface from a Part.BSplineSurface().
# len(knot_u) := nNodes_u + degree_u + 1
# len(knot_v) := nNodes_v + degree_v + 1
degree_u=3
degree_v=3
nNodes_u=6
nNodes_v=4
knot_u=[0,0,0,0,1.0/3.0,2.0/3.0,1,1,1,1]
knot_v=[0,0,0,0,1,1,1,1]
NURBS_Cubic_64_surf=Part.BSplineSurface()
NURBS_Cubic_64_surf.increaseDegree(degree_u,degree_v)
id=1
for i in range(0,len(knot_u)): #-1):
NURBS_Cubic_64_surf.insertUKnot(knot_u[i],id,0.0000001)
id=1
for i in range(0,len(knot_v)): #-1):
NURBS_Cubic_64_surf.insertVKnot(knot_v[i],id,0.0000001)
i=0
for jj in range(0,nNodes_v):
for ii in range(0,nNodes_u):
NURBS_Cubic_64_surf.setPole(ii+1,jj+1,grid_64[i][0],grid_64[i][1]);
i=i+1;
return NURBS_Cubic_64_surf
def blend_poly_2x4_1x6(poles_0,weights_0, poles_1, weights_1, scale_0, scale_1, scale_2, scale_3):
# blend two cubic bezier into a 6 point cubic NURBS. this function assumes poles_0 flow into poles_1 without checking.
#print ("weights_0 in blend_poly_2x4_1x6")
#print (weights_0)
WeightedPoles_0=[[poles_0[0],weights_0[0]], [poles_0[1],weights_0[1]], [poles_0[2],weights_0[2]], [poles_0[3],weights_0[3]]]
CubicCurve4_0= Bezier_Cubic_curve(WeightedPoles_0)
CubicCurve6_0=CubicCurve4_0
CubicCurve6_0.insertKnot(1.0/3.0) # add knots to convert bezier to 6P
CubicCurve6_0.insertKnot(2.0/3.0)
#print ("weights_1 in blend_poly_2x4_1x6")
#print (weights_1)
WeightedPoles_1=[[poles_1[0],weights_1[0]], [poles_1[1],weights_1[1]], [poles_1[2],weights_1[2]], [poles_1[3],weights_1[3]]]
CubicCurve4_1= Bezier_Cubic_curve(WeightedPoles_1) # checked good BSplineSurface object.
CubicCurve6_1=CubicCurve4_1
CubicCurve6_1.insertKnot(1.0/3.0) # add knots to convert bezier to 6P
CubicCurve6_1.insertKnot(2.0/3.0)
poles_6_0=CubicCurve6_0.getPoles()
weights_6_0=CubicCurve6_0.getWeights()
# if the weights are too similar to each other, the knot insertion can convert them all to 1
# this is bad for blends with arc, ellipses, etc.
#print ("weights_6_0 in blend_poly_2x4_1x6")
#print (weights_6_0)
poles_6_1=CubicCurve6_1.getPoles()
weights_6_1=CubicCurve6_1.getWeights()
# if the weights are too similar to each other, the knot insertion can convert them all to 1
# this is bad for blends with arcs, ellipses, etc.
#print ("weights_6_1 in blend_poly_2x4_1x6")
#print (weights_6_1)
# check original weights....set == at 1%?
# this code is behaving very strangely.
# it failed to reset on a strict comparison (< .001), but resets correctly on a loose comparison.
# the numbers under comparison were equal to 8 or more decimals???
if (((weights_0[0]-weights_0[1] / weights_0[0]).__pow__(2) < .1) and
((weights_0[0]-weights_0[2] / weights_0[0]).__pow__(2) < .1) and
((weights_0[0]-weights_0[3] / weights_0[0]).__pow__(2) < .1)):
a = weights_0[0]
#print ("a ", a)
weights_6_0 = [a,a,a,a,a,a]
print("reseting weights_0")
if (((weights_1[0]-weights_1[1] / weights_1[0]).__pow__(2) < .1) and
((weights_1[0]-weights_1[2] / weights_1[0]).__pow__(2) < .1) and
((weights_1[0]-weights_1[3] / weights_1[0]).__pow__(2) < .1)):
b = weights_1[0]
#print ("b ", b)
weights_6_1 = [b,b,b,b,b,b]
print("reseting weights_1")
p0=[poles_6_0[0],weights_6_0[0]]
p1=[poles_6_0[1],weights_6_0[1]]
p2=[poles_6_0[2],weights_6_0[2]]
p3=[poles_6_1[3],weights_6_1[3]] ###
p4=[poles_6_1[4],weights_6_1[4]] ###
p5=[poles_6_1[5],weights_6_1[5]] ###
corner='p01p10'
### calculate curvature components
## start point
if (equalVectors(poles_6_0[0],poles_6_0[3], .000001) == 0):
# the first curve is non degenerate
l0 = p1[0]-p0[0] # first control leg
tan0=Base.Vector(l0) # make clean copy
tan0.normalize() # unit tangent direction
l1=Base.Vector(tan0) # make clean copy
l1.multiply(tan0.dot(p2[0]-p1[0])) # scalar projection of second control leg along unit tangent
h1=(p2[0]-p1[0])-l1 # height of second control leg orthogonal to tangent
### scale first control leg
L0=Base.Vector(l0) # make clean copy
L0.multiply(scale_0) # apply tangent scale
p1_scl = [p0[0] + L0, p1[1]] # reposition second control point
### calc new heights for first inner control leg
H1 = Base.Vector(h1) # make clean copy
H1.multiply(scale_0.__pow__(2)) # apply height scale
# apply inner scale
L1 = p1_scl[0] - p0[0] # rescale to new tangent (scale_0 already applied)
L1 = L1.multiply(scale_1) # apply inner tangent scale
p2_scl = [p1_scl[0] + H1 + L1, p2[1]] # reposition third control point
else:
# the first curve is degenerate
# just pass the inputs forward
p1_scl = [p1[0], p1[1]]
p2_scl = [p2[0], p2[1]]
## end point
if (equalVectors(poles_6_0[0],poles_6_0[3], .000001) == 0):
# the first curve is non degenerate
l4 = p4[0]-p5[0] # last control leg
tan4=Base.Vector(l4) # make clean copy
tan4.normalize() # unit tangent direction
l3=Base.Vector(tan4) # make clean copy
l3.multiply(tan4.dot(p3[0]-p4[0])) # scalar projection of second to last control leg along unit tangent
h3=(p3[0]-p4[0])-l3 # height of second control leg orthogonal to tangent
### scale last control leg
L4=Base.Vector(l4) # make clean copy
L4.multiply(scale_3) # apply tangent scale
p4_scl = [p5[0] + L4, p4[1]] # reposition fifth control point
### calc new heights for last inner control leg
H3 = Base.Vector(h3) # make clean copy
H3.multiply(scale_3.__pow__(2)) # apply height scale
# apply inner scale
L3 = p4_scl[0] - p5[0] # rescale to new tangent (scale_3 already applied)
L3 = L3.multiply(scale_2) # apply inner tangent scale
p3_scl = [p4_scl[0] + H3 + L3, p3[1]] # reposition third control point
else:
# the second curve is degenerate
# just pass the inputs forward
p3_scl = [p3[0], p3[1]]
p4_scl = [p4[0], p4[1]]
poles=[p0[0], p1_scl[0], p2_scl[0], p3_scl[0], p4_scl[0], p5[0]]
# set the weights. No scaling at this point. No idea what happens if one of the input curve is an arc.
# it would probably be a mess, since the curvature formulas above do not incorporate weights yet.
# actually, it seems to work just fine, with both circle and ellipse arcs? curvature calc unaffected by weights?
weights = [p0[1], p1[1], p2[1], p3[1], p4[1], p5[1]]
WeightedPoles= [[poles[0],weights[0]], [poles[1],weights[1]], [poles[2],weights[2]], [poles[3],weights[3]], [poles[4],weights[4]], [poles[5],weights[5]]]
current_test = NURBS_Cubic_6P_curve(WeightedPoles)
# we need to return the scales so the function result is compatible with the
#'Fair' and 'G3' version of the blend function, which modify these values
# not a strict requirement
return [poles,weights, scale_1, scale_2]
def Cubic_Bezier_dCds(pole0, pole1, pole2, pole3):
# calculate the rate of change of curvature per unit length (chord)
# at the beginning of a cubic bezier curve defined by the given poles
# calculate start point curvature directly from poles
C0 = Cubic_Bezier_curvature(pole0[0], pole1[0], pole2[0])
if math.fabs(C0) < 1.0e-6:
C0= 0.0
# prepare cubic bezier object to subdivide
Curve = Bezier_Cubic_curve([pole0, pole1, pole2, pole3])
# setup refinement loop
t_seg = 0.05 # initial segmentation value
segment_degen = 'false'
tol= 0.01
error = 1.0
loop_count = 0
dCds_last = 'not_ready'
while (error > tol and loop_count < 100 and segment_degen != 'true'):
Curve.segment(0,t_seg)
Poles = Curve.getPoles()
# check start curvature after segmentation
C0_seg = Cubic_Bezier_curvature(Poles[0], Poles[1], Poles[2])
if math.fabs(C0_seg) < 1.0e-6:
C0_seg= 0.0
# if the start curvature changes dramatically after segmentation,
# the new values are invalid. not a valid test when C0 = 0.0 to begin with
if C0 != 0.0:
if math.fabs((C0_seg - C0)/C0) > 5*tol:
segment_degen = 'true'
print ('segmentation has collapsed the curve')
print ('C0', C0, 'C0_check', C0_seg)
print ('Cubic_Bezier_dCds step ', loop_count)
elif C0 == 0.0:
if math.fabs((C0_seg - C0)) > .00001:
segment_degen = 'true'
print ('segmentation has collapsed the curve')
print ('C0', C0, 'C0_check', C0_seg)
print ('Cubic_Bezier_dCds step ', loop_count)
# calculate curvature at the end of the current segment
Cs = Cubic_Bezier_curvature(Poles[3], Poles[2], Poles[1])
#if the start curvature and first cut curvature are equal, then dCds is 0
if math.fabs(C0-Cs) <= 1.0e-6 and loop_count == 0:
return 0.0
# calculate chord length of current segment
S = Base.Vector(Poles[3])-Base.Vector(Poles[0])
s = S.Length
dCds_seg = (Cs-C0)/s
#print ('step ', loop_count, ' dCds_seg ', dCds_seg)
if loop_count > 1:
error = math.fabs((dCds_seg - dCds_last)/dCds_last)
if segment_degen != 'true':
dCds_last = dCds_seg
t_seg = t_seg * 0.9
loop_count=loop_count + 1
#print 'step ', loop_count, ' dCds_seg ', dCds_seg, ' error ', error
if error > tol:
#print 'no dCds found within ', tol, ' Cubic_Bezier_dCds'
dCds = dCds_last
#print 'returning dCds = ', dCds, ' within ', error, ' Cubic_Bezier_dCds'
else:
dCds = dCds_seg
return dCds
def Cubic_6P_dCds(pole0, pole1, pole2, pole3, pole4, pole5):
# calculate the rate of change of curvature per unit length (chord)
# at the beginning of a cubic 6P curve defined by the given poles
# calculate start point curvature directly from poles.
C0 = Cubic_6P_curvature(pole0[0], pole1[0], pole2[0])
if math.fabs(C0) < 1.0e-5:
C0= 0.0
# prepare cubic 6P object to segment
Curve = NURBS_Cubic_6P_curve([pole0, pole1, pole2, pole3, pole4, pole5])
# cut the 6P below the first internal knot
Curve.segment(0,.25)
poles = Curve.getPoles()
weights = Curve.getWeights()
# rebuild the weighted poles
WeightedPoles = [[poles[0],weights[0]], [poles[1],weights[1]], [poles[2],weights[2]], [poles[3],weights[3]]]
# pass the weighted poles down to the Bezier dCds function
dCds = Cubic_Bezier_dCds(WeightedPoles[0], WeightedPoles[1], WeightedPoles[2], WeightedPoles[3])
return dCds
def blendG3_poly_2x4_1x6(poles_0,weights_0, poles_1, weights_1, scale_0, scale_1, scale_2, scale_3): # work in progress. complete mess
# blend two cubic bezier into a 6 point cubic NURBS.
# this function assumes poles_0 flow into poles_1 without checking.
# rebuild both bezier inputs from the poles and weights
WeightedPoles_0=[[poles_0[0],weights_0[0]], [poles_0[1],weights_0[1]], [poles_0[2],weights_0[2]], [poles_0[3],weights_0[3]]]
CubicCurve4_0= Bezier_Cubic_curve(WeightedPoles_0)
WeightedPoles_1=[[poles_1[0],weights_1[0]], [poles_1[1],weights_1[1]], [poles_1[2],weights_1[2]], [poles_1[3],weights_1[3]]]
CubicCurve4_1= Bezier_Cubic_curve(WeightedPoles_1)
# set end point dC/ds targets
C0 = Cubic_6P_curvature(WeightedPoles_0[0][0], WeightedPoles_0[1][0], WeightedPoles_0[2][0])
C1 = Cubic_6P_curvature(WeightedPoles_1[3][0], WeightedPoles_1[2][0], WeightedPoles_1[1][0])
DC = math.fabs(C0-C1)
dCds0 = Cubic_Bezier_dCds(WeightedPoles_0[0], WeightedPoles_0[1], WeightedPoles_0[2], WeightedPoles_0[3])
dCds1 = Cubic_Bezier_dCds(WeightedPoles_1[3], WeightedPoles_1[2], WeightedPoles_1[1], WeightedPoles_1[0])
DdCds = math.fabs(dCds0-dCds1)
# quick, cheap, and very incomplete symmetry test.
if DC < 1.0e-8 and DdCds < 1.0e-8:
symmetric = 1
else:
symmetric = 0
#print "dCds inputs: " "dCds0, ", dCds0, " dCds1, ", dCds1
if math.fabs(dCds0) < 5.0e-6:
dCds0 = 0.0
if math.fabs(dCds1) < 5.0e-6:
dCds1 = 0.0
print ("dCds targets: " "dCds0, ", dCds0, " dCds1, ", dCds1," C0, ", C0, " C1, ", C1, "symmetric: ", symmetric)
# convert 4P inputs to 6P
CubicCurve6_0=CubicCurve4_0
CubicCurve6_0.insertKnot(1.0/3.0) # add knots to convert bezier to 6P
CubicCurve6_0.insertKnot(2.0/3.0)
CubicCurve6_1=CubicCurve4_1
CubicCurve6_1.insertKnot(1.0/3.0) # add knots to convert bezier to 6P
CubicCurve6_1.insertKnot(2.0/3.0)
# extract poles and weights from 6Ps
poles_6_0=CubicCurve6_0.getPoles()
weights_6_0=CubicCurve6_0.getWeights()
poles_6_1=CubicCurve6_1.getPoles()
weights_6_1=CubicCurve6_1.getWeights()
# need to come back here and make sure fractional weights were not reset to 1 (in the case where an
# entire row is equal). this can screw up rational grids, even though it doesn't matter for curves.
# check Cubic_6P_dCds
WeightedPoles_6_0=[[poles_6_0[0],weights_6_0[0]],
[poles_6_0[1],weights_6_0[1]],
[poles_6_0[2],weights_6_0[2]],
[poles_6_0[3],weights_6_0[3]],
[poles_6_0[4],weights_6_0[4]],
[poles_6_0[5],weights_6_0[5]]]
dCds6_0 = Cubic_6P_dCds(WeightedPoles_6_0[0],
WeightedPoles_6_0[1],
WeightedPoles_6_0[2],
WeightedPoles_6_0[3],
WeightedPoles_6_0[4],
WeightedPoles_6_0[5])
WeightedPoles_6_1=[[poles_6_1[0],weights_6_1[0]],
[poles_6_1[1],weights_6_1[1]],
[poles_6_1[2],weights_6_1[2]],
[poles_6_1[3],weights_6_1[3]],
[poles_6_1[4],weights_6_1[4]],
[poles_6_1[5],weights_6_1[5]]]
dCds6_1 = Cubic_6P_dCds(WeightedPoles_6_1[5],
WeightedPoles_6_1[4],
WeightedPoles_6_1[3],
WeightedPoles_6_1[2],
WeightedPoles_6_1[1],
WeightedPoles_6_1[0])
print ("dCds 6P check: " "dCds6_0, ", dCds6_0, " dCds6_1, ", dCds6_1)
# compile the blend poly. this initial form is G2, but clumped towards the outer points.
p0=[poles_6_0[0],weights_6_0[0]]
p1=[poles_6_0[1],weights_6_0[1]]
p2=[poles_6_0[2],weights_6_0[2]]
p3=[poles_6_1[3],weights_6_1[3]]
p4=[poles_6_1[4],weights_6_1[4]]
p5=[poles_6_1[5],weights_6_1[5]]
### calculate curvature components
## start point
l0 = p1[0]-p0[0] # first control leg
tan0=Base.Vector(l0) # make clean copy
tan0.normalize() # unit tangent direction
l1=Base.Vector(tan0) # make clean copy
l1.multiply(tan0.dot(p2[0]-p1[0])) # scalar projection of second control leg along unit tangent
h1=(p2[0]-p1[0])-l1 # height of second control leg orthogonal to tangent
## end point
l4 = p4[0]-p5[0] # last control leg
tan4=Base.Vector(l4) # make clean copy
tan4.normalize() # unit tangent direction
l3=Base.Vector(tan4) # make clean copy
l3.multiply(tan4.dot(p3[0]-p4[0])) # scalar projection of second to last control leg along unit tangent
h3=(p3[0]-p4[0])-l3 # height of second control leg orthogonal to tangent
### scale first and last control legs
L0=Base.Vector(l0) # make clean copy
L0.multiply(scale_0) # apply tangent scale
p1_scl = [p0[0] + L0, p1[1]] # reposition second control point
L4=Base.Vector(l4) # make clean copy
L4.multiply(scale_3) # apply tangent scale
p4_scl = [p5[0] + L4, p4[1]] # reposition fifth control point
### calc new heights for inner control legs
H1 = Base.Vector(h1) # make clean copy
H1.multiply(scale_0.__pow__(2)) # apply height scale
H3 = Base.Vector(h3) # make clean copy
H3.multiply(scale_3.__pow__(2)) # apply height scale
# search loop initial parameters
scale_1i = 1.0
scale_2i = 1.0
error = 1.0
step_size = 0.1
dir_0_prev = 0.0
dir_1_prev = 0.0
nudge_prev = 'none'
streak_0_count = 0
streak_1_count = 0
#step_stage_complete = 0
loop_count = 0
tol= 1.0e-5
while (error > tol and loop_count < 200 ):
# reset for next iteration
L1 = p1_scl[0] - p0[0] # rescale to new tangent (scale_0 already applied)
L3 = p4_scl[0] - p5[0] # rescale to new tangent (scale_3 already applied)
# apply scales
L1 = L1.multiply(scale_1i) # apply inner tangent scale
p2_scl = [p1_scl[0] + H1 + L1, p2[1]] # reposition third control point
L3 = L3.multiply(scale_2i) # apply inner tangent scale
p3_scl = [p4_scl[0] + H3 + L3, p3[1]] # reposition third control point
# prepare poles and weights function output
poles=[p0[0], p1_scl[0], p2_scl[0], p3_scl[0], p4_scl[0], p5[0]]
weights = [p0[1], p1[1], p2[1], p3[1], p4[1], p5[1]]
# prepare weighted poles for curvature analysis
WeightedPoles_6_i = [[poles[0],weights[0]],
[poles[1],weights[1]],
[poles[2],weights[2]],
[poles[3],weights[3]],
[poles[4],weights[4]],
[poles[5],weights[5]]]
# check both end curvatures
dCds6_0i = Cubic_6P_dCds(WeightedPoles_6_i[0],
WeightedPoles_6_i[1],
WeightedPoles_6_i[2],
WeightedPoles_6_i[3],
WeightedPoles_6_i[4],
WeightedPoles_6_i[5])
dCds6_1i = Cubic_6P_dCds(WeightedPoles_6_i[5],
WeightedPoles_6_i[4],
WeightedPoles_6_i[3],
WeightedPoles_6_i[2],
WeightedPoles_6_i[1],
WeightedPoles_6_i[0])
# define current G3 error
# proportional in non-zero cases,? absolute if target is 0
# the proportional error is a problem.
# it prioritizes smaller errors near zero than large errors proportionally closer to the target.
# in practice, this causes divergent run away situations
if dCds0 != 0.0:
error_0 = (dCds6_0i - dCds6_0) * ( 1 + 1 / math.fabs(dCds0)) / 2
elif dCds0 == 0.0:
error_0 = (dCds6_0i - dCds6_0) * ( 1 + 1 / math.fabs(dCds6_0i)) / 2
else:
print ("dCds0, ", dCds0, "returned from Cubic_6P_dCds()")
if dCds1 != 0.0:
error_1 = (dCds6_1i - dCds6_1) * ( 1 + 1 / math.fabs(dCds1)) / 2
elif dCds1 == 0.0:
error_1 = (dCds6_1i - dCds6_1) * ( 1 + 1 / math.fabs(dCds6_1i)) / 2 # < this has caused div by 0 errors
else:
print ("dCds1, ", dCds1, "returned from Cubic_6P_dCds()")
# success criteria for normal exit
if math.fabs(error_0) <= tol and math.fabs(error_1) <= tol:
print ("final ", loop_count, ": ","scl[", scale_1i, ", ", scale_2i, "] dCds[", dCds6_0i, ", ", dCds6_1i,"] err[", error_0, ", ", error_1,"]")
return [poles,weights,scale_1i,scale_2i]
error = math.fabs(error_0) + math.fabs(error_1)
# establish bounds for the errors based on initial guess
if loop_count == 0:
error_0_max = error_0
error_1_max = error_1
# if the loop takes BOTH errors beyond the initial bounds, undo last action, and reduce step size.
if loop_count != 0 and math.fabs(error_0) > math.fabs(error_0_max) and math.fabs(error_1) > math.fabs(error_1_max):
#undo last action
if nudge_prev == 0:
scale_1i = scale_1i - dir_prev * step_size
elif nudge_prev == 1:
scale_2i = scale_2i - dir_prev * step_size
# reduce step size
step_size = step_size / 2.0
#print "div - error grew"
# determine the required adjustment direction for each endpoint.
if error_0 > tol / 2.0:
direction_0 = 1
elif error_0 < -tol / 2.0:
direction_0 = -1
else:
direction_0 = 0
if error_1 > tol / 2.0:
direction_1 = 1
elif error_1 < -tol / 2.0:
direction_1 = -1
else:
direction_1 = 0
# plan the next action
if math.fabs(error_0) >= math.fabs(error_1):
nudge = 0
dir = direction_0
if symmetric != 1:
streak_0_count = streak_0_count + 1
streak_1_count = 0
elif math.fabs(error_0) < math.fabs(error_1):
nudge = 1
dir = direction_1
if symmetric != 1:
streak_1_count = streak_1_count + 1
streak_0_count = 0
if symmetric == 1 and direction_0 == direction_1:
nudge = 2
direction_2 = direction_0
elif symmetric == 1 and direction_0 != direction_1:
nudge = 0
direction = 0
symmetric = 0
print ("symmetry assumption has broken down.")
# compare the next planned action to the last executed action
if nudge == nudge_prev and dir != dir_prev:
# if we are undoing the previous action, reduce step size
step_size = step_size / 2.0
#print "div - direct undo"
# execute planned action, unless it is the same action 5 times in a row.
if (nudge == 0 and streak_0_count <= 5) or streak_1_count > 4:
scale_1i = scale_1i + direction_0 * step_size
dir_prev = direction_0
nudge_prev = 0
streak_1_count = 0
elif (nudge == 1 and streak_1_count <= 5) or streak_0_count > 4:
scale_2i = scale_2i + direction_1 * step_size
dir_prev = direction_1
nudge_prev = 1
streak_0_count = 0
elif nudge == 2:
scale_1i = scale_1i + direction_0 * step_size
scale_2i = scale_2i + direction_1 * step_size
dir_prev = direction_2
nudge_prev = 2
upper_limit = 3.0
lower_limit = 0.75
if scale_1i > upper_limit:
scale_1i = upper_limit
if scale_2i > upper_limit:
scale_2i = upper_limit
if scale_1i < lower_limit:
scale_1i = lower_limit
if scale_2i < lower_limit:
scale_2i = lower_limit
# G3 loop message
print (loop_count, ": ","scl[", scale_1i, ", ", scale_2i, "] dCds[", dCds6_0i, ", ", dCds6_1i,"] err[", error_0, ", ", error_1,"] dir[", direction_0, ", ", direction_1,"]act[",nudge_prev, ", ", dir_prev,"] streaks [", streak_0_count, ", ", streak_1_count, "]")
if (scale_1i == upper_limit and scale_2i == upper_limit) or (scale_1i == lower_limit and scale_2i == lower_limit):
break
if (scale_1i == upper_limit and scale_2i == lower_limit) or (scale_1i == lower_limit and scale_2i == upper_limit):
break
loop_count=loop_count + 1
# G3 final message
print ("final ", loop_count, ": ","scl[", scale_1i, ", ", scale_2i, "] dCds[", dCds6_0i, ", ", dCds6_1i,"] err[", error_0, ", ", error_1,"]")
# make sure the final values have been applied
L1 = p1_scl[0] - p0[0] # rescale to new tangent (scale_0 already applied)
L3 = p4_scl[0] - p5[0] # rescale to new tangent (scale_3 already applied)
# apply scales
L1 = L1.multiply(scale_1i) # apply inner tangent scale
p2_scl = [p1_scl[0] + H1 + L1, p2[1]] # reposition third control point
L3 = L3.multiply(scale_2i) # apply inner tangent scale
p3_scl = [p4_scl[0] + H3 + L3, p3[1]] # reposition third control point
# prepare poles and weights function output
poles=[p0[0], p1_scl[0], p2_scl[0], p3_scl[0], p4_scl[0], p5[0]]
weights = [p0[1], p1[1], p2[1], p3[1], p4[1], p5[1]]
return [poles,weights,scale_1i,scale_2i]
def blendFair_poly_2x4_1x6(poles_0,weights_0, poles_1, weights_1, scale_0, scale_1, scale_2, scale_3): # work in progress. complete mess
# blend two cubic bezier into a 6 point cubic NURBS.
# this function assumes poles_0 flow into poles_1 without checking.
# rebuild both bezier inputs from the poles and weights
WeightedPoles_0=[[poles_0[0],weights_0[0]], [poles_0[1],weights_0[1]], [poles_0[2],weights_0[2]], [poles_0[3],weights_0[3]]]
CubicCurve4_0= Bezier_Cubic_curve(WeightedPoles_0)
WeightedPoles_1=[[poles_1[0],weights_1[0]], [poles_1[1],weights_1[1]], [poles_1[2],weights_1[2]], [poles_1[3],weights_1[3]]]
CubicCurve4_1= Bezier_Cubic_curve(WeightedPoles_1)
# set end point dC/ds targets
C0 = Cubic_6P_curvature(WeightedPoles_0[0][0], WeightedPoles_0[1][0], WeightedPoles_0[2][0])
C1 = Cubic_6P_curvature(WeightedPoles_1[3][0], WeightedPoles_1[2][0], WeightedPoles_1[1][0])
DC = math.fabs(C0-C1)
dCds0 = Cubic_Bezier_dCds(WeightedPoles_0[0], WeightedPoles_0[1], WeightedPoles_0[2], WeightedPoles_0[3])
dCds1 = Cubic_Bezier_dCds(WeightedPoles_1[3], WeightedPoles_1[2], WeightedPoles_1[1], WeightedPoles_1[0])
DdCds = math.fabs(dCds0-dCds1)
# quick, cheap, and very incomplete symmetry test.
if DC < 1.0e-8 and DdCds < 1.0e-8:
symmetric = 1
else:
symmetric = 0
#print "dCds inputs: " "dCds0, ", dCds0, " dCds1, ", dCds1
if math.fabs(dCds0) < 5.0e-6:
dCds0 = 0.0
if math.fabs(dCds1) < 5.0e-6:
dCds1 = 0.0
print ("dCds targets: " "dCds0, ", dCds0, " dCds1, ", dCds1," C0, ", C0, " C1, ", C1, "symmetric: ", symmetric)
# convert 4P inputs to 6P
CubicCurve6_0=CubicCurve4_0
CubicCurve6_0.insertKnot(1.0/3.0) # add knots to convert bezier to 6P
CubicCurve6_0.insertKnot(2.0/3.0)
CubicCurve6_1=CubicCurve4_1
CubicCurve6_1.insertKnot(1.0/3.0) # add knots to convert bezier to 6P
CubicCurve6_1.insertKnot(2.0/3.0)
# extract poles and weights from 6Ps
poles_6_0=CubicCurve6_0.getPoles()
weights_6_0=CubicCurve6_0.getWeights()
poles_6_1=CubicCurve6_1.getPoles()
weights_6_1=CubicCurve6_1.getWeights()
# check Cubic_6P_dCds
WeightedPoles_6_0=[[poles_6_0[0],weights_6_0[0]],
[poles_6_0[1],weights_6_0[1]],
[poles_6_0[2],weights_6_0[2]],
[poles_6_0[3],weights_6_0[3]],
[poles_6_0[4],weights_6_0[4]],
[poles_6_0[5],weights_6_0[5]]]
dCds6_0 = Cubic_6P_dCds(WeightedPoles_6_0[0],
WeightedPoles_6_0[1],
WeightedPoles_6_0[2],
WeightedPoles_6_0[3],
WeightedPoles_6_0[4],
WeightedPoles_6_0[5])
WeightedPoles_6_1=[[poles_6_1[0],weights_6_1[0]],
[poles_6_1[1],weights_6_1[1]],
[poles_6_1[2],weights_6_1[2]],
[poles_6_1[3],weights_6_1[3]],
[poles_6_1[4],weights_6_1[4]],
[poles_6_1[5],weights_6_1[5]]]
dCds6_1 = Cubic_6P_dCds(WeightedPoles_6_1[5],
WeightedPoles_6_1[4],
WeightedPoles_6_1[3],
WeightedPoles_6_1[2],
WeightedPoles_6_1[1],
WeightedPoles_6_1[0])
print ("dCds 6P check: " "dCds6_0, ", dCds6_0, " dCds6_1, ", dCds6_1)
# compile the blend poly. this initial form is G2, but clumped towards the outer points.
p0=[poles_6_0[0],weights_6_0[0]]
p1=[poles_6_0[1],weights_6_0[1]]
p2=[poles_6_0[2],weights_6_0[2]]
p3=[poles_6_1[3],weights_6_1[3]]
p4=[poles_6_1[4],weights_6_1[4]]
p5=[poles_6_1[5],weights_6_1[5]]
### calculate curvature components
## start point
l0 = p1[0]-p0[0] # first control leg
tan0=Base.Vector(l0) # make clean copy
tan0.normalize() # unit tangent direction
l1=Base.Vector(tan0) # make clean copy
l1.multiply(tan0.dot(p2[0]-p1[0])) # scalar projection of second control leg along unit tangent
h1=(p2[0]-p1[0])-l1 # height of second control leg orthogonal to tangent
## end point
l4 = p4[0]-p5[0] # last control leg
tan4=Base.Vector(l4) # make clean copy
tan4.normalize() # unit tangent direction
l3=Base.Vector(tan4) # make clean copy
l3.multiply(tan4.dot(p3[0]-p4[0])) # scalar projection of second to last control leg along unit tangent
h3=(p3[0]-p4[0])-l3 # height of second control leg orthogonal to tangent
### scale first and last control legs
L0=Base.Vector(l0) # make clean copy
L0.multiply(scale_0) # apply tangent scale
p1_scl = [p0[0] + L0, p1[1]] # reposition second control point
L4=Base.Vector(l4) # make clean copy
L4.multiply(scale_3) # apply tangent scale
p4_scl = [p5[0] + L4, p4[1]] # reposition fifth control point
### calc new heights for inner control legs
H1 = Base.Vector(h1) # make clean copy
H1.multiply(scale_0.__pow__(2)) # apply height scale
H3 = Base.Vector(h3) # make clean copy
H3.multiply(scale_3.__pow__(2)) # apply height scale
# output of the G2 portion
p0[0] # start point. unchanged
p1_scl # start tangent point. scaled
H1 # start inner tangent height. scaled
H3 # end inner tangent height. scaled
p4_scl[0] # end tangent point. scaled
p5[0] # end point. unchanged
# stuff involved in making a single fairing attempt and evaluating dCds to establish G3 error.
p0[0] # start point. unchanged
p1_scl[0] # start tangent point. unchanged
scale_1i # start inner tangent scale : independent variable in the attempt
H1 # start inner tangent height. unchanged
p2[1] # start inner tangent weight. unchanged
p3[1] # end inner tangent weight.
scale_2i # end inner tangent scale : independent variable in the attempt
H3 # end inner tangent height. unchanged
p4_scl[0] # end tangent point. unchanged
p5[0] # end point. unchanged
dCds0 # input poly start curvature derivative with respect to arclength (along +t). dependent variable in the attempt
dCds1 # input poly end curvature derivative with respect to arclength (along -t). dependent variable in the attempt