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 ) Standard_DomainError::Raise();
 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).
 539         SetXYZ(newtabP(k).XYZ()/newtabP2d->Value(k).X());
 540       tabWeights->SetValue(j,k,newtabP2d->Value(k).X());
 541     }
 542   }
 543
 544   for (j=1; j<=NbPoles2d; j++) {
 545     newtabP2d = new TColgp_HArray1OfPnt2d(1,NbVPoles);
 546     multC.Curve(NbUPoles+j,newtabP2d->ChangeArray1());
 547     seqPoles2d.Append(newtabP2d);
 548   }
 549
 550   done = Standard_True;
 551 }
 552
 553
 554 //=======================================================================
 555 //function : Perform
 556 //purpose  :
 557 //=======================================================================
 558
 559 void AppBlend_AppSurf::Perform(const Handle(TheLine)& Lin,
 560                                TheSectionGenerator& F,
 561                                const Standard_Integer NbMaxP)
 562 {
 563   done = Standard_False;
 564   if (Lin.IsNull()) {return;}
 565   Standard_Integer i,j,k;
 566   Standard_Integer NbUPoles,NbUKnots,NbPoles2d,NbVPoles;
 567   Standard_Boolean withderiv;
 568   AppParCurves_Constraint Cfirst=AppParCurves_NoConstraint,Clast=AppParCurves_NoConstraint;
 569
 570   Standard_Real mytol3d = 0.0, mytol2d = 0.0;
 571   gp_XYZ newDv;
 572
 573   seqPoles2d.Clear();
 574
 575   Standard_Integer NbPointTot = Lin->NbPoints();
 576
 577   F.GetShape(NbUPoles,NbUKnots,udeg,NbPoles2d);
 578
 579   tabUKnots  = new TColStd_HArray1OfReal (1,NbUKnots);
 580   tabUMults  = new TColStd_HArray1OfInteger (1,NbUKnots);
 581
 582   F.Knots(tabUKnots->ChangeArray1());
 583   F.Mults(tabUMults->ChangeArray1());
 584
 585   TColgp_Array1OfPnt tabAppP(1,NbUPoles);
 586   TColgp_Array1OfVec tabAppV(1,NbUPoles);
 587   Standard_Real X,Y,Z,DX,DY,DZ;
 588   X = Y = Z = RealLast();
 589   DX = DY = DZ = RealFirst();
 590
 591   TColgp_Array1OfPnt2d tabP2d(1,Max(1,NbPoles2d));
 592   TColgp_Array1OfVec2d tabV2d(1,Max(1,NbPoles2d));
 593   TColStd_Array1OfReal X2d(1,Max(1,NbPoles2d));X2d.Init(RealLast());
 594   TColStd_Array1OfReal Y2d(1,Max(1,NbPoles2d));Y2d.Init(RealLast());
 595   TColStd_Array1OfReal DX2d(1,Max(1,NbPoles2d));DX2d.Init(RealFirst());
 596   TColStd_Array1OfReal DY2d(1,Max(1,NbPoles2d));DY2d.Init(RealFirst());
 597
 598   TColStd_Array1OfReal tabW(1,NbUPoles),tabDW(1,NbUPoles);
 599
 600   TColgp_Array1OfPnt2d tabAppP2d(1,NbPoles2d+NbUPoles); // points2d + poids
 601   TColgp_Array1OfVec2d tabAppV2d(1,NbPoles2d+NbUPoles);
 602
 603   // On calcule les boites de chaque ligne (box for all lines)
 604   for(i = 1; i <= NbPointTot; i++){
 605     F.Section(Lin->Point(i),tabAppP,tabAppV,tabP2d,tabV2d,tabW,tabDW);
 606     Standard_Real x,y,z;
 607     for(j = 1; j <= NbUPoles; j++){
 608       tabAppP(j).Coord(x,y,z);
 609       if(x < X) X = x; if(x > DX) DX = x;
 610       if(y < Y) Y = y; if(y > DY) DY = y;
 611       if(z < Z) Z = z; if(z > DZ) DZ = z;
 612     }
 613     for(j = 1; j <= NbPoles2d; j++){
 614       tabP2d(j).Coord(x,y);
 615       if(x < X2d(j)) X2d(j) = x; if(x > DX2d(j)) DX2d(j) = x;
 616       if(y < Y2d(j)) Y2d(j) = y; if(y > DY2d(j)) DY2d(j) = y;
 617     }
 618   }
 619   // On calcule pour chaque ligne la transformation vers 0 1.
 620   Standard_Real seuil = 1000.*tol3d;
 621   Standard_Real seuil2d = 1000.*tol2d;
 622   if((DX - X) < seuil ){ DX = 1.; X = 0.; }
 623   else{ DX = 1./(DX - X); X *= -DX; }
 624   if((DY - Y) < seuil){ DY = 1.; Y = 0.; }
 625   else{ DY = 1./(DY - Y); Y *= -DY; }
 626   if((DZ - Z) < seuil){ DZ = 1.; Z = 0.; }
 627   else{ DZ = 1./(DZ - Z); Z *= -DZ; }
 628   for(j = 1; j <= NbPoles2d; j++){
 629     if((DX2d(j) - X2d(j)) < seuil2d){ DX2d(j) = 1.; X2d(j) = 0.; }
 630     else{ DX2d(j) = 1./(DX2d(j) - X2d(j)); X2d(j) *= -DX2d(j); }
 631     if((DY2d(j) - Y2d(j)) < seuil2d){ DY2d(j) = 1.; Y2d(j) = 0.; }
 632     else{ DY2d(j) = 1./(DY2d(j) - Y2d(j)); Y2d(j) *= -DY2d(j); }
 633   }
 634   if(!scal){
 635     DX = 1.; X = 0.;
 636     DY = 1.; Y = 0.;
 637     DZ = 1.; Z = 0.;
 638     for(j = 1; j <= NbPoles2d; j++){
 639       DX2d(j) = 1.; X2d(j) = 0.;
 640       DY2d(j) = 1.; Y2d(j) = 0.;
 641     }
 642   }
 643 //  modified by eap Thu Jan  3 14:45:22 2002 ___BEGIN___
 644   // Keep "inter-troncons" parameters, not only first and last
 645 //  Standard_Real Ufirst=0,Ulast=0;
 646   TColStd_SequenceOfReal aParamSeq;
 647    if (knownp) {
 648 //     Ufirst = F.Parameter(Lin->Point(1));
 649 //     Ulast = F.Parameter(Lin->Point(NbPointTot));
 650      aParamSeq.Append( F.Parameter (Lin->Point(1)) );
 651   }
 652 //  modified by EAP Thu Jan  3 14:45:41 2002 ___END___
 653
 654   Approx_MCurvesToBSpCurve concat;
 655
 656   //On calcule le nombre de troncons.
 657   Standard_Integer nbtronc = NbPointTot/NbMaxP;
 658   Standard_Integer reste = NbPointTot - (nbtronc * NbMaxP);
 659   // On regarde si il faut prendre un troncon de plus.
 660   Standard_Integer nmax = NbMaxP;
 661   if(nbtronc > 0 && reste > 0){
 662     nmax = NbPointTot/(nbtronc + 1);
 663     if(nmax > (2*NbMaxP)/3) {
 664       nbtronc++;
 665       reste = NbPointTot - (nbtronc * nmax);
 666     }
 667     else nmax = NbMaxP;
 668   }
 669   else if(nbtronc == 0){
 670     nbtronc = 1;
 671     nmax = reste;
 672     reste = 0;
 673   }
 674
 675   // Approximate each "troncon" with nb of Bezier's using AppDef_Compute
 676   // and concat them into BSpline with Approx_MCurvesToBSpCurve
 677
 678   TColStd_Array1OfInteger troncsize(1,nbtronc);
 679   TColStd_Array1OfInteger troncstart(1,nbtronc);
 680
 681   Standard_Integer rab = reste/nbtronc + 1;
 682   Standard_Integer start = 1;
 683   Standard_Integer itronc ;
 684   for( itronc = 1; itronc <= nbtronc; itronc++){
 685     troncstart(itronc) = start;
 686     Standard_Integer rabrab = Min(rab,reste);
 687     if(reste > 0){ reste -= rabrab; }
 688     troncsize(itronc) = nmax + rabrab + 1;
 689     start += (nmax + rabrab);
 690   }
 691   troncsize(nbtronc) = troncsize(nbtronc) - 1;
 692   for(itronc = 1; itronc <= nbtronc; itronc++){
 693     Standard_Integer NbPoint = troncsize(itronc);
 694     Standard_Integer StPoint = troncstart(itronc);
 695     AppDef_MultiPointConstraint multP;
 696     AppDef_MultiLine multL(NbPoint);
 697
 698     for (i=1; i<=NbPoint; i++) {
 699       Standard_Integer iLin = StPoint + i - 1;
 700       Standard_Real x,y,z;
 701       withderiv = F.Section(Lin->Point(iLin),tabAppP,tabAppV,tabP2d,tabV2d,
 702                             tabW,tabDW);
 703       if(AppBlend_GetContextApproxWithNoTgt()) withderiv = Standard_False;
 704
 705       for (j=1; j<=NbPoles2d; j++) {
 706         tabP2d(j).Coord(x,y);
 707         tabAppP2d(j).SetCoord(DX2d(j)*x+X2d(j),DY2d(j)*y+Y2d(j));
 708         if (withderiv) {
 709           tabV2d(j).Coord(x,y);
 710           tabAppV2d(j).SetCoord(DX2d(j)*x,DY2d(j)*y);
 711         }
 712       }
 713       for (j=1; j<=NbUPoles; j++) {
 714         // pour les courbes rationnelles il faut multiplier les poles par
 715         // leurs poids respectifs
 716         if (withderiv) {
 717           tabAppV2d(NbPoles2d+j).SetCoord(tabDW(j),0.);
 718           newDv.SetLinearForm(tabDW(j),tabAppP(j).XYZ(),tabW(j),tabAppV(j).XYZ());
 719           tabAppV(j).SetCoord(DX*newDv.X(),DY*newDv.Y(),DZ*newDv.Z());
 720         }
 721         tabAppP(j).SetXYZ(tabAppP(j).XYZ() * tabW(j));
 722         tabAppP2d(NbPoles2d+j).SetCoord(tabW(j),0.);
 723         tabAppP(j).Coord(x,y,z);
 724         tabAppP(j).SetCoord(DX*x+X,DY*y+Y,DZ*z+Z);
 725       }
 726       if (withderiv) {
 727         multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d,tabAppV,tabAppV2d);
 728         if(i == 1) Cfirst = AppParCurves_TangencyPoint;
 729         else if(i == NbPoint) Clast = AppParCurves_TangencyPoint;
 730       }
 731       else {
 732         multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d);
 733         if(i == 1) Cfirst = AppParCurves_PassPoint;
 734         else if(i == NbPoint) Clast = AppParCurves_PassPoint;
 735       }
 736       multL.SetValue(i,multP);
 737     }
 738
 739
 740   //IFV 04.06.07 occ13904
 741     if(NbPoint == 2) {
 742       dmin = 1;
 743       if(Cfirst == AppParCurves_PassPoint && Clast == AppParCurves_PassPoint) {
 744         dmax = 1;
 745       }
 746     }
 747
 748 //  modified by EAP Thu Jan  3 15:44:13 2002 ___BEGIN___
 749     Standard_Real Ufloc=0., Ulloc=0.;
 750     AppDef_Compute theapprox (dmin,dmax,tol3d,tol2d,nbit);
 751     if (knownp) {
 752       math_Vector theParams(1,NbPoint);
 753       // On recale les parametres entre 0 et 1.
 754       /*Standard_Real*/ Ufloc = F.Parameter(Lin->Point(StPoint));
 755       /*Standard_Real*/ Ulloc = F.Parameter(Lin->Point(StPoint+NbPoint-1));
 756 //  modified by EAP Thu Jan  3 15:45:17 2002 ___END___
 757       for (i=1; i <= NbPoint; i++) {
 758         Standard_Integer iLin = StPoint + i - 1;
 759         theParams(i) = (F.Parameter(Lin->Point(iLin))-Ufloc)/(Ulloc - Ufloc);
 760       }
 761       AppDef_Compute theAppDef1(theParams,dmin,dmax,tol3d,tol2d,nbit, Standard_True,Standard_True);
 762       theapprox = theAppDef1;
 763     }
 764     theapprox.SetConstraints(Cfirst,Clast);
 765     theapprox.Perform(multL);
 766
 767 //  modified by EAP Thu Jan  3 16:00:43 2002 ___BEGIN___
 768     // To know internal parameters if multicurve is approximated by several Bezier's
 769     TColStd_SequenceOfReal aPoleDistSeq;
 770     Standard_Real aWholeDist=0;
 771 //  modified by EAP Thu Jan  3 16:45:48 2002 ___END___
 772     Standard_Real TheTol3d, TheTol2d;
 773     for (Standard_Integer Index=1; Index<=theapprox.NbMultiCurves(); Index++) {
 774       AppParCurves_MultiCurve& mucu = theapprox.ChangeValue(Index);
 775       theapprox.Error(Index, TheTol3d, TheTol2d);
 776       mytol3d = Max(TheTol3d/DX, mytol3d);
 777       mytol3d = Max(TheTol3d/DY, mytol3d);
 778       mytol3d = Max(TheTol3d/DZ, mytol3d);
 779       for(j = 1; j <= NbUPoles; j++){
 780         mucu.Transform(j,
 781                        -X/DX,1./DX,
 782                        -Y/DY,1./DY,
 783                        -Z/DZ,1./DZ);
 784       }
 785       for(j = 1; j <= NbPoles2d; j++){
 786         mucu.Transform2d(j + NbUPoles,
 787                          -X2d(j)/DX2d(j),1./DX2d(j),
 788                          -Y2d(j)/DY2d(j),1./DY2d(j));
 789         mytol2d = Max(TheTol2d/DX2d(j), mytol2d);
 790         mytol2d = Max(TheTol2d/DY2d(j), mytol2d);
 791       }
 792       concat.Append(mucu);
 793
 794 //  modified by EAP Thu Jan  3 15:45:23 2002 ___BEGIN___
 795       if (knownp && theapprox.NbMultiCurves() > 1)
 796       {
 797         gp_Pnt aFirstPole = mucu.Pole(Index, 1);
 798         gp_Pnt aLastPole  = mucu.Pole(Index, mucu.NbPoles());
 799         aPoleDistSeq.Append (aFirstPole.Distance(aLastPole));
 800         aWholeDist += aPoleDistSeq.Last();
 801       }
 802     }
 803     if (knownp)
 804     {
 805       Standard_Integer iDist;
 806       Standard_Real iU = Ufloc;
 807       for (iDist=1; iDist<aPoleDistSeq.Length(); iDist++)
 808       {
 809         iU += aPoleDistSeq(iDist) / aWholeDist * (Ulloc - Ufloc);
 810         //cout << "Internal: " << iU << endl;
 811         aParamSeq.Append(iU);
 812       }
 813       aParamSeq.Append(Ulloc);
 814     }
 815 //  modified by EAP Thu Jan  3 15:45:27 2002 ___END___
 816   }
 817 #ifdef OCCT_DEBUG
 818   cout << "   Tolerances obtenues  --> 3d : "<< mytol3d << endl;
 819   cout << "                        --> 2d : "<< mytol2d << endl;
 820 #endif
 821   tol3dreached = mytol3d;
 822   tol2dreached = mytol2d;
 823   concat.Perform();
 824   const AppParCurves_MultiBSpCurve& multC = concat.Value();
 825   vdeg = multC.Degree();
 826   NbVPoles = multC.NbPoles();
 827
 828   tabPoles   = new TColgp_HArray2OfPnt (1,NbUPoles,1,NbVPoles);
 829   tabWeights = new TColStd_HArray2OfReal (1,NbUPoles,1,NbVPoles);
 830   tabVKnots  = new TColStd_HArray1OfReal (multC.Knots().Lower(),
 831                                           multC.Knots().Upper());
 832   tabVKnots->ChangeArray1() = multC.Knots();
 833
 834   if (knownp) {
 835 //  modified by EAP Fri Jan  4 12:07:30 2002 ___BEGIN___
 836     if (aParamSeq.Length() != tabVKnots->Length())
 837     {
 838       BSplCLib::Reparametrize(F.Parameter(Lin->Point(1)),
 839                               F.Parameter(Lin->Point(Lin->NbPoints())),
 840                               tabVKnots->ChangeArray1()
 841                               );
 842 #ifdef OCCT_DEBUG
 843       cout << "Warning: AppBlend_AppSurf::Perform(), bad length of aParamSeq: " <<
 844         aParamSeq.Length() << " instead of " << tabVKnots->Length() << endl;
 845 #endif
 846     }
 847     else
 848     {
 849       Standard_Integer iKnot, iTabKnot = tabVKnots->Lower();
 850       for (iKnot=1; iKnot<=aParamSeq.Length(); iKnot++, iTabKnot++)
 851       {
 852         //cout << "Replace " << tabVKnots->Value(iTabKnot) << " with " << aParamSeq(iKnot) << endl;
 853         tabVKnots->SetValue(iTabKnot, aParamSeq(iKnot));
 854       }
 855     }
 856 //  modified by EAP Fri Jan  4 12:07:35 2002 ___END___
 857   }
 858
 859   tabVMults  = new TColStd_HArray1OfInteger (multC.Multiplicities().Lower(),
 860                                              multC.Multiplicities().Upper());
 861   tabVMults->ChangeArray1() = multC.Multiplicities();
 862
 863
 864   TColgp_Array1OfPnt newtabP(1,NbVPoles);
 865   Handle(TColgp_HArray1OfPnt2d) newtabP2d =
 866     new TColgp_HArray1OfPnt2d(1,NbVPoles);
 867   for (j=1; j <=NbUPoles; j++) {
 868     multC.Curve(j,newtabP);
 869     multC.Curve(j+NbUPoles+NbPoles2d,newtabP2d->ChangeArray1());
 870     for (k=1; k<=NbVPoles; k++) {
 871       // pour les courbes rationnelles il faut maintenant diviser
 872       // les poles par leurs poids respectifs
 873       tabPoles->ChangeValue(j,k).
 874         SetXYZ(newtabP(k).XYZ()/newtabP2d->Value(k).X());
 875       tabWeights->SetValue(j,k,newtabP2d->Value(k).X());
 876     }
 877   }
 878
 879   for (j=1; j<=NbPoles2d; j++) {
 880     newtabP2d = new TColgp_HArray1OfPnt2d(1,NbVPoles);
 881     multC.Curve(NbUPoles+j,newtabP2d->ChangeArray1());
 882     seqPoles2d.Append(newtabP2d);
 883   }
 884
 885   done = Standard_True;
 886 }
 887
 888
 889 //=======================================================================
 890 //function : SurfShape
 891 //purpose  :
 892 //=======================================================================
 893
 894 void AppBlend_AppSurf::SurfShape (Standard_Integer& UDegree,
 895                                   Standard_Integer& VDegree,
 896                                   Standard_Integer& NbUPoles,
 897                                   Standard_Integer& NbVPoles,
 898                                   Standard_Integer& NbUKnots,
 899                                   Standard_Integer& NbVKnots) const
 900 {
 901   if (!done) {StdFail_NotDone::Raise();}
 902   UDegree  = udeg;
 903   VDegree  = vdeg;
 904   NbUPoles = tabPoles->ColLength();
 905   NbVPoles = tabPoles->RowLength();
 906   NbUKnots = tabUKnots->Length();
 907   NbVKnots = tabVKnots->Length();
 908 }
 909
 910
 911 void AppBlend_AppSurf::Surface(TColgp_Array2OfPnt& TPoles,
 912                                TColStd_Array2OfReal& TWeights,
 913                                TColStd_Array1OfReal& TUKnots,
 914                                TColStd_Array1OfReal& TVKnots,
 915                                TColStd_Array1OfInteger& TUMults,
 916                                TColStd_Array1OfInteger& TVMults) const
 917
 918 {
 919   if (!done) {StdFail_NotDone::Raise();}
 920   TPoles   = tabPoles->Array2();
 921   TWeights = tabWeights->Array2();
 922   TUKnots  = tabUKnots->Array1();
 923   TUMults  = tabUMults->Array1();
 924   TVKnots  = tabVKnots->Array1();
 925   TVMults  = tabVMults->Array1();
 926 }
 927
 928 //=======================================================================
 929 //function : Curves2dShape
 930 //purpose  :
 931 //=======================================================================
 932
 933 void AppBlend_AppSurf::Curves2dShape(Standard_Integer& Degree,
 934                                      Standard_Integer& NbPoles,
 935                                      Standard_Integer& NbKnots) const
 936 {
 937   if (!done) {StdFail_NotDone::Raise();}
 938   if (seqPoles2d.Length() == 0) {Standard_DomainError::Raise();}
 939   Degree = vdeg;
 940   NbPoles = tabPoles->ColLength();
 941   NbKnots = tabVKnots->Length();
 942 }
 943
 944 //=======================================================================
 945 //function : Curve2d
 946 //purpose  :
 947 //=======================================================================
 948
 949 void AppBlend_AppSurf::Curve2d(const Standard_Integer Index,
 950                                TColgp_Array1OfPnt2d& TPoles,
 951                                TColStd_Array1OfReal& TKnots,
 952                                TColStd_Array1OfInteger& TMults) const
 953 {
 954   if (!done) {StdFail_NotDone::Raise();}
 955   if (seqPoles2d.Length() == 0) {Standard_DomainError::Raise();}
 956   TPoles = seqPoles2d(Index)->Array1();
 957   TKnots  = tabVKnots->Array1();
 958   TMults  = tabVMults->Array1();
 959 }
 960
 961 //=======================================================================
 962 //function : TolCurveOnSurf
 963 //purpose  :
 964 //=======================================================================
 965
 966 Standard_Real AppBlend_AppSurf::TolCurveOnSurf(const Standard_Integer) const
 967 {
 968   return tol3dreached; //On ne s'embete pas !!
 969 }
 970
 971
 972