7792defa11878abb06823f44560f5cf9336cc3a7
[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 {
496   if (U2 < U1)
497     throw Standard_DomainError("Geom_BSplineCurve::Segment");
498   
499   Standard_Real NewU1, NewU2;
500   Standard_Real U,DU=0,aDDU=0;
501   Standard_Integer index;
502   Standard_Boolean wasPeriodic = periodic;
503
504   TColStd_Array1OfReal    Knots(1,2);
505   TColStd_Array1OfInteger Mults(1,2);
506
507   // define param ditance to keep (eap, Apr 18 2002, occ311)
508   if (periodic) {
509     Standard_Real Period = LastParameter() - FirstParameter();
510     DU = U2 - U1;
511     if (DU - Period > Precision::PConfusion())
512       throw Standard_DomainError("Geom_BSplineCurve::Segment");
513     if (DU > Period)
514       DU = Period;
515     aDDU = DU;
516   }
517
518   index = 0;
519   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
520                             U1,periodic,knots->Lower(),knots->Upper(),
521                             index,NewU1);
522   index = 0;
523   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
524                             U2,periodic,knots->Lower(),knots->Upper(),
525                             index,NewU2);
526
527   //-- DBB
528   Standard_Real aNu2 = NewU2;
529   //-- DBB
530
531
532   Knots( 1) = Min( NewU1, NewU2);
533   Knots( 2) = Max( NewU1, NewU2);
534   Mults( 1) = Mults( 2) = deg;
535
536   Standard_Real AbsUMax = Max(Abs(NewU1),Abs(NewU2));
537
538 //  Modified by Sergey KHROMOV - Fri Apr 11 12:15:40 2003 Begin
539   AbsUMax = Max(AbsUMax, Max(Abs(FirstParameter()),Abs(LastParameter())));
540 //  Modified by Sergey KHROMOV - Fri Apr 11 12:15:40 2003 End
541
542   Standard_Real Eps = 100. * Epsilon(AbsUMax);
543
544   InsertKnots( Knots, Mults, Eps);
545
546   if (periodic) { // set the origine at NewU1
547     index = 0;
548     BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
549                               U1,periodic,knots->Lower(),knots->Upper(),
550                               index,U);
551     // Test si l'insertion est Ok et decalage sinon. 
552     if ( Abs(knots->Value(index+1)-U) <= Eps) // <= pour etre homogene a InsertKnots
553       index++;
554     SetOrigin(index);
555     SetNotPeriodic();
556     NewU2 = NewU1 + DU;
557   }
558
559   // compute index1 and index2 to set the new knots and mults 
560   Standard_Integer index1 = 0, index2 = 0;
561   Standard_Integer FromU1 = knots->Lower();
562   Standard_Integer ToU2   = knots->Upper();
563   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
564                             NewU1,periodic,FromU1,ToU2,index1,U);
565   if ( Abs(knots->Value(index1+1)-U) <= Eps)
566     index1++;
567
568   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
569                             NewU2,periodic,FromU1,ToU2,index2,U);
570   if ( Abs(knots->Value(index2+1)-U) <= Eps || index2 == index1)
571     index2++;
572   
573   Standard_Integer nbknots = index2 - index1 + 1;
574
575   Handle(TColStd_HArray1OfReal) 
576     nknots = new TColStd_HArray1OfReal(1,nbknots);
577   Handle(TColStd_HArray1OfInteger) 
578     nmults = new TColStd_HArray1OfInteger(1,nbknots);
579
580   // to restore changed U1
581   if (DU > 0) // if was periodic
582     DU = NewU1 - U1;
583   
584   Standard_Integer i , k = 1;
585   for ( i = index1; i<= index2; i++) {
586     nknots->SetValue(k, knots->Value(i) - DU);
587     nmults->SetValue(k, mults->Value(i));
588     k++;
589   }
590   nmults->SetValue(      1, deg + 1);
591   nmults->SetValue(nbknots, deg + 1);
592
593
594   // compute index1 and index2 to set the new poles and weights
595   Standard_Integer pindex1 
596     = BSplCLib::PoleIndex(deg,index1,periodic,mults->Array1());
597   Standard_Integer pindex2 
598     = BSplCLib::PoleIndex(deg,index2,periodic,mults->Array1());
599
600   pindex1++;
601   pindex2 = Min( pindex2+1, poles->Length());
602
603   Standard_Integer nbpoles  = pindex2 - pindex1 + 1;
604
605   Handle(TColStd_HArray1OfReal) 
606     nweights = new TColStd_HArray1OfReal(1,nbpoles);
607   Handle(TColgp_HArray1OfPnt)
608     npoles = new TColgp_HArray1OfPnt(1,nbpoles);
609
610   k = 1;
611   if ( rational) {
612     nweights = new TColStd_HArray1OfReal( 1, nbpoles);
613     for ( i = pindex1; i <= pindex2; i++) {
614       npoles->SetValue(k, poles->Value(i));
615       nweights->SetValue(k, weights->Value(i));
616       k++;
617     }
618   }
619   else {
620     for ( i = pindex1; i <= pindex2; i++) {
621       npoles->SetValue(k, poles->Value(i));
622       k++;
623     }
624   }
625
626   //-- DBB
627   if( wasPeriodic ) {
628     nknots->ChangeValue(nknots->Lower()) = U1;
629     if( aNu2 < U2 ) {
630       nknots->ChangeValue(nknots->Upper()) = U1 + aDDU;
631     }
632   }
633   //-- DBB
634
635   knots = nknots;
636   mults = nmults;
637   poles = npoles;
638   if ( rational) 
639     weights = nweights;
640
641   maxderivinvok = 0;
642   UpdateKnots();
643 }
644
645 //=======================================================================
646 //function : SetKnot
647 //purpose  : 
648 //=======================================================================
649
650 void Geom_BSplineCurve::SetKnot
651 (const Standard_Integer Index,
652  const Standard_Real K)
653 {
654   if (Index < 1 || Index > knots->Length())     throw Standard_OutOfRange();
655   Standard_Real DK = Abs(Epsilon (K));
656   if (Index == 1) { 
657     if (K >= knots->Value(2) - DK) throw Standard_ConstructionError();
658   }
659   else if (Index == knots->Length()) {
660     if (K <= knots->Value (knots->Length()-1) + DK)  {
661       throw Standard_ConstructionError();
662     }
663   }
664   else {
665     if (K <= knots->Value(Index-1) + DK ||
666         K >= knots->Value(Index+1) - DK ) {
667       throw Standard_ConstructionError();
668     }
669   }
670   if (K != knots->Value (Index)) {
671     knots->SetValue (Index, K);
672     maxderivinvok = 0;
673     UpdateKnots();
674   }
675 }
676
677 //=======================================================================
678 //function : SetKnots
679 //purpose  : 
680 //=======================================================================
681
682 void Geom_BSplineCurve::SetKnots
683 (const TColStd_Array1OfReal& K)
684 {
685   CheckCurveData(poles->Array1(),K,mults->Array1(),deg,periodic);
686   knots->ChangeArray1() = K;
687   maxderivinvok = 0;
688   UpdateKnots();
689 }
690
691 //=======================================================================
692 //function : SetKnot
693 //purpose  : 
694 //=======================================================================
695
696 void Geom_BSplineCurve::SetKnot
697 (const Standard_Integer Index,
698  const Standard_Real K,
699  const Standard_Integer M)
700 {
701   IncreaseMultiplicity (Index, M);
702   SetKnot (Index, K);
703 }
704
705 //=======================================================================
706 //function : SetPeriodic
707 //purpose  : 
708 //=======================================================================
709
710 void Geom_BSplineCurve::SetPeriodic ()
711 {
712   Standard_Integer first = FirstUKnotIndex();
713   Standard_Integer last  = LastUKnotIndex();
714
715   Handle(TColStd_HArray1OfReal) tk = knots;
716   TColStd_Array1OfReal    cknots((knots->Array1())(first),first,last);
717   knots = new TColStd_HArray1OfReal(1,cknots.Length());
718   knots->ChangeArray1() = cknots;
719
720   Handle(TColStd_HArray1OfInteger) tm = mults;
721   TColStd_Array1OfInteger cmults((mults->Array1())(first),first,last);
722   cmults(first) = cmults(last) = Min(deg, Max( cmults(first), cmults(last)));
723   mults = new TColStd_HArray1OfInteger(1,cmults.Length());
724   mults->ChangeArray1() = cmults;
725
726   // compute new number of poles;
727   Standard_Integer nbp = BSplCLib::NbPoles(deg,Standard_True,cmults);
728   
729   Handle(TColgp_HArray1OfPnt) tp = poles;
730   TColgp_Array1OfPnt cpoles((poles->Array1())(1),1,nbp);
731   poles = new TColgp_HArray1OfPnt(1,nbp);
732   poles->ChangeArray1() = cpoles;
733   
734   if (rational) {
735     Handle(TColStd_HArray1OfReal) tw = weights;
736     TColStd_Array1OfReal cweights((weights->Array1())(1),1,nbp);
737     weights = new TColStd_HArray1OfReal(1,nbp);
738     weights->ChangeArray1() = cweights;
739   }
740
741   periodic = Standard_True;
742   
743   maxderivinvok = 0;
744   UpdateKnots();
745 }
746
747 //=======================================================================
748 //function : SetOrigin
749 //purpose  : 
750 //=======================================================================
751
752 void Geom_BSplineCurve::SetOrigin(const Standard_Integer Index)
753 {
754   if (!periodic)
755     throw Standard_NoSuchObject("Geom_BSplineCurve::SetOrigin");
756
757   Standard_Integer i,k;
758   Standard_Integer first = FirstUKnotIndex();
759   Standard_Integer last  = LastUKnotIndex();
760
761   if ((Index < first) || (Index > last))
762     throw Standard_DomainError("Geom_BSplineCurve::SetOrigin");
763
764   Standard_Integer nbknots = knots->Length();
765   Standard_Integer nbpoles = poles->Length();
766
767   Handle(TColStd_HArray1OfReal) nknots = 
768     new TColStd_HArray1OfReal(1,nbknots);
769   TColStd_Array1OfReal& newknots = nknots->ChangeArray1();
770
771   Handle(TColStd_HArray1OfInteger) nmults =
772     new TColStd_HArray1OfInteger(1,nbknots);
773   TColStd_Array1OfInteger& newmults = nmults->ChangeArray1();
774
775   // set the knots and mults
776   Standard_Real period = knots->Value(last) - knots->Value(first);
777   k = 1;
778   for ( i = Index; i <= last ; i++) {
779     newknots(k) = knots->Value(i);
780     newmults(k) = mults->Value(i);
781     k++;
782   }
783   for ( i = first+1; i <= Index; i++) {
784     newknots(k) = knots->Value(i) + period;
785     newmults(k) = mults->Value(i);
786     k++;
787   }
788
789   Standard_Integer index = 1;
790   for (i = first+1; i <= Index; i++) 
791     index += mults->Value(i);
792
793   // set the poles and weights
794   Handle(TColgp_HArray1OfPnt) npoles =
795     new TColgp_HArray1OfPnt(1,nbpoles);
796   Handle(TColStd_HArray1OfReal) nweights =
797     new TColStd_HArray1OfReal(1,nbpoles);
798   TColgp_Array1OfPnt   & newpoles   = npoles->ChangeArray1();
799   TColStd_Array1OfReal & newweights = nweights->ChangeArray1();
800   first = poles->Lower();
801   last  = poles->Upper();
802   if ( rational) {
803     k = 1;
804     for ( i = index; i <= last; i++) {
805       newpoles(k)   = poles->Value(i);
806       newweights(k) = weights->Value(i);
807       k++;
808     }
809     for ( i = first; i < index; i++) {
810       newpoles(k)   = poles->Value(i);
811       newweights(k) = weights->Value(i);
812       k++;
813     }
814   }
815   else {
816     k = 1;
817     for ( i = index; i <= last; i++) {
818       newpoles(k) = poles->Value(i);
819       k++;
820     }
821     for ( i = first; i < index; i++) {
822       newpoles(k) = poles->Value(i);
823       k++;
824     }
825   }
826
827   poles = npoles;
828   knots = nknots;
829   mults = nmults;
830   if (rational) 
831     weights = nweights;
832   maxderivinvok = 0;
833   UpdateKnots();
834 }
835
836 //=======================================================================
837 //function : SetOrigin
838 //purpose  : 
839 //=======================================================================
840
841 void Geom_BSplineCurve::SetOrigin(const Standard_Real U,
842                                   const Standard_Real Tol)
843 {
844   if (!periodic)
845     throw Standard_NoSuchObject("Geom_BSplineCurve::SetOrigin");
846   //U est il dans la period.
847   Standard_Real uf = FirstParameter(), ul = LastParameter();
848   Standard_Real u = U, period = ul - uf;
849   while (Tol < (uf-u)) u += period;
850   while (Tol > (ul-u)) u -= period;
851
852   if(Abs(U-u)>Tol) { //On reparametre la courbe
853     Standard_Real delta = U-u;
854     uf += delta;
855     ul += delta;
856     TColStd_Array1OfReal& kn = knots->ChangeArray1();
857     Standard_Integer fk = kn.Lower(), lk = kn.Upper();
858     for(Standard_Integer i = fk; i <= lk; i++){
859       kn.ChangeValue(i) += delta;
860     }
861     UpdateKnots();
862   }
863   if(Abs(U-uf)<Tol) return;
864   
865   TColStd_Array1OfReal& kn = knots->ChangeArray1();
866   Standard_Integer fk = kn.Lower(), lk = kn.Upper(),ik=0;
867   Standard_Real delta = RealLast();
868   for(Standard_Integer i = fk; i<= lk; i++){
869     Standard_Real dki = kn.Value(i)-U;
870     if(Abs(dki)<Abs(delta)){
871       ik = i;
872       delta = dki;
873     }
874   }
875   if(Abs(delta)>Tol){
876     InsertKnot(U);
877     if(delta < 0.) ik++;
878   }
879   SetOrigin(ik);
880 }
881
882 //=======================================================================
883 //function : SetNotPeriodic
884 //purpose  : 
885 //=======================================================================
886
887 void Geom_BSplineCurve::SetNotPeriodic () 
888
889   if ( periodic) {
890     Standard_Integer NbKnots, NbPoles;
891     BSplCLib::PrepareUnperiodize( deg, mults->Array1(),NbKnots,NbPoles);
892     
893     Handle(TColgp_HArray1OfPnt) npoles 
894       = new TColgp_HArray1OfPnt(1,NbPoles);
895     
896     Handle(TColStd_HArray1OfReal) nknots 
897       = new TColStd_HArray1OfReal(1,NbKnots);
898     
899     Handle(TColStd_HArray1OfInteger) nmults
900       = new TColStd_HArray1OfInteger(1,NbKnots);
901     
902     Handle(TColStd_HArray1OfReal) nweights;
903     if (IsRational())
904     {
905       nweights = new TColStd_HArray1OfReal(1,NbPoles);
906     }
907
908     BSplCLib::Unperiodize (deg,
909                            mults->Array1(), knots->Array1(), poles->Array1(),
910                            !nweights.IsNull() ? &weights->Array1() : BSplCLib::NoWeights(),
911                            nmults->ChangeArray1(), nknots->ChangeArray1(), npoles->ChangeArray1(),
912                            !nweights.IsNull() ? &nweights->ChangeArray1() : BSplCLib::NoWeights());
913     poles   = npoles;
914     weights = nweights;
915     mults   = nmults;
916     knots   = nknots;
917     periodic = Standard_False;
918     
919     maxderivinvok = 0;
920     UpdateKnots();
921   }
922 }
923
924 //=======================================================================
925 //function : SetPole
926 //purpose  : 
927 //=======================================================================
928
929 void Geom_BSplineCurve::SetPole
930 (const Standard_Integer Index,
931  const gp_Pnt& P)
932 {
933   if (Index < 1 || Index > poles->Length()) throw Standard_OutOfRange();
934   poles->SetValue (Index, P);
935   maxderivinvok = 0;
936 }
937
938 //=======================================================================
939 //function : SetPole
940 //purpose  : 
941 //=======================================================================
942
943 void Geom_BSplineCurve::SetPole
944 (const Standard_Integer Index,
945  const gp_Pnt& P,
946  const Standard_Real W)
947 {
948   SetPole(Index,P);
949   SetWeight(Index,W);
950 }
951
952 //=======================================================================
953 //function : SetWeight
954 //purpose  : 
955 //=======================================================================
956
957 void Geom_BSplineCurve::SetWeight
958 (const Standard_Integer Index,
959  const Standard_Real W)
960 {
961   if (Index < 1 || Index > poles->Length())   throw Standard_OutOfRange();
962
963   if (W <= gp::Resolution ())     throw Standard_ConstructionError();
964
965
966   Standard_Boolean rat = IsRational() || (Abs(W - 1.) > gp::Resolution());
967
968   if ( rat) {
969     if (rat && !IsRational()) {
970       weights = new TColStd_HArray1OfReal(1,poles->Length());
971       weights->Init(1.);
972     }
973     
974     weights->SetValue (Index, W);
975     
976     if (IsRational()) {
977       rat = Rational(weights->Array1());
978       if (!rat) weights.Nullify();
979     }
980     
981     rational = !weights.IsNull();
982   }
983   maxderivinvok = 0;
984 }
985
986 //=======================================================================
987 //function : MovePoint
988 //purpose  : 
989 //=======================================================================
990
991 void Geom_BSplineCurve::MovePoint(const Standard_Real U,
992                                   const gp_Pnt& P,
993                                   const Standard_Integer Index1,
994                                   const Standard_Integer Index2,
995                                   Standard_Integer& FirstModifiedPole,
996                                   Standard_Integer& LastmodifiedPole)
997 {
998   if (Index1 < 1 || Index1 > poles->Length() || 
999       Index2 < 1 || Index2 > poles->Length() || Index1 > Index2) {
1000     throw Standard_OutOfRange();
1001   }
1002   TColgp_Array1OfPnt npoles(1, poles->Length());
1003   gp_Pnt P0;
1004   D0(U, P0);
1005   gp_Vec Displ(P0, P);
1006   BSplCLib::MovePoint (U, Displ, Index1, Index2, deg, poles->Array1(),
1007                        rational ? &weights->Array1() : BSplCLib::NoWeights(), flatknots->Array1(),
1008                        FirstModifiedPole, LastmodifiedPole, npoles);
1009   if (FirstModifiedPole) {
1010     poles->ChangeArray1() = npoles;
1011     maxderivinvok = 0;
1012   }
1013 }
1014
1015 //=======================================================================
1016 //function : MovePointAndTangent
1017 //purpose  : 
1018 //=======================================================================
1019
1020 void Geom_BSplineCurve::MovePointAndTangent(const Standard_Real    U,
1021                                             const gp_Pnt&          P,
1022                                             const gp_Vec&          Tangent,
1023                                             const Standard_Real    Tolerance,
1024                                             const Standard_Integer StartingCondition,
1025                                             const Standard_Integer EndingCondition,
1026                                             Standard_Integer&      ErrorStatus) 
1027 {
1028   Standard_Integer ii ;
1029   if (IsPeriodic()) {
1030     //
1031     // for the time being do not deal with periodic curves
1032     //
1033     SetNotPeriodic() ;
1034   }
1035   TColgp_Array1OfPnt new_poles(1, poles->Length());
1036   gp_Pnt P0;
1037
1038
1039   gp_Vec delta_derivative;
1040   D1(U, P0,
1041      delta_derivative) ;
1042   gp_Vec delta(P0, P);
1043   for (ii = 1 ; ii <= 3 ; ii++) {
1044     delta_derivative.SetCoord(ii, Tangent.Coord(ii)-delta_derivative.Coord(ii));
1045   }
1046   BSplCLib::MovePointAndTangent(U,
1047                                 delta,
1048                                 delta_derivative,
1049                                 Tolerance,
1050                                 deg,
1051                                 StartingCondition,
1052                                 EndingCondition,
1053                                 poles->Array1(), 
1054                                 rational ? &weights->Array1() : BSplCLib::NoWeights(),
1055                                 flatknots->Array1(), 
1056                                 new_poles,
1057                                 ErrorStatus) ;
1058   if (!ErrorStatus) {
1059     poles->ChangeArray1() = new_poles;
1060     maxderivinvok = 0;
1061   }
1062 }
1063
1064 //=======================================================================
1065 //function : UpdateKnots
1066 //purpose  : 
1067 //=======================================================================
1068
1069 void Geom_BSplineCurve::UpdateKnots()
1070 {
1071   rational = !weights.IsNull();
1072
1073   Standard_Integer MaxKnotMult = 0;
1074   BSplCLib::KnotAnalysis(deg, 
1075                          periodic,
1076                          knots->Array1(), 
1077                          mults->Array1(), 
1078                          knotSet, MaxKnotMult);
1079   
1080   if (knotSet == GeomAbs_Uniform && !periodic)  {
1081     flatknots = knots;
1082   }
1083   else {
1084     flatknots = new TColStd_HArray1OfReal 
1085       (1, BSplCLib::KnotSequenceLength(mults->Array1(),deg,periodic));
1086
1087     BSplCLib::KnotSequence(knots->Array1(), 
1088                            mults->Array1(),
1089                            deg,periodic,
1090                            flatknots->ChangeArray1());
1091   }
1092   
1093   if (MaxKnotMult == 0)  smooth = GeomAbs_CN;
1094   else {
1095     switch (deg - MaxKnotMult) {
1096     case 0 :   smooth = GeomAbs_C0;   break;
1097     case 1 :   smooth = GeomAbs_C1;   break;
1098     case 2 :   smooth = GeomAbs_C2;   break;
1099     case 3 :   smooth = GeomAbs_C3;   break;
1100       default :  smooth = GeomAbs_C3;   break;
1101     }
1102   }
1103 }
1104
1105 //=======================================================================
1106 //function : Normalizes the parameters if the curve is periodic
1107 //purpose  : that is compute the cache so that it is valid
1108 //=======================================================================
1109
1110 void Geom_BSplineCurve::PeriodicNormalization(Standard_Real&  Parameter) const 
1111 {
1112   Standard_Real Period ;
1113
1114   if (periodic) {
1115     Period = flatknots->Value(flatknots->Upper() - deg) - flatknots->Value (deg + 1) ;
1116     while (Parameter > flatknots->Value(flatknots->Upper()-deg)) {
1117       Parameter -= Period ;
1118     }
1119     while (Parameter < flatknots->Value((deg + 1))) {
1120       Parameter +=  Period ;
1121     }
1122   }
1123 }
1124