src/GeomFill/GeomFill_Frenet.cxx

   1 // Created on: 1998-01-22
   2 // Created by: Philippe MANGIN/Roman BORISOV
   3 // Copyright (c) 1998-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17
  18 #include <Adaptor3d_HCurve.hxx>
  19 #include <Extrema_ExtPC.hxx>
  20 #include <GeomAbs_CurveType.hxx>
  21 #include <GeomFill_Frenet.hxx>
  22 #include <GeomFill_SnglrFunc.hxx>
  23 #include <GeomFill_TrihedronLaw.hxx>
  24 #include <GeomLib.hxx>
  25 #include <gp_Vec.hxx>
  26 #include <NCollection_Array1.hxx>
  27 #include <Precision.hxx>
  28 #include <Standard_ConstructionError.hxx>
  29 #include <Standard_OutOfRange.hxx>
  30 #include <Standard_Type.hxx>
  31 #include <TColgp_SequenceOfPnt2d.hxx>
  32 #include <TColStd_HArray1OfBoolean.hxx>
  33
  34 #include <algorithm>
  35 IMPLEMENT_STANDARD_RTTIEXT(GeomFill_Frenet,GeomFill_TrihedronLaw)
  36
  37 static const Standard_Real NullTol = 1.e-10;
  38 static const Standard_Real MaxSingular = 1.e-5;
  39
  40 static const Standard_Integer maxDerivOrder = 3;
  41
  42
  43 //=======================================================================
  44 //function : FDeriv
  45 //purpose  : computes (F/|F|)'
  46 //=======================================================================
  47 static gp_Vec FDeriv(const gp_Vec& F, const gp_Vec& DF)
  48 {
  49   Standard_Real Norma = F.Magnitude();
  50   gp_Vec Result = (DF - F*(F*DF)/(Norma*Norma))/Norma;
  51   return Result;
  52 }
  53
  54
  55 //=======================================================================
  56 //function : DDeriv
  57 //purpose  : computes (F/|F|)''
  58 //=======================================================================
  59 static gp_Vec DDeriv(const gp_Vec& F, const gp_Vec& DF, const gp_Vec& D2F)
  60 {
  61   Standard_Real Norma = F.Magnitude();
  62   gp_Vec Result = (D2F - 2*DF*(F*DF)/(Norma*Norma))/Norma -
  63      F*((DF.SquareMagnitude() + F*D2F
  64         - 3*(F*DF)*(F*DF)/(Norma*Norma))/(Norma*Norma*Norma));
  65   return Result;
  66 }
  67
  68 //=======================================================================
  69 //function : CosAngle
  70 //purpose  : Return a cosine between vectors theV1 and theV2.
  71 //=======================================================================
  72 static Standard_Real CosAngle(const gp_Vec& theV1, const gp_Vec& theV2)
  73   {
  74   const Standard_Real aTol = gp::Resolution();
  75
  76   const Standard_Real m1 = theV1.Magnitude(), m2 = theV2.Magnitude();
  77   if((m1 <= aTol) || (m2 <= aTol)) //Vectors are codirectional
  78     return 1.0;
  79
  80   Standard_Real aCAng = theV1.Dot(theV2)/(m1*m2);
  81
  82   if(aCAng > 1.0)
  83     aCAng = 1.0;
  84
  85   if(aCAng < -1.0)
  86     aCAng = -1.0;
  87
  88
  89   return aCAng;
  90   }
  91
  92 //=======================================================================
  93 //function : GeomFill_Frenet
  94 //purpose  :
  95 //=======================================================================
  96
  97 GeomFill_Frenet::GeomFill_Frenet()
  98 {
  99 }
 100
 101
 102 //=======================================================================
 103 //function : Copy
 104 //purpose  :
 105 //=======================================================================
 106
 107 Handle(GeomFill_TrihedronLaw) GeomFill_Frenet::Copy() const
 108 {
 109   Handle(GeomFill_Frenet) copy = new (GeomFill_Frenet)();
 110   if (!myCurve.IsNull()) copy->SetCurve(myCurve);
 111   return copy;
 112 }
 113
 114 //=======================================================================
 115 //function : SetCurve
 116 //purpose  :
 117 //=======================================================================
 118
 119  void GeomFill_Frenet::SetCurve(const Handle(Adaptor3d_HCurve)& C)
 120 {
 121   GeomFill_TrihedronLaw::SetCurve(C);
 122   if (! C.IsNull()) {
 123     GeomAbs_CurveType type;
 124     type = C->GetType();
 125     switch  (type) {
 126     case GeomAbs_Circle:
 127     case GeomAbs_Ellipse:
 128     case GeomAbs_Hyperbola:
 129     case GeomAbs_Parabola:
 130     case GeomAbs_Line:
 131       {
 132         // No probleme
 133         isSngl = Standard_False;
 134       }
 135      default :
 136        {
 137          // We have to search singulaties
 138          Init();
 139        }
 140     }
 141   }
 142 }
 143
 144 //=======================================================================
 145 //function : Init
 146 //purpose  :
 147 //=======================================================================
 148
 149  void GeomFill_Frenet::Init()
 150 {
 151   Standard_Integer i, j;
 152   GeomFill_SnglrFunc Func(myCurve);
 153   Standard_Real TolF = 1.0e-10, Tol = 10*TolF, Tol2 = Tol * Tol,
 154                 PTol = Precision::PConfusion();
 155
 156 // We want to determine if the curve has linear segments
 157   Standard_Integer NbIntC2 = myCurve->NbIntervals(GeomAbs_C2);
 158   Handle(TColStd_HArray1OfReal) myC2Disc =
 159     new TColStd_HArray1OfReal(1, NbIntC2 + 1);
 160   Handle(TColStd_HArray1OfBoolean) IsLin =
 161     new TColStd_HArray1OfBoolean(1, NbIntC2);
 162   Handle(TColStd_HArray1OfBoolean) IsConst =
 163     new TColStd_HArray1OfBoolean(1, NbIntC2);
 164   TColStd_Array1OfReal AveFunc(1,  NbIntC2);
 165   myCurve->Intervals(myC2Disc->ChangeArray1(), GeomAbs_C2);
 166   Standard_Integer NbControl = 10;
 167   Standard_Real Step, Average = 0, modulus;
 168   gp_Pnt C, C1;
 169   for(i = 1; i <= NbIntC2; i++) {
 170     Step = (myC2Disc->Value(i+1) - myC2Disc->Value(i))/NbControl;
 171     IsLin->ChangeValue(i) = Standard_True;
 172     IsConst->ChangeValue(i) = Standard_True;
 173     for(j = 1; j <= NbControl; j++) {
 174       Func.D0(myC2Disc->Value(i) + (j-1)*Step, C);
 175       if(j == 1) C1 = C;
 176       modulus = C.XYZ().Modulus();
 177       if(modulus > Tol) {
 178         IsLin->ChangeValue(i) = Standard_False;
 179       }
 180       Average += modulus;
 181
 182       if(IsConst->Value(i)) {
 183         if(Abs(C.X() - C1.X()) > Tol ||
 184            Abs(C.Y() - C1.Y()) > Tol ||
 185            Abs(C.Z() - C1.Z()) > Tol) {
 186           IsConst->ChangeValue(i) = Standard_False;
 187         }
 188       }
 189
 190     }
 191     AveFunc(i) = Average/NbControl;
 192   }
 193
 194 // Here we are looking for singularities
 195   TColStd_SequenceOfReal * SeqArray = new TColStd_SequenceOfReal[ NbIntC2 ];
 196   TColStd_SequenceOfReal SnglSeq;
 197 //  Standard_Real Value2, preValue=1.e200, t;
 198   Standard_Real Value2, t;
 199   Extrema_ExtPC Ext;
 200   gp_Pnt Origin(0, 0, 0);
 201
 202   for(i = 1; i <= NbIntC2; i++) {
 203     if (!IsLin->Value(i) && !IsConst->Value(i))
 204       {
 205         Func.SetRatio(1./AveFunc(i)); // Normalization
 206         Ext.Initialize(Func, myC2Disc->Value(i), myC2Disc->Value(i+1), TolF);
 207         Ext.Perform(Origin);
 208         if(Ext.IsDone() && Ext.NbExt() != 0)
 209           {
 210             for(j = 1; j <= Ext.NbExt(); j++)
 211               {
 212                 Value2 = Ext.SquareDistance(j);
 213                 if(Value2 < Tol2)
 214                   {
 215                     t = Ext.Point(j).Parameter();
 216                     SeqArray[i-1].Append(t);
 217                   }
 218               }
 219           }
 220         // sorting
 221         if(SeqArray[i-1].Length() != 0) {
 222           NCollection_Array1<Standard_Real> anArray( 1, SeqArray[i-1].Length() );
 223           for (j = 1; j <= anArray.Length(); j++)
 224             anArray(j) = SeqArray[i-1](j);
 225           std::sort (anArray.begin(), anArray.end());
 226           for (j = 1; j <= anArray.Length(); j++)
 227             SeqArray[i-1](j) = anArray(j);
 228         }
 229       }
 230   }
 231   //Filling SnglSeq by first sets of roots
 232   for(i = 0; i < NbIntC2; i++)
 233     for (j = 1; j <= SeqArray[i].Length(); j++)
 234       SnglSeq.Append( SeqArray[i](j) );
 235
 236   //Extrema works bad, need to pass second time
 237   for(i = 0; i < NbIntC2; i++)
 238     if (! SeqArray[i].IsEmpty())
 239       {
 240         SeqArray[i].Prepend( myC2Disc->Value(i+1) );
 241         SeqArray[i].Append( myC2Disc->Value(i+2) );
 242         Func.SetRatio(1./AveFunc(i+1)); // Normalization
 243         for (j = 1; j < SeqArray[i].Length(); j++)
 244           if (SeqArray[i](j+1) - SeqArray[i](j) > PTol)
 245             {
 246               Ext.Initialize(Func, SeqArray[i](j), SeqArray[i](j+1), TolF);
 247               Ext.Perform(Origin);
 248               if(Ext.IsDone())
 249                 {
 250                   for(Standard_Integer k = 1; k <= Ext.NbExt(); k++)
 251                     {
 252                       Value2 = Ext.SquareDistance(k);
 253                       if(Value2 < Tol2)
 254                         {
 255                           t = Ext.Point(k).Parameter();
 256                           if (t-SeqArray[i](j) > PTol && SeqArray[i](j+1)-t > PTol)
 257                             SnglSeq.Append(t);
 258                         }
 259                     }
 260                 }
 261             }
 262       }
 263
 264   delete [] SeqArray;
 265
 266   if(SnglSeq.Length() > 0) {
 267     // sorting
 268     NCollection_Array1<Standard_Real> anArray( 1, SnglSeq.Length() );
 269     for (i = 1; i <= SnglSeq.Length(); i++)
 270       anArray(i) = SnglSeq(i);
 271     std::sort (anArray.begin(), anArray.end());
 272     for (i = 1; i <= SnglSeq.Length(); i++)
 273       SnglSeq(i) = anArray(i);
 274
 275     // discard repeating elements
 276     Standard_Boolean found = Standard_True;
 277     j = 1;
 278     while (found)
 279       {
 280         found = Standard_False;
 281         for (i = j; i < SnglSeq.Length(); i++)
 282           if (SnglSeq(i+1) - SnglSeq(i) <= PTol)
 283             {
 284               SnglSeq.Remove(i+1);
 285               j = i;
 286               found = Standard_True;
 287               break;
 288             }
 289       }
 290
 291     mySngl = new TColStd_HArray1OfReal(1, SnglSeq.Length());
 292     for(i = 1; i <= mySngl->Length(); i++)
 293       mySngl->ChangeValue(i) = SnglSeq(i);
 294
 295 // computation of length of singular interval
 296     mySnglLen = new TColStd_HArray1OfReal(1, mySngl->Length());
 297     gp_Vec SnglDer, SnglDer2;
 298     Standard_Real norm;
 299     for(i = 1; i <= mySngl->Length(); i++) {
 300       Func.D2(mySngl->Value(i), C, SnglDer, SnglDer2);
 301       if ((norm = SnglDer.Magnitude()) > gp::Resolution())
 302         mySnglLen->ChangeValue(i) = Min(NullTol/norm, MaxSingular);
 303       else if ((norm = SnglDer2.Magnitude()) > gp::Resolution())
 304         mySnglLen->ChangeValue(i) = Min(Sqrt(2*NullTol/norm), MaxSingular);
 305       else mySnglLen->ChangeValue(i) = MaxSingular;
 306     }
 307 #ifdef OCCT_DEBUG
 308     for(i = 1; i <= mySngl->Length(); i++) {
 309       cout<<"Sngl("<<i<<") = "<<mySngl->Value(i)<<" Length = "<<mySnglLen->Value(i)<<endl;
 310     }
 311 #endif
 312     if(mySngl->Length() > 1) {
 313 // we have to merge singular points that have common parts of singular intervals
 314     TColgp_SequenceOfPnt2d tmpSeq;
 315     tmpSeq.Append(gp_Pnt2d(mySngl->Value(1), mySnglLen->Value(1)));
 316     Standard_Real U11, U12, U21, U22;
 317     for(i = 2; i<= mySngl->Length(); i++) {
 318       U12 = tmpSeq.Last().X() + tmpSeq.Last().Y();
 319       U21 = mySngl->Value(i) - mySnglLen->Value(i);
 320       if(U12 >= U21) {
 321         U11 = tmpSeq.Last().X() - tmpSeq.Last().Y();
 322         U22 = mySngl->Value(i) + mySnglLen->Value(i);
 323         tmpSeq.ChangeValue(tmpSeq.Length()) = gp_Pnt2d((U11 + U22)/2, (U22 - U11)/2);
 324       }
 325       else tmpSeq.Append(gp_Pnt2d(mySngl->Value(i), mySnglLen->Value(i)));
 326     }
 327     mySngl = new TColStd_HArray1OfReal(1, tmpSeq.Length());
 328     mySnglLen = new TColStd_HArray1OfReal(1, tmpSeq.Length());
 329     for(i = 1; i <= mySngl->Length(); i++) {
 330       mySngl->ChangeValue(i) = tmpSeq(i).X();
 331       mySnglLen->ChangeValue(i) = tmpSeq(i).Y();
 332     }
 333   }
 334 #ifdef OCCT_DEBUG
 335     cout<<"After merging"<<endl;
 336     for(i = 1; i <= mySngl->Length(); i++) {
 337       cout<<"Sngl("<<i<<") = "<<mySngl->Value(i)<<" Length = "<<mySnglLen->Value(i)<<endl;
 338     }
 339 #endif
 340     isSngl = Standard_True;
 341   }
 342   else isSngl = Standard_False;
 343 }
 344
 345 //=======================================================================
 346 //function :  RotateTrihedron
 347 //purpose  :  This function revolves the trihedron (which is determined of
 348 //            given "Tangent", "Normal" and "BiNormal" vectors)
 349 //            to coincide "Tangent" and "NewTangent" axes.
 350 //=======================================================================
 351 Standard_Boolean
 352           GeomFill_Frenet::RotateTrihedron( gp_Vec& Tangent,
 353                                             gp_Vec& Normal,
 354                                             gp_Vec& BiNormal,
 355                                             const gp_Vec& NewTangent) const
 356   {
 357   const Standard_Real anInfCOS = cos(Precision::Angular()); //0.99999995
 358   const Standard_Real aTol = gp::Resolution();
 359
 360   gp_Vec anAxis = Tangent.Crossed(NewTangent);
 361   const Standard_Real NT = anAxis.Magnitude();
 362   if(NT <= aTol)
 363     //No rotation required
 364     return Standard_True;
 365   else
 366     anAxis /= NT;   //Normalization
 367
 368   const Standard_Real aPx = anAxis.X(), aPy = anAxis.Y(), aPz = anAxis.Z();
 369   const Standard_Real aCAng = CosAngle(Tangent,NewTangent); //cosine
 370
 371   const Standard_Real anAddCAng = 1.0 - aCAng;
 372   const Standard_Real aSAng = sqrt(1.0 - aCAng*aCAng);  //sine
 373
 374 //According to rotate direction, sine of rotation angle might be
 375 //positive or negative.
 376 //We can research to choose necessary sign. But I think, it is more
 377 //effectively, to rotate "Tangent" vector in both direction. After that
 378 //we can choose necessary rotation direction in depend of results.
 379
 380   const gp_Vec aV11(anAddCAng*aPx*aPx+aCAng,      anAddCAng*aPx*aPy-aPz*aSAng,  anAddCAng*aPx*aPz+aPy*aSAng);
 381   const gp_Vec aV12(anAddCAng*aPx*aPx+aCAng,      anAddCAng*aPx*aPy+aPz*aSAng,  anAddCAng*aPx*aPz-aPy*aSAng);
 382   const gp_Vec aV21(anAddCAng*aPx*aPy+aPz*aSAng,  anAddCAng*aPy*aPy+aCAng,      anAddCAng*aPy*aPz-aPx*aSAng);
 383   const gp_Vec aV22(anAddCAng*aPx*aPy-aPz*aSAng,  anAddCAng*aPy*aPy+aCAng,      anAddCAng*aPy*aPz+aPx*aSAng);
 384   const gp_Vec aV31(anAddCAng*aPx*aPz-aPy*aSAng,  anAddCAng*aPy*aPz+aPx*aSAng,  anAddCAng*aPz*aPz+aCAng);
 385   const gp_Vec aV32(anAddCAng*aPx*aPz+aPy*aSAng,  anAddCAng*aPy*aPz-aPx*aSAng,  anAddCAng*aPz*aPz+aCAng);
 386
 387   gp_Vec aT1(Tangent.Dot(aV11), Tangent.Dot(aV21), Tangent.Dot(aV31));
 388   gp_Vec aT2(Tangent.Dot(aV12), Tangent.Dot(aV22), Tangent.Dot(aV32));
 389
 390   if(CosAngle(aT1,NewTangent) >= CosAngle(aT2,NewTangent))
 391     {
 392     Tangent = aT1;
 393     Normal = gp_Vec(Normal.Dot(aV11), Normal.Dot(aV21), Normal.Dot(aV31));
 394     BiNormal = gp_Vec(BiNormal.Dot(aV11), BiNormal.Dot(aV21), BiNormal.Dot(aV31));
 395     }
 396   else
 397     {
 398     Tangent = aT2;
 399     Normal = gp_Vec(Normal.Dot(aV12), Normal.Dot(aV22), Normal.Dot(aV32));
 400     BiNormal = gp_Vec(BiNormal.Dot(aV12), BiNormal.Dot(aV22), BiNormal.Dot(aV32));
 401     }
 402
 403   return CosAngle(Tangent,NewTangent) >= anInfCOS;
 404   }
 405
 406
 407 //=======================================================================
 408 //function : D0
 409 //purpose  :
 410 //=======================================================================
 411
 412  Standard_Boolean GeomFill_Frenet::D0(const Standard_Real theParam,
 413                                       gp_Vec& Tangent,
 414                                       gp_Vec& Normal,
 415                                       gp_Vec& BiNormal)
 416 {
 417   const Standard_Real aTol = gp::Resolution();
 418
 419   Standard_Real norm;
 420   Standard_Integer Index;
 421   Standard_Real Delta = 0.;
 422   if(IsSingular(theParam, Index))
 423     if (SingularD0(theParam, Index, Tangent, Normal, BiNormal, Delta))
 424       return Standard_True;
 425
 426   Standard_Real aParam = theParam + Delta;
 427   myTrimmed->D2(aParam, P, Tangent, BiNormal);
 428
 429   const Standard_Real DivisionFactor = 1.e-3;
 430   const Standard_Real anUinfium   = myTrimmed->FirstParameter();
 431   const Standard_Real anUsupremum = myTrimmed->LastParameter();
 432   const Standard_Real aDelta = (anUsupremum - anUinfium)*DivisionFactor;
 433   Standard_Real Ndu = Tangent.Magnitude();
 434
 435 //////////////////////////////////////////////////////////////////////////////////////////////////////
 436   if(Ndu <= aTol)
 437     {
 438     gp_Vec aTn;
 439 //Derivative is approximated by Taylor-series
 440
 441     Standard_Integer anIndex = 1; //Derivative order
 442     Standard_Boolean isDeriveFound = Standard_False;
 443
 444     do
 445       {
 446       aTn =  myTrimmed->DN(theParam,++anIndex);
 447       Ndu = aTn.Magnitude();
 448       isDeriveFound = Ndu > aTol;
 449       }
 450     while(!isDeriveFound && anIndex < maxDerivOrder);
 451
 452     if(isDeriveFound)
 453       {
 454       Standard_Real u;
 455
 456       if(theParam-anUinfium < aDelta)
 457         u = theParam+aDelta;
 458       else
 459         u = theParam-aDelta;
 460
 461       gp_Pnt P1, P2;
 462       myTrimmed->D0(Min(theParam, u),P1);
 463       myTrimmed->D0(Max(theParam, u),P2);
 464
 465       gp_Vec V1(P1,P2);
 466       Standard_Real aDirFactor = aTn.Dot(V1);
 467
 468       if(aDirFactor < 0.0)
 469         aTn = -aTn;
 470       }//if(IsDeriveFound)
 471     else
 472       {
 473 //Derivative is approximated by three points
 474
 475       gp_Pnt Ptemp(0.0,0.0,0.0); //(0,0,0)-coordinate
 476       gp_Pnt P1, P2, P3;
 477       Standard_Boolean IsParameterGrown;
 478
 479       if(theParam-anUinfium < 2*aDelta)
 480         {
 481         myTrimmed->D0(theParam,P1);
 482         myTrimmed->D0(theParam+aDelta,P2);
 483         myTrimmed->D0(theParam+2*aDelta,P3);
 484         IsParameterGrown = Standard_True;
 485         }
 486       else
 487         {
 488         myTrimmed->D0(theParam-2*aDelta,P1);
 489         myTrimmed->D0(theParam-aDelta,P2);
 490         myTrimmed->D0(theParam,P3);
 491         IsParameterGrown = Standard_False;
 492         }
 493
 494       gp_Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
 495
 496       if(IsParameterGrown)
 497         aTn=-3*V1+4*V2-V3;
 498       else
 499         aTn=V1-4*V2+3*V3;
 500       }//else of "if(IsDeriveFound)" condition
 501       Ndu = aTn.Magnitude();
 502       gp_Pnt Pt = P;
 503       Standard_Real dPar = 10.0*aDelta;
 504
 505 //Recursive calling is used for determine of trihedron for
 506 //point, which is near to given.
 507       if(theParam-anUinfium < dPar)
 508         {
 509         if(D0(aParam+dPar,Tangent,Normal,BiNormal) == Standard_False)
 510           return Standard_False;
 511         }
 512       else
 513         {
 514         if(D0(aParam-dPar,Tangent,Normal,BiNormal) == Standard_False)
 515           return Standard_False;
 516         }
 517
 518       P = Pt;
 519
 520       if(RotateTrihedron(Tangent,Normal,BiNormal,aTn) == Standard_False)
 521         {
 522 #ifdef OCCT_DEBUG
 523         cout << "Cannot coincide two tangents." << endl;
 524 #endif
 525         return Standard_False;
 526         }
 527     }//if(Ndu <= aTol)
 528   else
 529     {
 530     Tangent = Tangent/Ndu;
 531     BiNormal = Tangent.Crossed(BiNormal);
 532     norm = BiNormal.Magnitude();
 533     if (norm <= gp::Resolution())
 534       {
 535       gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent);
 536       BiNormal.SetXYZ(Axe.YDirection().XYZ());
 537       }
 538     else
 539       BiNormal.Normalize();
 540
 541     Normal = BiNormal;
 542     Normal.Cross(Tangent);
 543     }
 544
 545   return Standard_True;
 546 }
 547
 548 //=======================================================================
 549 //function : D1
 550 //purpose  :
 551 //=======================================================================
 552
 553  Standard_Boolean GeomFill_Frenet::D1(const Standard_Real Param,
 554                                       gp_Vec& Tangent,
 555                                       gp_Vec& DTangent,
 556                                       gp_Vec& Normal,
 557                                       gp_Vec& DNormal,
 558                                       gp_Vec& BiNormal,
 559                                       gp_Vec& DBiNormal)
 560 {
 561   Standard_Integer Index;
 562   Standard_Real Delta = 0.;
 563   if(IsSingular(Param, Index))
 564     if (SingularD1(Param, Index, Tangent, DTangent, Normal, DNormal, BiNormal, DBiNormal, Delta))
 565       return Standard_True;
 566
 567 //  Standard_Real Norma;
 568   Standard_Real theParam = Param + Delta;
 569   gp_Vec DC1, DC2, DC3;
 570   myTrimmed->D3(theParam, P, DC1, DC2, DC3);
 571   Tangent = DC1.Normalized();
 572
 573   //if (DC2.Magnitude() <= NullTol || Tangent.Crossed(DC2).Magnitude() <= NullTol) {
 574   if (Tangent.Crossed(DC2).Magnitude() <= gp::Resolution()) {
 575     gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent);
 576     Normal.SetXYZ(Axe.XDirection().XYZ());
 577     BiNormal.SetXYZ(Axe.YDirection().XYZ());
 578     DTangent.SetCoord(0,0,0);
 579     DNormal.SetCoord(0,0,0);
 580     DBiNormal.SetCoord(0,0,0);
 581     return Standard_True;
 582   }
 583   else
 584     BiNormal = Tangent.Crossed(DC2).Normalized();
 585
 586   Normal = BiNormal.Crossed(Tangent);
 587
 588   DTangent = FDeriv(DC1, DC2);
 589
 590   gp_Vec instead_DC1, instead_DC2;
 591   instead_DC1 = Tangent.Crossed(DC2);
 592   instead_DC2 = DTangent.Crossed(DC2) + Tangent.Crossed(DC3);
 593   DBiNormal = FDeriv(instead_DC1, instead_DC2);
 594
 595   DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
 596   return Standard_True;
 597 }
 598
 599 //=======================================================================
 600 //function : D2
 601 //purpose  :
 602 //=======================================================================
 603
 604  Standard_Boolean GeomFill_Frenet::D2(const Standard_Real Param,
 605                                       gp_Vec& Tangent,
 606                                       gp_Vec& DTangent,
 607                                       gp_Vec& D2Tangent,
 608                                       gp_Vec& Normal,
 609                                       gp_Vec& DNormal,
 610                                       gp_Vec& D2Normal,
 611                                       gp_Vec& BiNormal,
 612                                       gp_Vec& DBiNormal,
 613                                       gp_Vec& D2BiNormal)
 614 {
 615   Standard_Integer Index;
 616   Standard_Real Delta = 0.;
 617   if(IsSingular(Param, Index))
 618     if(SingularD2(Param, Index, Tangent, DTangent, D2Tangent,
 619                   Normal, DNormal, D2Normal,
 620                   BiNormal, DBiNormal, D2BiNormal,
 621                   Delta))
 622       return Standard_True;
 623
 624 //  Standard_Real Norma;
 625   Standard_Real theParam = Param + Delta;
 626   gp_Vec DC1, DC2, DC3, DC4;
 627   myTrimmed->D3(theParam, P, DC1, DC2, DC3);
 628   DC4 = myTrimmed->DN(theParam, 4);
 629
 630   Tangent = DC1.Normalized();
 631
 632   //if (DC2.Magnitude() <= NullTol || Tangent.Crossed(DC2).Magnitude() <= NullTol) {
 633   if (Tangent.Crossed(DC2).Magnitude() <= gp::Resolution()) {
 634     gp_Ax2 Axe (gp_Pnt(0,0,0), Tangent);
 635     Normal.SetXYZ(Axe.XDirection().XYZ());
 636     BiNormal.SetXYZ(Axe.YDirection().XYZ());
 637     DTangent.SetCoord(0,0,0);
 638     DNormal.SetCoord(0,0,0);
 639     DBiNormal.SetCoord(0,0,0);
 640     D2Tangent.SetCoord(0,0,0);
 641     D2Normal.SetCoord(0,0,0);
 642     D2BiNormal.SetCoord(0,0,0);
 643     return Standard_True;
 644   }
 645   else
 646     BiNormal = Tangent.Crossed(DC2).Normalized();
 647
 648   Normal = BiNormal.Crossed(Tangent);
 649
 650   DTangent = FDeriv(DC1, DC2);
 651   D2Tangent = DDeriv(DC1, DC2, DC3);
 652
 653   gp_Vec instead_DC1, instead_DC2, instead_DC3;
 654   instead_DC1 = Tangent.Crossed(DC2);
 655   instead_DC2 = DTangent.Crossed(DC2) + Tangent.Crossed(DC3);
 656   instead_DC3 = D2Tangent.Crossed(DC2) + 2*DTangent.Crossed(DC3) + Tangent.Crossed(DC4);
 657   DBiNormal = FDeriv(instead_DC1, instead_DC2);
 658   D2BiNormal = DDeriv(instead_DC1, instead_DC2, instead_DC3);
 659
 660   DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
 661
 662   D2Normal = D2BiNormal.Crossed(Tangent) + 2*DBiNormal.Crossed(DTangent) + BiNormal.Crossed(D2Tangent);
 663
 664   return Standard_True;
 665 }
 666
 667 //=======================================================================
 668 //function : NbIntervals
 669 //purpose  :
 670 //=======================================================================
 671
 672  Standard_Integer GeomFill_Frenet::NbIntervals(const GeomAbs_Shape S) const
 673 {
 674   GeomAbs_Shape tmpS = GeomAbs_C0;
 675   Standard_Integer NbTrimmed;
 676   switch (S) {
 677   case GeomAbs_C0: tmpS = GeomAbs_C2; break;
 678   case GeomAbs_C1: tmpS = GeomAbs_C3; break;
 679   case GeomAbs_C2:
 680   case GeomAbs_C3:
 681   case GeomAbs_CN: tmpS = GeomAbs_CN; break;
 682   default: throw Standard_OutOfRange();
 683   }
 684
 685   NbTrimmed = myCurve->NbIntervals(tmpS);
 686
 687   if (!isSngl) return NbTrimmed;
 688
 689   TColStd_Array1OfReal TrimInt(1, NbTrimmed + 1);
 690   myCurve->Intervals(TrimInt, tmpS);
 691
 692   TColStd_SequenceOfReal Fusion;
 693   GeomLib::FuseIntervals(TrimInt, mySngl->Array1(), Fusion);
 694
 695   return Fusion.Length() - 1;
 696 }
 697
 698 //=======================================================================
 699 //function : Intervals
 700 //purpose  :
 701 //=======================================================================
 702
 703  void GeomFill_Frenet::Intervals(TColStd_Array1OfReal& T,
 704                                  const GeomAbs_Shape S) const
 705 {
 706   GeomAbs_Shape tmpS = GeomAbs_C0;
 707   Standard_Integer NbTrimmed;
 708   switch (S) {
 709   case GeomAbs_C0: tmpS = GeomAbs_C2; break;
 710   case GeomAbs_C1: tmpS = GeomAbs_C3; break;
 711   case GeomAbs_C2:
 712   case GeomAbs_C3:
 713   case GeomAbs_CN: tmpS = GeomAbs_CN; break;
 714   default: throw Standard_OutOfRange();
 715   }
 716
 717   if (!isSngl) {
 718     myCurve->Intervals(T, tmpS);
 719     return;
 720   }
 721
 722   NbTrimmed = myCurve->NbIntervals(tmpS);
 723   TColStd_Array1OfReal TrimInt(1, NbTrimmed + 1);
 724   myCurve->Intervals(TrimInt, tmpS);
 725
 726   TColStd_SequenceOfReal Fusion;
 727   GeomLib::FuseIntervals(TrimInt, mySngl->Array1(), Fusion);
 728
 729   for (Standard_Integer i = 1; i <= Fusion.Length(); i++)
 730     T.ChangeValue(i) = Fusion.Value(i);
 731 }
 732
 733  void GeomFill_Frenet::GetAverageLaw(gp_Vec& ATangent,
 734                                      gp_Vec& ANormal,
 735                                      gp_Vec& ABiNormal)
 736 {
 737   Standard_Integer Num = 20; //order of digitalization
 738   gp_Vec T, N, BN;
 739   ATangent = gp_Vec(0, 0, 0);
 740   ANormal = gp_Vec(0, 0, 0);
 741   ABiNormal = gp_Vec(0, 0, 0);
 742   Standard_Real Step = (myTrimmed->LastParameter() -
 743                         myTrimmed->FirstParameter()) / Num;
 744   Standard_Real Param;
 745   for (Standard_Integer i = 0; i <= Num; i++) {
 746     Param = myTrimmed->FirstParameter() + i*Step;
 747     if (Param > myTrimmed->LastParameter()) Param = myTrimmed->LastParameter();
 748     D0(Param, T, N, BN);
 749     ATangent += T;
 750     ANormal += N;
 751     ABiNormal += BN;
 752   }
 753   ATangent /= Num + 1;
 754   ANormal /= Num + 1;
 755
 756   ATangent.Normalize();
 757   ABiNormal = ATangent.Crossed(ANormal).Normalized();
 758   ANormal = ABiNormal.Crossed(ATangent);
 759 }
 760
 761 //=======================================================================
 762 //function : IsConstant
 763 //purpose  :
 764 //=======================================================================
 765
 766  Standard_Boolean GeomFill_Frenet::IsConstant() const
 767 {
 768   return (myCurve->GetType() == GeomAbs_Line);
 769 }
 770
 771 //=======================================================================
 772 //function : IsOnlyBy3dCurve
 773 //purpose  :
 774 //=======================================================================
 775
 776  Standard_Boolean GeomFill_Frenet::IsOnlyBy3dCurve() const
 777 {
 778   return Standard_True;
 779 }
 780
 781 //=======================================================================
 782 //function : IsSingular
 783 //purpose  :
 784 //=======================================================================
 785
 786 Standard_Boolean GeomFill_Frenet::IsSingular(const Standard_Real U, Standard_Integer& Index) const
 787 {
 788   Standard_Integer i;
 789   if(!isSngl) return Standard_False;
 790   for(i = 1; i <= mySngl->Length(); i++) {
 791     if (Abs(U - mySngl->Value(i)) < mySnglLen->Value(i)) {
 792       Index = i;
 793       return Standard_True;
 794     }
 795   }
 796   return Standard_False;
 797 }
 798
 799 //=======================================================================
 800 //function : DoSingular
 801 //purpose  :
 802 //=======================================================================
 803
 804 Standard_Boolean GeomFill_Frenet::DoSingular(const Standard_Real U,
 805                                              const Standard_Integer Index,
 806                                              gp_Vec& Tangent,
 807                                              gp_Vec& BiNormal,
 808                                              Standard_Integer& n,
 809                                              Standard_Integer& k,
 810                                              Standard_Integer& TFlag,
 811                                              Standard_Integer& BNFlag,
 812                                              Standard_Real& Delta)
 813 {
 814   Standard_Integer i, MaxN = 20;
 815   Delta = 0.;
 816   Standard_Real h;
 817   h = 2*mySnglLen->Value(Index);
 818
 819   Standard_Real A, B;
 820   gp_Vec T, N, BN;
 821   TFlag = 1;
 822   BNFlag = 1;
 823   GetInterval(A, B);
 824   if (U >= (A + B)/2)
 825     h = -h;
 826   for(i = 1; i <= MaxN; i++) {
 827     Tangent = myTrimmed->DN(U, i);
 828     if(Tangent.Magnitude() > Precision::Confusion()) break;
 829   }
 830   if (i > MaxN) return Standard_False;
 831   Tangent.Normalize();
 832   n = i;
 833
 834   i++;
 835   for(; i <= MaxN; i++) {
 836     BiNormal = Tangent.Crossed(myTrimmed->DN(U, i));
 837     Standard_Real magn = BiNormal.Magnitude();
 838     if (magn > Precision::Confusion())
 839       {
 840         //modified by jgv, 12.08.03 for OCC605 ////
 841         gp_Vec NextBiNormal = Tangent.Crossed(myTrimmed->DN(U, i+1));
 842         if (NextBiNormal.Magnitude() > magn)
 843           {
 844             i++;
 845             BiNormal = NextBiNormal;
 846           }
 847         ///////////////////////////////////////////
 848         break;
 849       }
 850   }
 851   if (i > MaxN)
 852   {
 853     Delta = h;
 854     return Standard_False;
 855   }
 856
 857   BiNormal.Normalize();
 858   k = i;
 859
 860   D0(U + h, T, N, BN);
 861
 862   if(Tangent.Angle(T) > M_PI/2) TFlag = -1;
 863   if(BiNormal.Angle(BN) > M_PI/2) BNFlag = -1;
 864
 865   return Standard_True;
 866 }
 867
 868  Standard_Boolean GeomFill_Frenet::SingularD0(const Standard_Real Param,
 869                                               const Standard_Integer Index,
 870                                               gp_Vec& Tangent,
 871                                               gp_Vec& Normal,
 872                                               gp_Vec& BiNormal,
 873                                               Standard_Real& Delta)
 874 {
 875   Standard_Integer n, k, TFlag, BNFlag;
 876   if(!DoSingular(Param, Index, Tangent, BiNormal,
 877                  n, k, TFlag, BNFlag, Delta))
 878     return Standard_False;
 879
 880   Tangent *= TFlag;
 881   BiNormal *= BNFlag;
 882   Normal = BiNormal;
 883   Normal.Cross(Tangent);
 884
 885   return Standard_True;
 886 }
 887
 888  Standard_Boolean GeomFill_Frenet::SingularD1(const Standard_Real Param,
 889                                               const Standard_Integer Index,
 890                                               gp_Vec& Tangent,gp_Vec& DTangent,
 891                                               gp_Vec& Normal,gp_Vec& DNormal,
 892                                               gp_Vec& BiNormal,gp_Vec& DBiNormal,
 893                                               Standard_Real& Delta)
 894 {
 895   Standard_Integer n, k, TFlag, BNFlag;
 896   if(!DoSingular(Param, Index, Tangent, BiNormal, n, k, TFlag, BNFlag, Delta))
 897     return Standard_False;
 898
 899   gp_Vec F, DF, Dtmp;
 900   F = myTrimmed->DN(Param, n);
 901   DF = myTrimmed->DN(Param, n+1);
 902   DTangent = FDeriv(F, DF);
 903
 904   Dtmp = myTrimmed->DN(Param, k);
 905   F = Tangent.Crossed(Dtmp);
 906   DF = DTangent.Crossed(Dtmp) + Tangent.Crossed(myTrimmed->DN(Param, k+1));
 907   DBiNormal = FDeriv(F, DF);
 908
 909   if(TFlag < 0) {
 910     Tangent = -Tangent;
 911     DTangent = -DTangent;
 912   }
 913
 914   if(BNFlag < 0) {
 915     BiNormal = -BiNormal;
 916     DBiNormal = -DBiNormal;
 917   }
 918
 919   Normal = BiNormal.Crossed(Tangent);
 920   DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
 921
 922   return Standard_True;
 923 }
 924
 925  Standard_Boolean GeomFill_Frenet::SingularD2(const Standard_Real Param,
 926                                               const Standard_Integer Index,
 927                                               gp_Vec& Tangent,
 928                                               gp_Vec& DTangent,
 929                                               gp_Vec& D2Tangent,
 930                                               gp_Vec& Normal,
 931                                               gp_Vec& DNormal,
 932                                               gp_Vec& D2Normal,
 933                                               gp_Vec& BiNormal,
 934                                               gp_Vec& DBiNormal,
 935                                               gp_Vec& D2BiNormal,
 936                                               Standard_Real& Delta)
 937 {
 938   Standard_Integer n, k, TFlag, BNFlag;
 939   if(!DoSingular(Param, Index, Tangent, BiNormal, n, k, TFlag, BNFlag, Delta))
 940     return Standard_False;
 941
 942   gp_Vec F, DF, D2F, Dtmp1, Dtmp2;
 943   F = myTrimmed->DN(Param, n);
 944   DF = myTrimmed->DN(Param, n+1);
 945   D2F = myTrimmed->DN(Param, n+2);
 946   DTangent = FDeriv(F, DF);
 947   D2Tangent = DDeriv(F, DF, D2F);
 948
 949   Dtmp1 = myTrimmed->DN(Param, k);
 950   Dtmp2 = myTrimmed->DN(Param, k+1);
 951   F = Tangent.Crossed(Dtmp1);
 952   DF = DTangent.Crossed(Dtmp1) + Tangent.Crossed(Dtmp2);
 953   D2F = D2Tangent.Crossed(Dtmp1) + 2*DTangent.Crossed(Dtmp2) +
 954     Tangent.Crossed(myTrimmed->DN(Param, k+2));
 955   DBiNormal = FDeriv(F, DF);
 956   D2BiNormal = DDeriv(F, DF, D2F);
 957
 958   if(TFlag < 0) {
 959     Tangent = -Tangent;
 960     DTangent = -DTangent;
 961     D2Tangent = -D2Tangent;
 962   }
 963
 964   if(BNFlag < 0) {
 965     BiNormal = -BiNormal;
 966     DBiNormal = -DBiNormal;
 967     D2BiNormal = -D2BiNormal;
 968   }
 969
 970   Normal = BiNormal.Crossed(Tangent);
 971   DNormal = DBiNormal.Crossed(Tangent) + BiNormal.Crossed(DTangent);
 972   D2Normal = D2BiNormal.Crossed(Tangent) + 2*DBiNormal.Crossed(DTangent) + BiNormal.Crossed(D2Tangent);
 973
 974   return Standard_True;
 975 }