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