src/Approx/Approx_ComputeLine.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 <Approx_ParametrizationType.hxx>
  16 #include Approx_ParLeastSquareOfMyGradient_hxx
  17 #include <TColStd_Array1OfReal.hxx>
  18 #include <TColgp_Array1OfPnt.hxx>
  19 #include <TColgp_Array1OfPnt2d.hxx>
  20 #include <gp_Pnt.hxx>
  21 #include <gp_Pnt2d.hxx>
  22 #include <gp_Vec.hxx>
  23 #include <gp_Vec2d.hxx>
  24 #include <TColgp_Array1OfVec.hxx>
  25 #include <TColgp_Array1OfVec2d.hxx>
  26 #include <AppParCurves_Constraint.hxx>
  27 #include <AppParCurves_HArray1OfConstraintCouple.hxx>
  28 #include <AppParCurves_MultiPoint.hxx>
  29 #include <Precision.hxx>
  30 #include <math_IntegerVector.hxx>
  31 #include <math_Gauss.hxx>
  32 #include <math_Uzawa.hxx>
  33 #include <Approx_MCurvesToBSpCurve.hxx>
  34 #include <AppParCurves_ConstraintCouple.hxx>
  35
  36 #include <stdio.h>
  37
  38 #ifdef OCCT_DEBUG
  39 static Standard_Boolean mydebug = Standard_False;
  40
  41 #include <Geom_BezierCurve.hxx>
  42 #include <Geom2d_BezierCurve.hxx>
  43 #ifdef DRAW
  44 #include <DrawTrSurf.hxx>
  45 #include <Draw.hxx>
  46 #include <Draw_Appli.hxx>
  47 #endif
  48
  49 static void DUMP(const MultiLine& Line)
  50 {
  51   Standard_Integer i, j, nbP2d, nbP3d, firstP, lastP;
  52   gp_Pnt P1;
  53   gp_Pnt2d P12d;
  54
  55   firstP = LineTool::FirstPoint(Line);
  56   lastP  = LineTool::LastPoint(Line);
  57
  58   nbP3d = LineTool::NbP3d(Line);
  59   nbP2d = LineTool::NbP2d(Line);
  60   Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
  61   if (nbP3d == 0) mynbP3d = 1;
  62   if (nbP2d == 0) mynbP2d = 1;
  63
  64   TColgp_Array1OfPnt tabP(1, mynbP3d);
  65   TColgp_Array1OfPnt2d tabP2d(1, mynbP2d);
  66
  67   cout <<"DUMP de la MultiLine entre "<<firstP <<" et "<<lastP<<": "<<endl;
  68   for (i = firstP; i <= lastP; i++) {
  69     if (nbP3d != 0 && nbP2d != 0) LineTool::Value(Line, i, tabP, tabP2d);
  70     else if (nbP2d != 0)          LineTool::Value(Line, i, tabP2d);
  71     else if (nbP3d != 0)          LineTool::Value(Line, i, tabP);
  72
  73     cout << "point "<<i<<":"<< endl;
  74     for (j = 1; j <= nbP3d; j++) {
  75       P1 = tabP(j);
  76       cout <<P1.X()<<" "<<P1.Y()<<" "<<P1.Z()<<endl;
  77     }
  78     for (j = 1; j <= nbP2d; j++) {
  79       P12d = tabP2d(j);
  80       cout <<P12d.X()<<" "<<P12d.Y()<<endl;
  81     }
  82   }
  83
  84 }
  85
  86 static void DUMP(const AppParCurves_MultiCurve& C) {
  87   static Standard_Integer nbappel = 0;
  88   Standard_Integer i;
  89   Standard_Integer nbpoles = C.NbPoles();
  90
  91   Handle(Geom_BezierCurve) BSp;
  92   Handle(Geom2d_BezierCurve) BSp2d;
  93
  94   TColgp_Array1OfPnt tabPoles(1, nbpoles);
  95   TColgp_Array1OfPnt2d tabPoles2d(1, nbpoles);
  96   char solname[100];
  97
  98   nbappel++;
  99   for (i = 1; i <= C.NbCurves(); i++) {
 100     if (C.Dimension(i) == 3) {
 101       C.Curve(i, tabPoles);
 102       BSp = new Geom_BezierCurve(tabPoles);
 103       sprintf(solname, "%s%i%s_%i", "c", i, "3d", nbappel);
 104 #ifdef DRAW
 105       char* Temp = solname;
 106       DrawTrSurf::Set(Temp, BSp);
 107 //      DrawTrSurf::Set(solname, BSp);
 108 #endif
 109     }
 110     else {
 111       C.Curve(i, tabPoles2d);
 112       BSp2d = new Geom2d_BezierCurve(tabPoles2d);
 113       sprintf(solname, "%s%i%s_%i", "c", i, "2d", nbappel);
 114 #ifdef DRAW
 115       char* Temp = solname;
 116       DrawTrSurf::Set(Temp, BSp2d);
 117 //      DrawTrSurf::Set(solname, BSp2d);
 118 #endif
 119     }
 120   }
 121 #ifdef DRAW
 122   dout.Flush();
 123 #endif
 124 }
 125
 126 #endif
 127
 128 static Standard_Boolean CheckMultiCurve(const AppParCurves_MultiCurve& theMultiCurve,
 129                                         const MultiLine& theLine,
 130                                         const Standard_Integer theIndfirst,
 131                                         const Standard_Integer theIndlast,
 132                                         Standard_Integer& theIndbad)
 133 {
 134   const Standard_Integer nbp3d = LineTool::NbP3d(theLine);
 135   const Standard_Integer nbp2d = LineTool::NbP2d(theLine);
 136
 137   if (nbp3d > 1) //only simple cases are analysed
 138     return Standard_True;
 139
 140   const Standard_Real MinScalProd = -0.9;
 141   const Standard_Real SqTol3d = Precision::SquareConfusion();
 142
 143   theIndbad = 0;
 144   Standard_Integer indbads [4];
 145   indbads[1] = indbads[2] = indbads[3] = 0;
 146
 147   Standard_Integer NbCur = theMultiCurve.NbCurves();
 148   Standard_Boolean LoopFound = Standard_False;
 149
 150   Standard_Integer aNbP3d = Max(nbp3d, 1);
 151   Standard_Integer aNbP2d = Max(nbp2d, 1);
 152
 153   TColgp_Array1OfPnt tabP(1, aNbP3d);
 154   TColgp_Array1OfPnt2d tabP2d(1, aNbP2d);
 155
 156 #ifdef DRAW
 157   char* name = new char[100];
 158   Standard_Integer nbbc = 1;
 159   Standard_Integer indc = 1;
 160 #endif
 161   if (theMultiCurve.Dimension(1) == 3 /*myNbP3d == 1*/)
 162   {
 163     TColgp_Array1OfPnt aPoles(1, theMultiCurve.NbPoles());
 164     theMultiCurve.Curve(1, aPoles);
 165 #ifdef DRAW
 166     Handle(Geom_Curve) theBezier = new Geom_BezierCurve(aPoles);
 167     sprintf(name, "bc3d_%d_%d", indc, nbbc);
 168     DrawTrSurf::Set(name, theBezier);
 169 #endif
 170     gp_Vec FirstVec, SecondVec;
 171     Standard_Integer indp = 2;
 172     while (indp <= aPoles.Upper())
 173     {
 174       FirstVec = gp_Vec(aPoles(1), aPoles(indp++));
 175       Standard_Real aLength = FirstVec.Magnitude();
 176       if (aLength > gp::Resolution())
 177       {
 178         FirstVec /= aLength;
 179         break;
 180       }
 181     }
 182     gp_Pnt MidPnt = aPoles(indp-1);
 183     //for (Standard_Integer k = 3; k <= aPoles.Upper(); k++)
 184     while (indp <= aPoles.Upper())
 185     {
 186       SecondVec = gp_Vec(MidPnt, aPoles(indp));
 187       Standard_Real aLength = SecondVec.Magnitude();
 188       if (aLength <= gp::Resolution())
 189       {
 190         indp++;
 191         continue;
 192       }
 193       SecondVec /= aLength;
 194       Standard_Real ScalProd = FirstVec * SecondVec;
 195       if (ScalProd < MinScalProd)
 196       {
 197 #ifdef DRAW
 198         cout<<"ScalProd("<<indp-2<<","<<indp-1<<")-("<<indp-1<<","<<indp<<") = "<<ScalProd<<endl;
 199 #endif
 200         LoopFound = Standard_True;
 201         break;
 202       }
 203       FirstVec = SecondVec;
 204       MidPnt = aPoles(indp);
 205       indp++;
 206     }
 207     //Check: may be it is a real loop
 208     if (LoopFound)
 209     {
 210       for (Standard_Integer FirstInd = theIndfirst;
 211            FirstInd <= theIndlast - 2; FirstInd++)
 212       {
 213         LineTool::Value(theLine, FirstInd, tabP);
 214         gp_Pnt FirstPnt = tabP(1);
 215         for (Standard_Integer k = FirstInd+1; k < theIndlast; k++)
 216         {
 217           LineTool::Value(theLine, k, tabP);
 218           gp_Pnt Pnt1 = tabP(1);
 219           LineTool::Value(theLine, k+1, tabP);
 220           gp_Pnt Pnt2 = tabP(1);
 221           if (FirstPnt.SquareDistance(Pnt1) <= SqTol3d ||
 222               FirstPnt.SquareDistance(Pnt2) <= SqTol3d)
 223           {
 224             LoopFound = Standard_False;
 225             break;
 226           }
 227           gp_Vec Vec1(FirstPnt, Pnt1);
 228           Vec1.Normalize();
 229           gp_Vec Vec2(FirstPnt, Pnt2);
 230           Vec2.Normalize();
 231           Standard_Real ScalProd = Vec1 * Vec2;
 232           if (ScalProd < MinScalProd)
 233           {
 234             LoopFound = Standard_False;
 235             break;
 236           }
 237         }
 238         if (LoopFound == Standard_False)
 239           break;
 240       }
 241     }
 242     if (LoopFound)
 243     {
 244       //search <indbad>
 245       Standard_Real MaxSqDist = 0.;
 246       for (Standard_Integer k = theIndfirst+1; k <= theIndlast; k++)
 247       {
 248         LineTool::Value(theLine, k-1, tabP);
 249         gp_Pnt PrevPnt = tabP(1);
 250         LineTool::Value(theLine, k, tabP);
 251         gp_Pnt CurPnt  = tabP(1);
 252         Standard_Real aSqDist = PrevPnt.SquareDistance(CurPnt);
 253         if (aSqDist > MaxSqDist)
 254         {
 255           MaxSqDist = aSqDist;
 256           indbads[1] = k;
 257         }
 258       }
 259       for (Standard_Integer indcur = 2; indcur <= NbCur; indcur++)
 260       {
 261         MaxSqDist = 0.;
 262         for (Standard_Integer k = theIndfirst+1; k <= theIndlast; k++)
 263         {
 264           LineTool::Value(theLine, k-1, tabP2d);
 265           gp_Pnt2d PrevPnt = tabP2d(indcur-1);
 266           LineTool::Value(theLine, k, tabP2d);
 267           gp_Pnt2d CurPnt  = tabP2d(indcur-1);
 268           Standard_Real aSqDist = PrevPnt.SquareDistance(CurPnt);
 269           if (aSqDist > MaxSqDist)
 270           {
 271             MaxSqDist = aSqDist;
 272             indbads[indcur] = k;
 273           }
 274         }
 275       }
 276     }
 277   } //if (myNbP3d == 1)
 278   else //2d case
 279   {
 280     TColgp_Array1OfPnt2d aPoles2d(1, theMultiCurve.NbPoles());
 281     theMultiCurve.Curve(1, aPoles2d);
 282 #ifdef DRAW
 283     Handle(Geom2d_Curve) theBezier2d = new Geom2d_BezierCurve(aPoles2d);
 284     sprintf(name, "bc2d_%d_%d", indc, nbbc);
 285     DrawTrSurf::Set(name, theBezier2d);
 286 #endif
 287     gp_Vec2d FirstVec, SecondVec;
 288     FirstVec = gp_Vec2d(aPoles2d(1), aPoles2d(2));
 289     FirstVec.Normalize();
 290     gp_Pnt2d MidPnt = aPoles2d(2);
 291     for (Standard_Integer k = 3; k <= aPoles2d.Upper(); k++)
 292     {
 293       SecondVec = gp_Vec2d(MidPnt, aPoles2d(k));
 294       SecondVec.Normalize();
 295       Standard_Real ScalProd = FirstVec * SecondVec;
 296       if (ScalProd < MinScalProd)
 297       {
 298 #ifdef DRAW
 299         cout<<"ScalProd("<<k-2<<","<<k-1<<")-("<<k-1<<","<<k<<") = "<<ScalProd<<endl;
 300 #endif
 301         LoopFound = Standard_True;
 302         break;
 303       }
 304       FirstVec = SecondVec;
 305       MidPnt = aPoles2d(k);
 306     }
 307     //Check: may be it is a real loop
 308     if (LoopFound)
 309     {
 310       for (Standard_Integer FirstInd = theIndfirst;
 311            FirstInd <= theIndlast - 2; FirstInd++)
 312       {
 313         LineTool::Value(theLine, FirstInd, tabP2d);
 314         gp_Pnt2d FirstPnt = tabP2d(1);
 315         for (Standard_Integer k = FirstInd+1; k < theIndlast; k++)
 316         {
 317           LineTool::Value(theLine, k, tabP2d);
 318           gp_Pnt2d Pnt1 = tabP2d(1);
 319           LineTool::Value(theLine, k+1, tabP2d);
 320           gp_Pnt2d Pnt2 = tabP2d(1);
 321           if (FirstPnt.SquareDistance(Pnt1) <= SqTol3d ||
 322               FirstPnt.SquareDistance(Pnt2) <= SqTol3d)
 323           {
 324             LoopFound = Standard_False;
 325             break;
 326           }
 327           gp_Vec2d Vec1(FirstPnt, Pnt1);
 328           Vec1.Normalize();
 329           gp_Vec2d Vec2(FirstPnt, Pnt2);
 330           Vec2.Normalize();
 331           Standard_Real ScalProd = Vec1 * Vec2;
 332           if (ScalProd < MinScalProd)
 333           {
 334             LoopFound = Standard_False;
 335             break;
 336           }
 337         }
 338         if (LoopFound == Standard_False)
 339           break;
 340       }
 341     }
 342     if (LoopFound)
 343     {
 344       //search <indbad>
 345       for (Standard_Integer indcur = 1; indcur <= NbCur; indcur++)
 346       {
 347         Standard_Real MaxSqDist = 0.;
 348         for (Standard_Integer k = theIndfirst+1; k <= theIndlast; k++)
 349         {
 350           LineTool::Value(theLine, k-1, tabP2d);
 351           gp_Pnt2d PrevPnt = tabP2d(indcur);
 352           LineTool::Value(theLine, k, tabP2d);
 353           gp_Pnt2d CurPnt  = tabP2d(indcur);
 354           Standard_Real aSqDist = PrevPnt.SquareDistance(CurPnt);
 355           if (aSqDist > MaxSqDist)
 356           {
 357             MaxSqDist = aSqDist;
 358             indbads[indcur] = k;
 359           }
 360         }
 361       }
 362     }
 363   }
 364
 365   //Define <indbad>
 366   if (indbads[1] != 0 && indbads[2] != 0)
 367   {
 368     if (indbads[1] != indbads[2])
 369       LoopFound = Standard_False;
 370     else if (indbads[3] != 0 && indbads[1] != indbads[3])
 371       LoopFound = Standard_False;
 372   }
 373   if (LoopFound)
 374     theIndbad = indbads[1];
 375
 376   return (!LoopFound);
 377 }
 378
 379 void Approx_ComputeLine::FirstTangencyVector(const MultiLine&       Line,
 380                                              const Standard_Integer index,
 381                                              math_Vector&           V) const
 382 {
 383
 384   Standard_Integer i, j, nbP2d, nbP3d;
 385   nbP3d = LineTool::NbP3d(Line);
 386   nbP2d = LineTool::NbP2d(Line);
 387   Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
 388   if (nbP3d == 0) mynbP3d = 1;
 389   if (nbP2d == 0) mynbP2d = 1;
 390   Standard_Boolean Ok=Standard_False;
 391   TColgp_Array1OfVec TabV(1, mynbP3d);
 392   TColgp_Array1OfVec2d TabV2d(1, mynbP2d);
 393
 394   if (nbP3d != 0 && nbP2d != 0)
 395     Ok = LineTool::Tangency(Line, index, TabV, TabV2d);
 396   else if (nbP2d != 0)
 397     Ok = LineTool::Tangency(Line, index, TabV2d);
 398   else if (nbP3d != 0)
 399     Ok = LineTool::Tangency(Line, index, TabV);
 400
 401   if (Ok) {
 402     if (nbP3d != 0) {
 403       j = 1;
 404       for (i = TabV.Lower(); i <= TabV.Upper(); i++) {
 405         V(j)   = TabV(i).X();
 406         V(j+1) = TabV(i).Y();
 407         V(j+2) = TabV(i).Z();
 408         j += 3;
 409       }
 410     }
 411     if (nbP2d != 0) {
 412       j = nbP3d*3+1;
 413       for (i = TabV2d.Lower(); i <= TabV2d.Upper(); i++) {
 414         V(j)   = TabV2d(i).X();
 415         V(j+1) = TabV2d(i).Y();
 416         j += 2;
 417       }
 418     }
 419   }
 420   else {
 421
 422     // recherche d un vecteur tangent par construction d une parabole:
 423     AppParCurves_Constraint firstC, lastC;
 424     firstC = lastC = AppParCurves_PassPoint;
 425     Standard_Integer nbpoles = 3;
 426     math_Vector mypar(index, index+2);
 427     Parameters(Line, index, index+2, mypar);
 428     Approx_ParLeastSquareOfMyGradient
 429       LSQ(Line, index, index+2, firstC, lastC, mypar, nbpoles);
 430     AppParCurves_MultiCurve C = LSQ.BezierValue();
 431
 432     gp_Pnt myP;
 433     gp_Vec myV;
 434     gp_Pnt2d myP2d;
 435     gp_Vec2d myV2d;
 436     j = 1;
 437     for (i = 1; i <= nbP3d; i++) {
 438       C.D1(i, 0.0, myP, myV);
 439       V(j)   = myV.X();
 440       V(j+1) = myV.Y();
 441       V(j+2) = myV.Z();
 442       j += 3;
 443     }
 444     j = nbP3d*3+1;
 445     for (i = nbP3d+1; i <= nbP3d+nbP2d; i++) {
 446       C.D1(i, 0.0, myP2d, myV2d);
 447       V(j)   = myV2d.X();
 448       V(j+1) = myV2d.Y();
 449       j += 2;
 450     }
 451   }
 452 }
 453
 454
 455 void Approx_ComputeLine::LastTangencyVector(const MultiLine&       Line,
 456                                             const Standard_Integer index,
 457                                             math_Vector&           V) const
 458 {
 459   Standard_Integer i, j, nbP2d, nbP3d;
 460   nbP3d = LineTool::NbP3d(Line);
 461   nbP2d = LineTool::NbP2d(Line);
 462   Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
 463   if (nbP3d == 0) mynbP3d = 1;
 464   if (nbP2d == 0) mynbP2d = 1;
 465   Standard_Boolean Ok=Standard_False;
 466   TColgp_Array1OfVec TabV(1, mynbP3d);
 467   TColgp_Array1OfVec2d TabV2d(1, mynbP2d);
 468
 469
 470   if (nbP3d != 0 && nbP2d != 0)
 471     Ok = LineTool::Tangency(Line, index, TabV, TabV2d);
 472   else if (nbP2d != 0)
 473     Ok = LineTool::Tangency(Line, index, TabV2d);
 474   else if (nbP3d != 0)
 475     Ok = LineTool::Tangency(Line, index, TabV);
 476
 477   if (Ok) {
 478     if (nbP3d != 0) {
 479       j = 1;
 480       for (i = TabV.Lower(); i <= TabV.Upper(); i++) {
 481         V(j)   = TabV(i).X();
 482         V(j+1) = TabV(i).Y();
 483         V(j+2) = TabV(i).Z();
 484         j += 3;
 485       }
 486     }
 487     if (nbP2d != 0) {
 488       j = nbP3d*3+1;
 489       for (i = TabV2d.Lower(); i <= TabV2d.Upper(); i++) {
 490         V(j)   = TabV2d(i).X();
 491         V(j+1) = TabV2d(i).Y();
 492         j += 2;
 493       }
 494     }
 495   }
 496   else {
 497
 498     // recherche d un vecteur tangent par construction d une parabole:
 499     AppParCurves_Constraint firstC, lastC;
 500     firstC = lastC = AppParCurves_PassPoint;
 501     Standard_Integer nbpoles = 3;
 502     math_Vector mypar(index-2, index);
 503     Parameters(Line, index-2, index, mypar);
 504     Approx_ParLeastSquareOfMyGradient
 505       LSQ(Line, index-2, index, firstC, lastC, mypar, nbpoles);
 506     AppParCurves_MultiCurve C = LSQ.BezierValue();
 507
 508     gp_Pnt myP;
 509     gp_Vec myV;
 510     gp_Pnt2d myP2d;
 511     gp_Vec2d myV2d;
 512     j = 1;
 513     for (i = 1; i <= nbP3d; i++) {
 514       C.D1(i, 1.0, myP, myV);
 515       V(j)   = myV.X();
 516       V(j+1) = myV.Y();
 517       V(j+2) = myV.Z();
 518       j += 3;
 519     }
 520     j = nbP3d*3+1;
 521     for (i = nbP3d+1; i <= nbP3d+nbP2d; i++) {
 522       C.D1(i, 1.0, myP2d, myV2d);
 523       V(j)   = myV2d.X();
 524       V(j+1) = myV2d.Y();
 525       j += 2;
 526     }
 527   }
 528
 529 }
 530
 531
 532
 533 Standard_Real Approx_ComputeLine::
 534   SearchFirstLambda(const MultiLine&            Line,
 535                     const math_Vector&          TheParam,
 536                     const math_Vector&          V,
 537                     const Standard_Integer      index) const
 538 {
 539
 540   // dq/dw = lambda* V = (p2-p1)/(u2-u1)
 541
 542   Standard_Integer nbP2d, nbP3d;
 543   gp_Pnt P1, P2;
 544   gp_Pnt2d P12d, P22d;
 545   nbP3d = LineTool::NbP3d(Line);
 546   nbP2d = LineTool::NbP2d(Line);
 547   Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
 548   if (nbP3d == 0) mynbP3d = 1;
 549   if (nbP2d == 0) mynbP2d = 1;
 550   TColgp_Array1OfPnt tabP1(1, mynbP3d), tabP2(1, mynbP3d);
 551   TColgp_Array1OfPnt2d tabP12d(1, mynbP2d), tabP22d(1, mynbP2d);
 552
 553
 554   if (nbP3d != 0 && nbP2d != 0) LineTool::Value(Line, index, tabP1, tabP12d);
 555   else if (nbP2d != 0)          LineTool::Value(Line, index, tabP12d);
 556   else if (nbP3d != 0)          LineTool::Value(Line, index, tabP1);
 557
 558   if (nbP3d != 0 && nbP2d != 0) LineTool::Value(Line, index+1, tabP2, tabP22d);
 559   else if (nbP2d != 0)          LineTool::Value(Line, index+1, tabP22d);
 560   else if (nbP3d != 0)          LineTool::Value(Line, index+1, tabP2);
 561
 562
 563   Standard_Real U1 = TheParam(index), U2 = TheParam(index+1);
 564   Standard_Real lambda, S;
 565   Standard_Integer low = V.Lower();
 566
 567   if (nbP3d != 0) {
 568     P1 = tabP1(1);
 569     P2 = tabP2(1);
 570     gp_Vec P1P2(P1, P2), myV;
 571     myV.SetCoord(V(low), V(low+1), V(low+2));
 572     lambda = (P1P2.Magnitude())/(myV.Magnitude()*(U2-U1));
 573     S = (P1P2.Dot(myV)> 0.0) ? 1.0 : -1.0;
 574   }
 575   else {
 576     P12d = tabP12d(1);
 577     P22d = tabP22d(1);
 578     gp_Vec2d P1P2(P12d, P22d), myV;
 579     myV.SetCoord(V(low), V(low+1));
 580     lambda = (P1P2.Magnitude())/(myV.Magnitude()*(U2-U1));
 581     S = (P1P2.Dot(myV)> 0.0) ? 1.0 : -1.0;
 582   }
 583   return (S*lambda);
 584
 585 }
 586
 587
 588 Standard_Real Approx_ComputeLine::
 589  SearchLastLambda(const MultiLine&            Line,
 590                   const math_Vector&          TheParam,
 591                   const math_Vector&          V,
 592                   const Standard_Integer      index) const
 593 {
 594   // dq/dw = lambda* V = (p2-p1)/(u2-u1)
 595
 596   Standard_Integer nbP2d, nbP3d;
 597   gp_Pnt P1, P2;
 598   gp_Pnt2d P12d, P22d;
 599   nbP3d = LineTool::NbP3d(Line);
 600   nbP2d = LineTool::NbP2d(Line);
 601   Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
 602   if (nbP3d == 0) mynbP3d = 1;
 603   if (nbP2d == 0) mynbP2d = 1;
 604   TColgp_Array1OfPnt tabP(1, mynbP3d), tabP2(1, mynbP3d);
 605   TColgp_Array1OfPnt2d tabP2d(1, mynbP2d), tabP22d(1, mynbP2d);
 606
 607   if (nbP3d != 0 && nbP2d != 0) LineTool::Value(Line, index-1, tabP, tabP2d);
 608   else if (nbP2d != 0)          LineTool::Value(Line, index-1, tabP2d);
 609   else if (nbP3d != 0)          LineTool::Value(Line, index-1, tabP);
 610
 611   if (nbP3d != 0 && nbP2d != 0) LineTool::Value(Line, index, tabP2, tabP22d);
 612   else if (nbP2d != 0)          LineTool::Value(Line, index, tabP22d);
 613   else if (nbP3d != 0)          LineTool::Value(Line, index, tabP2);
 614
 615
 616   Standard_Real U1 = TheParam(index-1), U2 = TheParam(index);
 617   Standard_Real lambda, S;
 618   Standard_Integer low = V.Lower();
 619
 620   if (nbP3d != 0) {
 621     P1 = tabP(1);
 622     P2 = tabP2(1);
 623     gp_Vec P1P2(P1, P2), myV;
 624     myV.SetCoord(V(low), V(low+1), V(low+2));
 625     lambda = (P1P2.Magnitude())/(myV.Magnitude()*(U2-U1));
 626     S = (P1P2.Dot(myV)> 0.0) ? 1.0 : -1.0;
 627   }
 628   else {
 629     P12d = tabP2d(1);
 630     P22d = tabP22d(1);
 631     gp_Vec2d P1P2(P12d, P22d), myV;
 632     myV.SetCoord(V(low), V(low+1));
 633     lambda = (P1P2.Magnitude())/(myV.Magnitude()*(U2-U1));
 634     S = (P1P2.Dot(myV)> 0.0) ? 1.0 : -1.0;
 635   }
 636
 637   return (S*lambda);
 638 }
 639
 640
 641
 642 Approx_ComputeLine::Approx_ComputeLine
 643                     (const MultiLine& Line,
 644                      const math_Vector& Parameters,
 645                      const Standard_Integer degreemin,
 646                      const Standard_Integer degreemax,
 647                      const Standard_Real Tolerance3d,
 648                      const Standard_Real Tolerance2d,
 649                      const Standard_Integer NbIterations,
 650                      const Standard_Boolean cutting,
 651                      const Standard_Boolean Squares)
 652 : myMultiLineNb (0),
 653   myNbPlusOnePoint (0),
 654   myIsClear (Standard_False)
 655 {
 656   myfirstParam = new TColStd_HArray1OfReal(Parameters.Lower(),
 657                                            Parameters.Upper());
 658   for (Standard_Integer i = Parameters.Lower(); i <= Parameters.Upper(); i++) {
 659     myfirstParam->SetValue(i, Parameters(i));
 660   }
 661   myConstraints = new AppParCurves_HArray1OfConstraintCouple(1, 2);
 662   Par = Approx_IsoParametric;
 663   mydegremin = degreemin;
 664   mydegremax = degreemax;
 665   mytol3d = Tolerance3d;
 666   mytol2d = Tolerance2d;
 667   mysquares = Squares;
 668   mycut = cutting;
 669   myitermax = NbIterations;
 670   alldone = Standard_False;
 671   myfirstC = AppParCurves_TangencyPoint;
 672   mylastC = AppParCurves_TangencyPoint;
 673   Perform(Line);
 674 }
 675
 676
 677 Approx_ComputeLine::Approx_ComputeLine
 678                     (const math_Vector& Parameters,
 679                      const Standard_Integer degreemin,
 680                      const Standard_Integer degreemax,
 681                      const Standard_Real Tolerance3d,
 682                      const Standard_Real Tolerance2d,
 683                      const Standard_Integer NbIterations,
 684                      const Standard_Boolean cutting,
 685                      const Standard_Boolean Squares)
 686 : myMultiLineNb (0),
 687   myNbPlusOnePoint (0),
 688   myIsClear (Standard_False)
 689 {
 690   myfirstParam = new TColStd_HArray1OfReal(Parameters.Lower(),
 691                                            Parameters.Upper());
 692   for (Standard_Integer i = Parameters.Lower(); i <= Parameters.Upper(); i++) {
 693     myfirstParam->SetValue(i, Parameters(i));
 694   }
 695   myfirstC = AppParCurves_TangencyPoint;
 696   mylastC = AppParCurves_TangencyPoint;
 697   myConstraints = new AppParCurves_HArray1OfConstraintCouple(1, 2);
 698   Par = Approx_IsoParametric;
 699   mydegremin = degreemin;
 700   mydegremax = degreemax;
 701   mytol3d = Tolerance3d;
 702   mytol2d = Tolerance2d;
 703   mysquares = Squares;
 704   mycut = cutting;
 705   myitermax = NbIterations;
 706   alldone = Standard_False;
 707 }
 708
 709 Approx_ComputeLine::Approx_ComputeLine
 710                     (const Standard_Integer degreemin,
 711                      const Standard_Integer degreemax,
 712                      const Standard_Real Tolerance3d,
 713                      const Standard_Real Tolerance2d,
 714                      const Standard_Integer NbIterations,
 715                      const Standard_Boolean cutting,
 716                      const Approx_ParametrizationType parametrization,
 717                      const Standard_Boolean Squares)
 718 : myMultiLineNb (0),
 719   myNbPlusOnePoint (0),
 720   myIsClear (Standard_False)
 721 {
 722   myConstraints = new AppParCurves_HArray1OfConstraintCouple(1, 2);
 723   Par = parametrization;
 724   mydegremin = degreemin;
 725   mydegremax = degreemax;
 726   mytol3d = Tolerance3d;
 727   mytol2d = Tolerance2d;
 728   mysquares = Squares;
 729   mycut = cutting;
 730   myitermax = NbIterations;
 731   myfirstC = AppParCurves_TangencyPoint;
 732   mylastC = AppParCurves_TangencyPoint;
 733   alldone = Standard_False;
 734 }
 735
 736
 737 Approx_ComputeLine::Approx_ComputeLine
 738                     (const MultiLine& Line,
 739                      const Standard_Integer degreemin,
 740                      const Standard_Integer degreemax,
 741                      const Standard_Real Tolerance3d,
 742                      const Standard_Real Tolerance2d,
 743                      const Standard_Integer NbIterations,
 744                      const Standard_Boolean cutting,
 745                      const Approx_ParametrizationType parametrization,
 746                      const Standard_Boolean Squares)
 747 : myMultiLineNb (0),
 748   myNbPlusOnePoint (0),
 749   myIsClear (Standard_False)
 750 {
 751   myConstraints = new AppParCurves_HArray1OfConstraintCouple(1, 2);
 752   alldone = Standard_False;
 753   mydegremin = degreemin;
 754   mydegremax = degreemax;
 755   mytol3d = Tolerance3d;
 756   mytol2d = Tolerance2d;
 757   mysquares = Squares;
 758   mycut = cutting;
 759   myitermax = NbIterations;
 760   Par = parametrization;
 761   myfirstC = AppParCurves_TangencyPoint;
 762   mylastC = AppParCurves_TangencyPoint;
 763
 764   Perform(Line);
 765 }
 766
 767
 768
 769 void Approx_ComputeLine::Perform(const MultiLine& Line)
 770 {
 771 #ifdef OCCT_DEBUG
 772   if (mydebug) DUMP(Line);
 773 #endif
 774   if (!myIsClear)
 775   {
 776     myMultiCurves.Clear();
 777     myPar.Clear();
 778     Tolers3d.Clear();
 779     Tolers2d.Clear();
 780     myMultiLineNb = 0;
 781     //myNbPlusOnePoint = 0;
 782   }
 783   else myIsClear = Standard_False;
 784
 785   Standard_Integer i, nbp, Thefirstpt, Thelastpt, oldlastpt;
 786   Standard_Boolean Finish = Standard_False,
 787           begin = Standard_True, Ok = Standard_False,
 788           GoUp = Standard_False, Interpol;
 789   Standard_Real thetol3d, thetol2d;
 790   Approx_Status MyStatus;
 791 //  gp_Vec V13d, V23d;
 792 //  gp_Vec2d V2d;
 793   Thefirstpt = LineTool::FirstPoint(Line);
 794   Thelastpt  = LineTool::LastPoint(Line);
 795   Standard_Integer myfirstpt = Thefirstpt;
 796   Standard_Integer mylastpt = Thelastpt;
 797
 798   AppParCurves_ConstraintCouple myCouple1(myfirstpt, myfirstC);
 799   AppParCurves_ConstraintCouple myCouple2(mylastpt, mylastC);
 800   myConstraints->SetValue(1, myCouple1);
 801   myConstraints->SetValue(2, myCouple2);
 802
 803   math_Vector TheParam(Thefirstpt, Thelastpt);
 804
 805
 806   if (!mycut) {
 807     if(myfirstParam.IsNull()) {
 808       Parameters(Line, Thefirstpt, Thelastpt, TheParam);
 809     }
 810     else {
 811       for (i = myfirstParam->Lower(); i <= myfirstParam->Upper(); i++) {
 812         TheParam(i+Thefirstpt-1) = myfirstParam->Value(i);
 813       }
 814     }
 815     TheMultiCurve = AppParCurves_MultiCurve();
 816     MultiLine anOtherLine0;
 817     Standard_Boolean isOtherLine0Made = Standard_False;
 818     Standard_Integer indbad = 0;
 819     alldone = Compute(Line, myfirstpt, mylastpt, TheParam, thetol3d, thetol2d, indbad);
 820     if (indbad != 0)
 821     {
 822       isOtherLine0Made = LineTool::MakeMLOneMorePoint(Line, myfirstpt, mylastpt, indbad, anOtherLine0);
 823     }
 824     if (isOtherLine0Made)
 825     {
 826       myIsClear = Standard_True;
 827       //++myMultiLineNb;
 828       myNbPlusOnePoint++;
 829       Perform(anOtherLine0);
 830       alldone = Standard_True;
 831     }
 832     if(!alldone && TheMultiCurve.NbCurves() > 0) {
 833 #ifdef OCCT_DEBUG
 834       if (mydebug) DUMP(TheMultiCurve);
 835 #endif
 836       myMultiCurves.Append(TheMultiCurve);
 837       Tolers3d.Append(currenttol3d);
 838       Tolers2d.Append(currenttol2d);
 839       Standard_Integer mylen = mylastpt-myfirstpt+1;
 840       Standard_Integer myParLen = myParameters->Length();
 841       Standard_Integer aLen = (myParLen > mylen)? myParLen : mylen;
 842       Handle(TColStd_HArray1OfReal) ThePar =
 843         new TColStd_HArray1OfReal(myfirstpt, myfirstpt+aLen-1);
 844       for (i = 0; i < aLen; i++)
 845         ThePar->SetValue(myfirstpt+i, myParameters->Value(myParameters->Lower()+i));
 846       myPar.Append(ThePar);
 847     }
 848   }
 849   else {
 850     while (!Finish) {
 851       oldlastpt = mylastpt;
 852       // Gestion du decoupage de la multiline pour approximer:
 853       if(!begin) {
 854         if (!GoUp) {
 855           if (Ok) {
 856             // Calcul de la partie a approximer.
 857             myfirstpt = mylastpt;
 858             mylastpt  = Thelastpt;
 859             if (myfirstpt == Thelastpt) {
 860               Finish = Standard_True;
 861               alldone = Standard_True;
 862               return;
 863             }
 864           }
 865           else {
 866             nbp = mylastpt - myfirstpt + 1;
 867             MyStatus = LineTool::WhatStatus(Line, myfirstpt, mylastpt);
 868             if (MyStatus == Approx_NoPointsAdded && nbp <= mydegremax+1) {
 869               Interpol = ComputeCurve(Line, myfirstpt, mylastpt);
 870               if (Interpol) {
 871                 if (mylastpt == Thelastpt) {
 872                   Finish = Standard_True;
 873                   alldone = Standard_True;
 874                   return;
 875                 }
 876               }
 877             }
 878             mylastpt = Standard_Integer((myfirstpt + mylastpt)/2);
 879           }
 880         }
 881         GoUp = Standard_False;
 882       }
 883
 884       // Verification du nombre de points restants par rapport au degre
 885       // demande.
 886       // ==============================================================
 887       nbp = mylastpt - myfirstpt + 1;
 888       MyStatus = LineTool::WhatStatus(Line, myfirstpt, mylastpt);
 889       if (nbp <= mydegremax+5 ) {
 890         // Rajout necessaire de points si possible.
 891         // ========================================
 892         GoUp = Standard_False;
 893         Ok = Standard_True;
 894         if (MyStatus == Approx_PointsAdded) {
 895           // Appel recursif du decoupage:
 896           GoUp = Standard_True;
 897
 898           MultiLine anOtherLine1 = LineTool::MakeMLBetween(Line, myfirstpt, mylastpt, nbp-1);
 899
 900           Standard_Integer nbpdsotherligne = LineTool::FirstPoint (anOtherLine1) - LineTool::LastPoint (anOtherLine1);
 901
 902           //-- Si MakeML a echoue   on retourne une ligne vide
 903           if ((nbpdsotherligne == 0) || myMultiLineNb >= 3)
 904           {
 905             //-- cout<<" ** ApproxComputeLine MakeML Echec ** LBR lbr "<<endl;
 906             if (myfirstpt == mylastpt) break;  // Pour etre sur de ne pas
 907             // planter la station !!
 908             myCouple1.SetIndex(myfirstpt);
 909             myCouple2.SetIndex(mylastpt);
 910             myConstraints->SetValue(1, myCouple1);
 911             myConstraints->SetValue(2, myCouple2);
 912
 913             math_Vector Param(myfirstpt, mylastpt);
 914             Approx_ParametrizationType SavePar = Par;
 915             Par = Approx_IsoParametric;
 916             Parameters(Line, myfirstpt, mylastpt, Param);
 917             TheMultiCurve = AppParCurves_MultiCurve();
 918             MultiLine anOtherLine2;
 919             Standard_Boolean isOtherLine2Made = Standard_False;
 920             Standard_Integer indbad = 0;
 921             Ok = Compute(Line, myfirstpt, mylastpt, Param, thetol3d, thetol2d, indbad);
 922             if (indbad != 0)
 923             {
 924               isOtherLine2Made = LineTool::MakeMLOneMorePoint(Line, myfirstpt, mylastpt, indbad, anOtherLine2);
 925             }
 926             if (isOtherLine2Made)
 927             {
 928               myIsClear = Standard_True;
 929               //++myMultiLineNb;
 930               myNbPlusOnePoint++;
 931               Par = SavePar;
 932               Perform(anOtherLine2);
 933               Ok = Standard_True;
 934             }
 935
 936             if (!Ok) {
 937               Standard_Real tt3d = currenttol3d, tt2d = currenttol2d;
 938               Handle(TColStd_HArray1OfReal) saveParameters = myParameters;
 939               AppParCurves_MultiCurve saveMultiCurve = TheMultiCurve;
 940
 941               if(SavePar != Approx_IsoParametric)
 942                 Par = SavePar;
 943               else
 944                 Par = Approx_ChordLength;
 945
 946               Parameters(Line, myfirstpt, mylastpt, Param);
 947         isOtherLine2Made = Standard_False;
 948               indbad = 0;
 949               Ok = Compute(Line, myfirstpt, mylastpt, Param, thetol3d, thetol2d, indbad);
 950               if (indbad != 0)
 951               {
 952                 isOtherLine2Made = LineTool::MakeMLOneMorePoint (Line, myfirstpt, mylastpt, indbad, anOtherLine2);
 953               }
 954               if (isOtherLine2Made)
 955               {
 956                 myIsClear = Standard_True;
 957                 //++myMultiLineNb;
 958                 myNbPlusOnePoint++;
 959                 Perform (anOtherLine2);
 960                 Ok = Standard_True;
 961               }
 962
 963               if (!Ok && tt3d <= currenttol3d && tt2d <= currenttol2d) {
 964                 currenttol3d = tt3d; currenttol2d = tt2d;
 965                 myParameters = saveParameters;
 966                 TheMultiCurve = saveMultiCurve;
 967               }
 968             }
 969             Par = SavePar;
 970             if (myfirstpt == Thelastpt)
 971             {
 972               Finish = Standard_True;
 973               alldone = Standard_True;
 974               return;
 975             }
 976
 977             oldlastpt = mylastpt;
 978             if (!Ok) {
 979               tolreached = Standard_False;
 980               if (TheMultiCurve.NbCurves() == 0) {
 981                 myMultiCurves.Clear();
 982                 return;
 983               }
 984 #ifdef OCCT_DEBUG
 985               if (mydebug) DUMP(TheMultiCurve);
 986 #endif
 987               MultiLine anOtherLine3;
 988               Standard_Boolean isOtherLine3Made = Standard_False;
 989               Standard_Integer indbad2 = 0;
 990               if (!CheckMultiCurve(TheMultiCurve, Line,
 991                                    myfirstpt, mylastpt,
 992                                    indbad2))
 993               {
 994                 isOtherLine3Made = LineTool::MakeMLOneMorePoint (Line, myfirstpt, mylastpt, indbad2, anOtherLine3);
 995               }
 996               if (isOtherLine3Made)
 997               {
 998                 myIsClear = Standard_True;
 999                 //++myMultiLineNb;
1000                 myNbPlusOnePoint++;
1001                 Perform(anOtherLine3);
1002                 myfirstpt = mylastpt;
1003                 mylastpt = Thelastpt;
1004               }
1005               else
1006               {
1007                 myMultiCurves.Append(TheMultiCurve);
1008                 Tolers3d.Append(currenttol3d);
1009                 Tolers2d.Append(currenttol2d);
1010                 Standard_Integer mylen = oldlastpt-myfirstpt+1;
1011                 Standard_Integer myParLen = myParameters->Length();
1012                 Standard_Integer aLen = (myParLen > mylen)? myParLen : mylen;
1013                 Handle(TColStd_HArray1OfReal) ThePar =
1014                   new TColStd_HArray1OfReal(myfirstpt, myfirstpt+aLen-1);
1015                 for (i = 0; i < aLen; i++)
1016                   ThePar->SetValue(myfirstpt+i, myParameters->Value(myParameters->Lower()+i));
1017                 myPar.Append(ThePar);
1018               }
1019             }
1020             myfirstpt = oldlastpt;
1021             mylastpt = Thelastpt;
1022
1023           }
1024           else
1025           {
1026             myIsClear = Standard_True;
1027             ++myMultiLineNb;
1028             Perform(anOtherLine1);
1029             myfirstpt = mylastpt;
1030             mylastpt = Thelastpt;
1031           }
1032         }
1033
1034         if  (MyStatus == Approx_NoPointsAdded && !begin) {
1035           // On rend la meilleure approximation obtenue precedemment.
1036           // ========================================================
1037           GoUp = Standard_True;
1038           tolreached = Standard_False;
1039           if (TheMultiCurve.NbCurves() == 0) {
1040             myMultiCurves.Clear();
1041             return;
1042           }
1043 #ifdef OCCT_DEBUG
1044       if (mydebug) DUMP(TheMultiCurve);
1045 #endif
1046           myMultiCurves.Append(TheMultiCurve);
1047           Tolers3d.Append(currenttol3d);
1048           Tolers2d.Append(currenttol2d);
1049           Standard_Integer mylen = oldlastpt-myfirstpt+1;
1050           Standard_Integer myParLen = myParameters->Length();
1051           Standard_Integer aLen = (myParLen > mylen)? myParLen : mylen;
1052           Handle(TColStd_HArray1OfReal) ThePar =
1053             new TColStd_HArray1OfReal(myfirstpt, myfirstpt+aLen-1);
1054           for (i = 0; i < aLen; i++)
1055             ThePar->SetValue(myfirstpt+i, myParameters->Value(myParameters->Lower()+i));
1056           myPar.Append(ThePar);
1057
1058           myfirstpt = oldlastpt;
1059           mylastpt = Thelastpt;
1060         }
1061
1062         else if (MyStatus == Approx_NoApproximation) {
1063           // On ne fait pas d approximation entre myfirstpt et mylastpt.
1064           // ===========================================================
1065           // On stocke pour pouvoir en informer l utilisateur.
1066           GoUp = Standard_True;
1067           myfirstpt = mylastpt;
1068           mylastpt = Thelastpt;
1069         }
1070       }
1071
1072       if (myfirstpt == Thelastpt) {
1073         Finish = Standard_True;
1074         alldone = Standard_True;
1075         return;
1076       }
1077       if (!GoUp) {
1078         if (myfirstpt == mylastpt) break;  // Pour etre sur de ne pas
1079                                            // planter la station !!
1080         myCouple1.SetIndex(myfirstpt);
1081         myCouple2.SetIndex(mylastpt);
1082         myConstraints->SetValue(1, myCouple1);
1083         myConstraints->SetValue(2, myCouple2);
1084
1085         // Calcul des parametres sur ce nouvel intervalle.
1086         // On recupere les parametres initiaux lors du decoupage.
1087
1088         math_Vector Param(myfirstpt, mylastpt);
1089         if (begin) {
1090           if(myfirstParam.IsNull()) {
1091             Parameters(Line, myfirstpt, mylastpt, Param);
1092           }
1093           else {
1094             for (i = myfirstParam->Lower(); i <= myfirstParam->Upper(); i++) {
1095               Param(i) = myfirstParam->Value(i);
1096             }
1097             myfirstParam.Nullify();
1098           }
1099           TheParam = Param;
1100           begin = Standard_False;
1101         }
1102         else {
1103           Standard_Real pfirst = TheParam.Value(myfirstpt);
1104           Standard_Real plast = TheParam.Value(mylastpt);
1105           for (i = myfirstpt; i <= mylastpt; i++) {
1106             Param(i) = (TheParam.Value(i)-pfirst)/(plast-pfirst);
1107           }
1108         }
1109
1110         TheMultiCurve = AppParCurves_MultiCurve();
1111         MultiLine anOtherLine4;
1112         Standard_Boolean isOtherLine4Made = Standard_False;
1113         Standard_Integer indbad = 0;
1114         Ok = Compute(Line, myfirstpt, mylastpt, Param, thetol3d, thetol2d, indbad);
1115         if (indbad != 0)
1116         {
1117           isOtherLine4Made = LineTool::MakeMLOneMorePoint (Line, myfirstpt, mylastpt, indbad, anOtherLine4);
1118         }
1119         if (isOtherLine4Made)
1120         {
1121           myIsClear = Standard_True;
1122           //++myMultiLineNb;
1123           myNbPlusOnePoint++;
1124           Perform (anOtherLine4);
1125           Ok = Standard_True;
1126         }
1127         if (myfirstpt == Thelastpt)
1128         {
1129           Finish = Standard_True;
1130           alldone = Standard_True;
1131           return;
1132         }
1133       }
1134     }
1135   }
1136 }
1137
1138
1139
1140 const TColStd_Array1OfReal& Approx_ComputeLine::Parameters(const Standard_Integer Index) const
1141 {
1142   return (myPar.Value(Index))->Array1();
1143 }
1144
1145
1146 Standard_Integer Approx_ComputeLine::NbMultiCurves()const
1147 {
1148   return myMultiCurves.Length();
1149 }
1150
1151 AppParCurves_MultiCurve& Approx_ComputeLine::ChangeValue(const Standard_Integer Index)
1152 {
1153   return myMultiCurves.ChangeValue(Index);
1154 }
1155
1156
1157 const AppParCurves_MultiCurve& Approx_ComputeLine::Value(const Standard_Integer Index)
1158 const
1159 {
1160   return myMultiCurves.Value(Index);
1161 }
1162
1163
1164 const AppParCurves_MultiBSpCurve& Approx_ComputeLine::SplineValue()
1165 {
1166   Approx_MCurvesToBSpCurve Trans;
1167   Trans.Perform(myMultiCurves);
1168   myspline = Trans.Value();
1169   return myspline;
1170 }
1171
1172
1173
1174
1175
1176 void Approx_ComputeLine::Parameters(const MultiLine& Line,
1177                                const Standard_Integer firstP,
1178                                const Standard_Integer lastP,
1179                                math_Vector& TheParameters) const
1180 {
1181   Standard_Integer i, j, nbP2d, nbP3d;
1182   Standard_Real dist;
1183   gp_Pnt P1, P2;
1184   gp_Pnt2d P12d, P22d;
1185
1186   if (Par == Approx_ChordLength || Par == Approx_Centripetal) {
1187     nbP3d = LineTool::NbP3d(Line);
1188     nbP2d = LineTool::NbP2d(Line);
1189     Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
1190     if (nbP3d == 0) mynbP3d = 1;
1191     if (nbP2d == 0) mynbP2d = 1;
1192
1193     TheParameters(firstP) = 0.0;
1194     dist = 0.0;
1195     TColgp_Array1OfPnt tabP(1, mynbP3d);
1196     TColgp_Array1OfPnt tabPP(1, mynbP3d);
1197     TColgp_Array1OfPnt2d tabP2d(1, mynbP2d);
1198     TColgp_Array1OfPnt2d tabPP2d(1, mynbP2d);
1199
1200     for (i = firstP+1; i <= lastP; i++) {
1201       if (nbP3d != 0 && nbP2d != 0) LineTool::Value(Line, i-1, tabP, tabP2d);
1202       else if (nbP2d != 0)          LineTool::Value(Line, i-1, tabP2d);
1203       else if (nbP3d != 0)          LineTool::Value(Line, i-1, tabP);
1204
1205       if (nbP3d != 0 && nbP2d != 0) LineTool::Value(Line, i, tabPP, tabPP2d);
1206       else if (nbP2d != 0)          LineTool::Value(Line, i, tabPP2d);
1207       else if (nbP3d != 0)          LineTool::Value(Line, i, tabPP);
1208       dist = 0;
1209       for (j = 1; j <= nbP3d; j++) {
1210         P1 = tabP(j);
1211         P2 = tabPP(j);
1212         dist += P2.Distance(P1);
1213       }
1214       for (j = 1; j <= nbP2d; j++) {
1215         P12d = tabP2d(j);
1216         P22d = tabPP2d(j);
1217         dist += P22d.Distance(P12d);
1218       }
1219       if(Par == Approx_ChordLength)
1220         TheParameters(i) = TheParameters(i-1) + dist;
1221       else {// Par == Approx_Centripetal
1222         TheParameters(i) = TheParameters(i-1) + Sqrt(dist);
1223       }
1224     }
1225     for (i = firstP; i <= lastP; i++) TheParameters(i) /= TheParameters(lastP);
1226   }
1227   else {
1228     for (i = firstP; i <= lastP; i++) {
1229       TheParameters(i) = (Standard_Real(i)-firstP)/
1230                          (Standard_Real(lastP)-Standard_Real(firstP));
1231     }
1232   }
1233 }
1234
1235
1236 Standard_Boolean Approx_ComputeLine::Compute(const MultiLine& Line,
1237                                              const Standard_Integer fpt,
1238                                              const Standard_Integer lpt,
1239                                              math_Vector&     Para,
1240                                              Standard_Real&   TheTol3d,
1241                                              Standard_Real&   TheTol2d,
1242                                              Standard_Integer& indbad)
1243 {
1244   indbad = 0;
1245   Standard_Integer deg, i;
1246   Standard_Boolean mydone;
1247   Standard_Real Fv;
1248   Standard_Integer nbp = lpt-fpt+1;
1249
1250   math_Vector ParSav(Para.Lower(), Para.Upper());
1251   for (i = Para.Lower(); i <= Para.Upper(); i++) {
1252     ParSav(i) = Para(i);
1253   }
1254   Standard_Integer Mdegmax = mydegremax;
1255   if(nbp < Mdegmax+5 && mycut) {
1256     Mdegmax = nbp - 5;
1257   }
1258   if(Mdegmax < mydegremin)  {
1259     Mdegmax = mydegremin;
1260   }
1261
1262   currenttol3d = currenttol2d = RealLast();
1263   for (deg = Min(nbp-1,mydegremin); deg <= Mdegmax; deg++) {
1264     AppParCurves_MultiCurve mySCU(deg+1);
1265     if (mysquares) {
1266       Approx_ParLeastSquareOfMyGradient SQ(Line, fpt, lpt,
1267                                            myfirstC, mylastC, Para, deg+1);
1268       mydone = SQ.IsDone();
1269       mySCU = SQ.BezierValue();
1270       SQ.Error(Fv, TheTol3d, TheTol2d);
1271     }
1272     else {
1273       Approx_MyGradient GRAD(Line, fpt, lpt, myConstraints,
1274                              Para, deg, mytol3d, mytol2d, myitermax);
1275       mydone = GRAD.IsDone();
1276       mySCU = GRAD.Value();
1277       if (mySCU.NbCurves() == 0)
1278         continue;
1279       TheTol3d = GRAD.MaxError3d();
1280       TheTol2d = GRAD.MaxError2d();
1281     }
1282     Standard_Real uu1 = Para(Para.Lower()), uu2;
1283     Standard_Boolean restau = Standard_False;
1284     for ( i = Para.Lower()+1; i <= Para.Upper(); i++) {
1285       uu2 =  Para(i);
1286       if (uu2 <= uu1) {
1287         restau = Standard_True;
1288 //      cout << "restau = Standard_True" << endl;
1289         break;
1290       }
1291       uu1 = uu2;
1292     }
1293     if (restau) {
1294       for (i = Para.Lower(); i <= Para.Upper(); i++) {
1295         Para(i) = ParSav(i);
1296       }
1297     }
1298     if (mydone) {
1299       if (TheTol3d <= mytol3d && TheTol2d <= mytol2d) {
1300         // Stockage de la multicurve approximee.
1301         tolreached = Standard_True;
1302 #ifdef OCCT_DEBUG
1303         if (mydebug) DUMP(mySCU);
1304 #endif
1305         if (myNbPlusOnePoint != 0 &&
1306             !CheckMultiCurve(mySCU, Line,
1307                              fpt, lpt,
1308                              indbad))
1309         {
1310           return Standard_False;
1311         }
1312         else
1313         {
1314           myMultiCurves.Append(mySCU);
1315           // Stockage des parametres de la partie de MultiLine approximee:
1316           // A ameliorer !! (bq trop de recopies)
1317           Handle(TColStd_HArray1OfReal) ThePar =
1318             new TColStd_HArray1OfReal(Para.Lower(), Para.Upper());
1319           for (i = Para.Lower(); i <= Para.Upper(); i++) {
1320             ThePar->SetValue(i, Para(i));
1321           }
1322           myPar.Append(ThePar);
1323           Tolers3d.Append(TheTol3d);
1324           Tolers2d.Append(TheTol2d);
1325           return Standard_True;
1326         }
1327       }
1328     }
1329
1330     if (TheTol3d <= currenttol3d && TheTol2d <= currenttol2d) {
1331       TheMultiCurve = mySCU;
1332       currenttol3d = TheTol3d;
1333       currenttol2d = TheTol2d;
1334       myParameters = new TColStd_HArray1OfReal(Para.Lower(), Para.Upper());
1335       for (i = Para.Lower(); i <= Para.Upper(); i++) {
1336         myParameters->SetValue(i, Para(i));
1337       }
1338     }
1339
1340   }
1341
1342   return Standard_False;
1343 }
1344
1345
1346
1347
1348 Standard_Boolean  Approx_ComputeLine::ComputeCurve(const MultiLine& Line,
1349                                       const Standard_Integer firstpt,
1350                                       const Standard_Integer lastpt)
1351 {
1352   Standard_Integer i, j, nbP3d, nbP2d, deg;
1353   gp_Vec V13d, V23d;
1354   gp_Vec2d V12d, V22d;
1355   gp_Pnt P1, P2;
1356   gp_Pnt2d P12d, P22d;
1357   Standard_Boolean Tangent1, Tangent2, mydone= Standard_False;
1358 #ifdef OCCT_DEBUG
1359   Standard_Boolean Parallel;
1360 #endif
1361   Standard_Integer myfirstpt = firstpt, mylastpt = lastpt;
1362   Standard_Integer nbp = lastpt-firstpt+1, Kopt = 0;
1363   math_Vector Para(firstpt, lastpt);
1364
1365   Parameters(Line, firstpt, lastpt, Para);
1366
1367   nbP3d = LineTool::NbP3d(Line);
1368   nbP2d = LineTool::NbP2d(Line);
1369   Standard_Integer mynbP3d = nbP3d, mynbP2d = nbP2d;
1370   if (nbP3d == 0) mynbP3d = 1 ;
1371   if (nbP2d == 0) mynbP2d = 1 ;
1372
1373
1374   TColgp_Array1OfVec tabV1(1, mynbP3d), tabV2(1, mynbP3d);
1375   TColgp_Array1OfPnt tabP1(1, mynbP3d), tabP2(1, mynbP3d), tabP(1, mynbP3d);
1376   TColgp_Array1OfVec2d tabV12d(1, mynbP2d), tabV22d(1, mynbP2d);
1377   TColgp_Array1OfPnt2d tabP12d(1, mynbP2d), tabP22d(1, mynbP2d), tabP2d(1, mynbP2d);
1378
1379   if (nbP3d != 0 && nbP2d != 0) {
1380     LineTool::Value(Line, myfirstpt,tabP1,tabP12d);
1381     LineTool::Value(Line, mylastpt,tabP2,tabP22d);
1382     Tangent1 = LineTool::Tangency(Line, myfirstpt,tabV1,tabV12d);
1383     Tangent2 = LineTool::Tangency(Line, mylastpt,tabV2,tabV22d);
1384   }
1385   else if (nbP2d != 0) {
1386     LineTool::Value(Line, myfirstpt,tabP12d);
1387     LineTool::Value(Line, mylastpt,tabP22d);
1388     Tangent1 = LineTool::Tangency(Line, myfirstpt, tabV12d);
1389     Tangent2 = LineTool::Tangency(Line, mylastpt, tabV22d);
1390   }
1391   else {
1392     LineTool::Value(Line, myfirstpt,tabP1);
1393     LineTool::Value(Line, mylastpt,tabP2);
1394     Tangent1 = LineTool::Tangency(Line, myfirstpt, tabV1);
1395     Tangent2 = LineTool::Tangency(Line, mylastpt, tabV2);
1396   }
1397
1398   if (Tangent1) Kopt++;
1399   if (Tangent2) Kopt++;
1400
1401
1402   if (nbp == 2) {
1403     // S il n y a que 2 points, on verifie quand meme que les tangentes sont
1404     // alignees.
1405 #ifdef OCCT_DEBUG
1406     Parallel = Standard_True;
1407 #endif
1408     if (Tangent1) {
1409       for (i = 1; i <= nbP3d; i++) {
1410         gp_Vec PVec(tabP1(i), tabP2(i));
1411         V13d = tabV1(i);
1412         if (!PVec.IsParallel(V13d, Precision::Angular())) {
1413 #ifdef OCCT_DEBUG
1414           Parallel = Standard_False;
1415 #endif
1416           break;
1417         }
1418       }
1419       for (i = 1; i <= nbP2d; i++) {
1420         gp_Vec2d PVec2d(tabP12d(i), tabP22d(i));
1421         V12d = tabV12d(i);
1422         if (!PVec2d.IsParallel(V12d, Precision::Angular())) {
1423 #ifdef OCCT_DEBUG
1424           Parallel = Standard_False;
1425 #endif
1426           break;
1427         }
1428       }
1429     }
1430
1431     if (Tangent2) {
1432       for (i = 1; i <= nbP3d; i++) {
1433         gp_Vec PVec(tabP1(i), tabP2(i));
1434         V23d = tabV2(i);
1435         if (!PVec.IsParallel(V23d, Precision::Angular())) {
1436 #ifdef OCCT_DEBUG
1437           Parallel = Standard_False;
1438 #endif
1439           break;
1440         }
1441       }
1442       for (i = 1; i <= nbP2d; i++) {
1443         gp_Vec2d PVec2d(tabP12d(i), tabP22d(i));
1444         V22d = tabV22d(i);
1445         if (!PVec2d.IsParallel(V22d, Precision::Angular())) {
1446 #ifdef OCCT_DEBUG
1447           Parallel = Standard_False;
1448 #endif
1449           break;
1450         }
1451       }
1452     }
1453
1454 #ifdef OCCT_DEBUG
1455     if (!Parallel) {
1456       if (mydebug) cout <<"droite mais tangentes pas vraiment paralleles!!"<< endl;
1457     }
1458 #endif
1459     AppParCurves_MultiCurve mySCU(mydegremin+1);
1460     if (nbP3d != 0 && nbP2d != 0) {
1461       AppParCurves_MultiPoint MPole1(tabP1, tabP12d);
1462       AppParCurves_MultiPoint MPole2(tabP2, tabP22d);
1463       mySCU.SetValue(1, MPole1);
1464       mySCU.SetValue(mydegremin+1, MPole2);
1465       for (i = 2; i <= mydegremin; i++) {
1466         for (j = 1; j<= nbP3d; j++) {
1467           P1 = tabP1(j);
1468           P2 = tabP2(j);
1469           tabP(j).SetXYZ(P1.XYZ()+(i-1)*(P2.XYZ()-P1.XYZ())/mydegremin);
1470         }
1471         for (j = 1; j<= nbP2d; j++) {
1472           P12d = tabP12d(j);
1473           P22d = tabP22d(j);
1474           tabP2d(j).SetXY(P12d.XY()+(i-1)*(P22d.XY()-P12d.XY())/mydegremin);
1475         }
1476         AppParCurves_MultiPoint MPole(tabP, tabP2d);
1477         mySCU.SetValue(i, MPole);
1478       }
1479
1480     }
1481     else if (nbP3d != 0) {
1482       AppParCurves_MultiPoint MPole1(tabP1);
1483       AppParCurves_MultiPoint MPole2(tabP2);
1484       mySCU.SetValue(1, MPole1);
1485       mySCU.SetValue(mydegremin+1, MPole2);
1486       for (i = 2; i <= mydegremin; i++) {
1487         for (j = 1; j<= nbP3d; j++) {
1488           P1 = tabP1(j);
1489           P2 = tabP2(j);
1490           tabP(j).SetXYZ(P1.XYZ()+(i-1)*(P2.XYZ()-P1.XYZ())/mydegremin);
1491         }
1492         AppParCurves_MultiPoint MPole(tabP);
1493         mySCU.SetValue(i, MPole);
1494       }
1495     }
1496     else if (nbP2d != 0) {
1497       AppParCurves_MultiPoint MPole1(tabP12d);
1498       AppParCurves_MultiPoint MPole2(tabP22d);
1499       mySCU.SetValue(1, MPole1);
1500       mySCU.SetValue(mydegremin+1, MPole2);
1501       for (i = 2; i <= mydegremin; i++) {
1502         for (j = 1; j<= nbP2d; j++) {
1503           P12d = tabP12d(j);
1504           P22d = tabP22d(j);
1505           tabP2d(j).SetXY(P12d.XY()+(i-1)*(P22d.XY()-P12d.XY())/mydegremin);
1506         }
1507         AppParCurves_MultiPoint MPole(tabP2d);
1508         mySCU.SetValue(i, MPole);
1509       }
1510     }
1511     mydone = Standard_True;
1512     // Stockage de la multicurve approximee.
1513     tolreached = Standard_True;
1514 #ifdef OCCT_DEBUG
1515       if (mydebug) DUMP(mySCU);
1516 #endif
1517     myMultiCurves.Append(mySCU);
1518     Handle(TColStd_HArray1OfReal) ThePar = new TColStd_HArray1OfReal(Para.Lower(), Para.Upper());
1519     for (i = Para.Lower(); i <= Para.Upper(); i++) {
1520       ThePar->SetValue(i, Para(i));
1521     }
1522     myPar.Append(ThePar);
1523     Tolers3d.Append(Precision::Confusion());
1524     Tolers2d.Append(Precision::PConfusion());
1525     return mydone;
1526   }
1527
1528   // avec les tangentes.
1529   deg = nbp+1;
1530   AppParCurves_MultiCurve mySCU(deg+1);
1531   AppParCurves_Constraint Cons = AppParCurves_TangencyPoint;
1532   Standard_Real lambda1, lambda2;
1533   math_Vector V1(1, nbP3d*3+nbP2d*2);
1534   math_Vector V2(1, nbP3d*3+nbP2d*2);
1535   FirstTangencyVector(Line, myfirstpt, V1);
1536   lambda1 = SearchFirstLambda(Line, Para, V1, myfirstpt);
1537
1538   LastTangencyVector(Line, mylastpt, V2);
1539   lambda2 = SearchLastLambda(Line, Para, V2, mylastpt);
1540
1541   Approx_ParLeastSquareOfMyGradient
1542     LSQ(Line, myfirstpt, mylastpt,
1543         Cons, Cons, Para, deg+1);
1544
1545   lambda1 = lambda1/deg;
1546   lambda2 = lambda2/deg;
1547   LSQ.Perform(Para, V1, V2, lambda1, lambda2);
1548   mydone = LSQ.IsDone();
1549   mySCU = LSQ.BezierValue();
1550
1551   if (mydone) {
1552     Standard_Real Fv, TheTol3d, TheTol2d;
1553     LSQ.Error(Fv, TheTol3d, TheTol2d);
1554
1555     // Stockage de la multicurve approximee.
1556     tolreached = Standard_True;
1557 #ifdef OCCT_DEBUG
1558       if (mydebug) DUMP(mySCU);
1559 #endif
1560     myMultiCurves.Append(mySCU);
1561     Handle(TColStd_HArray1OfReal) ThePar =
1562       new TColStd_HArray1OfReal(Para.Lower(), Para.Upper());
1563     for (i = Para.Lower(); i <= Para.Upper(); i++) {
1564       ThePar->SetValue(i, Para(i));
1565     }
1566     myPar.Append(ThePar);
1567     Tolers3d.Append(TheTol3d);
1568     Tolers2d.Append(TheTol2d);
1569     return Standard_True;
1570   }
1571   return mydone;
1572 }
1573
1574
1575
1576 void Approx_ComputeLine::Init(const Standard_Integer degreemin,
1577                               const Standard_Integer degreemax,
1578                               const Standard_Real Tolerance3d,
1579                               const Standard_Real Tolerance2d,
1580                               const Standard_Integer NbIterations,
1581                               const Standard_Boolean cutting,
1582                               const Approx_ParametrizationType parametrization,
1583                               const Standard_Boolean Squares)
1584 {
1585   mydegremin = degreemin;
1586   mydegremax = degreemax;
1587   mytol3d = Tolerance3d;
1588   mytol2d = Tolerance2d;
1589   Par = parametrization;
1590   mysquares = Squares;
1591   mycut = cutting;
1592   myitermax = NbIterations;
1593 }
1594
1595
1596
1597 void Approx_ComputeLine::SetDegrees(const Standard_Integer degreemin,
1598                                     const Standard_Integer degreemax)
1599 {
1600   mydegremin = degreemin;
1601   mydegremax = degreemax;
1602 }
1603
1604
1605 void Approx_ComputeLine::SetTolerances(const Standard_Real Tolerance3d,
1606                                        const Standard_Real Tolerance2d)
1607 {
1608   mytol3d = Tolerance3d;
1609   mytol2d = Tolerance2d;
1610 }
1611
1612
1613 void Approx_ComputeLine::SetConstraints(const AppParCurves_Constraint FirstC,
1614                                         const AppParCurves_Constraint LastC)
1615 {
1616   myfirstC = FirstC;
1617   mylastC = LastC;
1618 }
1619
1620
1621
1622 Standard_Boolean Approx_ComputeLine::IsAllApproximated()
1623 const {
1624   return alldone;
1625 }
1626
1627 Standard_Boolean Approx_ComputeLine::IsToleranceReached()
1628 const {
1629   return tolreached;
1630 }
1631
1632 void Approx_ComputeLine::Error(const Standard_Integer Index,
1633                                Standard_Real& tol3d,
1634                                Standard_Real& tol2d) const
1635 {
1636   tol3d = Tolers3d.Value(Index);
1637   tol2d = Tolers2d.Value(Index);
1638 }
1639
1640 Approx_ParametrizationType Approx_ComputeLine::Parametrization() const
1641 {
1642   return Par;
1643 }