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