0031035: Coding - uninitialized class fields reported by Visual Studio Code Analysis
[occt.git] / src / GeomFill / GeomFill_SectionPlacement.cxx
1 // Created on: 1997-12-15
2 // Created by: Philippe MANGIN
3 // Copyright (c) 1997-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 <Bnd_Box.hxx>
20 #include <BndLib_Add3dCurve.hxx>
21 #include <Extrema_ExtCC.hxx>
22 #include <Extrema_POnCurv.hxx>
23 #include <Geom_BezierCurve.hxx>
24 #include <Geom_BSplineCurve.hxx>
25 #include <Geom_CartesianPoint.hxx>
26 #include <Geom_Conic.hxx>
27 #include <Geom_Curve.hxx>
28 #include <Geom_Geometry.hxx>
29 #include <Geom_Line.hxx>
30 #include <Geom_Plane.hxx>
31 #include <Geom_TrimmedCurve.hxx>
32 #include <GeomAbs_CurveType.hxx>
33 #include <GeomAdaptor_HSurface.hxx>
34 #include <GeomFill_LocationLaw.hxx>
35 #include <GeomFill_SectionPlacement.hxx>
36 #include <GeomLib.hxx>
37 #include <GeomLProp_CLProps.hxx>
38 #include <gp.hxx>
39 #include <gp_Ax2.hxx>
40 #include <gp_Ax3.hxx>
41 #include <gp_Dir.hxx>
42 #include <gp_Mat.hxx>
43 #include <gp_Pnt.hxx>
44 #include <gp_Trsf.hxx>
45 #include <gp_Vec.hxx>
46 #include <IntCurveSurface_HInter.hxx>
47 #include <IntCurveSurface_IntersectionPoint.hxx>
48 #include <Precision.hxx>
49 #include <StdFail_NotDone.hxx>
50 #include <TColgp_HArray1OfPnt.hxx>
51
52 //===============================================================
53 // Function :
54 // Purpose :
55 //===============================================================
56 static void Tangente(const Adaptor3d_Curve& Path,
57                      const Standard_Real Param,
58                      gp_Pnt& P,
59                      gp_Vec& Tang)
60 {
61   Path.D1(Param, P, Tang);
62   Standard_Real Norm = Tang.Magnitude();
63
64   for (Standard_Integer ii=2; (ii<12) && (Norm < Precision::Confusion()); 
65        ii++) {
66     Tang =  Path.DN(Param, ii);
67     Norm = Tang.Magnitude();
68   }
69
70   if (Norm > 100*gp::Resolution()) Tang /= Norm; 
71 }
72                    
73
74 static Standard_Real Penalite(const Standard_Real angle,
75                               const Standard_Real dist)
76 {
77   Standard_Real penal;
78   
79   if (dist < 1)
80     penal = Sqrt(dist);
81   else if (dist <2)
82     penal = Pow(dist, 2);
83   else
84     penal = dist + 2;
85
86   if (angle > 1.e-3) {
87     penal += 1./angle -2./M_PI;
88   }
89   else {
90     penal += 1.e3;
91   }
92
93   return penal;
94 }
95
96 static Standard_Real EvalAngle(const gp_Vec& V1, 
97                                const gp_Vec& V2)
98 {
99  Standard_Real angle;
100  angle = V1.Angle(V2);
101  if (angle > M_PI/2) angle = M_PI - angle;
102  return angle;
103 }
104
105 static Standard_Integer NbSamples(const Handle(Geom_Curve)& aCurve)
106 {
107   Standard_Real nbs = 100.; //on default
108
109   Handle(Geom_Curve) theCurve = aCurve;
110   if (aCurve->IsInstance(STANDARD_TYPE(Geom_TrimmedCurve)))
111     theCurve = (Handle(Geom_TrimmedCurve)::DownCast(aCurve))->BasisCurve();
112
113   if (theCurve->IsInstance(STANDARD_TYPE(Geom_Line)))
114     nbs = 1;
115   else if (theCurve->IsKind(STANDARD_TYPE(Geom_Conic)))
116     nbs = 4;
117   else if (theCurve->IsInstance(STANDARD_TYPE(Geom_BezierCurve)))
118     {
119       Handle(Geom_BezierCurve) BC = Handle(Geom_BezierCurve)::DownCast(theCurve);
120       nbs = 3 + BC->NbPoles();
121     }
122   else if (theCurve->IsInstance(STANDARD_TYPE(Geom_BSplineCurve)))
123     {
124       Handle(Geom_BSplineCurve) BC = Handle(Geom_BSplineCurve)::DownCast(theCurve);
125       nbs  = BC->NbKnots();
126       nbs *= BC->Degree();
127       Standard_Real ratio =
128         (aCurve->LastParameter() - aCurve->FirstParameter())/(BC->LastParameter() - BC->FirstParameter());
129       nbs *= ratio;
130       if(nbs < 4.0)
131         nbs = 4;
132     }
133
134   if (nbs > 300.)
135     nbs = 300;
136
137   return ((Standard_Integer)nbs);
138 }
139
140 //===============================================================
141 // Function :DistMini
142 // Purpose : Examine un extrema pour updater <Dist> & <Param>
143 //===============================================================
144 static void DistMini(const Extrema_ExtPC& Ext,
145                      const Adaptor3d_Curve& C,
146                      Standard_Real& Dist,
147                      Standard_Real& Param)
148 {
149   Standard_Real dist1, dist2;
150   Standard_Integer ii;
151   gp_Pnt P1, P2;
152   Standard_Real Dist2 = RealLast();
153
154   Ext.TrimmedSquareDistances(dist1, dist2, P1, P2);
155   if ( (dist1<Dist2) || (dist2<Dist2) ) {
156     if (dist1 < dist2) {
157       Dist2 = dist1;
158       Param  = C.FirstParameter();
159     }
160     else {
161       Dist2 = dist2;
162       Param  = C.LastParameter();
163     }
164   }
165
166   if (Ext.IsDone())
167     {
168       for (ii=1; ii<= Ext.NbExt(); ii++) {
169         if (Ext.SquareDistance(ii) < Dist2) {
170           Dist2 = Ext.SquareDistance(ii);
171           Param  =  Ext.Point(ii).Parameter();
172         }
173       }
174     }
175   Dist = sqrt (Dist2);
176 }
177                      
178
179 //===============================================================
180 // Function : Constructeur
181 // Purpose :
182 //===============================================================
183 GeomFill_SectionPlacement::
184 GeomFill_SectionPlacement(const Handle(GeomFill_LocationLaw)& L,
185               const Handle(Geom_Geometry)& Section) :
186               myLaw(L), /* myAdpSection(Section),  mySection(Section), */
187               SecParam(0.0), PathParam(0.0),
188               Dist(RealLast()), AngleMax(0.)
189 {
190
191   done   = Standard_False;
192   isplan = Standard_False;
193   myIsPoint = Standard_False;
194
195   if (Section->IsInstance(STANDARD_TYPE(Geom_CartesianPoint)))
196     {
197       myIsPoint = Standard_True;
198       myPoint   = Handle(Geom_CartesianPoint)::DownCast(Section)->Pnt();
199       isplan    = Standard_True;
200     }
201   else
202     {
203       Handle(Geom_Curve) CurveSection = Handle(Geom_Curve)::DownCast(Section);
204       myAdpSection.Load(CurveSection);
205       mySection = CurveSection;
206     }
207
208   Standard_Integer i, j, NbPoles=0;
209   Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
210
211   // Boite d'encombrement de la section pour en deduire le gabarit
212   Bnd_Box box;
213   if (myIsPoint)
214     box.Add(myPoint);
215   else
216     BndLib_Add3dCurve::Add(myAdpSection, 1.e-4, box);
217   box.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
218
219   Standard_Real DX = aXmax - aXmin ;
220   Standard_Real DY = aYmax - aYmin ;
221   Standard_Real DZ = aZmax - aZmin ;
222   Gabarit = Sqrt( DX*DX+ DY*DY + DZ*DZ )/2. ;
223
224   Gabarit += Precision::Confusion(); // Cas des toute petite
225   
226
227   // Initialisation de TheAxe pour les cas singulier
228   if (!myIsPoint)
229     {
230       gp_Pnt P;
231       gp_Vec V;
232       Tangente(myAdpSection, 
233                (myAdpSection.FirstParameter()+myAdpSection.LastParameter())/2,
234                P, V);
235       TheAxe.SetLocation(P);
236       TheAxe.SetDirection(V);
237  
238       // y a t'il un Plan moyen ?
239       GeomAbs_CurveType TheType = myAdpSection.GetType();
240       switch  (TheType) {
241       case GeomAbs_Circle:
242         {
243           isplan = Standard_True;
244           TheAxe =  myAdpSection.Circle().Axis();
245           break;
246         }
247       case GeomAbs_Ellipse:
248         {
249           isplan = Standard_True;
250           TheAxe =  myAdpSection.Ellipse().Axis();
251           break;
252         }
253       case GeomAbs_Hyperbola:
254         {
255           isplan = Standard_True;
256           TheAxe =  myAdpSection.Hyperbola().Axis();
257           break;
258         }
259       case GeomAbs_Parabola:
260         {
261           isplan = Standard_True;
262           TheAxe =  myAdpSection.Parabola().Axis();
263           break;
264         }
265       case GeomAbs_Line:
266         {
267           NbPoles = 0; // Pas de Plan !!
268           break;
269         }
270       case GeomAbs_BezierCurve:
271       case GeomAbs_BSplineCurve:
272         {
273           NbPoles = myAdpSection.NbPoles();
274           break;
275         }
276       default:
277         NbPoles = 21;
278       }
279       
280       
281       if (!isplan && NbPoles>2)
282         {
283           // Calcul d'un plan moyen.
284           Handle(TColgp_HArray1OfPnt) Pnts;
285           Standard_Real first = myAdpSection.FirstParameter();
286           Standard_Real last =  myAdpSection.LastParameter();
287           if (myAdpSection.IsPeriodic())
288           {
289             //Correct boundaries to avoid mistake of LocateU
290             Handle(Geom_Curve) aCurve = myAdpSection.Curve();
291             if (aCurve->IsInstance(STANDARD_TYPE(Geom_TrimmedCurve)))
292               aCurve = (Handle(Geom_TrimmedCurve)::DownCast(aCurve))->BasisCurve();
293             Standard_Real Ufirst = aCurve->FirstParameter();
294             Standard_Real aPeriod = aCurve->Period();
295             Standard_Real U1 = Ufirst + Floor((first - Ufirst)/aPeriod) * aPeriod;
296             Standard_Real U2 = U1 + aPeriod;
297             if (Abs(first - U1) <= Precision::PConfusion())
298               first = U1;
299             if (Abs(last - U2) <= Precision::PConfusion())
300               last = U2;
301           }
302           Standard_Real t, delta;
303           if (myAdpSection.GetType() == GeomAbs_BSplineCurve)
304             {
305               Handle(Geom_BSplineCurve) BC =
306                 Handle(Geom_BSplineCurve)::DownCast(myAdpSection.Curve());
307               Standard_Integer I1, I2, I3, I4;
308               BC->LocateU( first, Precision::Confusion(), I1, I2 );
309               BC->LocateU( last,  Precision::Confusion(), I3, I4 );
310               Standard_Integer NbKnots = I3 - I2 + 1;
311               
312               Standard_Integer NbLocalPnts = 10;
313               Standard_Integer NbPnts = (NbKnots-1) * NbLocalPnts;
314               if (I1 != I2)
315                 NbPnts += NbLocalPnts;
316               if (I3 != I4 && first < BC->Knot(I3))
317                 NbPnts += NbLocalPnts;
318               if (!myAdpSection.IsClosed())
319                 NbPnts++;
320               Pnts = new TColgp_HArray1OfPnt(1, NbPnts);
321               Standard_Integer nb = 1;
322               if (I1 != I2)
323                 {
324                   Standard_Real locallast = (BC->Knot(I2) < last)? BC->Knot(I2) : last;
325                   delta = (locallast - first) / NbLocalPnts;
326                   for (j = 0; j < NbLocalPnts; j++)
327                     {
328                       t = first + j*delta;
329                       Pnts->SetValue( nb++, myAdpSection.Value(t) );
330                     }
331                 }
332               for (i = I2; i < I3; i++)
333                 {
334                   t = BC->Knot(i);
335                   delta = (BC->Knot(i+1) - t) / NbLocalPnts;
336                   for (j = 0; j < NbLocalPnts; j++)
337                     {
338                       Pnts->SetValue( nb++, myAdpSection.Value(t) );
339                       t += delta;
340                     }
341                 }
342               if (I3 != I4 && first < BC->Knot(I3))
343                 {
344                   t = BC->Knot(I3);
345                   delta = (last - t) / NbLocalPnts;
346                   for (j = 0; j < NbLocalPnts; j++)
347                     {
348                       Pnts->SetValue( nb++, myAdpSection.Value(t) );
349                       t += delta;
350                     }
351                 }
352               if (!myAdpSection.IsClosed())
353                 Pnts->SetValue( nb, myAdpSection.Value(last) );
354             }
355           else // other type
356             {
357               Standard_Integer NbPnts = NbPoles-1;
358               if (!myAdpSection.IsClosed())
359                 NbPnts++;
360               Pnts = new TColgp_HArray1OfPnt(1, NbPnts);
361               delta = (last - first) / (NbPoles-1);
362               for (i = 0; i < NbPoles-1; i++)
363                 {
364                   t = first + i*delta;
365                   Pnts->SetValue( i+1, myAdpSection.Value(t) );
366                 }
367               if (!myAdpSection.IsClosed())
368                 Pnts->SetValue( NbPnts, myAdpSection.Value(last) );
369             }
370           
371           Standard_Boolean issing;
372           gp_Ax2 axe;
373           GeomLib::AxeOfInertia(Pnts->Array1(), axe, issing, Precision::Confusion());
374           if (!issing) {
375             isplan = Standard_True;
376             TheAxe.SetLocation(axe.Location());
377             TheAxe.SetDirection(axe.Direction());
378           }
379         }
380       
381       
382       myExt.Initialize(myAdpSection, 
383                        myAdpSection.FirstParameter(), 
384                        myAdpSection.LastParameter(), 
385                        Precision::Confusion());
386     }
387 }
388
389 //===============================================================
390 // Function :SetLocation
391 // Purpose :
392 //===============================================================
393 void GeomFill_SectionPlacement::
394 SetLocation(const Handle(GeomFill_LocationLaw)& L) 
395 {
396   myLaw = L;
397 }
398
399 //===============================================================
400 // Function : Perform
401 // Purpose : Le plus simple
402 //===============================================================
403 void GeomFill_SectionPlacement::Perform(const Standard_Real Tol) 
404 {
405   Handle(Adaptor3d_HCurve) Path;
406   Path = myLaw->GetCurve();
407   Perform(Path, Tol);
408 }
409
410 //===============================================================
411 // Function :Perform
412 // Purpose : Recherche automatique
413 //===============================================================
414 void GeomFill_SectionPlacement::Perform(const Handle(Adaptor3d_HCurve)& Path,
415                                         const Standard_Real Tol) 
416 {
417   Standard_Real IntTol = 1.e-5;
418   Standard_Real DistCenter = Precision::Infinite();
419
420   if (myIsPoint)
421     {
422       Extrema_ExtPC Projector(myPoint, Path->Curve(), Precision::Confusion());
423       DistMini( Projector, Path->Curve(), Dist, PathParam );
424       AngleMax = M_PI/2;
425     }
426   else
427     {
428       PathParam = Path->FirstParameter();
429       SecParam =  myAdpSection.FirstParameter();
430       
431       Standard_Real distaux, taux = 0.0, alpha;
432       gp_Pnt PonPath, PonSec, P;
433       gp_Vec VRef, dp1;
434       VRef.SetXYZ(TheAxe.Direction().XYZ());
435       
436       Tangente( Path->Curve(), PathParam, PonPath, dp1);
437       PonSec = myAdpSection.Value(SecParam);
438       Dist =  PonPath.Distance(PonSec);
439       if (Dist > Tol) { // On Cherche un meilleur point sur la section
440         myExt.Perform(PonPath);  
441         if ( myExt.IsDone() ) { 
442           DistMini(myExt, myAdpSection, Dist, SecParam);
443           PonSec = myAdpSection.Value(SecParam);
444         }
445       }
446       AngleMax = EvalAngle(VRef, dp1);
447       if (isplan) AngleMax = M_PI/2 - AngleMax;
448       
449       Standard_Boolean Trouve = Standard_False; 
450       Standard_Integer ii;
451       
452       if (isplan) {
453         // (1.1) Distances Point-Plan
454         Standard_Real DistPlan;
455         gp_Vec V1 (PonPath, TheAxe.Location());
456         DistPlan = Abs(V1.Dot(VRef));
457         if (DistPlan <= IntTol)
458           DistCenter = V1.Magnitude();
459         
460         gp_Pnt Plast = Path->Value( Path->LastParameter() );
461         V1.SetXYZ( TheAxe.Location().XYZ() - Plast.XYZ() );
462         DistPlan = Abs(V1.Dot(VRef));
463         if (DistPlan <= IntTol)
464           {
465             Standard_Real aDist = V1.Magnitude();
466             if (aDist < DistCenter)
467               {
468                 DistCenter = aDist;
469                 PonPath = Plast;
470                 PathParam = Path->LastParameter();
471               }
472           }
473         
474         // (1.2) Intersection Plan-courbe
475         gp_Ax3 axe (TheAxe.Location(), TheAxe.Direction());
476         Handle(Geom_Plane) plan = new (Geom_Plane)(axe);
477         Handle(GeomAdaptor_HSurface) adplan = 
478           new (GeomAdaptor_HSurface)(plan);
479         IntCurveSurface_HInter Intersector;
480         Intersector.Perform(Path, adplan);
481         if (Intersector.IsDone())
482           {
483             Standard_Real w;
484             Standard_Real aDist;
485             for (ii=1; ii<=Intersector.NbPoints(); ii++)
486               {
487                 w = Intersector.Point(ii).W();
488                 P = Path->Value( w );
489                 aDist = P.Distance( TheAxe.Location() );
490                 if (aDist < DistCenter)
491                   {
492                     DistCenter = aDist;
493                     PonPath = P;
494                     PathParam = w;
495                   }
496               }
497           }
498         if (!Intersector.IsDone() || Intersector.NbPoints() == 0)
499           {
500             Standard_Integer NbPnts = NbSamples( mySection );
501             TColgp_Array1OfPnt Pnts( 1, NbPnts+1 );
502             Standard_Real delta = (mySection->LastParameter()-mySection->FirstParameter())/NbPnts;
503             for (ii = 0; ii <= NbPnts; ii++)
504               Pnts(ii+1) = mySection->Value( mySection->FirstParameter() + ii*delta );
505             
506             gp_Pnt BaryCenter;
507             gp_Dir Xdir, Ydir;
508             Standard_Real Xgap, Ygap, Zgap;
509             GeomLib::Inertia( Pnts, BaryCenter, Xdir, Ydir, Xgap, Ygap, Zgap );
510             
511             gp_Pnt Pfirst = Path->Value( Path->FirstParameter() );
512             if (Pfirst.Distance(BaryCenter) < Plast.Distance(BaryCenter))
513               PathParam = Path->FirstParameter();
514             else
515               {
516                 PathParam = Path->LastParameter();
517                 Tangente( Path->Curve(), PathParam, PonPath, dp1);
518                 PonSec = myAdpSection.Value(SecParam);
519                 Dist =  PonPath.Distance(PonSec);
520                 if (Dist > Tol) { // On Cherche un meilleur point sur la section
521                   myExt.Perform(PonPath);  
522                   if ( myExt.IsDone() ) { 
523                     DistMini(myExt, myAdpSection, Dist, SecParam);
524                     PonSec = myAdpSection.Value(SecParam);
525                   }
526                 }
527                 AngleMax = EvalAngle(VRef, dp1);
528                 AngleMax = M_PI/2 - AngleMax;
529               }
530           }
531     
532         /*
533            // (1.1) Distances Point-Plan
534            Standard_Real DistPlan;
535            gp_Vec V1 (PonPath, TheAxe.Location());
536            DistPlan = Abs(V1.Dot(VRef));
537            
538            // On examine l'autre extremite
539            gp_Pnt P;
540            Tangente(Path->Curve(), Path->LastParameter(), P, dp1);
541            V1.SetXYZ(TheAxe.Location().XYZ()-P.XYZ());
542            if  (Abs(V1.Dot(VRef)) <= DistPlan ) { // On prend l'autre extremite
543            alpha =  M_PI/2 - EvalAngle(VRef, dp1);
544            distaux =  PonPath.Distance(PonSec);
545            if (distaux > Tol) {
546            myExt.Perform(P);  
547            if ( myExt.IsDone() ) 
548            DistMini(myExt, myAdpSection, distaux, taux);
549            }
550            else 
551            taux = SecParam;
552            
553            if (Choix(distaux, alpha)) {
554            Dist = distaux;
555            AngleMax = alpha;
556            PonPath = P;
557            PathParam = Path->LastParameter();
558            }
559            }
560            
561            // (1.2) Intersection Plan-courbe
562            gp_Ax3 axe (TheAxe.Location(), TheAxe.Direction());
563            Handle(Geom_Plane) plan = new (Geom_Plane)(axe);
564            Handle(GeomAdaptor_HSurface) adplan = 
565            new (GeomAdaptor_HSurface)(plan);
566            IntCurveSurface_HInter Intersector;
567            Intersector.Perform(Path, adplan);
568            if (Intersector.IsDone()) {
569            Standard_Real w;
570            gp_Vec V;
571            for (ii=1; ii<=Intersector.NbPoints(); ii++) {
572            w = Intersector.Point(ii).W();
573            //(1.3) test d'angle
574            Tangente( Path->Curve(), w, P, V);
575            alpha = M_PI/2 - EvalAngle(V, VRef);
576            //(1.4) Test de distance Point-Courbe
577            myExt.Perform(P);  
578            if ( myExt.IsDone() ) {
579            DistMini(myExt, myAdpSection, distaux, taux);
580            if (Choix(distaux, alpha)) {
581            Dist = distaux;
582            SecParam  = taux;
583            AngleMax = alpha;
584            PonPath = P;
585            PathParam = w;
586            PonSec = myAdpSection.Value(SecParam);
587            }
588            }
589            else {
590            distaux = P.Distance(PonSec);
591            if (Choix(distaux, alpha)) {
592            Dist = distaux;
593            AngleMax = alpha;
594            PonPath = P;
595            PathParam = w;
596            }
597            }
598            }
599            }
600     */
601 #ifdef OCCT_DEBUG
602         if (Intersector.NbPoints() == 0) {
603           Intersector.Dump();
604         } 
605 #endif
606       } 
607       
608       // Cas General
609       if (! isplan) {
610         // (2.1) Distance avec les extremites ...
611         myExt.Perform(PonPath);  
612         if ( myExt.IsDone() ) {
613           DistMini(myExt, myAdpSection, distaux, taux);
614           if (distaux < Dist) {
615             Dist = distaux;
616             SecParam  = taux;
617           } 
618         }
619         Trouve = (Dist <= Tol);
620         if (!Trouve) {
621           Tangente( Path->Curve(), Path->LastParameter(), P, dp1);
622           alpha = EvalAngle(VRef, dp1);
623           myExt.Perform(P);     
624           if ( myExt.IsDone() ) {
625             if ( myExt.IsDone() ) {
626               DistMini(myExt, myAdpSection, distaux, taux);
627               if (Choix(distaux, alpha)) {
628                 Dist = distaux;
629                 SecParam  = taux;
630                 AngleMax = alpha;
631                 PonPath = P;
632                 PathParam = Path->LastParameter();
633               }
634             }
635           }
636           Trouve = (Dist <= Tol); 
637         }
638         
639         // (2.2) Distance courbe-courbe
640         if (!Trouve) {
641           Extrema_ExtCC Ext(Path->Curve(), myAdpSection, 
642                             Path->FirstParameter(), Path->LastParameter(),
643                             myAdpSection.FirstParameter(), 
644                             myAdpSection.LastParameter(),
645                             Path->Resolution(Tol/100), 
646                             myAdpSection.Resolution(Tol/100));
647           if (Ext.IsDone()) {
648             Extrema_POnCurv P1, P2;
649             for (ii=1; ii<=Ext.NbExt(); ii++) {
650               distaux = sqrt (Ext.SquareDistance(ii));
651               Ext.Points(ii, P1, P2);
652               Tangente(Path->Curve(), P1.Parameter(), P, dp1);
653               alpha =  EvalAngle(VRef, dp1);
654               if (Choix(distaux, alpha)) {
655                 Trouve = Standard_True;
656                 Dist = distaux;
657                 PathParam = P1.Parameter();
658                 SecParam = P2.Parameter();
659                 PonSec =  P2.Value();
660                 PonPath = P;
661                 AngleMax = alpha;
662               }
663             }
664           }
665           if (!Trouve){ 
666             // Si l'on a toujours rien, on essai une distance point/path
667             // c'est la derniere chance.
668             Extrema_ExtPC PExt;
669             PExt.Initialize(Path->Curve(), 
670                             Path->FirstParameter(),  
671                             Path->LastParameter(),
672                             Precision::Confusion());
673             PExt.Perform(PonSec);
674             if ( PExt.IsDone() ) {
675               // modified for OCC13595  Tue Oct 17 15:00:08 2006.BEGIN
676               //              DistMini(PExt, myAdpSection, distaux, taux);
677               DistMini(PExt, Path->Curve(), distaux, taux);
678               // modified for OCC13595  Tue Oct 17 15:00:11 2006.END
679               Tangente(Path->Curve(), taux, P, dp1);
680               alpha = EvalAngle(VRef, dp1);
681               if (Choix(distaux, alpha)) {
682                 Dist = distaux;
683                 PonPath = P;
684                 AngleMax = alpha;
685                 PathParam = taux;
686               }
687             }
688           }
689         }
690       }
691     }
692
693   done = Standard_True;
694 }
695
696
697 //===============================================================
698 // Function : Perform
699 // Purpose : Calcul le placement pour un parametre donne.
700 //===============================================================
701 void GeomFill_SectionPlacement::Perform(const Standard_Real Param,
702                                         const Standard_Real Tol) 
703 {
704   done = Standard_True;
705   Handle(Adaptor3d_HCurve) Path;
706   Path = myLaw->GetCurve();
707
708   PathParam = Param;
709   if (myIsPoint)
710     {
711       gp_Pnt PonPath = Path->Value( PathParam );
712       Dist = PonPath.Distance(myPoint);
713       AngleMax = M_PI/2;
714     }
715   else
716     {
717       SecParam =  myAdpSection.FirstParameter();
718       
719       //  Standard_Real distaux, taux, alpha;
720       //  gp_Pnt PonPath, PonSec, P;
721       gp_Pnt PonPath, PonSec;
722       gp_Vec VRef, dp1;
723       VRef.SetXYZ(TheAxe.Direction().XYZ()); 
724       
725       Tangente( Path->Curve(), PathParam, PonPath, dp1);
726       PonSec = myAdpSection.Value(SecParam);
727       Dist =  PonPath.Distance(PonSec);
728       if (Dist > Tol) { // On Cherche un meilleur point sur la section
729         myExt.Perform(PonPath);  
730         if ( myExt.IsDone() ) { 
731           DistMini(myExt, myAdpSection, Dist, SecParam);
732           PonSec = myAdpSection.Value(SecParam);
733         }
734       }
735       AngleMax = EvalAngle(VRef, dp1);
736       if (isplan) AngleMax = M_PI/2 - AngleMax;
737     }
738
739   done = Standard_True;
740 }
741
742 //===============================================================
743 // Function :IsDone
744 // Purpose :
745 //===============================================================
746  Standard_Boolean GeomFill_SectionPlacement::IsDone() const
747 {
748   return done; 
749 }
750
751 //===============================================================
752 // Function : ParameterOnPath
753 // Purpose :
754 //===============================================================
755  Standard_Real GeomFill_SectionPlacement::ParameterOnPath() const
756 {
757   return PathParam;
758 }
759
760 //===============================================================
761 // Function : ParameterOnSection
762 // Purpose :
763 //===============================================================
764  Standard_Real GeomFill_SectionPlacement::ParameterOnSection() const
765 {
766  return SecParam;
767 }
768
769 //===============================================================
770 // Function : Distance
771 // Purpose :
772 //===============================================================
773  Standard_Real GeomFill_SectionPlacement::Distance() const
774 {
775   return Dist;
776 }
777
778 //===============================================================
779 // Function : Angle
780 // Purpose :
781 //===============================================================
782  Standard_Real GeomFill_SectionPlacement::Angle() const
783 {
784   return AngleMax;
785 }
786
787 //===============================================================
788 // Function : Transformation
789 // Purpose :
790 //===============================================================
791  gp_Trsf GeomFill_SectionPlacement::Transformation(
792          const Standard_Boolean WithTranslation,
793          const Standard_Boolean WithCorrection) const
794 {
795   gp_Vec V;
796   gp_Mat M;
797   gp_Dir DT, DN, D;
798 // modified by NIZHNY-MKK  Fri Oct 17 15:27:07 2003
799   gp_Pnt P(0., 0., 0.), PSection(0., 0., 0.);
800
801   // Calcul des reperes
802   myLaw->D0(PathParam, M, V);
803
804   P.SetXYZ(V.XYZ());
805   D.SetXYZ (M.Column(3));
806   DN.SetXYZ(M.Column(1));
807   gp_Ax3 Paxe(P, D, DN);
808
809
810   if (WithTranslation || WithCorrection) {
811     if (myIsPoint)
812       PSection = myPoint;
813     else
814       PSection = mySection->Value(SecParam);   
815   }
816   // modified by NIZHNY-MKK  Wed Oct  8 15:03:19 2003.BEGIN
817   gp_Trsf Rot;
818
819   if (WithCorrection && !myIsPoint) { 
820     Standard_Real angle;
821     
822     if (!isplan)
823       throw Standard_Failure("Illegal usage: can't rotate non-planar profile");
824     
825     gp_Dir ProfileNormal = TheAxe.Direction();
826     gp_Dir SpineStartDir = Paxe.Direction();
827
828     if (! ProfileNormal.IsParallel( SpineStartDir, Precision::Angular() )) {
829       gp_Dir DirAxeOfRotation = ProfileNormal ^ SpineStartDir;
830       angle = ProfileNormal.AngleWithRef( SpineStartDir, DirAxeOfRotation );
831       gp_Ax1 AxeOfRotation( TheAxe.Location(), DirAxeOfRotation );
832       Rot.SetRotation( AxeOfRotation, angle );
833     }
834     PSection.Transform(Rot);
835   }
836   // modified by NIZHNY-MKK  Wed Oct  8 15:03:21 2003.END
837
838   if (WithTranslation) { 
839     P.ChangeCoord().SetLinearForm(-1,PSection.XYZ(), 
840                                   V.XYZ());
841   }
842   else {
843     P.SetCoord(0., 0., 0.);
844   }
845
846   gp_Ax3 Saxe(P, gp::DZ(), gp::DX());
847
848   // Transfo
849   gp_Trsf Tf;
850   Tf.SetTransformation(Saxe, Paxe);
851
852   if (WithCorrection) {
853     // modified by NIZHNY-MKK  Fri Oct 17 15:26:36 2003.BEGIN
854     //     Standard_Real angle;
855     //     gp_Trsf Rot;
856     
857     //     if (!isplan)
858     //       throw Standard_Failure("Illegal usage: can't rotate non-planar profile");
859     
860     //     gp_Dir ProfileNormal = TheAxe.Direction();
861     //     gp_Dir SpineStartDir = Paxe.Direction();
862     //     if (! ProfileNormal.IsParallel( SpineStartDir, Precision::Angular() ))
863     //       {
864     //  gp_Dir DirAxeOfRotation = ProfileNormal ^ SpineStartDir;
865     //  angle = ProfileNormal.AngleWithRef( SpineStartDir, DirAxeOfRotation );
866     //  gp_Ax1 AxeOfRotation( TheAxe.Location(), DirAxeOfRotation );
867     //  Rot.SetRotation( AxeOfRotation, angle );
868     //       }
869     // modified by NIZHNY-MKK  Fri Oct 17 15:26:42 2003.END
870
871     Tf *= Rot;
872   }
873
874   return Tf;
875 }
876
877 //===============================================================
878 // Function : Section
879 // Purpose :
880 //===============================================================
881  Handle(Geom_Curve) GeomFill_SectionPlacement::
882  Section(const Standard_Boolean WithTranslation) const
883 {
884   Handle(Geom_Curve) TheSection =  
885     Handle(Geom_Curve)::DownCast(mySection->Copy());
886   TheSection->Transform(Transformation(WithTranslation, Standard_False));
887   return TheSection;
888 }
889
890
891 //===============================================================
892 // Function :
893 // Purpose :
894 //===============================================================
895  Handle(Geom_Curve) GeomFill_SectionPlacement::
896  ModifiedSection(const Standard_Boolean WithTranslation) const
897 {
898   Handle(Geom_Curve) TheSection =  
899     Handle(Geom_Curve)::DownCast(mySection->Copy());
900   TheSection->Transform(Transformation(WithTranslation, Standard_True));
901   return TheSection;
902 }
903
904 //===============================================================
905 // Function :SectionAxis
906 // Purpose :
907 //===============================================================
908  void GeomFill_SectionPlacement::SectionAxis(const gp_Mat& M,
909                                              gp_Vec& T,
910                                              gp_Vec& N, 
911                                              gp_Vec& BN) const
912 {
913   Standard_Real Eps = 1.e-10;
914   gp_Dir D;
915   gp_Vec PathNormal;
916   GeomLProp_CLProps CP(mySection, SecParam, 2, Eps);
917   if (CP.IsTangentDefined()) {
918     CP.Tangent(D);
919     T.SetXYZ(D.XYZ());
920     T.Normalize();
921     if (CP.Curvature() > Eps) {
922       CP.Normal(D);
923       N.SetXYZ(D.XYZ());
924     }
925     else { 
926       // Cas ambigu, on essai de recuperer la normal a la trajectoire
927       // Reste le probleme des points d'inflexions, qui n'est pas 
928       // bien traiter par LProp (pas de recuperation de la normal) : 
929       // A voir...
930       PathNormal.SetXYZ(M.Column(1));
931       PathNormal.Normalize();
932       BN = T ^ PathNormal;
933       if (BN.Magnitude() > Eps) {
934         BN.Normalize();
935         //N = BN ^ T;
936       }
937       N = BN ^ T;
938     }
939   }
940   else { // Cas indefinie, on prend le triedre 
941         // complet sur la trajectoire
942     T.SetXYZ(M.Column(3));
943     N.SetXYZ(M.Column(2));
944   }
945   BN = T ^ N;
946 }
947
948
949 //===============================================================
950 // Function :Choix
951 // Purpose : Decide si le couple (dist, angle) est "meilleur" que
952 // couple courrant.
953 //===============================================================
954  Standard_Boolean GeomFill_SectionPlacement::Choix(const Standard_Real dist,
955                                                    const Standard_Real  angle) const
956 {
957   Standard_Real evoldist, evolangle;
958   evoldist = dist - Dist;
959   evolangle = angle - AngleMax;
960   // (1) Si la gain en distance est > que le gabarit, on prend
961   if (evoldist < - Gabarit) 
962     return Standard_True;
963
964  //  (2) si l'ecart en distance est de l'ordre du gabarit
965   if (Abs(evoldist) < Gabarit) {
966 //  (2.1) si le gain en angle est important on garde
967     if  (evolangle > 0.5)
968        return Standard_True;
969  //  (2.2) si la variation d'angle est moderee on evalue
970  //  une fonction de penalite
971     if (Penalite(angle, dist/Gabarit) < Penalite(AngleMax, Dist/Gabarit) )
972       return Standard_True;
973   }
974
975   return Standard_False;
976 }
977
978
979
980
981
982
983
984