projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0031668: Visualization - WebGL sample doesn't work on Emscripten 1.39
[occt.git] / src / AdvApprox / AdvApprox_ApproxAFunction.cxx
CommitLineData
b311480e 1// Created on: 1995-05-29
2// Created by: Xavier BENVENISTE
3// Copyright (c) 1995-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 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
973c2be1 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.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
7fd59977 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
7fd59977 23
42cf5bc1 24#include <AdvApprox_ApproxAFunction.hxx>
25#include <AdvApprox_Cutting.hxx>
7fd59977 26#include <AdvApprox_DichoCutting.hxx>
42cf5bc1 27#include <AdvApprox_EvaluatorFunction.hxx>
7fd59977 28#include <AdvApprox_SimpleApprox.hxx>
7fd59977 29#include <BSplCLib.hxx>
42cf5bc1 30#include <Convert_CompPolynomialToPoles.hxx>
31#include <GeomAbs_Shape.hxx>
7fd59977 32#include <gp_Vec.hxx>
33#include <gp_Vec2d.hxx>
42cf5bc1 34#include <math_Vector.hxx>
35#include <PLib.hxx>
36#include <PLib_JacobiPolynomial.hxx>
37#include <Precision.hxx>
38#include <Standard_ConstructionError.hxx>
39#include <Standard_OutOfRange.hxx>
7fd59977 40#include <TColgp_Array1OfPnt.hxx>
41#include <TColgp_Array1OfPnt2d.hxx>
7fd59977 42#include <TColgp_HArray2OfPnt.hxx>
42cf5bc1 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>
7fd59977 49
0797d9d3 50#ifdef OCCT_DEBUG
7fd59977 51
52static Standard_Boolean AdvApprox_Debug = 0;
53
54//=====================================================
55static 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///===================================================
63{
64 Standard_Integer derive, OrdreDer, ier, NDIMEN = dimension;
65 Standard_Real* Ptr;
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);
68
69 if (h < 100*RealEpsilon()) {h = 100*RealEpsilon();}
70 ab[0] = Debut;
71 ab[1] = Fin;
72
73 for (OrdreDer=1, derive = 0;
74 OrdreDer <= Iordre; OrdreDer++, derive++) {
75 // Verif au debut
76 Ptr = (Standard_Real*) &V1.Value(1);
77 t = Debut+h;
78 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
79
80 Ptr = (Standard_Real*) &V2.Value(1);
81 t = Debut;
82 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
83
84 Diff = (V1 - V2)/h;
85
86 Ptr = (Standard_Real*) &Der.Value(1);
87 t = Debut;
88 Evaluator(&NDIMEN, ab, &t, &OrdreDer, Ptr, &ier);
89
90
91 if ( (Diff-Der).Norm() > eps * (Der.Norm()+1) ) {
04232180 92 std::cout << " Debug Ft au parametre t+ = " << t << std::endl;
93 std::cout << " Positionement sur la derive "<< OrdreDer
94 << " : " << Der << std::endl;
95 std::cout << " Erreur estime : " << (Der-Diff) << std::endl;
7fd59977 96 }
97
98 // Verif a la fin
99 Ptr = (Standard_Real*) &V1.Value(1);
100 t = Fin-h;
101 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
102
103 Ptr = (Standard_Real*) &V2.Value(1);
104 t = Fin;
105 Evaluator(&NDIMEN, ab, &t, &derive, Ptr, &ier);
106
107 Diff = (V2 - V1)/h;
108
109 Ptr = (Standard_Real*) &Der.Value(1);
110 t = Fin;
111 Evaluator(&NDIMEN, ab, &t, &OrdreDer, Ptr, &ier);
112
113
114 if ( (Diff-Der).Norm() > eps * (Der.Norm()+1) ) {
04232180 115 std::cout << " Debug Ft au parametre t- = " << t << std::endl;
116 std::cout << " Positionement sur la derive "<< OrdreDer
117 << " : " << Der << std::endl;
118 std::cout << " Erreur estime : " << (Der-Diff) << std::endl;
7fd59977 119 }
120 }
121}
122#endif
123
124//===================================================================
125static 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//====================================================================
140
141{
142 // Declaration
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);
151
152 gp_Vec V1,V2;
153 gp_Vec2d v1, v2;
154
155 TColStd_Array1OfReal Result(1, 2*(ContinuityOrder+1)*Dimension);
156 TColStd_Array1OfReal Prec(1, NbSpace), Suivant(1, NbSpace);
157
158
159 Res1 = (Standard_Real*) &(Result.ChangeValue(1));
160 Res2 = (Standard_Real*) &(Result.ChangeValue((ContinuityOrder+1)*Dimension+1));
161
162
163 //Init
164 Continuity.Init(0);
165 if (ContinuityOrder ==0) return;
166
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),
175 ContinuityOrder,
176 Deg1,
177 Dimension,
178 Coef1[0],
179 Res1[0]);
180 Standard_Integer Deg2 = NumCoeffPerCurve(NumCoeffPerCurve.Lower()+icurve) - 1;
181 PLib::EvalPolynomial(PolynomialIntervals(icurve+1, 1),
182 ContinuityOrder,
183 Deg2,
184 Dimension,
185 Coef2[0],
186 Res2[0]);
187
188 // Verif dans chaque sous espaces.
189 for (iordre=1; iordre<=ContinuityOrder && isCi; iordre++) {
190
191 // fixing a bug PRO18577
192
193 Standard_Real Toler = 1.0e-5;
194
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;
200
201 if( Abs(f1_divizor) < Toler ) // this is to avoid divizion by zero
202 //in this case fract1 = 5.14755758946803e-85
203 facteur1 = 0.0;
204 else{ fract1 = f1_dividend/f1_divizor;
205 facteur1 = Pow( fract1, iordre);
206 }
207 if( Abs(f2_divizor) < Toler )
208 //in this case fract2 = 6.77193633669143e-313
209 facteur2 = 0.0;
210 else{ fract2 = f2_dividend/f2_divizor;
211 facteur2 = Pow( fract2, iordre);
212 }
213 normal1 = Pow(f1_divizor, iordre);
214 normal2 = Pow(f2_divizor, iordre);
215
216
217 idim = 1;
218 Val1 = Res1+iordre*Dimension-1;
219 Val2 = Res2+iordre*Dimension-1;
220
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;
232 }
233 else {
234 Prec(idim) = 0;
235 Suivant(idim) = 0;
236 }
237 }
238
239 Val1+=Num1DSS;
240 Val2+=Num1DSS;
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]);
245 v1 *= facteur1;
246 v2 *= facteur2;
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;
254 }
255 else {
256 Prec(idim) = 0;
257 Suivant(idim) = 0;
258 }
259 }
260
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]);
265 V1 *= facteur1;
266 V2 *= facteur2;
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;
276 }
277 else {
278 Prec(idim) = 0;
279 Suivant(idim) = 0;
280 }
281 }
282 // Si c'est bon on update le tout
283 if (isCi) {
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);
289 }
290 }
291 }
292 }
293}
294
295
296//=======================================================================
297//function : ApproxFunction
298//
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.
302//
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 ".
332//C
333//C
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
342//C FONCNP.
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
355//C 0 = Tout est ok.
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.
363//C
364//
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).
369//
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.
375//C
376//C--> ATTENTION : Toute modification du calcul de NDJAC et NBPNT
377//C doit etre reportee dans le programme MAPNBP.
378//
379// HISTORIQUE DES MODIFICATIONS :
380// --------------------------------
381//
382// 20-08-1996 : PMN ; Creation d'apres la routine Fortran MAPFNC
383//======================================================================
384void 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)
422{
423// Standard_Real EpsPar = Precision::Confusion();
424 Standard_Integer IDIM, NUPIL,TheDeg;
0797d9d3 425#ifdef OCCT_DEBUG
7fd59977 426 Standard_Integer NDIMEN = TotalDimension;
427#endif
428 Standard_Boolean isCut = Standard_False;
429
430// Defintion des Tableaux C
431 Standard_Real * TABINT = (Standard_Real*) &IntervalsArray(1);
432
433
434 ErrorCode=0;
435 CoefficientArray.Init(0);
436
437//-------------------------- Verification des entrees ------------------
438
439 if ((MaxSegments<1)|| (Abs(Last-First)<1.e-9))
440 {ErrorCode=1;
441 return;}
442
443//--> La dimension totale doit etre la somme des dimensions des
444// sous-espaces
445 IDIM=0;
446 for ( Standard_Integer I=1; I<=TotalNumSS; I++) {IDIM += LocalDimension(I);}
447 if (IDIM != TotalDimension)
448 {ErrorCode=1;
449 return;}
450 GeomAbs_Shape Continuity=GeomAbs_C0;
451 switch (ContinuityOrder) {
452 case 0:
453 Continuity = GeomAbs_C0 ;
454 break ;
455 case 1:
456 Continuity = GeomAbs_C1 ;
457 break ;
458 case 2:
459 Continuity = GeomAbs_C2 ;
460 break ;
461 default:
9775fa61 462 throw Standard_ConstructionError();
7fd59977 463 }
464
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.
468
469
470 Standard_Integer NbGaussPoints, WorkDegree;
471
472 PLib::JacobiParameters(Continuity, NumMaxCoeffs-1, code_precis,
473 NbGaussPoints, WorkDegree);
474// NDJAC=WorkDegree;
475
476//------------------ Initialisation de la gestion des decoupes ---------
477
478 TABINT[0]=First;
479 TABINT[1]=Last;
480 NUPIL=1;
481 NumCurves=0;
482
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,
492 Continuity,
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 !
497
498//Portage HP le compilateur refuse le debranchement
499 goto_fin_de_boucle=Standard_False;
500
501//---- Calcul de la courbe d' approximation dans la base de Jacobi -----
502
503 Approx.Perform ( LocalDimension, LocalTolerancesArray,
504 TABINT[NumCurves], TABINT[NumCurves+1],
505 MaxDegree);
506 if (!Approx.IsDone()) {
507 ErrorCode = 1;
508 return;
509 }
510
511
512
513//---------- Calcul du degre de la courbe et de l' erreur max ----------
514
515 IDIM=0;
516 NumCoeffPerCurveArray(NumCurves + 1)=0;
517
518// L'erreur doit etre satisfaite sur tous les sous-espaces sinon, on decoupe
519
520 Standard_Boolean MaxMaxErr=Standard_True;
521 for ( IS=1; IS<=TotalNumSS; IS++)
522 { if (Approx.MaxError(IS) > LocalTolerancesArray(IS))
523 { MaxMaxErr = Standard_False;
524 break;
525 }
526 }
527
528 if (MaxMaxErr == Standard_True)
529 {
530 NumCurves++;
531// NumCoeffPerCurveArray(NumCurves) = Approx.Degree()+1;
532 }
533 else
534 {
535//-> ...sinon on essai de decouper l' intervalle courant en 2...
536 Standard_Real TMIL;
537 Standard_Boolean Large;
538
539 Large = CutTool.Value(TABINT[NumCurves], TABINT[NumCurves+1],
540 TMIL);
541
542 if (NUPIL<MaxSegments && Large) {
543
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];
551 }
552 IDEB[0] = TMIL;
553 NUPIL++;
554//--> ... et on recommence.
555//Portage HP le compilateur refuse le debranchement
556 goto_fin_de_boucle=Standard_True;
557// break;
558 }
559 else {
560//--> Si la pile est pleine...
561// pour l'instant ErrorCode=-1;
562 NumCurves++;
563// NumCoeffPerCurveArray(NumCurves) = Approx.Degree()+1;
564 }
565 }
566// IDIM += NDSES;
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));
573 }
574
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;
581
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;
588 i<f; i++, j++) {
589 CoefficientArray.SetValue(j, Coefficients.Value(i));
590 }
591
0797d9d3 592#ifdef OCCT_DEBUG
7fd59977 593 // Test des derives par difference finis
594 Standard_Integer IORDRE = ContinuityOrder;
595
596 if (IORDRE>0 && AdvApprox_Debug) {
597 MAPDBN(NDIMEN, TABINT[NumCurves-1],
598 TABINT[NumCurves], Evaluator, IORDRE);
599 }
600#endif
601//Portage HP le compilateur refuse le debranchement
602// fin_de_boucle:
603// {;} fin de la boucle while
604 }
605 return;
606}
607//=======================================================================
608//function : AdvApprox_ApproxAFunction
609//purpose : Constructeur avec Decoupe Dichotomique
610//=======================================================================
611
612AdvApprox_ApproxAFunction::
613AdvApprox_ApproxAFunction(const Standard_Integer Num1DSS,
614 const Standard_Integer Num2DSS,
615 const Standard_Integer Num3DSS,
857ffd5e 616 const Handle(TColStd_HArray1OfReal)& OneDTol,
617 const Handle(TColStd_HArray1OfReal)& TwoDTol,
618 const Handle(TColStd_HArray1OfReal)& ThreeDTol,
7fd59977 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),
628 myFirst(First),
629 myLast(Last),
630 myContinuity(Continuity),
631 myMaxDegree(MaxDeg),
632 myMaxSegments(MaxSeg),
633 myDone(Standard_False),
634 myHasResult(Standard_False),
635 myEvaluator((Standard_Address)&Func)
636{
637 AdvApprox_DichoCutting Cut;
638 Perform(Num1DSS, Num2DSS, Num3DSS, Cut);
639}
640
641AdvApprox_ApproxAFunction::
642AdvApprox_ApproxAFunction(const Standard_Integer Num1DSS,
643 const Standard_Integer Num2DSS,
644 const Standard_Integer Num3DSS,
857ffd5e 645 const Handle(TColStd_HArray1OfReal)& OneDTol,
646 const Handle(TColStd_HArray1OfReal)& TwoDTol,
647 const Handle(TColStd_HArray1OfReal)& ThreeDTol,
7fd59977 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),
658 myFirst(First),
659 myLast(Last),
660 myContinuity(Continuity),
661 myMaxDegree(MaxDeg),
662 myMaxSegments(MaxSeg),
663 myDone(Standard_False),
664 myHasResult(Standard_False),
665 myEvaluator((Standard_Address)&Func)
666{
667 Perform(Num1DSS, Num2DSS, Num3DSS, CutTool);
668}
669
670//=======================================================================
671//function : AdvApprox_ApproxAFunction::Perform
672//purpose : Make all the Work !!
673//=======================================================================
674void AdvApprox_ApproxAFunction::Perform(const Standard_Integer Num1DSS,
675 const Standard_Integer Num2DSS,
676 const Standard_Integer Num3DSS,
677 const AdvApprox_Cutting& CutTool)
678{
679 if (Num1DSS < 0 ||
680 Num2DSS < 0 ||
681 Num3DSS < 0 ||
682 Num1DSS + Num2DSS + Num3DSS <= 0 ||
683 myLast < myFirst ||
684 myMaxDegree < 1 ||
685 myMaxSegments < 0)
9775fa61 686 throw Standard_ConstructionError();
7fd59977 687 if (myMaxDegree > 14) {
688 myMaxDegree = 14 ;
689 }
690 //
691 // Input
692 //
693 myNumSubSpaces[0] = Num1DSS ;
694 myNumSubSpaces[1] = Num2DSS ;
695 myNumSubSpaces[2] = Num3DSS ;
696 Standard_Integer TotalNumSS =
697 Num1DSS + Num2DSS + Num3DSS,
698 ii,
699 jj,
700 kk,
701 index,
702 dim_index,
703 local_index;
704 Standard_Integer TotalDimension =
705 myNumSubSpaces[0] + 2 * myNumSubSpaces[1] + 3 * myNumSubSpaces[2] ;
706 Standard_Real error_value ;
707
708 Standard_Integer ContinuityOrder=0 ;
709 switch (myContinuity) {
710 case GeomAbs_C0:
711 ContinuityOrder = 0 ;
712 break ;
713 case GeomAbs_C1:
714 ContinuityOrder = 1 ;
715 break ;
716 case GeomAbs_C2:
717 ContinuityOrder = 2 ;
718 break ;
719 default:
9775fa61 720 throw Standard_ConstructionError();
7fd59977 721 }
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 ;
726 //
727 // WARNING : the polynomial coefficients are the
728 // taylor expansion of the polynomial at 0.0e0 !
729 //
730 ApproxStartEnd[0] = -1.0e0 ;
731 ApproxStartEnd[1] = 1.0e0 ;
732
733 TColStd_Array1OfInteger LocalDimension(1,TotalNumSS) ;
734
735
736 index = 1 ;
737 TColStd_Array1OfReal LocalTolerances(1,TotalNumSS) ;
738
739 for(jj = 1; jj <= myNumSubSpaces[0] ; jj++) {
740 LocalTolerances.SetValue(index,my1DTolerances->Value(jj)) ;
741 LocalDimension.SetValue(index,1) ;
742 index += 1 ;
743 }
744 for(jj = 1 ; jj <= myNumSubSpaces[1] ; jj++) {
745 LocalTolerances.SetValue(index,my2DTolerances->Value(jj)) ;
746 LocalDimension.SetValue(index,2) ;
747 index += 1 ;
748 }
749 for(jj = 1; jj <= myNumSubSpaces[2] ; jj++) {
750 LocalTolerances.SetValue(index,my3DTolerances->Value(jj)) ;
751 LocalDimension.SetValue(index,3) ;
752 index += 1 ;
753 }
754 //
755 // Ouput
756 //
757
758 Standard_Integer ErrorCode = 0,
759 NumCurves,
760 size =
761 myMaxSegments * NumMaxCoeffs * TotalDimension ;
762 Handle(TColStd_HArray1OfInteger) NumCoeffPerCurvePtr =
763 new TColStd_HArray1OfInteger (1,myMaxSegments) ;
764
765 Handle(TColStd_HArray1OfReal) LocalCoefficientsPtr =
766 new TColStd_HArray1OfReal(1,size) ;
767
768 Handle(TColStd_HArray1OfReal) IntervalsPtr =
769 new TColStd_HArray1OfReal (1,myMaxSegments+1) ;
770
771 TColStd_Array1OfReal ErrorMax(1,myMaxSegments * TotalNumSS) ;
772
773 TColStd_Array1OfReal AverageError(1,myMaxSegments * TotalNumSS) ;
774
32ca7a51 775 Approximation ( TotalDimension,
776 TotalNumSS, LocalDimension,
777 myFirst, myLast,
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
787 ErrorCode) ;
788
7fd59977 789 if (ErrorCode == 0 || ErrorCode == -1) {
790 //
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
793
794 TColStd_Array1OfInteger TabContinuity(1, NumCurves);
795 TColStd_Array2OfReal PolynomialIntervalsPtr (1,NumCurves,1,2);
796 for (ii = PolynomialIntervalsPtr.LowerRow() ;
797 ii <= PolynomialIntervalsPtr.UpperRow() ;
798 ii++) {
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]) ;
804 }
805
806 PrepareConvert(NumCurves, myMaxDegree, ContinuityOrder,
807 Num1DSS, Num2DSS, Num3DSS,
808 NumCoeffPerCurvePtr->Array1(),
809 LocalCoefficientsPtr->ChangeArray1(),
810 PolynomialIntervalsPtr, IntervalsPtr->Array1(),
811 LocalTolerances, ErrorMax,
812 TabContinuity);
813
814 Convert_CompPolynomialToPoles
815 AConverter(NumCurves,
816 TotalDimension,
817 myMaxDegree,
818 TabContinuity,
819 NumCoeffPerCurvePtr->Array1(),
820 LocalCoefficientsPtr->Array1(),
821 PolynomialIntervalsPtr,
822 IntervalsPtr->Array1());
823
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() ;
830 index = 0 ;
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)) ;
839 }
840 }
841
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) ;
847
848 }
849 my1DMaxError->SetValue(jj, error_value) ;
850 }
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) ;
856
857 }
858 error_value /= (Standard_Real) NumCurves ;
859 my1DAverageError->SetValue(jj, error_value) ;
860 }
861 index += myNumSubSpaces[0] ;
862 }
863
864 dim_index = index; //Pour le cas ou il n'y a pas de 2D
865
866 if (myNumSubSpaces[1] > 0) {
867 gp_Pnt2d Point2d ;
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++) {
876 Point2d.SetCoord(kk,
877 PolesPtr->Value(ii,local_index + kk)) ;
878 }
879 my2DPoles->SetValue(ii,jj, Point2d) ;
880 }
881 }
882
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) ;
888
889 }
890 my2DMaxError->SetValue(jj, error_value) ;
891 }
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) ;
897
898 }
899 error_value /= (Standard_Real) NumCurves ;
900 my2DAverageError->SetValue(jj, error_value) ;
901 }
902 index += myNumSubSpaces[1] ;
903 // Pour les poles il faut doubler le decalage :
904 dim_index = index + myNumSubSpaces[1];
905 }
906
907 if (myNumSubSpaces[2] > 0) {
908 gp_Pnt Point ;
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++) {
917 Point.SetCoord(kk,
918 PolesPtr->Value(ii,local_index + kk)) ;
919 }
920 my3DPoles->SetValue(ii,jj,Point) ;
921
922 }
923 }
924
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) ;
930
931 }
932 my3DMaxError->SetValue(jj, error_value) ;
933 }
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) ;
939
940 }
941 error_value /= (Standard_Real) NumCurves ;
942 my3DAverageError->SetValue(jj, error_value) ;
943 }
944 }
945 if (ErrorCode == 0) {
946 myDone = Standard_True ;
947 myHasResult = Standard_True ;
948 }
949 else if (ErrorCode == -1) {
950 myHasResult = Standard_True ;
951 }
952
953 }
954 }
955}
956
957
958//=======================================================================
959//function : Poles
960//purpose :
961//=======================================================================
962
963void AdvApprox_ApproxAFunction::Poles(const Standard_Integer Index,
964 TColgp_Array1OfPnt& P) const
965{
966 Standard_Integer ii ;
967 for (ii = P.Lower() ; ii <= P.Upper() ; ii++) {
968 P.SetValue(ii,my3DPoles->Value(ii,Index)) ;
969 }
970}
971
972
973//=======================================================================
974//function : NbPoles
975//purpose :
976//=======================================================================
977
978Standard_Integer AdvApprox_ApproxAFunction::NbPoles() const
979{ if (myDone || myHasResult) return BSplCLib::NbPoles(myDegree,
980 Standard_False,
981 myMults->Array1()) ;
982 return 0 ; }
983
984//=======================================================================
985//function : Poles
986//purpose :
987//=======================================================================
988
989void AdvApprox_ApproxAFunction::Poles2d(const Standard_Integer Index,
990 TColgp_Array1OfPnt2d& P) const
991{
992 Standard_Integer ii ;
993 for (ii = P.Lower() ; ii <= P.Upper() ; ii++) {
994 P.SetValue(ii,my2DPoles->Value(ii,Index)) ;
995 }
996}
997//=======================================================================
998//function : Poles
999//purpose :
1000//=======================================================================
1001
1002void AdvApprox_ApproxAFunction::Poles1d(const Standard_Integer Index,
1003 TColStd_Array1OfReal& P) const
1004{
1005 Standard_Integer ii ;
1006 for (ii = P.Lower() ; ii <= P.Upper() ; ii++) {
1007 P.SetValue(ii,my1DPoles->Value(ii,Index)) ;
1008 }
1009}
1010//=======================================================================
1011//function : MaxError
1012//purpose :
1013//=======================================================================
1014Handle(TColStd_HArray1OfReal)
1015 AdvApprox_ApproxAFunction::MaxError(const Standard_Integer D) const
1016
1017{
1018 Handle(TColStd_HArray1OfReal) EPtr ;
1019 if (D <= 0 ||
1020 D > 3) {
1021
9775fa61 1022 throw Standard_OutOfRange() ;
7fd59977 1023 }
1024 switch (D) {
1025 case 1:
1026 EPtr = my1DMaxError ;
1027 break ;
1028 case 2:
1029 EPtr = my2DMaxError ;
1030 break ;
1031 case 3:
1032 EPtr = my3DMaxError ;
1033 break ;
1034 }
1035 return EPtr ;
1036}
1037//=======================================================================
1038//function : MaxError
1039//purpose :
1040//=======================================================================
1041Standard_Real AdvApprox_ApproxAFunction::MaxError(const Standard_Integer D,
1042 const Standard_Integer Index) const
1043{
1044 Handle(TColStd_HArray1OfReal) EPtr =
1045 MaxError(D) ;
1046
1047 return EPtr->Value(Index) ;
1048}
1049
1050//=======================================================================
1051//function : AverageError
1052//purpose :
1053//=======================================================================
1054Handle(TColStd_HArray1OfReal)
1055 AdvApprox_ApproxAFunction::AverageError(const Standard_Integer D) const
1056
1057{
1058 Handle(TColStd_HArray1OfReal) EPtr ;
1059 if (D <= 0 ||
1060 D > 3) {
1061
9775fa61 1062 throw Standard_OutOfRange() ;
7fd59977 1063 }
1064 switch (D) {
1065 case 1:
1066 EPtr = my1DAverageError ;
1067 break ;
1068 case 2:
1069 EPtr = my2DAverageError ;
1070 break ;
1071 case 3:
1072 EPtr = my3DAverageError ;
1073 break ;
1074 }
1075 return EPtr ;
1076}
1077//=======================================================================
1078//function : AverageError
1079//purpose :
1080//=======================================================================
1081
1082Standard_Real AdvApprox_ApproxAFunction::AverageError(
1083 const Standard_Integer D,
1084 const Standard_Integer Index) const
1085
1086{
1087 Handle(TColStd_HArray1OfReal) EPtr =
1088 AverageError(D) ;
1089 return EPtr->Value(Index) ;
1090}
1091
1092void AdvApprox_ApproxAFunction::Dump(Standard_OStream& o) const
1093{
1094 Standard_Integer ii;
04232180 1095 o << "Dump of ApproxAFunction" << std::endl;
7fd59977 1096 if (myNumSubSpaces[0] > 0) {
04232180 1097 o << "Error(s) 1d = " << std::endl;
7fd59977 1098 for (ii=1; ii <= myNumSubSpaces[0]; ii++) {
04232180 1099 o << " " << MaxError(1, ii) << std::endl;
7fd59977 1100 }
1101 }
1102
1103 if (myNumSubSpaces[1] > 0) {
04232180 1104 o << "Error(s) 2d = " << std::endl;
7fd59977 1105 for (ii=1; ii <= myNumSubSpaces[1]; ii++) {
04232180 1106 o << " " << MaxError(2, ii) << std::endl;
7fd59977 1107 }
1108 }
1109
1110 if (myNumSubSpaces[2] > 0) {
04232180 1111 o << "Error(s) 3d = " << std::endl;
7fd59977 1112 for (ii=1; ii <= myNumSubSpaces[2]; ii++) {
04232180 1113 o << " " << MaxError(3, ii) << std::endl;
7fd59977 1114 }
1115 }
1116}
Open CASCADE Technology repository