| 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 |
| 5 | // |
| 6 | // This file is part of Open CASCADE Technology software library. |
| 7 | // |
| 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. |
| 13 | // |
| 14 | // Alternatively, this file may be used under the terms of Open CASCADE |
| 15 | // commercial license or contractual agreement. |
| 16 | |
| 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 |
| 23 | |
| 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> |
| 32 | #include <gp_Vec.hxx> |
| 33 | #include <gp_Vec2d.hxx> |
| 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> |
| 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> |
| 49 | |
| 50 | #ifdef OCCT_DEBUG |
| 51 | |
| 52 | static Standard_Boolean AdvApprox_Debug = 0; |
| 53 | |
| 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 | ///=================================================== |
| 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) ) { |
| 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; |
| 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) ) { |
| 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; |
| 119 | } |
| 120 | } |
| 121 | } |
| 122 | #endif |
| 123 | |
| 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 | //==================================================================== |
| 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 | //====================================================================== |
| 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) |
| 422 | { |
| 423 | // Standard_Real EpsPar = Precision::Confusion(); |
| 424 | Standard_Integer IDIM, NUPIL,TheDeg; |
| 425 | #ifdef OCCT_DEBUG |
| 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: |
| 462 | throw Standard_ConstructionError(); |
| 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 | |
| 592 | #ifdef OCCT_DEBUG |
| 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 | |
| 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), |
| 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 | |
| 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), |
| 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 | //======================================================================= |
| 674 | void 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) |
| 686 | throw Standard_ConstructionError(); |
| 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: |
| 720 | throw Standard_ConstructionError(); |
| 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 | |
| 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 | |
| 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 | |
| 963 | void 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 | |
| 978 | Standard_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 | |
| 989 | void 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 | |
| 1002 | void 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 | //======================================================================= |
| 1014 | Handle(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 | |
| 1022 | throw Standard_OutOfRange() ; |
| 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 | //======================================================================= |
| 1041 | Standard_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 | //======================================================================= |
| 1054 | Handle(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 | |
| 1062 | throw Standard_OutOfRange() ; |
| 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 | |
| 1082 | Standard_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 | |
| 1092 | void AdvApprox_ApproxAFunction::Dump(Standard_OStream& o) const |
| 1093 | { |
| 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; |
| 1100 | } |
| 1101 | } |
| 1102 | |
| 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; |
| 1107 | } |
| 1108 | } |
| 1109 | |
| 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; |
| 1114 | } |
| 1115 | } |
| 1116 | } |