projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0025763: SSP sample is unable to build the proper surface on the given cloud of points
[occt.git] / src / BSplCLib / BSplCLib_2.cxx
CommitLineData
b311480e 1// Created on: 1995-01-17
2// Created by: Laurent BOURESCHE
3// Copyright (c) 1995-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 8// This library is free software; you can redistribute it and/or modify it under
9// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1 10// by the Free Software Foundation, with special exception defined in the file
11// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12// distribution for complete text of the license and disclaimer of any warranty.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
b311480e 16
7fd59977 17// xab : modified 15-Mar-95 : added BuildBSplMatrix,FactorBandedMatrix,
18// EvalBsplineBasis,
19// EvalPolynomial : Horners method
20
7fd59977 21#include <Standard_Stream.hxx>
22
23#include <BSplCLib.hxx>
24#include <gp_Mat2d.hxx>
25#include <PLib.hxx>
26#include <Standard_ConstructionError.hxx>
27#include <TColStd_Array1OfReal.hxx>
28#include <TColStd_Array1OfInteger.hxx>
29#include <TColStd_HArray1OfReal.hxx>
30#include <TColStd_HArray1OfInteger.hxx>
31
32#include <math_Matrix.hxx>
33
34//=======================================================================
35//struct : BSplCLib_DataContainer
36//purpose: Auxiliary structure providing buffers for poles and knots used in
37// evaluation of bspline (allocated in the stack)
38//=======================================================================
39
40struct BSplCLib_DataContainer
41{
41194117 42 BSplCLib_DataContainer(Standard_Integer Degree)
7fd59977 43 {
105aae76 44 (void)Degree; // avoid compiler warning
41194117
K		 45 Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree() ||
46 BSplCLib::MaxDegree() > 25,
47 "BSplCLib: bspline degree is greater than maximum supported");
7fd59977 48 }
49
50 Standard_Real poles[2*(25+1)];
51 Standard_Real knots[2*25];
52 Standard_Real ders[4];
53};
54
55// methods for 1 dimensional BSplines
56
57//=======================================================================
58//function : BuildEval
59//purpose : builds the local array for evaluation
60//=======================================================================
61
62void BSplCLib::BuildEval(const Standard_Integer Degree,
63 const Standard_Integer Index,
64 const TColStd_Array1OfReal& Poles,
65 const TColStd_Array1OfReal& Weights,
66 Standard_Real& LP)
67{
68 Standard_Integer PLower = Poles.Lower();
69 Standard_Integer PUpper = Poles.Upper();
70 Standard_Integer i;
71 Standard_Integer ip = PLower + Index - 1;
72 Standard_Real w, *pole = &LP;
73 if (&Weights == NULL) {
74
75 for (i = 0; i <= Degree; i++) {
76 ip++;
77 if (ip > PUpper) ip = PLower;
78 pole[0] = Poles(ip);
79 pole += 1;
80 }
81 }
82 else {
83
84 for (i = 0; i <= Degree; i++) {
85 ip++;
86 if (ip > PUpper) ip = PLower;
87 pole[1] = w = Weights(ip);
88 pole[0] = Poles(ip) * w;
89 pole += 2;
90 }
91 }
92}
93
94//=======================================================================
95//function : PrepareEval
96//purpose : stores data for Eval in the local arrays
97// dc.poles and dc.knots
98//=======================================================================
99
100static void PrepareEval
101(Standard_Real& u,
102 Standard_Integer& index,
103 Standard_Integer& dim,
104 Standard_Boolean& rational,
105 const Standard_Integer Degree,
106 const Standard_Boolean Periodic,
107 const TColStd_Array1OfReal& Poles,
108 const TColStd_Array1OfReal& Weights,
109 const TColStd_Array1OfReal& Knots,
110 const TColStd_Array1OfInteger& Mults,
111 BSplCLib_DataContainer& dc)
112{
113 // Set the Index
114 BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
115
116 // make the knots
117 BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*dc.knots);
118 if (&Mults == NULL)
119 index -= Knots.Lower() + Degree;
120 else
121 index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults);
122
123 // check truly rational
124 rational = (&Weights != NULL);
125 if (rational) {
126 Standard_Integer WLower = Weights.Lower() + index;
127 rational = BSplCLib::IsRational(Weights, WLower, WLower + Degree);
128 }
129
130 // make the poles
131 if(rational) {
132 dim = 2;
133 BSplCLib::BuildEval(Degree, index, Poles, Weights , *dc.poles);
134 }
135 else {
136 dim = 1;
137 BSplCLib::BuildEval(Degree, index, Poles, BSplCLib::NoWeights(), *dc.poles);
138 }
139}
140
141//=======================================================================
142//function : D0
143//purpose :
144//=======================================================================
145
146void BSplCLib::D0
147(const Standard_Real U,
148 const Standard_Integer Index,
149 const Standard_Integer Degree,
150 const Standard_Boolean Periodic,
151 const TColStd_Array1OfReal& Poles,
152 const TColStd_Array1OfReal& Weights,
153 const TColStd_Array1OfReal& Knots,
154 const TColStd_Array1OfInteger& Mults,
155 Standard_Real& P)
156{
157 Standard_Integer dim,index = Index;
158 Standard_Real u = U;
159 Standard_Boolean rational;
160 BSplCLib_DataContainer dc(Degree);
161 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
162 BSplCLib::Eval(u,Degree,*dc.knots,dim,*dc.poles);
163 if (rational) P = dc.poles[0] / dc.poles[1];
164 else P = dc.poles[0];
165}
166
167//=======================================================================
168//function : D1
169//purpose :
170//=======================================================================
171
172void BSplCLib::D1
173(const Standard_Real U,
174 const Standard_Integer Index,
175 const Standard_Integer Degree,
176 const Standard_Boolean Periodic,
177 const TColStd_Array1OfReal& Poles,
178 const TColStd_Array1OfReal& Weights,
179 const TColStd_Array1OfReal& Knots,
180 const TColStd_Array1OfInteger& Mults,
181 Standard_Real& P,
182 Standard_Real& V)
183{
184 Standard_Integer dim,index = Index;
185 Standard_Real u = U;
186 Standard_Boolean rational;
187 BSplCLib_DataContainer dc(Degree);
188 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
189 BSplCLib::Bohm(u,Degree,1,*dc.knots,dim,*dc.poles);
190 Standard_Real *result = dc.poles;
191 if (rational) {
192 PLib::RationalDerivative(Degree,1,1,*dc.poles,*dc.ders);
193 result = dc.ders;
194 }
195 P = result[0];
196 V = result[1];
197}
198
199//=======================================================================
200//function : D2
201//purpose :
202//=======================================================================
203
204void BSplCLib::D2
205(const Standard_Real U,
206 const Standard_Integer Index,
207 const Standard_Integer Degree,
208 const Standard_Boolean Periodic,
209 const TColStd_Array1OfReal& Poles,
210 const TColStd_Array1OfReal& Weights,
211 const TColStd_Array1OfReal& Knots,
212 const TColStd_Array1OfInteger& Mults,
213 Standard_Real& P,
214 Standard_Real& V1,
215 Standard_Real& V2)
216{
217 Standard_Integer dim,index = Index;
218 Standard_Real u = U;
219 Standard_Boolean rational;
220 BSplCLib_DataContainer dc(Degree);
221 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
222 BSplCLib::Bohm(u,Degree,2,*dc.knots,dim,*dc.poles);
223 Standard_Real *result = dc.poles;
224 if (rational) {
225 PLib::RationalDerivative(Degree,2,1,*dc.poles,*dc.ders);
226 result = dc.ders;
227 }
228 P = result[0];
229 V1 = result[1];
230 if (!rational && (Degree < 2)) V2 = 0.;
231 else V2 = result[2];
232}
233
234//=======================================================================
235//function : D3
236//purpose :
237//=======================================================================
238
239void BSplCLib::D3
240(const Standard_Real U,
241 const Standard_Integer Index,
242 const Standard_Integer Degree,
243 const Standard_Boolean Periodic,
244 const TColStd_Array1OfReal& Poles,
245 const TColStd_Array1OfReal& Weights,
246 const TColStd_Array1OfReal& Knots,
247 const TColStd_Array1OfInteger& Mults,
248 Standard_Real& P,
249 Standard_Real& V1,
250 Standard_Real& V2,
251 Standard_Real& V3)
252{
253 Standard_Integer dim,index = Index;
254 Standard_Real u = U;
255 Standard_Boolean rational;
256 BSplCLib_DataContainer dc(Degree);
257 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
258 BSplCLib::Bohm(u,Degree,3,*dc.knots,dim,*dc.poles);
259 Standard_Real *result = dc.poles;
260 if (rational) {
261 PLib::RationalDerivative(Degree,3,1,*dc.poles,*dc.ders);
262 result = dc.ders;
263 }
264 P = result[0];
265 V1 = result[1];
266 if (!rational && (Degree < 2)) V2 = 0.;
267 else V2 = result[2];
268 if (!rational && (Degree < 3)) V3 = 0.;
269 else V3 = result[3];
270}
271
272//=======================================================================
273//function : DN
274//purpose :
275//=======================================================================
276
277void BSplCLib::DN
278(const Standard_Real U,
279 const Standard_Integer N,
280 const Standard_Integer Index,
281 const Standard_Integer Degree,
282 const Standard_Boolean Periodic,
283 const TColStd_Array1OfReal& Poles,
284 const TColStd_Array1OfReal& Weights,
285 const TColStd_Array1OfReal& Knots,
286 const TColStd_Array1OfInteger& Mults,
287 Standard_Real& VN)
288{
289 Standard_Integer dim,index = Index;
290 Standard_Real u = U;
291 Standard_Boolean rational;
292 BSplCLib_DataContainer dc(Degree);
293 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
294 BSplCLib::Bohm(u,Degree,N,*dc.knots,dim,*dc.poles);
295 if (rational) {
296 Standard_Real v;
297 PLib::RationalDerivative(Degree,N,1,*dc.poles,v,Standard_False);
298 VN = v;
299 }
300 else {
301 if (N > Degree) VN = 0.;
302 else VN = dc.poles[N];
303 }
304}
305
306//=======================================================================
307//function : Build BSpline Matrix
308//purpose : Builds the Bspline Matrix
309//=======================================================================
310
311Standard_Integer
312BSplCLib::BuildBSpMatrix(const TColStd_Array1OfReal& Parameters,
313 const TColStd_Array1OfInteger& ContactOrderArray,
314 const TColStd_Array1OfReal& FlatKnots,
315 const Standard_Integer Degree,
316 math_Matrix& Matrix,
317 Standard_Integer& UpperBandWidth,
318 Standard_Integer& LowerBandWidth)
319{
320 Standard_Integer ii,
321 jj,
322 Index,
323 ErrorCode,
324 ReturnCode = 0,
325 FirstNonZeroBsplineIndex,
326 BandWidth,
327 MaxOrder = 21,
328 Order ;
329
330 math_Matrix BSplineBasis(1, MaxOrder,
331 1, MaxOrder) ;
332
333 Order = Degree + 1 ;
334 UpperBandWidth = Degree ;
335 LowerBandWidth = Degree ;
336 BandWidth = UpperBandWidth + LowerBandWidth + 1 ;
337 if (Matrix.LowerRow() != Parameters.Lower() ||
338 Matrix.UpperRow() != Parameters.Upper() ||
339 Matrix.LowerCol() != 1 ||
340 Matrix.UpperCol() != BandWidth) {
341 ReturnCode = 1;
342 goto FINISH ;
343 }
344
345 for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
346 ErrorCode =
347 BSplCLib::EvalBsplineBasis(1,
348 ContactOrderArray(ii),
349 Order,
350 FlatKnots,
351 Parameters(ii),
352
353 FirstNonZeroBsplineIndex,
354 BSplineBasis) ;
355 if (ErrorCode != 0) {
356 ReturnCode = 2 ;
357 goto FINISH ;
358 }
359 Index = LowerBandWidth + 1 + FirstNonZeroBsplineIndex - ii ;
360
361 for (jj = 1 ; jj < Index ; jj++) {
362 Matrix.Value(ii,jj) = 0.0e0 ;
363 }
364
365 for (jj = 1 ; jj <= Order ; jj++) {
366 Matrix.Value(ii,Index) = BSplineBasis(ContactOrderArray(ii) + 1, jj) ;
367 Index += 1 ;
368 }
369
370 for (jj = Index ; jj <= BandWidth ; jj++) {
371 Matrix.Value(ii,jj) = 0.0e0 ;
372 }
373 }
374 FINISH : ;
375 return (ReturnCode) ;
376}
377
378//=======================================================================
379//function : Makes LU decompositiomn without Pivoting
380//purpose : Builds the Bspline Matrix
381//=======================================================================
382
383Standard_Integer
384BSplCLib::FactorBandedMatrix(math_Matrix& Matrix,
385 const Standard_Integer UpperBandWidth,
386 const Standard_Integer LowerBandWidth,
387 Standard_Integer& PivotIndexProblem)
388{
389 Standard_Integer ii,
390 jj,
391 kk,
392 Index,
393 MinIndex,
394 MaxIndex,
395 ReturnCode = 0,
396 BandWidth = UpperBandWidth + LowerBandWidth + 1 ;
397
398 Standard_Real Inverse ;
399 PivotIndexProblem = 0 ;
400
401 for (ii = Matrix.LowerRow() + 1 ; ii <= Matrix.UpperRow() ; ii++) {
402 MinIndex = ( LowerBandWidth - ii + 2 >= 1 ? LowerBandWidth - ii + 2 : 1) ;
403
404 for (jj = MinIndex ; jj <= LowerBandWidth ; jj++) {
405 Index = ii - LowerBandWidth + jj - 1 ;
406 Inverse = Matrix(Index ,LowerBandWidth + 1) ;
407 if (Abs(Inverse) > RealSmall()) {
408 Inverse = -1.0e0/Inverse ;
409 }
410 else {
411 ReturnCode = 1 ;
412 PivotIndexProblem = Index ;
413 goto FINISH ;
414 }
415 Matrix(ii,jj) = Matrix(ii,jj) * Inverse ;
416 MaxIndex = BandWidth + Index - ii ;
417
418 for (kk = jj + 1 ; kk <= MaxIndex ; kk++) {
419 Matrix(ii,kk) += Matrix(ii,jj) * Matrix(Index, kk + ii - Index) ;
420 }
421 }
422 }
423 FINISH :
424 return (ReturnCode) ;
425}
426
427//=======================================================================
428//function : Build BSpline Matrix
429//purpose : Builds the Bspline Matrix
430//=======================================================================
431
432Standard_Integer
433BSplCLib::EvalBsplineBasis
434//(const Standard_Integer Side, // = 1 rigth side, -1 left side
435(const Standard_Integer , // = 1 rigth side, -1 left side
436 const Standard_Integer DerivativeRequest,
437 const Standard_Integer Order,
438 const TColStd_Array1OfReal& FlatKnots,
439 const Standard_Real Parameter,
440 Standard_Integer& FirstNonZeroBsplineIndex,
ecbdb1b0 441 math_Matrix& BsplineBasis,
442 Standard_Boolean isPeriodic)
7fd59977 443{
444 // the matrix must have at least DerivativeRequest + 1
445 // row and Order columns
446 // the result are stored in the following way in
447 // the Bspline matrix
448 // Let i be the FirstNonZeroBsplineIndex and
449 // t be the parameter value, k the order of the
450 // knot vector, r the DerivativeRequest :
451 //
452 // B (t) B (t) B (t)
453 // i i+1 i+k-1
454 //
455 // (1) (1) (1)
456 // B (t) B (t) B (t)
457 // i i+1 i+k-1
458 //
459 //
460 //
461 //
462 // (r) (r) (r)
463 // B (t) B (t) B (t)
464 // i i+1 i+k-1
465 //
466 Standard_Integer
467 ReturnCode,
468 ii,
469 pp,
470 qq,
471 ss,
472 NumPoles,
473 LocalRequest ;
474// ,Index ;
475
476 Standard_Real NewParameter,
477 Inverse,
478 Factor,
479 LocalInverse,
480 Saved ;
481// , *FlatKnotsArray ;
482
483 ReturnCode = 0 ;
484 FirstNonZeroBsplineIndex = 0 ;
485 LocalRequest = DerivativeRequest ;
486 if (DerivativeRequest >= Order) {
487 LocalRequest = Order - 1 ;
488 }
489
490 if (BsplineBasis.LowerCol() != 1 ||
491 BsplineBasis.UpperCol() < Order ||
492 BsplineBasis.LowerRow() != 1 ||
493 BsplineBasis.UpperRow() <= LocalRequest) {
494 ReturnCode = 1;
495 goto FINISH ;
496 }
497 NumPoles = FlatKnots.Upper() - FlatKnots.Lower() + 1 - Order ;
498 BSplCLib::LocateParameter(Order - 1,
499 FlatKnots,
500 Parameter,
ecbdb1b0 501 isPeriodic,
7fd59977 502 Order,
503 NumPoles+1,
504 ii,
505 NewParameter) ;
506
507 FirstNonZeroBsplineIndex = ii - Order + 1 ;
508
509 BsplineBasis(1,1) = 1.0e0 ;
510 LocalRequest = DerivativeRequest ;
511 if (DerivativeRequest >= Order) {
512 LocalRequest = Order - 1 ;
513 }
514
515 for (qq = 2 ; qq <= Order - LocalRequest ; qq++) {
516 BsplineBasis(1,qq) = 0.0e0 ;
517
518 for (pp = 1 ; pp <= qq - 1 ; pp++) {
519 //
520 // this should be always invertible if ii is correctly computed
521 //
522 Factor = (Parameter - FlatKnots(ii - qq + pp + 1))
523 / (FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1)) ;
524 Saved = Factor * BsplineBasis(1,pp) ;
525 BsplineBasis(1,pp) *= (1.0e0 - Factor) ;
526 BsplineBasis(1,pp) += BsplineBasis(1,qq) ;
527 BsplineBasis(1,qq) = Saved ;
528 }
529 }
530
531 for (qq = Order - LocalRequest + 1 ; qq <= Order ; qq++) {
532
533 for (pp = 1 ; pp <= qq - 1 ; pp++) {
534 BsplineBasis(Order - qq + 2,pp) = BsplineBasis(1,pp) ;
535 }
536 BsplineBasis(1,qq) = 0.0e0 ;
537
538 for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) {
539 BsplineBasis(Order - ss + 2,qq) = 0.0e0 ;
540 }
541
542 for (pp = 1 ; pp <= qq - 1 ; pp++) {
543 Inverse = 1.0e0 / (FlatKnots(ii + pp) - FlatKnots(ii - qq + pp + 1)) ;
544 Factor = (Parameter - FlatKnots(ii - qq + pp + 1)) * Inverse ;
545 Saved = Factor * BsplineBasis(1,pp) ;
546 BsplineBasis(1,pp) *= (1.0e0 - Factor) ;
547 BsplineBasis(1,pp) += BsplineBasis(1,qq) ;
548 BsplineBasis(1,qq) = Saved ;
549 LocalInverse = (Standard_Real) (qq - 1) * Inverse ;
550
551 for (ss = Order - LocalRequest + 1 ; ss <= qq ; ss++) {
552 Saved = LocalInverse * BsplineBasis(Order - ss + 2, pp) ;
553 BsplineBasis(Order - ss + 2, pp) *= - LocalInverse ;
554 BsplineBasis(Order - ss + 2, pp) += BsplineBasis(Order - ss + 2,qq) ;
555 BsplineBasis(Order - ss + 2,qq) = Saved ;
556 }
557 }
558 }
559 FINISH :
560 return (ReturnCode) ;
561}
562
563//=======================================================================
564//function : MovePointAndTangent
565//purpose :
566//=======================================================================
567
568void BSplCLib::MovePointAndTangent(const Standard_Real U,
569 const Standard_Integer ArrayDimension,
570 Standard_Real &Delta,
571 Standard_Real &DeltaDerivatives,
572 const Standard_Real Tolerance,
573 const Standard_Integer Degree,
574 const Standard_Boolean Rational,
575 const Standard_Integer StartingCondition,
576 const Standard_Integer EndingCondition,
577 Standard_Real& Poles,
578 const TColStd_Array1OfReal& Weights,
579 const TColStd_Array1OfReal& FlatKnots,
580 Standard_Real& NewPoles,
581 Standard_Integer& ErrorStatus)
582{
583 Standard_Integer num_poles,
584 num_knots,
585 ii,
586 jj,
587 conditions,
588 start_num_poles,
589 end_num_poles,
590 index,
591 start_index,
592 end_index,
593 other_index,
594 type,
595 order ;
596
597 Standard_Real new_parameter,
598 value,
599 divide,
600 end_value,
601 start_value,
602 *poles_array,
603 *new_poles_array,
604 *delta_array,
605 *derivatives_array,
606 *weights_array ;
607
608 ErrorStatus = 0 ;
609 weights_array = NULL ;
610 if (Rational) {
611 weights_array = (Standard_Real *) &Weights(Weights.Lower()) ;
612 }
613
614 poles_array = &Poles ;
615 new_poles_array = &NewPoles ;
616 delta_array = &Delta ;
617 derivatives_array = &DeltaDerivatives ;
618 order = Degree + 1 ;
619 num_knots = FlatKnots.Length() ;
620 num_poles = num_knots - order ;
621 conditions = StartingCondition + EndingCondition + 4 ;
622 //
623 // check validity of input data
624 //
625 if (StartingCondition >= -1 &&
626 StartingCondition <= Degree &&
627 EndingCondition >= -1 &&
628 EndingCondition <= Degree &&
629 conditions <= num_poles) {
630 //
631 // check the parameter is within bounds
632 //
633 start_index = FlatKnots.Lower() + Degree ;
634 end_index = FlatKnots.Upper() - Degree ;
635 //
636 // check if there is enough room to move the poles
637 //
638 conditions = 1 ;
639 if (StartingCondition == -1) {
640 conditions = conditions && (FlatKnots(start_index) <= U) ;
641 }
642 else {
643 conditions = conditions && (FlatKnots(start_index) + Tolerance < U) ;
644 }
645 if (EndingCondition == -1) {
646 conditions = conditions && (FlatKnots(end_index) >= U) ;
647 }
648 else {
649 conditions = conditions && (FlatKnots(end_index) - Tolerance > U) ;
650 }
651
652 if (conditions) {
653 //
654 // build 2 auxialiary functions
655 //
656 TColStd_Array1OfReal schoenberg_points(1,num_poles) ;
657 TColStd_Array1OfReal first_function (1,num_poles) ;
658 TColStd_Array1OfReal second_function (1,num_poles) ;
659
660 BuildSchoenbergPoints(Degree,
661 FlatKnots,
662 schoenberg_points) ;
663 start_index = StartingCondition + 2 ;
664 end_index = num_poles - EndingCondition - 1 ;
665 LocateParameter(schoenberg_points,
666 U,
667 Standard_False,
668 start_index,
669 end_index,
670 index,
671 new_parameter,
672 0, 1) ;
673
674 if (index == start_index) {
675 other_index = index + 1 ;
676 }
677 else if (index == end_index) {
678 other_index = index -1 ;
679 }
680 else if (U - FlatKnots(index) < FlatKnots(index + 1) - U ) {
681 other_index = index - 1 ;
682 }
683 else {
684 other_index = index + 1 ;
685 }
686 type = 3 ;
687
688 start_num_poles = StartingCondition + 2 ;
689 end_num_poles = num_poles - EndingCondition - 1 ;
690 if (start_num_poles == 1) {
691 start_value = schoenberg_points(num_poles) - schoenberg_points(1) ;
692 start_value = schoenberg_points(1) - start_value ;
693 }
694 else {
695 start_value = schoenberg_points(start_num_poles - 1) ;
696 }
697 if (end_num_poles == num_poles) {
698 end_value = schoenberg_points(num_poles) - schoenberg_points(1) ;
699 end_value = schoenberg_points(num_poles) + end_value ;
700 }
701 else {
702 end_value = schoenberg_points(end_num_poles + 1) ;
703 }
704
705 for (ii = 1 ; ii < start_num_poles ; ii++) {
706 first_function(ii) = 0.e0 ;
707 second_function(ii) = 0.0e0 ;
708 }
709
710 for (ii = end_num_poles + 1 ; ii <= num_poles ; ii++) {
711 first_function(ii) = 0.e0 ;
712 second_function(ii) = 0.0e0 ;
713 }
714 divide = 1.0e0 / (schoenberg_points(index) - start_value) ;
715
716 for (ii = start_num_poles ; ii <= index ; ii++) {
717 value = schoenberg_points(ii) - start_value ;
718 value *= divide ;
719 first_function(ii) = 1.0e0 ;
720
721 for (jj = 0 ; jj < type ; jj++) {
722 first_function(ii) *= value ;
723 }
724 }
725 divide = 1.0e0 /(end_value - schoenberg_points(index)) ;
726
727 for (ii = index ; ii <= end_num_poles ; ii++) {
728 value = end_value - schoenberg_points(ii) ;
729 value *= divide ;
730 first_function(ii) = 1.0e0 ;
731
732 for (jj = 0 ; jj < type ; jj++) {
733 first_function(ii) *= value ;
734 }
735 }
736
737 divide = 1.0e0 / (schoenberg_points(other_index) - start_value) ;
738
739 for (ii = start_num_poles ; ii <= other_index ; ii++) {
740 value = schoenberg_points(ii) - start_value ;
741 value *= divide ;
742 second_function(ii) = 1.0e0 ;
743
744 for (jj = 0 ; jj < type ; jj++) {
745 second_function(ii) *= value ;
746 }
747 }
748 divide = 1.0e0/( end_value - schoenberg_points(other_index)) ;
749
750 for (ii = other_index ; ii <= end_num_poles ; ii++) {
751 value = end_value - schoenberg_points(ii) ;
752 value *= divide ;
753 second_function(ii) = 1.0e0 ;
754
755 for (jj = 0 ; jj < type ; jj++) {
756 second_function(ii) *= value ;
757 }
758 }
759
760 //
761 // compute the point and derivatives of both functions
762 //
763 Standard_Real results[2][2],
764 weights_results[2][2];
765 Standard_Integer extrap_mode[2],
766 derivative_request = 1,
767 dimension = 1 ;
768 Standard_Boolean periodic_flag = Standard_False ;
769
770 extrap_mode[0] = Degree ;
771 extrap_mode[1] = Degree ;
772 if (Rational) {
773 //
774 // evaluate in homogenised form
775 //
776 Eval(U,
777 periodic_flag,
778 derivative_request,
779 extrap_mode[0],
780 Degree,
781 FlatKnots,
782 dimension,
783 first_function(1),
784 weights_array[0],
785 results[0][0],
786 weights_results[0][0]) ;
787
788 Eval(U,
789 periodic_flag,
790 derivative_request,
791 extrap_mode[0],
792 Degree,
793 FlatKnots,
794 dimension,
795 second_function(1),
796 weights_array[0],
797 results[1][0],
798 weights_results[1][0]) ;
799 //
800 // compute the rational derivatives values
801 //
802
803 for (ii = 0 ; ii < 2 ; ii++) {
804 PLib::RationalDerivatives(1,
805 1,
806 results[ii][0],
807 weights_results[ii][0],
808 results[ii][0]) ;
809 }
810 }
811 else {
812 Eval(U,
813 Standard_False,
814 1,
815 extrap_mode[0],
816 Degree,
817 FlatKnots,
818 1,
819 first_function(1),
820 results[0][0]) ;
821
822 Eval(U,
823 Standard_False,
824 1,
825 extrap_mode[0],
826 Degree,
827 FlatKnots,
828 1,
829 second_function(1),
830 results[1][0]) ;
831 }
832 gp_Mat2d a_matrix ;
833
834 for (ii = 0 ; ii < 2 ; ii++) {
835
836 for (jj = 0 ; jj < 2 ; jj++) {
837 a_matrix.SetValue(ii+1,jj+1,results[ii][jj]) ;
838 }
839 }
840 a_matrix.Invert() ;
841 TColStd_Array1OfReal the_a_vector(0,ArrayDimension-1) ;
842 TColStd_Array1OfReal the_b_vector(0,ArrayDimension-1) ;
843
844 for( ii = 0 ; ii < ArrayDimension ; ii++) {
845 the_a_vector(ii) =
846 a_matrix.Value(1,1) * delta_array[ii] +
847 a_matrix.Value(2,1) * derivatives_array[ii] ;
848 the_b_vector(ii) =
849 a_matrix.Value(1,2) * delta_array[ii] +
850 a_matrix.Value(2,2) * derivatives_array[ii] ;
851 }
852 index = 0 ;
853
854 for (ii = 0 ; ii < num_poles ; ii++) {
855
856 for (jj = 0 ; jj < ArrayDimension ; jj++) {
857 new_poles_array[index] = poles_array[index] ;
858 new_poles_array[index] +=
859 first_function(ii+1) * the_a_vector(jj) ;
860 new_poles_array[index] +=
861 second_function(ii+1) * the_b_vector(jj) ;
862 index += 1 ;
863 }
864 }
865 }
866 else {
867 ErrorStatus = 1 ;
868 }
869 }
870 else {
871 ErrorStatus = 2 ;
872 }
873}
874
875//=======================================================================
876//function : FunctionMultiply
877//purpose :
878//=======================================================================
879
880void BSplCLib::FunctionMultiply
881(const BSplCLib_EvaluatorFunction & FunctionPtr,
882 const Standard_Integer BSplineDegree,
883 const TColStd_Array1OfReal & BSplineFlatKnots,
884 const Standard_Integer PolesDimension,
885 Standard_Real & Poles,
886 const TColStd_Array1OfReal & FlatKnots,
887 const Standard_Integer NewDegree,
888 Standard_Real & NewPoles,
889 Standard_Integer & Status)
890{
891 Standard_Integer ii,
892 jj,
893 index ;
894 Standard_Integer extrap_mode[2],
895 error_code,
896 num_new_poles,
897 derivative_request = 0 ;
898 Standard_Boolean periodic_flag = Standard_False ;
899 Standard_Real result,
900 start_end[2],
901 *array_of_poles,
902 *array_of_new_poles ;
903
904 array_of_poles = (Standard_Real *) &NewPoles ;
905 extrap_mode[0] =
906 extrap_mode[1] = BSplineDegree ;
907 num_new_poles =
908 FlatKnots.Length() - NewDegree - 1 ;
909 start_end[0] = FlatKnots(NewDegree+1) ;
910 start_end[1] = FlatKnots(num_new_poles+1) ;
911 TColStd_Array1OfReal parameters(1,num_new_poles) ;
912 TColStd_Array1OfInteger contact_order_array(1,num_new_poles) ;
913 TColStd_Array1OfReal new_poles_array(1, num_new_poles * PolesDimension) ;
914
915 array_of_new_poles =
916 (Standard_Real *) &new_poles_array(1) ;
917 BuildSchoenbergPoints(NewDegree,
918 FlatKnots,
919 parameters) ;
920 //
921 // on recadre sur les bornes
922 //
923 if (parameters(1) < start_end[0]) {
924 parameters(1) = start_end[0] ;
925 }
926 if (parameters(num_new_poles) > start_end[1]) {
927 parameters(num_new_poles) = start_end[1] ;
928 }
929 index = 0 ;
930
931 for (ii = 1 ; ii <= num_new_poles ; ii++) {
932 contact_order_array(ii) = 0 ;
41194117 933 FunctionPtr.Evaluate (contact_order_array(ii),
7fd59977 934 start_end,
935 parameters(ii),
936 result,
41194117 937 error_code);
7fd59977 938 if (error_code) {
939 Status = 1 ;
940 goto FINISH ;
941 }
942
943 Eval(parameters(ii),
944 periodic_flag,
945 derivative_request,
946 extrap_mode[0],
947 BSplineDegree,
948 BSplineFlatKnots,
949 PolesDimension,
950 Poles,
951 array_of_new_poles[index]) ;
952
953 for (jj = 0 ; jj < PolesDimension ; jj++) {
954 array_of_new_poles[index] *= result ;
955 index += 1 ;
956 }
957 }
958 Interpolate(NewDegree,
959 FlatKnots,
960 parameters,
961 contact_order_array,
962 PolesDimension,
963 array_of_new_poles[0],
964 Status) ;
965
966 for (ii = 0 ; ii < num_new_poles * PolesDimension ; ii++) {
967 array_of_poles[ii] = array_of_new_poles[ii] ;
968
969 }
970 FINISH :
971 ;
972}
973
974//=======================================================================
975// function : FunctionMultiply
976//purpose :
977//=======================================================================
978
979void BSplCLib::FunctionReparameterise
980(const BSplCLib_EvaluatorFunction & FunctionPtr,
981 const Standard_Integer BSplineDegree,
982 const TColStd_Array1OfReal & BSplineFlatKnots,
983 const Standard_Integer PolesDimension,
984 Standard_Real & Poles,
985 const TColStd_Array1OfReal & FlatKnots,
986 const Standard_Integer NewDegree,
987 Standard_Real & NewPoles,
988 Standard_Integer & Status)
989{
990 Standard_Integer ii,
991// jj,
992 index ;
993 Standard_Integer extrap_mode[2],
994 error_code,
995 num_new_poles,
996 derivative_request = 0 ;
997 Standard_Boolean periodic_flag = Standard_False ;
998 Standard_Real result,
999 start_end[2],
1000 *array_of_poles,
1001 *array_of_new_poles ;
1002
1003 array_of_poles = (Standard_Real *) &NewPoles ;
1004 extrap_mode[0] =
1005 extrap_mode[1] = BSplineDegree ;
1006 num_new_poles =
1007 FlatKnots.Length() - NewDegree - 1 ;
1008 start_end[0] = FlatKnots(NewDegree+1) ;
1009 start_end[1] = FlatKnots(num_new_poles+1) ;
1010 TColStd_Array1OfReal parameters(1,num_new_poles) ;
1011 TColStd_Array1OfInteger contact_order_array(1,num_new_poles) ;
1012 TColStd_Array1OfReal new_poles_array(1, num_new_poles * PolesDimension) ;
1013
1014 array_of_new_poles =
1015 (Standard_Real *) &new_poles_array(1) ;
1016 BuildSchoenbergPoints(NewDegree,
1017 FlatKnots,
1018 parameters) ;
1019 index = 0 ;
1020
1021 for (ii = 1 ; ii <= num_new_poles ; ii++) {
1022 contact_order_array(ii) = 0 ;
41194117 1023 FunctionPtr.Evaluate (contact_order_array(ii),
7fd59977 1024 start_end,
1025 parameters(ii),
1026 result,
41194117 1027 error_code);
7fd59977 1028 if (error_code) {
1029 Status = 1 ;
1030 goto FINISH ;
1031 }
1032
1033 Eval(result,
1034 periodic_flag,
1035 derivative_request,
1036 extrap_mode[0],
1037 BSplineDegree,
1038 BSplineFlatKnots,
1039 PolesDimension,
1040 Poles,
1041 array_of_new_poles[index]) ;
1042 index += PolesDimension ;
1043 }
1044 Interpolate(NewDegree,
1045 FlatKnots,
1046 parameters,
1047 contact_order_array,
1048 PolesDimension,
1049 array_of_new_poles[0],
1050 Status) ;
1051
1052 for (ii = 0 ; ii < num_new_poles * PolesDimension ; ii++) {
1053 array_of_poles[ii] = array_of_new_poles[ii] ;
1054
1055 }
1056 FINISH :
1057 ;
1058}
1059
1060//=======================================================================
1061//function : FunctionMultiply
1062//purpose :
1063//=======================================================================
1064
1065void BSplCLib::FunctionMultiply
1066(const BSplCLib_EvaluatorFunction & FunctionPtr,
1067 const Standard_Integer BSplineDegree,
1068 const TColStd_Array1OfReal & BSplineFlatKnots,
1069 const TColStd_Array1OfReal & Poles,
1070 const TColStd_Array1OfReal & FlatKnots,
1071 const Standard_Integer NewDegree,
1072 TColStd_Array1OfReal & NewPoles,
1073 Standard_Integer & Status)
1074{
1075 Standard_Integer num_bspline_poles =
1076 BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1077 Standard_Integer num_new_poles =
1078 FlatKnots.Length() - NewDegree - 1 ;
1079
1080 if (Poles.Length() != num_bspline_poles ||
1081 NewPoles.Length() != num_new_poles) {
1082 Standard_ConstructionError::Raise();
1083 }
1084 Standard_Real * array_of_poles =
1085 (Standard_Real *) &Poles(Poles.Lower()) ;
1086 Standard_Real * array_of_new_poles =
1087 (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1088 BSplCLib::FunctionMultiply(FunctionPtr,
1089 BSplineDegree,
1090 BSplineFlatKnots,
1091 1,
1092 array_of_poles[0],
1093 FlatKnots,
1094 NewDegree,
1095 array_of_new_poles[0],
1096 Status) ;
1097}
1098
1099//=======================================================================
1100//function : FunctionReparameterise
1101//purpose :
1102//=======================================================================
1103
1104void BSplCLib::FunctionReparameterise
1105(const BSplCLib_EvaluatorFunction & FunctionPtr,
1106 const Standard_Integer BSplineDegree,
1107 const TColStd_Array1OfReal & BSplineFlatKnots,
1108 const TColStd_Array1OfReal & Poles,
1109 const TColStd_Array1OfReal & FlatKnots,
1110 const Standard_Integer NewDegree,
1111 TColStd_Array1OfReal & NewPoles,
1112 Standard_Integer & Status)
1113{
1114 Standard_Integer num_bspline_poles =
1115 BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1116 Standard_Integer num_new_poles =
1117 FlatKnots.Length() - NewDegree - 1 ;
1118
1119 if (Poles.Length() != num_bspline_poles ||
1120 NewPoles.Length() != num_new_poles) {
1121 Standard_ConstructionError::Raise();
1122 }
1123 Standard_Real * array_of_poles =
1124 (Standard_Real *) &Poles(Poles.Lower()) ;
1125 Standard_Real * array_of_new_poles =
1126 (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1127 BSplCLib::FunctionReparameterise(
1128 FunctionPtr,
1129 BSplineDegree,
1130 BSplineFlatKnots,
1131 1,
1132 array_of_poles[0],
1133 FlatKnots,
1134 NewDegree,
1135 array_of_new_poles[0],
1136 Status) ;
1137}
1138
1139//=======================================================================
1140//function : MergeBSplineKnots
1141//purpose :
1142//=======================================================================
1143void BSplCLib::MergeBSplineKnots
1144(const Standard_Real Tolerance,
1145 const Standard_Real StartValue,
1146 const Standard_Real EndValue,
1147 const Standard_Integer Degree1,
1148 const TColStd_Array1OfReal & Knots1,
1149 const TColStd_Array1OfInteger & Mults1,
1150 const Standard_Integer Degree2,
1151 const TColStd_Array1OfReal & Knots2,
1152 const TColStd_Array1OfInteger & Mults2,
1153 Standard_Integer & NumPoles,
1154 Handle(TColStd_HArray1OfReal) & NewKnots,
1155 Handle(TColStd_HArray1OfInteger) & NewMults)
1156{
1157 Standard_Integer ii,
1158 jj,
1159 continuity,
1160 set_mults_flag,
1161 degree,
1162 index,
1163 num_knots ;
1164 if (StartValue < EndValue - Tolerance) {
1165 TColStd_Array1OfReal knots1(1,Knots1.Length()) ;
1166 TColStd_Array1OfReal knots2(1,Knots2.Length()) ;
1167 degree = Degree1 + Degree2 ;
1168 index = 1 ;
1169
1170 for (ii = Knots1.Lower() ; ii <= Knots1.Upper() ; ii++) {
1171 knots1(index) = Knots1(ii) ;
1172 index += 1 ;
1173 }
1174 index = 1 ;
1175
1176 for (ii = Knots2.Lower() ; ii <= Knots2.Upper() ; ii++) {
1177 knots2(index) = Knots2(ii) ;
1178 index += 1 ;
1179 }
1180 BSplCLib::Reparametrize(StartValue,
1181 EndValue,
1182 knots1) ;
1183
1184 BSplCLib::Reparametrize(StartValue,
1185 EndValue,
1186 knots2) ;
1187 num_knots = 0 ;
1188 jj = 1 ;
1189
1190 for (ii = 1 ; ii <= knots1.Length() ; ii++) {
1191
1192 while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) - Tolerance) {
1193 jj += 1 ;
1194 num_knots += 1 ;
1195 }
1196
1197 while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) + Tolerance) {
1198 jj += 1 ;
1199 }
1200 num_knots += 1 ;
1201 }
1202 NewKnots =
1203 new TColStd_HArray1OfReal(1,num_knots) ;
1204 NewMults =
1205 new TColStd_HArray1OfInteger(1,num_knots) ;
1206 num_knots = 1 ;
1207 jj = 1 ;
1208
1209 for (ii = 1 ; ii <= knots1.Length() ; ii++) {
1210
1211 while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) - Tolerance) {
1212 NewKnots->ChangeArray1()(num_knots) = knots2(jj) ;
1213 NewMults->ChangeArray1()(num_knots) = Mults2(jj) + Degree1 ;
1214 jj += 1 ;
1215 num_knots += 1 ;
1216 }
1217 set_mults_flag = 0 ;
1218
1219 while (jj <= knots2.Length() && knots2(jj) <= knots1(ii) + Tolerance) {
1220 continuity = Min(Degree1 - Mults1(ii), Degree2 - Mults2(jj)) ;
1221 set_mults_flag = 1 ;
1222 NewMults->ChangeArray1()(num_knots) = degree - continuity ;
1223 jj += 1 ;
1224 }
1225
1226 NewKnots->ChangeArray1()(num_knots) = knots1(ii) ;
1227 if (! set_mults_flag) {
1228 NewMults->ChangeArray1()(num_knots) = Mults1(ii) + Degree2 ;
1229 }
1230 num_knots += 1 ;
1231 }
1232 num_knots -= 1 ;
1233 NewMults->ChangeArray1()(1) = degree + 1 ;
1234 NewMults->ChangeArray1()(num_knots) = degree + 1 ;
1235 index = 0 ;
1236
1237 for (ii = 1 ; ii <= num_knots ; ii++) {
1238 index += NewMults->Value(ii) ;
1239 }
1240 NumPoles = index - degree - 1 ;
1241 }
1242}
1243
Open CASCADE Technology repository