0031035: Coding - uninitialized class fields reported by Visual Studio Code Analysis
[occt.git] / 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 : isSngl(Standard_False)
99 {
100 }
101
102
103 //=======================================================================
104 //function : Copy
105 //purpose  : 
106 //=======================================================================
107
108 Handle(GeomFill_TrihedronLaw) GeomFill_Frenet::Copy() const
109 {
110   Handle(GeomFill_Frenet) copy = new (GeomFill_Frenet)();
111   if (!myCurve.IsNull()) copy->SetCurve(myCurve);
112   return copy;
113 }
114
115 //=======================================================================
116 //function : SetCurve
117 //purpose  : 
118 //=======================================================================
119
120  void GeomFill_Frenet::SetCurve(const Handle(Adaptor3d_HCurve)& C) 
121 {
122   GeomFill_TrihedronLaw::SetCurve(C);
123   if (! C.IsNull()) {
124     GeomAbs_CurveType type;
125     type = C->GetType();
126     switch  (type) {
127     case GeomAbs_Circle:
128     case GeomAbs_Ellipse:
129     case GeomAbs_Hyperbola:
130     case GeomAbs_Parabola:
131     case GeomAbs_Line:
132       {
133         // No probleme
134         isSngl = Standard_False;
135         break;
136       }
137      default :
138        { 
139          // We have to search singulaties
140          Init();
141        }
142     }
143   }
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);
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);
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 }