// Created on: 1991-08-26 // Created by: JCV // Copyright (c) 1991-1999 Matra Datavision // Copyright (c) 1999-2014 OPEN CASCADE SAS // // This file is part of Open CASCADE Technology software library. // // This library is free software; you can redistribute it and/or modify it under // the terms of the GNU Lesser General Public License version 2.1 as published // by the Free Software Foundation, with special exception defined in the file // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT // distribution for complete text of the license and disclaimer of any warranty. // // Alternatively, this file may be used under the terms of Open CASCADE // commercial license or contractual agreement. // Modifed RLE Aug 93 - Complete rewrite // xab 21-Mar-95 implemented cache mecanism // pmn 25-09-96 Interpolation // jct 25-09-96 : Correction de l'alloc de LocalArray dans RationalDerivative. // pmn 07-10-96 : Correction de DN dans le cas rationnal. // pmn 06-02-97 : Correction des poids dans RationalDerivative. (PRO700) #include #include #include #include #include #include #include #include #include #include // for null derivatives static Standard_Real BSplSLib_zero[3] = {0.0, 0.0, 0.0}; //======================================================================= //struct : BSplCLib_DataContainer //purpose: Auxiliary structure providing buffers for poles and knots used in // evaluation of bspline (allocated in the stack) //======================================================================= struct BSplSLib_DataContainer { BSplSLib_DataContainer (Standard_Integer UDegree, Standard_Integer VDegree) { (void)UDegree; (void)VDegree; // just to avoid compiler warning in Release mode Standard_OutOfRange_Raise_if (UDegree > BSplCLib::MaxDegree() || VDegree > BSplCLib::MaxDegree() || BSplCLib::MaxDegree() > 25, "BSplSLib: bspline degree is greater than maximum supported"); } Standard_Real poles[4*(25+1)*(25+1)]; Standard_Real knots1[2*25]; Standard_Real knots2[2*25]; Standard_Real ders[48]; }; //************************************************************************** // Evaluation methods //************************************************************************** //======================================================================= //function : RationalDerivative //purpose : computes the rational derivatives when whe have the // the derivatives of the homogeneous numerator and the // the derivatives of the denominator //======================================================================= void BSplSLib::RationalDerivative(const Standard_Integer UDeg, const Standard_Integer VDeg, const Standard_Integer N, const Standard_Integer M, Standard_Real& HDerivatives, Standard_Real& RDerivatives, const Standard_Boolean All) { // // if All is True all derivatives are computed. if Not only // the requested N, M is computed // // Numerator(u,v) // let f(u,v) be a rational function = ------------------ // Denominator(u,v) // // // Let (N,M) the order of the derivatives we want : then since // we have : // // Numerator = f * Denominator // // we derive : // // (N,M) 1 ( (N M) (p q) (N -p M-q) ) // f = ------------ ( Numerator - SUM SUM a * f * Denominator ) // (0,0) ( p StoreDerivatives (All ? 0 : ii * 3); Standard_Real *RArray = (All ? &RDerivatives : (Standard_Real*)StoreDerivatives); NCollection_LocalArray StoreW (ii); Standard_Real *HomogeneousArray = &HDerivatives; Standard_Real denominator,Pii,Pip,Pjq; denominator = 1.0e0 / HomogeneousArray[3]; index_u = 0; index_u1 = 0; if (UDeg < N) MinN = UDeg; else MinN = N; if (VDeg < M) MinM = VDeg; else MinM = M; MinN1 = MinN + 1; MinM1 = MinM + 1; iiM1 = - M1; for (ii = 0 ; ii < MinN1 ; ii++) { iiM1 += M1; index_v = index_u; index_v1 = index_u1; index_w = iiM1; for (jj = 0 ; jj < MinM1 ; jj++) { RArray[index_v++] = HomogeneousArray[index_v1++]; RArray[index_v++] = HomogeneousArray[index_v1++]; RArray[index_v++] = HomogeneousArray[index_v1++]; StoreW[index_w++] = HomogeneousArray[index_v1++]; } for (jj = MinM1 ; jj < M1 ; jj++) { RArray[index_v++] = 0.; RArray[index_v++] = 0.; RArray[index_v++] = 0.; StoreW[index_w++] = 0.; } index_u1 += M4; index_u += M3; } index_v = MinN1 * M3; index_w = MinN1 * M1; for (ii = MinN1 ; ii < N1 ; ii++) { for (jj = 0 ; jj < M1 ; jj++) { RArray[index_v++] = 0.0e0; RArray[index_v++] = 0.0e0; RArray[index_v++] = 0.0e0; StoreW[index_w++] = 0.0e0; } } // --------------- Calculation ---------------- iiM1 = - M1; iiM3 = - M3; for (ii = 0 ; ii <= N ; ii++) { iiM1 += M1; iiM3 += M3; index1 = iiM3 - 3; jjM1 = iiM1; for (jj = 0 ; jj <= M ; jj++) { jjM1 ++; ppM1 = - M1; ppM3 = - M3; index1 += 3; for (pp = 0 ; pp < ii ; pp++) { ppM1 += M1; ppM3 += M3; index = ppM3; index2 = jjM1 - ppM1; Pip = PLib::Bin(ii,pp); for (qq = 0 ; qq <= jj ; qq++) { index2--; Pjq = Pip * PLib::Bin(jj,qq) * StoreW[index2]; RArray[index1] -= Pjq * RArray[index]; index++; index1++; RArray[index1] -= Pjq * RArray[index]; index++; index1++; RArray[index1] -= Pjq * RArray[index]; index++; index1 -= 2; } } index = iiM3; index2 = jj + 1; Pii = PLib::Bin(ii,ii); for (qq = 0 ; qq < jj ; qq++) { index2--; Pjq = Pii * PLib::Bin(jj,qq) * StoreW[index2]; RArray[index1] -= Pjq * RArray[index]; index++; index1++; RArray[index1] -= Pjq * RArray[index]; index++; index1++; RArray[index1] -= Pjq * RArray[index]; index++; index1 -= 2; } RArray[index1] *= denominator; index1++; RArray[index1] *= denominator; index1++; RArray[index1] *= denominator; index1 -= 2; } } if (!All) { RArray = &RDerivatives; index = N * M1 + M; index = (index << 1) + index; RArray[0] = StoreDerivatives[index]; index++; RArray[1] = StoreDerivatives[index]; index++; RArray[2] = StoreDerivatives[index]; } } //======================================================================= //function : PrepareEval //purpose : //======================================================================= // // PrepareEval : // // Prepare all data for computing points : // local arrays of knots // local array of poles (multiplied by the weights if rational) // // The first direction to compute (smaller degree) is returned // and the poles are stored according to this direction. static Standard_Boolean PrepareEval (const Standard_Real U, const Standard_Real V, const Standard_Integer Uindex, const Standard_Integer Vindex, const Standard_Integer UDegree, const Standard_Integer VDegree, const Standard_Boolean URat, const Standard_Boolean VRat, const Standard_Boolean UPer, const Standard_Boolean VPer, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, const TColStd_Array1OfReal& UKnots, const TColStd_Array1OfReal& VKnots, const TColStd_Array1OfInteger& UMults, const TColStd_Array1OfInteger& VMults, Standard_Real& u1, // first parameter to use Standard_Real& u2, // second parameter to use Standard_Integer& d1, // first degree Standard_Integer& d2, // second degree Standard_Boolean& rational, BSplSLib_DataContainer& dc) { rational = URat || VRat; Standard_Integer uindex = Uindex; Standard_Integer vindex = Vindex; Standard_Integer UKLower = UKnots.Lower(); Standard_Integer UKUpper = UKnots.Upper(); Standard_Integer VKLower = VKnots.Lower(); Standard_Integer VKUpper = VKnots.Upper(); if (UDegree <= VDegree) { // compute the indices if (uindex < UKLower || uindex > UKUpper) BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u1); else u1 = U; if (vindex < VKLower || vindex > VKUpper) BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u2); else u2 = V; // get the knots d1 = UDegree; d2 = VDegree; BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots1); BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots2); if (&UMults == NULL) uindex -= UKLower + UDegree; else uindex = BSplCLib::PoleIndex(UDegree,uindex,UPer,UMults); if (&VMults == NULL) vindex -= VKLower + VDegree; else vindex = BSplCLib::PoleIndex(VDegree,vindex,VPer,VMults); // get the poles Standard_Integer i,j,ip,jp; Standard_Real w, *pole = dc.poles; d1 = UDegree; d2 = VDegree; Standard_Integer PLowerRow = Poles.LowerRow(); Standard_Integer PUpperRow = Poles.UpperRow(); Standard_Integer PLowerCol = Poles.LowerCol(); Standard_Integer PUpperCol = Poles.UpperCol(); // verify if locally non rational if (rational) { rational = Standard_False; ip = PLowerRow + uindex; jp = PLowerCol + vindex; if(ip < PLowerRow) ip = PUpperRow; if(jp < PLowerCol) jp = PUpperCol; w = Weights.Value(ip,jp); Standard_Real eps = Epsilon(w); Standard_Real dw; for (i = 0; i <= UDegree && !rational; i++) { jp = PLowerCol + vindex; if(jp < PLowerCol) jp = PUpperCol; for (j = 0; j <= VDegree && !rational; j++) { dw = Weights.Value(ip,jp) - w; if (dw < 0) dw = - dw; rational = (dw > eps); jp++; if (jp > PUpperCol) jp = PLowerCol; } ip++; if (ip > PUpperRow) ip = PLowerRow; } } // copy the poles ip = PLowerRow + uindex; if(ip < PLowerRow) ip = PUpperRow; if (rational) { for (i = 0; i <= d1; i++) { jp = PLowerCol + vindex; if(jp < PLowerCol) jp = PUpperCol; for (j = 0; j <= d2; j++) { const gp_Pnt& P = Poles .Value(ip,jp); pole[3] = w = Weights.Value(ip,jp); pole[0] = P.X() * w; pole[1] = P.Y() * w; pole[2] = P.Z() * w; pole += 4; jp++; if (jp > PUpperCol) jp = PLowerCol; } ip++; if (ip > PUpperRow) ip = PLowerRow; } } else { for (i = 0; i <= d1; i++) { jp = PLowerCol + vindex; if(jp < PLowerCol) jp = PUpperCol; for (j = 0; j <= d2; j++) { const gp_Pnt& P = Poles.Value(ip,jp); pole[0] = P.X(); pole[1] = P.Y(); pole[2] = P.Z(); pole += 3; jp++; if (jp > PUpperCol) jp = PLowerCol; } ip++; if (ip > PUpperRow) ip = PLowerRow; } } return Standard_True; } else { // compute the indices if (uindex < UKLower || uindex > UKUpper) BSplCLib::LocateParameter(UDegree,UKnots,UMults,U,UPer,uindex,u2); else u2 = U; if (vindex < VKLower || vindex > VKUpper) BSplCLib::LocateParameter(VDegree,VKnots,VMults,V,VPer,vindex,u1); else u1 = V; // get the knots d2 = UDegree; d1 = VDegree; BSplCLib::BuildKnots(UDegree,uindex,UPer,UKnots,UMults,*dc.knots2); BSplCLib::BuildKnots(VDegree,vindex,VPer,VKnots,VMults,*dc.knots1); if (&UMults == NULL) uindex -= UKLower + UDegree; else uindex = BSplCLib::PoleIndex(UDegree,uindex,UPer,UMults); if (&VMults == NULL) vindex -= VKLower + VDegree; else vindex = BSplCLib::PoleIndex(VDegree,vindex,VPer,VMults); // get the poles Standard_Integer i,j,ip,jp; Standard_Real w, *pole = dc.poles; d1 = VDegree; d2 = UDegree; Standard_Integer PLowerRow = Poles.LowerRow(); Standard_Integer PUpperRow = Poles.UpperRow(); Standard_Integer PLowerCol = Poles.LowerCol(); Standard_Integer PUpperCol = Poles.UpperCol(); // verify if locally non rational if (rational) { rational = Standard_False; ip = PLowerRow + uindex; jp = PLowerCol + vindex; if(ip < PLowerRow) ip = PUpperRow; if(jp < PLowerCol) jp = PUpperCol; w = Weights.Value(ip,jp); Standard_Real eps = Epsilon(w); Standard_Real dw; for (i = 0; i <= UDegree && !rational; i++) { jp = PLowerCol + vindex; if(jp < PLowerCol) jp = PUpperCol; for (j = 0; j <= VDegree && !rational; j++) { dw = Weights.Value(ip,jp) - w; if (dw < 0) dw = - dw; rational = dw > eps; jp++; if (jp > PUpperCol) jp = PLowerCol; } ip++; if (ip > PUpperRow) ip = PLowerRow; } } // copy the poles jp = PLowerCol + vindex; if(jp < PLowerCol) jp = PUpperCol; if (rational) { for (i = 0; i <= d1; i++) { ip = PLowerRow + uindex; if(ip < PLowerRow) ip = PUpperRow; for (j = 0; j <= d2; j++) { const gp_Pnt& P = Poles.Value(ip,jp); pole[3] = w = Weights.Value(ip,jp); pole[0] = P.X() * w; pole[1] = P.Y() * w; pole[2] = P.Z() * w; pole += 4; ip++; if (ip > PUpperRow) ip = PLowerRow; } jp++; if (jp > PUpperCol) jp = PLowerCol; } } else { for (i = 0; i <= d1; i++) { ip = PLowerRow + uindex; if(ip < PLowerRow) ip = PUpperRow; for (j = 0; j <= d2; j++) { const gp_Pnt& P = Poles.Value(ip,jp); pole[0] = P.X(); pole[1] = P.Y(); pole[2] = P.Z(); pole += 3; ip++; if (ip > PUpperRow) ip = PLowerRow; } jp++; if (jp > PUpperCol) jp = PLowerCol; } } return Standard_False; } } //======================================================================= //function : D0 //purpose : //======================================================================= void BSplSLib::D0 (const Standard_Real U, const Standard_Real V, const Standard_Integer UIndex, const Standard_Integer VIndex, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, const TColStd_Array1OfReal& UKnots, const TColStd_Array1OfReal& VKnots, const TColStd_Array1OfInteger& UMults, const TColStd_Array1OfInteger& VMults, const Standard_Integer UDegree, const Standard_Integer VDegree, const Standard_Boolean URat, const Standard_Boolean VRat, const Standard_Boolean UPer, const Standard_Boolean VPer, gp_Pnt& P) { // Standard_Integer k ; Standard_Real W ; HomogeneousD0(U, V, UIndex, VIndex, Poles, Weights, UKnots, VKnots, UMults, VMults, UDegree, VDegree, URat, VRat, UPer, VPer, W, P) ; P.SetX(P.X() / W); P.SetY(P.Y() / W); P.SetZ(P.Z() / W); } //======================================================================= //function : D0 //purpose : //======================================================================= void BSplSLib::HomogeneousD0 (const Standard_Real U, const Standard_Real V, const Standard_Integer UIndex, const Standard_Integer VIndex, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, const TColStd_Array1OfReal& UKnots, const TColStd_Array1OfReal& VKnots, const TColStd_Array1OfInteger& UMults, const TColStd_Array1OfInteger& VMults, const Standard_Integer UDegree, const Standard_Integer VDegree, const Standard_Boolean URat, const Standard_Boolean VRat, const Standard_Boolean UPer, const Standard_Boolean VPer, Standard_Real & W, gp_Pnt& P) { Standard_Boolean rational; // Standard_Integer k,dim; Standard_Integer dim; Standard_Real u1,u2; Standard_Integer d1,d2; W = 1.0e0 ; BSplSLib_DataContainer dc (UDegree, VDegree); PrepareEval(U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer, Poles,Weights,UKnots,VKnots,UMults,VMults, u1,u2,d1,d2,rational,dc); if (rational) { dim = 4; BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles); BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles); W = dc.poles[3]; P.SetX(dc.poles[0]); P.SetY(dc.poles[1]); P.SetZ(dc.poles[2]); } else { dim = 3; BSplCLib::Eval(u1,d1,*dc.knots1,dim * (d2 + 1),*dc.poles); BSplCLib::Eval(u2,d2,*dc.knots2,dim,*dc.poles); P.SetX(dc.poles[0]); P.SetY(dc.poles[1]); P.SetZ(dc.poles[2]); } } //======================================================================= //function : D1 //purpose : //======================================================================= void BSplSLib::D1 (const Standard_Real U, const Standard_Real V, const Standard_Integer UIndex, const Standard_Integer VIndex, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, const TColStd_Array1OfReal& UKnots, const TColStd_Array1OfReal& VKnots, const TColStd_Array1OfInteger& UMults, const TColStd_Array1OfInteger& VMults, const Standard_Integer UDegree, const Standard_Integer VDegree, const Standard_Boolean URat, const Standard_Boolean VRat, const Standard_Boolean UPer, const Standard_Boolean VPer, gp_Pnt& P, gp_Vec& Vu, gp_Vec& Vv) { Standard_Boolean rational; // Standard_Integer k,dim,dim2; Standard_Integer dim,dim2; Standard_Real u1,u2; Standard_Integer d1,d2; Standard_Real *result, *resVu, *resVv; BSplSLib_DataContainer dc (UDegree, VDegree); if (PrepareEval (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer, Poles,Weights,UKnots,VKnots,UMults,VMults, u1,u2,d1,d2,rational,dc)) { if (rational) { dim = 4; dim2 = (d2 + 1) << 2; BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles); BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles); BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2)); BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders); result = dc.ders; resVu = result + 6; resVv = result + 3; } else { dim = 3; dim2 = d2 + 1; dim2 = (dim2 << 1) + dim2; BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles); BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles); BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2)); result = dc.poles; resVu = result + dim2; resVv = result + 3; } } else { if (rational) { dim = 4; dim2 = (d2 + 1) << 2; BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles); BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles); BSplCLib::Eval(u2,d2, *dc.knots2,dim ,*(dc.poles + dim2)); BSplSLib::RationalDerivative(d1,d2,1,1,*dc.poles,*dc.ders); result = dc.ders; resVu = result + 3; resVv = result + 6; } else { dim = 3; dim2 = d2 + 1; dim2 = (dim2 << 1) + dim2; BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim2,*dc.poles); BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*dc.poles); BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + dim2)); result = dc.poles; resVu = result + 3; resVv = result + dim2; } } P .SetX(result[0]); Vu.SetX(resVu [0]); Vv.SetX(resVv [0]); P .SetY(result[1]); Vu.SetY(resVu [1]); Vv.SetY(resVv [1]); P .SetZ(result[2]); Vu.SetZ(resVu [2]); Vv.SetZ(resVv [2]); } //======================================================================= //function : D1 //purpose : //======================================================================= void BSplSLib::HomogeneousD1 (const Standard_Real U, const Standard_Real V, const Standard_Integer UIndex, const Standard_Integer VIndex, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, const TColStd_Array1OfReal& UKnots, const TColStd_Array1OfReal& VKnots, const TColStd_Array1OfInteger& UMults, const TColStd_Array1OfInteger& VMults, const Standard_Integer UDegree, const Standard_Integer VDegree, const Standard_Boolean URat, const Standard_Boolean VRat, const Standard_Boolean UPer, const Standard_Boolean VPer, gp_Pnt& N, gp_Vec& Nu, gp_Vec& Nv, Standard_Real& D, Standard_Real& Du, Standard_Real& Dv) { Standard_Boolean rational; // Standard_Integer k,dim; Standard_Integer dim; Standard_Real u1,u2; Standard_Integer d1,d2; D = 1.0e0 ; Du = 0.0e0 ; Dv = 0.0e0 ; BSplSLib_DataContainer dc (UDegree, VDegree); Standard_Boolean ufirst = PrepareEval (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer, Poles,Weights,UKnots,VKnots,UMults,VMults, u1,u2,d1,d2,rational,dc); dim = rational ? 4 : 3; BSplCLib::Bohm(u1,d1,1,*dc.knots1,dim * (d2 + 1),*dc.poles); BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim,*dc.poles); BSplCLib::Eval(u2,d2,*dc.knots2,dim,*(dc.poles+dim*(d2+1))); Standard_Real *result, *resVu, *resVv; result = dc.poles; resVu = result + (ufirst ? dim*(d2+1) : dim); resVv = result + (ufirst ? dim : dim*(d2+1)); N .SetX(result[0]); Nu.SetX(resVu [0]); Nv.SetX(resVv [0]); N .SetY(result[1]); Nu.SetY(resVu [1]); Nv.SetY(resVv [1]); N .SetZ(result[2]); Nu.SetZ(resVu [2]); Nv.SetZ(resVv [2]); if (rational) { D = result[3]; Du = resVu [3]; Dv = resVv [3]; } } //======================================================================= //function : D2 //purpose : //======================================================================= void BSplSLib::D2 (const Standard_Real U, const Standard_Real V, const Standard_Integer UIndex, const Standard_Integer VIndex, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, const TColStd_Array1OfReal& UKnots, const TColStd_Array1OfReal& VKnots, const TColStd_Array1OfInteger& UMults, const TColStd_Array1OfInteger& VMults, const Standard_Integer UDegree, const Standard_Integer VDegree, const Standard_Boolean URat, const Standard_Boolean VRat, const Standard_Boolean UPer, const Standard_Boolean VPer, gp_Pnt& P, gp_Vec& Vu, gp_Vec& Vv, gp_Vec& Vuu, gp_Vec& Vvv, gp_Vec& Vuv) { Standard_Boolean rational; // Standard_Integer k,dim,dim2; Standard_Integer dim,dim2; Standard_Real u1,u2; Standard_Integer d1,d2; Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv; BSplSLib_DataContainer dc (UDegree, VDegree); if (PrepareEval (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer, Poles,Weights,UKnots,VKnots,UMults,VMults, u1,u2,d1,d2,rational,dc)) { if (rational) { dim = 4; dim2 = (d2 + 1) << 2; BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles); BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles); BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2)); if (d1 > 1) BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1))); BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders); result = dc.ders; resVu = result + 9; resVv = result + 3; resVuu = result + 18; resVvv = result + 6; resVuv = result + 12; } else { dim = 3; dim2 = d2 + 1; dim2 = (dim2 << 1) + dim2; BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles); BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles); BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2)); if (d1 > 1) BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1))); result = dc.poles; resVu = result + dim2; resVv = result + 3; if (UDegree <= 1) resVuu = BSplSLib_zero; else resVuu = result + (dim2 << 1); if (VDegree <= 1) resVvv = BSplSLib_zero; else resVvv = result + 6; resVuv = result + (d2 << 1) + d2 + 6; } } else { if (rational) { dim = 4; dim2 = (d2 + 1) << 2; BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles); BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles); BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2)); if (d1 > 1) BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1))); BSplSLib::RationalDerivative(d1,d2,2,2,*dc.poles,*dc.ders); result = dc.ders; resVu = result + 3; resVv = result + 9; resVuu = result + 6; resVvv = result + 18; resVuv = result + 12; } else { dim = 3; dim2 = d2 + 1; dim2 = (dim2 << 1) + dim2; BSplCLib::Bohm(u1,d1,2,*dc.knots1,dim2,*dc.poles); BSplCLib::Bohm(u2,d2,2,*dc.knots2,dim ,*dc.poles); BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + dim2)); if (d1 > 1) BSplCLib::Eval(u2,d2,*dc.knots2,dim ,*(dc.poles + (dim2 << 1))); result = dc.poles; resVu = result + 3; resVv = result + dim2; if (UDegree <= 1) resVuu = BSplSLib_zero; else resVuu = result + 6; if (VDegree <= 1) resVvv = BSplSLib_zero; else resVvv = result + (dim2 << 1); resVuv = result + (d2 << 1) + d2 + 6; } } P .SetX(result[0]); Vu .SetX(resVu [0]); Vv .SetX(resVv [0]); Vuu.SetX(resVuu[0]); Vvv.SetX(resVvv[0]); Vuv.SetX(resVuv[0]); P .SetY(result[1]); Vu .SetY(resVu [1]); Vv .SetY(resVv [1]); Vuu.SetY(resVuu[1]); Vvv.SetY(resVvv[1]); Vuv.SetY(resVuv[1]); P .SetZ(result[2]); Vu .SetZ(resVu [2]); Vv .SetZ(resVv [2]); Vuu.SetZ(resVuu[2]); Vvv.SetZ(resVvv[2]); Vuv.SetZ(resVuv[2]); } //======================================================================= //function : D3 //purpose : //======================================================================= void BSplSLib::D3 (const Standard_Real U, const Standard_Real V, const Standard_Integer UIndex, const Standard_Integer VIndex, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, const TColStd_Array1OfReal& UKnots, const TColStd_Array1OfReal& VKnots, const TColStd_Array1OfInteger& UMults, const TColStd_Array1OfInteger& VMults, const Standard_Integer UDegree, const Standard_Integer VDegree, const Standard_Boolean URat, const Standard_Boolean VRat, const Standard_Boolean UPer, const Standard_Boolean VPer, gp_Pnt& P, gp_Vec& Vu, gp_Vec& Vv, gp_Vec& Vuu, gp_Vec& Vvv, gp_Vec& Vuv, gp_Vec& Vuuu, gp_Vec& Vvvv, gp_Vec& Vuuv, gp_Vec& Vuvv) { Standard_Boolean rational; // Standard_Integer k,dim,dim2; Standard_Integer dim,dim2; Standard_Real u1,u2; Standard_Integer d1,d2; Standard_Real *result, *resVu, *resVv, *resVuu, *resVvv, *resVuv, *resVuuu, *resVvvv, *resVuuv, *resVuvv; BSplSLib_DataContainer dc (UDegree, VDegree); if (PrepareEval (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer, Poles,Weights,UKnots,VKnots,UMults,VMults, u1,u2,d1,d2,rational,dc)) { if (rational) { dim = 4; dim2 = (d2 + 1) << 2; BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles); BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles); BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2)); if (d1 > 1) BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1))); if (d1 > 2) BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2)); BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders); result = dc.ders; resVu = result + 12; resVv = result + 3; resVuu = result + 24; resVvv = result + 6; resVuv = result + 15; resVuuu = result + 36; resVvvv = result + 9; resVuuv = result + 27; resVuvv = result + 18; } else { dim = 3; dim2 = (d2 + 1); dim2 = (dim2 << 1) + dim2; BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles); BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles); BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2)); if (d1 > 1) BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1))); if (d1 > 2) BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2)); result = dc.poles; resVu = result + dim2; resVv = result + 3; if (UDegree <= 1) { resVuu = BSplSLib_zero; resVuuv = BSplSLib_zero; } else { resVuu = result + (dim2 << 1); resVuuv = result + (dim2 << 1) + 3; } if (VDegree <= 1) { resVvv = BSplSLib_zero; resVuvv = BSplSLib_zero; } else { resVvv = result + 6; resVuvv = result + dim2 + 6; } resVuv = result + (d2 << 1) + d2 + 6; if (UDegree <= 2) resVuuu = BSplSLib_zero; else resVuuu = result + (dim2 << 1) + dim2; if (VDegree <= 2) resVvvv = BSplSLib_zero; else resVvvv = result + 9; } } else { if (rational) { dim = 4; dim2 = (d2 + 1) << 2; BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles); BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles); BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2)); if (d1 > 1) BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1))); if (d1 > 2) BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2)); BSplSLib::RationalDerivative(d1,d2,3,3,*dc.poles,*dc.ders); result = dc.ders; resVu = result + 3; resVv = result + 12; resVuu = result + 6; resVvv = result + 24; resVuv = result + 15; resVuuu = result + 9; resVvvv = result + 36; resVuuv = result + 18; resVuvv = result + 27; } else { dim = 3; dim2 = (d2 + 1); dim2 = (dim2 << 1) + dim2; BSplCLib::Bohm (u1,d1,3,*dc.knots1,dim2,*dc.poles); BSplCLib::Bohm (u2,d2,3,*dc.knots2,dim ,*dc.poles); BSplCLib::Bohm (u2,d2,2,*dc.knots2,dim ,*(dc.poles + dim2)); if (d1 > 1) BSplCLib::Bohm(u2,d2,1,*dc.knots2,dim ,*(dc.poles + (dim2 << 1))); if (d1 > 2) BSplCLib::Eval(u2,d2 ,*dc.knots2,dim ,*(dc.poles + (dim2 << 1) + dim2)); result = dc.poles; resVu = result + 3; resVv = result + dim2; if (UDegree <= 1) { resVuu = BSplSLib_zero; resVuuv = BSplSLib_zero; } else { resVuu = result + 6; resVuuv = result + dim2 + 6; } if (VDegree <= 1) { resVvv = BSplSLib_zero; resVuvv = BSplSLib_zero; } else { resVvv = result + (dim2 << 1); resVuvv = result + (dim2 << 1) + 3; } resVuv = result + (d2 << 1) + d2 + 6; if (UDegree <= 2) resVuuu = BSplSLib_zero; else resVuuu = result + 9; if (VDegree <= 2) resVvvv = BSplSLib_zero; else resVvvv = result + (dim2 << 1) + dim2; } } P .SetX(result [0]); Vu .SetX(resVu [0]); Vv .SetX(resVv [0]); Vuu .SetX(resVuu [0]); Vvv .SetX(resVvv [0]); Vuv .SetX(resVuv [0]); Vuuu.SetX(resVuuu[0]); Vvvv.SetX(resVvvv[0]); Vuuv.SetX(resVuuv[0]); Vuvv.SetX(resVuvv[0]); P .SetY(result [1]); Vu .SetY(resVu [1]); Vv .SetY(resVv [1]); Vuu .SetY(resVuu [1]); Vvv .SetY(resVvv [1]); Vuv .SetY(resVuv [1]); Vuuu.SetY(resVuuu[1]); Vvvv.SetY(resVvvv[1]); Vuuv.SetY(resVuuv[1]); Vuvv.SetY(resVuvv[1]); P .SetZ(result [2]); Vu .SetZ(resVu [2]); Vv .SetZ(resVv [2]); Vuu .SetZ(resVuu [2]); Vvv .SetZ(resVvv [2]); Vuv .SetZ(resVuv [2]); Vuuu.SetZ(resVuuu[2]); Vvvv.SetZ(resVvvv[2]); Vuuv.SetZ(resVuuv[2]); Vuvv.SetZ(resVuvv[2]); } //======================================================================= //function : DN //purpose : //======================================================================= void BSplSLib::DN (const Standard_Real U, const Standard_Real V, const Standard_Integer Nu, const Standard_Integer Nv, const Standard_Integer UIndex, const Standard_Integer VIndex, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, const TColStd_Array1OfReal& UKnots, const TColStd_Array1OfReal& VKnots, const TColStd_Array1OfInteger& UMults, const TColStd_Array1OfInteger& VMults, const Standard_Integer UDegree, const Standard_Integer VDegree, const Standard_Boolean URat, const Standard_Boolean VRat, const Standard_Boolean UPer, const Standard_Boolean VPer, gp_Vec& Vn) { Standard_Boolean rational; Standard_Integer k,dim; Standard_Real u1,u2; Standard_Integer d1,d2; BSplSLib_DataContainer dc (UDegree, VDegree); Standard_Boolean ufirst = PrepareEval (U,V,UIndex,VIndex,UDegree,VDegree,URat,VRat,UPer,VPer, Poles,Weights,UKnots,VKnots,UMults,VMults, u1,u2,d1,d2,rational,dc); dim = rational ? 4 : 3; if (!rational) { if ((Nu > UDegree) || (Nv > VDegree)) { Vn.SetX(0.); Vn.SetY(0.); Vn.SetZ(0.); return; } } Standard_Integer n1 = ufirst ? Nu : Nv; Standard_Integer n2 = ufirst ? Nv : Nu; BSplCLib::Bohm(u1,d1,n1,*dc.knots1,dim * (d2 + 1),*dc.poles); for (k = 0; k <= Min(n1,d1); k++) BSplCLib::Bohm(u2,d2,n2,*dc.knots2,dim,*(dc.poles+k*dim*(d2+1))); Standard_Real *result; if (rational) { BSplSLib::RationalDerivative(d1,d2,n1,n2,*dc.poles,*dc.ders,Standard_False); result = dc.ders; // because Standard_False ci-dessus. } else { result = dc.poles + (n1 * (d2+1) + n2) * dim ; } Vn.SetX(result[0]); Vn.SetY(result[1]); Vn.SetZ(result[2]); } // // Surface modifications // // a surface is processed as a curve of curves. // i.e : a curve of parameter u where the current point is the set of poles // of the iso. // //======================================================================= //function : Iso //purpose : //======================================================================= void BSplSLib::Iso(const Standard_Real Param, const Standard_Boolean IsU, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger& Mults, const Standard_Integer Degree, const Standard_Boolean Periodic, TColgp_Array1OfPnt& CPoles, TColStd_Array1OfReal& CWeights) { Standard_Integer index = 0; Standard_Real u = Param; Standard_Boolean rational = &Weights != NULL; Standard_Integer dim = rational ? 4 : 3; // compute local knots NCollection_LocalArray locknots1 (2*Degree); BSplCLib::LocateParameter(Degree,Knots,Mults,u,Periodic,index,u); BSplCLib::BuildKnots(Degree,index,Periodic,Knots,Mults,*locknots1); if (&Mults == NULL) index -= Knots.Lower() + Degree; else index = BSplCLib::PoleIndex(Degree,index,Periodic,Mults); // copy the local poles // Standard_Integer f1,l1,f2,l2,i,j,k; Standard_Integer f1,l1,f2,l2,i,j; if (IsU) { f1 = Poles.LowerRow(); l1 = Poles.UpperRow(); f2 = Poles.LowerCol(); l2 = Poles.UpperCol(); } else { f1 = Poles.LowerCol(); l1 = Poles.UpperCol(); f2 = Poles.LowerRow(); l2 = Poles.UpperRow(); } NCollection_LocalArray locpoles ((Degree+1) * (l2-f2+1) * dim); Standard_Real w, *pole = locpoles; index += f1; for (i = 0; i <= Degree; i++) { for (j = f2; j <= l2; j++) { const gp_Pnt& P = IsU ? Poles(index,j) : Poles(j,index); if (rational) { pole[3] = w = IsU ? Weights(index,j) : Weights(j,index); pole[0] = P.X() * w; pole[1] = P.Y() * w; pole[2] = P.Z() * w; } else { pole[0] = P.X(); pole[1] = P.Y(); pole[2] = P.Z(); } pole += dim; } index++; if (index > l1) index = f1; } // compute the iso BSplCLib::Eval(u,Degree,*locknots1,(l2-f2+1)*dim,*locpoles); // get the result pole = locpoles; for (i = CPoles.Lower(); i <= CPoles.Upper(); i++) { gp_Pnt& P = CPoles(i); if (rational) { CWeights(i) = w = pole[3]; P.SetX( pole[0] / w); P.SetY( pole[1] / w); P.SetZ( pole[2] / w); } else { P.SetX( pole[0]); P.SetY( pole[1]); P.SetZ( pole[2]); } pole += dim; } // if the input is not rational but weights are wanted if (!rational && (&CWeights != NULL)) { for (i = CWeights.Lower(); i <= CWeights.Upper(); i++) CWeights(i) = 1.; } } //======================================================================= //function : Reverse //purpose : //======================================================================= void BSplSLib::Reverse( TColgp_Array2OfPnt& Poles, const Standard_Integer Last, const Standard_Boolean UDirection) { Standard_Integer i,j, l = Last; if ( UDirection) { l = Poles.LowerRow() + (l - Poles.LowerRow())%(Poles.ColLength()); TColgp_Array2OfPnt temp(0, Poles.ColLength()-1, Poles.LowerCol(), Poles.UpperCol()); for (i = Poles.LowerRow(); i <= l; i++) { for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) { temp(l-i,j) = Poles(i,j); } } for (i = l+1; i <= Poles.UpperRow(); i++) { for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) { temp(l+Poles.ColLength()-i,j) = Poles(i,j); } } for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) { for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) { Poles(i,j) = temp (i-Poles.LowerRow(),j); } } } else { l = Poles.LowerCol() + (l - Poles.LowerCol())%(Poles.RowLength()); TColgp_Array2OfPnt temp(Poles.LowerRow(), Poles.UpperRow(), 0, Poles.RowLength()-1); for (j = Poles.LowerCol(); j <= l; j++) { for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) { temp(i,l-j) = Poles(i,j); } } for (j = l+1; j <= Poles.UpperCol(); j++) { for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) { temp(i,l+Poles.RowLength()-j) = Poles(i,j); } } for (i = Poles.LowerRow(); i <= Poles.UpperRow(); i++) { for (j = Poles.LowerCol(); j <= Poles.UpperCol(); j++) { Poles(i,j) = temp (i,j-Poles.LowerCol()); } } } } //======================================================================= //function : Reverse //purpose : //======================================================================= void BSplSLib::Reverse( TColStd_Array2OfReal& Weights, const Standard_Integer Last, const Standard_Boolean UDirection) { Standard_Integer i,j, l = Last; if ( UDirection) { l = Weights.LowerRow() + (l - Weights.LowerRow())%(Weights.ColLength()); TColStd_Array2OfReal temp(0, Weights.ColLength()-1, Weights.LowerCol(), Weights.UpperCol()); for (i = Weights.LowerRow(); i <= l; i++) { for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) { temp(l-i,j) = Weights(i,j); } } for (i = l+1; i <= Weights.UpperRow(); i++) { for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) { temp(l+Weights.ColLength()-i,j) = Weights(i,j); } } for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) { for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) { Weights(i,j) = temp (i-Weights.LowerRow(),j); } } } else { l = Weights.LowerCol() + (l - Weights.LowerCol())%(Weights.RowLength()); TColStd_Array2OfReal temp(Weights.LowerRow(), Weights.UpperRow(), 0, Weights.RowLength()-1); for (j = Weights.LowerCol(); j <= l; j++) { for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) { temp(i,l-j) = Weights(i,j); } } for (j = l+1; j <= Weights.UpperCol(); j++) { for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) { temp(i,l+Weights.RowLength()-j) = Weights(i,j); } } for (i = Weights.LowerRow(); i <= Weights.UpperRow(); i++) { for (j = Weights.LowerCol(); j <= Weights.UpperCol(); j++) { Weights(i,j) = temp (i,j-Weights.LowerCol()); } } } } //======================================================================= //function : IsRational //purpose : //======================================================================= Standard_Boolean BSplSLib::IsRational (const TColStd_Array2OfReal& Weights, const Standard_Integer I1, const Standard_Integer I2, const Standard_Integer J1, const Standard_Integer J2, const Standard_Real Epsi) { Standard_Real eps = (Epsi > 0.0) ? Epsi : Epsilon(Weights(I1,I2)); Standard_Integer i,j; Standard_Integer fi = Weights.LowerRow(), li = Weights.ColLength(); Standard_Integer fj = Weights.LowerCol(), lj = Weights.RowLength(); for (i = I1 - fi; i < I2 - fi; i++) { for (j = J1 - fj; j < J2 - fj; j++) { if (Abs(Weights(fi+i%li,fj+j%lj)-Weights(fi+(i+1)%li,fj+j%lj))>eps) return Standard_True; } } return Standard_False; } //======================================================================= //function : SetPoles //purpose : //======================================================================= void BSplSLib::SetPoles(const TColgp_Array2OfPnt& Poles, TColStd_Array1OfReal& FP, const Standard_Boolean UDirection) { Standard_Integer i, j, l = FP.Lower(); Standard_Integer PLowerRow = Poles.LowerRow(); Standard_Integer PUpperRow = Poles.UpperRow(); Standard_Integer PLowerCol = Poles.LowerCol(); Standard_Integer PUpperCol = Poles.UpperCol(); if (UDirection) { for ( i = PLowerRow; i <= PUpperRow; i++) { for ( j = PLowerCol; j <= PUpperCol; j++) { const gp_Pnt& P = Poles.Value(i,j); FP(l) = P.X(); l++; FP(l) = P.Y(); l++; FP(l) = P.Z(); l++; } } } else { for ( j = PLowerCol; j <= PUpperCol; j++) { for ( i = PLowerRow; i <= PUpperRow; i++) { const gp_Pnt& P = Poles.Value(i,j); FP(l) = P.X(); l++; FP(l) = P.Y(); l++; FP(l) = P.Z(); l++; } } } } //======================================================================= //function : SetPoles //purpose : //======================================================================= void BSplSLib::SetPoles(const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, TColStd_Array1OfReal& FP, const Standard_Boolean UDirection) { Standard_Integer i, j, l = FP.Lower(); Standard_Integer PLowerRow = Poles.LowerRow(); Standard_Integer PUpperRow = Poles.UpperRow(); Standard_Integer PLowerCol = Poles.LowerCol(); Standard_Integer PUpperCol = Poles.UpperCol(); if (UDirection) { for ( i = PLowerRow; i <= PUpperRow; i++) { for ( j = PLowerCol; j <= PUpperCol; j++) { const gp_Pnt& P = Poles .Value(i,j); Standard_Real w = Weights.Value(i,j); FP(l) = P.X() * w; l++; FP(l) = P.Y() * w; l++; FP(l) = P.Z() * w; l++; FP(l) = w; l++; } } } else { for ( j = PLowerCol; j <= PUpperCol; j++) { for ( i = PLowerRow; i <= PUpperRow; i++) { const gp_Pnt& P = Poles .Value(i,j); Standard_Real w = Weights.Value(i,j); FP(l) = P.X() * w; l++; FP(l) = P.Y() * w; l++; FP(l) = P.Z() * w; l++; FP(l) = w; l++; } } } } //======================================================================= //function : GetPoles //purpose : //======================================================================= void BSplSLib::GetPoles(const TColStd_Array1OfReal& FP, TColgp_Array2OfPnt& Poles, const Standard_Boolean UDirection) { Standard_Integer i, j, l = FP.Lower(); Standard_Integer PLowerRow = Poles.LowerRow(); Standard_Integer PUpperRow = Poles.UpperRow(); Standard_Integer PLowerCol = Poles.LowerCol(); Standard_Integer PUpperCol = Poles.UpperCol(); if (UDirection) { for ( i = PLowerRow; i <= PUpperRow; i++) { for ( j = PLowerCol; j <= PUpperCol; j++) { gp_Pnt& P = Poles.ChangeValue(i,j); P.SetX(FP(l)); l++; P.SetY(FP(l)); l++; P.SetZ(FP(l)); l++; } } } else { for ( j = PLowerCol; j <= PUpperCol; j++) { for ( i = PLowerRow; i <= PUpperRow; i++) { gp_Pnt& P = Poles.ChangeValue(i,j); P.SetX(FP(l)); l++; P.SetY(FP(l)); l++; P.SetZ(FP(l)); l++; } } } } //======================================================================= //function : GetPoles //purpose : //======================================================================= void BSplSLib::GetPoles(const TColStd_Array1OfReal& FP, TColgp_Array2OfPnt& Poles, TColStd_Array2OfReal& Weights, const Standard_Boolean UDirection) { Standard_Integer i, j, l = FP.Lower(); Standard_Integer PLowerRow = Poles.LowerRow(); Standard_Integer PUpperRow = Poles.UpperRow(); Standard_Integer PLowerCol = Poles.LowerCol(); Standard_Integer PUpperCol = Poles.UpperCol(); if (UDirection) { for ( i = PLowerRow; i <= PUpperRow; i++) { for ( j = PLowerCol; j <= PUpperCol; j++) { Standard_Real w = FP( l + 3); Weights(i,j) = w; gp_Pnt& P = Poles.ChangeValue(i,j); P.SetX(FP(l) / w); l++; P.SetY(FP(l) / w); l++; P.SetZ(FP(l) / w); l++; l++; } } } else { for ( j = PLowerCol; j <= PUpperCol; j++) { for ( i = PLowerRow; i <= PUpperRow; i++) { Standard_Real w = FP( l + 3); Weights(i,j) = w; gp_Pnt& P = Poles.ChangeValue(i,j); P.SetX(FP(l) / w); l++; P.SetY(FP(l) / w); l++; P.SetZ(FP(l) / w); l++; l++; } } } } //======================================================================= //function : InsertKnots //purpose : //======================================================================= void BSplSLib::InsertKnots(const Standard_Boolean UDirection, const Standard_Integer Degree, const Standard_Boolean Periodic, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger& Mults, const TColStd_Array1OfReal& AddKnots, const TColStd_Array1OfInteger& AddMults, TColgp_Array2OfPnt& NewPoles, TColStd_Array2OfReal& NewWeights, TColStd_Array1OfReal& NewKnots, TColStd_Array1OfInteger& NewMults, const Standard_Real Epsilon, const Standard_Boolean Add ) { Standard_Boolean rational = &Weights != NULL; Standard_Integer dim = 3; if (rational) dim++; TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength()); TColStd_Array1OfReal newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength()); if (rational) SetPoles(Poles,Weights,poles,UDirection); else SetPoles(Poles,poles,UDirection); if (UDirection) { dim *= Poles.RowLength(); } else { dim *= Poles.ColLength(); } BSplCLib::InsertKnots(Degree,Periodic,dim,poles,Knots,Mults, AddKnots,AddMults,newpoles,NewKnots,NewMults, Epsilon,Add); if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection); else GetPoles(newpoles,NewPoles,UDirection); } //======================================================================= //function : RemoveKnot //purpose : //======================================================================= Standard_Boolean BSplSLib::RemoveKnot (const Standard_Boolean UDirection, const Standard_Integer Index, const Standard_Integer Mult, const Standard_Integer Degree, const Standard_Boolean Periodic, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger& Mults, TColgp_Array2OfPnt& NewPoles, TColStd_Array2OfReal& NewWeights, TColStd_Array1OfReal& NewKnots, TColStd_Array1OfInteger& NewMults, const Standard_Real Tolerance) { Standard_Boolean rational = &Weights != NULL; Standard_Integer dim = 3; if (rational) dim++; TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength()); TColStd_Array1OfReal newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength()); if (rational) SetPoles(Poles,Weights,poles,UDirection); else SetPoles(Poles,poles,UDirection); if (UDirection) { dim *= Poles.RowLength(); } else { dim *= Poles.ColLength(); } if ( !BSplCLib::RemoveKnot(Index,Mult,Degree,Periodic,dim, poles,Knots,Mults,newpoles,NewKnots,NewMults, Tolerance)) return Standard_False; if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection); else GetPoles(newpoles,NewPoles,UDirection); return Standard_True; } //======================================================================= //function : IncreaseDegree //purpose : //======================================================================= void BSplSLib::IncreaseDegree (const Standard_Boolean UDirection, const Standard_Integer Degree, const Standard_Integer NewDegree, const Standard_Boolean Periodic, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, const TColStd_Array1OfReal& Knots, const TColStd_Array1OfInteger& Mults, TColgp_Array2OfPnt& NewPoles, TColStd_Array2OfReal& NewWeights, TColStd_Array1OfReal& NewKnots, TColStd_Array1OfInteger& NewMults) { Standard_Boolean rational = &Weights != NULL; Standard_Integer dim = 3; if (rational) dim++; TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength()); TColStd_Array1OfReal newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength()); if (rational) SetPoles(Poles,Weights,poles,UDirection); else SetPoles(Poles,poles,UDirection); if (UDirection) { dim *= Poles.RowLength(); } else { dim *= Poles.ColLength(); } BSplCLib::IncreaseDegree(Degree,NewDegree,Periodic,dim,poles,Knots,Mults, newpoles,NewKnots,NewMults); if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection); else GetPoles(newpoles,NewPoles,UDirection); } //======================================================================= //function : Unperiodize //purpose : //======================================================================= void BSplSLib::Unperiodize (const Standard_Boolean UDirection, const Standard_Integer Degree, const TColStd_Array1OfInteger& Mults, const TColStd_Array1OfReal& Knots, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, TColStd_Array1OfInteger& NewMults, TColStd_Array1OfReal& NewKnots, TColgp_Array2OfPnt& NewPoles, TColStd_Array2OfReal& NewWeights) { Standard_Boolean rational = &Weights != NULL; Standard_Integer dim = 3; if (rational) dim++; TColStd_Array1OfReal poles( 1, dim*Poles.RowLength()*Poles.ColLength()); TColStd_Array1OfReal newpoles( 1, dim*NewPoles.RowLength()*NewPoles.ColLength()); if (rational) SetPoles(Poles,Weights,poles,UDirection); else SetPoles(Poles,poles,UDirection); if (UDirection) { dim *= Poles.RowLength(); } else { dim *= Poles.ColLength(); } BSplCLib::Unperiodize(Degree, dim, Mults, Knots, poles, NewMults, NewKnots, newpoles); if (rational) GetPoles(newpoles,NewPoles,NewWeights,UDirection); else GetPoles(newpoles,NewPoles,UDirection); } //======================================================================= //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 BSplSLib::BuildCache (const Standard_Real U, const Standard_Real V, const Standard_Real USpanDomain, const Standard_Real VSpanDomain, const Standard_Boolean UPeriodic, const Standard_Boolean VPeriodic, const Standard_Integer UDegree, const Standard_Integer VDegree, const Standard_Integer UIndex, const Standard_Integer VIndex, const TColStd_Array1OfReal& UFlatKnots, const TColStd_Array1OfReal& VFlatKnots, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, TColgp_Array2OfPnt& CachePoles, TColStd_Array2OfReal& CacheWeights) { Standard_Boolean rational,rational_u,rational_v,flag_u_or_v; Standard_Integer kk,d1,d1p1,d2,d2p1,ii,jj,iii,jjj,Index; Standard_Real u1,min_degree_domain,max_degree_domain,f,factor[2],u2; if (&Weights != NULL) rational_u = rational_v = Standard_True; else rational_u = rational_v = Standard_False; BSplSLib_DataContainer dc (UDegree, VDegree); flag_u_or_v = PrepareEval (U, V, UIndex, VIndex, UDegree, VDegree, rational_u, rational_v, UPeriodic, VPeriodic, Poles, Weights, UFlatKnots, VFlatKnots, (BSplCLib::NoMults()), (BSplCLib::NoMults()), u1, u2, d1, d2, rational, dc); d1p1 = d1 + 1; d2p1 = d2 + 1; if (rational) { BSplCLib::Bohm(u1,d1,d1,*dc.knots1,4 * d2p1,*dc.poles); for (kk = 0; kk <= d1 ; kk++) BSplCLib::Bohm(u2,d2,d2,*dc.knots2,4,*(dc.poles + kk * 4 * d2p1)); if (flag_u_or_v) { min_degree_domain = USpanDomain ; max_degree_domain = VSpanDomain ; } else { min_degree_domain = VSpanDomain ; max_degree_domain = USpanDomain ; } factor[0] = 1.0e0 ; for (ii = 0 ; ii <= d2 ; ii++) { iii = ii + 1; factor[1] = 1.0e0 ; for (jj = 0 ; jj <= d1 ; jj++) { jjj = jj + 1; Index = jj * d2p1 + ii ; Index = Index << 2; gp_Pnt& P = CachePoles(iii,jjj); f = factor[0] * factor[1]; P.SetX( f * dc.poles[Index]); Index++; P.SetY( f * dc.poles[Index]); Index++; P.SetZ( f * dc.poles[Index]); Index++; CacheWeights(iii ,jjj) = f * dc.poles[Index] ; factor[1] *= min_degree_domain / (Standard_Real) (jjj) ; } factor[0] *= max_degree_domain / (Standard_Real) (iii) ; } } else { BSplCLib::Bohm(u1,d1,d1,*dc.knots1,3 * d2p1,*dc.poles); for (kk = 0; kk <= d1 ; kk++) BSplCLib::Bohm(u2,d2,d2,*dc.knots2,3,*(dc.poles + kk * 3 * d2p1)); if (flag_u_or_v) { min_degree_domain = USpanDomain ; max_degree_domain = VSpanDomain ; } else { min_degree_domain = VSpanDomain ; max_degree_domain = USpanDomain ; } factor[0] = 1.0e0 ; for (ii = 0 ; ii <= d2 ; ii++) { iii = ii + 1; factor[1] = 1.0e0 ; for (jj = 0 ; jj <= d1 ; jj++) { jjj = jj + 1; Index = jj * d2p1 + ii ; Index = (Index << 1) + Index; gp_Pnt& P = CachePoles(iii,jjj); f = factor[0] * factor[1]; P.SetX( f * dc.poles[Index]); Index++; P.SetY( f * dc.poles[Index]); Index++; P.SetZ( f * dc.poles[Index]); factor[1] *= min_degree_domain / (Standard_Real) (jjj) ; } factor[0] *= max_degree_domain / (Standard_Real) (iii) ; } if (&Weights != NULL) { // // means that PrepareEval did found out that the surface was // locally polynomial but since the surface is constructed // with some weights we need to set the weight polynomial to constant // for (ii = 1 ; ii <= d2p1 ; ii++) { for (jj = 1 ; jj <= d1p1 ; jj++) { CacheWeights(ii,jj) = 0.0e0 ; } } CacheWeights(1,1) = 1.0e0 ; } } } //======================================================================= //function : CacheD0 //purpose : Evaluates the polynomial cache of the Bspline Curve // //======================================================================= void BSplSLib::CacheD0(const Standard_Real UParameter, const Standard_Real VParameter, const Standard_Integer UDegree, const Standard_Integer VDegree, const Standard_Real UCacheParameter, const Standard_Real VCacheParameter, const Standard_Real USpanLenght, const Standard_Real VSpanLenght, const TColgp_Array2OfPnt& PolesArray, const TColStd_Array2OfReal& WeightsArray, gp_Pnt& 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_Integer // ii, dimension, min_degree, max_degree ; Standard_Real new_parameter[2] , inverse ; Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ; Standard_Real * myPoint = (Standard_Real *) &aPoint ; if (UDegree <= VDegree) { min_degree = UDegree ; max_degree = VDegree ; new_parameter[1] = (UParameter - UCacheParameter) / USpanLenght ; new_parameter[0] = (VParameter - VCacheParameter) / VSpanLenght ; dimension = 3 * (UDegree + 1) ; } else { min_degree = VDegree ; max_degree = UDegree ; new_parameter[0] = (UParameter - UCacheParameter) / USpanLenght ; new_parameter[1] = (VParameter - VCacheParameter) / VSpanLenght ; dimension = 3 * (VDegree + 1) ; } NCollection_LocalArray locpoles(dimension); PLib::NoDerivativeEvalPolynomial(new_parameter[0], max_degree, dimension, max_degree*dimension, PArray[0], locpoles[0]) ; PLib::NoDerivativeEvalPolynomial(new_parameter[1], min_degree, 3, (min_degree << 1) + min_degree, locpoles[0], myPoint[0]) ; if (&WeightsArray != NULL) { dimension = min_degree + 1 ; Standard_Real * WArray = (Standard_Real *) &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ; PLib::NoDerivativeEvalPolynomial(new_parameter[0], max_degree, dimension, max_degree*dimension, WArray[0], locpoles[0]) ; PLib::NoDerivativeEvalPolynomial(new_parameter[1], min_degree, 1, min_degree, locpoles[0], inverse) ; inverse = 1.0e0/ inverse ; myPoint[0] *= inverse ; myPoint[1] *= inverse ; myPoint[2] *= inverse ; } } //======================================================================= //function : CacheD1 //purpose : Evaluates the polynomial cache of the Bspline Curve // //======================================================================= void BSplSLib::CacheD1(const Standard_Real UParameter, const Standard_Real VParameter, const Standard_Integer UDegree, const Standard_Integer VDegree, const Standard_Real UCacheParameter, const Standard_Real VCacheParameter, const Standard_Real USpanLenght, const Standard_Real VSpanLenght, const TColgp_Array2OfPnt& PolesArray, const TColStd_Array2OfReal& WeightsArray, gp_Pnt& aPoint, gp_Vec& aVecU, gp_Vec& aVecV) { // // 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, // jj, // kk, dimension, min_degree, max_degree ; Standard_Real inverse_min, inverse_max, new_parameter[2] ; Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ; Standard_Real local_poles_array[2][2][3], local_poles_and_weights_array[2][2][4], local_weights_array[2][2] ; Standard_Real * my_vec_min, * my_vec_max, * my_point ; my_point = (Standard_Real *) &aPoint ; // // initialize in case of rational evaluation // because RationalDerivative will use all // the coefficients // // if (&WeightsArray != NULL) { local_poles_array [0][0][0] = 0.0e0 ; local_poles_array [0][0][1] = 0.0e0 ; local_poles_array [0][0][2] = 0.0e0 ; local_weights_array [0][0] = 0.0e0 ; local_poles_and_weights_array[0][0][0] = 0.0e0 ; local_poles_and_weights_array[0][0][1] = 0.0e0 ; local_poles_and_weights_array[0][0][2] = 0.0e0 ; local_poles_and_weights_array[0][0][3] = 0.0e0 ; local_poles_array [0][1][0] = 0.0e0 ; local_poles_array [0][1][1] = 0.0e0 ; local_poles_array [0][1][2] = 0.0e0 ; local_weights_array [0][1] = 0.0e0 ; local_poles_and_weights_array[0][1][0] = 0.0e0 ; local_poles_and_weights_array[0][1][1] = 0.0e0 ; local_poles_and_weights_array[0][1][2] = 0.0e0 ; local_poles_and_weights_array[0][1][3] = 0.0e0 ; local_poles_array [1][0][0] = 0.0e0 ; local_poles_array [1][0][1] = 0.0e0 ; local_poles_array [1][0][2] = 0.0e0 ; local_weights_array [1][0] = 0.0e0 ; local_poles_and_weights_array[1][0][0] = 0.0e0 ; local_poles_and_weights_array[1][0][1] = 0.0e0 ; local_poles_and_weights_array[1][0][2] = 0.0e0 ; local_poles_and_weights_array[1][0][3] = 0.0e0 ; local_poles_array [1][1][0] = 0.0e0 ; local_poles_array [1][1][1] = 0.0e0 ; local_poles_array [1][1][2] = 0.0e0 ; local_weights_array [1][1] = 0.0e0 ; local_poles_and_weights_array[1][1][0] = 0.0e0 ; local_poles_and_weights_array[1][1][1] = 0.0e0 ; local_poles_and_weights_array[1][1][2] = 0.0e0 ; local_poles_and_weights_array[1][1][3] = 0.0e0 ; } if (UDegree <= VDegree) { min_degree = UDegree ; max_degree = VDegree ; inverse_min = 1.0e0/USpanLenght ; inverse_max = 1.0e0/VSpanLenght ; new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ; new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ; dimension = 3 * (UDegree + 1) ; my_vec_min = (Standard_Real *) &aVecU ; my_vec_max = (Standard_Real *) &aVecV ; } else { min_degree = VDegree ; max_degree = UDegree ; inverse_min = 1.0e0/VSpanLenght ; inverse_max = 1.0e0/USpanLenght ; new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ; new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ; dimension = 3 * (VDegree + 1) ; my_vec_min = (Standard_Real *) &aVecV ; my_vec_max = (Standard_Real *) &aVecU ; } NCollection_LocalArray locpoles (2 * dimension); PLib::EvalPolynomial(new_parameter[0], 1, max_degree, dimension, PArray[0], locpoles[0]) ; PLib::EvalPolynomial(new_parameter[1], 1, min_degree, 3, locpoles[0], local_poles_array[0][0][0]) ; PLib::NoDerivativeEvalPolynomial(new_parameter[1], min_degree, 3, (min_degree << 1) + min_degree, locpoles[dimension], local_poles_array[1][0][0]) ; if (&WeightsArray != NULL) { dimension = min_degree + 1 ; Standard_Real * WArray = (Standard_Real *) &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ; PLib::EvalPolynomial(new_parameter[0], 1, max_degree, dimension, WArray[0], locpoles[0]) ; PLib::EvalPolynomial(new_parameter[1], 1, min_degree, 1, locpoles[0], local_weights_array[0][0]) ; PLib::NoDerivativeEvalPolynomial(new_parameter[1], min_degree, 1, min_degree, locpoles[dimension], local_weights_array[1][0]) ; local_poles_and_weights_array[0][0][0] = local_poles_array [0][0][0] ; local_poles_and_weights_array[0][0][1] = local_poles_array [0][0][1] ; local_poles_and_weights_array[0][0][2] = local_poles_array [0][0][2] ; local_poles_and_weights_array[0][0][3] = local_weights_array[0][0] ; local_poles_and_weights_array[0][1][0] = local_poles_array [0][1][0] ; local_poles_and_weights_array[0][1][1] = local_poles_array [0][1][1] ; local_poles_and_weights_array[0][1][2] = local_poles_array [0][1][2] ; local_poles_and_weights_array[0][1][3] = local_weights_array[0][1] ; local_poles_and_weights_array[1][0][0] = local_poles_array [1][0][0] ; local_poles_and_weights_array[1][0][1] = local_poles_array [1][0][1] ; local_poles_and_weights_array[1][0][2] = local_poles_array [1][0][2] ; local_poles_and_weights_array[1][0][3] = local_weights_array[1][0] ; local_poles_and_weights_array[1][1][0] = local_poles_array [1][1][0] ; local_poles_and_weights_array[1][1][1] = local_poles_array [1][1][1] ; local_poles_and_weights_array[1][1][2] = local_poles_array [1][1][2] ; local_poles_and_weights_array[1][1][3] = local_weights_array[1][1] ; BSplSLib::RationalDerivative(1, 1, 1, 1, local_poles_and_weights_array[0][0][0], local_poles_array[0][0][0]) ; } my_point [0] = local_poles_array [0][0][0] ; my_vec_min[0] = inverse_min * local_poles_array[0][1][0] ; my_vec_max[0] = inverse_max * local_poles_array[1][0][0] ; my_point [1] = local_poles_array [0][0][1] ; my_vec_min[1] = inverse_min * local_poles_array[0][1][1] ; my_vec_max[1] = inverse_max * local_poles_array[1][0][1] ; my_point [2] = local_poles_array [0][0][2] ; my_vec_min[2] = inverse_min * local_poles_array[0][1][2] ; my_vec_max[2] = inverse_max * local_poles_array[1][0][2] ; } //======================================================================= //function : CacheD2 //purpose : Evaluates the polynomial cache of the Bspline Curve // //======================================================================= void BSplSLib::CacheD2(const Standard_Real UParameter, const Standard_Real VParameter, const Standard_Integer UDegree, const Standard_Integer VDegree, const Standard_Real UCacheParameter, const Standard_Real VCacheParameter, const Standard_Real USpanLenght, const Standard_Real VSpanLenght, const TColgp_Array2OfPnt& PolesArray, const TColStd_Array2OfReal& WeightsArray, gp_Pnt& aPoint, gp_Vec& aVecU, gp_Vec& aVecV, gp_Vec& aVecUU, gp_Vec& aVecUV, gp_Vec& aVecVV) { // // 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, // jj, kk, index, dimension, min_degree, max_degree ; Standard_Real inverse_min, inverse_max, new_parameter[2] ; Standard_Real * PArray = (Standard_Real *) &(PolesArray(PolesArray.LowerCol(), PolesArray.LowerRow())) ; Standard_Real local_poles_array[3][3][3], local_poles_and_weights_array[3][3][4], local_weights_array[3][3] ; Standard_Real * my_vec_min, * my_vec_max, * my_vec_min_min, * my_vec_max_max, * my_vec_min_max, * my_point ; my_point = (Standard_Real *) &aPoint ; // // initialize in case the min and max degree are less than 2 // local_poles_array[0][0][0] = 0.0e0 ; local_poles_array[0][0][1] = 0.0e0 ; local_poles_array[0][0][2] = 0.0e0 ; local_poles_array[0][1][0] = 0.0e0 ; local_poles_array[0][1][1] = 0.0e0 ; local_poles_array[0][1][2] = 0.0e0 ; local_poles_array[0][2][0] = 0.0e0 ; local_poles_array[0][2][1] = 0.0e0 ; local_poles_array[0][2][2] = 0.0e0 ; local_poles_array[1][0][0] = 0.0e0 ; local_poles_array[1][0][1] = 0.0e0 ; local_poles_array[1][0][2] = 0.0e0 ; local_poles_array[1][1][0] = 0.0e0 ; local_poles_array[1][1][1] = 0.0e0 ; local_poles_array[1][1][2] = 0.0e0 ; local_poles_array[1][2][0] = 0.0e0 ; local_poles_array[1][2][1] = 0.0e0 ; local_poles_array[1][2][2] = 0.0e0 ; local_poles_array[2][0][0] = 0.0e0 ; local_poles_array[2][0][1] = 0.0e0 ; local_poles_array[2][0][2] = 0.0e0 ; local_poles_array[2][1][0] = 0.0e0 ; local_poles_array[2][1][1] = 0.0e0 ; local_poles_array[2][1][2] = 0.0e0 ; local_poles_array[2][2][0] = 0.0e0 ; local_poles_array[2][2][1] = 0.0e0 ; local_poles_array[2][2][2] = 0.0e0 ; // // initialize in case of rational evaluation // because RationalDerivative will use all // the coefficients // // if (&WeightsArray != NULL) { local_poles_and_weights_array[0][0][0] = 0.0e0 ; local_poles_and_weights_array[0][0][1] = 0.0e0 ; local_poles_and_weights_array[0][0][2] = 0.0e0 ; local_poles_and_weights_array[0][1][0] = 0.0e0 ; local_poles_and_weights_array[0][1][1] = 0.0e0 ; local_poles_and_weights_array[0][1][2] = 0.0e0 ; local_poles_and_weights_array[0][2][0] = 0.0e0 ; local_poles_and_weights_array[0][2][1] = 0.0e0 ; local_poles_and_weights_array[0][2][2] = 0.0e0 ; local_poles_and_weights_array[1][0][0] = 0.0e0 ; local_poles_and_weights_array[1][0][1] = 0.0e0 ; local_poles_and_weights_array[1][0][2] = 0.0e0 ; local_poles_and_weights_array[1][1][0] = 0.0e0 ; local_poles_and_weights_array[1][1][1] = 0.0e0 ; local_poles_and_weights_array[1][1][2] = 0.0e0 ; local_poles_and_weights_array[1][2][0] = 0.0e0 ; local_poles_and_weights_array[1][2][1] = 0.0e0 ; local_poles_and_weights_array[1][2][2] = 0.0e0 ; local_poles_and_weights_array[2][0][0] = 0.0e0 ; local_poles_and_weights_array[2][0][1] = 0.0e0 ; local_poles_and_weights_array[2][0][2] = 0.0e0 ; local_poles_and_weights_array[2][1][0] = 0.0e0 ; local_poles_and_weights_array[2][1][1] = 0.0e0 ; local_poles_and_weights_array[2][1][2] = 0.0e0 ; local_poles_and_weights_array[2][2][0] = 0.0e0 ; local_poles_and_weights_array[2][2][1] = 0.0e0 ; local_poles_and_weights_array[2][2][2] = 0.0e0 ; local_poles_and_weights_array[0][0][3] = local_weights_array[0][0] = 0.0e0 ; local_poles_and_weights_array[0][1][3] = local_weights_array[0][1] = 0.0e0 ; local_poles_and_weights_array[0][2][3] = local_weights_array[0][2] = 0.0e0 ; local_poles_and_weights_array[1][0][3] = local_weights_array[1][0] = 0.0e0 ; local_poles_and_weights_array[1][1][3] = local_weights_array[1][1] = 0.0e0 ; local_poles_and_weights_array[1][2][3] = local_weights_array[1][2] = 0.0e0 ; local_poles_and_weights_array[2][0][3] = local_weights_array[2][0] = 0.0e0 ; local_poles_and_weights_array[2][1][3] = local_weights_array[2][1] = 0.0e0 ; local_poles_and_weights_array[2][2][3] = local_weights_array[2][2] = 0.0e0 ; } if (UDegree <= VDegree) { min_degree = UDegree ; max_degree = VDegree ; inverse_min = 1.0e0/USpanLenght ; inverse_max = 1.0e0/VSpanLenght ; new_parameter[0] = (VParameter - VCacheParameter) * inverse_max ; new_parameter[1] = (UParameter - UCacheParameter) * inverse_min ; dimension = 3 * (UDegree + 1) ; my_vec_min = (Standard_Real *) &aVecU ; my_vec_max = (Standard_Real *) &aVecV ; my_vec_min_min = (Standard_Real *) &aVecUU ; my_vec_min_max = (Standard_Real *) &aVecUV ; my_vec_max_max = (Standard_Real *) &aVecVV ; } else { min_degree = VDegree ; max_degree = UDegree ; inverse_min = 1.0e0/VSpanLenght ; inverse_max = 1.0e0/USpanLenght ; new_parameter[0] = (UParameter - UCacheParameter) * inverse_max ; new_parameter[1] = (VParameter - VCacheParameter) * inverse_min ; dimension = 3 * (VDegree + 1) ; my_vec_min = (Standard_Real *) &aVecV ; my_vec_max = (Standard_Real *) &aVecU ; my_vec_min_min = (Standard_Real *) &aVecVV ; my_vec_min_max = (Standard_Real *) &aVecUV ; my_vec_max_max = (Standard_Real *) &aVecUU ; } NCollection_LocalArray locpoles (3 * dimension); // // initialize in case min or max degree are less than 2 // Standard_Integer MinIndMax = 2; if ( max_degree < 2) MinIndMax = max_degree; Standard_Integer MinIndMin = 2; if ( min_degree < 2) MinIndMin = min_degree; index = MinIndMax * dimension ; for (ii = MinIndMax ; ii < 3 ; ii++) { for (kk = 0 ; kk < dimension ; kk++) { locpoles[index] = 0.0e0 ; index += 1 ; } } PLib::EvalPolynomial(new_parameter[0], MinIndMax, max_degree, dimension, PArray[0], locpoles[0]) ; PLib::EvalPolynomial(new_parameter[1], MinIndMin, min_degree, 3, locpoles[0], local_poles_array[0][0][0]) ; PLib::EvalPolynomial(new_parameter[1], 1, min_degree, 3, locpoles[dimension], local_poles_array[1][0][0]) ; PLib::NoDerivativeEvalPolynomial(new_parameter[1], min_degree, 3, (min_degree << 1) + min_degree, locpoles[dimension + dimension], local_poles_array[2][0][0]) ; if (&WeightsArray != NULL) { dimension = min_degree + 1 ; Standard_Real * WArray = (Standard_Real *) &WeightsArray(WeightsArray.LowerCol(),WeightsArray.LowerRow()) ; PLib::EvalPolynomial(new_parameter[0], MinIndMax, max_degree, dimension, WArray[0], locpoles[0]) ; PLib::EvalPolynomial(new_parameter[1], MinIndMin, min_degree, 1, locpoles[0], local_weights_array[0][0]) ; PLib::EvalPolynomial(new_parameter[1], 1, min_degree, 1, locpoles[dimension], local_weights_array[1][0]) ; PLib::NoDerivativeEvalPolynomial(new_parameter[1], min_degree, 1, min_degree, locpoles[dimension + dimension], local_weights_array[2][0]) ; local_poles_and_weights_array[0][0][0] = local_poles_array[0][0][0]; local_poles_and_weights_array[0][0][1] = local_poles_array[0][0][1]; local_poles_and_weights_array[0][0][2] = local_poles_array[0][0][2]; local_poles_and_weights_array[0][1][0] = local_poles_array[0][1][0]; local_poles_and_weights_array[0][1][1] = local_poles_array[0][1][1]; local_poles_and_weights_array[0][1][2] = local_poles_array[0][1][2]; local_poles_and_weights_array[0][2][0] = local_poles_array[0][2][0]; local_poles_and_weights_array[0][2][1] = local_poles_array[0][2][1]; local_poles_and_weights_array[0][2][2] = local_poles_array[0][2][2]; local_poles_and_weights_array[1][0][0] = local_poles_array[1][0][0]; local_poles_and_weights_array[1][0][1] = local_poles_array[1][0][1]; local_poles_and_weights_array[1][0][2] = local_poles_array[1][0][2]; local_poles_and_weights_array[1][1][0] = local_poles_array[1][1][0]; local_poles_and_weights_array[1][1][1] = local_poles_array[1][1][1]; local_poles_and_weights_array[1][1][2] = local_poles_array[1][1][2]; local_poles_and_weights_array[1][2][0] = local_poles_array[1][2][0]; local_poles_and_weights_array[1][2][1] = local_poles_array[1][2][1]; local_poles_and_weights_array[1][2][2] = local_poles_array[1][2][2]; local_poles_and_weights_array[2][0][0] = local_poles_array[2][0][0]; local_poles_and_weights_array[2][0][1] = local_poles_array[2][0][1]; local_poles_and_weights_array[2][0][2] = local_poles_array[2][0][2]; local_poles_and_weights_array[2][1][0] = local_poles_array[2][1][0]; local_poles_and_weights_array[2][1][1] = local_poles_array[2][1][1]; local_poles_and_weights_array[2][1][2] = local_poles_array[2][1][2]; local_poles_and_weights_array[2][2][0] = local_poles_array[2][2][0]; local_poles_and_weights_array[2][2][1] = local_poles_array[2][2][1]; local_poles_and_weights_array[2][2][2] = local_poles_array[2][2][2]; local_poles_and_weights_array[0][0][3] = local_weights_array[0][0]; local_poles_and_weights_array[0][1][3] = local_weights_array[0][1]; local_poles_and_weights_array[0][2][3] = local_weights_array[0][2]; local_poles_and_weights_array[1][0][3] = local_weights_array[1][0]; local_poles_and_weights_array[1][1][3] = local_weights_array[1][1]; local_poles_and_weights_array[1][2][3] = local_weights_array[1][2]; local_poles_and_weights_array[2][0][3] = local_weights_array[2][0]; local_poles_and_weights_array[2][1][3] = local_weights_array[2][1]; local_poles_and_weights_array[2][2][3] = local_weights_array[2][2]; BSplSLib::RationalDerivative(2, 2, 2, 2, local_poles_and_weights_array[0][0][0], local_poles_array[0][0][0]) ; } Standard_Real minmin = inverse_min * inverse_min; Standard_Real minmax = inverse_min * inverse_max; Standard_Real maxmax = inverse_max * inverse_max; my_point [0] = local_poles_array [0][0][0] ; my_vec_min [0] = inverse_min * local_poles_array[0][1][0] ; my_vec_max [0] = inverse_max * local_poles_array[1][0][0] ; my_vec_min_min[0] = minmin * local_poles_array [0][2][0] ; my_vec_min_max[0] = minmax * local_poles_array [1][1][0] ; my_vec_max_max[0] = maxmax * local_poles_array [2][0][0] ; my_point [1] = local_poles_array [0][0][1] ; my_vec_min [1] = inverse_min * local_poles_array[0][1][1] ; my_vec_max [1] = inverse_max * local_poles_array[1][0][1] ; my_vec_min_min[1] = minmin * local_poles_array [0][2][1] ; my_vec_min_max[1] = minmax * local_poles_array [1][1][1] ; my_vec_max_max[1] = maxmax * local_poles_array [2][0][1] ; my_point [2] = local_poles_array [0][0][2] ; my_vec_min [2] = inverse_min * local_poles_array[0][1][2] ; my_vec_max [2] = inverse_max * local_poles_array[1][0][2] ; my_vec_min_min[2] = minmin * local_poles_array [0][2][2] ; my_vec_min_max[2] = minmax * local_poles_array [1][1][2] ; my_vec_max_max[2] = maxmax * local_poles_array [2][0][2] ; } //======================================================================= //function : MovePoint //purpose : Find the new poles which allows an old point (with a // given u and v as parameters) to reach a new position //======================================================================= void BSplSLib::MovePoint (const Standard_Real U, const Standard_Real V, const gp_Vec& Displ, const Standard_Integer UIndex1, const Standard_Integer UIndex2, const Standard_Integer VIndex1, const Standard_Integer VIndex2, const Standard_Integer UDegree, const Standard_Integer VDegree, const Standard_Boolean Rational, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, const TColStd_Array1OfReal& UFlatKnots, const TColStd_Array1OfReal& VFlatKnots, Standard_Integer& UFirstIndex, Standard_Integer& ULastIndex, Standard_Integer& VFirstIndex, Standard_Integer& VLastIndex, TColgp_Array2OfPnt& NewPoles) { // calculate the UBSplineBasis in the parameter U Standard_Integer UFirstNonZeroBsplineIndex; math_Matrix UBSplineBasis(1, 1, 1, UDegree+1); Standard_Integer ErrorCod1 = BSplCLib::EvalBsplineBasis(1, 0, UDegree+1, UFlatKnots, U, UFirstNonZeroBsplineIndex, UBSplineBasis); // calculate the VBSplineBasis in the parameter V Standard_Integer VFirstNonZeroBsplineIndex; math_Matrix VBSplineBasis(1, 1, 1, VDegree+1); Standard_Integer ErrorCod2 = BSplCLib::EvalBsplineBasis(1, 0, VDegree+1, VFlatKnots, V, VFirstNonZeroBsplineIndex, VBSplineBasis); if (ErrorCod1 || ErrorCod2) { UFirstIndex = 0; ULastIndex = 0; VFirstIndex = 0; VLastIndex = 0; return; } // find the span which is predominant for parameter U UFirstIndex = UFirstNonZeroBsplineIndex; ULastIndex = UFirstNonZeroBsplineIndex + UDegree ; if (UFirstIndex < UIndex1) UFirstIndex = UIndex1; if (ULastIndex > UIndex2) ULastIndex = UIndex2; Standard_Real maxValue = 0.0; Standard_Integer i, ukk1=0, ukk2; for (i = UFirstIndex-UFirstNonZeroBsplineIndex+1; i <= ULastIndex-UFirstNonZeroBsplineIndex+1; i++) { if (UBSplineBasis(1,i) > maxValue) { ukk1 = i + UFirstNonZeroBsplineIndex - 1; maxValue = UBSplineBasis(1,i); } } // find a ukk2 if symetriy ukk2 = ukk1; i = ukk1 - UFirstNonZeroBsplineIndex + 2; if ((ukk1+1) <= ULastIndex) { if (Abs(UBSplineBasis(1, ukk1-UFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) { ukk2 = ukk1+1; } } // find the span which is predominant for parameter V VFirstIndex = VFirstNonZeroBsplineIndex; VLastIndex = VFirstNonZeroBsplineIndex + VDegree ; if (VFirstIndex < VIndex1) VFirstIndex = VIndex1; if (VLastIndex > VIndex2) VLastIndex = VIndex2; maxValue = 0.0; Standard_Integer j, vkk1=0, vkk2; for (j = VFirstIndex-VFirstNonZeroBsplineIndex+1; j <= VLastIndex-VFirstNonZeroBsplineIndex+1; j++) { if (VBSplineBasis(1,j) > maxValue) { vkk1 = j + VFirstNonZeroBsplineIndex - 1; maxValue = VBSplineBasis(1,j); } } // find a vkk2 if symetriy vkk2 = vkk1; j = vkk1 - VFirstNonZeroBsplineIndex + 2; if ((vkk1+1) <= VLastIndex) { if (Abs(VBSplineBasis(1, vkk1-VFirstNonZeroBsplineIndex+2) - maxValue) < 1.e-10) { vkk2 = vkk1+1; } } // compute the vector of displacement Standard_Real D1 = 0.0; Standard_Real D2 = 0.0; Standard_Real hN, Coef, DvalU, DvalV; Standard_Integer ii, jj; for (i = 1; i <= UDegree+1; i++) { ii = i + UFirstNonZeroBsplineIndex - 1; if (ii < ukk1) { DvalU = ukk1-ii; } else if (ii > ukk2) { DvalU = ii - ukk2; } else { DvalU = 0.0; } for (j = 1; j <= VDegree+1; j++) { jj = j + VFirstNonZeroBsplineIndex - 1; if (Rational) { hN = Weights(ii, jj)*UBSplineBasis(1, i)*VBSplineBasis(1,j); D2 += hN; } else { hN = UBSplineBasis(1, i)*VBSplineBasis(1,j); } if (ii >= UFirstIndex && ii <= ULastIndex && jj >= VFirstIndex && jj <= VLastIndex) { if (jj < vkk1) { DvalV = vkk1-jj; } else if (jj > vkk2) { DvalV = jj - vkk2; } else { DvalV = 0.0; } D1 += 1./(DvalU + DvalV + 1.) * hN; } } } if (Rational) { Coef = D2/D1; } else { Coef = 1./D1; } // compute the new poles for (i=Poles.LowerRow(); i<=Poles.UpperRow(); i++) { if (i < ukk1) { DvalU = ukk1-i; } else if (i > ukk2) { DvalU = i - ukk2; } else { DvalU = 0.0; } for (j=Poles.LowerCol(); j<=Poles.UpperCol(); j++) { if (i >= UFirstIndex && i <= ULastIndex && j >= VFirstIndex && j <= VLastIndex) { if (j < vkk1) { DvalV = vkk1-j; } else if (j > vkk2) { DvalV = j - vkk2; } else { DvalV = 0.0; } NewPoles(i,j) = Poles(i,j).Translated((Coef/(DvalU + DvalV + 1.))*Displ); } else { NewPoles(i,j) = Poles(i,j); } } } } //======================================================================= // function : Resolution // purpose : this computes an estimate for the maximum of the // partial derivatives both in U and in V // // // The calculation resembles at the calculation of curves with // additional index for the control point. Let Si,j be the // control points for ls surface and Di,j the weights. // The checking of upper bounds for the partial derivatives // will be omitted and Su is the next upper bound in the polynomial case : // // // // | Si,j - Si-1,j | // d * Max | ------------- | // i = 2,n | ti+d - ti | // i=1.m // // // and in the rational case : // // // // Di,j * (Si,j - Sk,j) - Di-1,j * (Si-1,j - Sk,j) // Max Max d * ----------------------------------------------- // k=1,n i dans Rj ti+d - ti // j=1,m // ---------------------------------------------------------------------- // // Min Di,j // i=1,n // j=1,m // // // // with Rj = {j-d, ...., j+d+d+1}. // // //======================================================================= void BSplSLib::Resolution(const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, const TColStd_Array1OfReal& UKnots, const TColStd_Array1OfReal& VKnots, const TColStd_Array1OfInteger& UMults, const TColStd_Array1OfInteger& VMults, const Standard_Integer UDegree, const Standard_Integer VDegree, const Standard_Boolean URational, const Standard_Boolean VRational, const Standard_Boolean UPeriodic, const Standard_Boolean VPeriodic, const Standard_Real Tolerance3D, Standard_Real& UTolerance, Standard_Real& VTolerance) { Standard_Real Wij,Wmj,Wji,Wjm; Standard_Real Xij,Xmj,Xji,Xjm,Xpq,Xqp; Standard_Real Yij,Ymj,Yji,Yjm,Ypq,Yqp; Standard_Real Zij,Zmj,Zji,Zjm,Zpq,Zqp; Standard_Real factor,value,min,min_weights=0,inverse,max_derivative[2]; max_derivative[0] = max_derivative[1] = 0.0e0 ; Standard_Integer PRowLength, PColLength; Standard_Integer ii,jj,pp,qq,ii_index,jj_index,pp_index,qq_index; Standard_Integer ii_minus,upper[2],lower[2],poles_length[2]; Standard_Integer num_poles[2],num_flat_knots[2]; num_flat_knots[0] = BSplCLib::KnotSequenceLength(UMults, UDegree, UPeriodic) ; num_flat_knots[1] = BSplCLib::KnotSequenceLength(VMults, VDegree, VPeriodic) ; TColStd_Array1OfReal flat_knots_in_u(1,num_flat_knots[0]) ; TColStd_Array1OfReal flat_knots_in_v(1,num_flat_knots[1]) ; BSplCLib::KnotSequence(UKnots, UMults, UDegree, UPeriodic, flat_knots_in_u) ; BSplCLib::KnotSequence(VKnots, VMults, VDegree, VPeriodic, flat_knots_in_v) ; PRowLength = Poles.RowLength(); PColLength = Poles.ColLength(); if (URational || VRational) { Standard_Integer Wsize = PRowLength * PColLength; const Standard_Real * WG = &Weights(Weights.LowerRow(),Weights.LowerCol()); min_weights = WG[0]; for (ii = 1 ; ii < Wsize ; ii++) { min = WG[ii]; if (min_weights > min) min_weights = min; } } Standard_Integer UD1 = UDegree + 1; Standard_Integer VD1 = VDegree + 1; num_poles[0] = num_flat_knots[0] - UD1; num_poles[1] = num_flat_knots[1] - VD1; poles_length[0] = PColLength; poles_length[1] = PRowLength; if(URational) { Standard_Integer UD2 = UDegree << 1; Standard_Integer VD2 = VDegree << 1; for (ii = 2 ; ii <= num_poles[0] ; ii++) { ii_index = (ii - 1) % poles_length[0] + 1 ; ii_minus = (ii - 2) % poles_length[0] + 1 ; inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ; inverse = 1.0e0 / inverse ; lower[0] = ii - UD1; if (lower[0] < 1) lower[0] = 1; upper[0] = ii + UD2 + 1; if (upper[0] > num_poles[0]) upper[0] = num_poles[0]; for ( jj = 1 ; jj <= num_poles[1] ; jj++) { jj_index = (jj - 1) % poles_length[1] + 1 ; lower[1] = jj - VD1; if (lower[1] < 1) lower[1] = 1; upper[1] = jj + VD2 + 1; if (upper[1] > num_poles[1]) upper[1] = num_poles[1]; const gp_Pnt& Pij = Poles .Value(ii_index,jj_index); Wij = Weights.Value(ii_index,jj_index); const gp_Pnt& Pmj = Poles .Value(ii_minus,jj_index); Wmj = Weights.Value(ii_minus,jj_index); Xij = Pij.X(); Yij = Pij.Y(); Zij = Pij.Z(); Xmj = Pmj.X(); Ymj = Pmj.Y(); Zmj = Pmj.Z(); for (pp = lower[0] ; pp <= upper[0] ; pp++) { pp_index = (pp - 1) % poles_length[0] + 1 ; for (qq = lower[1] ; qq <= upper[1] ; qq++) { value = 0.0e0 ; qq_index = (qq - 1) % poles_length[1] + 1 ; const gp_Pnt& Ppq = Poles.Value(pp_index,qq_index); Xpq = Ppq.X(); Ypq = Ppq.Y(); Zpq = Ppq.Z(); factor = (Xpq - Xij) * Wij; factor -= (Xpq - Xmj) * Wmj; if (factor < 0) factor = - factor; value += factor ; factor = (Ypq - Yij) * Wij; factor -= (Ypq - Ymj) * Wmj; if (factor < 0) factor = - factor; value += factor ; factor = (Zpq - Zij) * Wij; factor -= (Zpq - Zmj) * Wmj; if (factor < 0) factor = - factor; value += factor ; value *= inverse ; if (max_derivative[0] < value) max_derivative[0] = value ; } } } } max_derivative[0] /= min_weights ; } else { for (ii = 2 ; ii <= num_poles[0] ; ii++) { ii_index = (ii - 1) % poles_length[0] + 1 ; ii_minus = (ii - 2) % poles_length[0] + 1 ; inverse = flat_knots_in_u(ii + UDegree) - flat_knots_in_u(ii) ; inverse = 1.0e0 / inverse ; for ( jj = 1 ; jj <= num_poles[1] ; jj++) { jj_index = (jj - 1) % poles_length[1] + 1 ; value = 0.0e0 ; const gp_Pnt& Pij = Poles.Value(ii_index,jj_index); const gp_Pnt& Pmj = Poles.Value(ii_minus,jj_index); factor = Pij.X() - Pmj.X(); if (factor < 0) factor = - factor; value += factor; factor = Pij.Y() - Pmj.Y(); if (factor < 0) factor = - factor; value += factor; factor = Pij.Z() - Pmj.Z(); if (factor < 0) factor = - factor; value += factor; value *= inverse ; if (max_derivative[0] < value) max_derivative[0] = value ; } } } max_derivative[0] *= UDegree ; if(VRational) { Standard_Integer UD2 = UDegree << 1; Standard_Integer VD2 = VDegree << 1; for (ii = 2 ; ii <= num_poles[1] ; ii++) { ii_index = (ii - 1) % poles_length[1] + 1 ; ii_minus = (ii - 2) % poles_length[1] + 1 ; inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ; inverse = 1.0e0 / inverse ; lower[0] = ii - VD1; if (lower[0] < 1) lower[0] = 1; upper[0] = ii + VD2 + 1; if (upper[0] > num_poles[1]) upper[0] = num_poles[1]; for ( jj = 1 ; jj <= num_poles[0] ; jj++) { jj_index = (jj - 1) % poles_length[0] + 1 ; lower[1] = jj - UD1; if (lower[1] < 1) lower[1] = 1; upper[1] = jj + UD2 + 1; if (upper[1] > num_poles[0]) upper[1] = num_poles[0]; const gp_Pnt& Pji = Poles .Value(jj_index,ii_index); Wji = Weights.Value(jj_index,ii_index); const gp_Pnt& Pjm = Poles .Value(jj_index,ii_minus); Wjm = Weights.Value(jj_index,ii_minus); Xji = Pji.X(); Yji = Pji.Y(); Zji = Pji.Z(); Xjm = Pjm.X(); Yjm = Pjm.Y(); Zjm = Pjm.Z(); for (pp = lower[1] ; pp <= upper[1] ; pp++) { pp_index = (pp - 1) % poles_length[1] + 1 ; for (qq = lower[0] ; qq <= upper[0] ; qq++) { value = 0.0e0 ; qq_index = (qq - 1) % poles_length[0] + 1 ; const gp_Pnt& Pqp = Poles.Value(qq_index,pp_index); Xqp = Pqp.X(); Yqp = Pqp.Y(); Zqp = Pqp.Z(); factor = (Xqp - Xji) * Wji; factor -= (Xqp - Xjm) * Wjm; if (factor < 0) factor = - factor; value += factor ; factor = (Yqp - Yji) * Wji; factor -= (Yqp - Yjm) * Wjm; if (factor < 0) factor = - factor; value += factor ; factor = (Zqp - Zji) * Wji; factor -= (Zqp - Zjm) * Wjm; if (factor < 0) factor = - factor; value += factor ; value *= inverse ; if (max_derivative[1] < value) max_derivative[1] = value ; } } } } max_derivative[1] /= min_weights ; } else { for (ii = 2 ; ii <= num_poles[1] ; ii++) { ii_index = (ii - 1) % poles_length[1] + 1 ; ii_minus = (ii - 2) % poles_length[1] + 1 ; inverse = flat_knots_in_v(ii + VDegree) - flat_knots_in_v(ii) ; inverse = 1.0e0 / inverse ; for ( jj = 1 ; jj <= num_poles[0] ; jj++) { jj_index = (jj - 1) % poles_length[0] + 1 ; value = 0.0e0 ; const gp_Pnt& Pji = Poles.Value(jj_index,ii_index); const gp_Pnt& Pjm = Poles.Value(jj_index,ii_minus); factor = Pji.X() - Pjm.X() ; if (factor < 0) factor = - factor; value += factor; factor = Pji.Y() - Pjm.Y() ; if (factor < 0) factor = - factor; value += factor; factor = Pji.Z() - Pjm.Z() ; if (factor < 0) factor = - factor; value += factor; value *= inverse ; if (max_derivative[1] < value) max_derivative[1] = value ; } } } max_derivative[1] *= VDegree ; max_derivative[0] *= M_SQRT2 ; max_derivative[1] *= M_SQRT2 ; if(max_derivative[0] && max_derivative[1]) { UTolerance = Tolerance3D / max_derivative[0] ; VTolerance = Tolerance3D / max_derivative[1] ; } else { UTolerance=VTolerance=0.0; #ifdef OCCT_DEBUG cout<<"ElSLib.cxx : maxderivative = 0.0 "<Init(0); for (ii=1; ii <= VLength; ii++) { for (jj=1, ll=1; jj<= ULength; jj++, ll+=4) { Points(ii,ll) = Poles(jj, ii).X(); Points(ii,ll+1) = Poles(jj, ii).Y(); Points(ii,ll+2) = Poles(jj, ii).Z(); Points(ii,ll+3) = Weights(jj,ii) ; } } // interpolation of iso u poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ; BSplCLib::Interpolate(VDegree, VFlatKnots, VParameters, ContactOrder->Array1(), dimension, poles_array[0], InversionProblem) ; if (InversionProblem != 0) return; // extraction of iso v dimension = VLength*4; TColStd_Array2OfReal IsoPoles(1, ULength, 1, dimension); ContactOrder = new (TColStd_HArray1OfInteger)(1, ULength); ContactOrder->Init(0); poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ; for (ii=1, kk=1; ii <= ULength; ii++, kk+=4) { for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) { IsoPoles (ii,ll) = Points(jj, kk); IsoPoles (ii,ll+1) = Points(jj, kk+1); IsoPoles (ii,ll+2) = Points(jj, kk+2); IsoPoles (ii,ll+3) = Points(jj, kk+3); } } // interpolation of iso v BSplCLib::Interpolate(UDegree, UFlatKnots, UParameters, ContactOrder->Array1(), dimension, poles_array[0], InversionProblem); // return results for (ii=1; ii <= ULength; ii++) { for (jj=1, ll=1; jj<= VLength; jj++, ll+=4) { gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2)); Poles.SetValue(ii, jj, Pnt); Weights.SetValue(ii,jj,IsoPoles(ii,ll+3)) ; } } } //======================================================================= //function : Interpolate //purpose : //======================================================================= void BSplSLib::Interpolate(const Standard_Integer UDegree, const Standard_Integer VDegree, const TColStd_Array1OfReal& UFlatKnots, const TColStd_Array1OfReal& VFlatKnots, const TColStd_Array1OfReal& UParameters, const TColStd_Array1OfReal& VParameters, TColgp_Array2OfPnt& Poles, Standard_Integer& InversionProblem) { Standard_Integer ii, jj, ll, kk, dimension; Standard_Integer ULength = UParameters.Length(); Standard_Integer VLength = VParameters.Length(); Standard_Real * poles_array; // extraction of iso u dimension = 3*ULength; TColStd_Array2OfReal Points(1, VLength, 1, dimension); Handle(TColStd_HArray1OfInteger) ContactOrder = new (TColStd_HArray1OfInteger)(1, VLength); ContactOrder->Init(0); for (ii=1; ii <= VLength; ii++) { for (jj=1, ll=1; jj<= ULength; jj++, ll+=3) { Points(ii,ll) = Poles(jj, ii).X(); Points(ii,ll+1) = Poles(jj, ii).Y(); Points(ii,ll+2) = Poles(jj, ii).Z(); } } // interpolation of iso u poles_array = (Standard_Real *) &Points.ChangeValue(1,1) ; BSplCLib::Interpolate(VDegree, VFlatKnots, VParameters, ContactOrder->Array1(), dimension, poles_array[0], InversionProblem) ; if (InversionProblem != 0) return; // extraction of iso v dimension = VLength*3; TColStd_Array2OfReal IsoPoles(1, ULength, 1, dimension); ContactOrder = new (TColStd_HArray1OfInteger)(1, ULength); ContactOrder->Init(0); poles_array = (Standard_Real *) &IsoPoles.ChangeValue(1,1) ; for (ii=1, kk=1; ii <= ULength; ii++, kk+=3) { for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) { IsoPoles (ii,ll) = Points(jj, kk); IsoPoles (ii,ll+1) = Points(jj, kk+1); IsoPoles (ii,ll+2) = Points(jj, kk+2); } } // interpolation of iso v BSplCLib::Interpolate(UDegree, UFlatKnots, UParameters, ContactOrder->Array1(), dimension, poles_array[0], InversionProblem); // return results for (ii=1; ii <= ULength; ii++) { for (jj=1, ll=1; jj<= VLength; jj++, ll+=3) { gp_Pnt Pnt(IsoPoles(ii,ll), IsoPoles(ii,ll+1), IsoPoles(ii,ll+2)); Poles.SetValue(ii, jj, Pnt); } } } //======================================================================= //function : FunctionMultiply //purpose : //======================================================================= void BSplSLib::FunctionMultiply (const BSplSLib_EvaluatorFunction& Function, const Standard_Integer UBSplineDegree, const Standard_Integer VBSplineDegree, const TColStd_Array1OfReal& UBSplineKnots, const TColStd_Array1OfReal& VBSplineKnots, const TColStd_Array1OfInteger & UMults, const TColStd_Array1OfInteger & VMults, const TColgp_Array2OfPnt& Poles, const TColStd_Array2OfReal& Weights, const TColStd_Array1OfReal& UFlatKnots, const TColStd_Array1OfReal& VFlatKnots, const Standard_Integer UNewDegree, const Standard_Integer VNewDegree, TColgp_Array2OfPnt& NewNumerator, TColStd_Array2OfReal& NewDenominator, Standard_Integer& Status) { Standard_Integer num_uparameters, // ii,jj,kk, ii,jj, error_code, num_vparameters ; Standard_Real result ; num_uparameters = UFlatKnots.Length() - UNewDegree - 1 ; num_vparameters = VFlatKnots.Length() - VNewDegree - 1 ; TColStd_Array1OfReal UParameters(1,num_uparameters) ; TColStd_Array1OfReal VParameters(1,num_vparameters) ; if ((NewNumerator.ColLength() == num_uparameters) && (NewNumerator.RowLength() == num_vparameters) && (NewDenominator.ColLength() == num_uparameters) && (NewDenominator.RowLength() == num_vparameters)) { BSplCLib::BuildSchoenbergPoints(UNewDegree, UFlatKnots, UParameters) ; BSplCLib::BuildSchoenbergPoints(VNewDegree, VFlatKnots, VParameters) ; for (ii = 1 ; ii <= num_uparameters ; ii++) { for (jj = 1 ; jj <= num_vparameters ; jj++) { HomogeneousD0(UParameters(ii), VParameters(jj), 0, 0, Poles, Weights, UBSplineKnots, VBSplineKnots, UMults, VMults, UBSplineDegree, VBSplineDegree, Standard_True, Standard_True, Standard_False, Standard_False, NewDenominator(ii,jj), NewNumerator(ii,jj)) ; Function.Evaluate (0, UParameters(ii), VParameters(jj), result, error_code) ; if (error_code) { Standard_ConstructionError::Raise(); } gp_Pnt& P = NewNumerator(ii,jj); P.SetX(P.X() * result); P.SetY(P.Y() * result); P.SetZ(P.Z() * result); NewDenominator(ii,jj) *= result ; } } Interpolate(UNewDegree, VNewDegree, UFlatKnots, VFlatKnots, UParameters, VParameters, NewNumerator, NewDenominator, Status) ; } else { Standard_ConstructionError::Raise(); } }