src/AppBlend/AppBlend_AppSurf.gxx

   1 // Copyright (c) 1995-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 #include <AppDef_MultiLine.hxx>
  16 #include <AppDef_MultiPointConstraint.hxx>
  17 #include <AppParCurves_MultiBSpCurve.hxx>
  18 #include <AppParCurves_MultiCurve.hxx>
  19 #include <AppDef_BSplineCompute.hxx>
  20 #include <AppDef_Compute.hxx>
  21 #include <AppParCurves_Constraint.hxx>
  22 #include <Approx_MCurvesToBSpCurve.hxx>
  23 #include <TColgp_Array1OfPnt.hxx>
  24 #include <TColgp_Array1OfPnt2d.hxx>
  25 #include <TColgp_Array1OfVec.hxx>
  26 #include <TColgp_Array1OfVec2d.hxx>
  27 #include <gp_Vec.hxx>
  28 #include <gp_Vec2d.hxx>
  29 #include <gp_Pnt.hxx>
  30 #include <gp_Pnt2d.hxx>
  31 #include <math_Vector.hxx>
  32 #include <BSplCLib.hxx>
  33
  34 #include <StdFail_NotDone.hxx>
  35 #include <AppParCurves_HArray1OfConstraintCouple.hxx>
  36 #include <AppDef_Variational.hxx>
  37
  38 static   Standard_Boolean scal = 1;
  39 #ifdef DEB
  40 extern Standard_Boolean AppBlend_GetContextSplineApprox();
  41 extern Standard_Boolean AppBlend_GetContextApproxWithNoTgt();
  42 #endif
  43
  44 //  modified by EAP (Edward AGAPOV) Fri Jan 4 2002, bug OCC9
  45 //  --- keep pipe parametrized like path
  46
  47
  48 //=======================================================================
  49 //function : AppBlend_AppSurf
  50 //purpose  :
  51 //=======================================================================
  52
  53 AppBlend_AppSurf::AppBlend_AppSurf ():done(Standard_False) {}
  54
  55
  56 //=======================================================================
  57 //function : AppBlend_AppSurf
  58 //purpose  :
  59 //=======================================================================
  60
  61 AppBlend_AppSurf::AppBlend_AppSurf (const Standard_Integer Degmin,
  62                                     const Standard_Integer Degmax,
  63                                     const Standard_Real Tol3d,
  64                                     const Standard_Real Tol2d,
  65                                     const Standard_Integer NbIt,
  66                                     const Standard_Boolean KnownParameters):
  67        done(Standard_False),dmin(Degmin),dmax(Degmax),
  68        tol3d(Tol3d),tol2d(Tol2d),nbit(NbIt),knownp(KnownParameters)
  69 {
  70   continuity = GeomAbs_C2;
  71   paramtype = Approx_ChordLength;
  72   critweights[0]=0.4;
  73   critweights[1]=0.2;
  74   critweights[2]=0.4;
  75 }
  76
  77 //=======================================================================
  78 //function : Init
  79 //purpose  :
  80 //=======================================================================
  81
  82 void AppBlend_AppSurf::Init (const Standard_Integer Degmin,
  83                              const Standard_Integer Degmax,
  84                              const Standard_Real Tol3d,
  85                              const Standard_Real Tol2d,
  86                              const Standard_Integer NbIt,
  87                              const Standard_Boolean KnownParameters)
  88 {
  89   done  = Standard_False;
  90   dmin  = Degmin;
  91   dmax  = Degmax;
  92   tol3d = Tol3d;
  93   tol2d = Tol2d;
  94   nbit  = NbIt;
  95   knownp = KnownParameters;
  96   continuity = GeomAbs_C2;
  97   paramtype = Approx_ChordLength;
  98   critweights[0]=0.4;
  99   critweights[1]=0.2;
 100   critweights[2]=0.4;
 101 }
 102
 103 //=======================================================================
 104 //function : CriteriumWeight
 105 //purpose  : returns the Weights associed  to the criterium used in
 106 //           the  optimization.
 107 //=======================================================================
 108 //
 109 void AppBlend_AppSurf::CriteriumWeight(Standard_Real& W1, Standard_Real& W2, Standard_Real& W3) const
 110 {
 111   W1 = critweights[0];
 112   W2 = critweights[1];
 113   W3 = critweights[2] ;
 114 }
 115 //=======================================================================
 116 //function : SetCriteriumWeight
 117 //purpose  :
 118 //=======================================================================
 119
 120 void AppBlend_AppSurf::SetCriteriumWeight(const Standard_Real W1, const Standard_Real W2, const Standard_Real W3)
 121 {
 122   if (W1 < 0 || W2 < 0 || W3 < 0 ) Standard_DomainError::Raise();
 123   critweights[0] = W1;
 124   critweights[1] = W2;
 125   critweights[2] = W3;
 126 }
 127 //=======================================================================
 128 //function : SetContinuity
 129 //purpose  :
 130 //=======================================================================
 131
 132 void AppBlend_AppSurf::SetContinuity (const GeomAbs_Shape TheCont)
 133 {
 134   continuity = TheCont;
 135 }
 136
 137 //=======================================================================
 138 //function : Continuity
 139 //purpose  :
 140 //=======================================================================
 141
 142 GeomAbs_Shape AppBlend_AppSurf::Continuity () const
 143 {
 144   return continuity;
 145 }
 146
 147 //=======================================================================
 148 //function : SetParType
 149 //purpose  :
 150 //=======================================================================
 151
 152 void AppBlend_AppSurf::SetParType (const Approx_ParametrizationType ParType)
 153 {
 154   paramtype = ParType;
 155 }
 156
 157 //=======================================================================
 158 //function : ParType
 159 //purpose  :
 160 //=======================================================================
 161
 162 Approx_ParametrizationType AppBlend_AppSurf::ParType () const
 163 {
 164   return paramtype;
 165 }
 166
 167
 168 //=======================================================================
 169 //function : Perform
 170 //purpose  :
 171 //=======================================================================
 172
 173 void AppBlend_AppSurf::Perform(const Handle(TheLine)& Lin,
 174                                TheSectionGenerator& F,
 175                                const Standard_Boolean SpApprox)
 176
 177 {
 178   InternalPerform(Lin, F, SpApprox, Standard_False);
 179 }
 180
 181 //=======================================================================
 182 //function : PerformSmoothing
 183 //purpose  :
 184 //=======================================================================
 185
 186 void AppBlend_AppSurf::PerformSmoothing(const Handle(TheLine)& Lin,
 187                                           TheSectionGenerator& F)
 188
 189 {
 190   InternalPerform(Lin, F, Standard_True, Standard_True);
 191 }
 192
 193 //=======================================================================
 194 //function : InternalPerform
 195 //purpose  :
 196 //=======================================================================
 197
 198 void AppBlend_AppSurf::InternalPerform(const Handle(TheLine)& Lin,
 199                                        TheSectionGenerator& F,
 200                                        const Standard_Boolean SpApprox,
 201                                        const Standard_Boolean UseSmoothing)
 202
 203 {
 204   done = Standard_False;
 205   if (Lin.IsNull()) {return;}
 206   Standard_Integer i,j,k,NbPoint;
 207   Standard_Integer NbUPoles,NbUKnots,NbPoles2d,NbVPoles;
 208   Standard_Boolean withderiv;
 209   AppParCurves_Constraint Cfirst,Clast;
 210
 211   Standard_Real mytol3d,mytol2d;
 212   gp_XYZ newDv;
 213
 214   seqPoles2d.Clear();
 215
 216   NbPoint=Lin->NbPoints();
 217   AppDef_MultiPointConstraint multP;
 218   AppDef_MultiLine multL(NbPoint);
 219
 220   F.GetShape(NbUPoles,NbUKnots,udeg,NbPoles2d);
 221
 222   tabUKnots  = new TColStd_HArray1OfReal (1,NbUKnots);
 223   tabUMults  = new TColStd_HArray1OfInteger (1,NbUKnots);
 224
 225   F.Knots(tabUKnots->ChangeArray1());
 226   F.Mults(tabUMults->ChangeArray1());
 227
 228   TColgp_Array1OfPnt tabAppP(1,NbUPoles);
 229   TColgp_Array1OfVec tabAppV(1,NbUPoles);
 230
 231   TColgp_Array1OfPnt2d tabP2d(1,Max(1,NbPoles2d));
 232   TColgp_Array1OfVec2d tabV2d(1,Max(1,NbPoles2d));
 233
 234   TColStd_Array1OfReal tabW(1,NbUPoles),tabDW(1,NbUPoles);
 235
 236   TColgp_Array1OfPnt2d tabAppP2d(1,NbPoles2d+NbUPoles); // points2d + poids
 237   TColgp_Array1OfVec2d tabAppV2d(1,NbPoles2d+NbUPoles);
 238
 239
 240   AppParCurves_MultiBSpCurve multC;
 241
 242 //  Standard_Boolean SpApprox = Standard_False;
 243
 244   withderiv = F.Section(Lin->Point(1),tabAppP,tabAppV,tabP2d,tabV2d,
 245                         tabW,tabDW);
 246
 247 #if DEB
 248 //  SpApprox = AppBlend_GetContextSplineApprox();
 249   if(AppBlend_GetContextApproxWithNoTgt()) withderiv = Standard_False;
 250 #endif
 251
 252
 253   for (j=1; j<=NbPoles2d; j++) {
 254     tabAppP2d(j) = tabP2d(j);
 255     if (withderiv) {
 256       tabAppV2d(j) = tabV2d(j);
 257     }
 258   }
 259   for (j=1; j<=NbUPoles; j++) {
 260     // pour les courbes rationnelles il faut multiplier les poles par
 261     // leurs poids respectifs
 262     if (withderiv) {
 263       tabAppV2d(NbPoles2d+j).SetCoord(tabDW(j),0.);
 264       newDv.SetLinearForm(tabDW(j),tabAppP(j).XYZ(),tabW(j),tabAppV(j).XYZ());
 265       tabAppV(j).SetXYZ(newDv);
 266     }
 267     tabAppP(j).SetXYZ(tabAppP(j).XYZ() * tabW(j));
 268     tabAppP2d(NbPoles2d+j).SetCoord(tabW(j),0.);
 269   }
 270
 271   if (withderiv) {
 272     multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d,tabAppV,tabAppV2d);
 273     Cfirst = AppParCurves_TangencyPoint;
 274   }
 275   else {
 276     multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d);
 277     Cfirst = AppParCurves_PassPoint;
 278   }
 279   multL.SetValue(1,multP);
 280
 281   for (i=2; i<=NbPoint-1; i++) {
 282     if (SpApprox) {
 283       F.Section(Lin->Point(i),tabAppP,tabP2d,tabW);
 284       for (j=1; j<=NbPoles2d; j++) {
 285         tabAppP2d(j) = tabP2d(j);
 286       }
 287       for (j=1; j<=NbUPoles; j++) {
 288         // pour les courbes rationnelles il faut multiplier les poles par
 289         // leurs poids respectifs
 290         tabAppP(j).SetXYZ(tabAppP(j).XYZ() * tabW(j));
 291         tabAppP2d(NbPoles2d+j).SetCoord(tabW(j),0.);
 292       }
 293       multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d);
 294       multL.SetValue(i,multP);
 295     }
 296 // ***********************
 297     else {
 298       withderiv = F.Section(Lin->Point(i),tabAppP,tabAppV,tabP2d,tabV2d,
 299                             tabW,tabDW);
 300 #if DEB
 301       if(AppBlend_GetContextApproxWithNoTgt()) withderiv = Standard_False;
 302 #endif
 303
 304       for (j=1; j<=NbPoles2d; j++) {
 305         tabAppP2d(j) = tabP2d(j);
 306         if (withderiv) {
 307           tabAppV2d(j) = tabV2d(j);
 308         }
 309       }
 310       for (j=1; j<=NbUPoles; j++) {
 311         // pour les courbes rationnelles il faut multiplier les poles par
 312         // leurs poids respectifs
 313         if (withderiv) {
 314           tabAppV2d(NbPoles2d+j).SetCoord(tabDW(j),0.);
 315           newDv.SetLinearForm(tabDW(j),tabAppP(j).XYZ(),tabW(j),tabAppV(j).XYZ());
 316           tabAppV(j).SetXYZ(newDv);
 317         }
 318         tabAppP(j).SetXYZ(tabAppP(j).XYZ() * tabW(j));
 319         tabAppP2d(NbPoles2d+j).SetCoord(tabW(j),0.);
 320       }
 321       if (withderiv) {
 322         multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d,tabAppV,tabAppV2d);
 323       }
 324       else {
 325         multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d);
 326       }
 327       multL.SetValue(i,multP);
 328     }
 329 // ******************************
 330   }
 331
 332   withderiv = F.Section(Lin->Point(NbPoint),tabAppP,tabAppV,tabP2d,tabV2d,
 333                         tabW,tabDW);
 334 #if DEB
 335       if(AppBlend_GetContextApproxWithNoTgt()) withderiv = Standard_False;
 336 #endif
 337
 338   for (j=1; j<=NbPoles2d; j++) {
 339     tabAppP2d(j) = tabP2d(j);
 340     if (withderiv) {
 341       tabAppV2d(j) = tabV2d(j);
 342     }
 343   }
 344   for (j=1; j<=NbUPoles; j++) {
 345     // pour les courbes rationnelles il faut multiplier les poles par
 346     // leurs poids respectifs
 347     if (withderiv) {
 348       tabAppV2d(NbPoles2d+j).SetCoord(tabDW(j),0.);
 349       newDv.SetLinearForm(tabDW(j),tabAppP(j).XYZ(),tabW(j),tabAppV(j).XYZ());
 350       tabAppV(j).SetXYZ(newDv);
 351     }
 352     tabAppP(j).SetXYZ(tabAppP(j).XYZ() * tabW(j));
 353     tabAppP2d(NbPoles2d+j).SetCoord(tabW(j),0.);
 354   }
 355
 356   if (withderiv) {
 357     multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d,tabAppV,tabAppV2d);
 358     Clast = AppParCurves_TangencyPoint;
 359   }
 360   else {
 361     multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d);
 362     Clast = AppParCurves_PassPoint;
 363   }
 364   multL.SetValue(NbPoint,multP);
 365
 366   //IFV 04.06.07 occ13904
 367   if(NbPoint == 2) {
 368     dmin = 1;
 369     if(Cfirst == AppParCurves_PassPoint && Clast == AppParCurves_PassPoint) {
 370       dmax = 1;
 371     }
 372   }
 373
 374
 375   if (!SpApprox) {
 376     AppDef_Compute theapprox (dmin,dmax,tol3d,tol2d,nbit, Standard_True, paramtype);
 377     if (knownp) {
 378       math_Vector theParams(1,NbPoint);
 379
 380       // On recale les parametres entre 0 et 1.
 381       theParams(1) = 0.;
 382       theParams(NbPoint) = 1.;
 383       Standard_Real Uf = F.Parameter(Lin->Point(1));
 384       Standard_Real Ul = F.Parameter(Lin->Point(NbPoint))-Uf;
 385       for (i=2; i<NbPoint; i++) {
 386         theParams(i) = (F.Parameter(Lin->Point(i))-Uf)/Ul;;
 387       }
 388       AppDef_Compute theAppDef(theParams,dmin,dmax,tol3d,tol2d,nbit,
 389                                  Standard_True, Standard_True);
 390       theapprox = theAppDef;
 391     }
 392     theapprox.SetConstraints(Cfirst,Clast);
 393     theapprox.Perform(multL);
 394
 395     Standard_Real TheTol3d, TheTol2d;
 396     mytol3d = mytol2d = 0.0;
 397     for (Standard_Integer Index=1; Index<=theapprox.NbMultiCurves(); Index++) {
 398       theapprox.Error(Index, TheTol3d, TheTol2d);
 399       mytol3d = Max(TheTol3d, mytol3d);
 400       mytol2d = Max(TheTol2d, mytol2d);
 401     }
 402 #ifdef DEB
 403     cout << " Tolerances obtenues  --> 3d : "<< mytol3d << endl;
 404     cout << "                      --> 2d : "<< mytol2d << endl;
 405 #endif
 406     multC = theapprox.SplineValue();
 407   }
 408
 409   else {
 410     if(!UseSmoothing) {
 411       Standard_Boolean UseSquares = Standard_False;
 412       if(nbit == 0) UseSquares = Standard_True;
 413       AppDef_BSplineCompute theapprox (dmin,dmax,tol3d,tol2d,nbit,Standard_True, paramtype,
 414                                        UseSquares);
 415       if(continuity == GeomAbs_C0) {
 416         theapprox.SetContinuity(0);
 417       }
 418       if(continuity == GeomAbs_C1) {
 419         theapprox.SetContinuity(1);
 420       }
 421       else if(continuity == GeomAbs_C2) {
 422         theapprox.SetContinuity(2);
 423       }
 424       else {
 425         theapprox.SetContinuity(3);
 426       }
 427
 428       theapprox.SetConstraints(Cfirst,Clast);
 429
 430       if (knownp) {
 431         math_Vector theParams(1,NbPoint);
 432         // On recale les parametres entre 0 et 1.
 433         theParams(1) = 0.;
 434         theParams(NbPoint) = 1.;
 435         Standard_Real Uf = F.Parameter(Lin->Point(1));
 436         Standard_Real Ul = F.Parameter(Lin->Point(NbPoint))-Uf;
 437         for (i=2; i<NbPoint; i++) {
 438           theParams(i) = (F.Parameter(Lin->Point(i))-Uf)/Ul;;
 439         }
 440
 441         theapprox.Init(dmin,dmax,tol3d,tol2d,nbit,Standard_True,
 442                        Approx_IsoParametric,Standard_True);
 443         theapprox.SetParameters(theParams);
 444       }
 445       theapprox.Perform(multL);
 446       theapprox.Error(mytol3d,mytol2d);
 447 #ifdef DEB
 448       cout << " Tolerances obtenues  --> 3d : "<< mytol3d << endl;
 449       cout << "                      --> 2d : "<< mytol2d << endl;
 450 #endif
 451       tol3dreached = mytol3d;
 452       tol2dreached = mytol2d;
 453       multC = theapprox.Value();
 454     }
 455     else {
 456       //Variational algo
 457       Handle(AppParCurves_HArray1OfConstraintCouple) TABofCC =
 458         new AppParCurves_HArray1OfConstraintCouple(1, NbPoint);
 459       AppParCurves_Constraint  Constraint=AppParCurves_NoConstraint;
 460
 461       for(i = 1; i <= NbPoint; ++i) {
 462         AppParCurves_ConstraintCouple ACC(i,Constraint);
 463         TABofCC->SetValue(i,ACC);
 464       }
 465
 466       TABofCC->ChangeValue(1).SetConstraint(Cfirst);
 467       TABofCC->ChangeValue(NbPoint).SetConstraint(Clast);
 468
 469       AppDef_Variational Variation(multL, 1, NbPoint, TABofCC);
 470
 471 //===================================
 472       Standard_Integer theMaxSegments = 1000;
 473       Standard_Boolean theWithMinMax = Standard_False;
 474       Standard_Boolean theWithCutting = Standard_True;
 475 //===================================
 476
 477       Variation.SetMaxDegree(dmax);
 478       Variation.SetContinuity(continuity);
 479       Variation.SetMaxSegment(theMaxSegments);
 480
 481       Variation.SetTolerance(tol3d);
 482       Variation.SetWithMinMax(theWithMinMax);
 483       Variation.SetWithCutting(theWithCutting);
 484       Variation.SetNbIterations(nbit);
 485
 486       Variation.SetCriteriumWeight(critweights[0], critweights[1], critweights[2]);
 487
 488       if(!Variation.IsCreated()) {
 489         return;
 490       }
 491
 492       if(Variation.IsOverConstrained()) {
 493         return;
 494       }
 495
 496       try {
 497         Variation.Approximate();
 498       }
 499       catch (Standard_Failure) {
 500         return;
 501       }
 502
 503       if(!Variation.IsDone()) {
 504         return;
 505       }
 506
 507       mytol3d = Variation.MaxError();
 508       mytol2d = 0.;
 509 #ifdef DEB
 510       cout << " Tolerances obtenues  --> 3d : "<< mytol3d << endl;
 511       cout << "                      --> 2d : "<< mytol2d << endl;
 512 #endif
 513       tol3dreached = mytol3d;
 514       tol2dreached = mytol2d;
 515       multC = Variation.Value();
 516     }
 517   }
 518
 519   vdeg = multC.Degree();
 520   NbVPoles = multC.NbPoles();
 521
 522   tabPoles   = new TColgp_HArray2OfPnt (1,NbUPoles,1,NbVPoles);
 523   tabWeights = new TColStd_HArray2OfReal (1,NbUPoles,1,NbVPoles);
 524   tabVKnots  = new TColStd_HArray1OfReal (multC.Knots().Lower(),
 525                                           multC.Knots().Upper());
 526   tabVKnots->ChangeArray1() = multC.Knots();
 527
 528   if (knownp && !UseSmoothing) {
 529     BSplCLib::Reparametrize(F.Parameter(Lin->Point(1)),
 530                             F.Parameter(Lin->Point(NbPoint)),
 531                             tabVKnots->ChangeArray1());
 532   }
 533
 534   tabVMults  = new TColStd_HArray1OfInteger (multC.Multiplicities().Lower(),
 535                                              multC.Multiplicities().Upper());
 536   tabVMults->ChangeArray1() = multC.Multiplicities();
 537
 538
 539   TColgp_Array1OfPnt newtabP(1,NbVPoles);
 540   Handle(TColgp_HArray1OfPnt2d) newtabP2d =
 541     new TColgp_HArray1OfPnt2d(1,NbVPoles);
 542   for (j=1; j <=NbUPoles; j++) {
 543     multC.Curve(j,newtabP);
 544     multC.Curve(j+NbUPoles+NbPoles2d,newtabP2d->ChangeArray1());
 545     for (k=1; k<=NbVPoles; k++) {
 546       // pour les courbes rationnelles il faut maintenant diviser
 547       // les poles par leurs poids respectifs
 548       tabPoles->ChangeValue(j,k).
 549         SetXYZ(newtabP(k).XYZ()/newtabP2d->Value(k).X());
 550       tabWeights->SetValue(j,k,newtabP2d->Value(k).X());
 551     }
 552   }
 553
 554   for (j=1; j<=NbPoles2d; j++) {
 555     newtabP2d = new TColgp_HArray1OfPnt2d(1,NbVPoles);
 556     multC.Curve(NbUPoles+j,newtabP2d->ChangeArray1());
 557     seqPoles2d.Append(newtabP2d);
 558   }
 559
 560   done = Standard_True;
 561 }
 562
 563
 564 //=======================================================================
 565 //function : Perform
 566 //purpose  :
 567 //=======================================================================
 568
 569 void AppBlend_AppSurf::Perform(const Handle(TheLine)& Lin,
 570                                TheSectionGenerator& F,
 571                                const Standard_Integer NbMaxP)
 572 {
 573   done = Standard_False;
 574   if (Lin.IsNull()) {return;}
 575   Standard_Integer i,j,k;
 576   Standard_Integer NbUPoles,NbUKnots,NbPoles2d,NbVPoles;
 577   Standard_Boolean withderiv;
 578   AppParCurves_Constraint Cfirst=AppParCurves_NoConstraint,Clast=AppParCurves_NoConstraint;
 579
 580   Standard_Real mytol3d = 0.0, mytol2d = 0.0;
 581   gp_XYZ newDv;
 582
 583   seqPoles2d.Clear();
 584
 585   Standard_Integer NbPointTot = Lin->NbPoints();
 586
 587   F.GetShape(NbUPoles,NbUKnots,udeg,NbPoles2d);
 588
 589   tabUKnots  = new TColStd_HArray1OfReal (1,NbUKnots);
 590   tabUMults  = new TColStd_HArray1OfInteger (1,NbUKnots);
 591
 592   F.Knots(tabUKnots->ChangeArray1());
 593   F.Mults(tabUMults->ChangeArray1());
 594
 595   TColgp_Array1OfPnt tabAppP(1,NbUPoles);
 596   TColgp_Array1OfVec tabAppV(1,NbUPoles);
 597   Standard_Real X,Y,Z,DX,DY,DZ;
 598   X = Y = Z = RealLast();
 599   DX = DY = DZ = RealFirst();
 600
 601   TColgp_Array1OfPnt2d tabP2d(1,Max(1,NbPoles2d));
 602   TColgp_Array1OfVec2d tabV2d(1,Max(1,NbPoles2d));
 603   TColStd_Array1OfReal X2d(1,Max(1,NbPoles2d));X2d.Init(RealLast());
 604   TColStd_Array1OfReal Y2d(1,Max(1,NbPoles2d));Y2d.Init(RealLast());
 605   TColStd_Array1OfReal DX2d(1,Max(1,NbPoles2d));DX2d.Init(RealFirst());
 606   TColStd_Array1OfReal DY2d(1,Max(1,NbPoles2d));DY2d.Init(RealFirst());
 607
 608   TColStd_Array1OfReal tabW(1,NbUPoles),tabDW(1,NbUPoles);
 609
 610   TColgp_Array1OfPnt2d tabAppP2d(1,NbPoles2d+NbUPoles); // points2d + poids
 611   TColgp_Array1OfVec2d tabAppV2d(1,NbPoles2d+NbUPoles);
 612
 613   // On calcule les boites de chaque ligne (box for all lines)
 614   for(i = 1; i <= NbPointTot; i++){
 615     F.Section(Lin->Point(i),tabAppP,tabAppV,tabP2d,tabV2d,tabW,tabDW);
 616     Standard_Real x,y,z;
 617     for(j = 1; j <= NbUPoles; j++){
 618       tabAppP(j).Coord(x,y,z);
 619       if(x < X) X = x; if(x > DX) DX = x;
 620       if(y < Y) Y = y; if(y > DY) DY = y;
 621       if(z < Z) Z = z; if(z > DZ) DZ = z;
 622     }
 623     for(j = 1; j <= NbPoles2d; j++){
 624       tabP2d(j).Coord(x,y);
 625       if(x < X2d(j)) X2d(j) = x; if(x > DX2d(j)) DX2d(j) = x;
 626       if(y < Y2d(j)) Y2d(j) = y; if(y > DY2d(j)) DY2d(j) = y;
 627     }
 628   }
 629   // On calcule pour chaque ligne la transformation vers 0 1.
 630   Standard_Real seuil = 1000.*tol3d;
 631   Standard_Real seuil2d = 1000.*tol2d;
 632   if((DX - X) < seuil ){ DX = 1.; X = 0.; }
 633   else{ DX = 1./(DX - X); X *= -DX; }
 634   if((DY - Y) < seuil){ DY = 1.; Y = 0.; }
 635   else{ DY = 1./(DY - Y); Y *= -DY; }
 636   if((DZ - Z) < seuil){ DZ = 1.; Z = 0.; }
 637   else{ DZ = 1./(DZ - Z); Z *= -DZ; }
 638   for(j = 1; j <= NbPoles2d; j++){
 639     if((DX2d(j) - X2d(j)) < seuil2d){ DX2d(j) = 1.; X2d(j) = 0.; }
 640     else{ DX2d(j) = 1./(DX2d(j) - X2d(j)); X2d(j) *= -DX2d(j); }
 641     if((DY2d(j) - Y2d(j)) < seuil2d){ DY2d(j) = 1.; Y2d(j) = 0.; }
 642     else{ DY2d(j) = 1./(DY2d(j) - Y2d(j)); Y2d(j) *= -DY2d(j); }
 643   }
 644   if(!scal){
 645     DX = 1.; X = 0.;
 646     DY = 1.; Y = 0.;
 647     DZ = 1.; Z = 0.;
 648     for(j = 1; j <= NbPoles2d; j++){
 649       DX2d(j) = 1.; X2d(j) = 0.;
 650       DY2d(j) = 1.; Y2d(j) = 0.;
 651     }
 652   }
 653 //  modified by eap Thu Jan  3 14:45:22 2002 ___BEGIN___
 654   // Keep "inter-troncons" parameters, not only first and last
 655 //  Standard_Real Ufirst=0,Ulast=0;
 656   TColStd_SequenceOfReal aParamSeq;
 657    if (knownp) {
 658 //     Ufirst = F.Parameter(Lin->Point(1));
 659 //     Ulast = F.Parameter(Lin->Point(NbPointTot));
 660      aParamSeq.Append( F.Parameter (Lin->Point(1)) );
 661   }
 662 //  modified by EAP Thu Jan  3 14:45:41 2002 ___END___
 663
 664   Approx_MCurvesToBSpCurve concat;
 665
 666   //On calcule le nombre de troncons.
 667   Standard_Integer nbtronc = NbPointTot/NbMaxP;
 668   Standard_Integer reste = NbPointTot - (nbtronc * NbMaxP);
 669   // On regarde si il faut prendre un troncon de plus.
 670   Standard_Integer nmax = NbMaxP;
 671   if(nbtronc > 0 && reste > 0){
 672     nmax = NbPointTot/(nbtronc + 1);
 673     if(nmax > (2*NbMaxP)/3) {
 674       nbtronc++;
 675       reste = NbPointTot - (nbtronc * nmax);
 676     }
 677     else nmax = NbMaxP;
 678   }
 679   else if(nbtronc == 0){
 680     nbtronc = 1;
 681     nmax = reste;
 682     reste = 0;
 683   }
 684
 685   // Approximate each "troncon" with nb of Bezier's using AppDef_Compute
 686   // and concat them into BSpline with Approx_MCurvesToBSpCurve
 687
 688   TColStd_Array1OfInteger troncsize(1,nbtronc);
 689   TColStd_Array1OfInteger troncstart(1,nbtronc);
 690
 691   Standard_Integer rab = reste/nbtronc + 1;
 692   Standard_Integer start = 1;
 693   Standard_Integer itronc ;
 694   for( itronc = 1; itronc <= nbtronc; itronc++){
 695     troncstart(itronc) = start;
 696     Standard_Integer rabrab = Min(rab,reste);
 697     if(reste > 0){ reste -= rabrab; }
 698     troncsize(itronc) = nmax + rabrab + 1;
 699     start += (nmax + rabrab);
 700   }
 701   troncsize(nbtronc) = troncsize(nbtronc) - 1;
 702   for(itronc = 1; itronc <= nbtronc; itronc++){
 703     Standard_Integer NbPoint = troncsize(itronc);
 704     Standard_Integer StPoint = troncstart(itronc);
 705     AppDef_MultiPointConstraint multP;
 706     AppDef_MultiLine multL(NbPoint);
 707
 708     for (i=1; i<=NbPoint; i++) {
 709       Standard_Integer iLin = StPoint + i - 1;
 710       Standard_Real x,y,z;
 711       withderiv = F.Section(Lin->Point(iLin),tabAppP,tabAppV,tabP2d,tabV2d,
 712                             tabW,tabDW);
 713 #if DEB
 714       if(AppBlend_GetContextApproxWithNoTgt()) withderiv = Standard_False;
 715 #endif
 716
 717       for (j=1; j<=NbPoles2d; j++) {
 718         tabP2d(j).Coord(x,y);
 719         tabAppP2d(j).SetCoord(DX2d(j)*x+X2d(j),DY2d(j)*y+Y2d(j));
 720         if (withderiv) {
 721           tabV2d(j).Coord(x,y);
 722           tabAppV2d(j).SetCoord(DX2d(j)*x,DY2d(j)*y);
 723         }
 724       }
 725       for (j=1; j<=NbUPoles; j++) {
 726         // pour les courbes rationnelles il faut multiplier les poles par
 727         // leurs poids respectifs
 728         if (withderiv) {
 729           tabAppV2d(NbPoles2d+j).SetCoord(tabDW(j),0.);
 730           newDv.SetLinearForm(tabDW(j),tabAppP(j).XYZ(),tabW(j),tabAppV(j).XYZ());
 731           tabAppV(j).SetCoord(DX*newDv.X(),DY*newDv.Y(),DZ*newDv.Z());
 732         }
 733         tabAppP(j).SetXYZ(tabAppP(j).XYZ() * tabW(j));
 734         tabAppP2d(NbPoles2d+j).SetCoord(tabW(j),0.);
 735         tabAppP(j).Coord(x,y,z);
 736         tabAppP(j).SetCoord(DX*x+X,DY*y+Y,DZ*z+Z);
 737       }
 738       if (withderiv) {
 739         multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d,tabAppV,tabAppV2d);
 740         if(i == 1) Cfirst = AppParCurves_TangencyPoint;
 741         else if(i == NbPoint) Clast = AppParCurves_TangencyPoint;
 742       }
 743       else {
 744         multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d);
 745         if(i == 1) Cfirst = AppParCurves_PassPoint;
 746         else if(i == NbPoint) Clast = AppParCurves_PassPoint;
 747       }
 748       multL.SetValue(i,multP);
 749     }
 750
 751
 752   //IFV 04.06.07 occ13904
 753     if(NbPoint == 2) {
 754       dmin = 1;
 755       if(Cfirst == AppParCurves_PassPoint && Clast == AppParCurves_PassPoint) {
 756         dmax = 1;
 757       }
 758     }
 759
 760 //  modified by EAP Thu Jan  3 15:44:13 2002 ___BEGIN___
 761     Standard_Real Ufloc=0., Ulloc=0.;
 762     AppDef_Compute theapprox (dmin,dmax,tol3d,tol2d,nbit);
 763     if (knownp) {
 764       math_Vector theParams(1,NbPoint);
 765       // On recale les parametres entre 0 et 1.
 766       /*Standard_Real*/ Ufloc = F.Parameter(Lin->Point(StPoint));
 767       /*Standard_Real*/ Ulloc = F.Parameter(Lin->Point(StPoint+NbPoint-1));
 768 //  modified by EAP Thu Jan  3 15:45:17 2002 ___END___
 769       for (i=1; i <= NbPoint; i++) {
 770         Standard_Integer iLin = StPoint + i - 1;
 771         theParams(i) = (F.Parameter(Lin->Point(iLin))-Ufloc)/(Ulloc - Ufloc);
 772       }
 773       AppDef_Compute theAppDef1(theParams,dmin,dmax,tol3d,tol2d,nbit, Standard_True,Standard_True);
 774       theapprox = theAppDef1;
 775     }
 776     theapprox.SetConstraints(Cfirst,Clast);
 777     theapprox.Perform(multL);
 778
 779 //  modified by EAP Thu Jan  3 16:00:43 2002 ___BEGIN___
 780     // To know internal parameters if multicurve is approximated by several Bezier's
 781     TColStd_SequenceOfReal aPoleDistSeq;
 782     Standard_Real aWholeDist=0;
 783 //  modified by EAP Thu Jan  3 16:45:48 2002 ___END___
 784     Standard_Real TheTol3d, TheTol2d;
 785     for (Standard_Integer Index=1; Index<=theapprox.NbMultiCurves(); Index++) {
 786       AppParCurves_MultiCurve& mucu = theapprox.ChangeValue(Index);
 787       theapprox.Error(Index, TheTol3d, TheTol2d);
 788       mytol3d = Max(TheTol3d/DX, mytol3d);
 789       mytol3d = Max(TheTol3d/DY, mytol3d);
 790       mytol3d = Max(TheTol3d/DZ, mytol3d);
 791       for(j = 1; j <= NbUPoles; j++){
 792         mucu.Transform(j,
 793                        -X/DX,1./DX,
 794                        -Y/DY,1./DY,
 795                        -Z/DZ,1./DZ);
 796       }
 797       for(j = 1; j <= NbPoles2d; j++){
 798         mucu.Transform2d(j + NbUPoles,
 799                          -X2d(j)/DX2d(j),1./DX2d(j),
 800                          -Y2d(j)/DY2d(j),1./DY2d(j));
 801         mytol2d = Max(TheTol2d/DX2d(j), mytol2d);
 802         mytol2d = Max(TheTol2d/DY2d(j), mytol2d);
 803       }
 804       concat.Append(mucu);
 805
 806 //  modified by EAP Thu Jan  3 15:45:23 2002 ___BEGIN___
 807       if (knownp && theapprox.NbMultiCurves() > 1)
 808       {
 809         gp_Pnt aFirstPole = mucu.Pole(Index, 1);
 810         gp_Pnt aLastPole  = mucu.Pole(Index, mucu.NbPoles());
 811         aPoleDistSeq.Append (aFirstPole.Distance(aLastPole));
 812         aWholeDist += aPoleDistSeq.Last();
 813       }
 814     }
 815     if (knownp)
 816     {
 817       Standard_Integer iDist;
 818       Standard_Real iU = Ufloc;
 819       for (iDist=1; iDist<aPoleDistSeq.Length(); iDist++)
 820       {
 821         iU += aPoleDistSeq(iDist) / aWholeDist * (Ulloc - Ufloc);
 822         //cout << "Internal: " << iU << endl;
 823         aParamSeq.Append(iU);
 824       }
 825       aParamSeq.Append(Ulloc);
 826     }
 827 //  modified by EAP Thu Jan  3 15:45:27 2002 ___END___
 828   }
 829 #ifdef DEB
 830   cout << "   Tolerances obtenues  --> 3d : "<< mytol3d << endl;
 831   cout << "                        --> 2d : "<< mytol2d << endl;
 832 #endif
 833   tol3dreached = mytol3d;
 834   tol2dreached = mytol2d;
 835   concat.Perform();
 836   const AppParCurves_MultiBSpCurve& multC = concat.Value();
 837   vdeg = multC.Degree();
 838   NbVPoles = multC.NbPoles();
 839
 840   tabPoles   = new TColgp_HArray2OfPnt (1,NbUPoles,1,NbVPoles);
 841   tabWeights = new TColStd_HArray2OfReal (1,NbUPoles,1,NbVPoles);
 842   tabVKnots  = new TColStd_HArray1OfReal (multC.Knots().Lower(),
 843                                           multC.Knots().Upper());
 844   tabVKnots->ChangeArray1() = multC.Knots();
 845
 846   if (knownp) {
 847 //  modified by EAP Fri Jan  4 12:07:30 2002 ___BEGIN___
 848     if (aParamSeq.Length() != tabVKnots->Length())
 849     {
 850       BSplCLib::Reparametrize(F.Parameter(Lin->Point(1)),
 851                               F.Parameter(Lin->Point(Lin->NbPoints())),
 852                               tabVKnots->ChangeArray1()
 853                               );
 854 #ifdef DEB
 855       cout << "Warning: AppBlend_AppSurf::Perform(), bad length of aParamSeq: " <<
 856         aParamSeq.Length() << " instead of " << tabVKnots->Length() << endl;
 857 #endif
 858     }
 859     else
 860     {
 861       Standard_Integer iKnot, iTabKnot = tabVKnots->Lower();
 862       for (iKnot=1; iKnot<=aParamSeq.Length(); iKnot++, iTabKnot++)
 863       {
 864         //cout << "Replace " << tabVKnots->Value(iTabKnot) << " with " << aParamSeq(iKnot) << endl;
 865         tabVKnots->SetValue(iTabKnot, aParamSeq(iKnot));
 866       }
 867     }
 868 //  modified by EAP Fri Jan  4 12:07:35 2002 ___END___
 869   }
 870
 871   tabVMults  = new TColStd_HArray1OfInteger (multC.Multiplicities().Lower(),
 872                                              multC.Multiplicities().Upper());
 873   tabVMults->ChangeArray1() = multC.Multiplicities();
 874
 875
 876   TColgp_Array1OfPnt newtabP(1,NbVPoles);
 877   Handle(TColgp_HArray1OfPnt2d) newtabP2d =
 878     new TColgp_HArray1OfPnt2d(1,NbVPoles);
 879   for (j=1; j <=NbUPoles; j++) {
 880     multC.Curve(j,newtabP);
 881     multC.Curve(j+NbUPoles+NbPoles2d,newtabP2d->ChangeArray1());
 882     for (k=1; k<=NbVPoles; k++) {
 883       // pour les courbes rationnelles il faut maintenant diviser
 884       // les poles par leurs poids respectifs
 885       tabPoles->ChangeValue(j,k).
 886         SetXYZ(newtabP(k).XYZ()/newtabP2d->Value(k).X());
 887       tabWeights->SetValue(j,k,newtabP2d->Value(k).X());
 888     }
 889   }
 890
 891   for (j=1; j<=NbPoles2d; j++) {
 892     newtabP2d = new TColgp_HArray1OfPnt2d(1,NbVPoles);
 893     multC.Curve(NbUPoles+j,newtabP2d->ChangeArray1());
 894     seqPoles2d.Append(newtabP2d);
 895   }
 896
 897   done = Standard_True;
 898 }
 899
 900
 901 //=======================================================================
 902 //function : SurfShape
 903 //purpose  :
 904 //=======================================================================
 905
 906 void AppBlend_AppSurf::SurfShape (Standard_Integer& UDegree,
 907                                   Standard_Integer& VDegree,
 908                                   Standard_Integer& NbUPoles,
 909                                   Standard_Integer& NbVPoles,
 910                                   Standard_Integer& NbUKnots,
 911                                   Standard_Integer& NbVKnots) const
 912 {
 913   if (!done) {StdFail_NotDone::Raise();}
 914   UDegree  = udeg;
 915   VDegree  = vdeg;
 916   NbUPoles = tabPoles->ColLength();
 917   NbVPoles = tabPoles->RowLength();
 918   NbUKnots = tabUKnots->Length();
 919   NbVKnots = tabVKnots->Length();
 920 }
 921
 922
 923 void AppBlend_AppSurf::Surface(TColgp_Array2OfPnt& TPoles,
 924                                TColStd_Array2OfReal& TWeights,
 925                                TColStd_Array1OfReal& TUKnots,
 926                                TColStd_Array1OfReal& TVKnots,
 927                                TColStd_Array1OfInteger& TUMults,
 928                                TColStd_Array1OfInteger& TVMults) const
 929
 930 {
 931   if (!done) {StdFail_NotDone::Raise();}
 932   TPoles   = tabPoles->Array2();
 933   TWeights = tabWeights->Array2();
 934   TUKnots  = tabUKnots->Array1();
 935   TUMults  = tabUMults->Array1();
 936   TVKnots  = tabVKnots->Array1();
 937   TVMults  = tabVMults->Array1();
 938 }
 939
 940 //=======================================================================
 941 //function : Curves2dShape
 942 //purpose  :
 943 //=======================================================================
 944
 945 void AppBlend_AppSurf::Curves2dShape(Standard_Integer& Degree,
 946                                      Standard_Integer& NbPoles,
 947                                      Standard_Integer& NbKnots) const
 948 {
 949   if (!done) {StdFail_NotDone::Raise();}
 950   if (seqPoles2d.Length() == 0) {Standard_DomainError::Raise();}
 951   Degree = vdeg;
 952   NbPoles = tabPoles->ColLength();
 953   NbKnots = tabVKnots->Length();
 954 }
 955
 956 //=======================================================================
 957 //function : Curve2d
 958 //purpose  :
 959 //=======================================================================
 960
 961 void AppBlend_AppSurf::Curve2d(const Standard_Integer Index,
 962                                TColgp_Array1OfPnt2d& TPoles,
 963                                TColStd_Array1OfReal& TKnots,
 964                                TColStd_Array1OfInteger& TMults) const
 965 {
 966   if (!done) {StdFail_NotDone::Raise();}
 967   if (seqPoles2d.Length() == 0) {Standard_DomainError::Raise();}
 968   TPoles = seqPoles2d(Index)->Array1();
 969   TKnots  = tabVKnots->Array1();
 970   TMults  = tabVMults->Array1();
 971 }
 972
 973 //=======================================================================
 974 //function : TolCurveOnSurf
 975 //purpose  :
 976 //=======================================================================
 977
 978 Standard_Real AppBlend_AppSurf::TolCurveOnSurf(const Standard_Integer) const
 979 {
 980   return tol3dreached; //On ne s'embete pas !!
 981 }
 982
 983
 984