1 // Created on: 1995-05-29
2 // Created by: Xavier BENVENISTE
3 // Copyright (c) 1995-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
22 // Modified PMN 22/05/1996 : Pb de recopie des poles
23 // Modified PMN 20/08/1996 : Introduction de la decoupe parametrable.
24 // (plus verif des derives en debug)
25 // Modified PMN 20/08/1996 : Linearisation de F(t) => Les sous Espaces
26 // facilement approchable ont de petites erreurs.
27 // Modified PMN 15/04/1997 : Gestion fine de la continuite aux lieux de decoupes
29 #include <math_Vector.hxx>
30 #include <AdvApprox_ApproxAFunction.ixx>
31 #include <AdvApprox_EvaluatorFunction.hxx>
32 #include <AdvApprox_DichoCutting.hxx>
33 #include <AdvApprox_SimpleApprox.hxx>
35 #include <Precision.hxx>
36 #include <GeomAbs_Shape.hxx>
37 #include <Convert_CompPolynomialToPoles.hxx>
38 #include <BSplCLib.hxx>
40 #include <TColStd_HArray1OfReal.hxx>
41 #include <TColStd_HArray2OfReal.hxx>
42 #include <TColStd_HArray1OfInteger.hxx>
43 #include <TColStd_Array1OfReal.hxx>
44 #include <TColStd_Array1OfInteger.hxx>
47 #include <gp_Vec2d.hxx>
49 #include <TColgp_Array1OfPnt2d.hxx>
50 #include <TColgp_Array1OfPnt.hxx>
51 #include <TColgp_Array1OfPnt2d.hxx>
52 #include <TColgp_Array1OfPnt.hxx>
53 #include <TColgp_HArray2OfPnt2d.hxx>
54 #include <TColgp_HArray2OfPnt.hxx>
57 #include <Standard_ConstructionError.hxx>
58 #include <Standard_OutOfRange.hxx>
61 #include <PLib_JacobiPolynomial.hxx>
66 static Standard_Boolean AdvApprox_Debug = 0;
68 //=====================================================
69 static void MAPDBN(const Standard_Integer dimension,
70 const Standard_Real Debut,
71 const Standard_Real Fin,
72 AdvApprox_EvaluatorFunction& Evaluator,
73 const Standard_Integer Iordre)
74 // Objet : Controle par difference finis, des derives
75 // Warning : En mode Debug, uniquement
76 ///===================================================
78 Standard_Integer derive, OrdreDer, ier, NDIMEN = dimension;
80 Standard_Real h = 1.e-4*(Fin-Debut+1.e-3), eps = 1.e-3, t, ab[2];
81 math_Vector V1(1,NDIMEN), V2(1,NDIMEN), Diff(1,NDIMEN), Der(1,NDIMEN);
83 if (h < 100*RealEpsilon()) {h = 100*RealEpsilon();}
87 for (OrdreDer=1, derive = 0;
88 OrdreDer <= Iordre; OrdreDer++, derive++) {
90 Ptr = (Standard_Real*) &V1.Value(1);
92 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
94 Ptr = (Standard_Real*) &V2.Value(1);
96 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
100 Ptr = (Standard_Real*) &Der.Value(1);
102 Evaluator(&NDIMEN, ab, &t, &OrdreDer, Ptr, &ier);
105 if ( (Diff-Der).Norm() > eps * (Der.Norm()+1) ) {
106 cout << " Debug Ft au parametre t+ = " << t << endl;
107 cout << " Positionement sur la derive "<< OrdreDer
108 << " : " << Der << endl;
109 cout << " Erreur estime : " << (Der-Diff) << endl;
113 Ptr = (Standard_Real*) &V1.Value(1);
115 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
117 Ptr = (Standard_Real*) &V2.Value(1);
119 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
123 Ptr = (Standard_Real*) &Der.Value(1);
125 Evaluator(&NDIMEN, ab, &t, &OrdreDer, Ptr, &ier);
128 if ( (Diff-Der).Norm() > eps * (Der.Norm()+1) ) {
129 cout << " Debug Ft au parametre t- = " << t << endl;
130 cout << " Positionement sur la derive "<< OrdreDer
131 << " : " << Der << endl;
132 cout << " Erreur estime : " << (Der-Diff) << endl;
138 //===================================================================
139 static void PrepareConvert(const Standard_Integer NumCurves,
140 const Standard_Integer MaxDegree,
141 const Standard_Integer ContinuityOrder,
142 const Standard_Integer Num1DSS,
143 const Standard_Integer Num2DSS,
144 const Standard_Integer Num3DSS,
145 const TColStd_Array1OfInteger& NumCoeffPerCurve,
146 TColStd_Array1OfReal& Coefficients,
147 const TColStd_Array2OfReal& PolynomialIntervals,
148 const TColStd_Array1OfReal& TrueIntervals,
149 const TColStd_Array1OfReal& LocalTolerance,
150 TColStd_Array1OfReal& ErrorMax,
151 TColStd_Array1OfInteger& Continuity)
152 // Pour determiner les continuites locales
153 //====================================================================
157 Standard_Boolean isCi;
158 Standard_Integer icurve, idim, iordre, ii,
159 Dimension=Num1DSS + 2*Num2DSS + 3*Num3DSS,
160 NbSpace = Num1DSS + Num2DSS + Num3DSS;
161 Standard_Real diff, moy, facteur1, facteur2, normal1, normal2, eps;
162 Standard_Real *Res1, *Res2, *Val1, *Val2;
163 Standard_Real *Coef1, *Coef2;
164 Standard_Integer RealDegree = Max(MaxDegree + 1, 2 * ContinuityOrder + 2);
169 TColStd_Array1OfReal Result(1, 2*(ContinuityOrder+1)*Dimension);
170 TColStd_Array1OfReal Prec(1, NbSpace), Suivant(1, NbSpace);
173 Res1 = (Standard_Real*) &(Result.ChangeValue(1));
174 Res2 = (Standard_Real*) &(Result.ChangeValue((ContinuityOrder+1)*Dimension+1));
179 if (ContinuityOrder ==0) return;
181 for (icurve=1; icurve<NumCurves; icurve++) {
182 // Init et positionement au noeud
183 isCi = Standard_True;
184 Coef1 = (Standard_Real*) &(Coefficients.Value(
185 (icurve-1) * Dimension * RealDegree + Coefficients.Lower()));
186 Coef2 = Coef1 + Dimension * RealDegree;
187 Standard_Integer Deg1 = NumCoeffPerCurve(NumCoeffPerCurve.Lower()+icurve-1) - 1;
188 PLib::EvalPolynomial(PolynomialIntervals(icurve, 2),
194 Standard_Integer Deg2 = NumCoeffPerCurve(NumCoeffPerCurve.Lower()+icurve) - 1;
195 PLib::EvalPolynomial(PolynomialIntervals(icurve+1, 1),
202 // Verif dans chaque sous espaces.
203 for (iordre=1; iordre<=ContinuityOrder && isCi; iordre++) {
205 // fixing a bug PRO18577
207 Standard_Real Toler = 1.0e-5;
209 Standard_Real f1_dividend = PolynomialIntervals(icurve,2)-PolynomialIntervals(icurve,1);
210 Standard_Real f2_dividend = PolynomialIntervals(icurve+1,2)-PolynomialIntervals(icurve+1,1);
211 Standard_Real f1_divizor = TrueIntervals(icurve+1)-TrueIntervals(icurve);
212 Standard_Real f2_divizor = TrueIntervals(icurve+2)-TrueIntervals(icurve+1);
213 Standard_Real fract1, fract2;
215 if( Abs(f1_divizor) < Toler ) // this is to avoid divizion by zero
216 //in this case fract1 = 5.14755758946803e-85
218 else{ fract1 = f1_dividend/f1_divizor;
219 facteur1 = Pow( fract1, iordre);
221 if( Abs(f2_divizor) < Toler )
222 //in this case fract2 = 6.77193633669143e-313
224 else{ fract2 = f2_dividend/f2_divizor;
225 facteur2 = Pow( fract2, iordre);
227 normal1 = Pow(f1_divizor, iordre);
228 normal2 = Pow(f2_divizor, iordre);
232 Val1 = Res1+iordre*Dimension-1;
233 Val2 = Res2+iordre*Dimension-1;
235 for (ii=1; ii<=Num1DSS && isCi; ii++, idim++) {
236 eps = LocalTolerance(idim)*0.01;
237 diff = Abs (Val1[ii]*facteur1 - Val2[ii]*facteur2);
238 moy = Abs (Val1[ii]*facteur1 + Val2[ii]*facteur2);
239 // Un premier controle sur la valeur relative
240 if (diff > moy*1.e-9) {
241 Prec(idim) = diff*normal1;
242 Suivant(idim) = diff*normal2;
243 // Et un second sur le majorant d'erreur engendree
244 if (Prec(idim)>eps ||
245 Suivant(idim)>eps) isCi=Standard_False;
255 for (ii=1; ii<=Num2DSS && isCi; ii++, idim++, Val1+=2,Val2+=2) {
256 eps = LocalTolerance(idim)*0.01;
257 v1.SetCoord(Val1[1],Val1[2]);
258 v2.SetCoord(Val2[1],Val2[2]);
261 diff = Abs(v1.X() - v2.X()) + Abs(v1.Y() - v2.Y());
262 moy = Abs(v1.X() + v2.X()) + Abs(v1.Y() + v2.Y());
263 if (diff > moy*1.e-9) {
264 Prec(idim) = diff*normal1;
265 Suivant(idim) = diff*normal2;
266 if (Prec(idim)>eps ||
267 Suivant(idim)>eps) isCi=Standard_False;
275 for (ii=1; ii<=Num3DSS && isCi; ii++, idim++, Val1+=3,Val2+=3) {
276 eps = LocalTolerance(idim)*0.01;
277 V1.SetCoord(Val1[1], Val1[2], Val1[3]);
278 V2.SetCoord(Val2[1], Val2[2], Val2[3]);
281 diff = Abs(V1.X() - V2.X()) + Abs(V1.Y() - V2.Y()) +
282 Abs(V1.Z() - V2.Z());
283 moy = Abs(V1.X() + V2.X()) + Abs(V1.Y() + V2.Y()) +
284 Abs(V1.Z() + V2.Z());
285 if (diff > moy*1.e-9) {
286 Prec(idim) = diff*normal1;
287 Suivant(idim) = diff*normal2;
288 if (Prec(idim)>eps ||
289 Suivant(idim)>eps) isCi=Standard_False;
296 // Si c'est bon on update le tout
298 Continuity(icurve+1) = iordre;
299 Standard_Integer index = (icurve-1)*NbSpace+1;
300 for (ii=index, idim=1; idim<=NbSpace; ii++,idim++) {
301 ErrorMax(ii) += Prec(idim);
302 ErrorMax(ii+NbSpace) += Suivant(idim);
310 //=======================================================================
311 //function : ApproxFunction
313 //purpose : Approximation d' UNE fonction non polynomiale (dans l'espace de
314 // dimension NDIMEN) par N courbes polynomiales, par la methode
315 // des moindres carres. Le parametre de la fonction est conserve.
317 // ARGUMENTS D'ENTREE :
318 //C ------------------
319 //C NDIMEN : Dimension totale de l' espace (somme des dimensions
320 //C des sous-espaces).
321 //C NBSESP : Nombre de sous-espaces "independants".
322 //C NDIMSE : Table des dimensions des sous-espaces.
323 //C ABFONC : Bornes de l' intervalle (a,b) de definition de la
324 //C fonction a approcher.
325 //C FONCNP : Fonction externe de positionnement sur la fonction non
326 //C polynomiale a approcher.
327 //C IORDRE : Ordre de contrainte aux extremites de chaque courbe
328 //C polynomiale cree :
329 //C -1 = pas de contraintes,
330 //C 0 = contraintes de passage aux bornes (i.e. C0),
331 //C 1 = C0 + contraintes de derivees 1eres (i.e. C1),
332 //C 2 = C1 + contraintes de derivees 2ndes (i.e. C2).
333 //C NBCRMX : Nombre maxi de courbes polynomiales a calculer.
334 //C NCFLIM : Nombre LIMITE de coeff des "courbes" polynomiales
335 //C d' approximation. Doit etre superieur ou egal a
336 //C 2*IORDRE + 2 et inferieur ou egal a 50 et a NCOFMX).
337 //C Limite a 14 apres l'appel a VRIRAZ pour eviter le bruit.
338 //C EPSAPR : Table des erreurs d' approximations souhaitees
339 //C sous-espace par sous-espace.
340 //C ICODEO : Code d' init. des parametres de discretisation.
341 //C (choix de NBPNT et de NDJAC).
342 //C = 1 calcul rapide avec precision moyenne sur les coeff.
343 //C = 2 calcul rapide avec meilleure precision " " ".
344 //C = 3 calcul un peu plus lent avec bonne precision ".
345 //C = 4 calcul lent avec la meilleure precision possible ".
348 //C ARGUMENTS DE SORTIE :
349 //C -------------------
350 //C NBCRBE : Nombre de courbes polynomiales creees.
351 //C NCFTAB : Table des nombres de coeff. significatifs des NBCRBE
352 //C "courbes" calculees.
353 //C CRBTAB : Tableau des coeff. des "courbes" polynomiales calculees.
354 //C Doit etre dimensionne a CRBTAB(NDIMEN,NCOFMX,NBCRMX).
355 //C TABINT : Table des NBCRBE + 1 bornes des intervalles de decoupe de
357 //C ERRMAX : Table des erreurs (sous-espace par sous espace)
358 //C MAXIMALES commises dans l' approximation de FONCNP par
359 //C par les NBCRBE courbes polynomiales.
360 //C ERRMOY : Table des erreurs (sous-espace par sous espace)
361 //C MOYENNES commises dans l' approximation de FONCNP par
362 //C par les NBCRBE courbes polynomiales.
363 //C IERCOD : Code d' erreur :
364 //C -1 = ERRMAX > EPSAPR pour au moins un des sous-espace.
365 //C On a alors NCRBMX courbes calculees, certaines de
366 //C degre mathematique min(NCFLIM,14)-1 ou la precision
367 //C demandee n' est pas atteinte.
368 //C -2 = Pb dans le mode DEBUG
370 //C 1 = Pb d' incoherence des entrees.
371 //C 10 = Pb de calcul de l' interpolation des contraintes.
372 //C 13 = Pb dans l' allocation dynamique.
373 //C 33 = Pb dans la recuperation des donnees du block data
374 //C des coeff. d' integration par la methode de GAUSS.
375 //C >100 Pb dans l' evaluation de FONCNP, le code d' erreur
376 //C renvoye est egal au code d' erreur de FONCNP + 100.
379 // --> La i-eme courbe polynomiale creee correspond a l' approximation
380 //C de FONCNP sur l' intervalle (TABINT(i),TABINT(i+1)). TABINT(i)
381 // est une suite STRICTEMENT monotone avec TABINT(1)=ABFONC(1) et
382 // TABINT(NBCRBE+1)=ABFONC(2).
384 // --> Si IERCOD=-1, la precision demandee n' est pas atteinte (ERRMAX
385 //C est superieur a EPSAPR sur au moins l' un des sous-espaces), mais
386 // on donne le meilleur resultat possible pour min(NCFLIM,14),NBCRMX
387 // et EPSAPR choisis par l' utilisateur Dans ce cas (ainsi que pour
388 //C IERCOD=0), on a une solution.
390 //C--> ATTENTION : Toute modification du calcul de NDJAC et NBPNT
391 //C doit etre reportee dans le programme MAPNBP.
393 // HISTORIQUE DES MODIFICATIONS :
394 // --------------------------------
396 // 20-08-1996 : PMN ; Creation d'apres la routine Fortran MAPFNC
397 //======================================================================
398 void AdvApprox_ApproxAFunction::Approximation(
399 const Standard_Integer TotalDimension,
400 // Dimension totale de l' espace
401 // (somme des dimensions des sous-espaces).
402 const Standard_Integer TotalNumSS,//Nombre de sous-espaces "independants".
403 const TColStd_Array1OfInteger& LocalDimension,//dimensions des sous-espaces.
404 const Standard_Real First,
405 const Standard_Real Last,
406 // Intervalle (a,b) de definition de la fonction a approcher.
407 AdvApprox_EvaluatorFunction& Evaluator,
408 // Fonction externe de positionnement sur la fonction a approcher.
409 const AdvApprox_Cutting& CutTool,
410 // Maniere de decouper en 2 un intervalle.
411 const Standard_Integer ContinuityOrder,
412 // ordre de continutie a respecter (-1, 0, 1, 2)
413 const Standard_Integer NumMaxCoeffs,
414 // Nombre LIMITE de coeff des "courbes" polynomiales
415 // d' approximation. Doit etre superieur ou egal a
416 // Doit etre superieur ou egal a 2*ContinuityOrder + 2
417 // Limite a 14 pour eviter le bruit.
418 const Standard_Integer MaxSegments,
419 //Nombre maxi de courbes polynomiales a calculer.
420 const TColStd_Array1OfReal& LocalTolerancesArray,
421 //Tolerances d'approximation souhaitees sous-espace par sous-espace.
422 const Standard_Integer code_precis,
423 //Code d' init. des parametres de discretisation.
424 //C (choix de NBPNT et de NDJAC).
425 //C = 1 calcul rapide avec precision moyenne sur les coeff.
426 //C = 2 calcul rapide avec meilleure precision " " ".
427 //C = 3 calcul un peu plus lent avec bonne precision ".
428 //C = 4 calcul lent avec la meilleure precision possible ".
429 Standard_Integer& NumCurves,
430 TColStd_Array1OfInteger& NumCoeffPerCurveArray,
431 TColStd_Array1OfReal& CoefficientArray,
432 TColStd_Array1OfReal& IntervalsArray,
433 TColStd_Array1OfReal& ErrorMaxArray,
434 TColStd_Array1OfReal& AverageErrorArray,
435 Standard_Integer& ErrorCode)
437 // Standard_Real EpsPar = Precision::Confusion();
438 Standard_Integer IDIM, NUPIL,TheDeg;
440 Standard_Integer NDIMEN = TotalDimension;
442 Standard_Boolean isCut = Standard_False;
444 // Defintion des Tableaux C
445 Standard_Real * TABINT = (Standard_Real*) &IntervalsArray(1);
449 CoefficientArray.Init(0);
451 //-------------------------- Verification des entrees ------------------
453 if ((MaxSegments<1)|| (Abs(Last-First)<1.e-9))
457 //--> La dimension totale doit etre la somme des dimensions des
460 for ( Standard_Integer I=1; I<=TotalNumSS; I++) {IDIM += LocalDimension(I);}
461 if (IDIM != TotalDimension)
464 GeomAbs_Shape Continuity=GeomAbs_C0;
465 switch (ContinuityOrder) {
467 Continuity = GeomAbs_C0 ;
470 Continuity = GeomAbs_C1 ;
473 Continuity = GeomAbs_C2 ;
476 Standard_ConstructionError::Raise();
479 //--------------------- Choix du nombre de points ----------------------
480 //---------- et du degre de lissage dans la base orthogonale -----------
481 //--> NDJAC est le degre de "travail" dans la base orthogonale.
484 Standard_Integer NbGaussPoints, WorkDegree;
486 PLib::JacobiParameters(Continuity, NumMaxCoeffs-1, code_precis,
487 NbGaussPoints, WorkDegree);
490 //------------------ Initialisation de la gestion des decoupes ---------
497 //C**********************************************************************
498 //C APPROXIMATION AVEC DECOUPES
499 //C**********************************************************************
500 Handle(PLib_JacobiPolynomial) JacobiBase = new (PLib_JacobiPolynomial) (WorkDegree, Continuity);
501 //Portage HP le compilateur refuse le debranchement
502 Standard_Integer IS ;
503 Standard_Boolean goto_fin_de_boucle;
504 Standard_Integer MaxDegree = NumMaxCoeffs-1;
505 AdvApprox_SimpleApprox Approx (TotalDimension, TotalNumSS,
507 WorkDegree, NbGaussPoints,
508 JacobiBase, Evaluator);
509 while ((NUPIL-NumCurves)!=0) {
510 //--> Lorsque l' on a atteint le haut de la pile, c' est fini !
512 //Portage HP le compilateur refuse le debranchement
513 goto_fin_de_boucle=Standard_False;
515 //---- Calcul de la courbe d' approximation dans la base de Jacobi -----
517 Approx.Perform ( LocalDimension, LocalTolerancesArray,
518 TABINT[NumCurves], TABINT[NumCurves+1],
520 if (!Approx.IsDone()) {
527 //---------- Calcul du degre de la courbe et de l' erreur max ----------
530 NumCoeffPerCurveArray(NumCurves + 1)=0;
532 // L'erreur doit etre satisfaite sur tous les sous-espaces sinon, on decoupe
534 Standard_Boolean MaxMaxErr=Standard_True;
535 for ( IS=1; IS<=TotalNumSS; IS++)
536 { if (Approx.MaxError(IS) > LocalTolerancesArray(IS))
537 { MaxMaxErr = Standard_False;
542 if (MaxMaxErr == Standard_True)
545 // NumCoeffPerCurveArray(NumCurves) = Approx.Degree()+1;
549 //-> ...sinon on essai de decouper l' intervalle courant en 2...
551 Standard_Boolean Large;
553 Large = CutTool.Value(TABINT[NumCurves], TABINT[NumCurves+1],
556 if (NUPIL<MaxSegments && Large) {
558 // if (IS!=1) {NumCurves--;}
559 isCut = Standard_True; // Ca y est, on le sait !
560 Standard_Real * IDEB = TABINT+NumCurves+1;
561 Standard_Real * IDEB1 = IDEB+1;
562 Standard_Integer ILONG= NUPIL-NumCurves-1;
563 for (Standard_Integer iI=ILONG; iI>=0; iI--) {
564 IDEB1[iI] = IDEB[iI];
568 //--> ... et on recommence.
569 //Portage HP le compilateur refuse le debranchement
570 goto_fin_de_boucle=Standard_True;
574 //--> Si la pile est pleine...
575 // pour l'instant ErrorCode=-1;
577 // NumCoeffPerCurveArray(NumCurves) = Approx.Degree()+1;
581 //Portage HP le compilateur refuse le debranchement
582 if (goto_fin_de_boucle) continue;
583 for (IS=1; IS<=TotalNumSS; IS++) {
584 ErrorMaxArray.SetValue(IS+(NumCurves-1)*TotalNumSS,Approx. MaxError(IS));
585 //---> ...et calcul de l' erreur moyenne.
586 AverageErrorArray.SetValue(IS+(NumCurves-1)*TotalNumSS,Approx.AverageError(IS));
589 Handle(TColStd_HArray1OfReal) HJacCoeff = Approx.Coefficients();
590 TheDeg = Approx.Degree();
591 if (isCut && (TheDeg < 2*ContinuityOrder+1) )
592 // Pour ne pas bruiter les derives aux bouts, et garder une continuite
593 // correcte sur la BSpline resultat.
594 TheDeg = 2*ContinuityOrder+1;
596 NumCoeffPerCurveArray(NumCurves) = TheDeg+1;
597 TColStd_Array1OfReal Coefficients(0,(TheDeg+1)*TotalDimension-1);
598 JacobiBase->ToCoefficients (TotalDimension, TheDeg,
599 HJacCoeff->Array1(), Coefficients);
600 Standard_Integer i,j, f = (TheDeg+1)*TotalDimension;
601 for (i=0,j=(NumCurves-1)*TotalDimension*NumMaxCoeffs+1;
603 CoefficientArray.SetValue(j, Coefficients.Value(i));
607 // Test des derives par difference finis
608 Standard_Integer IORDRE = ContinuityOrder;
610 if (IORDRE>0 && AdvApprox_Debug) {
611 MAPDBN(NDIMEN, TABINT[NumCurves-1],
612 TABINT[NumCurves], Evaluator, IORDRE);
615 //Portage HP le compilateur refuse le debranchement
617 // {;} fin de la boucle while
621 //=======================================================================
622 //function : AdvApprox_ApproxAFunction
623 //purpose : Constructeur avec Decoupe Dichotomique
624 //=======================================================================
626 AdvApprox_ApproxAFunction::
627 AdvApprox_ApproxAFunction(const Standard_Integer Num1DSS,
628 const Standard_Integer Num2DSS,
629 const Standard_Integer Num3DSS,
630 const Handle_TColStd_HArray1OfReal& OneDTol,
631 const Handle_TColStd_HArray1OfReal& TwoDTol,
632 const Handle_TColStd_HArray1OfReal& ThreeDTol,
633 const Standard_Real First,
634 const Standard_Real Last,
635 const GeomAbs_Shape Continuity,
636 const Standard_Integer MaxDeg,
637 const Standard_Integer MaxSeg,
638 const AdvApprox_EvaluatorFunction& Func) :
639 my1DTolerances(OneDTol),
640 my2DTolerances(TwoDTol),
641 my3DTolerances(ThreeDTol),
644 myContinuity(Continuity),
646 myMaxSegments(MaxSeg),
647 myDone(Standard_False),
648 myHasResult(Standard_False),
649 myEvaluator((Standard_Address)&Func)
651 AdvApprox_DichoCutting Cut;
652 Perform(Num1DSS, Num2DSS, Num3DSS, Cut);
655 AdvApprox_ApproxAFunction::
656 AdvApprox_ApproxAFunction(const Standard_Integer Num1DSS,
657 const Standard_Integer Num2DSS,
658 const Standard_Integer Num3DSS,
659 const Handle_TColStd_HArray1OfReal& OneDTol,
660 const Handle_TColStd_HArray1OfReal& TwoDTol,
661 const Handle_TColStd_HArray1OfReal& ThreeDTol,
662 const Standard_Real First,
663 const Standard_Real Last,
664 const GeomAbs_Shape Continuity,
665 const Standard_Integer MaxDeg,
666 const Standard_Integer MaxSeg,
667 const AdvApprox_EvaluatorFunction& Func,
668 const AdvApprox_Cutting& CutTool) :
669 my1DTolerances(OneDTol),
670 my2DTolerances(TwoDTol),
671 my3DTolerances(ThreeDTol),
674 myContinuity(Continuity),
676 myMaxSegments(MaxSeg),
677 myDone(Standard_False),
678 myHasResult(Standard_False),
679 myEvaluator((Standard_Address)&Func)
681 Perform(Num1DSS, Num2DSS, Num3DSS, CutTool);
684 //=======================================================================
685 //function : AdvApprox_ApproxAFunction::Perform
686 //purpose : Make all the Work !!
687 //=======================================================================
688 void AdvApprox_ApproxAFunction::Perform(const Standard_Integer Num1DSS,
689 const Standard_Integer Num2DSS,
690 const Standard_Integer Num3DSS,
691 const AdvApprox_Cutting& CutTool)
696 Num1DSS + Num2DSS + Num3DSS <= 0 ||
700 Standard_ConstructionError::Raise();
701 if (myMaxDegree > 14) {
707 myNumSubSpaces[0] = Num1DSS ;
708 myNumSubSpaces[1] = Num2DSS ;
709 myNumSubSpaces[2] = Num3DSS ;
710 Standard_Integer TotalNumSS =
711 Num1DSS + Num2DSS + Num3DSS,
718 Standard_Integer TotalDimension =
719 myNumSubSpaces[0] + 2 * myNumSubSpaces[1] + 3 * myNumSubSpaces[2] ;
720 Standard_Real error_value ;
722 Standard_Integer ContinuityOrder=0 ;
723 switch (myContinuity) {
725 ContinuityOrder = 0 ;
728 ContinuityOrder = 1 ;
731 ContinuityOrder = 2 ;
734 Standard_ConstructionError::Raise();
736 Standard_Real ApproxStartEnd[2] ;
737 Standard_Integer NumMaxCoeffs = Max(myMaxDegree + 1, 2 * ContinuityOrder + 2);
738 myMaxDegree = NumMaxCoeffs - 1;
739 Standard_Integer code_precis = 1 ;
741 // WARNING : the polynomial coefficients are the
742 // taylor expansion of the polynomial at 0.0e0 !
744 ApproxStartEnd[0] = -1.0e0 ;
745 ApproxStartEnd[1] = 1.0e0 ;
747 TColStd_Array1OfInteger LocalDimension(1,TotalNumSS) ;
751 TColStd_Array1OfReal LocalTolerances(1,TotalNumSS) ;
753 for(jj = 1; jj <= myNumSubSpaces[0] ; jj++) {
754 LocalTolerances.SetValue(index,my1DTolerances->Value(jj)) ;
755 LocalDimension.SetValue(index,1) ;
758 for(jj = 1 ; jj <= myNumSubSpaces[1] ; jj++) {
759 LocalTolerances.SetValue(index,my2DTolerances->Value(jj)) ;
760 LocalDimension.SetValue(index,2) ;
763 for(jj = 1; jj <= myNumSubSpaces[2] ; jj++) {
764 LocalTolerances.SetValue(index,my3DTolerances->Value(jj)) ;
765 LocalDimension.SetValue(index,3) ;
772 Standard_Integer ErrorCode = 0,
775 myMaxSegments * NumMaxCoeffs * TotalDimension ;
776 Handle(TColStd_HArray1OfInteger) NumCoeffPerCurvePtr =
777 new TColStd_HArray1OfInteger (1,myMaxSegments) ;
779 Handle(TColStd_HArray1OfReal) LocalCoefficientsPtr =
780 new TColStd_HArray1OfReal(1,size) ;
782 Handle(TColStd_HArray1OfReal) IntervalsPtr =
783 new TColStd_HArray1OfReal (1,myMaxSegments+1) ;
785 TColStd_Array1OfReal ErrorMax(1,myMaxSegments * TotalNumSS) ;
787 TColStd_Array1OfReal AverageError(1,myMaxSegments * TotalNumSS) ;
789 Approximation ( TotalDimension,
790 TotalNumSS, LocalDimension,
792 *(AdvApprox_EvaluatorFunction*)myEvaluator,
793 CutTool, ContinuityOrder, NumMaxCoeffs,
794 myMaxSegments, LocalTolerances, code_precis,
795 NumCurves, // Nombre de courbe en sortie
796 NumCoeffPerCurvePtr->ChangeArray1(), // Nbre de coeff par courbe
797 LocalCoefficientsPtr->ChangeArray1(), // Les Coeffs solutions
798 IntervalsPtr->ChangeArray1(), // La Table des decoupes
799 ErrorMax, // Majoration de l'erreur
800 AverageError, // Erreur moyenne constatee
803 if (ErrorCode == 0 || ErrorCode == -1) {
805 // si tout est OK ou bien on a un resultat dont l une des erreurs max est
806 // plus grande que la tolerance demandee
808 TColStd_Array1OfInteger TabContinuity(1, NumCurves);
809 TColStd_Array2OfReal PolynomialIntervalsPtr (1,NumCurves,1,2);
810 for (ii = PolynomialIntervalsPtr.LowerRow() ;
811 ii <= PolynomialIntervalsPtr.UpperRow() ;
813 // On force un degre 1 minimum (PRO5474)
814 NumCoeffPerCurvePtr->ChangeValue(ii) =
815 Max(2, NumCoeffPerCurvePtr->Value(ii));
816 PolynomialIntervalsPtr.SetValue(ii,1,ApproxStartEnd[0]) ;
817 PolynomialIntervalsPtr.SetValue(ii,2,ApproxStartEnd[1]) ;
820 PrepareConvert(NumCurves, myMaxDegree, ContinuityOrder,
821 Num1DSS, Num2DSS, Num3DSS,
822 NumCoeffPerCurvePtr->Array1(),
823 LocalCoefficientsPtr->ChangeArray1(),
824 PolynomialIntervalsPtr, IntervalsPtr->Array1(),
825 LocalTolerances, ErrorMax,
828 Convert_CompPolynomialToPoles
829 AConverter(NumCurves,
833 NumCoeffPerCurvePtr->Array1(),
834 LocalCoefficientsPtr->Array1(),
835 PolynomialIntervalsPtr,
836 IntervalsPtr->Array1());
838 if (AConverter.IsDone()) {
839 Handle(TColStd_HArray2OfReal) PolesPtr ;
840 AConverter.Poles(PolesPtr) ;
841 AConverter.Knots(myKnots) ;
842 AConverter.Multiplicities(myMults) ;
843 myDegree = AConverter.Degree() ;
845 if (myNumSubSpaces[0] > 0) {
846 my1DPoles = new TColStd_HArray2OfReal(1,PolesPtr->ColLength(),
847 1,myNumSubSpaces[0]) ;
848 my1DMaxError = new TColStd_HArray1OfReal(1,myNumSubSpaces[0]) ;
849 my1DAverageError = new TColStd_HArray1OfReal(1,myNumSubSpaces[0]) ;
850 for (ii = 1 ; ii <= PolesPtr->ColLength() ; ii++) {
851 for (jj = 1 ; jj <= myNumSubSpaces[0] ; jj++) {
852 my1DPoles->SetValue(ii,jj, PolesPtr->Value(ii,jj)) ;
856 for (jj = 1 ; jj <= myNumSubSpaces[0] ; jj++) {
857 error_value = 0.0e0 ;
858 for (ii = 1 ; ii <= NumCurves ; ii++) {
859 local_index = (ii - 1) * TotalNumSS ;
860 error_value = Max(ErrorMax(local_index + jj),error_value) ;
863 my1DMaxError->SetValue(jj, error_value) ;
865 for (jj = 1 ; jj <= myNumSubSpaces[0] ; jj++) {
866 error_value = 0.0e0 ;
867 for (ii = 1 ; ii <= NumCurves ; ii++) {
868 local_index = (ii - 1) * TotalNumSS ;
869 error_value += AverageError(local_index + jj) ;
872 error_value /= (Standard_Real) NumCurves ;
873 my1DAverageError->SetValue(jj, error_value) ;
875 index += myNumSubSpaces[0] ;
878 dim_index = index; //Pour le cas ou il n'y a pas de 2D
880 if (myNumSubSpaces[1] > 0) {
882 my2DPoles = new TColgp_HArray2OfPnt2d(1,PolesPtr->ColLength(),
883 1,myNumSubSpaces[1]) ;
884 my2DMaxError = new TColStd_HArray1OfReal(1,myNumSubSpaces[1]) ;
885 my2DAverageError = new TColStd_HArray1OfReal(1,myNumSubSpaces[1]) ;
886 for (ii = 1 ; ii <= PolesPtr->ColLength() ; ii++) {
887 for (jj = 1 ; jj <= myNumSubSpaces[1] ; jj++) {
888 local_index = index + (jj-1) * 2 ;
889 for (kk = 1; kk <= 2 ; kk++) {
891 PolesPtr->Value(ii,local_index + kk)) ;
893 my2DPoles->SetValue(ii,jj, Point2d) ;
897 for (jj = 1 ; jj <= myNumSubSpaces[1] ; jj++) {
898 error_value = 0.0e0 ;
899 for (ii = 1 ; ii <= NumCurves ; ii++) {
900 local_index = (ii - 1) * TotalNumSS + index ;
901 error_value = Max(ErrorMax(local_index + jj),error_value) ;
904 my2DMaxError->SetValue(jj, error_value) ;
906 for (jj = 1 ; jj <= myNumSubSpaces[1] ; jj++) {
907 error_value = 0.0e0 ;
908 for (ii = 1 ; ii <= NumCurves ; ii++) {
909 local_index = (ii - 1) * TotalNumSS + index ;
910 error_value += AverageError(local_index + jj) ;
913 error_value /= (Standard_Real) NumCurves ;
914 my2DAverageError->SetValue(jj, error_value) ;
916 index += myNumSubSpaces[1] ;
917 // Pour les poles il faut doubler le decalage :
918 dim_index = index + myNumSubSpaces[1];
921 if (myNumSubSpaces[2] > 0) {
923 my3DPoles = new TColgp_HArray2OfPnt(1,PolesPtr->ColLength(),
924 1,myNumSubSpaces[2]) ;
925 my3DMaxError = new TColStd_HArray1OfReal(1,myNumSubSpaces[2]) ;
926 my3DAverageError = new TColStd_HArray1OfReal(1,myNumSubSpaces[2]) ;
927 for (ii = 1 ; ii <= PolesPtr->ColLength() ; ii++) {
928 for (jj = 1 ; jj <= myNumSubSpaces[2] ; jj++) {
929 local_index = dim_index + (jj-1) * 3 ;
930 for (kk = 1; kk <= 3 ; kk++) {
932 PolesPtr->Value(ii,local_index + kk)) ;
934 my3DPoles->SetValue(ii,jj,Point) ;
939 for (jj = 1 ; jj <= myNumSubSpaces[2] ; jj++) {
940 error_value = 0.0e0 ;
941 for (ii = 1 ; ii <= NumCurves ; ii++) {
942 local_index = (ii - 1) * TotalNumSS + index ;
943 error_value = Max(ErrorMax(local_index + jj),error_value) ;
946 my3DMaxError->SetValue(jj, error_value) ;
948 for (jj = 1 ; jj <= myNumSubSpaces[2] ; jj++) {
949 error_value = 0.0e0 ;
950 for (ii = 1 ; ii <= NumCurves ; ii++) {
951 local_index = (ii - 1) * TotalNumSS + index ;
952 error_value += AverageError(local_index + jj) ;
955 error_value /= (Standard_Real) NumCurves ;
956 my3DAverageError->SetValue(jj, error_value) ;
959 if (ErrorCode == 0) {
960 myDone = Standard_True ;
961 myHasResult = Standard_True ;
963 else if (ErrorCode == -1) {
964 myHasResult = Standard_True ;
972 //=======================================================================
975 //=======================================================================
977 void AdvApprox_ApproxAFunction::Poles(const Standard_Integer Index,
978 TColgp_Array1OfPnt& P) const
980 Standard_Integer ii ;
981 for (ii = P.Lower() ; ii <= P.Upper() ; ii++) {
982 P.SetValue(ii,my3DPoles->Value(ii,Index)) ;
987 //=======================================================================
990 //=======================================================================
992 Standard_Integer AdvApprox_ApproxAFunction::NbPoles() const
993 { if (myDone || myHasResult) return BSplCLib::NbPoles(myDegree,
998 //=======================================================================
1001 //=======================================================================
1003 void AdvApprox_ApproxAFunction::Poles2d(const Standard_Integer Index,
1004 TColgp_Array1OfPnt2d& P) const
1006 Standard_Integer ii ;
1007 for (ii = P.Lower() ; ii <= P.Upper() ; ii++) {
1008 P.SetValue(ii,my2DPoles->Value(ii,Index)) ;
1011 //=======================================================================
1014 //=======================================================================
1016 void AdvApprox_ApproxAFunction::Poles1d(const Standard_Integer Index,
1017 TColStd_Array1OfReal& P) const
1019 Standard_Integer ii ;
1020 for (ii = P.Lower() ; ii <= P.Upper() ; ii++) {
1021 P.SetValue(ii,my1DPoles->Value(ii,Index)) ;
1024 //=======================================================================
1025 //function : MaxError
1027 //=======================================================================
1028 Handle(TColStd_HArray1OfReal)
1029 AdvApprox_ApproxAFunction::MaxError(const Standard_Integer D) const
1032 Handle(TColStd_HArray1OfReal) EPtr ;
1036 Standard_OutOfRange::Raise() ;
1040 EPtr = my1DMaxError ;
1043 EPtr = my2DMaxError ;
1046 EPtr = my3DMaxError ;
1051 //=======================================================================
1052 //function : MaxError
1054 //=======================================================================
1055 Standard_Real AdvApprox_ApproxAFunction::MaxError(const Standard_Integer D,
1056 const Standard_Integer Index) const
1058 Handle(TColStd_HArray1OfReal) EPtr =
1061 return EPtr->Value(Index) ;
1064 //=======================================================================
1065 //function : AverageError
1067 //=======================================================================
1068 Handle(TColStd_HArray1OfReal)
1069 AdvApprox_ApproxAFunction::AverageError(const Standard_Integer D) const
1072 Handle(TColStd_HArray1OfReal) EPtr ;
1076 Standard_OutOfRange::Raise() ;
1080 EPtr = my1DAverageError ;
1083 EPtr = my2DAverageError ;
1086 EPtr = my3DAverageError ;
1091 //=======================================================================
1092 //function : AverageError
1094 //=======================================================================
1096 Standard_Real AdvApprox_ApproxAFunction::AverageError(
1097 const Standard_Integer D,
1098 const Standard_Integer Index) const
1101 Handle(TColStd_HArray1OfReal) EPtr =
1103 return EPtr->Value(Index) ;
1106 void AdvApprox_ApproxAFunction::Dump(Standard_OStream& o) const
1108 Standard_Integer ii;
1109 o << "Dump of ApproxAFunction" << endl;
1110 if (myNumSubSpaces[0] > 0) {
1111 o << "Error(s) 1d = " << endl;
1112 for (ii=1; ii <= myNumSubSpaces[0]; ii++) {
1113 o << " " << MaxError(1, ii) << endl;
1117 if (myNumSubSpaces[1] > 0) {
1118 o << "Error(s) 2d = " << endl;
1119 for (ii=1; ii <= myNumSubSpaces[1]; ii++) {
1120 o << " " << MaxError(2, ii) << endl;
1124 if (myNumSubSpaces[2] > 0) {
1125 o << "Error(s) 3d = " << endl;
1126 for (ii=1; ii <= myNumSubSpaces[2]; ii++) {
1127 o << " " << MaxError(3, ii) << endl;