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