src/GeomFill/GeomFill_NSections.cxx

   1 // Created on: 1998-12-14
   2 // Created by: Joelle CHAUVET
   3 // Copyright (c) 1998-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:    Fri Jan  8 15:47:20 1999
  18 //              enfin un calcul exact pour D1 et D2
  19 //              le calcul par differ. finies est garde dans verifD1 et verifD2
  20 // Modified:    Mon Jan 18 11:06:46 1999
  21 //              mise au point de D1, D2 et IsConstant
  22
  23 #include <BSplCLib.hxx>
  24 #include <Convert_ParameterisationType.hxx>
  25 #include <GCPnts_AbscissaPoint.hxx>
  26 #include <Geom_BSplineCurve.hxx>
  27 #include <Geom_BSplineSurface.hxx>
  28 #include <Geom_Circle.hxx>
  29 #include <Geom_Curve.hxx>
  30 #include <Geom_Geometry.hxx>
  31 #include <Geom_Surface.hxx>
  32 #include <Geom_TrimmedCurve.hxx>
  33 #include <GeomAdaptor_Curve.hxx>
  34 #include <GeomAdaptor_Surface.hxx>
  35 #include <GeomConvert.hxx>
  36 #include <GeomFill_AppSurf.hxx>
  37 #include <GeomFill_Line.hxx>
  38 #include <GeomFill_NSections.hxx>
  39 #include <GeomFill_SectionGenerator.hxx>
  40 #include <gp_Circ.hxx>
  41 #include <gp_Lin.hxx>
  42 #include <gp_Pnt.hxx>
  43 #include <Precision.hxx>
  44 #include <Standard_OutOfRange.hxx>
  45 #include <Standard_Type.hxx>
  46 #include <TColGeom_Array1OfCurve.hxx>
  47 #include <TColgp_Array2OfPnt.hxx>
  48 #include <TColStd_Array1OfInteger.hxx>
  49 #include <TColStd_Array1OfReal.hxx>
  50
  51 #include <stdio.h>
  52 IMPLEMENT_STANDARD_RTTIEXT(GeomFill_NSections,GeomFill_SectionLaw)
  53
  54 #ifdef OCCT_DEBUG
  55 # ifdef DRAW
  56 #  include <DrawTrSurf.hxx>
  57 #include <Geom_Curve.hxx>
  58 # endif
  59 static Standard_Boolean Affich = 0;
  60 static Standard_Integer NbSurf = 0;
  61 #endif
  62
  63 #ifdef OCCT_DEBUG
  64 // verification des fonctions de derivation D1 et D2 par differences finies
  65 static Standard_Boolean verifD1(const TColgp_Array1OfPnt& P1,
  66                          const TColStd_Array1OfReal& W1,
  67                          const TColgp_Array1OfPnt& P2,
  68                          const TColStd_Array1OfReal& W2,
  69                          const TColgp_Array1OfVec& DPoles,
  70                          const TColStd_Array1OfReal& DWeights,
  71                          const Standard_Real pTol,
  72                          const Standard_Real wTol,
  73                          const Standard_Real pas)
  74 {
  75   Standard_Boolean ok = Standard_True;
  76   Standard_Integer ii, L =  P1.Length();
  77   Standard_Real dw;
  78   gp_Vec dP;
  79   for (ii=1; ii<=L; ii++) {
  80     dw = (W2(ii)-W1(ii)) / pas;
  81     if (Abs(dw-DWeights(ii))>wTol) {
  82       if (Affich) {
  83         std::cout<<"erreur dans la derivee 1ere du poids pour l'indice "<<ii<<std::endl;
  84         std::cout<<"par diff finies : "<<dw<<std::endl;
  85         std::cout<<"resultat obtenu : "<<DWeights(ii)<<std::endl;
  86       }
  87       ok = Standard_False;
  88     }
  89     dP.SetXYZ( (P2(ii).XYZ()- P1(ii).XYZ()) /pas );
  90     gp_Vec diff = dP - DPoles(ii);
  91     if (diff.Magnitude()>pTol) {
  92       if (Affich) {
  93         std::cout<<"erreur dans la derivee 1ere du pole pour l'indice "<<ii<<std::endl;
  94         std::cout<<"par diff finies : ("<<dP.X()
  95                               <<" "<<dP.Y()
  96                               <<" "<<dP.Z()<<")"<<std::endl;
  97         std::cout<<"resultat obtenu : ("<<DPoles(ii).X()
  98                               <<" "<<DPoles(ii).Y()
  99                               <<" "<<DPoles(ii).Z()<<")"<<std::endl;
 100       }
 101       ok = Standard_False;
 102     }
 103   }
 104   return ok;
 105 }
 106
 107 static Standard_Boolean verifD2(const TColgp_Array1OfVec& DP1,
 108                          const TColStd_Array1OfReal& DW1,
 109                          const TColgp_Array1OfVec& DP2,
 110                          const TColStd_Array1OfReal& DW2,
 111                          const TColgp_Array1OfVec& D2Poles,
 112                          const TColStd_Array1OfReal& D2Weights,
 113                          const Standard_Real pTol,
 114                          const Standard_Real wTol,
 115                          const Standard_Real pas)
 116 {
 117   Standard_Boolean ok = Standard_True;
 118   Standard_Integer ii, L =  DP1.Length();
 119   Standard_Real d2w;
 120   gp_Vec d2P;
 121   for (ii=1; ii<=L; ii++) {
 122     Standard_Real dw1 = DW1(ii), dw2 = DW2(ii);
 123     d2w = (dw2-dw1) / pas;
 124     if (Abs(d2w-D2Weights(ii))>wTol) {
 125       if (Affich) {
 126         std::cout<<"erreur dans la derivee 2nde du poids pour l'indice "<<ii<<std::endl;
 127         std::cout<<"par diff finies : "<<d2w<<std::endl;
 128         std::cout<<"resultat obtenu : "<<D2Weights(ii)<<std::endl;
 129       }
 130       ok = Standard_False;
 131     }
 132     d2P.SetXYZ( (DP2(ii).XYZ()- DP1(ii).XYZ()) /pas );
 133     gp_Vec diff = d2P - D2Poles(ii);
 134     if (diff.Magnitude()>pTol) {
 135       if (Affich) {
 136         std::cout<<"erreur dans la derivee 2nde du pole pour l'indice "<<ii<<std::endl;
 137         std::cout<<"par diff finies : ("<<d2P.X()
 138                               <<" "<<d2P.Y()
 139                               <<" "<<d2P.Z()<<")"<<std::endl;
 140         std::cout<<"resultat obtenu : ("<<D2Poles(ii).X()
 141                               <<" "<<D2Poles(ii).Y()
 142                               <<" "<<D2Poles(ii).Z()<<")"<<std::endl;
 143       }
 144       ok = Standard_False;
 145     }
 146   }
 147   return ok;
 148 }
 149 #endif
 150
 151 // fonction d'evaluation des poles et des poids de mySurface pour D1 et D2
 152 static void ResultEval(const Handle(Geom_BSplineSurface)& surf,
 153                        const Standard_Real V,
 154                        const Standard_Integer deriv,
 155                        TColStd_Array1OfReal& Result)
 156 {
 157   Standard_Boolean rational = surf->IsVRational() ;
 158   Standard_Integer gap = 3;
 159   if ( rational ) gap++;
 160   Standard_Integer Cdeg = surf->VDegree(),
 161   Cdim = surf->NbUPoles() * gap,
 162   NbP = surf->NbVPoles();
 163
 164   //  les noeuds plats
 165   Standard_Integer Ksize = NbP + Cdeg + 1;
 166   TColStd_Array1OfReal FKnots(1,Ksize);
 167   surf->VKnotSequence(FKnots);
 168
 169   //  les poles
 170   Standard_Integer Psize = Cdim * NbP;
 171   TColStd_Array1OfReal SPoles(1,Psize);
 172   Standard_Integer ii, jj, ipole=1;
 173   for (jj=1;jj<=NbP;jj++) {
 174     for (ii=1;ii<=surf->NbUPoles();ii++) {
 175       SPoles(ipole) = surf->Pole(ii,jj).X();
 176       SPoles(ipole+1) = surf->Pole(ii,jj).Y();
 177       SPoles(ipole+2) = surf->Pole(ii,jj).Z();
 178       if (rational) {
 179         SPoles(ipole+3) = surf->Weight(ii,jj);
 180         SPoles(ipole) *= SPoles(ipole+3);
 181         SPoles(ipole+1) *= SPoles(ipole+3);
 182         SPoles(ipole+2) *= SPoles(ipole+3);
 183       }
 184       ipole+=gap;
 185     }
 186   }
 187   Standard_Real * Padr = (Standard_Real *) &SPoles(1);
 188
 189   Standard_Boolean periodic_flag = Standard_False ;
 190   Standard_Integer extrap_mode[2];
 191   extrap_mode[0] = extrap_mode[1] = Cdeg;
 192   TColStd_Array1OfReal  EvalBS(1, Cdim * (deriv+1)) ;
 193   Standard_Real * Eadr = (Standard_Real *) &EvalBS(1) ;
 194   BSplCLib::Eval(V,periodic_flag,deriv,extrap_mode[0],
 195                  Cdeg,FKnots,Cdim,*Padr,*Eadr);
 196
 197   for (ii=1;ii<=Cdim;ii++) {
 198     Result(ii) = EvalBS(ii+deriv*Cdim);
 199   }
 200 }
 201
 202
 203 //=======================================================================
 204 //function : GeomFill_NSections
 205 //purpose  :
 206 //=======================================================================
 207
 208 GeomFill_NSections::GeomFill_NSections(const TColGeom_SequenceOfCurve& NC)
 209 {
 210   mySections = NC;
 211   UFirst = 0.;
 212   ULast = 1.;
 213   VFirst = 0.;
 214   VLast = 1.;
 215   myRefSurf.Nullify();
 216   ComputeSurface();
 217 }
 218
 219 //=======================================================================
 220 //function : GeomFill_NSections
 221 //purpose  :
 222 //=======================================================================
 223
 224 GeomFill_NSections::GeomFill_NSections(const TColGeom_SequenceOfCurve& NC,
 225                                        const TColStd_SequenceOfReal& NP)
 226 {
 227   mySections = NC;
 228   myParams = NP;
 229   UFirst = 0.;
 230   ULast = 1.;
 231   VFirst = 0.;
 232   VLast = 1.;
 233   myRefSurf.Nullify();
 234   ComputeSurface();
 235 }
 236
 237 //=======================================================================
 238 //function : GeomFill_NSections
 239 //purpose  :
 240 //=======================================================================
 241
 242 GeomFill_NSections::GeomFill_NSections(const TColGeom_SequenceOfCurve& NC,
 243                                        const TColStd_SequenceOfReal& NP,
 244                                        const Standard_Real UF,
 245                                        const Standard_Real UL,
 246                                        const Standard_Real VF,
 247                                        const Standard_Real VL)
 248 {
 249   mySections = NC;
 250   myParams = NP;
 251   UFirst = UF;
 252   ULast = UL;
 253   VFirst = VF;
 254   VLast = VL;
 255   myRefSurf.Nullify();
 256   ComputeSurface();
 257 }
 258
 259 //=======================================================================
 260 //function : GeomFill_NSections
 261 //purpose  :
 262 //=======================================================================
 263
 264 GeomFill_NSections::GeomFill_NSections(const TColGeom_SequenceOfCurve& NC,
 265                                        const GeomFill_SequenceOfTrsf& Trsfs,
 266                                        const TColStd_SequenceOfReal& NP,
 267                                        const Standard_Real UF,
 268                                        const Standard_Real UL,
 269                                        const Standard_Real VF,
 270                                        const Standard_Real VL,
 271                                        const Handle(Geom_BSplineSurface)& Surf)
 272 {
 273   mySections = NC;
 274   myTrsfs = Trsfs;
 275   myParams = NP;
 276   UFirst = UF;
 277   ULast = UL;
 278   VFirst = VF;
 279   VLast = VL;
 280   myRefSurf = Surf;
 281   ComputeSurface();
 282 }
 283
 284 //=======================================================
 285 // Purpose :D0
 286 //=======================================================
 287  Standard_Boolean GeomFill_NSections::D0(const Standard_Real V,
 288                                               TColgp_Array1OfPnt& Poles,
 289                                               TColStd_Array1OfReal& Weights)
 290 {
 291   if (mySurface.IsNull()) {
 292     return Standard_False;
 293   }
 294   else {
 295     Handle(Geom_BSplineCurve) Curve
 296       = Handle(Geom_BSplineCurve)::DownCast(mySurface->VIso( V, Standard_False ));
 297     TColgp_Array1OfPnt poles(1,mySurface->NbUPoles());
 298     TColStd_Array1OfReal weights(1,mySurface->NbUPoles());
 299     Curve->Poles(poles);
 300     Curve->Weights(weights);
 301     Standard_Integer ii, L =  Poles.Length();
 302     for (ii=1; ii<=L; ii++) {
 303       Poles(ii).SetXYZ(poles(ii).XYZ());
 304       Weights(ii) = weights(ii);
 305     }
 306   }
 307   return Standard_True;
 308 }
 309
 310 //=======================================================
 311 // Purpose :D1
 312 //=======================================================
 313  Standard_Boolean GeomFill_NSections::D1(const Standard_Real V,
 314                                               TColgp_Array1OfPnt& Poles,
 315                                               TColgp_Array1OfVec& DPoles,
 316                                               TColStd_Array1OfReal& Weights,
 317                                               TColStd_Array1OfReal& DWeights)
 318 {
 319   if (mySurface.IsNull() ) return Standard_False;
 320
 321   Standard_Boolean ok = D0(V,Poles,Weights);
 322   if (!ok) return Standard_False;
 323
 324   Standard_Integer L =  Poles.Length(), derivative_request = 1;
 325   Standard_Boolean rational = mySurface->IsVRational() ;
 326   Standard_Integer gap = 3;
 327   if (rational) gap++;
 328
 329   Standard_Integer dimResult = mySurface->NbUPoles() * gap;
 330   Handle(Geom_BSplineSurface) surf_deper;
 331   if (mySurface->IsVPeriodic()) {
 332     surf_deper = Handle(Geom_BSplineSurface)::DownCast(mySurface->Copy());
 333     surf_deper->SetVNotPeriodic();
 334     dimResult = surf_deper->NbUPoles() * gap;
 335   }
 336   TColStd_Array1OfReal Result(1,dimResult);
 337   if (mySurface->IsVPeriodic()) {
 338     ResultEval(surf_deper,V,derivative_request,Result);
 339   }
 340   else {
 341     ResultEval(mySurface,V,derivative_request,Result);
 342   }
 343
 344
 345   Standard_Real ww, EpsW = 10*Precision::PConfusion();
 346   Standard_Boolean NullWeight = Standard_False;
 347   if (!rational) DWeights.Init(0.);
 348   Standard_Integer indice = 1, ii;
 349
 350   //  recopie des poles du resultat sous forme de points 3D et de poids
 351   for (ii=1; ii<=L && (!NullWeight) ; ii++) {
 352     DPoles(ii).SetX( Result(indice) );
 353     DPoles(ii).SetY( Result(indice+1) );
 354     DPoles(ii).SetZ( Result(indice+2) );
 355     if (rational) {
 356       ww = Weights(ii);
 357       if (ww < EpsW) {
 358         NullWeight = Standard_True;
 359       }
 360       else {
 361         DWeights(ii) = Result(indice+3);
 362         DPoles(ii)
 363           .SetXYZ( ( DPoles(ii).XYZ()-DWeights(ii)*Poles(ii).Coord() ) / ww );
 364       }
 365     }
 366     indice += gap;
 367   }
 368   if (NullWeight) return Standard_False;
 369
 370   // verif par diff finies sous debug sauf pour les surfaces periodiques
 371 #ifdef OCCT_DEBUG
 372   if (!mySurface->IsVPeriodic()) {
 373     Standard_Real pas = 1.e-6, wTol = 1.e-4, pTol = 1.e-3;
 374     Standard_Real V1,V2;
 375     Standard_Boolean ok1,ok2;
 376     TColStd_Array1OfReal W1(1,L),W2(1,L);
 377     TColgp_Array1OfPnt P1(1,L),P2(1,L);
 378     gp_Pnt nul(0.,0.,0.);
 379     W1.Init(0.);
 380     W2.Init(0.);
 381     P1.Init(nul);
 382     P2.Init(nul);
 383
 384     V1 = V;
 385     V2 = V+pas;
 386     ok1 = D0(V1,P1,W1);
 387     ok2 = D0(V2,P2,W2);
 388     if (!ok1 || !ok2) std::cout<<"probleme en D0"<<std::endl;
 389     Standard_Boolean check = verifD1(P1,W1,P2,W2,DPoles,DWeights,pTol,wTol,pas);
 390     if (!check) std::cout<<"D1 incorrecte en V = "<<V<<std::endl;
 391   }
 392 #endif
 393
 394   return Standard_True;
 395 }
 396
 397 //=======================================================
 398 // Purpose :D2
 399 //=======================================================
 400  Standard_Boolean GeomFill_NSections::D2(const Standard_Real V,
 401                                               TColgp_Array1OfPnt& Poles,
 402                                               TColgp_Array1OfVec& DPoles,
 403                                               TColgp_Array1OfVec& D2Poles,
 404                                               TColStd_Array1OfReal& Weights,
 405                                               TColStd_Array1OfReal& DWeights,
 406                                               TColStd_Array1OfReal& D2Weights)
 407 {
 408   if (mySurface.IsNull() ) return Standard_False;
 409
 410   // pb dans BSplCLib::Eval() pour les surfaces rationnelles de degre 1
 411   // si l'ordre de derivation est egal a 2.
 412   if (mySurface->VDegree()<2) return Standard_False;
 413
 414   Standard_Boolean ok = D1(V,Poles,DPoles,Weights,DWeights);
 415   if (!ok) return Standard_False;
 416
 417   Standard_Integer L =  Poles.Length(), derivative_request = 2;
 418   Standard_Boolean rational = mySurface->IsVRational() ;
 419   Standard_Integer gap = 3;
 420   if (rational) gap++;
 421
 422   Standard_Integer dimResult = mySurface->NbUPoles() * gap;
 423   Handle(Geom_BSplineSurface) surf_deper;
 424   if (mySurface->IsVPeriodic()) {
 425     surf_deper = Handle(Geom_BSplineSurface)::DownCast(mySurface->Copy());
 426     surf_deper->SetVNotPeriodic();
 427     dimResult = surf_deper->NbUPoles() * gap;
 428   }
 429   TColStd_Array1OfReal Result(1,dimResult);
 430   if (mySurface->IsVPeriodic()) {
 431     ResultEval(surf_deper,V,derivative_request,Result);
 432   }
 433   else {
 434     ResultEval(mySurface,V,derivative_request,Result);
 435   }
 436
 437
 438   Standard_Real ww, EpsW = 10*Precision::PConfusion();
 439   Standard_Boolean NullWeight = Standard_False;
 440   if (!rational) D2Weights.Init(0.);
 441   Standard_Integer indice = 1, ii;
 442
 443   //  recopie des poles du resultat sous forme de points 3D et de poids
 444   for (ii=1; ii<=L && (!NullWeight) ; ii++) {
 445     D2Poles(ii).SetX( Result(indice) );
 446     D2Poles(ii).SetY( Result(indice+1) );
 447     D2Poles(ii).SetZ( Result(indice+2) );
 448     if (rational) {
 449       ww = Weights(ii);
 450       if (ww < EpsW) {
 451         NullWeight = Standard_True;
 452       }
 453       else {
 454         D2Weights(ii) = Result(indice+3);
 455         D2Poles(ii)
 456           .SetXYZ( ( D2Poles(ii).XYZ() - D2Weights(ii)*Poles(ii).Coord()
 457                     - 2*DWeights(ii)*DPoles(ii).XYZ() ) / ww );
 458       }
 459     }
 460     indice += gap;
 461   }
 462   if (NullWeight) return Standard_False;
 463
 464   // verif par diff finies sous debug sauf pour les surfaces periodiques
 465 #ifdef OCCT_DEBUG
 466   if (!mySurface->IsVPeriodic()) {
 467     Standard_Real V1,V2;
 468     Standard_Boolean ok1,ok2;
 469     Standard_Real pas = 1.e-6, wTol = 1.e-4, pTol = 1.e-3;
 470     TColStd_Array1OfReal W1(1,L),W2(1,L),DW1(1,L),DW2(1,L);
 471     TColgp_Array1OfPnt P1(1,L),P2(1,L);
 472     TColgp_Array1OfVec DP1(1,L),DP2(1,L);
 473     gp_Pnt nul(0.,0.,0.);
 474     gp_Vec Vnul(0.,0.,0.);
 475     W1.Init(0.);
 476     W2.Init(0.);
 477     DW1.Init(0.);
 478     DW2.Init(0.);
 479     P1.Init(nul);
 480     P2.Init(nul);
 481     DP1.Init(Vnul);
 482     DP2.Init(Vnul);
 483
 484     V1 = V;
 485     V2 = V+pas;
 486     ok1 = D1(V1,P1,DP1,W1,DW1);
 487     ok2 = D1(V2,P2,DP2,W2,DW2);
 488     if (!ok1 || !ok2) std::cout<<"probleme en D0 ou en D1"<<std::endl;
 489     Standard_Boolean check = verifD2(DP1,DW1,DP2,DW2,D2Poles,D2Weights,pTol,wTol,pas);
 490     if (!check) std::cout<<"D2 incorrecte en V = "<<V<<std::endl;
 491   }
 492 #endif
 493
 494   return Standard_True;
 495 }
 496
 497 //=======================================================
 498 // Purpose :BSplineSurface()
 499 //=======================================================
 500  Handle(Geom_BSplineSurface)
 501      GeomFill_NSections::BSplineSurface() const
 502 {
 503   return mySurface;
 504 }
 505
 506
 507 //=======================================================
 508 // Purpose :SetSurface()
 509 //=======================================================
 510  void GeomFill_NSections::SetSurface(const Handle(Geom_BSplineSurface)& RefSurf)
 511 {
 512   myRefSurf = RefSurf;
 513 }
 514
 515 //=======================================================
 516 // Purpose :ComputeSurface()
 517 //=======================================================
 518  void GeomFill_NSections::ComputeSurface()
 519 {
 520
 521   Handle(Geom_BSplineSurface) BS;
 522   if (myRefSurf.IsNull()) {
 523
 524     Standard_Real myPres3d = 1.e-06;
 525     Standard_Integer i,j,jdeb=1,jfin=mySections.Length();
 526
 527     if (jfin <= jdeb)
 528     {
 529       //We will not be able to create surface from single curve.
 530       return;
 531     }
 532
 533     GeomFill_SectionGenerator section;
 534     Handle(Geom_BSplineSurface) surface;
 535
 536     for (j=jdeb; j<=jfin; j++) {
 537
 538         // read the j-th curve
 539         Handle(Geom_Curve) curv = mySections(j);
 540
 541         // transformation to BSpline reparametrized to [UFirst,ULast]
 542         Handle(Geom_BSplineCurve) curvBS = Handle(Geom_BSplineCurve)::DownCast (curv);
 543         if (curvBS.IsNull())
 544         {
 545           curvBS = GeomConvert::CurveToBSplineCurve (curv, Convert_QuasiAngular);
 546         }
 547
 548         TColStd_Array1OfReal BSK(1,curvBS->NbKnots());
 549         curvBS->Knots(BSK);
 550         BSplCLib::Reparametrize(UFirst,ULast,BSK);
 551         curvBS->SetKnots(BSK);
 552
 553         section.AddCurve(curvBS);
 554     }
 555
 556     /*
 557     if (s2Point) {
 558       curv =  mySections(jfin+1);
 559       first =  curv->FirstParameter();
 560       last = curv->LastParameter();
 561       TColgp_Array1OfPnt Extremities(1,2);
 562       Extremities(1) = curv->Value(first);
 563       Extremities(2) = curv->Value(last);
 564       TColStd_Array1OfReal Bounds(1,2);
 565       Bounds(1) = UFirst;
 566       Bounds(2) = ULast;
 567       Standard_Real Deg = 1;
 568       TColStd_Array1OfInteger Mult(1,2);
 569       Mult(1) = (Standard_Integer ) Deg+1;
 570       Mult(2) = (Standard_Integer ) Deg+1;
 571       Handle(Geom_BSplineCurve) BSPoint
 572         = new Geom_BSplineCurve(Extremities,Bounds,Mult,(Standard_Integer ) Deg);
 573       section.AddCurve(BSPoint);
 574     }*/
 575
 576     Standard_Integer Nbcurves = mySections.Length();
 577     Standard_Integer Nbpar = myParams.Length();
 578     if (Nbpar > 0)
 579     {
 580       Handle(TColStd_HArray1OfReal) HPar
 581         = new TColStd_HArray1OfReal(1, Nbpar);
 582       for (i = 1; i <= Nbpar; i++) {
 583         HPar->SetValue(i, myParams(i));
 584       }
 585       section.SetParam(HPar);
 586     }
 587     section.Perform(Precision::PConfusion());
 588
 589     Handle(GeomFill_Line) line = new GeomFill_Line(Nbcurves);
 590     Standard_Integer nbIt = 0, degmin = 2, degmax = 6;
 591     Standard_Boolean knownP = Nbpar > 0;
 592     GeomFill_AppSurf anApprox(degmin, degmax, myPres3d, myPres3d, nbIt, knownP);
 593     anApprox.SetContinuity(GeomAbs_C1);
 594     Standard_Boolean SpApprox = Standard_True;
 595     anApprox.Perform(line, section, SpApprox);
 596
 597     BS =
 598       new Geom_BSplineSurface(anApprox.SurfPoles(), anApprox.SurfWeights(),
 599                               anApprox.SurfUKnots(), anApprox.SurfVKnots(),
 600                               anApprox.SurfUMults(), anApprox.SurfVMults(),
 601                               anApprox.UDegree(), anApprox.VDegree());
 602   }
 603
 604   else {
 605
 606     // segmentation de myRefSurf
 607     Standard_Real Ui1, Ui2, V0, V1;
 608     BS = Handle(Geom_BSplineSurface)::DownCast(myRefSurf->Copy());
 609     Ui1 = UFirst;
 610     Ui2 = ULast;
 611     Standard_Integer i1, i2;
 612     myRefSurf->LocateU( Ui1, Precision::PConfusion(), i1, i2 );
 613     if (Abs(Ui1 - myRefSurf->UKnot(i1)) <= Precision::PConfusion())
 614       Ui1 = myRefSurf->UKnot(i1);
 615     if (Abs(Ui1 - myRefSurf->UKnot(i2)) <= Precision::PConfusion())
 616       Ui1 = myRefSurf->UKnot(i2);
 617     myRefSurf->LocateU( Ui2, Precision::PConfusion(), i1, i2 );
 618     if (Abs(Ui2 - myRefSurf->UKnot(i1)) <= Precision::PConfusion())
 619       Ui2 = myRefSurf->UKnot(i1);
 620     if (Abs(Ui2 - myRefSurf->UKnot(i2)) <= Precision::PConfusion())
 621       Ui2 = myRefSurf->UKnot(i2);
 622     V0  = myRefSurf->VKnot(myRefSurf->FirstVKnotIndex());
 623     V1  = myRefSurf->VKnot(myRefSurf->LastVKnotIndex());
 624     BS->CheckAndSegment(Ui1,Ui2,V0,V1);
 625   }
 626   mySurface = BS;
 627   // On augmente le degre pour que le positionnement D2 soit correct
 628   if (mySurface->VDegree()<2) {
 629     mySurface->IncreaseDegree(mySurface->UDegree(),2);
 630   }
 631 #ifdef OCCT_DEBUG
 632   NbSurf++;
 633   if (Affich) {
 634 #ifdef DRAW
 635     char name[256];
 636     sprintf(name,"NS_Surf_%d",NbSurf);
 637     DrawTrSurf::Set(name,BS);
 638     std::cout<<std::endl<<"RESULTAT de ComputeSurface : NS_Surf_"<<NbSurf<<std::endl<<std::endl;
 639 #endif
 640   }
 641 #endif
 642 }
 643
 644 //=======================================================
 645 // Purpose :SectionShape
 646 //=======================================================
 647  void GeomFill_NSections::SectionShape(Standard_Integer& NbPoles,
 648                                             Standard_Integer& NbKnots,
 649                                             Standard_Integer& Degree) const
 650 {
 651   if (mySurface.IsNull())
 652     return;
 653
 654   NbPoles = mySurface->NbUPoles();
 655   NbKnots = mySurface->NbUKnots();
 656   Degree  = mySurface->UDegree();
 657 }
 658
 659 //=======================================================
 660 // Purpose :Knots
 661 //=======================================================
 662  void GeomFill_NSections::Knots(TColStd_Array1OfReal& TKnots) const
 663 {
 664   if (!mySurface.IsNull())
 665     mySurface->UKnots(TKnots);
 666 }
 667
 668 //=======================================================
 669 // Purpose :Mults
 670 //=======================================================
 671  void GeomFill_NSections::Mults(TColStd_Array1OfInteger& TMults) const
 672 {
 673   if (!mySurface.IsNull())
 674     mySurface->UMultiplicities(TMults);
 675 }
 676
 677
 678 //=======================================================
 679 // Purpose :IsRational
 680 //=======================================================
 681  Standard_Boolean GeomFill_NSections::IsRational() const
 682 {
 683   if (!mySurface.IsNull())
 684     return mySurface->IsURational();
 685
 686   return Standard_False;
 687 }
 688
 689 //=======================================================
 690 // Purpose :IsUPeriodic
 691 //=======================================================
 692  Standard_Boolean GeomFill_NSections::IsUPeriodic() const
 693 {
 694   if (!mySurface.IsNull())
 695     return  mySurface->IsUPeriodic();
 696
 697   return Standard_False;
 698 }
 699
 700 //=======================================================
 701 // Purpose :IsVPeriodic
 702 //=======================================================
 703  Standard_Boolean GeomFill_NSections::IsVPeriodic() const
 704 {
 705   if (!mySurface.IsNull())
 706     return  mySurface->IsVPeriodic();
 707
 708   return Standard_False;
 709 }
 710
 711 //=======================================================
 712 // Purpose :NbIntervals
 713 //=======================================================
 714  Standard_Integer GeomFill_NSections::NbIntervals(const GeomAbs_Shape S) const
 715 {
 716   if (mySurface.IsNull())
 717     return 0;
 718
 719   GeomAdaptor_Surface AdS(mySurface);
 720   return AdS.NbVIntervals(S);
 721 }
 722
 723
 724 //=======================================================
 725 // Purpose :Intervals
 726 //=======================================================
 727  void GeomFill_NSections::Intervals(TColStd_Array1OfReal& T,
 728                                          const GeomAbs_Shape S) const
 729 {
 730   if (mySurface.IsNull())
 731     return;
 732
 733   GeomAdaptor_Surface AdS(mySurface);
 734   AdS.VIntervals(T,S);
 735 }
 736
 737
 738 //=======================================================
 739 // Purpose : SetInterval
 740 //=======================================================
 741 // void GeomFill_NSections::SetInterval(const Standard_Real F,
 742  void GeomFill_NSections::SetInterval(const Standard_Real ,
 743 //                                           const Standard_Real L)
 744                                            const Standard_Real )
 745 {
 746   // rien a faire : mySurface est supposee Cn en V
 747 }
 748
 749 //=======================================================
 750 // Purpose : GetInterval
 751 //=======================================================
 752  void GeomFill_NSections::GetInterval(Standard_Real& F,
 753                                            Standard_Real& L) const
 754 {
 755   F = VFirst;
 756   L = VLast;
 757 }
 758
 759 //=======================================================
 760 // Purpose : GetDomain
 761 //=======================================================
 762  void GeomFill_NSections::GetDomain(Standard_Real& F,
 763                                          Standard_Real& L) const
 764 {
 765   F = VFirst;
 766   L = VLast;
 767 }
 768
 769 //=======================================================
 770 // Purpose : GetTolerance
 771 //=======================================================
 772  void GeomFill_NSections::GetTolerance(const Standard_Real BoundTol,
 773                                             const Standard_Real SurfTol,
 774 //                                            const Standard_Real AngleTol,
 775                                             const Standard_Real ,
 776                                             TColStd_Array1OfReal& Tol3d) const
 777 {
 778   Tol3d.Init(SurfTol);
 779   if (BoundTol<SurfTol) {
 780     Tol3d(Tol3d.Lower()) = BoundTol;
 781     Tol3d(Tol3d.Upper()) = BoundTol;
 782   }
 783 }
 784
 785 //=======================================================
 786 // Purpose : BarycentreOfSurf
 787 //=======================================================
 788  gp_Pnt GeomFill_NSections::BarycentreOfSurf() const
 789 {
 790   gp_Pnt P, Bary;
 791   Bary.SetCoord(0., 0., 0.);
 792
 793   if (mySurface.IsNull())
 794     return Bary;
 795
 796   Standard_Integer ii,jj;
 797   Standard_Real U0, U1, V0, V1;
 798   mySurface->Bounds(U0,U1,V0,V1);
 799   Standard_Real V = V0, DeltaV = ( V1 - V0 ) / 20;
 800   Standard_Real U = U0, DeltaU = ( U1 - U0 ) / 20;
 801   for(jj=0;jj<=20;jj++,V+=DeltaV) {
 802     for (ii=0 ; ii <=20; ii++, U+=DeltaU) {
 803       P = mySurface->Value(U,V);
 804       Bary.ChangeCoord() += P.XYZ();
 805     }
 806   }
 807
 808   Bary.ChangeCoord() /= (21*21);
 809   return Bary;
 810 }
 811
 812
 813 //=======================================================
 814 // Purpose : MaximalSection
 815 //=======================================================
 816  Standard_Real GeomFill_NSections::MaximalSection() const
 817 {
 818   Standard_Real L, Lmax=0.;
 819   Standard_Integer ii;
 820   for (ii=1; ii <=mySections.Length(); ii++) {
 821     GeomAdaptor_Curve AC (mySections(ii));
 822     L = GCPnts_AbscissaPoint::Length(AC);
 823     if (L>Lmax) Lmax = L;
 824   }
 825   return Lmax;
 826 }
 827
 828
 829 //=======================================================
 830 // Purpose : GetMinimalWeight
 831 //=======================================================
 832 void GeomFill_NSections::GetMinimalWeight(TColStd_Array1OfReal& Weights) const
 833 {
 834   if (mySurface.IsNull())
 835     return;
 836
 837   if (mySurface->IsURational()) {
 838     Standard_Integer NbU = mySurface->NbUPoles(),
 839                      NbV = mySurface->NbVPoles();
 840     TColStd_Array2OfReal WSurf(1,NbU,1,NbV);
 841     mySurface->Weights(WSurf);
 842     Standard_Integer i,j;
 843     for (i=1;i<=NbU;i++) {
 844       Standard_Real min = WSurf(i,1);
 845       for (j=2;j<=NbV;j++) {
 846         if (min> WSurf(i,j)) min = WSurf(i,j);
 847       }
 848       Weights.SetValue(i,min);
 849     }
 850   }
 851   else {
 852     Weights.Init(1);
 853   }
 854
 855 }
 856
 857
 858 //=======================================================
 859 // Purpose : IsConstant
 860 //=======================================================
 861  Standard_Boolean GeomFill_NSections::IsConstant(Standard_Real& Error) const
 862 {
 863   // on se limite a 2 sections
 864   Standard_Boolean isconst = (mySections.Length()==2);
 865   Standard_Real Err = 0.;
 866
 867   if (isconst) {
 868     GeomAdaptor_Curve AC1(mySections(1));
 869     GeomAbs_CurveType CType = AC1.GetType();
 870     GeomAdaptor_Curve AC2(mySections(2));
 871     // les sections doivent avoir le meme type
 872     isconst = ( AC2.GetType() == CType);
 873
 874     if (isconst) {
 875       if (CType == GeomAbs_Circle) {
 876         gp_Circ C1 = AC1.Circle();
 877         gp_Circ C2 = AC2.Circle();
 878         Standard_Real Tol = 1.e-7;
 879         Standard_Boolean samedir, samerad, samepos;
 880         samedir = (C1.Axis().IsParallel(C2.Axis(),1.e-4));
 881         samerad = (Abs(C1.Radius()-C2.Radius())<Tol);
 882         samepos = (C1.Location().Distance(C2.Location())<Tol);
 883         if (!samepos) {
 884           gp_Ax1 D(C1.Location(),gp_Vec(C1.Location(),C2.Location()));
 885           samepos = (C1.Axis().IsParallel(D,1.e-4));
 886         }
 887         isconst = samedir && samerad && samepos;
 888       }
 889       else if (CType == GeomAbs_Line) {
 890         gp_Lin L1 = AC1.Line();
 891         gp_Lin L2 = AC2.Line();
 892         Standard_Real Tol = 1.e-7;
 893         Standard_Boolean samedir, samelength, samepos;
 894         samedir = (L1.Direction().IsParallel(L2.Direction(),1.e-4));
 895         gp_Pnt P11 = AC1.Value(AC1.FirstParameter()),
 896                P12 = AC1.Value(AC1.LastParameter()),
 897                P21 = AC2.Value(AC2.FirstParameter()),
 898                P22 = AC2.Value(AC2.LastParameter());
 899         samelength = (Abs(P11.Distance(P12)-P21.Distance(P22))<Tol);
 900         // l'ecart entre les 2 sections ne compte pas
 901         samepos = ( ( P11.Distance(P21)<Tol && P12.Distance(P22)<Tol )
 902                      || ( P12.Distance(P21)<Tol && P11.Distance(P22)<Tol ) );
 903         //samepos = Standard_True;
 904         isconst = samedir && samelength && samepos;
 905       }
 906       else {
 907         isconst = Standard_False;
 908       }
 909     }
 910
 911   }
 912
 913   Error = Err;
 914   return isconst;
 915 }
 916
 917
 918 //=======================================================
 919 // Purpose : ConstantSection
 920 //=======================================================
 921  Handle(Geom_Curve) GeomFill_NSections::ConstantSection() const
 922 {
 923 //  Standard_Real Err;
 924 //  if (!IsConstant(Err)) throw StdFail_NotDone("The Law is not Constant!");
 925   Handle(Geom_Curve) C;
 926   C = Handle(Geom_Curve)::DownCast( mySections(1)->Copy());
 927   return C;
 928 }
 929
 930
 931 //=======================================================
 932 // Purpose : IsConicalLaw
 933 //=======================================================
 934  Standard_Boolean GeomFill_NSections::IsConicalLaw(Standard_Real& Error) const
 935 {
 936   Standard_Boolean isconic = (mySections.Length()==2);
 937   Standard_Real Err = 0.;
 938   if (isconic) {
 939     GeomAdaptor_Curve AC1(mySections(1));
 940     GeomAdaptor_Curve AC2(mySections(2));
 941     isconic = ( AC1.GetType() == GeomAbs_Circle )
 942                      &&  ( AC2.GetType() == GeomAbs_Circle ) ;
 943     if (isconic) {
 944       gp_Circ C1 = AC1.Circle();
 945       if (!myTrsfs.IsEmpty())
 946         C1.Transform(myTrsfs(1).Inverted());
 947       gp_Circ C2 = AC2.Circle();
 948       if (!myTrsfs.IsEmpty())
 949         C2.Transform(myTrsfs(2).Inverted());
 950       Standard_Real Tol = 1.e-7;
 951       //Standard_Boolean samedir, linearrad, sameaxis;
 952       isconic = (C1.Axis().IsParallel(C2.Axis(),1.e-4));
 953       //  pour 2 sections, la variation du rayon est forcement lineaire
 954       //linearrad = Standard_True;
 955       //  formule plus generale pour 3 sections au moins
 956       //  Standard_Real param0 = C2.Radius()*myParams(1) - C1.Radius()*myParams(2);
 957       //  param0 = param0 / (C2.Radius()-C1.Radius()) ;
 958       //  linearrad = ( Abs( C3.Radius()*myParams(1)-C1.Radius()*myParams(3)
 959       //                          - param0*(C3.Radius()-C1.Radius()) ) < Tol);
 960       if (isconic)
 961       {
 962         gp_Lin Line1(C1.Axis());
 963         isconic = (Line1.Distance(C2.Location()) < Tol);
 964         /*
 965         sameaxis = (C1.Location().Distance(C2.Location())<Tol);
 966         if (!sameaxis) {
 967         gp_Ax1 D(C1.Location(),gp_Vec(C1.Location(),C2.Location()));
 968         sameaxis = (C1.Axis().IsParallel(D,1.e-4));
 969         }
 970         isconic = samedir && linearrad && sameaxis;
 971         */
 972         if (isconic)
 973         {
 974           //// Modified by jgv, 18.02.2009 for OCC20866 ////
 975           Standard_Real first1 = AC1.FirstParameter(), last1 = AC1.LastParameter();
 976           Standard_Real first2 = AC2.FirstParameter(), last2 = AC2.LastParameter();
 977           isconic = (Abs(first1-first2) <= Precision::PConfusion() &&
 978                      Abs(last1-last2)   <= Precision::PConfusion());
 979           //////////////////////////////////////////////////
 980         }
 981       }
 982     }
 983   }
 984
 985   Error = Err;
 986   return isconic;
 987 }
 988
 989
 990 //=======================================================
 991 // Purpose : CirclSection
 992 //=======================================================
 993  Handle(Geom_Curve) GeomFill_NSections::CirclSection(const Standard_Real V) const
 994 {
 995   Standard_Real Err;
 996   if (!IsConicalLaw(Err)) throw StdFail_NotDone("The Law is not Conical!");
 997
 998   GeomAdaptor_Curve AC1(mySections(1));
 999   GeomAdaptor_Curve AC2(mySections(mySections.Length()));
1000   gp_Circ C1 = AC1.Circle();
1001   gp_Circ C2 = AC2.Circle();
1002
1003   Standard_Real p1 = myParams(1), p2 = myParams(myParams.Length());
1004   Standard_Real radius = ( C2.Radius() - C1.Radius() ) * (V - p1) / (p2 - p1)
1005                                   + C1.Radius();
1006
1007   C1.SetRadius(radius);
1008   Handle(Geom_Curve) C = new (Geom_Circle) (C1);
1009
1010   const Standard_Real aParF = AC1.FirstParameter();
1011   const Standard_Real aParL = AC1.LastParameter();
1012   const Standard_Real aPeriod = AC1.IsPeriodic() ? AC1.Period() : 0.0;
1013
1014   if ((aPeriod == 0.0) || (Abs(aParL - aParF - aPeriod) > Precision::PConfusion()))
1015   {
1016     Handle(Geom_Curve) Cbis = new Geom_TrimmedCurve(C, aParF, aParL);
1017     C = Cbis;
1018   }
1019   return C;
1020 }