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