projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0024428: Implementation of LGPL license
[occt.git] / src / Hermit / Hermit.cxx
CommitLineData
b311480e 1// Created on: 1997-01-15
2// Created by: Stagiaire Francois DUMONT
3// Copyright (c) 1997-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
973c2be1 8// This library is free software; you can redistribute it and / or modify it
9// under the terms of the GNU Lesser General Public 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.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
17#include <Hermit.ixx>
18#include <Geom_BSplineCurve.hxx>
19#include <Geom2d_BSplineCurve.hxx>
20#include <BSplCLib.hxx>
21#include <TColStd_Array1OfReal.hxx>
22#include <TColStd_HArray1OfReal.hxx>
23#include <TColStd_Array1OfInteger.hxx>
24#include <TColStd_HArray1OfInteger.hxx>
25#include <PLib.hxx>
26#include <Standard_Boolean.hxx>
27#include <gp_Pnt2d.hxx>
28#include <TColgp_Array1OfPnt2d.hxx>
29#include <Standard_DimensionError.hxx>
30#include <Standard_Real.hxx>
31#include <TCollection_CompareOfReal.hxx>
32#include <SortTools_QuickSortOfReal.hxx>
33#include <Precision.hxx>
34
35//=======================================================================
36//function : HermiteCoeff
37//purpose : calculate the Hermite coefficients of degree 3 from BS and
38// store them in TAB(4 coefficients)
39//=======================================================================
40
41static void HermiteCoeff(const Handle(Geom_BSplineCurve)& BS,
42 TColStd_Array1OfReal& TAB)
43
44{
45 TColStd_Array1OfReal Knots(1,BS->NbKnots());
46 TColStd_Array1OfReal Weights(1,BS->NbPoles());
47 TColStd_Array1OfInteger Mults(1,BS->NbKnots());
48 Standard_Integer Degree,Index0,Index1; // denominateur value for u=0 & u=1
49 Standard_Real Denom0,Denom1, // denominator value for u=0 & u=1
50 Deriv0,Deriv1 ; // derivative denominator value for u=0 & 1
51 Standard_Boolean Periodic;
52
53 BS->Knots(Knots);
54 BSplCLib::Reparametrize(0.0,1.0,Knots); //affinity on the nodal vector
55 BS->Weights(Weights);
56 BS->Multiplicities(Mults);
57 Degree = BS->Degree();
58 Periodic = BS->IsPeriodic();
59 Index0 = BS->FirstUKnotIndex();
60 Index1 = BS->LastUKnotIndex()-1;
61
62 BSplCLib::D1(0.0,Index0,Degree,Periodic,Weights,BSplCLib::NoWeights(),Knots,Mults,Denom0,Deriv0);
63 BSplCLib::D1(1.0,Index1,Degree,Periodic,Weights,BSplCLib::NoWeights(),Knots,Mults,Denom1,Deriv1);
64 TAB(0) = 1/Denom0; //Hermit coefficients
65 TAB(1) = -Deriv0/(Denom0*Denom0);
66 TAB(2) = -Deriv1/(Denom1*Denom1);
67 TAB(3) = 1/Denom1;
68
69}
70
71//=======================================================================
72//function : HermiteCoeff
73//purpose : calculate the Hermite coefficients of degree 3 from BS and
74// store them in TAB(4 coefficients)
75//=======================================================================
76
77static void HermiteCoeff(const Handle(Geom2d_BSplineCurve)& BS,
78 TColStd_Array1OfReal& TAB)
79
80{
81 TColStd_Array1OfReal Knots(1,BS->NbKnots());
82 TColStd_Array1OfReal Weights(1,BS->NbPoles());
83 TColStd_Array1OfInteger Mults(1,BS->NbKnots());
84 Standard_Integer Degree,Index0,Index1;
85 Standard_Real Denom0,Denom1, // denominateur value for u=0 & u=1
86 Deriv0,Deriv1 ; // denominator value for u=0 & u=1
87 Standard_Boolean Periodic; // derivative denominatur value for u=0 & 1
88
89 BS->Knots(Knots);
90 BSplCLib::Reparametrize(0.0,1.0,Knots); //affinity on the nodal vector
91 BS->Weights(Weights);
92 BS->Multiplicities(Mults);
93 Degree = BS->Degree();
94 Periodic = BS->IsPeriodic();
95 Index0 = BS->FirstUKnotIndex();
96 Index1 = BS->LastUKnotIndex()-1;
97
98 BSplCLib::D1(0.0,Index0,Degree,Periodic,Weights,BSplCLib::NoWeights(),Knots,Mults,Denom0,Deriv0);
99 BSplCLib::D1(1.0,Index1,Degree,Periodic,Weights,BSplCLib::NoWeights(),Knots,Mults,Denom1,Deriv1);
100 TAB(0) = 1/Denom0; //Hermit coefficients
101 TAB(1) = -Deriv0/(Denom0*Denom0);
102 TAB(2) = -Deriv1/(Denom1*Denom1);
103 TAB(3) = 1/Denom1;
104
105}
106
107//=======================================================================
108//function : SignDenom
109//purpose : give the sign of Herm(0) True=Positive
110//=======================================================================
111
112static Standard_Boolean SignDenom(const TColgp_Array1OfPnt2d& Poles)
113
114{
115 Standard_Boolean Result;
116
117 if (Poles(0).Y()<0)
118 Result=Standard_False;
119 else Result=Standard_True;
120 return Result;
121}
122
123//=======================================================================
124//function : Polemax
125//purpose : give the max and the min of the Poles (by their index)
126//=======================================================================
127
128
129static void Polemax(const TColgp_Array1OfPnt2d& Poles,
130 Standard_Integer& min,
131 Standard_Integer& max)
132
133{
134// Standard_Integer i,index=0;
135 Standard_Integer i;
136 Standard_Real Max,Min; //intermediate value of max and min ordinates
137 min=0;max=0; //initialisation of the indices
138
139 Min=Poles(0).Y(); //initialisation of the intermediate value
140 Max=Poles(0).Y();
141 for (i=1;i<=(Poles.Length()-1);i++){
142 if (Poles(i).Y()<Min){
143 Min=Poles(i).Y();
144 min=i;
145 }
146 if (Poles(i).Y()>Max){
147 Max=Poles(i).Y();
148 max=i;
149 }
150 }
151}
152
153//=======================================================================
154//function : PolyTest
155//purpose : give the knots U4 and U5 to insert to a(u)
156//=======================================================================
157
158
159static void PolyTest(const TColStd_Array1OfReal& Herm,
160 const Handle(Geom_BSplineCurve)& BS,
161 Standard_Real& U4,
162 Standard_Real& U5,
163 Standard_Integer& boucle,
164 const Standard_Real TolPoles,
165// const Standard_Real TolKnots,
166 const Standard_Real ,
167 const Standard_Real Ux,
168 const Standard_Real Uy)
169
170{
171 Standard_Integer index,i,
172 I1=0,I2=0,I3=0,I4=0; //knots index
173 TColgp_Array1OfPnt2d Polesinit(0,3) ;
174 Handle(TColStd_HArray1OfReal) Knots; //array of the BSpline knots + the ones inserted
175 Standard_Integer cas=0,mark=0,dercas=0, //loop marks
176 min,max; //Pole min and max indices
177 Standard_Real Us1,Us2,a; //bounderies value of the knots to be inserted
178// gp_Pnt2d P ;
179 TCollection_CompareOfReal Comp;
180
181 U4=0.0;U5=1.0; //default value
182 if (Ux!=1.0){
183 BS->LocateU(Ux,0.0,I1,I2); //localization of the inserted knots
184 if (Uy!=0.0)
185 BS->LocateU(Uy,0.0,I3,I4);
186 }
187
188 if (I1==I2) //definition and filling of the
189 if((I3==I4)||(I3==0)){ //array of knots
190 Knots=new TColStd_HArray1OfReal(1,BS->NbKnots());
191 for (i=1;i<=BS->NbKnots();i++)
192 Knots->SetValue(i,BS->Knot(i));
193 }
194 else{
195 Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()+1);
196 for (i=1;i<=BS->NbKnots();i++)
197 Knots->SetValue(i,BS->Knot(i));
198 Knots->SetValue(BS->NbKnots()+1,Uy);
199 }
200 else{
201 if((I3==I4)||(I3==0)){
202 Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()+1);
203 for (i=1;i<=BS->NbKnots();i++)
204 Knots->SetValue(i,BS->Knot(i));
205 Knots->SetValue(BS->NbKnots()+1,Ux);
206 }
207 else{
208 Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()+2);
209 for (i=1;i<=BS->NbKnots();i++)
210 Knots->SetValue(i,BS->Knot(i));
211 Knots->SetValue(BS->NbKnots()+1,Ux);
212 Knots->SetValue(BS->NbKnots()+2,Uy);
213 }
214 }
215
216 TColStd_Array1OfReal knots(1,Knots->Length());
217 knots=Knots->ChangeArray1();
218 SortTools_QuickSortOfReal::Sort(knots,Comp); //sort of the array of knots
219
220 Polesinit(0).SetCoord(0.0,Herm(0)); //poles of the Hermite polynome in the BSpline form
221 Polesinit(1).SetCoord(0.0,Herm(0)+Herm(1)/3.0);
222 Polesinit(2).SetCoord(0.0,Herm(3)-Herm(2)/3.0);
223 Polesinit(3).SetCoord(0.0,Herm(3));
224
225 //loop to check the tolerances on poles
226 if (TolPoles!=0.0){
227 Polemax(Polesinit,min,max);
228 Standard_Real Polemin=Polesinit(min).Y();
229 Standard_Real Polemax=Polesinit(max).Y();
230 if (((Polemax)>=((1/TolPoles)*Polemin))||((Polemin==0.0)&&(Polemax>=(1/TolPoles)))){
231 if (Polesinit(0).Y()>=(1/TolPoles)*Polesinit(3).Y()||Polesinit(0).Y()<=TolPoles*Polesinit(3).Y())
eafb234b 232 Standard_DimensionError::Raise("Hermit Impossible Tolerance");
233 if ((max==0)||(max==3))
234 {
235 for (i=0;i<=3;i++)
236 Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-TolPoles*Polemax));
237 }
238 else if ((max==1)||(max==2)) {
239 if ((min==0)||(min==3))
240 {
241 for (i=0;i<=3;i++)
242 Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-(1/TolPoles)*Polemin));
243 }
244 else{
245 if ((TolPoles*Polemax<Polesinit(0).Y())&&(TolPoles*Polemax<Polesinit(3).Y())){
246 for (i=0;i<=3;i++)
247 Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-TolPoles*Polemax));
248 mark=1;
249 }
250 if ((1/TolPoles*Polemin>Polesinit(0).Y())&&(1/TolPoles*Polemin>Polesinit(3).Y())&&(mark==0)){
251 for (i=0;i<=3;i++)
252 Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-1/TolPoles*Polemin));
253 mark=1;
254 }
255 if (mark==0){
256 Standard_Real Pole0,Pole3;
257 Pole0=Polesinit(0).Y();
258 Pole3=Polesinit(3).Y();
259 if (Pole0<3){
260 a=Log10(Pole3/Pole0);
261 if (boucle==2)
262 {
263 for (i=0;i<=3;i++)
264 Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole3*(Pow(10.0,(-0.5*Log10(TolPoles)-a/2.0)))));
265 }
266 if (boucle==1)
267 {
268 for (i=0;i<=3;i++)
269 Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole0*(Pow(10.0,(a/2.0+0.5*Log10(TolPoles))))));
270 dercas=1;
271 }
272 }
273 if (Pole0>Pole3)
274 {
275 a=Log10(Pole0/Pole3);
276 if (boucle==2)
277 {
278 for (i=0;i<=3;i++)
279 Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole0*(Pow(10.0,(-0.5*Log10(TolPoles)-a/2.0)))));
280 }
281 if (boucle==1)
282 {
283 for (i=0;i<=3;i++)
284 Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole3*(Pow(10.0,(a/2.0+0.5*Log10(TolPoles))))));
285 dercas=1;
286 }
287 }
288 }
289 }
290 }
7fd59977 291 }
eafb234b 292 } // end of the loop
7fd59977 293
eafb234b 294 if (!SignDenom(Polesinit)) //invertion of the polynome sign
295 {
7fd59977 296 for (index=0;index<=3;index++)
297 Polesinit(index).SetCoord(0.0,-Polesinit(index).Y());
eafb234b 298 }
299
300 // loop of positivity
301 if ((Polesinit(1).Y()<0.0)&&(Polesinit(2).Y()>=0.0))
302 {
7fd59977 303 Us1=Polesinit(0).Y()/(Polesinit(0).Y()-Polesinit(1).Y());
304 if (boucle==2)
305 Us1=Us1*knots(2);
306 if (boucle==1)
307 if (Ux!=0.0)
308 Us1=Us1*Ux;
309 BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1);
310 if (I1<2)
311 U4=Us1;
312 else
313 U4=knots(I1);
314 }
315
eafb234b 316 if ((Polesinit(1).Y()>=0.0)&&(Polesinit(2).Y()<0.0))
317 {
7fd59977 318 Us2=Polesinit(2).Y()/(Polesinit(2).Y()-Polesinit(3).Y());
319 if (boucle==2)
320 Us2=knots(knots.Length()-1)+Us2*(1-knots(knots.Length()-1));
321 if (boucle==1)
322 if (Ux!=0.0)
323 Us2=Uy+Us2*(1-Uy);
324 BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I1,Us2);
325 if (I1>=(knots.Length()-1))
326 U5=Us2;
327 else
328 U5=knots(I1+1);
329 }
330
331 if (dercas==1)
332 boucle++;
333
334 if ((Polesinit(1).Y()<0.0)&&(Polesinit(2).Y()<0.0)){
335 Us1=Polesinit(0).Y()/(Polesinit(0).Y()-Polesinit(1).Y());
336 Us2=Polesinit(2).Y()/(Polesinit(2).Y()-Polesinit(3).Y());
337 if (boucle!=0)
338 if (Ux!=0.0){
339 Us1=Us1*Ux;
340 Us2=Uy+Us2*(1-Uy);
341 }
342 if (Us2<=Us1){
343 BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1);
344 if (knots(I1)>=Us2) //insertion of one knot for the two poles
345 U4=knots(I1);
346 else{
347 if (I1>=2){ //insertion to the left and
348 U4=knots(I1); //to the right without a new knot
349 BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I3,Us2);
350 if (I3<(BS->NbKnots()-1)){
351 U5=knots(I3+1);
352 cas=1;
353 }
354 }
355 if(cas==0) //insertion of only one new knot
356 U4=(Us1+Us2)/2;
357 }
358 }
359 else{ //insertion of two knots
360 BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1);
361 if (I1>=2)
362 U4=knots(I1);
363 else
364 U4=Us1;
365 BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I3,Us2);
366 if (I3<(BS->NbKnots()-1))
367 U5=knots(I3+1);
368 else
369 U5=Us2;
370 }
371 }
372}
373
374//=======================================================================
375//function : PolyTest
376//purpose : give the knots U4 and U5 to insert to a(u)
377//=======================================================================
378
379
380static void PolyTest(const TColStd_Array1OfReal& Herm,
381 const Handle(Geom2d_BSplineCurve)& BS,
382 Standard_Real& U4,
383 Standard_Real& U5,
384 Standard_Integer& boucle,
385 const Standard_Real TolPoles,
386// const Standard_Real TolKnots,
387 const Standard_Real ,
388 const Standard_Real Ux,
389 const Standard_Real Uy)
390
391{
392 Standard_Integer index,i,
393 I1=0,I2=0,I3=0,I4=0; //knots index
394 TColgp_Array1OfPnt2d Polesinit(0,3) ;
395 Handle(TColStd_HArray1OfReal) Knots; //array of the BSpline knots + the ones inserted
396 Standard_Integer cas=0,mark=0,dercas=0, //loop marks
397 min,max; //Pole min and max indices
398 Standard_Real Us1,Us2,a; //bounderies value of the knots to be inserted
399// gp_Pnt2d P ;
400 TCollection_CompareOfReal Comp;
401
402 U4=0.0;U5=1.0; //default value
403 if (Ux!=1.0){
404 BS->LocateU(Ux,0.0,I1,I2); //localization of the inserted knots
405 if (Uy!=0.0)
406 BS->LocateU(Uy,0.0,I3,I4);
407 }
408
409 if (I1==I2) //definition and filling of the
eafb234b 410 {
7fd59977 411 if((I3==I4)||(I3==0)){ //array of knots
412 Knots=new TColStd_HArray1OfReal(1,BS->NbKnots());
413 for (i=1;i<=BS->NbKnots();i++)
414 Knots->SetValue(i,BS->Knot(i));
415 }
416 else{
417 Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()+1);
418 for (i=1;i<=BS->NbKnots();i++)
419 Knots->SetValue(i,BS->Knot(i));
420 Knots->SetValue(BS->NbKnots()+1,Uy);
421 }
eafb234b 422 }
423 else
424 {
7fd59977 425 if((I3==I4)||(I3==0)){
426 Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()+1);
427 for (i=1;i<=BS->NbKnots();i++)
428 Knots->SetValue(i,BS->Knot(i));
429 Knots->SetValue(BS->NbKnots()+1,Ux);
430 }
431 else{
432 Knots=new TColStd_HArray1OfReal(1,BS->NbKnots()+2);
433 for (i=1;i<=BS->NbKnots();i++)
434 Knots->SetValue(i,BS->Knot(i));
435 Knots->SetValue(BS->NbKnots()+1,Ux);
436 Knots->SetValue(BS->NbKnots()+2,Uy);
437 }
438 }
439
440 TColStd_Array1OfReal knots(1,Knots->Length());
441 knots=Knots->ChangeArray1();
442 SortTools_QuickSortOfReal::Sort(knots,Comp); //sort of the array of knots
443
444 Polesinit(0).SetCoord(0.0,Herm(0)); //poles of the Hermite polynome in the BSpline form
445 Polesinit(1).SetCoord(0.0,Herm(0)+Herm(1)/3.0);
446 Polesinit(2).SetCoord(0.0,Herm(3)-Herm(2)/3.0);
447 Polesinit(3).SetCoord(0.0,Herm(3));
448
eafb234b 449 // loop to check the tolerances on poles
450 if (TolPoles!=0.0)
451 {
7fd59977 452 Polemax(Polesinit,min,max);
453 Standard_Real Polemin=Polesinit(min).Y();
454 Standard_Real Polemax=Polesinit(max).Y();
eafb234b 455 if (((Polemax)>=((1/TolPoles)*Polemin))||((Polemin==0.0)&&(Polemax>=(1/TolPoles))))
456 {
7fd59977 457 if (Polesinit(0).Y()>=(1/TolPoles)*Polesinit(3).Y()||Polesinit(0).Y()<=TolPoles*Polesinit(3).Y())
eafb234b 458 Standard_DimensionError::Raise("Hermit Impossible Tolerance");
459 if ((max==0)||(max==3))
460 {
461 for (i=0;i<=3;i++)
462 Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-TolPoles*Polemax));
463 }
464 else if ((max==1)||(max==2))
465 {
466 if ((min==0)||(min==3))
467 {
468 for (i=0;i<=3;i++)
469 Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-(1/TolPoles)*Polemin));
470 }
471 else
472 {
473 if ((TolPoles*Polemax<Polesinit(0).Y())&&(TolPoles*Polemax<Polesinit(3).Y()))
474 {
475 for (i=0;i<=3;i++)
476 Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-TolPoles*Polemax));
477 mark=1;
478 }
479
480 if ((1/TolPoles*Polemin>Polesinit(0).Y())&&(1/TolPoles*Polemin>Polesinit(3).Y())&&(mark==0))
481 {
482 for (i=0;i<=3;i++)
483 Polesinit(i).SetCoord(0.0,(Polesinit(i).Y()-1/TolPoles*Polemin));
484 mark=1;
485 }
486 if (mark==0)
487 {
488 Standard_Real Pole0,Pole3;
489 Pole0=Polesinit(0).Y();
490 Pole3=Polesinit(3).Y();
491 if (Pole0<3)
492 {
493 a=Log10(Pole3/Pole0);
494 if (boucle==2)
495 {
496 for (i=0;i<=3;i++)
497 Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole3*(Pow(10.0,(-0.5*Log10(TolPoles)-a/2.0)))));
498 }
499 if (boucle==1)
500 {
501 for (i=0;i<=3;i++)
502 Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole0*(Pow(10.0,(a/2.0+0.5*Log10(TolPoles))))));
503 dercas=1;
504 }
505 }
506 if (Pole0>Pole3)
507 {
508 a=Log10(Pole0/Pole3);
509 if (boucle==2)
510 {
511 for (i=0;i<=3;i++)
512 Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole0*(Pow(10.0,(-0.5*Log10(TolPoles)-a/2.0)))));
513 }
514 else if (boucle==1)
515 {
516 for (i=0;i<=3;i++)
517 Polesinit(i).SetCoord(0.0, Polesinit(i).Y()-(Pole3*(Pow(10.0,(a/2.0+0.5*Log10(TolPoles))))));
518 dercas=1;
519 }
520 }
521 }
522 }
523 }
7fd59977 524 }
eafb234b 525 } // end of the loop
7fd59977 526
eafb234b 527 if (!SignDenom(Polesinit)) // invertion of the polynome sign
528 {
7fd59977 529 for (index=0;index<=3;index++)
530 Polesinit(index).SetCoord(0.0,-Polesinit(index).Y());
eafb234b 531 }
532
533 // boucle de positivite
534 if ((Polesinit(1).Y()<0.0)&&(Polesinit(2).Y()>=0.0))
535 {
7fd59977 536 Us1=Polesinit(0).Y()/(Polesinit(0).Y()-Polesinit(1).Y());
537 if (boucle==2)
538 Us1=Us1*knots(2);
539 if (boucle==1)
540 if (Ux!=0.0)
541 Us1=Us1*Ux;
542 BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1);
543 if (I1<2)
544 U4=Us1;
545 else
546 U4=knots(I1);
547 }
548
eafb234b 549 if ((Polesinit(1).Y()>=0.0)&&(Polesinit(2).Y()<0.0))
550 {
7fd59977 551 Us2=Polesinit(2).Y()/(Polesinit(2).Y()-Polesinit(3).Y());
552 if (boucle==2)
553 Us2=knots(knots.Length()-1)+Us2*(1-knots(knots.Length()-1));
554 if (boucle==1)
555 if (Ux!=0.0)
556 Us2=Uy+Us2*(1-Uy);
557 BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I1,Us2);
558 if (I1>=(knots.Length()-1))
559 U5=Us2;
560 else
561 U5=knots(I1+1);
562 }
563
564 if (dercas==1)
565 boucle++;
566
567 if ((Polesinit(1).Y()<0.0)&&(Polesinit(2).Y()<0.0)){
568 Us1=Polesinit(0).Y()/(Polesinit(0).Y()-Polesinit(1).Y());
569 Us2=Polesinit(2).Y()/(Polesinit(2).Y()-Polesinit(3).Y());
570 if (boucle!=0)
571 if (Ux!=0.0){
572 Us1=Us1*Ux;
573 Us2=Uy+Us2*(1-Uy);
574 }
575 if (Us2<=Us1){
576 BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1);
577 if (knots(I1)>=Us2) //insertion of one knot for the two poles
578 U4=knots(I1);
579 else{
580 if (I1>=2){ //insertion to the left and
581 U4=knots(I1); //to the right without a new knot
582 BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I3,Us2);
583 if (I3<(BS->NbKnots()-1)){
584 U5=knots(I3+1);
585 cas=1;
586 }
587 }
588 if(cas==0) //insertion of only one new knot
589 U4=(Us1+Us2)/2;
590 }
591 }
592 else{ //insertion of two knots
593 BSplCLib::LocateParameter(3,knots,Us1,Standard_False,1,knots.Length(),I1,Us1);
594 if (I1>=2)
595 U4=knots(I1);
596 else
597 U4=Us1;
598 BSplCLib::LocateParameter(3,knots,Us2,Standard_False,1,knots.Length(),I3,Us2);
599 if (I3<(BS->NbKnots()-1))
600 U5=knots(I3+1);
601 else
602 U5=Us2;
603 }
604 }
605}
606
607//=======================================================================
608//function : InsertKnots
609//purpose : insert the knots in BS knot sequence if they are not null.
610//=======================================================================
611
612static void InsertKnots(Handle(Geom2d_BSplineCurve)& BS,
613 const Standard_Real U4,
614 const Standard_Real U5)
615
616{
617 if (U4!=0.0) //insertion of :0 knot if U4=0
618 BS->InsertKnot(U4); // 1 knot if U4=U5
619 if ((U5!=1.0)&&(U5!=U4)) // 2 knots otherwise
620 BS->InsertKnot(U5);
621
622}
623
624//=======================================================================
625//function : MovePoles
626//purpose : move the poles above the x axis
627//=======================================================================
628
629static void MovePoles(Handle(Geom2d_BSplineCurve)& BS)
630
631{
632 gp_Pnt2d P ;
633// Standard_Integer i,index;
634 Standard_Integer i;
635
636 for (i=3;i<=(BS->NbPoles()-2);i++){
637 P.SetCoord(1,(BS->Pole(i).Coord(1))); //raising of the no constrained poles to
638 P.SetCoord(2,(BS->Pole(1).Coord(2))); //the first pole level
639 BS->SetPole(i,P);
640 }
641}
642
643//=======================================================================
644//function : Solution
645//purpose :
646//=======================================================================
647
648Handle(Geom2d_BSplineCurve) Hermit::Solution(const Handle(Geom_BSplineCurve)& BS,
649 const Standard_Real TolPoles,
650 const Standard_Real TolKnots)
651
652{
653 TColStd_Array1OfReal Herm(0,3);
654 Standard_Real Upos1=0.0, Upos2=1.0, //positivity knots
655 Ux=0.0, Uy=1.0,
656 Utol1=0.0, Utol2=1.0, //tolerance knots
657 Uint1=0.0, Uint2=1.0; //tolerance knots for the first loop
658 Standard_Integer boucle=1; //loop mark
659 TColStd_Array1OfReal Knots(1,2);
660 TColStd_Array1OfInteger Multiplicities(1,2);
661 TColgp_Array1OfPnt2d Poles(1,4);
662 Standard_Integer zeroboucle = 0 ;
663 HermiteCoeff(BS,Herm); //computation of the Hermite coefficient
664
665 Poles(1).SetCoord(0.0,Herm(0)); //poles of the Hermite polynome in the BSpline form
666 Poles(2).SetCoord(0.0,Herm(0)+Herm(1)/3.0);
667 Poles(3).SetCoord(0.0,Herm(3)-Herm(2)/3.0);
668 Poles(4).SetCoord(0.0,Herm(3));
669 Knots(1)=0.0;
670 Knots(2)=1.0;
671 Multiplicities(1)=4;
672 Multiplicities(2)=4;
673
674 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(Poles,Knots,Multiplicities,3);//creation of the basic
675 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(Poles,Knots,Multiplicities,3);//BSpline without modif
676
677 PolyTest(Herm,BS,Upos1,Upos2,zeroboucle,Precision::Confusion(),Precision::Confusion(),1.0,0.0);//computation of the positivity knots
678 InsertKnots(BS2,Upos1,Upos2); //and insertion
679
680 if (Upos1!=0.0)
681 if (Upos2!=1.0){
682 Ux=Min(Upos1,Upos2);
683 Uy=Max(Upos1,Upos2);
684 }
685 else{
686 Ux=Upos1;
687 Uy=Upos1;
688 }
689 else{
690 Ux=Upos2;
691 Uy=Upos2;
692 }
693
694 Herm(0)=BS2->Pole(1).Y(); //computation of the Hermite coefficient on the
695 Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y()); //positive BSpline
696 Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y());
697 Herm(3)=BS2->Pole(BS2->NbPoles()).Y();
698
699 PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Ux,Uy); //computation of the tolerance knots
700 InsertKnots(BS2,Utol1,Utol2); //and insertion
701
702 if (boucle==2){ //insertion of two knots
703 Herm(0)=BS2->Pole(1).Y();
704 Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y());
705 Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y());
706 Herm(3)=BS2->Pole(BS2->NbPoles()).Y();
707 if (Utol1==0.0){
708 Uint2=Utol2;
709 PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint2,0.0);
710 }
711 else{
712 Uint1=Utol1;
713 PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint1,0.0);
714 }
715 InsertKnots(BS2,Utol1,Utol2);
716 }
717 if ((BS2->Knot(2)<TolKnots)||(BS2->Knot(BS2->NbKnots()-1)>(1-TolKnots))) //checking of the knots tolerance
718 Standard_DimensionError::Raise("Hermit Impossible Tolerance");
719 else{
720 if ((Upos2==1.0)&&(Utol2==1.0)&&(Uint2==1.0)) //test on the final inserted knots
721 InsertKnots(BS1,BS2->Knot(2),1.0);
722 else{
723 if ((Upos1==0.0)&&(Utol1==0.0)&&(Uint1==0.0))
724 InsertKnots(BS1,BS2->Knot(BS2->NbKnots()-1),1.0);
725 else
726 InsertKnots(BS1,BS2->Knot(BS2->NbKnots()-1),BS2->Knot(2));
727 }
728 MovePoles(BS1); //relocation of the no-contrained knots
729 }
730 return BS1;
731}
732
733
734//=======================================================================
735// Solution
736//=======================================================================
737
738Handle(Geom2d_BSplineCurve) Hermit::Solution(const Handle(Geom2d_BSplineCurve)& BS,
739 const Standard_Real TolPoles,
740 const Standard_Real TolKnots)
741
742{
743 TColStd_Array1OfReal Herm(0,3);
744 Standard_Real Upos1=0.0, Upos2=1.0, //positivity knots
745 Ux=0.0, Uy=1.0,
746 Utol1=0.0, Utol2=1.0, //tolerance knots
747 Uint1=0.0, Uint2=1.0; //tolerance knots for the first loop
748 Standard_Integer boucle=1; //loop mark
749 TColStd_Array1OfReal Knots(1,2);
750 TColStd_Array1OfInteger Multiplicities(1,2);
751 TColgp_Array1OfPnt2d Poles(1,4);
752 Standard_Integer zeroboucle = 0 ;
753 HermiteCoeff(BS,Herm); //computation of the Hermite coefficient
754
755 Poles(1).SetCoord(0.0,Herm(0)); //poles of the Hermite polynome in the BSpline form
756 Poles(2).SetCoord(0.0,Herm(0)+Herm(1)/3.0);
757 Poles(3).SetCoord(0.0,Herm(3)-Herm(2)/3.0);
758 Poles(4).SetCoord(0.0,Herm(3));
759 Knots(1)=0.0;
760 Knots(2)=1.0;
761 Multiplicities(1)=4;
762 Multiplicities(2)=4;
763
764 Handle(Geom2d_BSplineCurve) BS1=new Geom2d_BSplineCurve(Poles,Knots,Multiplicities,3);//creation of the basic
765 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(Poles,Knots,Multiplicities,3);//BSpline without modif
766
767 PolyTest(Herm,BS,Upos1,Upos2,zeroboucle,Precision::Confusion(),Precision::Confusion(),1.0,0.0);//computation of the positivity knots
768 InsertKnots(BS2,Upos1,Upos2); //and insertion
769
770 if (Upos1!=0.0)
771 if (Upos2!=1.0){
772 Ux=Min(Upos1,Upos2);
773 Uy=Max(Upos1,Upos2);
774 }
775 else{
776 Ux=Upos1;
777 Uy=Upos1;
778 }
779 else{
780 Ux=Upos2;
781 Uy=Upos2;
782 }
783
784 Herm(0)=BS2->Pole(1).Y(); //computation of the Hermite coefficient on the
785 Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y()); //positive BSpline
786 Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y());
787 Herm(3)=BS2->Pole(BS2->NbPoles()).Y();
788
789 PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Ux,Uy); //computation of the tolerance knots
790 InsertKnots(BS2,Utol1,Utol2); //and insertion
791
792 if (boucle==2){ //insertion of two knots
793 Herm(0)=BS2->Pole(1).Y();
794 Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y());
795 Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y());
796 Herm(3)=BS2->Pole(BS2->NbPoles()).Y();
797 if (Utol1==0.0){
798 Uint2=Utol2;
799 PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint2,0.0);
800 }
801 else{
802 Uint1=Utol1;
803 PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint1,0.0);
804 }
805 InsertKnots(BS2,Utol1,Utol2);
806 }
807 if ((BS2->Knot(2)<TolKnots)||(BS2->Knot(BS2->NbKnots()-1)>(1-TolKnots))) //checking of the knots tolerance
808 Standard_DimensionError::Raise("Hermit Impossible Tolerance");
809 else{
810 if ((Upos2==1.0)&&(Utol2==1.0)&&(Uint2==1.0)) //test on the final inserted knots
811 InsertKnots(BS1,BS2->Knot(2),1.0);
812 else{
813 if ((Upos1==0.0)&&(Utol1==0.0)&&(Uint1==0.0))
814 InsertKnots(BS1,BS2->Knot(BS2->NbKnots()-1),1.0);
815 else
816 InsertKnots(BS1,BS2->Knot(BS2->NbKnots()-1),BS2->Knot(2));
817 }
818 MovePoles(BS1); //relocation of the no-contrained knots
819 }
820 return BS1;
821}
822
823//=======================================================================
824//function : Solutionbis
825//purpose :
826//=======================================================================
827
828void Hermit::Solutionbis(const Handle(Geom_BSplineCurve)& BS,
829 Standard_Real & Knotmin,
830 Standard_Real & Knotmax,
831 const Standard_Real TolPoles,
832 const Standard_Real TolKnots)
833
834{
835 TColStd_Array1OfReal Herm(0,3);
836 Standard_Real Upos1=0.0, Upos2=1.0, //positivity knots
837 Ux=0.0, Uy=1.0,
838 Utol1=0.0, Utol2=1.0, //tolerance knots
839 Uint1=0.0, Uint2=1.0; //tolerance knots for the first loop
840 Standard_Integer boucle=1; //loop mark
841 TColStd_Array1OfReal Knots(1,2);
842 TColStd_Array1OfInteger Multiplicities(1,2);
843 TColgp_Array1OfPnt2d Poles(1,4);
844 Standard_Integer zeroboucle = 0 ;
845 HermiteCoeff(BS,Herm); //computation of the Hermite coefficient
846
847 Poles(1).SetCoord(0.0,Herm(0)); //poles of the Hermite polynome in the BSpline form
848 Poles(2).SetCoord(0.0,Herm(0)+Herm(1)/3.0);
849 Poles(3).SetCoord(0.0,Herm(3)-Herm(2)/3.0);
850 Poles(4).SetCoord(0.0,Herm(3));
851 Knots(1)=0.0;
852 Knots(2)=1.0;
853 Multiplicities(1)=4;
854 Multiplicities(2)=4;
855
856 Handle(Geom2d_BSplineCurve) BS2=new Geom2d_BSplineCurve(Poles,Knots,Multiplicities,3);//creation of the basic
857 //BSpline without modif
858
859 PolyTest(Herm,BS,Upos1,Upos2,zeroboucle,Precision::Confusion(),Precision::Confusion(),1.0,0.0);//computation of the positivity knots
860 InsertKnots(BS2,Upos1,Upos2); //and insertion
861
862 if (Upos1!=0.0)
863 if (Upos2!=1.0){
864 Ux=Min(Upos1,Upos2);
865 Uy=Max(Upos1,Upos2);
866 }
867 else{
868 Ux=Upos1;
869 Uy=Upos1;
870 }
871 else{
872 Ux=Upos2;
873 Uy=Upos2;
874 }
875
876 Herm(0)=BS2->Pole(1).Y(); //computation of the Hermite coefficient on the
877 Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y()); //positive BSpline
878 Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y());
879 Herm(3)=BS2->Pole(BS2->NbPoles()).Y();
880
881 PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Ux,Uy); //computation of the tolerance knots
882 InsertKnots(BS2,Utol1,Utol2); //and insertion
883
884 if (boucle==2){ //insertion of two knots
885 Herm(0)=BS2->Pole(1).Y();
886 Herm(1)=3*(BS2->Pole(2).Y()-BS2->Pole(1).Y());
887 Herm(2)=3*(BS2->Pole(BS2->NbPoles()).Y()-BS2->Pole(BS2->NbPoles()-1).Y());
888 Herm(3)=BS2->Pole(BS2->NbPoles()).Y();
889 if (Utol1==0.0){
890 Uint2=Utol2;
891 PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint2,0.0);
892 }
893 else{
894 Uint1=Utol1;
895 PolyTest(Herm,BS,Utol1,Utol2,boucle,TolPoles,TolKnots,Uint1,0.0);
896 }
897 InsertKnots(BS2,Utol1,Utol2);
898 }
899 if ((BS2->Knot(2)<TolKnots)||(BS2->Knot(BS2->NbKnots()-1)>(1-TolKnots))) //checking of the knots tolerance
900 Standard_DimensionError::Raise("Hermit Impossible Tolerance");
901 else{
902 if ((Upos2==1.0)&&(Utol2==1.0)&&(Uint2==1.0)) //test on the final inserted knots
903 Knotmin=BS2->Knot(2);
904 else{
905 if ((Upos1==0.0)&&(Utol1==0.0)&&(Uint1==0.0))
906 Knotmax=BS2->Knot(BS2->NbKnots()-1);
907 else{
908 Knotmin=BS2->Knot(2);
909 Knotmax=BS2->Knot(BS2->NbKnots()-1);
910 }
911 }
912 }
913}
914
915
916
Open CASCADE Technology repository