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