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