1229c61c5c909028c4ffbcc0f10b9c2bffe9e284
[occt.git] / src / Geom2d / Geom2d_BSplineCurve.cxx
1 // Created on: 1993-03-25
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 Geom2d_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, Init methods, Remove typedefs
25 //  14-Mar-96 : xab implemented MovePointAndTangent
26
27 //SAMTECH Jan 2002 : add text to Raise()
28
29 #define No_Standard_OutOfRange
30
31
32 #include <BSplCLib.hxx>
33 #include <BSplCLib_KnotDistribution.hxx>
34 #include <BSplCLib_MultDistribution.hxx>
35 #include <Geom2d_BSplineCurve.hxx>
36 #include <Geom2d_Geometry.hxx>
37 #include <Geom2d_UndefinedDerivative.hxx>
38 #include <gp.hxx>
39 #include <gp_Pnt2d.hxx>
40 #include <gp_Trsf2d.hxx>
41 #include <gp_Vec2d.hxx>
42 #include <Precision.hxx>
43 #include <Standard_ConstructionError.hxx>
44 #include <Standard_DimensionError.hxx>
45 #include <Standard_DomainError.hxx>
46 #include <Standard_NoSuchObject.hxx>
47 #include <Standard_NotImplemented.hxx>
48 #include <Standard_OutOfRange.hxx>
49 #include <Standard_RangeError.hxx>
50 #include <Standard_Type.hxx>
51
52 IMPLEMENT_STANDARD_RTTIEXT(Geom2d_BSplineCurve,Geom2d_BoundedCurve)
53
54 //=======================================================================
55 //function : CheckCurveData
56 //purpose  : Internal use only
57 //=======================================================================
58 static void CheckCurveData
59 (const TColgp_Array1OfPnt2d&         CPoles,
60  const TColStd_Array1OfReal&       CKnots,
61  const TColStd_Array1OfInteger&    CMults,
62  const Standard_Integer            Degree,
63  const Standard_Boolean            Periodic)
64 {
65   if (Degree < 1 || Degree > Geom2d_BSplineCurve::MaxDegree()) {
66     throw Standard_ConstructionError("BSpline curve : invalid degree");
67   }
68   
69   if (CPoles.Length() < 2)                throw Standard_ConstructionError("BSpline curve : at least 2 poles required");
70   if (CKnots.Length() != CMults.Length()) throw Standard_ConstructionError("BSpline curve : Knot and Mult array size mismatch");
71   
72   for (Standard_Integer I = CKnots.Lower(); I < CKnots.Upper(); I++) {
73     if (CKnots (I+1) - CKnots (I) <= Epsilon (Abs(CKnots (I)))) {
74       throw Standard_ConstructionError("BSpline curve : Knots interval values too close");
75     }
76   }
77   
78   if (CPoles.Length() != BSplCLib::NbPoles(Degree,Periodic,CMults))
79     throw Standard_ConstructionError("BSpline curve : # Poles and degree mismatch");
80 }
81
82 //=======================================================================
83 //function : Rational
84 //purpose  : check rationality of an array of weights
85 //=======================================================================
86
87 static Standard_Boolean Rational(const TColStd_Array1OfReal& W)
88 {
89   Standard_Integer i, n = W.Length();
90   Standard_Boolean rat = Standard_False;
91   for (i = 1; i < n; i++) {
92     rat =  Abs(W(i) - W(i+1)) > gp::Resolution();
93     if (rat) break;
94   }
95   return rat;
96 }
97
98 //=======================================================================
99 //function : Copy
100 //purpose  : 
101 //=======================================================================
102
103 Handle(Geom2d_Geometry) Geom2d_BSplineCurve::Copy() const
104 {
105   Handle(Geom2d_BSplineCurve) C;
106   if (IsRational()) 
107     C = new Geom2d_BSplineCurve(poles->Array1(),
108                                 weights->Array1(),
109                                 knots->Array1(),
110                                 mults->Array1(),
111                                 deg,periodic);
112   else
113     C = new Geom2d_BSplineCurve(poles->Array1(),
114                                 knots->Array1(),
115                                 mults->Array1(),
116                                 deg,periodic);
117   return C;
118 }
119
120 //=======================================================================
121 //function : Geom2d_BSplineCurve
122 //purpose  : 
123 //=======================================================================
124
125 Geom2d_BSplineCurve::Geom2d_BSplineCurve
126 (const TColgp_Array1OfPnt2d&     Poles,
127  const TColStd_Array1OfReal&     Knots,
128  const TColStd_Array1OfInteger&  Mults,
129  const Standard_Integer          Degree,
130  const Standard_Boolean          Periodic) :
131  rational(Standard_False),
132  periodic(Periodic),
133  deg(Degree),
134  maxderivinvok(Standard_False)
135 {
136   // check
137   
138   CheckCurveData(Poles,
139                  Knots,
140                  Mults,
141                  Degree,
142                  Periodic);
143
144   // copy arrays
145
146   poles =  new TColgp_HArray1OfPnt2d(1,Poles.Length());
147   poles->ChangeArray1() = Poles;
148
149   knots = new TColStd_HArray1OfReal(1,Knots.Length());
150   knots->ChangeArray1() = Knots;
151
152   mults = new TColStd_HArray1OfInteger(1,Mults.Length());
153   mults->ChangeArray1() = Mults;
154
155   UpdateKnots();
156 }
157
158 //=======================================================================
159 //function : Geom2d_BSplineCurve
160 //purpose  : 
161 //=======================================================================
162
163 Geom2d_BSplineCurve::Geom2d_BSplineCurve
164 (const TColgp_Array1OfPnt2d&      Poles,
165  const TColStd_Array1OfReal&    Weights,
166  const TColStd_Array1OfReal&    Knots,
167  const TColStd_Array1OfInteger& Mults,
168  const Standard_Integer         Degree,
169  const Standard_Boolean         Periodic)  :
170  rational(Standard_True),
171  periodic(Periodic),
172  deg(Degree),
173  maxderivinvok(Standard_False)
174
175 {
176
177   // check
178   
179   CheckCurveData(Poles,
180                  Knots,
181                  Mults,
182                  Degree,
183                  Periodic);
184
185   if (Weights.Length() != Poles.Length())
186     throw Standard_ConstructionError("Geom2d_BSplineCurve :Weights and Poles array size mismatch");
187
188   Standard_Integer i;
189   for (i = Weights.Lower(); i <= Weights.Upper(); i++) {
190     if (Weights(i) <= gp::Resolution()) {
191       throw Standard_ConstructionError("Geom2d_BSplineCurve: Weights values too small");
192     }
193   }
194   
195   // check really rational
196   rational = Rational(Weights);
197   
198   // copy arrays
199   
200   poles =  new TColgp_HArray1OfPnt2d(1,Poles.Length());
201   poles->ChangeArray1() = Poles;
202   if (rational) {
203     weights =  new TColStd_HArray1OfReal(1,Weights.Length());
204     weights->ChangeArray1() = Weights;
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 }
215
216 //=======================================================================
217 //function : MaxDegree
218 //purpose  : 
219 //=======================================================================
220
221 Standard_Integer Geom2d_BSplineCurve::MaxDegree () 
222
223   return BSplCLib::MaxDegree(); 
224 }
225
226 //=======================================================================
227 //function : IncreaseDegree
228 //purpose  : 
229 //=======================================================================
230
231 void Geom2d_BSplineCurve::IncreaseDegree
232 (const Standard_Integer Degree)
233 {
234   if (Degree == deg) return;
235
236   if (Degree < deg || Degree > Geom2d_BSplineCurve::MaxDegree()) {
237     throw Standard_ConstructionError("BSpline curve : IncreaseDegree : bad degree value");
238   }
239
240   Standard_Integer FromK1 = FirstUKnotIndex ();
241   Standard_Integer ToK2   = LastUKnotIndex  ();
242   
243   Standard_Integer Step   = Degree - deg;
244   
245   Handle(TColgp_HArray1OfPnt2d) npoles = new
246     TColgp_HArray1OfPnt2d(1,poles->Length() + Step * (ToK2-FromK1));
247
248   Standard_Integer nbknots = BSplCLib::IncreaseDegreeCountKnots
249     (deg,Degree,periodic,mults->Array1());
250
251   Handle(TColStd_HArray1OfReal) nknots = 
252     new TColStd_HArray1OfReal(1,nbknots);
253
254   Handle(TColStd_HArray1OfInteger) nmults = 
255     new TColStd_HArray1OfInteger(1,nbknots);
256   
257   Handle(TColStd_HArray1OfReal) nweights;
258   
259   if (IsRational()) {
260     nweights = new TColStd_HArray1OfReal(1,npoles->Upper());
261   }
262
263   BSplCLib::IncreaseDegree (deg, Degree, periodic,
264                             poles->Array1(), !nweights.IsNull() ? &weights->Array1() : BSplCLib::NoWeights(),
265                             knots->Array1(), mults->Array1(),
266                             npoles->ChangeArray1(), !nweights.IsNull() ? &nweights->ChangeArray1() : BSplCLib::NoWeights(),
267                             nknots->ChangeArray1(), nmults->ChangeArray1());
268
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 Geom2d_BSplineCurve::IncreaseMultiplicity
283 (const Standard_Integer Index,
284  const Standard_Integer M)
285 {
286   TColStd_Array1OfReal k(1,1);
287   k(1) = knots->Value(Index);
288   TColStd_Array1OfInteger m(1,1);
289   m(1) = M - mults->Value(Index);
290   InsertKnots(k,m,Epsilon(1.),Standard_True);
291 }
292
293 //=======================================================================
294 //function : IncreaseMultiplicity
295 //purpose  : 
296 //=======================================================================
297
298 void Geom2d_BSplineCurve::IncreaseMultiplicity
299 (const Standard_Integer I1,
300  const Standard_Integer I2,
301  const Standard_Integer M)
302 {
303   Handle(TColStd_HArray1OfReal) tk = knots;
304   TColStd_Array1OfReal k((knots->Array1())(I1),I1,I2);
305   TColStd_Array1OfInteger m(I1,I2);
306   Standard_Integer i;
307   for (i = I1; i <= I2; i++)
308     m(i) = M - mults->Value(i);
309   InsertKnots(k,m,Epsilon(1.),Standard_True);
310 }
311
312 //=======================================================================
313 //function : IncrementMultiplicity
314 //purpose  : 
315 //=======================================================================
316
317 void Geom2d_BSplineCurve::IncrementMultiplicity
318 (const Standard_Integer I1,
319  const Standard_Integer I2,
320  const Standard_Integer Step)
321 {
322   Handle(TColStd_HArray1OfReal) tk = knots;
323   TColStd_Array1OfReal    k((knots->Array1())(I1),I1,I2);
324   TColStd_Array1OfInteger m(I1,I2);
325   m.Init(Step);
326   InsertKnots(k,m,Epsilon(1.),Standard_True);
327 }
328
329 //=======================================================================
330 //function : InsertKnot
331 //purpose  : 
332 //=======================================================================
333
334 void Geom2d_BSplineCurve::InsertKnot
335 (const Standard_Real U, 
336  const Standard_Integer M, 
337  const Standard_Real ParametricTolerance)
338 {
339   TColStd_Array1OfReal k(1,1);
340   k(1) = U;
341   TColStd_Array1OfInteger m(1,1);
342   m(1) = M;
343   InsertKnots(k,m,ParametricTolerance);
344 }
345
346 //=======================================================================
347 //function : InsertKnots
348 //purpose  : 
349 //=======================================================================
350 void  Geom2d_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("Geom2d_BSplineCurve::InsertKnots");
362
363   if (nbpoles == poles->Length()) return;
364
365   Handle(TColgp_HArray1OfPnt2d)      npoles = new TColgp_HArray1OfPnt2d(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  Geom2d_BSplineCurve::RemoveKnot
400 (const Standard_Integer Index,
401  const Standard_Integer M, 
402  const Standard_Real Tolerance)
403 {
404   if (M < 0) return Standard_True;
405
406   Standard_Integer I1  = FirstUKnotIndex ();
407   Standard_Integer I2  = LastUKnotIndex  ();
408
409   if (Index < I1 || Index > I2)  {
410     throw Standard_OutOfRange("BSpline curve : RemoveKnot : index out of range");
411   }
412   
413   const TColgp_Array1OfPnt2d & oldpoles   = poles->Array1();
414   
415   Standard_Integer step = mults->Value(Index) - M;
416   if (step <= 0) return Standard_True;
417   
418   Handle(TColgp_HArray1OfPnt2d) npoles =
419     new TColgp_HArray1OfPnt2d(1,oldpoles.Length()-step);
420   
421   Handle(TColStd_HArray1OfReal)    nknots  = knots;
422   Handle(TColStd_HArray1OfInteger) nmults  = mults;
423   
424   if (M == 0) {
425     nknots = new TColStd_HArray1OfReal(1,knots->Length()-1);
426     nmults = new TColStd_HArray1OfInteger(1,knots->Length()-1);
427   }
428   
429   Handle(TColStd_HArray1OfReal) nweights;
430   if (IsRational())
431   {
432     nweights = new TColStd_HArray1OfReal(1,npoles->Length());
433   }
434
435   if (!BSplCLib::RemoveKnot (Index, M, deg, periodic,
436                              poles->Array1(), !nweights.IsNull() ? &weights->Array1() : BSplCLib::NoWeights(),
437                              knots->Array1(),mults->Array1(),
438                              npoles->ChangeArray1(), !nweights.IsNull() ? &nweights->ChangeArray1() : BSplCLib::NoWeights(),
439                              nknots->ChangeArray1(),nmults->ChangeArray1(),
440                              Tolerance))
441   {
442     return Standard_False;
443   }
444
445   weights = nweights;
446   poles = npoles;
447   knots = nknots;
448   mults = nmults;
449   
450   UpdateKnots();
451   maxderivinvok = 0;
452   return Standard_True;
453 }
454
455 //=======================================================================
456 //function : InsertPoleAfter
457 //purpose  : 
458 //=======================================================================
459
460 void Geom2d_BSplineCurve::InsertPoleAfter
461 (const Standard_Integer Index,
462  const gp_Pnt2d& P,
463  const Standard_Real Weight)
464 {
465   if (Index < 0 || Index > poles->Length())  throw Standard_OutOfRange("BSpline curve : InsertPoleAfter: Index and #pole mismatch");
466
467   if (Weight <= gp::Resolution())     throw Standard_ConstructionError("BSpline curve : InsertPoleAfter: Weight too small");
468
469   if (knotSet == GeomAbs_NonUniform || knotSet == GeomAbs_PiecewiseBezier) {
470     throw Standard_ConstructionError("BSpline curve : InsertPoleAfter : bad knotSet type");
471   }
472
473   const TColStd_Array1OfReal& cknots   = knots->Array1();
474   Standard_Integer nbknots = cknots.Length();
475
476   Handle(TColStd_HArray1OfReal) nknots =  
477     new TColStd_HArray1OfReal(1,nbknots+1);
478
479   TColStd_Array1OfReal& newknots = nknots->ChangeArray1();
480
481   Standard_Integer i;
482   for (i = 1; i < nbknots; i++) {
483     newknots (i) = cknots(i);
484   }
485   
486   newknots (nbknots+1) = 2 * newknots (nbknots) - newknots(nbknots-1);
487   
488   Handle(TColStd_HArray1OfInteger) nmults =
489     new TColStd_HArray1OfInteger(1,nbknots+1);
490
491   TColStd_Array1OfInteger& newmults     = nmults->ChangeArray1();
492   const TColStd_Array1OfInteger& cmults = mults->Array1();
493
494   for (i = 2; i <= nbknots; i++) newmults (i) = 1;
495   newmults (1)         =  cmults(1);
496   newmults (nbknots+1) =  cmults(nbknots+1);
497   
498   const TColgp_Array1OfPnt2d& cpoles = poles->Array1();
499   Standard_Integer nbpoles = cpoles.Length();
500   Handle(TColgp_HArray1OfPnt2d) npoles = 
501     new TColgp_HArray1OfPnt2d(1, nbpoles+1);
502   TColgp_Array1OfPnt2d& newpoles = npoles->ChangeArray1();
503
504   // insert the pole
505
506   for (i = 1; i <= Index; i++)
507     newpoles(i) = cpoles(i);
508
509   newpoles(Index+1) = P;
510
511   for (i = Index+1; i <= nbpoles; i++)
512     newpoles(i+1) = cpoles(i);
513
514   // Insert the weight
515
516   Handle(TColStd_HArray1OfReal) nweights;
517   Standard_Boolean rat = IsRational() || Abs(Weight-1.) > gp::Resolution();
518   
519   if (rat) {
520     nweights = new TColStd_HArray1OfReal(1,nbpoles+1);
521     TColStd_Array1OfReal& newweights = nweights->ChangeArray1();
522
523     for (i = 1; i <= Index; i++)
524       if (IsRational())
525         newweights(i) = weights->Value(i);
526       else
527         newweights(i) = 1.;
528     
529     newweights(Index+1) = Weight;
530
531     for (i = Index+1; i <= nbpoles; i++)
532       if (IsRational())
533         newweights(i+1) = weights->Value(i);
534       else
535         newweights(i+1) = 1.;
536   }
537   
538   poles   = npoles;
539   weights = nweights;
540   knots   = nknots;
541   mults   = nmults;
542   maxderivinvok = 0;
543   UpdateKnots();
544 }
545
546 //=======================================================================
547 //function : InsertPoleBefore
548 //purpose  : 
549 //=======================================================================
550
551 void Geom2d_BSplineCurve::InsertPoleBefore
552 (const Standard_Integer Index,
553  const gp_Pnt2d& P,
554  const Standard_Real Weight)
555 {
556   InsertPoleAfter(Index-1,P,Weight);
557 }
558
559 //=======================================================================
560 //function : RemovePole
561 //purpose  : 
562 //=======================================================================
563
564 void Geom2d_BSplineCurve::RemovePole
565 (const Standard_Integer Index)
566 {
567   if (Index < 1 || Index > poles->Length())  throw Standard_OutOfRange("BSpline curve :RemovePole : Index and #pole mismatch");
568
569   if (poles->Length() <= 2)           throw Standard_ConstructionError("BSpline curve : RemovePole : #pole is already minimum");
570
571   if (knotSet == GeomAbs_NonUniform || knotSet == GeomAbs_PiecewiseBezier) 
572     throw Standard_ConstructionError("BSpline curve : RemovePole: bad knotSet type");
573
574   Standard_Integer i;
575   Handle(TColStd_HArray1OfReal) nknots =
576     new TColStd_HArray1OfReal(1,knots->Length()-1);
577   TColStd_Array1OfReal& newknots = nknots->ChangeArray1();
578
579   Handle(TColStd_HArray1OfInteger) nmults =
580     new TColStd_HArray1OfInteger(1,mults->Length()-1);
581   TColStd_Array1OfInteger& newmults = nmults->ChangeArray1();
582
583   for (i = 1; i < newknots.Length(); i++) {
584     newknots (i) = knots->Value (i);
585     newmults (i) = 1;
586   }
587   newmults(1) = mults->Value(1);
588   newknots(newknots.Upper()) = knots->Value (knots->Upper());
589   newmults(newmults.Upper()) = mults->Value (mults->Upper());
590
591
592   Handle(TColgp_HArray1OfPnt2d) npoles =
593     new TColgp_HArray1OfPnt2d(1, poles->Upper()-1);
594   TColgp_Array1OfPnt2d& newpoles = npoles->ChangeArray1();
595
596   for (i = 1; i < Index; i++)
597     newpoles(i) = poles->Value(i);
598   for (i = Index; i < newpoles.Length(); i++)
599     newpoles(i) = poles->Value(i+1);
600
601   Handle(TColStd_HArray1OfReal) nweights;
602   if (IsRational()) {
603     nweights = new TColStd_HArray1OfReal(1,newpoles.Length());
604     TColStd_Array1OfReal& newweights = nweights->ChangeArray1();
605     for (i = 1; i < Index; i++)
606       newweights(i) = weights->Value(i);
607     for (i = Index; i < newweights.Length(); i++)
608       newweights(i) = weights->Value(i+1);
609   }
610
611   poles   = npoles;
612   weights = nweights;
613   knots   = nknots;
614   mults   = nmults;
615   UpdateKnots();
616 }
617
618 //=======================================================================
619 //function : Reverse
620 //purpose  : 
621 //=======================================================================
622
623 void Geom2d_BSplineCurve::Reverse ()
624
625   BSplCLib::Reverse(knots->ChangeArray1());
626   BSplCLib::Reverse(mults->ChangeArray1());
627   Standard_Integer last;
628   if (periodic)
629     last = flatknots->Upper() - deg - 1;
630   else
631     last = poles->Upper();
632   BSplCLib::Reverse(poles->ChangeArray1(),last);
633   if (rational)
634     BSplCLib::Reverse(weights->ChangeArray1(),last);
635   UpdateKnots();
636 }
637
638 //=======================================================================
639 //function : ReversedParameter
640 //purpose  : 
641 //=======================================================================
642
643 Standard_Real Geom2d_BSplineCurve::ReversedParameter( const Standard_Real U) const
644 {
645   return (FirstParameter() + LastParameter() - U);
646 }
647
648 //=======================================================================
649 //function : Segment
650 //purpose  : 
651 //=======================================================================
652 void Geom2d_BSplineCurve::Segment(const Standard_Real aU1,
653                                   const Standard_Real aU2,
654                                   const Standard_Real theTolerance)
655 {
656   if (aU2 < aU1)
657     throw Standard_DomainError("Geom2d_BSplineCurve::Segment");
658   //
659   Standard_Real AbsUMax = Max(Abs(FirstParameter()),Abs(LastParameter()));
660   Standard_Real Eps = Max (Epsilon(AbsUMax), theTolerance);
661   Standard_Real NewU1, NewU2;
662   Standard_Real U, DU=0;
663   Standard_Integer i, k, index;
664   //
665   //f
666   // Checking the input bounds aUj (j=1,2). 
667   // For the case when aUj==knot(i), 
668   // in order to prevent the insertion of a new knot that will be too closed 
669   // to the existing knot,  
670   // we assign Uj=knot(i)
671   Standard_Integer n1, n2;
672   Standard_Real U1, U2;
673   //
674   U1=aU1;
675   U2=aU2;
676   n1=knots->Lower();
677   n2=knots->Upper();
678   for (i=n1; i<=n2; ++i) {
679     U=knots->Value(i);
680     if (Abs(U-aU1)<=Eps) {
681       U1=U;
682     }
683     else if (Abs(U-aU2)<=Eps) {
684       U2=U;
685     }
686   }
687   // Henceforward we use U1, U2 as bounds of the segment
688   //t
689   // 
690   TColStd_Array1OfReal    Knots(1,2);
691   TColStd_Array1OfInteger Mults(1,2);
692   //
693   // define param ditance to keep (eap, Apr 18 2002, occ311)
694   if (periodic) {
695     Standard_Real Period = LastParameter() - FirstParameter();
696     DU = U2 - U1;
697     if (DU - Period > Precision::PConfusion())
698       throw Standard_DomainError("Geom2d_BSplineCurve::Segment");
699     if (DU > Period)
700       DU = Period;
701   }
702   //
703   index = 0;
704   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
705                             U1,periodic,knots->Lower(),knots->Upper(),
706                             index,NewU1);
707   index = 0;
708   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
709                             U2,periodic,knots->Lower(),knots->Upper(),
710                             index,NewU2);
711   Knots(1) = Min( NewU1, NewU2);
712   Knots(2) = Max( NewU1, NewU2);
713   Mults(1) = Mults( 2) = deg;
714   InsertKnots(Knots, Mults, Eps);
715   
716   if (periodic) { // set the origine at NewU1
717     index = 0;
718     BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
719                               U1,periodic,knots->Lower(),knots->Upper(),
720                               index,U);
721     // Eps = Epsilon(knots->Value(index+1));
722     if ( Abs(knots->Value(index+1)-U) <= Eps) {
723       index++;
724     }
725     SetOrigin(index);
726     SetNotPeriodic();
727     NewU2 = NewU1 + DU;
728   }
729
730   // compute index1 and index2 to set the new knots and mults 
731   Standard_Integer index1 = 0, index2 = 0;
732   Standard_Integer FromU1 = knots->Lower();
733   Standard_Integer ToU2   = knots->Upper();
734   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
735                             NewU1,periodic,FromU1,ToU2,index1,U);
736   BSplCLib::LocateParameter(deg,knots->Array1(),mults->Array1(),
737                             NewU2,periodic,FromU1,ToU2,index2,U);
738   // Eps = Epsilon(knots->Value(index2+1));
739   if ( Abs(knots->Value(index2+1)-U) <= Eps){
740     index2++;
741   }
742   
743   Standard_Integer nbknots = index2 - index1 + 1;
744
745   Handle(TColStd_HArray1OfReal) 
746     nknots = new TColStd_HArray1OfReal(1,nbknots);
747   Handle(TColStd_HArray1OfInteger) 
748     nmults = new TColStd_HArray1OfInteger(1,nbknots);
749
750   // to restore changed U1
751   if (DU > 0) {// if was periodic
752     DU = NewU1 - U1;
753   }
754   //
755   k = 1;
756   //
757   for ( i = index1; i<= index2; i++) {
758     nknots->SetValue(k, knots->Value(i) - DU);
759     nmults->SetValue(k, mults->Value(i));
760     k++;
761   }
762   nmults->SetValue(      1, deg + 1);
763   nmults->SetValue(nbknots, deg + 1);
764
765
766   // compute index1 and index2 to set the new poles and weights
767   Standard_Integer pindex1 
768     = BSplCLib::PoleIndex(deg,index1,periodic,mults->Array1());
769   Standard_Integer pindex2 
770     = BSplCLib::PoleIndex(deg,index2,periodic,mults->Array1());
771
772   pindex1++;
773   pindex2 = Min( pindex2+1, poles->Length());
774
775   Standard_Integer nbpoles  = pindex2 - pindex1 + 1;
776
777   Handle(TColStd_HArray1OfReal) 
778     nweights = new TColStd_HArray1OfReal(1,nbpoles);
779   Handle(TColgp_HArray1OfPnt2d)
780     npoles = new TColgp_HArray1OfPnt2d(1,nbpoles);
781
782   k = 1;
783   if ( rational) {
784     nweights = new TColStd_HArray1OfReal( 1, nbpoles);
785     for ( i = pindex1; i <= pindex2; i++) {
786       npoles->SetValue(k, poles->Value(i));
787       nweights->SetValue(k, weights->Value(i));
788       k++;
789     }
790   }
791   else {
792     for ( i = pindex1; i <= pindex2; i++) {
793       npoles->SetValue(k, poles->Value(i));
794       k++;
795     }
796   }
797
798   knots = nknots;
799   mults = nmults;
800   poles = npoles;
801   if (rational){ 
802     weights = nweights;
803   }
804   UpdateKnots();
805 }
806
807 //=======================================================================
808 //function : SetKnot
809 //purpose  : 
810 //=======================================================================
811
812 void Geom2d_BSplineCurve::SetKnot
813 (const Standard_Integer Index,
814  const Standard_Real K)
815 {
816   if (Index < 1 || Index > knots->Length())     throw Standard_OutOfRange("BSpline curve : SetKnot:  Index and #pole mismatch");
817   Standard_Real DK = Abs(Epsilon (K));
818   if (Index == 1) { 
819     if (K >= knots->Value(2) - DK) throw Standard_ConstructionError("BSpline curve :SetKnot :K out of range");
820   }
821   else if (Index == knots->Length()) {
822     if (K <= knots->Value (knots->Length()-1) + DK)  {
823       throw Standard_ConstructionError("BSpline curve : SetKnot : K out of range");
824     }
825   }
826   else {
827     if (K <= knots->Value(Index-1) + DK ||
828         K >= knots->Value(Index+1) - DK ) {
829       throw Standard_ConstructionError("BSpline curve : SetKnot: K out of range");
830     }
831   }
832   if (K != knots->Value (Index)) {
833     knots->SetValue (Index, K);
834     maxderivinvok = 0;
835     UpdateKnots();
836   }
837 }
838
839 //=======================================================================
840 //function : SetKnots
841 //purpose  : 
842 //=======================================================================
843
844 void Geom2d_BSplineCurve::SetKnots
845 (const TColStd_Array1OfReal& K)
846 {
847   CheckCurveData(poles->Array1(),K,mults->Array1(),deg,periodic);
848   knots->ChangeArray1() = K;
849   maxderivinvok = 0;
850   UpdateKnots();
851 }
852
853 //=======================================================================
854 //function : SetKnot
855 //purpose  : 
856 //=======================================================================
857
858 void Geom2d_BSplineCurve::SetKnot
859 (const Standard_Integer Index,
860  const Standard_Real K,
861  const Standard_Integer M)
862 {
863   IncreaseMultiplicity (Index, M);
864   SetKnot (Index, K);
865 }
866
867 //=======================================================================
868 //function : SetPeriodic
869 //purpose  : 
870 //=======================================================================
871
872 void Geom2d_BSplineCurve::SetPeriodic ()
873 {
874   Standard_Integer first = FirstUKnotIndex();
875   Standard_Integer last  = LastUKnotIndex();
876
877   Handle(TColStd_HArray1OfReal) tk = knots;
878   TColStd_Array1OfReal    cknots((knots->Array1())(first),first,last);
879   knots = new TColStd_HArray1OfReal(1,cknots.Length());
880   knots->ChangeArray1() = cknots;
881
882   Handle(TColStd_HArray1OfInteger) tm = mults;
883   TColStd_Array1OfInteger cmults((mults->Array1())(first),first,last);
884   cmults(first) = cmults(last) = Min(deg, Max( cmults(first), cmults(last)));
885   mults = new TColStd_HArray1OfInteger(1,cmults.Length());
886   mults->ChangeArray1() = cmults;
887
888   // compute new number of poles;
889   Standard_Integer nbp = BSplCLib::NbPoles(deg,Standard_True,cmults);
890   
891   Handle(TColgp_HArray1OfPnt2d) tp = poles;
892   TColgp_Array1OfPnt2d cpoles((poles->Array1())(1),1,nbp);
893   poles = new TColgp_HArray1OfPnt2d(1,nbp);
894   poles->ChangeArray1() = cpoles;
895   
896   if (rational) {
897     Handle(TColStd_HArray1OfReal) tw = weights;
898     TColStd_Array1OfReal cweights((weights->Array1())(1),1,nbp);
899     weights = new TColStd_HArray1OfReal(1,nbp);
900     weights->ChangeArray1() = cweights;
901   }
902
903   periodic = Standard_True;
904   
905   maxderivinvok = 0;
906   UpdateKnots();
907 }
908
909 //=======================================================================
910 //function : SetOrigin
911 //purpose  : 
912 //=======================================================================
913
914 void Geom2d_BSplineCurve::SetOrigin(const Standard_Integer Index)
915 {
916   if (!periodic)
917     throw Standard_NoSuchObject("Geom2d_BSplineCurve::SetOrigin");
918   Standard_Integer i,k;
919   Standard_Integer first = FirstUKnotIndex();
920   Standard_Integer last  = LastUKnotIndex();
921
922   if ((Index < first) || (Index > last))
923   throw Standard_DomainError("Geom2d_BSplineCurve::SetOrigin");
924
925   Standard_Integer nbknots = knots->Length();
926   Standard_Integer nbpoles = poles->Length();
927
928   Handle(TColStd_HArray1OfReal) nknots = 
929     new TColStd_HArray1OfReal(1,nbknots);
930   TColStd_Array1OfReal& newknots = nknots->ChangeArray1();
931
932   Handle(TColStd_HArray1OfInteger) nmults =
933     new TColStd_HArray1OfInteger(1,nbknots);
934   TColStd_Array1OfInteger& newmults = nmults->ChangeArray1();
935
936   // set the knots and mults
937   Standard_Real period = knots->Value(last) - knots->Value(first);
938   k = 1;
939   for ( i = Index; i <= last ; i++) {
940     newknots(k) = knots->Value(i);
941     newmults(k) = mults->Value(i);
942     k++;
943   }
944   for ( i = first+1; i <= Index; i++) {
945     newknots(k) = knots->Value(i) + period;
946     newmults(k) = mults->Value(i);
947     k++;
948   }
949
950   Standard_Integer index = 1;
951   for (i = first+1; i <= Index; i++) 
952     index += mults->Value(i);
953
954   // set the poles and weights
955   Handle(TColgp_HArray1OfPnt2d) npoles =
956     new TColgp_HArray1OfPnt2d(1,nbpoles);
957   Handle(TColStd_HArray1OfReal) nweights =
958     new TColStd_HArray1OfReal(1,nbpoles);
959   TColgp_Array1OfPnt2d & newpoles   = npoles->ChangeArray1();
960   TColStd_Array1OfReal & newweights = nweights->ChangeArray1();
961   first = poles->Lower();
962   last  = poles->Upper();
963   if ( rational) {
964     k = 1;
965     for ( i = index; i <= last; i++) {
966       newpoles(k)   = poles->Value(i);
967       newweights(k) = weights->Value(i);
968       k++;
969     }
970     for ( i = first; i < index; i++) {
971       newpoles(k)   = poles->Value(i);
972       newweights(k) = weights->Value(i);
973       k++;
974     }
975   }
976   else {
977     k = 1;
978     for ( i = index; i <= last; i++) {
979       newpoles(k) = poles->Value(i);
980       k++;
981     }
982     for ( i = first; i < index; i++) {
983       newpoles(k) = poles->Value(i);
984       k++;
985     }
986   }
987
988   poles = npoles;
989   knots = nknots;
990   mults = nmults;
991   if (rational) 
992     weights = nweights;
993   maxderivinvok = 0;
994   UpdateKnots();
995 }
996
997 //=======================================================================
998 //function : SetNotPeriodic
999 //purpose  : 
1000 //=======================================================================
1001
1002 void Geom2d_BSplineCurve::SetNotPeriodic () 
1003
1004   if (periodic) {
1005     Standard_Integer NbKnots, NbPoles;
1006     BSplCLib::PrepareUnperiodize( deg, mults->Array1(),NbKnots,NbPoles);
1007     
1008     Handle(TColgp_HArray1OfPnt2d) npoles 
1009       = new TColgp_HArray1OfPnt2d(1,NbPoles);
1010     
1011     Handle(TColStd_HArray1OfReal) nknots 
1012       = new TColStd_HArray1OfReal(1,NbKnots);
1013     
1014     Handle(TColStd_HArray1OfInteger) nmults
1015       = new TColStd_HArray1OfInteger(1,NbKnots);
1016     
1017     Handle(TColStd_HArray1OfReal) nweights;
1018     if (IsRational()) {
1019       nweights = new TColStd_HArray1OfReal(1,NbPoles);
1020     }
1021
1022     BSplCLib::Unperiodize (deg,
1023                            mults->Array1(), knots->Array1(), poles->Array1(),
1024                            !nweights.IsNull() ? &weights->Array1() : BSplCLib::NoWeights(),
1025                            nmults->ChangeArray1(), nknots->ChangeArray1(), npoles->ChangeArray1(),
1026                            !nweights.IsNull() ? &nweights->ChangeArray1() : BSplCLib::NoWeights());
1027     poles   = npoles;
1028     weights = nweights;
1029     mults   = nmults;
1030     knots   = nknots;
1031     periodic = Standard_False;
1032     maxderivinvok = 0;
1033     UpdateKnots();
1034   }
1035 }
1036
1037 //=======================================================================
1038 //function : SetPole
1039 //purpose  : 
1040 //=======================================================================
1041
1042 void Geom2d_BSplineCurve::SetPole
1043 (const Standard_Integer Index,
1044  const gp_Pnt2d& P)
1045 {
1046   if (Index < 1 || Index > poles->Length()) throw Standard_OutOfRange("BSpline curve : SetPole : index and #pole mismatch");
1047   poles->SetValue (Index, P);
1048   maxderivinvok = 0;
1049 }
1050
1051 //=======================================================================
1052 //function : SetPole
1053 //purpose  : 
1054 //=======================================================================
1055
1056 void Geom2d_BSplineCurve::SetPole
1057 (const Standard_Integer Index,
1058  const gp_Pnt2d& P,
1059  const Standard_Real W)
1060 {
1061   SetPole(Index,P);
1062   SetWeight(Index,W);
1063 }
1064
1065 //=======================================================================
1066 //function : SetWeight
1067 //purpose  : 
1068 //=======================================================================
1069
1070 void Geom2d_BSplineCurve::SetWeight
1071 (const Standard_Integer Index,
1072  const Standard_Real W)
1073 {
1074   if (Index < 1 || Index > poles->Length())   throw Standard_OutOfRange("BSpline curve : SetWeight: Index and #pole mismatch");
1075
1076   if (W <= gp::Resolution ())     throw Standard_ConstructionError("BSpline curve : SetWeight: Weight too small");
1077
1078
1079   Standard_Boolean rat = IsRational() || (Abs(W - 1.) > gp::Resolution());
1080
1081   if ( rat) { 
1082     if (rat && !IsRational()) {
1083       weights = new TColStd_HArray1OfReal(1,poles->Length());
1084       weights->Init(1.);
1085     }
1086     
1087     weights->SetValue (Index, W);
1088     
1089     if (IsRational()) {
1090       rat = Rational(weights->Array1());
1091       if (!rat) weights.Nullify();
1092     }
1093     
1094     rational = !weights.IsNull();
1095   }
1096   
1097   maxderivinvok = 0;
1098 }
1099
1100 //=======================================================================
1101 //function : MovePoint
1102 //purpose  : 
1103 //=======================================================================
1104
1105 void Geom2d_BSplineCurve::MovePoint(const Standard_Real U,
1106                                     const gp_Pnt2d& P,
1107                                     const Standard_Integer Index1,
1108                                     const Standard_Integer Index2,
1109                                     Standard_Integer& FirstModifiedPole,
1110                                     Standard_Integer& LastmodifiedPole)
1111 {
1112   if (Index1 < 1 || Index1 > poles->Length() || 
1113       Index2 < 1 || Index2 > poles->Length() || Index1 > Index2) {
1114     throw Standard_OutOfRange("BSpline curve :  MovePoint: Index and #pole mismatch");
1115   }
1116   TColgp_Array1OfPnt2d npoles(1, poles->Length());
1117   gp_Pnt2d P0;
1118   D0(U, P0);
1119   gp_Vec2d Displ(P0, P);
1120
1121   BSplCLib::MovePoint (U, Displ, Index1, Index2, deg, poles->Array1(),
1122                        rational ? &weights->Array1() : BSplCLib::NoWeights(), flatknots->Array1(),
1123                        FirstModifiedPole, LastmodifiedPole, npoles);
1124   if (FirstModifiedPole) {
1125     poles->ChangeArray1() = npoles;
1126     maxderivinvok = 0;
1127   }
1128 }
1129
1130 //=======================================================================
1131 //function : MovePointAndTangent
1132 //purpose  : 
1133 //=======================================================================
1134
1135 void Geom2d_BSplineCurve::
1136 MovePointAndTangent(const Standard_Real    U,
1137                     const gp_Pnt2d&        P,
1138                     const gp_Vec2d&        Tangent,
1139                     const Standard_Real    Tolerance,
1140                     const Standard_Integer StartingCondition,
1141                     const Standard_Integer EndingCondition,
1142                     Standard_Integer&      ErrorStatus) 
1143 {
1144   Standard_Integer ii ;
1145   if (IsPeriodic()) {
1146     //
1147     // for the time being do not deal with periodic curves
1148     //
1149     SetNotPeriodic() ;
1150   }
1151   TColgp_Array1OfPnt2d new_poles(1, poles->Length());
1152   gp_Pnt2d P0;
1153   
1154
1155   gp_Vec2d delta_derivative;
1156   D1(U, P0,
1157      delta_derivative) ;
1158   gp_Vec2d delta(P0, P);
1159   for (ii = 1 ; ii <= 2 ; ii++) {
1160     delta_derivative.SetCoord(ii, 
1161                               Tangent.Coord(ii)- delta_derivative.Coord(ii)) ;
1162   }
1163
1164   BSplCLib::MovePointAndTangent (U,
1165                                  delta,
1166                                  delta_derivative,
1167                                  Tolerance,
1168                                  deg,
1169                                  StartingCondition,
1170                                  EndingCondition,
1171                                  poles->Array1(),
1172                                  rational ? &weights->Array1() : BSplCLib::NoWeights(),
1173                                  flatknots->Array1(),
1174                                  new_poles,
1175                                  ErrorStatus);
1176   if (!ErrorStatus) {
1177     poles->ChangeArray1() = new_poles;
1178     maxderivinvok = 0;
1179   }
1180 }
1181
1182 //=======================================================================
1183 //function : UpdateKnots
1184 //purpose  : 
1185 //=======================================================================
1186
1187 void Geom2d_BSplineCurve::UpdateKnots()
1188 {
1189
1190   rational = !weights.IsNull();
1191
1192   Standard_Integer MaxKnotMult = 0;
1193   BSplCLib::KnotAnalysis (deg,
1194                 periodic,
1195                 knots->Array1(), 
1196                 mults->Array1(), 
1197                 knotSet, MaxKnotMult);
1198   
1199   if (knotSet == GeomAbs_Uniform && !periodic)  {
1200     flatknots = knots;
1201   }
1202   else {
1203     flatknots = new TColStd_HArray1OfReal 
1204       (1, BSplCLib::KnotSequenceLength(mults->Array1(),deg,periodic));
1205
1206     BSplCLib::KnotSequence (knots->Array1(), 
1207                             mults->Array1(),
1208                             deg,periodic,
1209                             flatknots->ChangeArray1());
1210   }
1211   
1212   if (MaxKnotMult == 0)  smooth = GeomAbs_CN;
1213   else {
1214     switch (deg - MaxKnotMult) {
1215     case 0 :   smooth = GeomAbs_C0;   break;
1216     case 1 :   smooth = GeomAbs_C1;   break;
1217     case 2 :   smooth = GeomAbs_C2;   break;
1218     case 3 :   smooth = GeomAbs_C3;   break;
1219       default :  smooth = GeomAbs_C3;   break;
1220     }
1221   }
1222 }
1223
1224 //=======================================================================
1225 //function : Normalizes the parameters if the curve is periodic
1226 //purpose  : that is compute the cache so that it is valid
1227 //=======================================================================
1228
1229 void Geom2d_BSplineCurve::PeriodicNormalization(Standard_Real&  Parameter) const 
1230 {
1231   Standard_Real Period ;
1232
1233   if (periodic) {
1234     Period = flatknots->Value(flatknots->Upper() - deg) - flatknots->Value (deg + 1) ;
1235     while (Parameter > flatknots->Value(flatknots->Upper()-deg)) {
1236       Parameter -= Period ;
1237     }
1238     while (Parameter < flatknots->Value((deg + 1))) {
1239       Parameter +=  Period ;
1240     }
1241   }
1242 }
1243