4019dacce7f61300c4c19ed6c00c68ec0c1a106c
[occt.git] / src / Geom / Geom_BSplineCurve.cxx
1 // Created on: 1993-03-09
2 // Created by: JCV
3 // Copyright (c) 1993-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
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
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.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 //Avril 1991 : constructeurs + methodes de lecture.
18 //Mai 1991   : revue des specifs + debut de realisation des classes tool =>
19 //             implementation des methodes Set et calcul du point courant.
20 //Juillet 1991 : voir egalement File Geom_BSplineCurve_1.cxx
21 //Juin    1992 : mise a plat des valeurs nodales - amelioration des
22 //               performances sur calcul du point courant
23
24 //RLE Aug 1993  Remove Swaps, Remove typedefs, Update BSplCLib
25 //              debug periodic, IncreaseDegree
26 // 21-Mar-95 : xab implemented cache
27 // 14-Mar-96 : xab implemented MovePointAndTangent 
28 // 13-Oct-96 : pmn Bug dans SetPeriodic (PRO6088) et Segment (PRO6250)
29
30 #define No_Standard_OutOfRange
31
32 #include <Geom_BSplineCurve.ixx>
33 #include <gp.hxx>
34 #include <ElCLib.hxx>
35 #include <BSplCLib.hxx>
36 #include <BSplCLib_KnotDistribution.hxx>
37 #include <BSplCLib_MultDistribution.hxx>
38 #include <Standard_NotImplemented.hxx>
39 #include <Standard_ConstructionError.hxx>
40 #include <Standard_OutOfRange.hxx>
41 #include <Standard_Real.hxx>
42
43 //=======================================================================
44 //function : CheckCurveData
45 //purpose  : Internal use only
46 //=======================================================================
47
48 static void CheckCurveData
49 (const TColgp_Array1OfPnt&         CPoles,
50  const TColStd_Array1OfReal&       CKnots,
51  const TColStd_Array1OfInteger&    CMults,
52  const Standard_Integer            Degree,
53  const Standard_Boolean            Periodic)
54 {
55   if (Degree < 1 || Degree > Geom_BSplineCurve::MaxDegree()) {
56     Standard_ConstructionError::Raise();  
57   }
58   
59   if (CPoles.Length() < 2)                Standard_ConstructionError::Raise();
60   if (CKnots.Length() != CMults.Length()) Standard_ConstructionError::Raise();
61   
62   for (Standard_Integer I = CKnots.Lower(); I < CKnots.Upper(); I++) {
63     if (CKnots (I+1) - CKnots (I) <= Epsilon (Abs(CKnots (I)))) {
64       Standard_ConstructionError::Raise();
65     }
66   }
67   
68   if (CPoles.Length() != BSplCLib::NbPoles(Degree,Periodic,CMults))
69     Standard_ConstructionError::Raise();
70 }
71
72 //=======================================================================
73 //function : Rational
74 //purpose  : check rationality of an array of weights
75 //=======================================================================
76
77 static Standard_Boolean Rational(const TColStd_Array1OfReal& W)
78 {
79   Standard_Integer i, n = W.Length();
80   Standard_Boolean rat = Standard_False;
81   for (i = 1; i < n; i++) {
82     rat =  Abs(W(i) - W(i+1)) > gp::Resolution();
83     if (rat) break;
84   }
85   return rat;
86 }
87
88 //=======================================================================
89 //function : Copy
90 //purpose  : 
91 //=======================================================================
92
93 Handle(Geom_Geometry) Geom_BSplineCurve::Copy() const
94 {
95   Handle(Geom_BSplineCurve) C;
96   if (IsRational()) 
97     C = new Geom_BSplineCurve(poles->Array1(),
98                               weights->Array1(),
99                               knots->Array1(),
100                               mults->Array1(),
101                               deg,periodic);
102   else
103     C = new Geom_BSplineCurve(poles->Array1(),
104                               knots->Array1(),
105                               mults->Array1(),
106                               deg,periodic);
107   return C;
108 }
109
110 //=======================================================================
111 //function : Geom_BSplineCurve
112 //purpose  : 
113 //=======================================================================
114
115 Geom_BSplineCurve::Geom_BSplineCurve
116 (const TColgp_Array1OfPnt&       Poles,
117  const TColStd_Array1OfReal&     Knots,
118  const TColStd_Array1OfInteger&  Mults,
119  const Standard_Integer          Degree,
120  const Standard_Boolean          Periodic) :
121  rational(Standard_False),
122  periodic(Periodic),
123  deg(Degree),
124  maxderivinvok(Standard_False)
125 {
126   // check
127   
128   CheckCurveData (Poles,
129                   Knots,
130                   Mults,
131                   Degree,
132                   Periodic);
133
134
135   // copy arrays
136
137   poles =  new TColgp_HArray1OfPnt(1,Poles.Length());
138   poles->ChangeArray1() = Poles;
139   
140
141   knots = new TColStd_HArray1OfReal(1,Knots.Length());
142   knots->ChangeArray1() = Knots;
143
144   mults = new TColStd_HArray1OfInteger(1,Mults.Length());
145   mults->ChangeArray1() = Mults;
146
147   UpdateKnots();
148   cachepoles = new TColgp_HArray1OfPnt(1,Degree + 1);
149   parametercache = 0.0e0 ;
150   spanlenghtcache = 0.0e0 ;
151   spanindexcache = 0 ;
152
153 }
154
155 //=======================================================================
156 //function : Geom_BSplineCurve
157 //purpose  : 
158 //=======================================================================
159
160 Geom_BSplineCurve::Geom_BSplineCurve
161 (const TColgp_Array1OfPnt&      Poles,
162  const TColStd_Array1OfReal&    Weights,
163  const TColStd_Array1OfReal&    Knots,
164  const TColStd_Array1OfInteger& Mults,
165  const Standard_Integer         Degree,
166  const Standard_Boolean         Periodic,
167  const Standard_Boolean         CheckRational)  :
168  rational(Standard_True),
169  periodic(Periodic),
170  deg(Degree),
171  maxderivinvok(Standard_False)
172
173 {
174
175   // check
176   
177   CheckCurveData (Poles,
178                   Knots,
179                   Mults,
180                   Degree,
181                   Periodic);
182
183   if (Weights.Length() != Poles.Length())
184     Standard_ConstructionError::Raise("Geom_BSplineCurve");
185
186   Standard_Integer i;
187   for (i = Weights.Lower(); i <= Weights.Upper(); i++) {
188     if (Weights(i) <= gp::Resolution())  
189       Standard_ConstructionError::Raise("Geom_BSplineCurve");
190   }
191   
192   // check really rational
193   if (CheckRational)
194     rational = Rational(Weights);
195   
196   // copy arrays
197   
198   poles =  new TColgp_HArray1OfPnt(1,Poles.Length());
199   poles->ChangeArray1() = Poles;
200   cachepoles = new TColgp_HArray1OfPnt(1,Degree + 1);
201   if (rational) {
202     weights =  new TColStd_HArray1OfReal(1,Weights.Length());
203     weights->ChangeArray1() = Weights;
204     cacheweights  = new TColStd_HArray1OfReal(1,Degree + 1);
205   }
206
207   knots = new TColStd_HArray1OfReal(1,Knots.Length());
208   knots->ChangeArray1() = Knots;
209
210   mults = new TColStd_HArray1OfInteger(1,Mults.Length());
211   mults->ChangeArray1() = Mults;
212
213   UpdateKnots();
214   parametercache = 0.0e0 ;
215   spanlenghtcache = 0.0e0 ;
216   spanindexcache = 0 ;
217 }
218
219 //=======================================================================
220 //function : MaxDegree
221 //purpose  : 
222 //=======================================================================
223
224 Standard_Integer Geom_BSplineCurve::MaxDegree () 
225
226   return BSplCLib::MaxDegree(); 
227 }
228
229 //=======================================================================
230 //function : IncreaseDegree
231 //purpose  : 
232 //=======================================================================
233
234 void Geom_BSplineCurve::IncreaseDegree  (const Standard_Integer Degree)
235 {
236   if (Degree == deg) return;
237   
238   if (Degree < deg || Degree > Geom_BSplineCurve::MaxDegree()) {
239     Standard_ConstructionError::Raise();
240   }
241   Standard_Integer FromK1 = FirstUKnotIndex ();
242   Standard_Integer ToK2   = LastUKnotIndex  ();
243   
244   Standard_Integer Step   = Degree - deg;
245   
246   Handle(TColgp_HArray1OfPnt) npoles = new
247     TColgp_HArray1OfPnt(1,poles->Length() + Step * (ToK2-FromK1));
248
249   Standard_Integer nbknots = BSplCLib::IncreaseDegreeCountKnots
250     (deg,Degree,periodic,mults->Array1());
251
252   Handle(TColStd_HArray1OfReal) nknots = 
253     new TColStd_HArray1OfReal(1,nbknots);
254
255   Handle(TColStd_HArray1OfInteger) nmults = 
256     new TColStd_HArray1OfInteger(1,nbknots);
257
258   Handle(TColStd_HArray1OfReal) nweights;
259   
260   if (IsRational()) {
261     
262     nweights = new TColStd_HArray1OfReal(1,npoles->Upper());
263     
264     BSplCLib::IncreaseDegree
265       (deg,Degree, periodic,
266        poles->Array1(),weights->Array1(),
267        knots->Array1(),mults->Array1(),
268        npoles->ChangeArray1(),nweights->ChangeArray1(),
269        nknots->ChangeArray1(),nmults->ChangeArray1());
270   }
271   else {
272     BSplCLib::IncreaseDegree
273       (deg,Degree, periodic,
274        poles->Array1(),BSplCLib::NoWeights(),
275        knots->Array1(),mults->Array1(),
276        npoles->ChangeArray1(),
277        *((TColStd_Array1OfReal*) NULL),
278        nknots->ChangeArray1(),nmults->ChangeArray1());
279   }
280
281   deg     = Degree;
282   poles   = npoles;
283   weights = nweights;
284   knots   = nknots;
285   mults   = nmults;
286   UpdateKnots();
287   
288 }
289
290 //=======================================================================
291 //function : IncreaseMultiplicity
292 //purpose  : 
293 //=======================================================================
294
295 void Geom_BSplineCurve::IncreaseMultiplicity  (const Standard_Integer Index,
296                                                const Standard_Integer M)
297 {
298   TColStd_Array1OfReal k(1,1);
299   k(1) = knots->Value(Index);
300   TColStd_Array1OfInteger m(1,1);
301   m(1) = M - mults->Value(Index);
302   InsertKnots(k,m,Epsilon(1.),Standard_True);
303 }
304
305 //=======================================================================
306 //function : IncreaseMultiplicity
307 //purpose  : 
308 //=======================================================================
309
310 void Geom_BSplineCurve::IncreaseMultiplicity  (const Standard_Integer I1,
311                                                const Standard_Integer I2,
312                                                const Standard_Integer M)
313 {
314   Handle(TColStd_HArray1OfReal)  tk = knots;
315   TColStd_Array1OfReal k((knots->Array1())(I1),I1,I2);
316   TColStd_Array1OfInteger m(I1,I2);
317   Standard_Integer i;
318   for (i = I1; i <= I2; i++)
319     m(i) = M - mults->Value(i);
320   InsertKnots(k,m,Epsilon(1.),Standard_True);
321 }
322
323 //=======================================================================
324 //function : IncrementMultiplicity
325 //purpose  : 
326 //=======================================================================
327
328 void Geom_BSplineCurve::IncrementMultiplicity
329 (const Standard_Integer I1,
330  const Standard_Integer I2,
331  const Standard_Integer Step)
332 {
333   Handle(TColStd_HArray1OfReal) tk = knots;
334   TColStd_Array1OfReal    k((knots->Array1())(I1),I1,I2);
335   TColStd_Array1OfInteger m(I1,I2) ;
336   m.Init(Step);
337   InsertKnots(k,m,Epsilon(1.),Standard_True);
338 }
339
340 //=======================================================================
341 //function : InsertKnot
342 //purpose  : 
343 //=======================================================================
344
345 void Geom_BSplineCurve::InsertKnot
346 (const Standard_Real U, 
347  const Standard_Integer M, 
348  const Standard_Real ParametricTolerance,
349  const Standard_Boolean Add)
350 {
351   TColStd_Array1OfReal k(1,1);
352   k(1) = U;
353   TColStd_Array1OfInteger m(1,1);
354   m(1) = M;
355   InsertKnots(k,m,ParametricTolerance,Add);
356 }
357
358 //=======================================================================
359 //function : InsertKnots
360 //purpose  : 
361 //=======================================================================
362
363 void  Geom_BSplineCurve::InsertKnots(const TColStd_Array1OfReal& Knots, 
364                                      const TColStd_Array1OfInteger& Mults,
365                                      const Standard_Real Epsilon,
366                                      const Standard_Boolean Add)
367 {
368   // Check and compute new sizes
369   Standard_Integer nbpoles,nbknots;
370
371   if (!BSplCLib::PrepareInsertKnots(deg,periodic,
372                                     knots->Array1(),mults->Array1(),
373                                     Knots,Mults,nbpoles,nbknots,Epsilon,Add))
374     Standard_ConstructionError::Raise("Geom_BSplineCurve::InsertKnots");
375
376   if (nbpoles == poles->Length()) return;
377
378   Handle(TColgp_HArray1OfPnt)      npoles = new TColgp_HArray1OfPnt(1,nbpoles);
379   Handle(TColStd_HArray1OfReal)    nknots = knots;
380   Handle(TColStd_HArray1OfInteger) nmults = mults;
381
382   if (nbknots != knots->Length()) {
383     nknots = new TColStd_HArray1OfReal(1,nbknots);
384     nmults = new TColStd_HArray1OfInteger(1,nbknots);
385   }
386
387   if (rational) {
388     Handle(TColStd_HArray1OfReal) nweights = 
389       new TColStd_HArray1OfReal(1,nbpoles);
390     BSplCLib::InsertKnots(deg,periodic,
391                           poles->Array1(), weights->Array1(),
392                           knots->Array1(), mults->Array1(),
393                           Knots, Mults,
394                           npoles->ChangeArray1(), nweights->ChangeArray1(),
395                           nknots->ChangeArray1(), nmults->ChangeArray1(),
396                           Epsilon, Add);
397     weights = nweights;
398   }
399   else {
400     BSplCLib::InsertKnots(deg,periodic,
401                           poles->Array1(), BSplCLib::NoWeights(),
402                           knots->Array1(), mults->Array1(),
403                           Knots, Mults,
404                           npoles->ChangeArray1(),
405                           *((TColStd_Array1OfReal*) NULL),
406                           nknots->ChangeArray1(), nmults->ChangeArray1(),
407                           Epsilon, Add);
408   }
409
410   poles = npoles;
411   knots = nknots;
412   mults = nmults;
413   UpdateKnots();
414
415 }
416
417 //=======================================================================
418 //function : RemoveKnot
419 //purpose  : 
420 //=======================================================================
421
422 Standard_Boolean  Geom_BSplineCurve::RemoveKnot(const Standard_Integer Index,
423                                                 const Standard_Integer M, 
424                                                 const Standard_Real Tolerance)
425 {
426   if (M < 0) return Standard_True;
427
428   Standard_Integer I1  = FirstUKnotIndex ();
429   Standard_Integer I2  = LastUKnotIndex  ();
430
431   if ( !periodic && (Index <= I1 || Index >= I2) ) {
432     Standard_OutOfRange::Raise();
433   }
434   else if ( periodic  && (Index < I1 || Index > I2)) {
435     Standard_OutOfRange::Raise();
436   }
437   
438   const TColgp_Array1OfPnt   & oldpoles   = poles->Array1();
439
440   Standard_Integer step = mults->Value(Index) - M;
441   if (step <= 0) return Standard_True;
442
443   Handle(TColgp_HArray1OfPnt) npoles =
444     new TColgp_HArray1OfPnt(1,oldpoles.Length()-step);
445
446   Handle(TColStd_HArray1OfReal)    nknots  = knots;
447   Handle(TColStd_HArray1OfInteger) nmults  = mults;
448
449   if (M == 0) {
450     nknots = new TColStd_HArray1OfReal(1,knots->Length()-1);
451     nmults = new TColStd_HArray1OfInteger(1,knots->Length()-1);
452   }
453
454   if (IsRational()) {
455     Handle(TColStd_HArray1OfReal) nweights = 
456       new TColStd_HArray1OfReal(1,npoles->Length());
457     if (!BSplCLib::RemoveKnot
458         (Index, M, deg, periodic,
459          poles->Array1(),weights->Array1(), 
460          knots->Array1(),mults->Array1(),
461          npoles->ChangeArray1(), nweights->ChangeArray1(),
462          nknots->ChangeArray1(),nmults->ChangeArray1(),
463          Tolerance))
464       return Standard_False;
465     weights = nweights;
466   }
467   else {
468     if (!BSplCLib::RemoveKnot
469         (Index, M, deg, periodic,
470          poles->Array1(), BSplCLib::NoWeights(),
471          knots->Array1(),mults->Array1(),
472          npoles->ChangeArray1(),
473          *((TColStd_Array1OfReal*) NULL),
474          nknots->ChangeArray1(),nmults->ChangeArray1(),
475          Tolerance))
476       return Standard_False;
477   }
478   
479   poles = npoles;
480   knots = nknots;
481   mults = nmults;
482
483   UpdateKnots();
484   maxderivinvok = 0;
485   return Standard_True;
486 }
487
488 //=======================================================================
489 //function : Reverse
490 //purpose  : 
491 //=======================================================================
492
493 void Geom_BSplineCurve::Reverse ()
494
495   BSplCLib::Reverse(knots->ChangeArray1());
496   BSplCLib::Reverse(mults->ChangeArray1());
497   Standard_Integer last;
498   if (periodic)
499     last = flatknots->Upper() - deg - 1;
500   else
501     last = poles->Upper();
502   BSplCLib::Reverse(poles->ChangeArray1(),last);
503   if (rational)
504     BSplCLib::Reverse(weights->ChangeArray1(),last);
505   UpdateKnots();
506 }
507
508 //=======================================================================
509 //function : ReversedParameter
510 //purpose  : 
511 //=======================================================================
512
513 Standard_Real Geom_BSplineCurve::ReversedParameter
514 (const Standard_Real U) const
515 {
516   return (FirstParameter() + LastParameter() - U);
517 }
518
519 //=======================================================================
520 //function : Segment
521 //purpose  : 
522 //=======================================================================
523
524 void Geom_BSplineCurve::Segment(const Standard_Real U1,
525                                 const Standard_Real U2)
526 {
527   Standard_DomainError_Raise_if ( U2 < U1,
528                                  "Geom_BSplineCurve::Segment");
529   
530   Standard_Real NewU1, NewU2;
531   Standard_Real U,DU=0,aDDU=0;
532   Standard_Integer index;
533   Standard_Boolean wasPeriodic = periodic;
534
535   TColStd_Array1OfReal    Knots(1,2);
536   TColStd_Array1OfInteger Mults(1,2);
537
538   // define param ditance to keep (eap, Apr 18 2002, occ311)
539   if (periodic) {
540     Standard_Real Period = LastParameter() - FirstParameter();
541     DU = U2 - U1;
542     while (DU > Period)
543       DU -= Period;
544     if (DU <= Epsilon(Period))
545       DU = Period;
546     aDDU = DU;
547   }
548
549   index = 0;
550   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
551                             U1,periodic,knots->Lower(),knots->Upper(),
552                             index,NewU1);
553   index = 0;
554   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
555                             U2,periodic,knots->Lower(),knots->Upper(),
556                             index,NewU2);
557
558   //-- DBB
559   Standard_Real aNu2 = NewU2;
560   //-- DBB
561
562
563   Knots( 1) = Min( NewU1, NewU2);
564   Knots( 2) = Max( NewU1, NewU2);
565   Mults( 1) = Mults( 2) = deg;
566
567   Standard_Real AbsUMax = Max(Abs(NewU1),Abs(NewU2));
568
569 //  Modified by Sergey KHROMOV - Fri Apr 11 12:15:40 2003 Begin
570   AbsUMax = Max(AbsUMax, Max(Abs(FirstParameter()),Abs(LastParameter())));
571 //  Modified by Sergey KHROMOV - Fri Apr 11 12:15:40 2003 End
572
573   Standard_Real Eps = 100. * Epsilon(AbsUMax);
574
575   InsertKnots( Knots, Mults, Eps);
576
577   if (periodic) { // set the origine at NewU1
578     index = 0;
579     BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
580                               U1,periodic,knots->Lower(),knots->Upper(),
581                               index,U);
582     // Test si l'insertion est Ok et decalage sinon. 
583     if ( Abs(knots->Value(index+1)-U) <= Eps) // <= pour etre homogene a InsertKnots
584       index++;
585     SetOrigin(index);
586     SetNotPeriodic();
587     NewU2 = NewU1 + DU;
588   }
589
590   // compute index1 and index2 to set the new knots and mults 
591   Standard_Integer index1 = 0, index2 = 0;
592   Standard_Integer FromU1 = knots->Lower();
593   Standard_Integer ToU2   = knots->Upper();
594   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
595                             NewU1,periodic,FromU1,ToU2,index1,U);
596   if ( Abs(knots->Value(index1+1)-U) <= Eps)
597     index1++;
598
599   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
600                             NewU2,periodic,FromU1,ToU2,index2,U);
601   if ( Abs(knots->Value(index2+1)-U) <= Eps)
602     index2++;
603   
604   Standard_Integer nbknots = index2 - index1 + 1;
605
606   Handle(TColStd_HArray1OfReal) 
607     nknots = new TColStd_HArray1OfReal(1,nbknots);
608   Handle(TColStd_HArray1OfInteger) 
609     nmults = new TColStd_HArray1OfInteger(1,nbknots);
610
611   // to restore changed U1
612   if (DU > 0) // if was periodic
613     DU = NewU1 - U1;
614   
615   Standard_Integer i , k = 1;
616   for ( i = index1; i<= index2; i++) {
617     nknots->SetValue(k, knots->Value(i) - DU);
618     nmults->SetValue(k, mults->Value(i));
619     k++;
620   }
621   nmults->SetValue(      1, deg + 1);
622   nmults->SetValue(nbknots, deg + 1);
623
624
625   // compute index1 and index2 to set the new poles and weights
626   Standard_Integer pindex1 
627     = BSplCLib::PoleIndex(deg,index1,periodic,mults->Array1());
628   Standard_Integer pindex2 
629     = BSplCLib::PoleIndex(deg,index2,periodic,mults->Array1());
630
631   pindex1++;
632   pindex2 = Min( pindex2+1, poles->Length());
633
634   Standard_Integer nbpoles  = pindex2 - pindex1 + 1;
635
636   Handle(TColStd_HArray1OfReal) 
637     nweights = new TColStd_HArray1OfReal(1,nbpoles);
638   Handle(TColgp_HArray1OfPnt)
639     npoles = new TColgp_HArray1OfPnt(1,nbpoles);
640
641   k = 1;
642   if ( rational) {
643     nweights = new TColStd_HArray1OfReal( 1, nbpoles);
644     for ( i = pindex1; i <= pindex2; i++) {
645       npoles->SetValue(k, poles->Value(i));
646       nweights->SetValue(k, weights->Value(i));
647       k++;
648     }
649   }
650   else {
651     for ( i = pindex1; i <= pindex2; i++) {
652       npoles->SetValue(k, poles->Value(i));
653       k++;
654     }
655   }
656
657   //-- DBB
658   if( wasPeriodic ) {
659     nknots->ChangeValue(nknots->Lower()) = U1;
660     if( aNu2 < U2 ) {
661       nknots->ChangeValue(nknots->Upper()) = U1 + aDDU;
662     }
663   }
664   //-- DBB
665
666   knots = nknots;
667   mults = nmults;
668   poles = npoles;
669   if ( rational) 
670     weights = nweights;
671
672   maxderivinvok = 0;
673   UpdateKnots();
674 }
675
676 //=======================================================================
677 //function : SetKnot
678 //purpose  : 
679 //=======================================================================
680
681 void Geom_BSplineCurve::SetKnot
682 (const Standard_Integer Index,
683  const Standard_Real K)
684 {
685   if (Index < 1 || Index > knots->Length())     Standard_OutOfRange::Raise();
686   Standard_Real DK = Abs(Epsilon (K));
687   if (Index == 1) { 
688     if (K >= knots->Value(2) - DK) Standard_ConstructionError::Raise();
689   }
690   else if (Index == knots->Length()) {
691     if (K <= knots->Value (knots->Length()-1) + DK)  {
692       Standard_ConstructionError::Raise();
693     }
694   }
695   else {
696     if (K <= knots->Value(Index-1) + DK ||
697         K >= knots->Value(Index+1) - DK ) {
698       Standard_ConstructionError::Raise();
699     }
700   }
701   if (K != knots->Value (Index)) {
702     knots->SetValue (Index, K);
703     maxderivinvok = 0;
704     UpdateKnots();
705   }
706 }
707
708 //=======================================================================
709 //function : SetKnots
710 //purpose  : 
711 //=======================================================================
712
713 void Geom_BSplineCurve::SetKnots
714 (const TColStd_Array1OfReal& K)
715 {
716   CheckCurveData(poles->Array1(),K,mults->Array1(),deg,periodic);
717   knots->ChangeArray1() = K;
718   maxderivinvok = 0;
719   UpdateKnots();
720 }
721
722 //=======================================================================
723 //function : SetKnot
724 //purpose  : 
725 //=======================================================================
726
727 void Geom_BSplineCurve::SetKnot
728 (const Standard_Integer Index,
729  const Standard_Real K,
730  const Standard_Integer M)
731 {
732   IncreaseMultiplicity (Index, M);
733   SetKnot (Index, K);
734 }
735
736 //=======================================================================
737 //function : SetPeriodic
738 //purpose  : 
739 //=======================================================================
740
741 void Geom_BSplineCurve::SetPeriodic ()
742 {
743   Standard_Integer first = FirstUKnotIndex();
744   Standard_Integer last  = LastUKnotIndex();
745
746   Handle(TColStd_HArray1OfReal) tk = knots;
747   TColStd_Array1OfReal    cknots((knots->Array1())(first),first,last);
748   knots = new TColStd_HArray1OfReal(1,cknots.Length());
749   knots->ChangeArray1() = cknots;
750
751   Handle(TColStd_HArray1OfInteger) tm = mults;
752   TColStd_Array1OfInteger cmults((mults->Array1())(first),first,last);
753   cmults(first) = cmults(last) = Min(deg, Max( cmults(first), cmults(last)));
754   mults = new TColStd_HArray1OfInteger(1,cmults.Length());
755   mults->ChangeArray1() = cmults;
756
757   // compute new number of poles;
758   Standard_Integer nbp = BSplCLib::NbPoles(deg,Standard_True,cmults);
759   
760   Handle(TColgp_HArray1OfPnt) tp = poles;
761   TColgp_Array1OfPnt cpoles((poles->Array1())(1),1,nbp);
762   poles = new TColgp_HArray1OfPnt(1,nbp);
763   poles->ChangeArray1() = cpoles;
764   
765   if (rational) {
766     Handle(TColStd_HArray1OfReal) tw = weights;
767     TColStd_Array1OfReal cweights((weights->Array1())(1),1,nbp);
768     weights = new TColStd_HArray1OfReal(1,nbp);
769     weights->ChangeArray1() = cweights;
770   }
771
772   periodic = Standard_True;
773   
774   maxderivinvok = 0;
775   UpdateKnots();
776 }
777
778 //=======================================================================
779 //function : SetOrigin
780 //purpose  : 
781 //=======================================================================
782
783 void Geom_BSplineCurve::SetOrigin(const Standard_Integer Index)
784 {
785   Standard_NoSuchObject_Raise_if( !periodic,
786                                  "Geom_BSplineCurve::SetOrigin");
787   Standard_Integer i,k;
788   Standard_Integer first = FirstUKnotIndex();
789   Standard_Integer last  = LastUKnotIndex();
790
791   Standard_DomainError_Raise_if( (Index < first) || (Index > last),
792                                 "Geom_BSplineCurve::SetOrigine");
793
794   Standard_Integer nbknots = knots->Length();
795   Standard_Integer nbpoles = poles->Length();
796
797   Handle(TColStd_HArray1OfReal) nknots = 
798     new TColStd_HArray1OfReal(1,nbknots);
799   TColStd_Array1OfReal& newknots = nknots->ChangeArray1();
800
801   Handle(TColStd_HArray1OfInteger) nmults =
802     new TColStd_HArray1OfInteger(1,nbknots);
803   TColStd_Array1OfInteger& newmults = nmults->ChangeArray1();
804
805   // set the knots and mults
806   Standard_Real period = knots->Value(last) - knots->Value(first);
807   k = 1;
808   for ( i = Index; i <= last ; i++) {
809     newknots(k) = knots->Value(i);
810     newmults(k) = mults->Value(i);
811     k++;
812   }
813   for ( i = first+1; i <= Index; i++) {
814     newknots(k) = knots->Value(i) + period;
815     newmults(k) = mults->Value(i);
816     k++;
817   }
818
819   Standard_Integer index = 1;
820   for (i = first+1; i <= Index; i++) 
821     index += mults->Value(i);
822
823   // set the poles and weights
824   Handle(TColgp_HArray1OfPnt) npoles =
825     new TColgp_HArray1OfPnt(1,nbpoles);
826   Handle(TColStd_HArray1OfReal) nweights =
827     new TColStd_HArray1OfReal(1,nbpoles);
828   TColgp_Array1OfPnt   & newpoles   = npoles->ChangeArray1();
829   TColStd_Array1OfReal & newweights = nweights->ChangeArray1();
830   first = poles->Lower();
831   last  = poles->Upper();
832   if ( rational) {
833     k = 1;
834     for ( i = index; i <= last; i++) {
835       newpoles(k)   = poles->Value(i);
836       newweights(k) = weights->Value(i);
837       k++;
838     }
839     for ( i = first; i < index; i++) {
840       newpoles(k)   = poles->Value(i);
841       newweights(k) = weights->Value(i);
842       k++;
843     }
844   }
845   else {
846     k = 1;
847     for ( i = index; i <= last; i++) {
848       newpoles(k) = poles->Value(i);
849       k++;
850     }
851     for ( i = first; i < index; i++) {
852       newpoles(k) = poles->Value(i);
853       k++;
854     }
855   }
856
857   poles = npoles;
858   knots = nknots;
859   mults = nmults;
860   if (rational) 
861     weights = nweights;
862   maxderivinvok = 0;
863   UpdateKnots();
864 }
865
866 //=======================================================================
867 //function : SetOrigin
868 //purpose  : 
869 //=======================================================================
870
871 void Geom_BSplineCurve::SetOrigin(const Standard_Real U,
872                                   const Standard_Real Tol)
873 {
874   Standard_NoSuchObject_Raise_if( !periodic,
875                                  "Geom_BSplineCurve::SetOrigin");
876   //U est il dans la period.
877   Standard_Real uf = FirstParameter(), ul = LastParameter();
878   Standard_Real u = U, period = ul - uf;
879   while (Tol < (uf-u)) u += period;
880   while (Tol > (ul-u)) u -= period;
881
882   if(Abs(U-u)>Tol) { //On reparametre la courbe
883     Standard_Real delta = U-u;
884     uf += delta;
885     ul += delta;
886     TColStd_Array1OfReal& kn = knots->ChangeArray1();
887     Standard_Integer fk = kn.Lower(), lk = kn.Upper();
888     for(Standard_Integer i = fk; i <= lk; i++){
889       kn.ChangeValue(i) += delta;
890     }
891     UpdateKnots();
892   }
893   if(Abs(U-uf)<Tol) return;
894   
895   TColStd_Array1OfReal& kn = knots->ChangeArray1();
896   Standard_Integer fk = kn.Lower(), lk = kn.Upper(),ik=0;
897   Standard_Real delta = RealLast();
898   for(Standard_Integer i = fk; i<= lk; i++){
899     Standard_Real dki = kn.Value(i)-U;
900     if(Abs(dki)<Abs(delta)){
901       ik = i;
902       delta = dki;
903     }
904   }
905   if(Abs(delta)>Tol){
906     InsertKnot(U);
907     if(delta < 0.) ik++;
908   }
909   SetOrigin(ik);
910 }
911
912 //=======================================================================
913 //function : SetNotPeriodic
914 //purpose  : 
915 //=======================================================================
916
917 void Geom_BSplineCurve::SetNotPeriodic () 
918
919   if ( periodic) {
920     Standard_Integer NbKnots, NbPoles;
921     BSplCLib::PrepareUnperiodize( deg, mults->Array1(),NbKnots,NbPoles);
922     
923     Handle(TColgp_HArray1OfPnt) npoles 
924       = new TColgp_HArray1OfPnt(1,NbPoles);
925     
926     Handle(TColStd_HArray1OfReal) nknots 
927       = new TColStd_HArray1OfReal(1,NbKnots);
928     
929     Handle(TColStd_HArray1OfInteger) nmults
930       = new TColStd_HArray1OfInteger(1,NbKnots);
931     
932     Handle(TColStd_HArray1OfReal) nweights;
933     
934     if (IsRational()) {
935       
936       nweights = new TColStd_HArray1OfReal(1,NbPoles);
937       
938       BSplCLib::Unperiodize
939         (deg,mults->Array1(),knots->Array1(),poles->Array1(),
940          weights->Array1(),nmults->ChangeArray1(),
941          nknots->ChangeArray1(),npoles->ChangeArray1(),
942          nweights->ChangeArray1());
943       
944     }
945     else {
946       
947       BSplCLib::Unperiodize
948         (deg,mults->Array1(),knots->Array1(),poles->Array1(),
949          BSplCLib::NoWeights(),nmults->ChangeArray1(),
950          nknots->ChangeArray1(),npoles->ChangeArray1(),
951          *((TColStd_Array1OfReal*) NULL));
952       
953     }
954     poles   = npoles;
955     weights = nweights;
956     mults   = nmults;
957     knots   = nknots;
958     periodic = Standard_False;
959     
960     maxderivinvok = 0;
961     UpdateKnots();
962   }
963 }
964
965 //=======================================================================
966 //function : SetPole
967 //purpose  : 
968 //=======================================================================
969
970 void Geom_BSplineCurve::SetPole
971 (const Standard_Integer Index,
972  const gp_Pnt& P)
973 {
974   if (Index < 1 || Index > poles->Length()) Standard_OutOfRange::Raise();
975   poles->SetValue (Index, P);
976   maxderivinvok = 0;
977   InvalidateCache() ;
978 }
979
980 //=======================================================================
981 //function : SetPole
982 //purpose  : 
983 //=======================================================================
984
985 void Geom_BSplineCurve::SetPole
986 (const Standard_Integer Index,
987  const gp_Pnt& P,
988  const Standard_Real W)
989 {
990   SetPole(Index,P);
991   SetWeight(Index,W);
992 }
993
994 //=======================================================================
995 //function : SetWeight
996 //purpose  : 
997 //=======================================================================
998
999 void Geom_BSplineCurve::SetWeight
1000 (const Standard_Integer Index,
1001  const Standard_Real W)
1002 {
1003   if (Index < 1 || Index > poles->Length())   Standard_OutOfRange::Raise();
1004
1005   if (W <= gp::Resolution ())     Standard_ConstructionError::Raise();
1006
1007
1008   Standard_Boolean rat = IsRational() || (Abs(W - 1.) > gp::Resolution());
1009
1010   if ( rat) {
1011     if (rat && !IsRational()) {
1012       weights = new TColStd_HArray1OfReal(1,poles->Length());
1013       weights->Init(1.);
1014     }
1015     
1016     weights->SetValue (Index, W);
1017     
1018     if (IsRational()) {
1019       rat = Rational(weights->Array1());
1020       if (!rat) weights.Nullify();
1021     }
1022     
1023     rational = !weights.IsNull();
1024   }
1025   maxderivinvok = 0;
1026   InvalidateCache() ;
1027 }
1028
1029 //=======================================================================
1030 //function : MovePoint
1031 //purpose  : 
1032 //=======================================================================
1033
1034 void Geom_BSplineCurve::MovePoint(const Standard_Real U,
1035                                   const gp_Pnt& P,
1036                                   const Standard_Integer Index1,
1037                                   const Standard_Integer Index2,
1038                                   Standard_Integer& FirstModifiedPole,
1039                                   Standard_Integer& LastmodifiedPole)
1040 {
1041   if (Index1 < 1 || Index1 > poles->Length() || 
1042       Index2 < 1 || Index2 > poles->Length() || Index1 > Index2) {
1043     Standard_OutOfRange::Raise();
1044   }
1045   TColgp_Array1OfPnt npoles(1, poles->Length());
1046   gp_Pnt P0;
1047   D0(U, P0);
1048   gp_Vec Displ(P0, P);
1049   BSplCLib::MovePoint(U, Displ, Index1, Index2, deg, rational, poles->Array1(), 
1050                       weights->Array1(), flatknots->Array1(), 
1051                       FirstModifiedPole, LastmodifiedPole, npoles);
1052   if (FirstModifiedPole) {
1053     poles->ChangeArray1() = npoles;
1054     maxderivinvok = 0;
1055     InvalidateCache() ;
1056   }
1057 }
1058
1059 //=======================================================================
1060 //function : MovePointAndTangent
1061 //purpose  : 
1062 //=======================================================================
1063
1064 void Geom_BSplineCurve::
1065 MovePointAndTangent(const Standard_Real U,
1066                     const gp_Pnt&       P,
1067                     const gp_Vec&       Tangent,
1068                     const Standard_Real    Tolerance,
1069                     const Standard_Integer StartingCondition,
1070                     const Standard_Integer EndingCondition,
1071                     Standard_Integer&      ErrorStatus) 
1072 {
1073   Standard_Integer ii ;
1074   if (IsPeriodic()) {
1075     //
1076     // for the time being do not deal with periodic curves
1077     //
1078     SetNotPeriodic() ;
1079   }
1080   TColgp_Array1OfPnt new_poles(1, poles->Length());
1081   gp_Pnt P0;
1082
1083
1084   gp_Vec delta_derivative;
1085   D1(U, P0,
1086      delta_derivative) ;
1087   gp_Vec delta(P0, P);
1088   for (ii = 1 ; ii <= 3 ; ii++) {
1089     delta_derivative.SetCoord(ii, 
1090                               Tangent.Coord(ii)- delta_derivative.Coord(ii)) ;
1091   }
1092   BSplCLib::MovePointAndTangent(U,
1093                                 delta,
1094                                 delta_derivative,
1095                                 Tolerance,
1096                                 deg,
1097                                 rational,
1098                                 StartingCondition,
1099                                 EndingCondition,
1100                                 poles->Array1(), 
1101                                 weights->Array1(), 
1102                                 flatknots->Array1(), 
1103                                 new_poles,
1104                                 ErrorStatus) ;
1105   if (!ErrorStatus) {
1106     poles->ChangeArray1() = new_poles;
1107     maxderivinvok = 0;
1108     InvalidateCache() ;
1109   }
1110 }
1111
1112 //=======================================================================
1113 //function : UpdateKnots
1114 //purpose  : 
1115 //=======================================================================
1116
1117 void Geom_BSplineCurve::UpdateKnots()
1118 {
1119   rational = !weights.IsNull();
1120
1121   Standard_Integer MaxKnotMult = 0;
1122   BSplCLib::KnotAnalysis (deg, 
1123                 periodic,
1124                 knots->Array1(), 
1125                 mults->Array1(), 
1126                 knotSet, MaxKnotMult);
1127   
1128   if (knotSet == GeomAbs_Uniform && !periodic)  {
1129     flatknots = knots;
1130   }
1131   else {
1132     flatknots = new TColStd_HArray1OfReal 
1133       (1, BSplCLib::KnotSequenceLength(mults->Array1(),deg,periodic));
1134
1135     BSplCLib::KnotSequence (knots->Array1(), 
1136                             mults->Array1(),
1137                             deg,periodic,
1138                             flatknots->ChangeArray1());
1139   }
1140   
1141   if (MaxKnotMult == 0)  smooth = GeomAbs_CN;
1142   else {
1143     switch (deg - MaxKnotMult) {
1144     case 0 :   smooth = GeomAbs_C0;   break;
1145     case 1 :   smooth = GeomAbs_C1;   break;
1146     case 2 :   smooth = GeomAbs_C2;   break;
1147     case 3 :   smooth = GeomAbs_C3;   break;
1148       default :  smooth = GeomAbs_C3;   break;
1149     }
1150   }
1151   InvalidateCache() ;
1152 }
1153
1154 //=======================================================================
1155 //function : Invalidate the Cache
1156 //purpose  : as the name says
1157 //=======================================================================
1158
1159 void Geom_BSplineCurve::InvalidateCache() 
1160 {
1161   validcache = 0 ;
1162 }
1163
1164 //=======================================================================
1165 //function : check if the Cache is valid
1166 //purpose  : as the name says
1167 //=======================================================================
1168
1169 Standard_Boolean Geom_BSplineCurve::IsCacheValid
1170 (const Standard_Real  U)  const 
1171 {
1172   //Roman Lygin 26.12.08, performance improvements
1173   //1. avoided using NewParameter = (U - parametercache) / spanlenghtcache
1174   //to check against [0, 1), as division is CPU consuming
1175   //2. minimized use of if, as branching is also CPU consuming
1176   Standard_Real aDelta = U - parametercache;
1177
1178   return ( validcache &&
1179       (aDelta >= 0.0e0) &&
1180       ((aDelta < spanlenghtcache) || (spanindexcache == flatknots->Upper() - deg)) );
1181 }
1182
1183 //=======================================================================
1184 //function : Normalizes the parameters if the curve is periodic
1185 //purpose  : that is compute the cache so that it is valid
1186 //=======================================================================
1187
1188 void Geom_BSplineCurve::PeriodicNormalization(Standard_Real&  Parameter) const 
1189 {
1190   Standard_Real Period ;
1191
1192   if (periodic) {
1193     Period = flatknots->Value(flatknots->Upper() - deg) - flatknots->Value (deg + 1) ;
1194     while (Parameter > flatknots->Value(flatknots->Upper()-deg)) {
1195       Parameter -= Period ;
1196     }
1197     while (Parameter < flatknots->Value((deg + 1))) {
1198       Parameter +=  Period ;
1199     }
1200   }
1201 }
1202
1203 //=======================================================================
1204 //function : Validate the Cache
1205 //purpose  : that is compute the cache so that it is valid
1206 //=======================================================================
1207
1208 void Geom_BSplineCurve::ValidateCache(const Standard_Real  Parameter) 
1209 {
1210   Standard_Real NewParameter ;
1211   Standard_Integer LocalIndex = 0 ;
1212   //
1213   // check if the degree did not change
1214   //
1215   if (cachepoles->Upper() < deg + 1) {
1216     cachepoles = new TColgp_HArray1OfPnt(1,deg + 1);
1217     if (rational) {
1218       cacheweights  = new TColStd_HArray1OfReal(1,deg + 1);
1219     }
1220   }
1221   BSplCLib::LocateParameter(deg,
1222                             (flatknots->Array1()),
1223                             (BSplCLib::NoMults()),
1224                             Parameter,
1225                             periodic,
1226                             LocalIndex,
1227                             NewParameter);
1228   spanindexcache = LocalIndex ;
1229   if (Parameter == flatknots->Value(LocalIndex + 1)) {
1230     
1231     LocalIndex += 1 ;
1232     parametercache = flatknots->Value(LocalIndex) ;
1233     if (LocalIndex == flatknots->Upper() - deg) {
1234       //
1235       // for the last span if the parameter is outside of 
1236       // the domain of the curve than use the last knot
1237       // and normalize with the last span Still set the
1238       // spanindexcache to flatknots->Upper() - deg so that
1239       // the IsCacheValid will know for sure we are extending
1240       // the Bspline 
1241       //
1242       
1243       spanlenghtcache = flatknots->Value(LocalIndex - 1) - parametercache ;
1244     }
1245     else {
1246       spanlenghtcache = flatknots->Value(LocalIndex + 1) - parametercache ;
1247     }
1248   }
1249   else {
1250     parametercache = flatknots->Value(LocalIndex) ;
1251     spanlenghtcache = flatknots->Value(LocalIndex + 1) - parametercache ;
1252   }
1253   
1254   if  (rational) {
1255     BSplCLib::BuildCache(parametercache,
1256                          spanlenghtcache,
1257                          periodic,
1258                          deg,
1259                          (flatknots->Array1()),
1260                          poles->Array1(),
1261                          weights->Array1(),
1262                          cachepoles->ChangeArray1(),
1263                          cacheweights->ChangeArray1()) ;
1264   }
1265   else {
1266     BSplCLib::BuildCache(parametercache,
1267                          spanlenghtcache,
1268                          periodic,
1269                          deg,
1270                          (flatknots->Array1()),
1271                          poles->Array1(),
1272                          *((TColStd_Array1OfReal*) NULL),
1273                          cachepoles->ChangeArray1(),
1274                          *((TColStd_Array1OfReal*) NULL)) ;
1275   }
1276   validcache = 1 ;
1277 }
1278