// File: BSplCLib_3.cxx
// Created: Fri May 19 18:33:07 1995
// Author: Xavier BENVENISTE
// <xab@nonox>
#include <BSplCLib.ixx>
#include <Standard_NotImplemented.hxx>
// BSpline Curve in 3d space
// ***************************
#define No_Standard_RangeError
#define No_Standard_OutOfRange
#define Dimension_gen 3
#define Array1OfPoints TColgp_Array1OfPnt
#define Point gp_Pnt
#define Vector gp_Vec
#define PointToCoords(carr,pnt,op) \
(carr)[0] = (pnt).X() op, \
(carr)[1] = (pnt).Y() op, \
(carr)[2] = (pnt).Z() op
#define CoordsToPoint(pnt,carr,op) \
(pnt).SetX ((carr)[0] op), \
(pnt).SetY ((carr)[1] op), \
(pnt).SetZ ((carr)[2] op)
#define NullifyPoint(pnt) \
(pnt).SetCoord (0.,0.,0.)
#define NullifyCoords(carr) \
(carr)[0] = (carr)[1] = (carr)[2] = 0.
#define ModifyCoords(carr,op) \
(carr)[0] op, \
(carr)[1] op, \
(carr)[2] op
#define CopyCoords(carr,carr2) \
(carr)[0] = (carr2)[0], \
(carr)[1] = (carr2)[1], \
(carr)[2] = (carr2)[2]
#define BSplCLib_DataContainer BSplCLib_DataContainer_3d
#include <BSplCLib_CurveComputation.gxx>
// xab : modified 15-Mar-94 added EvalBSplineMatrix, FactorBandedMatrix, SolveBandedSystem
// Eval, BuildCache, cacheD0, cacheD1, cacheD2, cacheD3, LocalMatrix,
// EvalPolynomial, RationalDerivatives.
// xab : 22-Mar-94 : fixed rational problem in BuildCache
// xab : 12-Mar-96 : added MovePointAndTangent
//
#include <TColStd_Array1OfInteger.hxx>
#include <TColStd_Array1OfReal.hxx>
#include <gp_Vec2d.hxx>
#include <Standard_ConstructionError.hxx>
#include <PLib.hxx>
#include <math_Matrix.hxx>
#define No_Standard_RangeError
#define No_Standard_OutOfRange
//=======================================================================
//struct : BSplCLib_DataContainer
//purpose: Auxiliary structure providing buffers for poles and knots used in
// evaluation of bspline (allocated in the stack)
//=======================================================================
struct BSplCLib_DataContainer
{
BSplCLib_DataContainer(Standard_Integer Degree)
{
if ( Degree > BSplCLib::MaxDegree() || BSplCLib::MaxDegree() > 25 )
Standard_OutOfRange::Raise ("BSplCLib: bspline degree is greater than maximum supported");
}
Standard_Real poles[(25+1)*(Dimension_gen+1)];
Standard_Real knots[2*25];
Standard_Real ders[Dimension_gen*4];
};
//=======================================================================
//function : Reverse
//purpose :
//=======================================================================
void BSplCLib::Reverse(Array1OfPoints& Poles,
const Standard_Integer L)
{
Standard_Integer i, l = L;
l = Poles.Lower() + (l-Poles.Lower()) % (Poles.Upper()-Poles.Lower()+1);
Array1OfPoints temp(0,Poles.Length()-1);
for (i = Poles.Lower(); i <= l; i++)
temp(l-i) = Poles(i);
for (i = l+1; i <= Poles.Upper(); i++)
temp(l-Poles.Lower()+Poles.Upper()-i+1) = Poles(i);
for (i = Poles.Lower(); i <= Poles.Upper(); i++)
Poles(i) = temp(i-Poles.Lower());
}
//
// CURVES MODIFICATIONS
//
//=======================================================================
//function : RemoveKnot
//purpose :
//=======================================================================
Standard_Boolean BSplCLib::RemoveKnot
(const Standard_Integer Index,
const Standard_Integer Mult,
const Standard_Integer Degree,
const Standard_Boolean Periodic,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
const TColStd_Array1OfReal& Knots,
const TColStd_Array1OfInteger& Mults,
Array1OfPoints& NewPoles,
TColStd_Array1OfReal& NewWeights,
TColStd_Array1OfReal& NewKnots,
TColStd_Array1OfInteger& NewMults,
const Standard_Real Tolerance)
{
Standard_Boolean rational = &Weights != NULL;
Standard_Integer dim;
dim = Dimension_gen;
if (rational) dim++;
TColStd_Array1OfReal poles(1,dim*Poles.Length());
TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
if (rational) PLib::SetPoles(Poles,Weights,poles);
else PLib::SetPoles(Poles,poles);
if (!RemoveKnot(Index,Mult,Degree,Periodic,dim,
poles,Knots,Mults,newpoles,NewKnots,NewMults,Tolerance))
return Standard_False;
if (rational) PLib::GetPoles(newpoles,NewPoles,NewWeights);
else PLib::GetPoles(newpoles,NewPoles);
return Standard_True;
}
//=======================================================================
//function : InsertKnots
//purpose : insert an array of knots and multiplicities
//=======================================================================
void BSplCLib::InsertKnots
(const Standard_Integer Degree,
const Standard_Boolean Periodic,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
const TColStd_Array1OfReal& Knots,
const TColStd_Array1OfInteger& Mults,
const TColStd_Array1OfReal& AddKnots,
const TColStd_Array1OfInteger& AddMults,
Array1OfPoints& NewPoles,
TColStd_Array1OfReal& NewWeights,
TColStd_Array1OfReal& NewKnots,
TColStd_Array1OfInteger& NewMults,
const Standard_Real Epsilon,
const Standard_Boolean Add)
{
Standard_Boolean rational = &Weights != NULL;
Standard_Integer dim;
dim = Dimension_gen;
if (rational) dim++;
TColStd_Array1OfReal poles(1,dim*Poles.Length());
TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
if (rational) PLib::SetPoles(Poles,Weights,poles);
else PLib::SetPoles(Poles,poles);
InsertKnots(Degree,Periodic,dim,poles,Knots,Mults,
AddKnots,AddMults,newpoles,NewKnots,NewMults,Epsilon,Add);
if (rational) PLib::GetPoles(newpoles,NewPoles,NewWeights);
else PLib::GetPoles(newpoles,NewPoles);
}
//=======================================================================
//function : InsertKnot
//purpose :
//=======================================================================
void BSplCLib::InsertKnot(const Standard_Integer ,
const Standard_Real U,
const Standard_Integer UMult,
const Standard_Integer Degree,
const Standard_Boolean Periodic,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
const TColStd_Array1OfReal& Knots,
const TColStd_Array1OfInteger& Mults,
Array1OfPoints& NewPoles,
TColStd_Array1OfReal& NewWeights)
{
TColStd_Array1OfReal k(1,1);
k(1) = U;
TColStd_Array1OfInteger m(1,1);
m(1) = UMult;
TColStd_Array1OfReal nk(1,Knots.Length()+1);
TColStd_Array1OfInteger nm(1,Knots.Length()+1);
InsertKnots(Degree,Periodic,Poles,Weights,Knots,Mults,
k,m,NewPoles,NewWeights,nk,nm,Epsilon(U));
}
//=======================================================================
//function : RaiseMultiplicity
//purpose :
//=======================================================================
void BSplCLib::RaiseMultiplicity(const Standard_Integer KnotIndex,
const Standard_Integer Mult,
const Standard_Integer Degree,
const Standard_Boolean Periodic,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
const TColStd_Array1OfReal& Knots,
const TColStd_Array1OfInteger& Mults,
Array1OfPoints& NewPoles,
TColStd_Array1OfReal& NewWeights)
{
TColStd_Array1OfReal k(1,1);
k(1) = Knots(KnotIndex);
TColStd_Array1OfInteger m(1,1);
m(1) = Mult - Mults(KnotIndex);
TColStd_Array1OfReal nk(1,Knots.Length());
TColStd_Array1OfInteger nm(1,Knots.Length());
InsertKnots(Degree,Periodic,Poles,Weights,Knots,Mults,
k,m,NewPoles,NewWeights,nk,nm,Epsilon(k(1)));
}
//=======================================================================
//function : IncreaseDegree
//purpose :
//=======================================================================
void BSplCLib::IncreaseDegree
(const Standard_Integer Degree,
const Standard_Integer NewDegree,
const Standard_Boolean Periodic,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
const TColStd_Array1OfReal& Knots,
const TColStd_Array1OfInteger& Mults,
Array1OfPoints& NewPoles,
TColStd_Array1OfReal& NewWeights,
TColStd_Array1OfReal& NewKnots,
TColStd_Array1OfInteger& NewMults)
{
Standard_Boolean rational = &Weights != NULL;
Standard_Integer dim;
dim = Dimension_gen;
if (rational) dim++;
TColStd_Array1OfReal poles(1,dim*Poles.Length());
TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
if (rational) PLib::SetPoles(Poles,Weights,poles);
else PLib::SetPoles(Poles,poles);
IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults,
newpoles,NewKnots,NewMults);
if (rational) PLib::GetPoles(newpoles,NewPoles,NewWeights);
else PLib::GetPoles(newpoles,NewPoles);
}
//=======================================================================
//function : Unperiodize
//purpose :
//=======================================================================
void BSplCLib::Unperiodize
(const Standard_Integer Degree,
const TColStd_Array1OfInteger& Mults,
const TColStd_Array1OfReal& Knots,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
TColStd_Array1OfInteger& NewMults,
TColStd_Array1OfReal& NewKnots,
Array1OfPoints& NewPoles,
TColStd_Array1OfReal& NewWeights)
{
Standard_Boolean rational = &Weights != NULL;
Standard_Integer dim;
dim = Dimension_gen;
if (rational) dim++;
TColStd_Array1OfReal poles(1,dim*Poles.Length());
TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
if (rational) PLib::SetPoles(Poles,Weights,poles);
else PLib::SetPoles(Poles,poles);
Unperiodize(Degree, dim, Mults, Knots, poles,
NewMults,NewKnots, newpoles);
if (rational) PLib::GetPoles(newpoles,NewPoles,NewWeights);
else PLib::GetPoles(newpoles,NewPoles);
}
//=======================================================================
//function : Trimming
//purpose :
//=======================================================================
void BSplCLib::Trimming(const Standard_Integer Degree,
const Standard_Boolean Periodic,
const TColStd_Array1OfReal& Knots,
const TColStd_Array1OfInteger& Mults,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
const Standard_Real U1,
const Standard_Real U2,
TColStd_Array1OfReal& NewKnots,
TColStd_Array1OfInteger& NewMults,
Array1OfPoints& NewPoles,
TColStd_Array1OfReal& NewWeights)
{
Standard_Boolean rational = &Weights != NULL;
Standard_Integer dim;
dim = Dimension_gen;
if (rational) dim++;
TColStd_Array1OfReal poles(1,dim*Poles.Length());
TColStd_Array1OfReal newpoles(1,dim*NewPoles.Length());
if (rational) PLib::SetPoles(Poles,Weights,poles);
else PLib::SetPoles(Poles,poles);
Trimming(Degree, Periodic, dim, Knots, Mults, poles, U1, U2,
NewKnots, NewMults, newpoles);
if (rational) PLib::GetPoles(newpoles,NewPoles,NewWeights);
else PLib::GetPoles(newpoles,NewPoles);
}
//--------------------------------------------------------------------------
//ELEMENTARY COMPUTATIONS
//--------------------------------------------------------------------------
//
// All the following methods are methods of low level used in BSplCLib
// but also in BSplSLib (but not in Geom)
// At the creation time of this package there is no possibility
// to declare private methods of package and to declare friend
// methods of package. It could interesting to declare that BSplSLib
// is a package friend to BSplCLib because it uses all the basis methods
// of BSplCLib.
//--------------------------------------------------------------------------
//=======================================================================
//function : BuildEval
//purpose : builds the local array for evaluation
//=======================================================================
void BSplCLib::BuildEval(const Standard_Integer Degree,
const Standard_Integer Index,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
Standard_Real& LP)
{
Standard_Real w, *pole = &LP;
Standard_Integer PLower = Poles.Lower();
Standard_Integer PUpper = Poles.Upper();
Standard_Integer i;
Standard_Integer ip = PLower + Index - 1;
if (&Weights == NULL) {
for (i = 0; i <= Degree; i++) {
ip++;
if (ip > PUpper) ip = PLower;
const Point& P = Poles(ip);
PointToCoords (pole,P,+0);
pole += Dimension_gen;
}
}
else {
for (i = 0; i <= Degree; i++) {
ip++;
if (ip > PUpper) ip = PLower;
const Point& P = Poles(ip);
pole[Dimension_gen] = w = Weights(ip);
PointToCoords (pole, P, * w);
pole += Dimension_gen + 1;
}
}
}
//=======================================================================
//function : PrepareEval
//purpose : stores data for Eval in the local arrays
// dc.poles and dc.knots
//=======================================================================
static void PrepareEval
(Standard_Real& u,
Standard_Integer& index,
Standard_Integer& dim,
Standard_Boolean& rational,
const Standard_Integer Degree,
const Standard_Boolean Periodic,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
const TColStd_Array1OfReal& Knots,
const TColStd_Array1OfInteger& Mults,
BSplCLib_DataContainer& dc)
{
// Set the Index
BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u);
// make the knots
BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*dc.knots);
if (&Mults == NULL)
index -= Knots.Lower() + Degree;
else
index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults);
// check truly rational
rational = (&Weights != NULL);
if (rational) {
Standard_Integer WLower = Weights.Lower() + index;
rational = BSplCLib::IsRational(Weights, WLower, WLower + Degree);
}
// make the poles
if (rational) {
dim = Dimension_gen + 1;
BSplCLib::BuildEval(Degree, index, Poles, Weights , *dc.poles);
}
else {
dim = Dimension_gen;
BSplCLib::BuildEval(Degree, index, Poles, BSplCLib::NoWeights(), *dc.poles);
}
}
//=======================================================================
//function : D0
//purpose :
//=======================================================================
void BSplCLib::D0
(const Standard_Real U,
const Standard_Integer Index,
const Standard_Integer Degree,
const Standard_Boolean Periodic,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
const TColStd_Array1OfReal& Knots,
const TColStd_Array1OfInteger& Mults,
Point& P)
{
// Standard_Integer k,dim,index = Index;
Standard_Integer dim,index = Index;
Standard_Real u = U;
Standard_Boolean rational;
BSplCLib_DataContainer dc(Degree);
PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
BSplCLib::Eval(u,Degree,*dc.knots,dim,*dc.poles);
if (rational) {
Standard_Real w = dc.poles[Dimension_gen];
CoordsToPoint (P, dc.poles, / w);
}
else
CoordsToPoint (P, dc.poles, );
}
//=======================================================================
//function : D1
//purpose :
//=======================================================================
void BSplCLib::D1
(const Standard_Real U,
const Standard_Integer Index,
const Standard_Integer Degree,
const Standard_Boolean Periodic,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
const TColStd_Array1OfReal& Knots,
const TColStd_Array1OfInteger& Mults,
Point& P,
Vector& V)
{
Standard_Integer dim,index = Index;
Standard_Real u = U;
Standard_Boolean rational;
BSplCLib_DataContainer dc(Degree);
PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
BSplCLib::Bohm(u,Degree,1,*dc.knots,dim,*dc.poles);
Standard_Real *result = dc.poles;
if (rational) {
PLib::RationalDerivative(Degree,1,Dimension_gen,*dc.poles,*dc.ders);
result = dc.ders;
}
CoordsToPoint (P, result, );
CoordsToPoint (V, result + Dimension_gen, );
}
//=======================================================================
//function : D2
//purpose :
//=======================================================================
void BSplCLib::D2
(const Standard_Real U,
const Standard_Integer Index,
const Standard_Integer Degree,
const Standard_Boolean Periodic,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
const TColStd_Array1OfReal& Knots,
const TColStd_Array1OfInteger& Mults,
Point& P,
Vector& V1,
Vector& V2)
{
Standard_Integer dim,index = Index;
Standard_Real u = U;
Standard_Boolean rational;
BSplCLib_DataContainer dc(Degree);
PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
BSplCLib::Bohm(u,Degree,2,*dc.knots,dim,*dc.poles);
Standard_Real *result = dc.poles;
if (rational) {
PLib::RationalDerivative(Degree,2,Dimension_gen,*dc.poles,*dc.ders);
result = dc.ders;
}
CoordsToPoint (P, result, );
CoordsToPoint (V1, result + Dimension_gen, );
if (!rational && (Degree < 2))
NullifyPoint (V2);
else
CoordsToPoint (V2, result + 2 * Dimension_gen, );
}
//=======================================================================
//function : D3
//purpose :
//=======================================================================
void BSplCLib::D3
(const Standard_Real U,
const Standard_Integer Index,
const Standard_Integer Degree,
const Standard_Boolean Periodic,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
const TColStd_Array1OfReal& Knots,
const TColStd_Array1OfInteger& Mults,
Point& P,
Vector& V1,
Vector& V2,
Vector& V3)
{
Standard_Integer dim,index = Index;
Standard_Real u = U;
Standard_Boolean rational;
BSplCLib_DataContainer dc(Degree);
PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
BSplCLib::Bohm(u,Degree,3,*dc.knots,dim,*dc.poles);
Standard_Real *result = dc.poles;
if (rational) {
PLib::RationalDerivative(Degree,3,Dimension_gen,*dc.poles,*dc.ders);
result = dc.ders;
}
CoordsToPoint (P, result, );
CoordsToPoint (V1, result + Dimension_gen, );
if (!rational && (Degree < 2))
NullifyPoint (V2);
else
CoordsToPoint (V2, result + 2 * Dimension_gen, );
if (!rational && (Degree < 3))
NullifyPoint (V3);
else
CoordsToPoint (V3, result + 3 * Dimension_gen, );
}
//=======================================================================
//function : DN
//purpose :
//=======================================================================
void BSplCLib::DN
(const Standard_Real U,
const Standard_Integer N,
const Standard_Integer Index,
const Standard_Integer Degree,
const Standard_Boolean Periodic,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
const TColStd_Array1OfReal& Knots,
const TColStd_Array1OfInteger& Mults,
Vector& VN)
{
Standard_Integer dim,index = Index;
Standard_Real u = U;
Standard_Boolean rational;
BSplCLib_DataContainer dc(Degree);
PrepareEval(u,index,dim,rational,Degree,Periodic,Poles,Weights,Knots,Mults,dc);
BSplCLib::Bohm(u,Degree,N,*dc.knots,dim,*dc.poles);
if (rational) {
Standard_Real v[Dimension_gen];
PLib::RationalDerivative(Degree,N,Dimension_gen,*dc.poles,v[0],Standard_False);
CoordsToPoint (VN, v, );
}
else {
if (N > Degree)
NullifyPoint (VN);
else {
Standard_Real *DN = dc.poles + N * Dimension_gen;
CoordsToPoint (VN, DN, );
}
}
}
//=======================================================================
//function : Solves a LU factored Matrix
//purpose :
//=======================================================================
Standard_Integer
BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
const Standard_Integer UpperBandWidth,
const Standard_Integer LowerBandWidth,
Array1OfPoints& PolesArray)
{
Standard_Real *PArray ;
PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
return BSplCLib::SolveBandedSystem(Matrix,
UpperBandWidth,
LowerBandWidth,
Dimension_gen,
PArray[0]) ;
}
//=======================================================================
//function : Solves a LU factored Matrix
//purpose : if HomogeneousFlag is 1 then the input and the output
// will be homogeneous that is no division or multiplication
// by weigths will happen. On the contrary if HomogeneousFlag
// is 0 then the poles will be multiplied first by the weights
// and after interpolation they will be devided by the weights
//=======================================================================
Standard_Integer
BSplCLib::SolveBandedSystem(const math_Matrix& Matrix,
const Standard_Integer UpperBandWidth,
const Standard_Integer LowerBandWidth,
const Standard_Boolean HomogeneousFlag,
Array1OfPoints& PolesArray,
TColStd_Array1OfReal& WeightsArray)
{
Standard_Real *PArray,
* WArray ;
PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
return BSplCLib::SolveBandedSystem(Matrix,
UpperBandWidth,
LowerBandWidth,
HomogeneousFlag,
Dimension_gen,
PArray[0],
WArray[0]) ;
}
//=======================================================================
//function : Evaluates a Bspline function : uses the ExtrapMode
//purpose : the function is extrapolated using the Taylor expansion
// of degree ExtrapMode[0] to the left and the Taylor
// expansion of degree ExtrapMode[1] to the right
// if the HomogeneousFlag == True than the Poles are supposed
// to be stored homogeneously and the result will also be homogeneous
// Valid only if Weights
//=======================================================================
void BSplCLib::Eval(const Standard_Real Parameter,
const Standard_Boolean PeriodicFlag,
const Standard_Boolean HomogeneousFlag,
Standard_Integer& ExtrapMode,
const Standard_Integer Degree,
const TColStd_Array1OfReal& FlatKnots,
const Array1OfPoints& PolesArray,
const TColStd_Array1OfReal& WeightsArray,
Point& aPoint,
Standard_Real& aWeight)
{
Standard_Real Inverse,
P[Dimension_gen],
*PArray,
*WArray ;
Standard_Integer kk ;
PArray = (Standard_Real *) &PolesArray(PolesArray.Lower()) ;
WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
if (HomogeneousFlag) {
BSplCLib::Eval(Parameter,
PeriodicFlag,
0,
ExtrapMode,
Degree,
FlatKnots,
Dimension_gen,
PArray[0],
P[0]) ;
BSplCLib::Eval(Parameter,
PeriodicFlag,
0,
ExtrapMode,
Degree,
FlatKnots,
1,
WArray[0],
aWeight) ;
}
else {
BSplCLib::Eval(Parameter,
PeriodicFlag,
0,
ExtrapMode,
Degree,
FlatKnots,
Dimension_gen,
PArray[0],
WArray[0],
P[0],
aWeight) ;
Inverse = 1.0e0 / aWeight ;
for (kk = 0 ; kk < Dimension_gen ; kk++) {
P[kk] *= Inverse ;
}
}
for (kk = 0 ; kk < Dimension_gen ; kk++)
aPoint.SetCoord(kk + 1, P[kk]) ;
}
//=======================================================================
//function : CacheD0
//purpose : Evaluates the polynomial cache of the Bspline Curve
//
//=======================================================================
void BSplCLib::CacheD0(const Standard_Real Parameter,
const Standard_Integer Degree,
const Standard_Real CacheParameter,
const Standard_Real SpanLenght,
const Array1OfPoints& PolesArray,
const TColStd_Array1OfReal& WeightsArray,
Point& aPoint)
{
//
// the CacheParameter is where the cache polynomial was evaluated in homogeneous
// form
// the SpanLenght is the normalizing factor so that the polynomial is between
// 0 and 1
Standard_Real NewParameter, Inverse;
Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
Standard_Real * myPoint = (Standard_Real *) &aPoint ;
NewParameter = (Parameter - CacheParameter) / SpanLenght ;
PLib::NoDerivativeEvalPolynomial(NewParameter,
Degree,
Dimension_gen,
Degree * Dimension_gen,
PArray[0],
myPoint[0]) ;
if (&WeightsArray != NULL) {
Standard_Real *
WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
PLib::NoDerivativeEvalPolynomial(NewParameter,
Degree,
1,
Degree,
WArray[0],
Inverse) ;
Inverse = 1.0e0 / Inverse;
ModifyCoords (myPoint, *= Inverse);
}
}
//=======================================================================
//function : CacheD1
//purpose : Evaluates the polynomial cache of the Bspline Curve
//
//=======================================================================
void BSplCLib::CacheD1(const Standard_Real Parameter,
const Standard_Integer Degree,
const Standard_Real CacheParameter,
const Standard_Real SpanLenght,
const Array1OfPoints& PolesArray,
const TColStd_Array1OfReal& WeightsArray,
Point& aPoint,
Vector& aVector)
{
//
// the CacheParameter is where the cache polynomial was evaluated in homogeneous
// form
// the SpanLenght is the normalizing factor so that the polynomial is between
// 0 and 1
Standard_Real LocalPDerivatives[Dimension_gen << 1] ;
// Standard_Real LocalWDerivatives[2], NewParameter, Inverse ;
Standard_Real LocalWDerivatives[2], NewParameter ;
Standard_Real *
PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
Standard_Real *
myPoint = (Standard_Real *) &aPoint ;
Standard_Real *
myVector = (Standard_Real *) &aVector ;
NewParameter = (Parameter - CacheParameter) / SpanLenght ;
PLib::EvalPolynomial(NewParameter,
1,
Degree,
Dimension_gen,
PArray[0],
LocalPDerivatives[0]) ;
//
// unormalize derivatives since those are computed normalized
//
ModifyCoords (LocalPDerivatives + Dimension_gen, /= SpanLenght);
if (&WeightsArray != NULL) {
Standard_Real *
WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
PLib::EvalPolynomial(NewParameter,
1,
Degree,
1,
WArray[0],
LocalWDerivatives[0]) ;
//
// unormalize the result since the polynomial stored in the cache
// is normalized between 0 and 1
//
LocalWDerivatives[1] /= SpanLenght ;
PLib::RationalDerivatives(1,
Dimension_gen,
LocalPDerivatives[0],
LocalWDerivatives[0],
LocalPDerivatives[0]) ;
}
CopyCoords (myPoint, LocalPDerivatives);
CopyCoords (myVector, LocalPDerivatives + Dimension_gen);
}
//=======================================================================
//function : CacheD2
//purpose : Evaluates the polynomial cache of the Bspline Curve
//
//=======================================================================
void BSplCLib::CacheD2(const Standard_Real Parameter,
const Standard_Integer Degree,
const Standard_Real CacheParameter,
const Standard_Real SpanLenght,
const Array1OfPoints& PolesArray,
const TColStd_Array1OfReal& WeightsArray,
Point& aPoint,
Vector& aVector1,
Vector& aVector2)
{
//
// the CacheParameter is where the cache polynomial was evaluated in homogeneous
// form
// the SpanLenght is the normalizing factor so that the polynomial is between
// 0 and 1
Standard_Integer ii,Index,EndIndex;
Standard_Real LocalPDerivatives[(Dimension_gen << 1) + Dimension_gen] ;
// Standard_Real LocalWDerivatives[3], NewParameter, Factor, Inverse ;
Standard_Real LocalWDerivatives[3], NewParameter, Factor ;
Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
Standard_Real * myPoint = (Standard_Real *) &aPoint ;
Standard_Real * myVector1 = (Standard_Real *) &aVector1 ;
Standard_Real * myVector2 = (Standard_Real *) &aVector2 ;
NewParameter = (Parameter - CacheParameter) / SpanLenght ;
PLib::EvalPolynomial(NewParameter,
2,
Degree,
Dimension_gen,
PArray[0],
LocalPDerivatives[0]) ;
//
// unormalize derivatives since those are computed normalized
//
Factor = 1.0e0/ SpanLenght ;
Index = Dimension_gen ;
EndIndex = Min (2, Degree) ;
for (ii = 1 ; ii <= EndIndex ; ii++) {
ModifyCoords (LocalPDerivatives + Index, *= Factor);
Factor /= SpanLenght ;
Index += Dimension_gen;
}
Index = (Degree + 1) * Dimension_gen;
for (ii = Degree ; ii < 2 ; ii++) {
NullifyCoords (LocalPDerivatives + Index);
Index += Dimension_gen;
}
if (&WeightsArray != NULL) {
Standard_Real *
WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
PLib::EvalPolynomial(NewParameter,
2,
Degree,
1,
WArray[0],
LocalWDerivatives[0]) ;
for (ii = Degree + 1 ; ii <= 2 ; ii++) {
LocalWDerivatives[ii] = 0.0e0 ;
}
//
// unormalize the result since the polynomial stored in the cache
// is normalized between 0 and 1
//
Factor = 1.0e0 / SpanLenght ;
for (ii = 1 ; ii <= EndIndex ; ii++) {
LocalWDerivatives[ii] *= Factor ;
Factor /= SpanLenght ;
}
PLib::RationalDerivatives(2,
Dimension_gen,
LocalPDerivatives[0],
LocalWDerivatives[0],
LocalPDerivatives[0]) ;
}
CopyCoords (myPoint, LocalPDerivatives);
CopyCoords (myVector1, LocalPDerivatives + Dimension_gen);
CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2);
}
//=======================================================================
//function : CacheD3
//purpose : Evaluates the polynomial cache of the Bspline Curve
//
//=======================================================================
void BSplCLib::CacheD3(const Standard_Real Parameter,
const Standard_Integer Degree,
const Standard_Real CacheParameter,
const Standard_Real SpanLenght,
const Array1OfPoints& PolesArray,
const TColStd_Array1OfReal& WeightsArray,
Point& aPoint,
Vector& aVector1,
Vector& aVector2,
Vector& aVector3)
{
//
// the CacheParameter is where the cache polynomial was evaluated in homogeneous
// form
// the SpanLenght is the normalizing factor so that the polynomial is between
// 0 and 1
Standard_Integer ii, Index, EndIndex;
Standard_Real LocalPDerivatives[Dimension_gen << 2] ;
// Standard_Real LocalWDerivatives[4], Factor, NewParameter, Inverse ;
Standard_Real LocalWDerivatives[4], Factor, NewParameter ;
Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.Lower())) ;
Standard_Real * myPoint = (Standard_Real *) &aPoint ;
Standard_Real * myVector1 = (Standard_Real *) &aVector1 ;
Standard_Real * myVector2 = (Standard_Real *) &aVector2 ;
Standard_Real * myVector3 = (Standard_Real *) &aVector3 ;
NewParameter = (Parameter - CacheParameter) / SpanLenght ;
PLib::EvalPolynomial(NewParameter,
3,
Degree,
Dimension_gen,
PArray[0],
LocalPDerivatives[0]) ;
Index = (Degree + 1) * Dimension_gen;
for (ii = Degree ; ii < 3 ; ii++) {
NullifyCoords (LocalPDerivatives + Index);
Index += Dimension_gen;
}
//
// unormalize derivatives since those are computed normalized
//
Factor = 1.0e0 / SpanLenght ;
Index = Dimension_gen ;
EndIndex = Min(3,Degree) ;
for (ii = 1 ; ii <= EndIndex ; ii++) {
ModifyCoords (LocalPDerivatives + Index, *= Factor);
Factor /= SpanLenght;
Index += Dimension_gen;
}
if (&WeightsArray != NULL) {
Standard_Real *
WArray = (Standard_Real *) &WeightsArray(WeightsArray.Lower()) ;
PLib::EvalPolynomial(NewParameter,
3,
Degree,
1,
WArray[0],
LocalWDerivatives[0]) ;
//
// unormalize the result since the polynomial stored in the cache
// is normalized between 0 and 1
//
Factor = 1.0e0 / SpanLenght ;
for (ii = 1 ; ii <= EndIndex ; ii++) {
LocalWDerivatives[ii] *= Factor ;
Factor /= SpanLenght ;
}
for (ii = (Degree + 1) ; ii <= 3 ; ii++) {
LocalWDerivatives[ii] = 0.0e0 ;
}
PLib::RationalDerivatives(3,
Dimension_gen,
LocalPDerivatives[0],
LocalWDerivatives[0],
LocalPDerivatives[0]) ;
}
CopyCoords (myPoint, LocalPDerivatives);
CopyCoords (myVector1, LocalPDerivatives + Dimension_gen);
CopyCoords (myVector2, LocalPDerivatives + Dimension_gen * 2);
CopyCoords (myVector3, LocalPDerivatives + Dimension_gen * 3);
}
//=======================================================================
//function : BuildCache
//purpose : Stores theTaylor expansion normalized between 0,1 in the
// Cache : in case of a rational function the Poles are
// stored in homogeneous form
//=======================================================================
void BSplCLib::BuildCache
(const Standard_Real U,
const Standard_Real SpanDomain,
const Standard_Boolean Periodic,
const Standard_Integer Degree,
const TColStd_Array1OfReal& FlatKnots,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
Array1OfPoints& CachePoles,
TColStd_Array1OfReal& CacheWeights)
{
Standard_Integer ii,
Dimension,
LocalIndex,
index = 0 ;
Standard_Real u = U,
LocalValue ;
Standard_Boolean rational;
BSplCLib_DataContainer dc(Degree);
PrepareEval(u,
index,
Dimension,
rational,
Degree,
Periodic,
Poles,
Weights,
FlatKnots,
(BSplCLib::NoMults()),
dc);
//
// Watch out : PrepareEval checks if locally the Bspline is polynomial
// therefore rational flag could be set to False even if there are
// Weights. The Dimension is set accordingly.
//
BSplCLib::Bohm(u,Degree,Degree,*dc.knots,Dimension,*dc.poles);
LocalValue = 1.0e0 ;
LocalIndex = 0 ;
if (rational) {
for (ii = 1 ; ii <= Degree + 1 ; ii++) {
CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue);
LocalIndex += Dimension_gen + 1;
LocalValue *= SpanDomain / (Standard_Real) ii ;
}
LocalIndex = Dimension_gen;
LocalValue = 1.0e0 ;
for (ii = 1 ; ii <= Degree + 1 ; ii++) {
CacheWeights(ii) = dc.poles[LocalIndex] * LocalValue ;
LocalIndex += Dimension_gen + 1;
LocalValue *= SpanDomain / (Standard_Real) ii ;
}
}
else {
for (ii = 1 ; ii <= Degree + 1 ; ii++) {
CoordsToPoint (CachePoles(ii), dc.poles + LocalIndex, * LocalValue);
LocalIndex += Dimension_gen;
LocalValue *= SpanDomain / (Standard_Real) ii ;
}
if (&Weights != NULL) {
for (ii = 1 ; ii <= Degree + 1 ; ii++)
CacheWeights(ii) = 0.0e0 ;
CacheWeights(1) = 1.0e0 ;
}
}
}
//=======================================================================
//function : Interpolate
//purpose :
//=======================================================================
void BSplCLib::Interpolate(const Standard_Integer Degree,
const TColStd_Array1OfReal& FlatKnots,
const TColStd_Array1OfReal& Parameters,
const TColStd_Array1OfInteger& ContactOrderArray,
Array1OfPoints& Poles,
Standard_Integer& InversionProblem)
{
Standard_Real *PArray ;
PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
BSplCLib::Interpolate(Degree,
FlatKnots,
Parameters,
ContactOrderArray,
Dimension_gen,
PArray[0],
InversionProblem) ;
}
//=======================================================================
//function : Interpolate
//purpose :
//=======================================================================
void BSplCLib::Interpolate(const Standard_Integer Degree,
const TColStd_Array1OfReal& FlatKnots,
const TColStd_Array1OfReal& Parameters,
const TColStd_Array1OfInteger& ContactOrderArray,
Array1OfPoints& Poles,
TColStd_Array1OfReal& Weights,
Standard_Integer& InversionProblem)
{
Standard_Real *PArray,
* WArray ;
PArray = (Standard_Real *) &Poles(Poles.Lower()) ;
WArray = (Standard_Real *) &Weights(Weights.Lower()) ;
BSplCLib::Interpolate(Degree,
FlatKnots,
Parameters,
ContactOrderArray,
Dimension_gen,
PArray[0],
WArray[0],
InversionProblem) ;
}
//=======================================================================
//function : MovePoint
//purpose : Find the new poles which allows an old point (with a
// given u as parameter) to reach a new position
//=======================================================================
void BSplCLib::MovePoint (const Standard_Real U,
const Vector& Displ,
const Standard_Integer Index1,
const Standard_Integer Index2,
const Standard_Integer Degree,
const Standard_Boolean Rational,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
const TColStd_Array1OfReal& FlatKnots,
Standard_Integer& FirstIndex,
Standard_Integer& LastIndex,
Array1OfPoints& NewPoles)
{
// calculate the BSplineBasis in the parameter U
Standard_Integer FirstNonZeroBsplineIndex;
math_Matrix BSplineBasis(1, 1,
1, Degree+1);
Standard_Integer ErrorCode =
BSplCLib::EvalBsplineBasis(1,
0,
Degree+1,
FlatKnots,
U,
FirstNonZeroBsplineIndex,
BSplineBasis);
if (ErrorCode != 0) {
FirstIndex = 0;
LastIndex = 0;
for (Standard_Integer i=Poles.Lower(); i<=Poles.Upper(); i++) {
NewPoles(i) = Poles(i);
}
return;
}
// find the span which is predominant for parameter U
FirstIndex = FirstNonZeroBsplineIndex;
LastIndex = FirstNonZeroBsplineIndex + Degree ;
if (FirstIndex < Index1) FirstIndex = Index1;
if (LastIndex > Index2) LastIndex = Index2;
Standard_Real maxValue = 0.0;
Standard_Integer i, kk1=0, kk2, ii;
for (i = FirstIndex-FirstNonZeroBsplineIndex+1;
i <= LastIndex-FirstNonZeroBsplineIndex+1; i++) {
if (BSplineBasis(1,i) > maxValue) {
kk1 = i + FirstNonZeroBsplineIndex - 1;
maxValue = BSplineBasis(1,i);
}
}
// find a kk2 if symetriy
kk2 = kk1;
i = kk1 - FirstNonZeroBsplineIndex + 2;
if ((kk1+1) <= LastIndex) {
if (Abs(BSplineBasis(1, kk1-FirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) {
kk2 = kk1+1;
}
}
// compute the vector of displacement
Standard_Real D1 = 0.0;
Standard_Real D2 = 0.0;
Standard_Real hN, Coef, Dval;
for (i = 1; i <= Degree+1; i++) {
ii = i + FirstNonZeroBsplineIndex - 1;
if (Rational) {
hN = Weights(ii)*BSplineBasis(1, i);
D2 += hN;
}
else {
hN = BSplineBasis(1, i);
}
if (ii >= FirstIndex && ii <= LastIndex) {
if (ii < kk1) {
Dval = kk1-ii;
}
else if (ii > kk2) {
Dval = ii - kk2;
}
else {
Dval = 0.0;
}
D1 += 1./(Dval + 1.) * hN;
}
}
if (Rational) {
Coef = D2/D1;
}
else {
Coef = 1./D1;
}
// compute the new poles
for (i=Poles.Lower(); i<=Poles.Upper(); i++) {
if (i >= FirstIndex && i <= LastIndex) {
if (i < kk1) {
Dval = kk1-i;
}
else if (i > kk2) {
Dval = i - kk2;
}
else {
Dval = 0.0;
}
NewPoles(i) = Poles(i).Translated((Coef/(Dval+1.))*Displ);
}
else {
NewPoles(i) = Poles(i);
}
}
}
//=======================================================================
//function : MovePoint
//purpose : Find the new poles which allows an old point (with a
// given u as parameter) to reach a new position
//=======================================================================
//=======================================================================
//function : MovePointAndTangent
//purpose :
//=======================================================================
void BSplCLib::MovePointAndTangent (const Standard_Real U,
const Vector& Delta,
const Vector& DeltaDerivatives,
const Standard_Real Tolerance,
const Standard_Integer Degree,
const Standard_Boolean Rational,
const Standard_Integer StartingCondition,
const Standard_Integer EndingCondition,
const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
const TColStd_Array1OfReal& FlatKnots,
Array1OfPoints& NewPoles,
Standard_Integer & ErrorStatus)
{
Standard_Real *delta_array,
*delta_derivative_array,
*poles_array,
*new_poles_array ;
Standard_Integer num_poles ;
num_poles = Poles.Length() ;
if (NewPoles.Length() != num_poles) {
Standard_ConstructionError::Raise();
}
delta_array = (Standard_Real *) &Delta ;
delta_derivative_array = (Standard_Real *)&DeltaDerivatives ;
poles_array = (Standard_Real *) &Poles(Poles.Lower()) ;
new_poles_array = (Standard_Real *) &NewPoles(NewPoles.Lower()) ;
MovePointAndTangent(U,
Dimension_gen,
delta_array[0],
delta_derivative_array[0],
Tolerance,
Degree,
Rational,
StartingCondition,
EndingCondition,
poles_array[0],
Weights,
FlatKnots,
new_poles_array[0],
ErrorStatus) ;
}
//=======================================================================
//function : Resolution
//purpose :
//=======================================================================
void BSplCLib::Resolution(const Array1OfPoints& Poles,
const TColStd_Array1OfReal& Weights,
const Standard_Integer NumPoles,
const TColStd_Array1OfReal& FlatKnots,
const Standard_Integer Degree,
const Standard_Real Tolerance3D,
Standard_Real& UTolerance)
{
Standard_Real *PolesArray ;
PolesArray = (Standard_Real *) &Poles(Poles.Lower()) ;
BSplCLib::Resolution(PolesArray[0],
Dimension_gen,
NumPoles,
Weights,
FlatKnots,
Degree,
Tolerance3D,
UTolerance) ;
}
//=======================================================================
//function : FunctionMultiply
//purpose :
//=======================================================================
void BSplCLib::FunctionMultiply
(const BSplCLib_EvaluatorFunction & FunctionPtr,
const Standard_Integer BSplineDegree,
const TColStd_Array1OfReal & BSplineFlatKnots,
const Array1OfPoints & Poles,
const TColStd_Array1OfReal & FlatKnots,
const Standard_Integer NewDegree,
Array1OfPoints & NewPoles,
Standard_Integer & Status)
{
Standard_Integer num_bspline_poles =
BSplineFlatKnots.Length() - BSplineDegree - 1 ;
Standard_Integer num_new_poles =
FlatKnots.Length() - NewDegree - 1 ;
if (Poles.Length() != num_bspline_poles ||
NewPoles.Length() != num_new_poles) {
Standard_ConstructionError::Raise();
}
Standard_Real * array_of_poles =
(Standard_Real *) &Poles(Poles.Lower()) ;
Standard_Real * array_of_new_poles =
(Standard_Real *) &NewPoles(NewPoles.Lower()) ;
BSplCLib::FunctionMultiply(FunctionPtr,
BSplineDegree,
BSplineFlatKnots,
Dimension_gen,
array_of_poles[0],
FlatKnots,
NewDegree,
array_of_new_poles[0],
Status) ;
}
//=======================================================================
//function : FunctionReparameterise
//purpose :
//=======================================================================
void BSplCLib::FunctionReparameterise
(const BSplCLib_EvaluatorFunction & FunctionPtr,
const Standard_Integer BSplineDegree,
const TColStd_Array1OfReal & BSplineFlatKnots,
const Array1OfPoints & Poles,
const TColStd_Array1OfReal & FlatKnots,
const Standard_Integer NewDegree,
Array1OfPoints & NewPoles,
Standard_Integer & Status)
{
Standard_Integer num_bspline_poles =
BSplineFlatKnots.Length() - BSplineDegree - 1 ;
Standard_Integer num_new_poles =
FlatKnots.Length() - NewDegree - 1 ;
if (Poles.Length() != num_bspline_poles ||
NewPoles.Length() != num_new_poles) {
Standard_ConstructionError::Raise();
}
Standard_Real * array_of_poles =
(Standard_Real *) &Poles(Poles.Lower()) ;
Standard_Real * array_of_new_poles =
(Standard_Real *) &NewPoles(NewPoles.Lower()) ;
BSplCLib::FunctionReparameterise(FunctionPtr,
BSplineDegree,
BSplineFlatKnots,
Dimension_gen,
array_of_poles[0],
FlatKnots,
NewDegree,
array_of_new_poles[0],
Status) ;
}