Warning, /include/opencascade/BSplCLib_CurveComputation.gxx is written in an unsupported language. File is not indexed.
0001 // Copyright (c) 1995-1999 Matra Datavision
0002 // Copyright (c) 1999-2014 OPEN CASCADE SAS
0003 //
0004 // This file is part of Open CASCADE Technology software library.
0005 //
0006 // This library is free software; you can redistribute it and/or modify it under
0007 // the terms of the GNU Lesser General Public License version 2.1 as published
0008 // by the Free Software Foundation, with special exception defined in the file
0009 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0010 // distribution for complete text of the license and disclaimer of any warranty.
0011 //
0012 // Alternatively, this file may be used under the terms of Open CASCADE
0013 // commercial license or contractual agreement.
0014
0015 // xab : modified 15-Mar-94 added EvalBSplineMatrix, FactorBandedMatrix, SolveBandedSystem
0016 // Eval, BuildCache, cacheD0, cacheD1, cacheD2, cacheD3, LocalMatrix,
0017 // EvalPolynomial, RationalDerivatives.
0018 // xab : 22-Mar-94 : fixed rational problem in BuildCache
0019 // xab : 12-Mar-96 : added MovePointAndTangent
0020 #include <TColStd_Array1OfInteger.hxx>
0021 #include <TColStd_Array1OfReal.hxx>
0022 #include <TColStd_Array2OfReal.hxx>
0023 #include <gp_Vec2d.hxx>
0024 #include <Standard_ConstructionError.hxx>
0025 #include <PLib.hxx>
0026 #include <math_Matrix.hxx>
0027
0028 //=======================================================================
0029 //struct : BSplCLib_DataContainer
0030 //purpose: Auxiliary structure providing buffers for poles and knots used in
0031 // evaluation of bspline (allocated in the stack)
0032 //=======================================================================
0033
0034 struct BSplCLib_DataContainer
0035 {
0036 BSplCLib_DataContainer(Standard_Integer Degree)
0037 {
0038 (void)Degree; // avoid compiler warning
0039 Standard_OutOfRange_Raise_if (Degree > BSplCLib::MaxDegree() ||
0040 BSplCLib::MaxDegree() > 25,
0041 "BSplCLib: bspline degree is greater than maximum supported");
0042 }
0043
0044 Standard_Real poles[(25+1)*(Dimension_gen+1)];
0045 Standard_Real knots[2*25];
0046 Standard_Real ders[Dimension_gen*4];
0047 };
0048
0049 //=======================================================================
0050 //function : Reverse
0051 //purpose :
0052 //=======================================================================
0053
0054 void BSplCLib::Reverse(Array1OfPoints& Poles,
0055 const Standard_Integer L)
0056 {
0057 Standard_Integer i, l = L;
0058 l = Poles.Lower() + (l-Poles.Lower()) % (Poles.Upper()-Poles.Lower()+1);
0059
0060 Array1OfPoints temp(0,Poles.Length()-1);
0061
0062 for (i = Poles.Lower(); i <= l; i++)
0063 temp(l-i) = Poles(i);
0064
0065 for (i = l+1; i <= Poles.Upper(); i++)
0066 temp(l-Poles.Lower()+Poles.Upper()-i+1) = Poles(i);
0067
0068 for (i = Poles.Lower(); i <= Poles.Upper(); i++)
0069 Poles(i) = temp(i-Poles.Lower());
0070 }
0071
0072 //
0073 // CURVES MODIFICATIONS
0074 //
0075
0076 //=======================================================================
0077 //function : RemoveKnot
0078 //purpose :
0079 //=======================================================================
0080
0081 Standard_Boolean BSplCLib::RemoveKnot
0082 (const Standard_Integer Index,
0083 const Standard_Integer Mult,
0084 const Standard_Integer Degree,
0085 const Standard_Boolean Periodic,
0086 const Array1OfPoints& Poles,
0087 const TColStd_Array1OfReal* Weights,
0088 const TColStd_Array1OfReal& Knots,
0089 const TColStd_Array1OfInteger& Mults,
0090 Array1OfPoints& NewPoles,
0091 TColStd_Array1OfReal* NewWeights,
0092 TColStd_Array1OfReal& NewKnots,
0093 TColStd_Array1OfInteger& NewMults,
0094 const Standard_Real Tolerance)
0095 {
0096 Standard_Boolean rational = Weights != NULL;
0097 Standard_Integer dim;
0098 dim = Dimension_gen;
0099 if (rational) dim++;
0100
0101 TColStd_Array1OfReal poles(1,dim*Poles.Length());
0102 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
0103
0104 if (rational) PLib::SetPoles(Poles,*Weights,poles);
0105 else PLib::SetPoles(Poles,poles);
0106
0107 if (!RemoveKnot(Index,Mult,Degree,Periodic,dim,
0108 poles,Knots,Mults,newpoles,NewKnots,NewMults,Tolerance))
0109 return Standard_False;
0110
0111 if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
0112 else PLib::GetPoles(newpoles,NewPoles);
0113 return Standard_True;
0114 }
0115
0116 //=======================================================================
0117 //function : InsertKnots
0118 //purpose : insert an array of knots and multiplicities
0119 //=======================================================================
0120
0121 void BSplCLib::InsertKnots
0122 (const Standard_Integer Degree,
0123 const Standard_Boolean Periodic,
0124 const Array1OfPoints& Poles,
0125 const TColStd_Array1OfReal* Weights,
0126 const TColStd_Array1OfReal& Knots,
0127 const TColStd_Array1OfInteger& Mults,
0128 const TColStd_Array1OfReal& AddKnots,
0129 const TColStd_Array1OfInteger* AddMults,
0130 Array1OfPoints& NewPoles,
0131 TColStd_Array1OfReal* NewWeights,
0132 TColStd_Array1OfReal& NewKnots,
0133 TColStd_Array1OfInteger& NewMults,
0134 const Standard_Real Epsilon,
0135 const Standard_Boolean Add)
0136 {
0137 Standard_Boolean rational = Weights != NULL;
0138 Standard_Integer dim;
0139 dim = Dimension_gen;
0140 if (rational) dim++;
0141
0142 TColStd_Array1OfReal poles(1,dim*Poles.Length());
0143 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
0144
0145 if (rational) PLib::SetPoles(Poles,*Weights,poles);
0146 else PLib::SetPoles(Poles,poles);
0147
0148 InsertKnots(Degree,Periodic,dim,poles,Knots,Mults,
0149 AddKnots,AddMults,newpoles,NewKnots,NewMults,Epsilon,Add);
0150
0151 if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
0152 else PLib::GetPoles(newpoles,NewPoles);
0153 }
0154
0155 //=======================================================================
0156 //function : InsertKnot
0157 //purpose :
0158 //=======================================================================
0159
0160 void BSplCLib::InsertKnot(const Standard_Integer ,
0161 const Standard_Real U,
0162 const Standard_Integer UMult,
0163 const Standard_Integer Degree,
0164 const Standard_Boolean Periodic,
0165 const Array1OfPoints& Poles,
0166 const TColStd_Array1OfReal* Weights,
0167 const TColStd_Array1OfReal& Knots,
0168 const TColStd_Array1OfInteger& Mults,
0169 Array1OfPoints& NewPoles,
0170 TColStd_Array1OfReal* NewWeights)
0171 {
0172 TColStd_Array1OfReal k(1,1);
0173 k(1) = U;
0174 TColStd_Array1OfInteger m(1,1);
0175 m(1) = UMult;
0176 TColStd_Array1OfReal nk(1,Knots.Length()+1);
0177 TColStd_Array1OfInteger nm(1,Knots.Length()+1);
0178 InsertKnots(Degree,Periodic,Poles,Weights,Knots,Mults,
0179 k,&m,NewPoles,NewWeights,nk,nm,Epsilon(U));
0180 }
0181
0182 //=======================================================================
0183 //function : RaiseMultiplicity
0184 //purpose :
0185 //=======================================================================
0186
0187 void BSplCLib::RaiseMultiplicity(const Standard_Integer KnotIndex,
0188 const Standard_Integer Mult,
0189 const Standard_Integer Degree,
0190 const Standard_Boolean Periodic,
0191 const Array1OfPoints& Poles,
0192 const TColStd_Array1OfReal* Weights,
0193 const TColStd_Array1OfReal& Knots,
0194 const TColStd_Array1OfInteger& Mults,
0195 Array1OfPoints& NewPoles,
0196 TColStd_Array1OfReal* NewWeights)
0197 {
0198 TColStd_Array1OfReal k(1,1);
0199 k(1) = Knots(KnotIndex);
0200 TColStd_Array1OfInteger m(1,1);
0201 m(1) = Mult - Mults(KnotIndex);
0202 TColStd_Array1OfReal nk(1,Knots.Length());
0203 TColStd_Array1OfInteger nm(1,Knots.Length());
0204 InsertKnots(Degree,Periodic,Poles,Weights,Knots,Mults,
0205 k,&m,NewPoles,NewWeights,nk,nm,Epsilon(k(1)));
0206 }
0207
0208 //=======================================================================
0209 //function : IncreaseDegree
0210 //purpose :
0211 //=======================================================================
0212
0213 void BSplCLib::IncreaseDegree
0214 (const Standard_Integer Degree,
0215 const Standard_Integer NewDegree,
0216 const Standard_Boolean Periodic,
0217 const Array1OfPoints& Poles,
0218 const TColStd_Array1OfReal* Weights,
0219 const TColStd_Array1OfReal& Knots,
0220 const TColStd_Array1OfInteger& Mults,
0221 Array1OfPoints& NewPoles,
0222 TColStd_Array1OfReal* NewWeights,
0223 TColStd_Array1OfReal& NewKnots,
0224 TColStd_Array1OfInteger& NewMults)
0225 {
0226 Standard_Boolean rational = Weights != NULL;
0227 Standard_Integer dim;
0228 dim = Dimension_gen;
0229 if (rational) dim++;
0230
0231 TColStd_Array1OfReal poles(1,dim*Poles.Length());
0232 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
0233
0234 if (rational) PLib::SetPoles(Poles,*Weights,poles);
0235 else PLib::SetPoles(Poles,poles);
0236
0237 IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults,
0238 newpoles,NewKnots,NewMults);
0239
0240 if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
0241 else PLib::GetPoles(newpoles,NewPoles);
0242 }
0243
0244 //=======================================================================
0245 //function : Unperiodize
0246 //purpose :
0247 //=======================================================================
0248
0249 void BSplCLib::Unperiodize
0250 (const Standard_Integer Degree,
0251 const TColStd_Array1OfInteger& Mults,
0252 const TColStd_Array1OfReal& Knots,
0253 const Array1OfPoints& Poles,
0254 const TColStd_Array1OfReal* Weights,
0255 TColStd_Array1OfInteger& NewMults,
0256 TColStd_Array1OfReal& NewKnots,
0257 Array1OfPoints& NewPoles,
0258 TColStd_Array1OfReal* NewWeights)
0259 {
0260 Standard_Boolean rational = Weights != NULL;
0261 Standard_Integer dim;
0262 dim = Dimension_gen;
0263 if (rational) dim++;
0264
0265 TColStd_Array1OfReal poles(1,dim*Poles.Length());
0266 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
0267
0268 if (rational) PLib::SetPoles(Poles,*Weights,poles);
0269 else PLib::SetPoles(Poles,poles);
0270
0271 Unperiodize(Degree, dim, Mults, Knots, poles,
0272 NewMults,NewKnots, newpoles);
0273
0274 if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
0275 else PLib::GetPoles(newpoles,NewPoles);
0276 }
0277
0278 //=======================================================================
0279 //function : Trimming
0280 //purpose :
0281 //=======================================================================
0282
0283 void BSplCLib::Trimming(const Standard_Integer Degree,
0284 const Standard_Boolean Periodic,
0285 const TColStd_Array1OfReal& Knots,
0286 const TColStd_Array1OfInteger& Mults,
0287 const Array1OfPoints& Poles,
0288 const TColStd_Array1OfReal* Weights,
0289 const Standard_Real U1,
0290 const Standard_Real U2,
0291 TColStd_Array1OfReal& NewKnots,
0292 TColStd_Array1OfInteger& NewMults,
0293 Array1OfPoints& NewPoles,
0294 TColStd_Array1OfReal* NewWeights)
0295 {
0296 Standard_Boolean rational = Weights != NULL;
0297 Standard_Integer dim;
0298 dim = Dimension_gen;
0299 if (rational) dim++;
0300
0301 TColStd_Array1OfReal poles(1,dim*Poles.Length());
0302 TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
0303
0304 if (rational) PLib::SetPoles(Poles,*Weights,poles);
0305 else PLib::SetPoles(Poles,poles);
0306
0307 Trimming(Degree, Periodic, dim, Knots, Mults, poles, U1, U2,
0308 NewKnots, NewMults, newpoles);
0309
0310 if (rational) PLib::GetPoles(newpoles,NewPoles,*NewWeights);
0311 else PLib::GetPoles(newpoles,NewPoles);
0312 }
0313
0314 //--------------------------------------------------------------------------
0315 //ELEMENTARY COMPUTATIONS
0316 //--------------------------------------------------------------------------
0317 //
0318 // All the following methods are methods of low level used in BSplCLib
0319 // but also in BSplSLib (but not in Geom)
0320 // At the creation time of this package there is no possibility
0321 // to declare private methods of package and to declare friend
0322 // methods of package. It could interesting to declare that BSplSLib
0323 // is a package friend to BSplCLib because it uses all the basis methods
0324 // of BSplCLib.
0325 //--------------------------------------------------------------------------
0326
0327 //=======================================================================
0328 //function : BuildEval
0329 //purpose : builds the local array for evaluation
0330 //=======================================================================
0331
0332 void BSplCLib::BuildEval(const Standard_Integer Degree,
0333 const Standard_Integer Index,
0334 const Array1OfPoints& Poles,
0335 const TColStd_Array1OfReal* Weights,
0336 Standard_Real& LP)
0337 {
0338 Standard_Real w, *pole = &LP;
0339 Standard_Integer PLower = Poles.Lower();
0340 Standard_Integer PUpper = Poles.Upper();
0341 Standard_Integer i;
0342 Standard_Integer ip = PLower + Index - 1;
0343 if (Weights == NULL) {
0344 for (i = 0; i <= Degree; i++) {
0345 ip++;
0346 if (ip > PUpper) ip = PLower;
0347 const Point& P = Poles(ip);
0348 PointToCoords (pole,P,+0);
0349 pole += Dimension_gen;
0350 }
0351 }
0352 else {
0353 for (i = 0; i <= Degree; i++) {
0354 ip++;
0355 if (ip > PUpper) ip = PLower;
0356 const Point& P = Poles(ip);
0357 pole[Dimension_gen] = w = (*Weights)(ip);
0358 PointToCoords (pole, P, * w);
0359 pole += Dimension_gen + 1;
0360 }
0361 }
0362 }
0363
0364 //=======================================================================
0365 //function : PrepareEval
0366 //purpose : stores data for Eval in the local arrays
0367 // dc.poles and dc.knots
0368 //=======================================================================
0369
0370 static void PrepareEval
0371 (Standard_Real& u,
0372 Standard_Integer& index,
0373 Standard_Integer& dim,
0374 Standard_Boolean& rational,
0375 const Standard_Integer Degree,
0376 const Standard_Boolean Periodic,
0377 const Array1OfPoints& Poles,
0378 const TColStd_Array1OfReal* Weights,
0379 const TColStd_Array1OfReal& Knots,
0380 const TColStd_Array1OfInteger* Mults,
0381 BSplCLib_DataContainer& dc)
0382 {
0383 // Set the Index
0384 BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
0385
0386 // make the knots
0387 BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*dc.knots);
0388 if (Mults == NULL)
0389 index -= Knots.Lower() + Degree;
0390 else
0391 index = BSplCLib::PoleIndex(Degree,index,Periodic,*Mults);
0392
0393 // check truly rational
0394 rational = (Weights != NULL);
0395 if (rational) {
0396 Standard_Integer WLower = Weights->Lower() + index;
0397 rational = BSplCLib::IsRational(*Weights, WLower, WLower + Degree);
0398 }
0399
0400 // make the poles
0401 if (rational) {
0402 dim = Dimension_gen + 1;
0403 BSplCLib::BuildEval(Degree, index, Poles, Weights , *dc.poles);
0404 }
0405 else {
0406 dim = Dimension_gen;
0407 BSplCLib::BuildEval(Degree, index, Poles, BSplCLib::NoWeights(), *dc.poles);
0408 }
0409 }
0410
0411 //=======================================================================
0412 //function : D0
0413 //purpose :
0414 //=======================================================================
0415
0416 void BSplCLib::D0
0417 (const Standard_Real U,
0418 const Standard_Integer Index,
0419 const Standard_Integer Degree,
0420 const Standard_Boolean Periodic,
0421 const Array1OfPoints& Poles,
0422 const TColStd_Array1OfReal* Weights,
0423 const TColStd_Array1OfReal& Knots,
0424 const TColStd_Array1OfInteger* Mults,
0425 Point& P)
0426 {
0427 // Standard_Integer k,dim,index = Index;
0428 Standard_Integer dim,index = Index;
0429 Standard_Real u = U;
0430 Standard_Boolean rational;
0431 BSplCLib_DataContainer dc(Degree);
0432 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
0433 BSplCLib::Eval(u,Degree,*dc.knots,dim,*dc.poles);
0434
0435 if (rational) {
0436 Standard_Real w = dc.poles[Dimension_gen];
0437 CoordsToPoint (P, dc.poles, / w);
0438 }
0439 else
0440 CoordsToPoint (P, dc.poles, );
0441 }
0442
0443 //=======================================================================
0444 //function : D1
0445 //purpose :
0446 //=======================================================================
0447
0448 void BSplCLib::D1
0449 (const Standard_Real U,
0450 const Standard_Integer Index,
0451 const Standard_Integer Degree,
0452 const Standard_Boolean Periodic,
0453 const Array1OfPoints& Poles,
0454 const TColStd_Array1OfReal* Weights,
0455 const TColStd_Array1OfReal& Knots,
0456 const TColStd_Array1OfInteger* Mults,
0457 Point& P,
0458 Vector& V)
0459 {
0460 Standard_Integer dim,index = Index;
0461 Standard_Real u = U;
0462 Standard_Boolean rational;
0463 BSplCLib_DataContainer dc(Degree);
0464 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
0465 BSplCLib::Bohm(u,Degree,1,*dc.knots,dim,*dc.poles);
0466 Standard_Real *result = dc.poles;
0467 if (rational) {
0468 PLib::RationalDerivative(Degree,1,Dimension_gen,*dc.poles,*dc.ders);
0469 result = dc.ders;
0470 }
0471
0472 CoordsToPoint (P, result, );
0473 CoordsToPoint (V, result + Dimension_gen, );
0474 }
0475
0476 //=======================================================================
0477 //function : D2
0478 //purpose :
0479 //=======================================================================
0480
0481 void BSplCLib::D2
0482 (const Standard_Real U,
0483 const Standard_Integer Index,
0484 const Standard_Integer Degree,
0485 const Standard_Boolean Periodic,
0486 const Array1OfPoints& Poles,
0487 const TColStd_Array1OfReal* Weights,
0488 const TColStd_Array1OfReal& Knots,
0489 const TColStd_Array1OfInteger* Mults,
0490 Point& P,
0491 Vector& V1,
0492 Vector& V2)
0493 {
0494 Standard_Integer dim,index = Index;
0495 Standard_Real u = U;
0496 Standard_Boolean rational;
0497 BSplCLib_DataContainer dc(Degree);
0498 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
0499 BSplCLib::Bohm(u,Degree,2,*dc.knots,dim,*dc.poles);
0500 Standard_Real *result = dc.poles;
0501 if (rational) {
0502 PLib::RationalDerivative(Degree,2,Dimension_gen,*dc.poles,*dc.ders);
0503 result = dc.ders;
0504 }
0505
0506 CoordsToPoint (P, result, );
0507 CoordsToPoint (V1, result + Dimension_gen, );
0508 if (!rational && (Degree < 2))
0509 NullifyPoint (V2);
0510 else
0511 CoordsToPoint (V2, result + 2 * Dimension_gen, );
0512 }
0513
0514 //=======================================================================
0515 //function : D3
0516 //purpose :
0517 //=======================================================================
0518
0519 void BSplCLib::D3
0520 (const Standard_Real U,
0521 const Standard_Integer Index,
0522 const Standard_Integer Degree,
0523 const Standard_Boolean Periodic,
0524 const Array1OfPoints& Poles,
0525 const TColStd_Array1OfReal* Weights,
0526 const TColStd_Array1OfReal& Knots,
0527 const TColStd_Array1OfInteger* Mults,
0528 Point& P,
0529 Vector& V1,
0530 Vector& V2,
0531 Vector& V3)
0532 {
0533 Standard_Integer dim,index = Index;
0534 Standard_Real u = U;
0535 Standard_Boolean rational;
0536 BSplCLib_DataContainer dc(Degree);
0537 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
0538 BSplCLib::Bohm(u,Degree,3,*dc.knots,dim,*dc.poles);
0539 Standard_Real *result = dc.poles;
0540 if (rational) {
0541 PLib::RationalDerivative(Degree,3,Dimension_gen,*dc.poles,*dc.ders);
0542 result = dc.ders;
0543 }
0544
0545 CoordsToPoint (P, result, );
0546 CoordsToPoint (V1, result + Dimension_gen, );
0547 if (!rational && (Degree < 2))
0548 NullifyPoint (V2);
0549 else
0550 CoordsToPoint (V2, result + 2 * Dimension_gen, );
0551 if (!rational && (Degree < 3))
0552 NullifyPoint (V3);
0553 else
0554 CoordsToPoint (V3, result + 3 * Dimension_gen, );
0555 }
0556
0557 //=======================================================================
0558 //function : DN
0559 //purpose :
0560 //=======================================================================
0561
0562 void BSplCLib::DN
0563 (const Standard_Real U,
0564 const Standard_Integer N,
0565 const Standard_Integer Index,
0566 const Standard_Integer Degree,
0567 const Standard_Boolean Periodic,
0568 const Array1OfPoints& Poles,
0569 const TColStd_Array1OfReal* Weights,
0570 const TColStd_Array1OfReal& Knots,
0571 const TColStd_Array1OfInteger* Mults,
0572 Vector& VN)
0573 {
0574 Standard_Integer dim,index = Index;
0575 Standard_Real u = U;
0576 Standard_Boolean rational;
0577 BSplCLib_DataContainer dc(Degree);
0578 PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
0579 BSplCLib::Bohm(u,Degree,N,*dc.knots,dim,*dc.poles);
0580
0581 if (rational) {
0582 Standard_Real v[Dimension_gen];
0583 PLib::RationalDerivative(Degree,N,Dimension_gen,*dc.poles,v[0],Standard_False);
0584 CoordsToPoint (VN, v, );
0585 }
0586 else {
0587 if (N > Degree)
0588 NullifyPoint (VN);
0589 else {
0590 Standard_Real *DN = dc.poles + N * Dimension_gen;
0591 CoordsToPoint (VN, DN, );
0592 }
0593 }
0594 }
0595
0596 //=======================================================================
0597 //function : Solves a LU factored Matrix
0598 //purpose :
0599 //=======================================================================
0600
0601 Standard_Integer
0602 BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
0603 const Standard_Integer UpperBandWidth,
0604 const Standard_Integer LowerBandWidth,
0605 Array1OfPoints& PolesArray)
0606 {
0607 Standard_Real *PArray ;
0608 PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
0609
0610 return BSplCLib::SolveBandedSystem(Matrix,
0611 UpperBandWidth,
0612 LowerBandWidth,
0613 Dimension_gen,
0614 PArray[0]) ;
0615 }
0616
0617 //=======================================================================
0618 //function : Solves a LU factored Matrix
0619 //purpose : if HomogeneousFlag is 1 then the input and the output
0620 // will be homogeneous that is no division or multiplication
0621 // by weigths will happen. On the contrary if HomogeneousFlag
0622 // is 0 then the poles will be multiplied first by the weights
0623 // and after interpolation they will be divided by the weights
0624 //=======================================================================
0625
0626 Standard_Integer
0627 BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
0628 const Standard_Integer UpperBandWidth,
0629 const Standard_Integer LowerBandWidth,
0630 const Standard_Boolean HomogeneousFlag,
0631 Array1OfPoints& PolesArray,
0632 TColStd_Array1OfReal& WeightsArray)
0633 {
0634 Standard_Real *PArray,
0635 * WArray ;
0636 PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
0637 WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
0638 return BSplCLib::SolveBandedSystem(Matrix,
0639 UpperBandWidth,
0640 LowerBandWidth,
0641 HomogeneousFlag,
0642 Dimension_gen,
0643 PArray[0],
0644 WArray[0]) ;
0645 }
0646
0647 //=======================================================================
0648 //function : Evaluates a Bspline function : uses the ExtrapMode
0649 //purpose : the function is extrapolated using the Taylor expansion
0650 // of degree ExtrapMode[0] to the left and the Taylor
0651 // expansion of degree ExtrapMode[1] to the right
0652 // if the HomogeneousFlag == True than the Poles are supposed
0653 // to be stored homogeneously and the result will also be homogeneous
0654 // Valid only if Weights
0655 //=======================================================================
0656 void BSplCLib::Eval(const Standard_Real Parameter,
0657 const Standard_Boolean PeriodicFlag,
0658 const Standard_Boolean HomogeneousFlag,
0659 Standard_Integer& ExtrapMode,
0660 const Standard_Integer Degree,
0661 const TColStd_Array1OfReal& FlatKnots,
0662 const Array1OfPoints& PolesArray,
0663 const TColStd_Array1OfReal& WeightsArray,
0664 Point& aPoint,
0665 Standard_Real& aWeight)
0666 {
0667 Standard_Real Inverse,
0668 P[Dimension_gen],
0669 *PArray,
0670 *WArray ;
0671 Standard_Integer kk ;
0672 PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
0673 WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
0674 if (HomogeneousFlag) {
0675 BSplCLib::Eval(Parameter,
0676 PeriodicFlag,
0677 0,
0678 ExtrapMode,
0679 Degree,
0680 FlatKnots,
0681 Dimension_gen,
0682 PArray[0],
0683 P[0]) ;
0684 BSplCLib::Eval(Parameter,
0685 PeriodicFlag,
0686 0,
0687 ExtrapMode,
0688 Degree,
0689 FlatKnots,
0690 1,
0691 WArray[0],
0692 aWeight) ;
0693 }
0694 else {
0695 BSplCLib::Eval(Parameter,
0696 PeriodicFlag,
0697 0,
0698 ExtrapMode,
0699 Degree,
0700 FlatKnots,
0701 Dimension_gen,
0702 PArray[0],
0703 WArray[0],
0704 P[0],
0705 aWeight) ;
0706 Inverse = 1.0e0 / aWeight ;
0707
0708 for (kk = 0 ; kk < Dimension_gen ; kk++) {
0709 P[kk] *= Inverse ;
0710 }
0711 }
0712
0713 for (kk = 0 ; kk < Dimension_gen ; kk++)
0714 aPoint.SetCoord(kk + 1, P[kk]) ;
0715 }
0716
0717 //=======================================================================
0718 //function : CacheD0
0719 //purpose : Evaluates the polynomial cache of the Bspline Curve
0720 //
0721 //=======================================================================
0722 void BSplCLib::CacheD0(const Standard_Real Parameter,
0723 const Standard_Integer Degree,
0724 const Standard_Real CacheParameter,
0725 const Standard_Real SpanLenght,
0726 const Array1OfPoints& PolesArray,
0727 const TColStd_Array1OfReal* WeightsArray,
0728 Point& aPoint)
0729 {
0730 //
0731 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
0732 // form
0733 // the SpanLenght is the normalizing factor so that the polynomial is between
0734 // 0 and 1
0735 Standard_Real NewParameter, Inverse;
0736 Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
0737 Standard_Real * myPoint = (Standard_Real *) &aPoint ;
0738 NewParameter = (Parameter - CacheParameter) / SpanLenght ;
0739 PLib::NoDerivativeEvalPolynomial(NewParameter,
0740 Degree,
0741 Dimension_gen,
0742 Degree * Dimension_gen,
0743 PArray[0],
0744 myPoint[0]) ;
0745 if (WeightsArray != NULL) {
0746 const TColStd_Array1OfReal& refWeights = *WeightsArray;
0747 Standard_Real *
0748 WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
0749 PLib::NoDerivativeEvalPolynomial(NewParameter,
0750 Degree,
0751 1,
0752 Degree,
0753 WArray[0],
0754 Inverse) ;
0755
0756 Inverse = 1.0e0 / Inverse;
0757 ModifyCoords (myPoint, *= Inverse);
0758 }
0759 }
0760
0761 //=======================================================================
0762 //function : CacheD1
0763 //purpose : Evaluates the polynomial cache of the Bspline Curve
0764 //
0765 //=======================================================================
0766 void BSplCLib::CacheD1(const Standard_Real Parameter,
0767 const Standard_Integer Degree,
0768 const Standard_Real CacheParameter,
0769 const Standard_Real SpanLenght,
0770 const Array1OfPoints& PolesArray,
0771 const TColStd_Array1OfReal* WeightsArray,
0772 Point& aPoint,
0773 Vector& aVector)
0774 {
0775 //
0776 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
0777 // form
0778 // the SpanLenght is the normalizing factor so that the polynomial is between
0779 // 0 and 1
0780 Standard_Real LocalPDerivatives[Dimension_gen << 1] ;
0781 // Standard_Real LocalWDerivatives[2], NewParameter, Inverse ;
0782 Standard_Real LocalWDerivatives[2], NewParameter ;
0783
0784 Standard_Real *
0785 PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
0786 Standard_Real *
0787 myPoint = (Standard_Real *) &aPoint ;
0788 Standard_Real *
0789 myVector = (Standard_Real *) &aVector ;
0790 NewParameter = (Parameter - CacheParameter) / SpanLenght ;
0791 PLib::EvalPolynomial(NewParameter,
0792 1,
0793 Degree,
0794 Dimension_gen,
0795 PArray[0],
0796 LocalPDerivatives[0]) ;
0797 //
0798 // unormalize derivatives since those are computed normalized
0799 //
0800
0801 ModifyCoords (LocalPDerivatives + Dimension_gen, /= SpanLenght);
0802
0803 if (WeightsArray != NULL) {
0804 const TColStd_Array1OfReal& refWeights = *WeightsArray;
0805 Standard_Real *
0806 WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
0807 PLib::EvalPolynomial(NewParameter,
0808 1,
0809 Degree,
0810 1,
0811 WArray[0],
0812 LocalWDerivatives[0]) ;
0813 //
0814 // unormalize the result since the polynomial stored in the cache
0815 // is normalized between 0 and 1
0816 //
0817 LocalWDerivatives[1] /= SpanLenght ;
0818
0819 PLib::RationalDerivatives(1,
0820 Dimension_gen,
0821 LocalPDerivatives[0],
0822 LocalWDerivatives[0],
0823 LocalPDerivatives[0]) ;
0824 }
0825
0826 CopyCoords (myPoint, LocalPDerivatives);
0827 CopyCoords (myVector, LocalPDerivatives + Dimension_gen);
0828 }
0829
0830 //=======================================================================
0831 //function : CacheD2
0832 //purpose : Evaluates the polynomial cache of the Bspline Curve
0833 //
0834 //=======================================================================
0835 void BSplCLib::CacheD2(const Standard_Real Parameter,
0836 const Standard_Integer Degree,
0837 const Standard_Real CacheParameter,
0838 const Standard_Real SpanLenght,
0839 const Array1OfPoints& PolesArray,
0840 const TColStd_Array1OfReal* WeightsArray,
0841 Point& aPoint,
0842 Vector& aVector1,
0843 Vector& aVector2)
0844 {
0845 //
0846 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
0847 // form
0848 // the SpanLenght is the normalizing factor so that the polynomial is between
0849 // 0 and 1
0850 Standard_Integer ii,Index,EndIndex;
0851 Standard_Real LocalPDerivatives[(Dimension_gen << 1) + Dimension_gen] ;
0852 // Standard_Real LocalWDerivatives[3], NewParameter, Factor, Inverse ;
0853 Standard_Real LocalWDerivatives[3], NewParameter, Factor ;
0854 Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
0855 Standard_Real * myPoint = (Standard_Real *) &aPoint ;
0856 Standard_Real * myVector1 = (Standard_Real *) &aVector1 ;
0857 Standard_Real * myVector2 = (Standard_Real *) &aVector2 ;
0858 NewParameter = (Parameter - CacheParameter) / SpanLenght ;
0859 PLib::EvalPolynomial(NewParameter,
0860 2,
0861 Degree,
0862 Dimension_gen,
0863 PArray[0],
0864 LocalPDerivatives[0]) ;
0865 //
0866 // unormalize derivatives since those are computed normalized
0867 //
0868 Factor = 1.0e0/ SpanLenght ;
0869 Index = Dimension_gen ;
0870 EndIndex = Min (2, Degree) ;
0871
0872 for (ii = 1 ; ii <= EndIndex ; ii++) {
0873 ModifyCoords (LocalPDerivatives + Index, *= Factor);
0874 Factor /= SpanLenght ;
0875 Index += Dimension_gen;
0876 }
0877
0878 Index = (Degree + 1) * Dimension_gen;
0879 for (ii = Degree ; ii < 2 ; ii++) {
0880 NullifyCoords (LocalPDerivatives + Index);
0881 Index += Dimension_gen;
0882 }
0883
0884 if (WeightsArray != NULL) {
0885 const TColStd_Array1OfReal& refWeights = *WeightsArray;
0886 Standard_Real *
0887 WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
0888
0889 PLib::EvalPolynomial(NewParameter,
0890 2,
0891 Degree,
0892 1,
0893 WArray[0],
0894 LocalWDerivatives[0]) ;
0895
0896 for (ii = Degree + 1 ; ii <= 2 ; ii++) {
0897 LocalWDerivatives[ii] = 0.0e0 ;
0898 }
0899 //
0900 // unormalize the result since the polynomial stored in the cache
0901 // is normalized between 0 and 1
0902 //
0903 Factor = 1.0e0 / SpanLenght ;
0904
0905 for (ii = 1 ; ii <= EndIndex ; ii++) {
0906 LocalWDerivatives[ii] *= Factor ;
0907 Factor /= SpanLenght ;
0908 }
0909 PLib::RationalDerivatives(2,
0910 Dimension_gen,
0911 LocalPDerivatives[0],
0912 LocalWDerivatives[0],
0913 LocalPDerivatives[0]) ;
0914 }
0915
0916 CopyCoords (myPoint, LocalPDerivatives);
0917 CopyCoords (myVector1, LocalPDerivatives + Dimension_gen);
0918 CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2);
0919 }
0920
0921 //=======================================================================
0922 //function : CacheD3
0923 //purpose : Evaluates the polynomial cache of the Bspline Curve
0924 //
0925 //=======================================================================
0926 void BSplCLib::CacheD3(const Standard_Real Parameter,
0927 const Standard_Integer Degree,
0928 const Standard_Real CacheParameter,
0929 const Standard_Real SpanLenght,
0930 const Array1OfPoints& PolesArray,
0931 const TColStd_Array1OfReal* WeightsArray,
0932 Point& aPoint,
0933 Vector& aVector1,
0934 Vector& aVector2,
0935 Vector& aVector3)
0936 {
0937 //
0938 // the CacheParameter is where the cache polynomial was evaluated in homogeneous
0939 // form
0940 // the SpanLenght is the normalizing factor so that the polynomial is between
0941 // 0 and 1
0942 Standard_Integer ii, Index, EndIndex;
0943 Standard_Real LocalPDerivatives[Dimension_gen << 2] ;
0944 // Standard_Real LocalWDerivatives[4], Factor, NewParameter, Inverse ;
0945 Standard_Real LocalWDerivatives[4], Factor, NewParameter ;
0946 Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
0947 Standard_Real * myPoint = (Standard_Real *) &aPoint ;
0948 Standard_Real * myVector1 = (Standard_Real *) &aVector1 ;
0949 Standard_Real * myVector2 = (Standard_Real *) &aVector2 ;
0950 Standard_Real * myVector3 = (Standard_Real *) &aVector3 ;
0951 NewParameter = (Parameter - CacheParameter) / SpanLenght ;
0952 PLib::EvalPolynomial(NewParameter,
0953 3,
0954 Degree,
0955 Dimension_gen,
0956 PArray[0],
0957 LocalPDerivatives[0]) ;
0958
0959 Index = (Degree + 1) * Dimension_gen;
0960 for (ii = Degree ; ii < 3 ; ii++) {
0961 NullifyCoords (LocalPDerivatives + Index);
0962 Index += Dimension_gen;
0963 }
0964
0965 //
0966 // unormalize derivatives since those are computed normalized
0967 //
0968 Factor = 1.0e0 / SpanLenght ;
0969 Index = Dimension_gen ;
0970 EndIndex = Min(3,Degree) ;
0971
0972 for (ii = 1 ; ii <= EndIndex ; ii++) {
0973 ModifyCoords (LocalPDerivatives + Index, *= Factor);
0974 Factor /= SpanLenght;
0975 Index += Dimension_gen;
0976 }
0977
0978 if (WeightsArray != NULL) {
0979 const TColStd_Array1OfReal& refWeights = *WeightsArray;
0980 Standard_Real *
0981 WArray = (Standard_Real *) &refWeights(refWeights.Lower()) ;
0982
0983 PLib::EvalPolynomial(NewParameter,
0984 3,
0985 Degree,
0986 1,
0987 WArray[0],
0988 LocalWDerivatives[0]) ;
0989 //
0990 // unormalize the result since the polynomial stored in the cache
0991 // is normalized between 0 and 1
0992 //
0993 Factor = 1.0e0 / SpanLenght ;
0994
0995 for (ii = 1 ; ii <= EndIndex ; ii++) {
0996 LocalWDerivatives[ii] *= Factor ;
0997 Factor /= SpanLenght ;
0998 }
0999
1000 for (ii = (Degree + 1) ; ii <= 3 ; ii++) {
1001 LocalWDerivatives[ii] = 0.0e0 ;
1002 }
1003 PLib::RationalDerivatives(3,
1004 Dimension_gen,
1005 LocalPDerivatives[0],
1006 LocalWDerivatives[0],
1007 LocalPDerivatives[0]) ;
1008 }
1009
1010 CopyCoords (myPoint, LocalPDerivatives);
1011 CopyCoords (myVector1, LocalPDerivatives + Dimension_gen);
1012 CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2);
1013 CopyCoords (myVector3, LocalPDerivatives + Dimension_gen * 3);
1014 }
1015
1016 //=======================================================================
1017 //function : BuildCache
1018 //purpose : Stores theTaylor expansion normalized between 0,1 in the
1019 // Cache : in case of a rational function the Poles are
1020 // stored in homogeneous form
1021 //=======================================================================
1022
1023 void BSplCLib::BuildCache
1024 (const Standard_Real U,
1025 const Standard_Real SpanDomain,
1026 const Standard_Boolean Periodic,
1027 const Standard_Integer Degree,
1028 const TColStd_Array1OfReal& FlatKnots,
1029 const Array1OfPoints& Poles,
1030 const TColStd_Array1OfReal* Weights,
1031 Array1OfPoints& CachePoles,
1032 TColStd_Array1OfReal* CacheWeights)
1033 {
1034 Standard_Integer ii,
1035 Dimension,
1036 LocalIndex,
1037 index = 0 ;
1038 Standard_Real u = U,
1039 LocalValue ;
1040 Standard_Boolean rational;
1041
1042 BSplCLib_DataContainer dc(Degree);
1043 PrepareEval(u,
1044 index,
1045 Dimension,
1046 rational,
1047 Degree,
1048 Periodic,
1049 Poles,
1050 Weights,
1051 FlatKnots,
1052 (BSplCLib::NoMults()),
1053 dc);
1054 //
1055 // Watch out : PrepareEval checks if locally the Bspline is polynomial
1056 // therefore rational flag could be set to False even if there are
1057 // Weights. The Dimension is set accordingly.
1058 //
1059
1060 BSplCLib::Bohm(u,Degree,Degree,*dc.knots,Dimension,*dc.poles);
1061
1062 LocalValue = 1.0e0 ;
1063 LocalIndex = 0 ;
1064
1065 if (rational) {
1066
1067 for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1068 CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue);
1069 LocalIndex += Dimension_gen + 1;
1070 LocalValue *= SpanDomain / (Standard_Real) ii ;
1071 }
1072
1073 LocalIndex = Dimension_gen;
1074 LocalValue = 1.0e0 ;
1075 for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1076 (*CacheWeights)(ii) = dc.poles[LocalIndex] * LocalValue ;
1077 LocalIndex += Dimension_gen + 1;
1078 LocalValue *= SpanDomain / (Standard_Real) ii ;
1079 }
1080 }
1081 else {
1082
1083 for (ii = 1 ; ii <= Degree + 1 ; ii++) {
1084 CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue);
1085 LocalIndex += Dimension_gen;
1086 LocalValue *= SpanDomain / (Standard_Real) ii ;
1087 }
1088
1089 if (Weights != NULL) {
1090 for (ii = 1 ; ii <= Degree + 1 ; ii++)
1091 (*CacheWeights)(ii) = 0.0e0 ;
1092 (*CacheWeights)(1) = 1.0e0 ;
1093 }
1094 }
1095 }
1096
1097 void BSplCLib::BuildCache(const Standard_Real theParameter,
1098 const Standard_Real theSpanDomain,
1099 const Standard_Boolean thePeriodicFlag,
1100 const Standard_Integer theDegree,
1101 const Standard_Integer theSpanIndex,
1102 const TColStd_Array1OfReal& theFlatKnots,
1103 const Array1OfPoints& thePoles,
1104 const TColStd_Array1OfReal* theWeights,
1105 TColStd_Array2OfReal& theCacheArray)
1106 {
1107 Standard_Real aParam = theParameter;
1108 Standard_Integer anIndex = theSpanIndex;
1109 Standard_Integer aDimension;
1110 Standard_Boolean isRational;
1111
1112 BSplCLib_DataContainer dc(theDegree);
1113 PrepareEval(aParam,
1114 anIndex,
1115 aDimension,
1116 isRational,
1117 theDegree,
1118 thePeriodicFlag,
1119 thePoles,
1120 theWeights,
1121 theFlatKnots,
1122 (BSplCLib::NoMults()),
1123 dc);
1124 //
1125 // Watch out : PrepareEval checks if locally the Bspline is polynomial
1126 // therefore rational flag could be set to False even if there are
1127 // Weights. The Dimension is set accordingly.
1128 //
1129
1130 Standard_Integer aCacheShift = // helps to store weights when PrepareEval did not found that the curve is locally polynomial
1131 (theWeights != NULL && !isRational) ? aDimension + 1 : aDimension;
1132
1133 BSplCLib::Bohm(aParam, theDegree, theDegree, *dc.knots, aDimension, *dc.poles);
1134
1135 Standard_Real aCoeff = 1.0;
1136 Standard_Real* aCache = (Standard_Real *) &(theCacheArray(theCacheArray.LowerRow(), theCacheArray.LowerCol()));
1137 Standard_Real* aPolyCoeffs = dc.poles;
1138
1139 for (Standard_Integer i = 0; i <= theDegree; i++)
1140 {
1141 for (Standard_Integer j = 0; j< aDimension; j++)
1142 aCache[j] = aPolyCoeffs[j] * aCoeff;
1143 aCoeff *= theSpanDomain / (i+1);
1144 aPolyCoeffs += aDimension;
1145 aCache += aDimension;
1146 if (aCacheShift > aDimension)
1147 {
1148 aCache[0] = 0.0;
1149 aCache++;
1150 }
1151 }
1152
1153 if (aCacheShift > aDimension)
1154 theCacheArray.SetValue(theCacheArray.LowerRow(), theCacheArray.LowerCol() + aCacheShift - 1, 1.0);
1155 }
1156
1157 //=======================================================================
1158 //function : Interpolate
1159 //purpose :
1160 //=======================================================================
1161
1162 void BSplCLib::Interpolate(const Standard_Integer Degree,
1163 const TColStd_Array1OfReal& FlatKnots,
1164 const TColStd_Array1OfReal& Parameters,
1165 const TColStd_Array1OfInteger& ContactOrderArray,
1166 Array1OfPoints& Poles,
1167 Standard_Integer& InversionProblem)
1168
1169 {
1170 Standard_Real *PArray ;
1171
1172 PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1173
1174 BSplCLib::Interpolate(Degree,
1175 FlatKnots,
1176 Parameters,
1177 ContactOrderArray,
1178 Dimension_gen,
1179 PArray[0],
1180 InversionProblem) ;
1181 }
1182
1183 //=======================================================================
1184 //function : Interpolate
1185 //purpose :
1186 //=======================================================================
1187
1188 void BSplCLib::Interpolate(const Standard_Integer Degree,
1189 const TColStd_Array1OfReal& FlatKnots,
1190 const TColStd_Array1OfReal& Parameters,
1191 const TColStd_Array1OfInteger& ContactOrderArray,
1192 Array1OfPoints& Poles,
1193 TColStd_Array1OfReal& Weights,
1194 Standard_Integer& InversionProblem)
1195 {
1196 Standard_Real *PArray,
1197 * WArray ;
1198 PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1199 WArray = (Standard_Real *) &Weights(Weights.Lower()) ;
1200 BSplCLib::Interpolate(Degree,
1201 FlatKnots,
1202 Parameters,
1203 ContactOrderArray,
1204 Dimension_gen,
1205 PArray[0],
1206 WArray[0],
1207 InversionProblem) ;
1208 }
1209
1210 //=======================================================================
1211 //function : MovePoint
1212 //purpose : Find the new poles which allows an old point (with a
1213 // given u as parameter) to reach a new position
1214 //=======================================================================
1215
1216 void BSplCLib::MovePoint (const Standard_Real U,
1217 const Vector& Displ,
1218 const Standard_Integer Index1,
1219 const Standard_Integer Index2,
1220 const Standard_Integer Degree,
1221 const Array1OfPoints& Poles,
1222 const TColStd_Array1OfReal* Weights,
1223 const TColStd_Array1OfReal& FlatKnots,
1224 Standard_Integer& FirstIndex,
1225 Standard_Integer& LastIndex,
1226 Array1OfPoints& NewPoles)
1227 {
1228 // calculate the BSplineBasis in the parameter U
1229 Standard_Integer FirstNonZeroBsplineIndex;
1230 math_Matrix BSplineBasis(1, 1,
1231 1, Degree+1);
1232 Standard_Integer ErrorCode =
1233 BSplCLib::EvalBsplineBasis(0,
1234 Degree+1,
1235 FlatKnots,
1236 U,
1237 FirstNonZeroBsplineIndex,
1238 BSplineBasis);
1239 if (ErrorCode != 0) {
1240 FirstIndex = 0;
1241 LastIndex = 0;
1242
1243 for (Standard_Integer i=Poles.Lower(); i<=Poles.Upper(); i++) {
1244 NewPoles(i) = Poles(i);
1245 }
1246 return;
1247 }
1248
1249 // find the span which is predominant for parameter U
1250 FirstIndex = FirstNonZeroBsplineIndex;
1251 LastIndex = FirstNonZeroBsplineIndex + Degree ;
1252 if (FirstIndex < Index1) FirstIndex = Index1;
1253 if (LastIndex > Index2) LastIndex = Index2;
1254
1255 Standard_Real maxValue = 0.0;
1256 Standard_Integer i, kk1=0, kk2, ii;
1257
1258 for (i = FirstIndex-FirstNonZeroBsplineIndex+1;
1259 i <= LastIndex-FirstNonZeroBsplineIndex+1; i++) {
1260 if (BSplineBasis(1,i) > maxValue) {
1261 kk1 = i + FirstNonZeroBsplineIndex - 1;
1262 maxValue = BSplineBasis(1,i);
1263 }
1264 }
1265
1266 // find a kk2 if symetriy
1267 kk2 = kk1;
1268 i = kk1 - FirstNonZeroBsplineIndex + 2;
1269 if ((kk1+1) <= LastIndex) {
1270 if (Abs(BSplineBasis(1, kk1-FirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
1271 kk2 = kk1+1;
1272 }
1273 }
1274
1275 // compute the vector of displacement
1276 Standard_Real D1 = 0.0;
1277 Standard_Real D2 = 0.0;
1278 Standard_Real hN, Coef, Dval;
1279
1280 for (i = 1; i <= Degree+1; i++) {
1281 ii = i + FirstNonZeroBsplineIndex - 1;
1282 if (Weights != NULL) {
1283 hN = Weights->Value(ii)*BSplineBasis(1, i);
1284 D2 += hN;
1285 }
1286 else {
1287 hN = BSplineBasis(1, i);
1288 }
1289 if (ii >= FirstIndex && ii <= LastIndex) {
1290 if (ii < kk1) {
1291 Dval = kk1-ii;
1292 }
1293 else if (ii > kk2) {
1294 Dval = ii - kk2;
1295 }
1296 else {
1297 Dval = 0.0;
1298 }
1299 D1 += 1./(Dval + 1.) * hN;
1300 }
1301 }
1302
1303 if (Weights != NULL) {
1304 Coef = D2/D1;
1305 }
1306 else {
1307 Coef = 1./D1;
1308 }
1309
1310 // compute the new poles
1311
1312 for (i=Poles.Lower(); i<=Poles.Upper(); i++) {
1313 if (i >= FirstIndex && i <= LastIndex) {
1314 if (i < kk1) {
1315 Dval = kk1-i;
1316 }
1317 else if (i > kk2) {
1318 Dval = i - kk2;
1319 }
1320 else {
1321 Dval = 0.0;
1322 }
1323 NewPoles(i) = Poles(i).Translated((Coef/(Dval+1.))*Displ);
1324 }
1325 else {
1326 NewPoles(i) = Poles(i);
1327 }
1328 }
1329 }
1330
1331 //=======================================================================
1332 //function : MovePoint
1333 //purpose : Find the new poles which allows an old point (with a
1334 // given u as parameter) to reach a new position
1335 //=======================================================================
1336
1337 //=======================================================================
1338 //function : MovePointAndTangent
1339 //purpose :
1340 //=======================================================================
1341
1342 void BSplCLib::MovePointAndTangent (const Standard_Real U,
1343 const Vector& Delta,
1344 const Vector& DeltaDerivatives,
1345 const Standard_Real Tolerance,
1346 const Standard_Integer Degree,
1347 const Standard_Integer StartingCondition,
1348 const Standard_Integer EndingCondition,
1349 const Array1OfPoints& Poles,
1350 const TColStd_Array1OfReal* Weights,
1351 const TColStd_Array1OfReal& FlatKnots,
1352 Array1OfPoints& NewPoles,
1353 Standard_Integer & ErrorStatus)
1354 {
1355 Standard_Real *delta_array,
1356 *delta_derivative_array,
1357 *poles_array,
1358 *new_poles_array ;
1359
1360 Standard_Integer num_poles ;
1361 num_poles = Poles.Length() ;
1362
1363 if (NewPoles.Length() != num_poles) {
1364 throw Standard_ConstructionError();
1365 }
1366 delta_array = (Standard_Real *) &Delta ;
1367 delta_derivative_array = (Standard_Real *)&DeltaDerivatives ;
1368 poles_array = (Standard_Real *) &Poles(Poles.Lower()) ;
1369
1370 new_poles_array = (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1371 MovePointAndTangent(U,
1372 Dimension_gen,
1373 delta_array[0],
1374 delta_derivative_array[0],
1375 Tolerance,
1376 Degree,
1377 StartingCondition,
1378 EndingCondition,
1379 poles_array[0],
1380 Weights,
1381 FlatKnots,
1382 new_poles_array[0],
1383 ErrorStatus) ;
1384 }
1385
1386 //=======================================================================
1387 //function : Resolution
1388 //purpose :
1389 //=======================================================================
1390
1391 void BSplCLib::Resolution(const Array1OfPoints& Poles,
1392 const TColStd_Array1OfReal* Weights,
1393 const Standard_Integer NumPoles,
1394 const TColStd_Array1OfReal& FlatKnots,
1395 const Standard_Integer Degree,
1396 const Standard_Real Tolerance3D,
1397 Standard_Real& UTolerance)
1398 {
1399 Standard_Real *PolesArray ;
1400 PolesArray = (Standard_Real *) &Poles(Poles.Lower()) ;
1401 BSplCLib::Resolution(PolesArray[0],
1402 Dimension_gen,
1403 NumPoles,
1404 Weights,
1405 FlatKnots,
1406 Degree,
1407 Tolerance3D,
1408 UTolerance) ;
1409 }
1410
1411 //=======================================================================
1412 //function : FunctionMultiply
1413 //purpose :
1414 //=======================================================================
1415
1416 void BSplCLib::FunctionMultiply
1417 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1418 const Standard_Integer BSplineDegree,
1419 const TColStd_Array1OfReal & BSplineFlatKnots,
1420 const Array1OfPoints & Poles,
1421 const TColStd_Array1OfReal & FlatKnots,
1422 const Standard_Integer NewDegree,
1423 Array1OfPoints & NewPoles,
1424 Standard_Integer & theStatus)
1425 {
1426 Standard_Integer num_bspline_poles =
1427 BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1428 Standard_Integer num_new_poles =
1429 FlatKnots.Length() - NewDegree - 1 ;
1430
1431 if (Poles.Length() != num_bspline_poles ||
1432 NewPoles.Length() != num_new_poles) {
1433 throw Standard_ConstructionError();
1434 }
1435 Standard_Real * array_of_poles =
1436 (Standard_Real *) &Poles(Poles.Lower()) ;
1437 Standard_Real * array_of_new_poles =
1438 (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1439 BSplCLib::FunctionMultiply(FunctionPtr,
1440 BSplineDegree,
1441 BSplineFlatKnots,
1442 Dimension_gen,
1443 array_of_poles[0],
1444 FlatKnots,
1445 NewDegree,
1446 array_of_new_poles[0],
1447 theStatus);
1448 }
1449
1450 //=======================================================================
1451 //function : FunctionReparameterise
1452 //purpose :
1453 //=======================================================================
1454
1455 void BSplCLib::FunctionReparameterise
1456 (const BSplCLib_EvaluatorFunction & FunctionPtr,
1457 const Standard_Integer BSplineDegree,
1458 const TColStd_Array1OfReal & BSplineFlatKnots,
1459 const Array1OfPoints & Poles,
1460 const TColStd_Array1OfReal & FlatKnots,
1461 const Standard_Integer NewDegree,
1462 Array1OfPoints & NewPoles,
1463 Standard_Integer & theStatus)
1464 {
1465 Standard_Integer num_bspline_poles =
1466 BSplineFlatKnots.Length() - BSplineDegree - 1 ;
1467 Standard_Integer num_new_poles =
1468 FlatKnots.Length() - NewDegree - 1 ;
1469
1470 if (Poles.Length() != num_bspline_poles ||
1471 NewPoles.Length() != num_new_poles) {
1472 throw Standard_ConstructionError();
1473 }
1474 Standard_Real * array_of_poles =
1475 (Standard_Real *) &Poles(Poles.Lower()) ;
1476 Standard_Real * array_of_new_poles =
1477 (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
1478 BSplCLib::FunctionReparameterise(FunctionPtr,
1479 BSplineDegree,
1480 BSplineFlatKnots,
1481 Dimension_gen,
1482 array_of_poles[0],
1483 FlatKnots,
1484 NewDegree,
1485 array_of_new_poles[0],
1486 theStatus);
1487 }
