1 // Created on: 1995-05-29
2 // Created by: Xavier BENVENISTE
3 // Copyright (c) 1995-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
17 // Modified PMN 22/05/1996 : Pb de recopie des poles
18 // Modified PMN 20/08/1996 : Introduction de la decoupe parametrable.
19 // (plus verif des derives en debug)
20 // Modified PMN 20/08/1996 : Linearisation de F(t) => Les sous Espaces
21 // facilement approchable ont de petites erreurs.
22 // Modified PMN 15/04/1997 : Gestion fine de la continuite aux lieux de decoupes
24 #include <AdvApprox_ApproxAFunction.hxx>
25 #include <AdvApprox_Cutting.hxx>
26 #include <AdvApprox_DichoCutting.hxx>
27 #include <AdvApprox_EvaluatorFunction.hxx>
28 #include <AdvApprox_SimpleApprox.hxx>
29 #include <BSplCLib.hxx>
30 #include <Convert_CompPolynomialToPoles.hxx>
31 #include <GeomAbs_Shape.hxx>
33 #include <gp_Vec2d.hxx>
34 #include <math_Vector.hxx>
36 #include <PLib_JacobiPolynomial.hxx>
37 #include <Precision.hxx>
38 #include <Standard_ConstructionError.hxx>
39 #include <Standard_OutOfRange.hxx>
40 #include <TColgp_Array1OfPnt.hxx>
41 #include <TColgp_Array1OfPnt2d.hxx>
42 #include <TColgp_HArray2OfPnt.hxx>
43 #include <TColgp_HArray2OfPnt2d.hxx>
44 #include <TColStd_Array1OfInteger.hxx>
45 #include <TColStd_Array1OfReal.hxx>
46 #include <TColStd_HArray1OfInteger.hxx>
47 #include <TColStd_HArray1OfReal.hxx>
48 #include <TColStd_HArray2OfReal.hxx>
52 static Standard_Boolean AdvApprox_Debug = 0;
54 //=====================================================
55 static void MAPDBN(const Standard_Integer dimension,
56 const Standard_Real Debut,
57 const Standard_Real Fin,
58 AdvApprox_EvaluatorFunction& Evaluator,
59 const Standard_Integer Iordre)
60 // Objet : Controle par difference finis, des derives
61 // Warning : En mode Debug, uniquement
62 ///===================================================
64 Standard_Integer derive, OrdreDer, ier, NDIMEN = dimension;
66 Standard_Real h = 1.e-4*(Fin-Debut+1.e-3), eps = 1.e-3, t, ab[2];
67 math_Vector V1(1,NDIMEN), V2(1,NDIMEN), Diff(1,NDIMEN), Der(1,NDIMEN);
69 if (h < 100*RealEpsilon()) {h = 100*RealEpsilon();}
73 for (OrdreDer=1, derive = 0;
74 OrdreDer <= Iordre; OrdreDer++, derive++) {
76 Ptr = (Standard_Real*) &V1.Value(1);
78 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
80 Ptr = (Standard_Real*) &V2.Value(1);
82 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
86 Ptr = (Standard_Real*) &Der.Value(1);
88 Evaluator(&NDIMEN, ab, &t, &OrdreDer, Ptr, &ier);
91 if ( (Diff-Der).Norm() > eps * (Der.Norm()+1) ) {
92 cout << " Debug Ft au parametre t+ = " << t << endl;
93 cout << " Positionement sur la derive "<< OrdreDer
94 << " : " << Der << endl;
95 cout << " Erreur estime : " << (Der-Diff) << endl;
99 Ptr = (Standard_Real*) &V1.Value(1);
101 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
103 Ptr = (Standard_Real*) &V2.Value(1);
105 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
109 Ptr = (Standard_Real*) &Der.Value(1);
111 Evaluator(&NDIMEN, ab, &t, &OrdreDer, Ptr, &ier);
114 if ( (Diff-Der).Norm() > eps * (Der.Norm()+1) ) {
115 cout << " Debug Ft au parametre t- = " << t << endl;
116 cout << " Positionement sur la derive "<< OrdreDer
117 << " : " << Der << endl;
118 cout << " Erreur estime : " << (Der-Diff) << endl;
124 //===================================================================
125 static void PrepareConvert(const Standard_Integer NumCurves,
126 const Standard_Integer MaxDegree,
127 const Standard_Integer ContinuityOrder,
128 const Standard_Integer Num1DSS,
129 const Standard_Integer Num2DSS,
130 const Standard_Integer Num3DSS,
131 const TColStd_Array1OfInteger& NumCoeffPerCurve,
132 TColStd_Array1OfReal& Coefficients,
133 const TColStd_Array2OfReal& PolynomialIntervals,
134 const TColStd_Array1OfReal& TrueIntervals,
135 const TColStd_Array1OfReal& LocalTolerance,
136 TColStd_Array1OfReal& ErrorMax,
137 TColStd_Array1OfInteger& Continuity)
138 // Pour determiner les continuites locales
139 //====================================================================
143 Standard_Boolean isCi;
144 Standard_Integer icurve, idim, iordre, ii,
145 Dimension=Num1DSS + 2*Num2DSS + 3*Num3DSS,
146 NbSpace = Num1DSS + Num2DSS + Num3DSS;
147 Standard_Real diff, moy, facteur1, facteur2, normal1, normal2, eps;
148 Standard_Real *Res1, *Res2, *Val1, *Val2;
149 Standard_Real *Coef1, *Coef2;
150 Standard_Integer RealDegree = Max(MaxDegree + 1, 2 * ContinuityOrder + 2);
155 TColStd_Array1OfReal Result(1, 2*(ContinuityOrder+1)*Dimension);
156 TColStd_Array1OfReal Prec(1, NbSpace), Suivant(1, NbSpace);
159 Res1 = (Standard_Real*) &(Result.ChangeValue(1));
160 Res2 = (Standard_Real*) &(Result.ChangeValue((ContinuityOrder+1)*Dimension+1));
165 if (ContinuityOrder ==0) return;
167 for (icurve=1; icurve<NumCurves; icurve++) {
168 // Init et positionement au noeud
169 isCi = Standard_True;
170 Coef1 = (Standard_Real*) &(Coefficients.Value(
171 (icurve-1) * Dimension * RealDegree + Coefficients.Lower()));
172 Coef2 = Coef1 + Dimension * RealDegree;
173 Standard_Integer Deg1 = NumCoeffPerCurve(NumCoeffPerCurve.Lower()+icurve-1) - 1;
174 PLib::EvalPolynomial(PolynomialIntervals(icurve, 2),
180 Standard_Integer Deg2 = NumCoeffPerCurve(NumCoeffPerCurve.Lower()+icurve) - 1;
181 PLib::EvalPolynomial(PolynomialIntervals(icurve+1, 1),
188 // Verif dans chaque sous espaces.
189 for (iordre=1; iordre<=ContinuityOrder && isCi; iordre++) {
191 // fixing a bug PRO18577
193 Standard_Real Toler = 1.0e-5;
195 Standard_Real f1_dividend = PolynomialIntervals(icurve,2)-PolynomialIntervals(icurve,1);
196 Standard_Real f2_dividend = PolynomialIntervals(icurve+1,2)-PolynomialIntervals(icurve+1,1);
197 Standard_Real f1_divizor = TrueIntervals(icurve+1)-TrueIntervals(icurve);
198 Standard_Real f2_divizor = TrueIntervals(icurve+2)-TrueIntervals(icurve+1);
199 Standard_Real fract1, fract2;
201 if( Abs(f1_divizor) < Toler ) // this is to avoid divizion by zero
202 //in this case fract1 = 5.14755758946803e-85
204 else{ fract1 = f1_dividend/f1_divizor;
205 facteur1 = Pow( fract1, iordre);
207 if( Abs(f2_divizor) < Toler )
208 //in this case fract2 = 6.77193633669143e-313
210 else{ fract2 = f2_dividend/f2_divizor;
211 facteur2 = Pow( fract2, iordre);
213 normal1 = Pow(f1_divizor, iordre);
214 normal2 = Pow(f2_divizor, iordre);
218 Val1 = Res1+iordre*Dimension-1;
219 Val2 = Res2+iordre*Dimension-1;
221 for (ii=1; ii<=Num1DSS && isCi; ii++, idim++) {
222 eps = LocalTolerance(idim)*0.01;
223 diff = Abs (Val1[ii]*facteur1 - Val2[ii]*facteur2);
224 moy = Abs (Val1[ii]*facteur1 + Val2[ii]*facteur2);
225 // Un premier controle sur la valeur relative
226 if (diff > moy*1.e-9) {
227 Prec(idim) = diff*normal1;
228 Suivant(idim) = diff*normal2;
229 // Et un second sur le majorant d'erreur engendree
230 if (Prec(idim)>eps ||
231 Suivant(idim)>eps) isCi=Standard_False;
241 for (ii=1; ii<=Num2DSS && isCi; ii++, idim++, Val1+=2,Val2+=2) {
242 eps = LocalTolerance(idim)*0.01;
243 v1.SetCoord(Val1[1],Val1[2]);
244 v2.SetCoord(Val2[1],Val2[2]);
247 diff = Abs(v1.X() - v2.X()) + Abs(v1.Y() - v2.Y());
248 moy = Abs(v1.X() + v2.X()) + Abs(v1.Y() + v2.Y());
249 if (diff > moy*1.e-9) {
250 Prec(idim) = diff*normal1;
251 Suivant(idim) = diff*normal2;
252 if (Prec(idim)>eps ||
253 Suivant(idim)>eps) isCi=Standard_False;
261 for (ii=1; ii<=Num3DSS && isCi; ii++, idim++, Val1+=3,Val2+=3) {
262 eps = LocalTolerance(idim)*0.01;
263 V1.SetCoord(Val1[1], Val1[2], Val1[3]);
264 V2.SetCoord(Val2[1], Val2[2], Val2[3]);
267 diff = Abs(V1.X() - V2.X()) + Abs(V1.Y() - V2.Y()) +
268 Abs(V1.Z() - V2.Z());
269 moy = Abs(V1.X() + V2.X()) + Abs(V1.Y() + V2.Y()) +
270 Abs(V1.Z() + V2.Z());
271 if (diff > moy*1.e-9) {
272 Prec(idim) = diff*normal1;
273 Suivant(idim) = diff*normal2;
274 if (Prec(idim)>eps ||
275 Suivant(idim)>eps) isCi=Standard_False;
282 // Si c'est bon on update le tout
284 Continuity(icurve+1) = iordre;
285 Standard_Integer index = (icurve-1)*NbSpace+1;
286 for (ii=index, idim=1; idim<=NbSpace; ii++,idim++) {
287 ErrorMax(ii) += Prec(idim);
288 ErrorMax(ii+NbSpace) += Suivant(idim);
296 //=======================================================================
297 //function : ApproxFunction
299 //purpose : Approximation d' UNE fonction non polynomiale (dans l'espace de
300 // dimension NDIMEN) par N courbes polynomiales, par la methode
301 // des moindres carres. Le parametre de la fonction est conserve.
303 // ARGUMENTS D'ENTREE :
304 //C ------------------
305 //C NDIMEN : Dimension totale de l' espace (somme des dimensions
306 //C des sous-espaces).
307 //C NBSESP : Nombre de sous-espaces "independants".
308 //C NDIMSE : Table des dimensions des sous-espaces.
309 //C ABFONC : Bornes de l' intervalle (a,b) de definition de la
310 //C fonction a approcher.
311 //C FONCNP : Fonction externe de positionnement sur la fonction non
312 //C polynomiale a approcher.
313 //C IORDRE : Ordre de contrainte aux extremites de chaque courbe
314 //C polynomiale cree :
315 //C -1 = pas de contraintes,
316 //C 0 = contraintes de passage aux bornes (i.e. C0),
317 //C 1 = C0 + contraintes de derivees 1eres (i.e. C1),
318 //C 2 = C1 + contraintes de derivees 2ndes (i.e. C2).
319 //C NBCRMX : Nombre maxi de courbes polynomiales a calculer.
320 //C NCFLIM : Nombre LIMITE de coeff des "courbes" polynomiales
321 //C d' approximation. Doit etre superieur ou egal a
322 //C 2*IORDRE + 2 et inferieur ou egal a 50 et a NCOFMX).
323 //C Limite a 14 apres l'appel a VRIRAZ pour eviter le bruit.
324 //C EPSAPR : Table des erreurs d' approximations souhaitees
325 //C sous-espace par sous-espace.
326 //C ICODEO : Code d' init. des parametres de discretisation.
327 //C (choix de NBPNT et de NDJAC).
328 //C = 1 calcul rapide avec precision moyenne sur les coeff.
329 //C = 2 calcul rapide avec meilleure precision " " ".
330 //C = 3 calcul un peu plus lent avec bonne precision ".
331 //C = 4 calcul lent avec la meilleure precision possible ".
334 //C ARGUMENTS DE SORTIE :
335 //C -------------------
336 //C NBCRBE : Nombre de courbes polynomiales creees.
337 //C NCFTAB : Table des nombres de coeff. significatifs des NBCRBE
338 //C "courbes" calculees.
339 //C CRBTAB : Tableau des coeff. des "courbes" polynomiales calculees.
340 //C Doit etre dimensionne a CRBTAB(NDIMEN,NCOFMX,NBCRMX).
341 //C TABINT : Table des NBCRBE + 1 bornes des intervalles de decoupe de
343 //C ERRMAX : Table des erreurs (sous-espace par sous espace)
344 //C MAXIMALES commises dans l' approximation de FONCNP par
345 //C par les NBCRBE courbes polynomiales.
346 //C ERRMOY : Table des erreurs (sous-espace par sous espace)
347 //C MOYENNES commises dans l' approximation de FONCNP par
348 //C par les NBCRBE courbes polynomiales.
349 //C IERCOD : Code d' erreur :
350 //C -1 = ERRMAX > EPSAPR pour au moins un des sous-espace.
351 //C On a alors NCRBMX courbes calculees, certaines de
352 //C degre mathematique min(NCFLIM,14)-1 ou la precision
353 //C demandee n' est pas atteinte.
354 //C -2 = Pb dans le mode DEBUG
356 //C 1 = Pb d' incoherence des entrees.
357 //C 10 = Pb de calcul de l' interpolation des contraintes.
358 //C 13 = Pb dans l' allocation dynamique.
359 //C 33 = Pb dans la recuperation des donnees du block data
360 //C des coeff. d' integration par la methode de GAUSS.
361 //C >100 Pb dans l' evaluation de FONCNP, le code d' erreur
362 //C renvoye est egal au code d' erreur de FONCNP + 100.
365 // --> La i-eme courbe polynomiale creee correspond a l' approximation
366 //C de FONCNP sur l' intervalle (TABINT(i),TABINT(i+1)). TABINT(i)
367 // est une suite STRICTEMENT monotone avec TABINT(1)=ABFONC(1) et
368 // TABINT(NBCRBE+1)=ABFONC(2).
370 // --> Si IERCOD=-1, la precision demandee n' est pas atteinte (ERRMAX
371 //C est superieur a EPSAPR sur au moins l' un des sous-espaces), mais
372 // on donne le meilleur resultat possible pour min(NCFLIM,14),NBCRMX
373 // et EPSAPR choisis par l' utilisateur Dans ce cas (ainsi que pour
374 //C IERCOD=0), on a une solution.
376 //C--> ATTENTION : Toute modification du calcul de NDJAC et NBPNT
377 //C doit etre reportee dans le programme MAPNBP.
379 // HISTORIQUE DES MODIFICATIONS :
380 // --------------------------------
382 // 20-08-1996 : PMN ; Creation d'apres la routine Fortran MAPFNC
383 //======================================================================
384 void AdvApprox_ApproxAFunction::Approximation(
385 const Standard_Integer TotalDimension,
386 // Dimension totale de l' espace
387 // (somme des dimensions des sous-espaces).
388 const Standard_Integer TotalNumSS,//Nombre de sous-espaces "independants".
389 const TColStd_Array1OfInteger& LocalDimension,//dimensions des sous-espaces.
390 const Standard_Real First,
391 const Standard_Real Last,
392 // Intervalle (a,b) de definition de la fonction a approcher.
393 AdvApprox_EvaluatorFunction& Evaluator,
394 // Fonction externe de positionnement sur la fonction a approcher.
395 const AdvApprox_Cutting& CutTool,
396 // Maniere de decouper en 2 un intervalle.
397 const Standard_Integer ContinuityOrder,
398 // ordre de continutie a respecter (-1, 0, 1, 2)
399 const Standard_Integer NumMaxCoeffs,
400 // Nombre LIMITE de coeff des "courbes" polynomiales
401 // d' approximation. Doit etre superieur ou egal a
402 // Doit etre superieur ou egal a 2*ContinuityOrder + 2
403 // Limite a 14 pour eviter le bruit.
404 const Standard_Integer MaxSegments,
405 //Nombre maxi de courbes polynomiales a calculer.
406 const TColStd_Array1OfReal& LocalTolerancesArray,
407 //Tolerances d'approximation souhaitees sous-espace par sous-espace.
408 const Standard_Integer code_precis,
409 //Code d' init. des parametres de discretisation.
410 //C (choix de NBPNT et de NDJAC).
411 //C = 1 calcul rapide avec precision moyenne sur les coeff.
412 //C = 2 calcul rapide avec meilleure precision " " ".
413 //C = 3 calcul un peu plus lent avec bonne precision ".
414 //C = 4 calcul lent avec la meilleure precision possible ".
415 Standard_Integer& NumCurves,
416 TColStd_Array1OfInteger& NumCoeffPerCurveArray,
417 TColStd_Array1OfReal& CoefficientArray,
418 TColStd_Array1OfReal& IntervalsArray,
419 TColStd_Array1OfReal& ErrorMaxArray,
420 TColStd_Array1OfReal& AverageErrorArray,
421 Standard_Integer& ErrorCode)
423 // Standard_Real EpsPar = Precision::Confusion();
424 Standard_Integer IDIM, NUPIL,TheDeg;
426 Standard_Integer NDIMEN = TotalDimension;
428 Standard_Boolean isCut = Standard_False;
430 // Defintion des Tableaux C
431 Standard_Real * TABINT = (Standard_Real*) &IntervalsArray(1);
435 CoefficientArray.Init(0);
437 //-------------------------- Verification des entrees ------------------
439 if ((MaxSegments<1)|| (Abs(Last-First)<1.e-9))
443 //--> La dimension totale doit etre la somme des dimensions des
446 for ( Standard_Integer I=1; I<=TotalNumSS; I++) {IDIM += LocalDimension(I);}
447 if (IDIM != TotalDimension)
450 GeomAbs_Shape Continuity=GeomAbs_C0;
451 switch (ContinuityOrder) {
453 Continuity = GeomAbs_C0 ;
456 Continuity = GeomAbs_C1 ;
459 Continuity = GeomAbs_C2 ;
462 Standard_ConstructionError::Raise();
465 //--------------------- Choix du nombre de points ----------------------
466 //---------- et du degre de lissage dans la base orthogonale -----------
467 //--> NDJAC est le degre de "travail" dans la base orthogonale.
470 Standard_Integer NbGaussPoints, WorkDegree;
472 PLib::JacobiParameters(Continuity, NumMaxCoeffs-1, code_precis,
473 NbGaussPoints, WorkDegree);
476 //------------------ Initialisation de la gestion des decoupes ---------
483 //C**********************************************************************
484 //C APPROXIMATION AVEC DECOUPES
485 //C**********************************************************************
486 Handle(PLib_JacobiPolynomial) JacobiBase = new (PLib_JacobiPolynomial) (WorkDegree, Continuity);
487 //Portage HP le compilateur refuse le debranchement
488 Standard_Integer IS ;
489 Standard_Boolean goto_fin_de_boucle;
490 Standard_Integer MaxDegree = NumMaxCoeffs-1;
491 AdvApprox_SimpleApprox Approx (TotalDimension, TotalNumSS,
493 WorkDegree, NbGaussPoints,
494 JacobiBase, Evaluator);
495 while ((NUPIL-NumCurves)!=0) {
496 //--> Lorsque l' on a atteint le haut de la pile, c' est fini !
498 //Portage HP le compilateur refuse le debranchement
499 goto_fin_de_boucle=Standard_False;
501 //---- Calcul de la courbe d' approximation dans la base de Jacobi -----
503 Approx.Perform ( LocalDimension, LocalTolerancesArray,
504 TABINT[NumCurves], TABINT[NumCurves+1],
506 if (!Approx.IsDone()) {
513 //---------- Calcul du degre de la courbe et de l' erreur max ----------
516 NumCoeffPerCurveArray(NumCurves + 1)=0;
518 // L'erreur doit etre satisfaite sur tous les sous-espaces sinon, on decoupe
520 Standard_Boolean MaxMaxErr=Standard_True;
521 for ( IS=1; IS<=TotalNumSS; IS++)
522 { if (Approx.MaxError(IS) > LocalTolerancesArray(IS))
523 { MaxMaxErr = Standard_False;
528 if (MaxMaxErr == Standard_True)
531 // NumCoeffPerCurveArray(NumCurves) = Approx.Degree()+1;
535 //-> ...sinon on essai de decouper l' intervalle courant en 2...
537 Standard_Boolean Large;
539 Large = CutTool.Value(TABINT[NumCurves], TABINT[NumCurves+1],
542 if (NUPIL<MaxSegments && Large) {
544 // if (IS!=1) {NumCurves--;}
545 isCut = Standard_True; // Ca y est, on le sait !
546 Standard_Real * IDEB = TABINT+NumCurves+1;
547 Standard_Real * IDEB1 = IDEB+1;
548 Standard_Integer ILONG= NUPIL-NumCurves-1;
549 for (Standard_Integer iI=ILONG; iI>=0; iI--) {
550 IDEB1[iI] = IDEB[iI];
554 //--> ... et on recommence.
555 //Portage HP le compilateur refuse le debranchement
556 goto_fin_de_boucle=Standard_True;
560 //--> Si la pile est pleine...
561 // pour l'instant ErrorCode=-1;
563 // NumCoeffPerCurveArray(NumCurves) = Approx.Degree()+1;
567 //Portage HP le compilateur refuse le debranchement
568 if (goto_fin_de_boucle) continue;
569 for (IS=1; IS<=TotalNumSS; IS++) {
570 ErrorMaxArray.SetValue(IS+(NumCurves-1)*TotalNumSS,Approx. MaxError(IS));
571 //---> ...et calcul de l' erreur moyenne.
572 AverageErrorArray.SetValue(IS+(NumCurves-1)*TotalNumSS,Approx.AverageError(IS));
575 Handle(TColStd_HArray1OfReal) HJacCoeff = Approx.Coefficients();
576 TheDeg = Approx.Degree();
577 if (isCut && (TheDeg < 2*ContinuityOrder+1) )
578 // Pour ne pas bruiter les derives aux bouts, et garder une continuite
579 // correcte sur la BSpline resultat.
580 TheDeg = 2*ContinuityOrder+1;
582 NumCoeffPerCurveArray(NumCurves) = TheDeg+1;
583 TColStd_Array1OfReal Coefficients(0,(TheDeg+1)*TotalDimension-1);
584 JacobiBase->ToCoefficients (TotalDimension, TheDeg,
585 HJacCoeff->Array1(), Coefficients);
586 Standard_Integer i,j, f = (TheDeg+1)*TotalDimension;
587 for (i=0,j=(NumCurves-1)*TotalDimension*NumMaxCoeffs+1;
589 CoefficientArray.SetValue(j, Coefficients.Value(i));
593 // Test des derives par difference finis
594 Standard_Integer IORDRE = ContinuityOrder;
596 if (IORDRE>0 && AdvApprox_Debug) {
597 MAPDBN(NDIMEN, TABINT[NumCurves-1],
598 TABINT[NumCurves], Evaluator, IORDRE);
601 //Portage HP le compilateur refuse le debranchement
603 // {;} fin de la boucle while
607 //=======================================================================
608 //function : AdvApprox_ApproxAFunction
609 //purpose : Constructeur avec Decoupe Dichotomique
610 //=======================================================================
612 AdvApprox_ApproxAFunction::
613 AdvApprox_ApproxAFunction(const Standard_Integer Num1DSS,
614 const Standard_Integer Num2DSS,
615 const Standard_Integer Num3DSS,
616 const Handle(TColStd_HArray1OfReal)& OneDTol,
617 const Handle(TColStd_HArray1OfReal)& TwoDTol,
618 const Handle(TColStd_HArray1OfReal)& ThreeDTol,
619 const Standard_Real First,
620 const Standard_Real Last,
621 const GeomAbs_Shape Continuity,
622 const Standard_Integer MaxDeg,
623 const Standard_Integer MaxSeg,
624 const AdvApprox_EvaluatorFunction& Func) :
625 my1DTolerances(OneDTol),
626 my2DTolerances(TwoDTol),
627 my3DTolerances(ThreeDTol),
630 myContinuity(Continuity),
632 myMaxSegments(MaxSeg),
633 myDone(Standard_False),
634 myHasResult(Standard_False),
635 myEvaluator((Standard_Address)&Func)
637 AdvApprox_DichoCutting Cut;
638 Perform(Num1DSS, Num2DSS, Num3DSS, Cut);
641 AdvApprox_ApproxAFunction::
642 AdvApprox_ApproxAFunction(const Standard_Integer Num1DSS,
643 const Standard_Integer Num2DSS,
644 const Standard_Integer Num3DSS,
645 const Handle(TColStd_HArray1OfReal)& OneDTol,
646 const Handle(TColStd_HArray1OfReal)& TwoDTol,
647 const Handle(TColStd_HArray1OfReal)& ThreeDTol,
648 const Standard_Real First,
649 const Standard_Real Last,
650 const GeomAbs_Shape Continuity,
651 const Standard_Integer MaxDeg,
652 const Standard_Integer MaxSeg,
653 const AdvApprox_EvaluatorFunction& Func,
654 const AdvApprox_Cutting& CutTool) :
655 my1DTolerances(OneDTol),
656 my2DTolerances(TwoDTol),
657 my3DTolerances(ThreeDTol),
660 myContinuity(Continuity),
662 myMaxSegments(MaxSeg),
663 myDone(Standard_False),
664 myHasResult(Standard_False),
665 myEvaluator((Standard_Address)&Func)
667 Perform(Num1DSS, Num2DSS, Num3DSS, CutTool);
670 //=======================================================================
671 //function : AdvApprox_ApproxAFunction::Perform
672 //purpose : Make all the Work !!
673 //=======================================================================
674 void AdvApprox_ApproxAFunction::Perform(const Standard_Integer Num1DSS,
675 const Standard_Integer Num2DSS,
676 const Standard_Integer Num3DSS,
677 const AdvApprox_Cutting& CutTool)
682 Num1DSS + Num2DSS + Num3DSS <= 0 ||
686 Standard_ConstructionError::Raise();
687 if (myMaxDegree > 14) {
693 myNumSubSpaces[0] = Num1DSS ;
694 myNumSubSpaces[1] = Num2DSS ;
695 myNumSubSpaces[2] = Num3DSS ;
696 Standard_Integer TotalNumSS =
697 Num1DSS + Num2DSS + Num3DSS,
704 Standard_Integer TotalDimension =
705 myNumSubSpaces[0] + 2 * myNumSubSpaces[1] + 3 * myNumSubSpaces[2] ;
706 Standard_Real error_value ;
708 Standard_Integer ContinuityOrder=0 ;
709 switch (myContinuity) {
711 ContinuityOrder = 0 ;
714 ContinuityOrder = 1 ;
717 ContinuityOrder = 2 ;
720 Standard_ConstructionError::Raise();
722 Standard_Real ApproxStartEnd[2] ;
723 Standard_Integer NumMaxCoeffs = Max(myMaxDegree + 1, 2 * ContinuityOrder + 2);
724 myMaxDegree = NumMaxCoeffs - 1;
725 Standard_Integer code_precis = 1 ;
727 // WARNING : the polynomial coefficients are the
728 // taylor expansion of the polynomial at 0.0e0 !
730 ApproxStartEnd[0] = -1.0e0 ;
731 ApproxStartEnd[1] = 1.0e0 ;
733 TColStd_Array1OfInteger LocalDimension(1,TotalNumSS) ;
737 TColStd_Array1OfReal LocalTolerances(1,TotalNumSS) ;
739 for(jj = 1; jj <= myNumSubSpaces[0] ; jj++) {
740 LocalTolerances.SetValue(index,my1DTolerances->Value(jj)) ;
741 LocalDimension.SetValue(index,1) ;
744 for(jj = 1 ; jj <= myNumSubSpaces[1] ; jj++) {
745 LocalTolerances.SetValue(index,my2DTolerances->Value(jj)) ;
746 LocalDimension.SetValue(index,2) ;
749 for(jj = 1; jj <= myNumSubSpaces[2] ; jj++) {
750 LocalTolerances.SetValue(index,my3DTolerances->Value(jj)) ;
751 LocalDimension.SetValue(index,3) ;
758 Standard_Integer ErrorCode = 0,
761 myMaxSegments * NumMaxCoeffs * TotalDimension ;
762 Handle(TColStd_HArray1OfInteger) NumCoeffPerCurvePtr =
763 new TColStd_HArray1OfInteger (1,myMaxSegments) ;
765 Handle(TColStd_HArray1OfReal) LocalCoefficientsPtr =
766 new TColStd_HArray1OfReal(1,size) ;
768 Handle(TColStd_HArray1OfReal) IntervalsPtr =
769 new TColStd_HArray1OfReal (1,myMaxSegments+1) ;
771 TColStd_Array1OfReal ErrorMax(1,myMaxSegments * TotalNumSS) ;
773 TColStd_Array1OfReal AverageError(1,myMaxSegments * TotalNumSS) ;
775 Approximation ( TotalDimension,
776 TotalNumSS, LocalDimension,
778 *(AdvApprox_EvaluatorFunction*)myEvaluator,
779 CutTool, ContinuityOrder, NumMaxCoeffs,
780 myMaxSegments, LocalTolerances, code_precis,
781 NumCurves, // Nombre de courbe en sortie
782 NumCoeffPerCurvePtr->ChangeArray1(), // Nbre de coeff par courbe
783 LocalCoefficientsPtr->ChangeArray1(), // Les Coeffs solutions
784 IntervalsPtr->ChangeArray1(), // La Table des decoupes
785 ErrorMax, // Majoration de l'erreur
786 AverageError, // Erreur moyenne constatee
789 if (ErrorCode == 0 || ErrorCode == -1) {
791 // si tout est OK ou bien on a un resultat dont l une des erreurs max est
792 // plus grande que la tolerance demandee
794 TColStd_Array1OfInteger TabContinuity(1, NumCurves);
795 TColStd_Array2OfReal PolynomialIntervalsPtr (1,NumCurves,1,2);
796 for (ii = PolynomialIntervalsPtr.LowerRow() ;
797 ii <= PolynomialIntervalsPtr.UpperRow() ;
799 // On force un degre 1 minimum (PRO5474)
800 NumCoeffPerCurvePtr->ChangeValue(ii) =
801 Max(2, NumCoeffPerCurvePtr->Value(ii));
802 PolynomialIntervalsPtr.SetValue(ii,1,ApproxStartEnd[0]) ;
803 PolynomialIntervalsPtr.SetValue(ii,2,ApproxStartEnd[1]) ;
806 PrepareConvert(NumCurves, myMaxDegree, ContinuityOrder,
807 Num1DSS, Num2DSS, Num3DSS,
808 NumCoeffPerCurvePtr->Array1(),
809 LocalCoefficientsPtr->ChangeArray1(),
810 PolynomialIntervalsPtr, IntervalsPtr->Array1(),
811 LocalTolerances, ErrorMax,
814 Convert_CompPolynomialToPoles
815 AConverter(NumCurves,
819 NumCoeffPerCurvePtr->Array1(),
820 LocalCoefficientsPtr->Array1(),
821 PolynomialIntervalsPtr,
822 IntervalsPtr->Array1());
824 if (AConverter.IsDone()) {
825 Handle(TColStd_HArray2OfReal) PolesPtr ;
826 AConverter.Poles(PolesPtr) ;
827 AConverter.Knots(myKnots) ;
828 AConverter.Multiplicities(myMults) ;
829 myDegree = AConverter.Degree() ;
831 if (myNumSubSpaces[0] > 0) {
832 my1DPoles = new TColStd_HArray2OfReal(1,PolesPtr->ColLength(),
833 1,myNumSubSpaces[0]) ;
834 my1DMaxError = new TColStd_HArray1OfReal(1,myNumSubSpaces[0]) ;
835 my1DAverageError = new TColStd_HArray1OfReal(1,myNumSubSpaces[0]) ;
836 for (ii = 1 ; ii <= PolesPtr->ColLength() ; ii++) {
837 for (jj = 1 ; jj <= myNumSubSpaces[0] ; jj++) {
838 my1DPoles->SetValue(ii,jj, PolesPtr->Value(ii,jj)) ;
842 for (jj = 1 ; jj <= myNumSubSpaces[0] ; jj++) {
843 error_value = 0.0e0 ;
844 for (ii = 1 ; ii <= NumCurves ; ii++) {
845 local_index = (ii - 1) * TotalNumSS ;
846 error_value = Max(ErrorMax(local_index + jj),error_value) ;
849 my1DMaxError->SetValue(jj, error_value) ;
851 for (jj = 1 ; jj <= myNumSubSpaces[0] ; jj++) {
852 error_value = 0.0e0 ;
853 for (ii = 1 ; ii <= NumCurves ; ii++) {
854 local_index = (ii - 1) * TotalNumSS ;
855 error_value += AverageError(local_index + jj) ;
858 error_value /= (Standard_Real) NumCurves ;
859 my1DAverageError->SetValue(jj, error_value) ;
861 index += myNumSubSpaces[0] ;
864 dim_index = index; //Pour le cas ou il n'y a pas de 2D
866 if (myNumSubSpaces[1] > 0) {
868 my2DPoles = new TColgp_HArray2OfPnt2d(1,PolesPtr->ColLength(),
869 1,myNumSubSpaces[1]) ;
870 my2DMaxError = new TColStd_HArray1OfReal(1,myNumSubSpaces[1]) ;
871 my2DAverageError = new TColStd_HArray1OfReal(1,myNumSubSpaces[1]) ;
872 for (ii = 1 ; ii <= PolesPtr->ColLength() ; ii++) {
873 for (jj = 1 ; jj <= myNumSubSpaces[1] ; jj++) {
874 local_index = index + (jj-1) * 2 ;
875 for (kk = 1; kk <= 2 ; kk++) {
877 PolesPtr->Value(ii,local_index + kk)) ;
879 my2DPoles->SetValue(ii,jj, Point2d) ;
883 for (jj = 1 ; jj <= myNumSubSpaces[1] ; jj++) {
884 error_value = 0.0e0 ;
885 for (ii = 1 ; ii <= NumCurves ; ii++) {
886 local_index = (ii - 1) * TotalNumSS + index ;
887 error_value = Max(ErrorMax(local_index + jj),error_value) ;
890 my2DMaxError->SetValue(jj, error_value) ;
892 for (jj = 1 ; jj <= myNumSubSpaces[1] ; jj++) {
893 error_value = 0.0e0 ;
894 for (ii = 1 ; ii <= NumCurves ; ii++) {
895 local_index = (ii - 1) * TotalNumSS + index ;
896 error_value += AverageError(local_index + jj) ;
899 error_value /= (Standard_Real) NumCurves ;
900 my2DAverageError->SetValue(jj, error_value) ;
902 index += myNumSubSpaces[1] ;
903 // Pour les poles il faut doubler le decalage :
904 dim_index = index + myNumSubSpaces[1];
907 if (myNumSubSpaces[2] > 0) {
909 my3DPoles = new TColgp_HArray2OfPnt(1,PolesPtr->ColLength(),
910 1,myNumSubSpaces[2]) ;
911 my3DMaxError = new TColStd_HArray1OfReal(1,myNumSubSpaces[2]) ;
912 my3DAverageError = new TColStd_HArray1OfReal(1,myNumSubSpaces[2]) ;
913 for (ii = 1 ; ii <= PolesPtr->ColLength() ; ii++) {
914 for (jj = 1 ; jj <= myNumSubSpaces[2] ; jj++) {
915 local_index = dim_index + (jj-1) * 3 ;
916 for (kk = 1; kk <= 3 ; kk++) {
918 PolesPtr->Value(ii,local_index + kk)) ;
920 my3DPoles->SetValue(ii,jj,Point) ;
925 for (jj = 1 ; jj <= myNumSubSpaces[2] ; jj++) {
926 error_value = 0.0e0 ;
927 for (ii = 1 ; ii <= NumCurves ; ii++) {
928 local_index = (ii - 1) * TotalNumSS + index ;
929 error_value = Max(ErrorMax(local_index + jj),error_value) ;
932 my3DMaxError->SetValue(jj, error_value) ;
934 for (jj = 1 ; jj <= myNumSubSpaces[2] ; jj++) {
935 error_value = 0.0e0 ;
936 for (ii = 1 ; ii <= NumCurves ; ii++) {
937 local_index = (ii - 1) * TotalNumSS + index ;
938 error_value += AverageError(local_index + jj) ;
941 error_value /= (Standard_Real) NumCurves ;
942 my3DAverageError->SetValue(jj, error_value) ;
945 if (ErrorCode == 0) {
946 myDone = Standard_True ;
947 myHasResult = Standard_True ;
949 else if (ErrorCode == -1) {
950 myHasResult = Standard_True ;
958 //=======================================================================
961 //=======================================================================
963 void AdvApprox_ApproxAFunction::Poles(const Standard_Integer Index,
964 TColgp_Array1OfPnt& P) const
966 Standard_Integer ii ;
967 for (ii = P.Lower() ; ii <= P.Upper() ; ii++) {
968 P.SetValue(ii,my3DPoles->Value(ii,Index)) ;
973 //=======================================================================
976 //=======================================================================
978 Standard_Integer AdvApprox_ApproxAFunction::NbPoles() const
979 { if (myDone || myHasResult) return BSplCLib::NbPoles(myDegree,
984 //=======================================================================
987 //=======================================================================
989 void AdvApprox_ApproxAFunction::Poles2d(const Standard_Integer Index,
990 TColgp_Array1OfPnt2d& P) const
992 Standard_Integer ii ;
993 for (ii = P.Lower() ; ii <= P.Upper() ; ii++) {
994 P.SetValue(ii,my2DPoles->Value(ii,Index)) ;
997 //=======================================================================
1000 //=======================================================================
1002 void AdvApprox_ApproxAFunction::Poles1d(const Standard_Integer Index,
1003 TColStd_Array1OfReal& P) const
1005 Standard_Integer ii ;
1006 for (ii = P.Lower() ; ii <= P.Upper() ; ii++) {
1007 P.SetValue(ii,my1DPoles->Value(ii,Index)) ;
1010 //=======================================================================
1011 //function : MaxError
1013 //=======================================================================
1014 Handle(TColStd_HArray1OfReal)
1015 AdvApprox_ApproxAFunction::MaxError(const Standard_Integer D) const
1018 Handle(TColStd_HArray1OfReal) EPtr ;
1022 Standard_OutOfRange::Raise() ;
1026 EPtr = my1DMaxError ;
1029 EPtr = my2DMaxError ;
1032 EPtr = my3DMaxError ;
1037 //=======================================================================
1038 //function : MaxError
1040 //=======================================================================
1041 Standard_Real AdvApprox_ApproxAFunction::MaxError(const Standard_Integer D,
1042 const Standard_Integer Index) const
1044 Handle(TColStd_HArray1OfReal) EPtr =
1047 return EPtr->Value(Index) ;
1050 //=======================================================================
1051 //function : AverageError
1053 //=======================================================================
1054 Handle(TColStd_HArray1OfReal)
1055 AdvApprox_ApproxAFunction::AverageError(const Standard_Integer D) const
1058 Handle(TColStd_HArray1OfReal) EPtr ;
1062 Standard_OutOfRange::Raise() ;
1066 EPtr = my1DAverageError ;
1069 EPtr = my2DAverageError ;
1072 EPtr = my3DAverageError ;
1077 //=======================================================================
1078 //function : AverageError
1080 //=======================================================================
1082 Standard_Real AdvApprox_ApproxAFunction::AverageError(
1083 const Standard_Integer D,
1084 const Standard_Integer Index) const
1087 Handle(TColStd_HArray1OfReal) EPtr =
1089 return EPtr->Value(Index) ;
1092 void AdvApprox_ApproxAFunction::Dump(Standard_OStream& o) const
1094 Standard_Integer ii;
1095 o << "Dump of ApproxAFunction" << endl;
1096 if (myNumSubSpaces[0] > 0) {
1097 o << "Error(s) 1d = " << endl;
1098 for (ii=1; ii <= myNumSubSpaces[0]; ii++) {
1099 o << " " << MaxError(1, ii) << endl;
1103 if (myNumSubSpaces[1] > 0) {
1104 o << "Error(s) 2d = " << endl;
1105 for (ii=1; ii <= myNumSubSpaces[1]; ii++) {
1106 o << " " << MaxError(2, ii) << endl;
1110 if (myNumSubSpaces[2] > 0) {
1111 o << "Error(s) 3d = " << endl;
1112 for (ii=1; ii <= myNumSubSpaces[2]; ii++) {
1113 o << " " << MaxError(3, ii) << endl;