src/GeomConvert/GeomConvert_ApproxSurface.cxx

   1 // Copyright (c) 1997-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   6 // This library is free software; you can redistribute it and/or modify it under
   7 // the terms of the GNU Lesser General Public License version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15
  16 #include <Adaptor3d_HSurface.hxx>
  17 #include <AdvApp2Var_ApproxAFunc2Var.hxx>
  18 #include <AdvApprox_PrefAndRec.hxx>
  19 #include <Geom_BSplineSurface.hxx>
  20 #include <Geom_Surface.hxx>
  21 #include <GeomAdaptor_HSurface.hxx>
  22 #include <GeomConvert_ApproxSurface.hxx>
  23 #include <Precision.hxx>
  24 #include <Standard_OutOfRange.hxx>
  25 #include <Standard_Real.hxx>
  26 #include <TColStd_HArray1OfReal.hxx>
  27 #include <TColStd_HArray2OfReal.hxx>
  28
  29 class GeomConvert_ApproxSurface_Eval : public AdvApp2Var_EvaluatorFunc2Var
  30 {
  31
  32 public:
  33
  34   GeomConvert_ApproxSurface_Eval (const Handle(Adaptor3d_HSurface)& theAdaptor)
  35   : myAdaptor (theAdaptor) {}
  36
  37   virtual void Evaluate (Standard_Integer* theDimension,
  38                          Standard_Real*    theUStartEnd,
  39                          Standard_Real*    theVStartEnd,
  40                          Standard_Integer* theFavorIso,
  41                          Standard_Real*    theConstParam,
  42                          Standard_Integer* theNbParams,
  43                          Standard_Real*    theParameters,
  44                          Standard_Integer* theUOrder,
  45                          Standard_Integer* theVOrder,
  46                          Standard_Real*    theResult,
  47                          Standard_Integer* theErrorCode) const;
  48
  49 private:
  50
  51   mutable Handle(Adaptor3d_HSurface) myAdaptor;
  52
  53 };
  54
  55
  56 void GeomConvert_ApproxSurface_Eval::Evaluate (Standard_Integer * Dimension,
  57                            // Dimension
  58                            Standard_Real    * UStartEnd,
  59                            // StartEnd[2] in U
  60                            Standard_Real    * VStartEnd,
  61                            // StartEnd[2] in V
  62                            Standard_Integer * FavorIso,
  63                            // Choice of constante, 1 for U, 2 for V
  64                            Standard_Real    * ConstParam,
  65                            // Value of constant parameter
  66                            Standard_Integer * NbParams,
  67                            // Number of parameters N
  68                            Standard_Real    * Parameters,
  69                            // Values of parameters,
  70                            Standard_Integer * UOrder,
  71                            // Derivative Request in U
  72                            Standard_Integer * VOrder,
  73                            // Derivative Request in V
  74                            Standard_Real    * Result,
  75                            // Result[Dimension,N]
  76                            Standard_Integer * ErrorCode) const
  77                            // Error Code
  78 {
  79   *ErrorCode = 0;
  80 //  Standard_Integer idim;
  81   Standard_Integer jpar;
  82   Standard_Real Upar,Vpar;
  83
  84 // Dimension incorrecte
  85   if (*Dimension!=3) {
  86     *ErrorCode = 1;
  87   }
  88
  89 // Parametres incorrects
  90 /* if (*FavorIso==1) {
  91     Upar = *ConstParam;
  92     if (( Upar < UStartEnd[0] ) || ( Upar > UStartEnd[1] )) {
  93       *ErrorCode = 2;
  94     }
  95     for (jpar=1;jpar<=*NbParams;jpar++) {
  96       Vpar = Parameters[jpar-1];
  97       if (( Vpar < VStartEnd[0] ) || ( Vpar > VStartEnd[1] )) {
  98         *ErrorCode = 2;
  99       }
 100     }
 101  }
 102  else {
 103     Vpar = *ConstParam;
 104     if (( Vpar < VStartEnd[0] ) || ( Vpar > VStartEnd[1] )) {
 105       *ErrorCode = 2;
 106     }
 107     for (jpar=1;jpar<=*NbParams;jpar++) {
 108       Upar = Parameters[jpar-1];
 109       if (( Upar < UStartEnd[0] ) || ( Upar > UStartEnd[1] )) {
 110         *ErrorCode = 2;
 111       }
 112     }
 113  }*/
 114
 115 // Initialisation
 116
 117   myAdaptor = myAdaptor->UTrim (UStartEnd[0], UStartEnd[1], Precision::PConfusion());
 118   myAdaptor = myAdaptor->VTrim (VStartEnd[0], VStartEnd[1], Precision::PConfusion());
 119 /*
 120   for (idim=1;idim<=*Dimension;idim++) {
 121     for (jpar=1;jpar<=*NbParams;jpar++) {
 122       Result[idim-1+(jpar-1)*(*Dimension)] = 0.;
 123     }
 124   }*/
 125
 126
 127  Standard_Integer Order = *UOrder + *VOrder;
 128  gp_Pnt pnt;
 129  gp_Vec vect, v1, v2, v3, v4, v5, v6, v7, v8, v9;
 130
 131  if (*FavorIso==1) {
 132   Upar = *ConstParam;
 133   switch (Order) {
 134   case 0 :
 135     for (jpar=1;jpar<=*NbParams;jpar++) {
 136         Vpar = Parameters[jpar-1];
 137         pnt = myAdaptor->Value (Upar, Vpar);
 138         Result[(jpar-1)*(*Dimension)] = pnt.X();
 139         Result[1+(jpar-1)*(*Dimension)] = pnt.Y();
 140         Result[2+(jpar-1)*(*Dimension)] = pnt.Z();
 141     }
 142     break;
 143   case 1 :
 144     for (jpar=1;jpar<=*NbParams;jpar++) {
 145         Vpar = Parameters[jpar-1];
 146         myAdaptor->D1 (Upar, Vpar, pnt, v1, v2);
 147         if (*UOrder==1) {
 148           Result[(jpar-1)*(*Dimension)] = v1.X();
 149           Result[1+(jpar-1)*(*Dimension)] = v1.Y();
 150           Result[2+(jpar-1)*(*Dimension)] = v1.Z();
 151         }
 152         else {
 153           Result[(jpar-1)*(*Dimension)] = v2.X();
 154           Result[1+(jpar-1)*(*Dimension)] = v2.Y();
 155           Result[2+(jpar-1)*(*Dimension)] = v2.Z();
 156         }
 157     }
 158     break;
 159   case 2 :
 160     for (jpar=1;jpar<=*NbParams;jpar++) {
 161         Vpar = Parameters[jpar-1];
 162         myAdaptor->D2 (Upar, Vpar, pnt, v1, v2, v3, v4, v5);
 163         if (*UOrder==2) {
 164           Result[(jpar-1)*(*Dimension)] = v3.X();
 165           Result[1+(jpar-1)*(*Dimension)] = v3.Y();
 166           Result[2+(jpar-1)*(*Dimension)] = v3.Z();
 167         }
 168         else if (*UOrder==1) {
 169           Result[(jpar-1)*(*Dimension)] = v5.X();
 170           Result[1+(jpar-1)*(*Dimension)] = v5.Y();
 171           Result[2+(jpar-1)*(*Dimension)] = v5.Z();
 172         }
 173         else if (*UOrder==0) {
 174           Result[(jpar-1)*(*Dimension)] = v4.X();
 175           Result[1+(jpar-1)*(*Dimension)] = v4.Y();
 176           Result[2+(jpar-1)*(*Dimension)] = v4.Z();
 177         }
 178     }
 179     break;
 180   case 3 :
 181     for (jpar=1;jpar<=*NbParams;jpar++) {
 182         Vpar = Parameters[jpar-1];
 183         myAdaptor->D3 (Upar, Vpar, pnt, v1, v2, v3, v4, v5, v6, v7, v8, v9);
 184         if (*UOrder==2) {
 185           Result[(jpar-1)*(*Dimension)] = v8.X();
 186           Result[1+(jpar-1)*(*Dimension)] = v8.Y();
 187           Result[2+(jpar-1)*(*Dimension)] = v8.Z();
 188         }
 189         else if (*UOrder==1) {
 190           Result[(jpar-1)*(*Dimension)] = v9.X();
 191           Result[1+(jpar-1)*(*Dimension)] = v9.Y();
 192           Result[2+(jpar-1)*(*Dimension)] = v9.Z();
 193         }
 194     }
 195     break;
 196   case 4 :
 197     for (jpar=1;jpar<=*NbParams;jpar++) {
 198         Vpar = Parameters[jpar-1];
 199         vect = myAdaptor->DN (Upar, Vpar, *UOrder, *VOrder);
 200         Result[(jpar-1)*(*Dimension)] = vect.X();
 201         Result[1+(jpar-1)*(*Dimension)] = vect.Y();
 202         Result[2+(jpar-1)*(*Dimension)] = vect.Z();
 203     }
 204     break;
 205   }
 206  }
 207  else {
 208   Vpar = *ConstParam;
 209   switch (Order) {
 210   case 0 :
 211     for (jpar=1;jpar<=*NbParams;jpar++) {
 212         Upar = Parameters[jpar-1];
 213         pnt = myAdaptor->Value (Upar, Vpar);
 214         Result[(jpar-1)*(*Dimension)] = pnt.X();
 215         Result[1+(jpar-1)*(*Dimension)] = pnt.Y();
 216         Result[2+(jpar-1)*(*Dimension)] = pnt.Z();
 217     }
 218     break;
 219   case 1 :
 220     for (jpar=1;jpar<=*NbParams;jpar++) {
 221         Upar = Parameters[jpar-1];
 222         myAdaptor->D1 (Upar, Vpar, pnt, v1, v2);
 223         if (*UOrder==1) {
 224           Result[(jpar-1)*(*Dimension)] = v1.X();
 225           Result[1+(jpar-1)*(*Dimension)] = v1.Y();
 226           Result[2+(jpar-1)*(*Dimension)] = v1.Z();
 227         }
 228         else {
 229           Result[(jpar-1)*(*Dimension)] = v2.X();
 230           Result[1+(jpar-1)*(*Dimension)] = v2.Y();
 231           Result[2+(jpar-1)*(*Dimension)] = v2.Z();
 232         }
 233     }
 234     break;
 235   case 2 :
 236     for (jpar=1;jpar<=*NbParams;jpar++) {
 237         Upar = Parameters[jpar-1];
 238         myAdaptor->D2 (Upar, Vpar, pnt, v1, v2, v3, v4, v5);
 239         if (*UOrder==2) {
 240           Result[(jpar-1)*(*Dimension)] = v3.X();
 241           Result[1+(jpar-1)*(*Dimension)] = v3.Y();
 242           Result[2+(jpar-1)*(*Dimension)] = v3.Z();
 243         }
 244         else if (*UOrder==1) {
 245           Result[(jpar-1)*(*Dimension)] = v5.X();
 246           Result[1+(jpar-1)*(*Dimension)] = v5.Y();
 247           Result[2+(jpar-1)*(*Dimension)] = v5.Z();
 248         }
 249         else if (*UOrder==0) {
 250           Result[(jpar-1)*(*Dimension)] = v4.X();
 251           Result[1+(jpar-1)*(*Dimension)] = v4.Y();
 252           Result[2+(jpar-1)*(*Dimension)] = v4.Z();
 253         }
 254     }
 255     break;
 256   case 3 :
 257     for (jpar=1;jpar<=*NbParams;jpar++) {
 258         Upar = Parameters[jpar-1];
 259         myAdaptor->D3 (Upar, Vpar, pnt, v1, v2, v3, v4, v5, v6, v7, v8, v9);
 260         if (*UOrder==2) {
 261           Result[(jpar-1)*(*Dimension)] = v8.X();
 262           Result[1+(jpar-1)*(*Dimension)] = v8.Y();
 263           Result[2+(jpar-1)*(*Dimension)] = v8.Z();
 264         }
 265         else if (*UOrder==1) {
 266           Result[(jpar-1)*(*Dimension)] = v9.X();
 267           Result[1+(jpar-1)*(*Dimension)] = v9.Y();
 268           Result[2+(jpar-1)*(*Dimension)] = v9.Z();
 269         }
 270     }
 271     break;
 272   case 4 :
 273     for (jpar=1;jpar<=*NbParams;jpar++) {
 274         Upar = Parameters[jpar-1];
 275         vect = myAdaptor->DN (Upar, Vpar, *UOrder, *VOrder);
 276         Result[(jpar-1)*(*Dimension)] = vect.X();
 277         Result[1+(jpar-1)*(*Dimension)] = vect.Y();
 278         Result[2+(jpar-1)*(*Dimension)] = vect.Z();
 279     }
 280     break;
 281   }
 282  }
 283
 284 }
 285
 286
 287 //=======================================================================
 288 //function : GeomConvert_ApproxSurface
 289 //purpose  :
 290 //=======================================================================
 291
 292 GeomConvert_ApproxSurface::GeomConvert_ApproxSurface(const Handle(Geom_Surface)& Surf,
 293                                                      const Standard_Real Tol3d,
 294                                                      const GeomAbs_Shape UContinuity,
 295                                                      const GeomAbs_Shape VContinuity,
 296                                                      const Standard_Integer MaxDegU,
 297                                                      const Standard_Integer MaxDegV,
 298                                                      const Standard_Integer MaxSegments,
 299                                                      const Standard_Integer PrecisCode)
 300 {
 301   Handle(Adaptor3d_HSurface) aSurfAdaptor = new GeomAdaptor_HSurface (Surf);
 302   Approximate(aSurfAdaptor, Tol3d, UContinuity, VContinuity, MaxDegU, MaxDegV, MaxSegments, PrecisCode);
 303 }
 304
 305 GeomConvert_ApproxSurface::GeomConvert_ApproxSurface(const Handle(Adaptor3d_HSurface)& Surf,
 306                                                      const Standard_Real Tol3d,
 307                                                      const GeomAbs_Shape UContinuity,
 308                                                      const GeomAbs_Shape VContinuity,
 309                                                      const Standard_Integer MaxDegU,
 310                                                      const Standard_Integer MaxDegV,
 311                                                      const Standard_Integer MaxSegments,
 312                                                      const Standard_Integer PrecisCode)
 313 {
 314   Approximate(Surf, Tol3d, UContinuity, VContinuity, MaxDegU, MaxDegV, MaxSegments, PrecisCode);
 315 }
 316
 317 void GeomConvert_ApproxSurface::Approximate(const Handle(Adaptor3d_HSurface)& theSurf,
 318                                             const Standard_Real theTol3d,
 319                                             const GeomAbs_Shape theUContinuity,
 320                                             const GeomAbs_Shape theVContinuity,
 321                                             const Standard_Integer theMaxDegU,
 322                                             const Standard_Integer theMaxDegV,
 323                                             const Standard_Integer theMaxSegments,
 324                                             const Standard_Integer thePrecisCode)
 325 {
 326   Standard_Real U0 = theSurf->FirstUParameter();
 327   Standard_Real U1 = theSurf->LastUParameter();
 328   Standard_Real V0 = theSurf->FirstVParameter();
 329   Standard_Real V1 = theSurf->LastVParameter();
 330
 331 // " Init des nombres de sous-espaces et des tolerances"
 332   Standard_Integer nb1 = 0, nb2 = 0, nb3 = 1;
 333   Handle(TColStd_HArray1OfReal) nul1 = new TColStd_HArray1OfReal(1,1);
 334   nul1->SetValue(1, 0.);
 335   Handle(TColStd_HArray2OfReal) nul2 = new TColStd_HArray2OfReal(1,1,1,4);
 336   nul2->SetValue(1, 1, 0.);
 337   nul2->SetValue(1, 2, 0.);
 338   nul2->SetValue(1, 3, 0.);
 339   nul2->SetValue(1, 4, 0.);
 340   Handle(TColStd_HArray1OfReal) eps3D = new TColStd_HArray1OfReal(1,1);
 341   eps3D->SetValue(1, theTol3d);
 342   Handle(TColStd_HArray2OfReal) epsfr = new TColStd_HArray2OfReal(1,1,1,4);
 343   epsfr->SetValue(1, 1, theTol3d);
 344   epsfr->SetValue(1, 2, theTol3d);
 345   epsfr->SetValue(1, 3, theTol3d);
 346   epsfr->SetValue(1, 4, theTol3d);
 347
 348 // " Init du type d'iso"
 349   GeomAbs_IsoType IsoType = GeomAbs_IsoV;
 350   Standard_Integer NbDec;
 351
 352   NbDec = theSurf->NbUIntervals(GeomAbs_C2);
 353   TColStd_Array1OfReal UDec_C2(1, NbDec+1);
 354   theSurf->UIntervals(UDec_C2, GeomAbs_C2);
 355   NbDec = theSurf->NbVIntervals(GeomAbs_C2);
 356   TColStd_Array1OfReal VDec_C2(1, NbDec+1);
 357   theSurf->VIntervals(VDec_C2, GeomAbs_C2);
 358
 359   NbDec = theSurf->NbUIntervals(GeomAbs_C3);
 360   TColStd_Array1OfReal UDec_C3(1, NbDec+1);
 361   theSurf->UIntervals(UDec_C3, GeomAbs_C3);
 362
 363   NbDec = theSurf->NbVIntervals(GeomAbs_C3);
 364   TColStd_Array1OfReal VDec_C3(1, NbDec+1);
 365   theSurf->VIntervals(VDec_C3, GeomAbs_C3);
 366   // Approximation avec decoupe preferentiel
 367   // aux lieux de discontinuitees C2
 368   AdvApprox_PrefAndRec pUDec(UDec_C2,UDec_C3);
 369   AdvApprox_PrefAndRec pVDec(VDec_C2,VDec_C3);
 370
 371 //POP pour WNT
 372   GeomConvert_ApproxSurface_Eval ev (theSurf);
 373   AdvApp2Var_ApproxAFunc2Var approx(nb1, nb2, nb3,
 374                                     nul1,nul1,eps3D,
 375                                     nul2,nul2,epsfr,
 376                                     U0,U1,V0,V1,
 377                                     IsoType,theUContinuity,theVContinuity,thePrecisCode,
 378                                     theMaxDegU,theMaxDegV,theMaxSegments,ev,
 379                                     pUDec,pVDec);
 380
 381   myMaxError  = approx.MaxError(3,1);
 382   myBSplSurf  = approx.Surface(1);
 383   myIsDone    = approx.IsDone();
 384   myHasResult = approx.HasResult();
 385 }
 386
 387
 388 //=======================================================================
 389 //function : Surface
 390 //purpose  :
 391 //=======================================================================
 392
 393  Handle(Geom_BSplineSurface) GeomConvert_ApproxSurface::Surface() const
 394 {
 395   return myBSplSurf;
 396 }
 397
 398 //=======================================================================
 399 //function : IsDone
 400 //purpose  :
 401 //=======================================================================
 402
 403  Standard_Boolean GeomConvert_ApproxSurface::IsDone() const
 404 {
 405   return myIsDone;
 406 }
 407
 408 //=======================================================================
 409 //function : HasResult
 410 //purpose  :
 411 //=======================================================================
 412
 413  Standard_Boolean GeomConvert_ApproxSurface::HasResult() const
 414 {
 415   return myHasResult;
 416 }
 417
 418 //=======================================================================
 419 //function : MaxError
 420 //purpose  :
 421 //=======================================================================
 422
 423  Standard_Real GeomConvert_ApproxSurface::MaxError() const
 424 {
 425   return myMaxError;
 426 }
 427 //=======================================================================
 428 //function : Dump
 429 //purpose  :
 430 //=======================================================================
 431
 432 void GeomConvert_ApproxSurface::Dump(Standard_OStream& o) const
 433 {
 434   o<<endl;
 435   if (!myHasResult) { o<<"No result"<<endl; }
 436   else {
 437     o<<"Result max error :"<< myMaxError <<endl;
 438     }
 439   o<<endl;
 440 }
 441
 442
 443
 444
 445
 446
 447